Flattening axial intensity oscillations of a diffracted Bessel beam through a cardioid-like hole

We present a new feasible way to flatten the axial intensity oscillations for diffraction of a finite-sized Bessel beam, through designing a cardioid-like hole. The boundary formula of the cardioid-like hole is given analytically. Numerical results by the complete Rayleigh-Sommerfeld method reveal that the Bessel beam propagates stably in a considerably long axial range, after passing through the cardioid-like hole. Compared with the gradually absorbing apodization technique in previous papers, in this paper a hard truncation of the incident Bessel beam is employed at the cardioid-like hole edges. The proposed hard apodization technique takes two advantages in suppressing the axial intensity oscillations, i.e., easier implementation and higher accuracy. It is expected to have practical applications in laser machining, light sectioning, or optical trapping.


Introduction
The nondiffracting beam exhibits a strict propagating invariance property in free space with good axial intensity uniformity as well as stable transverse resolution [1]. Therefore, it attracts intensive research interest due to its wide applications in laser machining [2], light sectioning [3], interferometry [4], optical trapping [5,6], microscopy [7] and so on. Durnin et al. have demonstrated mathematically that the three-dimensional infinite-sized Bessel beams are exact solutions to the scalar wave equations [8,9]. However, it is physically impossible to generate an infinite-sized Bessel beam, in that it is not square integrable and contains infinite energy. Approximate Bessel beams were created in experiments by using an annular slit in the back focal plane of a lens [10][11][12], holographic techniques [12][13][14][15], axicons [16,17], or other ways [18]. In many previous papers, it was reported that the Bessel beam diffracted from a finitesized circular hole presented prominent axial intensity oscillations [19][20][21]. A sudden cut-off of the incident field leads to the strong diffraction effect at the hole edges. This problem was generally solved by covering a gradually absorbing mask. Almost complete suppression of axial intensity oscillations was achieved through apodization of the incident Bessel beam with flattened Gaussian or trigonometric apodized functions [20,22,23].
Although the gradually absorbing apodization technology presents a flat axial intensity profile theoretically, it faces two problems in practical applications. Firstly, the amplitude modulation is usually realized by using the spatial light modulator (SLM) [24][25][26][27][28], but the transmittance is difficult to be adjusted exactly the same as the expected value. If each pixel on the SLM brings about a transmittance error, the obtained intensity distribution will deviate much away from the ideal one. Secondly, it would be difficult that both the total size and spatial resolution of the SLM meet requirements simultaneously so that the apodized field drops very smoothly around the aperture edges. Based on these two considerations, a question is raised. Can we flatten the axial intensity oscillation of a finite-sized Bessel beam by a hard apodization (HA) technique? Here the HA technique means that there only exists a transmitted hole on the input plane, with other parts being opaque. Gray-level transmittance is no longer used for relieving the experimental difficulty and heightening the achieved intensity accuracy. For comparison, we refer to the gradually absorbing apodization technology as the soft apodization (SA) technique.
This paper is organized as follows. In Section 2, the design principles are described in detail, and the boundary formula of the transmitted hole is deduced analytically. For demonstrating the validity of our proposed strategy, in Section 3 the parameters are given and wave propagating behaviors are simulated by the complete Rayleigh-Sommerfeld method [29]. Physical explanations are also delivered. Section 3 consists of two subsections. In subsection 3.1, axial intensity distributions are calculated, and the intensity uniformity is characterized. In subsection 3.2, intensities on several cross-sectional planes within the beam propagation range are evaluated. A brief conclusion is drawn in Section 4 with some discussions.

Boundary formula of the transmitted hole for hard apodization
Let's first consider the diffraction of a Bessel beam from a finite-sized circular hole. The input plane is the xy plane, which situates at z 0 = 0 m. The incident Bessel beam propagates along the z-axis, whose amplitude distribution satisfies the J 0 (β r) function, where J 0 is the zeroth order Bessel function; β is the propagation constant; r = x 2 + y 2 denotes the radial distance from an arbitrary source point (x, y) to the origin on the input plane.
In the transmitted region (z > 0), the field at any point (x ′ , y ′ , z) is calculated by the complete Rayleigh-Sommerfeld method as [29] E(x ′ , y ′ , z) = kz i2π where k = 2π/λ represents the wave vector and λ represents the incident wavelength; Ω denotes integral region of the circular hole; ρ = (x ′ − x) 2 + (y ′ − y) 2 + z 2 is the distance from the source point (x, y, 0) on the input plane to the observation point (x ′ , y ′ , z). E 0 (x, y, 0) = J 0 (β r) = J 0 (β x 2 + y 2 ) stands for the incident field. The intensity at the observation point is given by I(x ′ , y ′ , z) = |E(x ′ , y ′ , z)| 2 , where |E| represents magnitude of the complex number E. On assuming x ′ = y ′ = 0, we will obtain the axial intensity distribution of the Bessel beam diffracted from the circular hole. In this case, the axial field can be simplified in polar coordinates as follows: where R is the radius of the circular hole; E 0 (r, 0) = J 0 (β r) represents the incident Bessel beam on the input plane z=0; r = x 2 + y 2 denotes the radial distance on the input plane. It is well known that strong axial intensity oscillations are encountered due to a hard truncation of the incident field. For suppressing the axial intensity oscillations, in Ref. [20] Cox et al. introduced a trigonometric apodized function as follows where T (r) represents the amplitude transmittance function on the input plane; ε is a smoothing parameter ranging from 0 to 1. Under this SA condition, the transmittance gradually decreases from 1 to 0 within an annular region r ∈ [εR, R], as displayed in Fig. 1(a). The green and pink dashed curves depict two circles with radii of εR and R, respectively. It is seen from Fig.  1(a) that the color gradually becomes from red to blue when the radius increases from εR to R. Correspondingly, the transmittance monotonically decreases from 1 to 0. For calculating the transmitted axial intensity, Eq. (2) is used where the input field now becomes E 0 (r, 0) = J 0 (β r)T (r). By using this SA technique, they achieved a flat axial intensity distribution of the diffracted Bessel beam, as seen in Figs. 4(a) and 4(b) in Ref. [20]. In the following, we will discuss our strategy. From Eq. (2), the field at any axial point (0, 0, z) is a diffraction superposition of all the source points on the input plane. Since there is only one integral variable r in Eq. (2), the optical system is rotationally symmetric. Namely, in the input plane each source point on a definite circle with the same radius r has the same incident field E 0 (r, 0) = J 0 (β r), contributing equally to a given axial observation point. Therefore, we can make such a equivalent change. For instance, in Fig. 1(a), if the transmittance T (r 0 ) equals 0.5 for some radius r 0 , we may let the incident field passes one half while the other half is blocked at radius r 0 . It is easy to understand that we should get the same axial intensity distributions in both cases. For other transmittance values, the only difference lies in the filling ratio of the transmitted part. Now we turn to explore the general shape that the transmitted hole should be like. In Fig.  1(b), the shadowed region represents the equivalent transmitted part for a given radius r within r ∈ [εR, R]. The filling ratio F(r) should fulfill the following equation After mathematical transformation, we may write Eq. (4) explicitly as follows All the points satisfying Eq. (5) form a closed curve, as drawn by the cardioid-like cyan solid curve in Fig. 1(c). From the above analysis, if the incident Bessel beam passes through a cardioid-like hole while it is prevented outside, we can expect the same flat axial intensity distribution as that in Ref. [20]. Figure 1(c) shows the transmittance distribution on the input plane under the proposed HA condition, indicating the transmittance being 1 and 0 inside and outside the region encircled by the cardioid-like curve, respectively.

Propagating behaviors of the apodized Bessel beams
In this section, we will perform numerical simulations to prove the validity of the proposed HA technique. Parameters are chosen as follows. For the incident Bessel beam, the propagation constant is β = 10 4 m −1 with wavelength of λ = 0.5 µm. The circular hole has a radius of R = 50 mm.  In this subsection, the axial intensity distribution of the diffracted Bessel beam is investigated. Firstly, by assuming E 0 (r, 0) to be J 0 (β r) in Eq. (2), we can obtain the axial intensity distribution of the diffracted Bessel beam from a finite-sized circular hole without apodization, as shown by the black curve in Fig. 2(a). It is seen from Fig. 2(a) that the axial intensity oscillates with a relative error of about ±7%. Secondly, axial intensity profiles are flattened by the SA and HA techniques to prove their equivalence, as shown by the blue solid and red dashed curves in Fig. 2(b), respectively. In both cases, we choose the same smoothing parameter ε 0 as 0.5. Under the SA condition, the incident field on the input plane in Eq. (2) is given by E s 0 (r, 0) = J 0 (β r) × T (r), where T (r) is given by Eq. (3). The integral region Ω 1 is the circular hole. Under the HA condition, since the cardioid-like hole is rotationally asymmetric, axial intensity should be calculated from Eq. (1), where the incident field on the input plane is E 0 (x, y, 0) = J 0 (β x 2 + y 2 ) and the integral region Ω 2 is the cardioid-like hole determined by Eq. (5). It is demonstrated in Fig. 2(b) that both curves overlap, namely, we have achieved the same flat axial intensity distributions of the diffracted Bessel beam by the previous SA and the proposed HA techniques. Thirdly, through optimizing the smoothing parameter ε in Eq. (5), we have computed the diffracted axial intensity distributions of the Bessel beam enclosed by the cardioid-like curve, as drawn in Fig. 2(c). The smoothing parameter ε ranges from 0.7 to 0.9 for every 0.05. The blue, green, red, cyan, and black curves correspond to different smoothing parameters of 0.7, 0.75, 0.8 ,0.85, and 0.9, respectively. With the increase of the smoothing parameter, the propagation distance is lengthened while the axial intensity uniformity is weakened, as observed in Fig. 2(c). It can be interpreted. When the smoothing parameter becomes larger, the transmittance drops more rapidly from 1 to 0 near the hole edges. Therefore, the diffraction effect is stronger and the axial intensity oscillations appear to be more prominent. When the smoothing parameter ε reaches 1, it degenerates to the case in Fig. 2(a). On taking both the axial intensity uniformity and the propagation distance into account, the smoothing parameter is compromised as 0.8 in the following.

The axial intensity distribution of the apodized Bessel beam
Next, the axial intensity uniformity is quantitatively characterized when the smoothing parameter is 0.8. The flat axial intensity region is defined as the axial domain with intensity relative error less than ±0.01%. Numerical results reveal that the flat axial intensity region extends out to z f = 32.01 m. The axial range expands to 39.15 and 46.42 m when the restriction of intensity relative error is relaxed to ±0.1% and ±1%, respectively. Then, the axial intensity oscillations are gradually strengthened. After the axial intensity reaches its peak value of 1.037 at z f = 48.44m, it drops rapidly towards 0. When the propagation distance is longer than 61.46 m, the axial intensity drops below 1%. The numerical simulations agree well with the results obtained from geometrical optics. Since the Bessel beam can be considered as superposition of plane waves whose wave vectors lie on a cone with an angle θ with respect to the z-axis [10,20], where β = k sin(θ ) = 2π sin(θ )/λ . In geometrical optics, the effective axial beam range is given by z max = R/ tan(θ ) = 62.83 m.

Intensity distributions on several cross-sectional planes for the hard apodized Bessel beam
Transverse position x (mm) Transverse position y (mm) Since the cardioid-like curve in Fig. 1(c) is rotationally asymmetric, the field on crosssectional planes should also be asymmetric. In this case, the transmitted intensity distributions are calculated from Eq. (1) and the integral region Ω is the cardioid-like hole. To see how asymmetric the field is, figure 3 displays the regional intensity deviation on three crosssectional planes. Figure 3(a) presents the intensity deviation ∆I 1 = I 1 (x 1 , y 1 ) − I 0 (x, y), where I 1 (x 1 , y 1 ) represents the intensity distribution on the lateral plane z 1 = 10 m and I 0 (x, y) = |J 0 (β x 2 + y 2 )| 2 is the intensity of the incident Bessel beam on the input plane z 0 = 0 m. Figures 3(b) and 3(c) are the same as Fig. 3(a) except for z 2 = 20 m, and z 3 = 30 m, respectively. A cardioid-like shape of intensity deviation is observed in Fig. 3, which suggests that asymmetry of the cardioid-like hole leads to the intensity asymmetry. With the propagation of Bessel beam, the size of the cardioid-like shape is gradually decreased. The intensity deviation is in the magnitude of 10 −3 from the color bars, indicating the intensity relative error is less than 1%. Accordingly, the diffracted Bessel beam preserves a relatively stable propagation property in a considerably long axial range, although somewhat asymmetric. To see the intensity asymmetry property more clearly, we extract the intensity distributions along the x-axis and the y-axis on the above three cross-sectional planes, as shown in Figs. 4(a), 4(b), and 4(c). The blue solid and red dashed curves correspond to intensities along the x-axis and the y-axis, respectively. It is seen from Figs. 4(a) to 4(c) that the intensity profiles on both axes almost overlap. In order to magnify the intensity difference along the x-axis and the y-axis, we calculate the intensity deviation ∆I x = I i (x i , 0) − |J 0 (β |x|)| 2 and ∆I y = I i (0, y i ) − |J 0 (β |y|)| 2 on the three cross-sectional planes, as shown by the blue and red curves in Fig.  4(d) from top to bottom. I i (x i , 0) and I i (0, y i ) (i=1,2,3) represent the intensities along the x-axis and the y-axis, respectively. The intensity deviation ∆I x or ∆I y presents the difference between the intensity at the observation point and that at the corresponding point on the input plane, indicating the propagation invariance property. It is seen from Fig. 4(d) that the red curves (∆I y ) are symmetric about the origin point while the blue curves (∆I x ) are asymmetric about the origin point. The reason is that the cardioid-like hole is symmetric about the x-axis but asymmetric about the y-axis. Since the maximum intensity deviation is less than ±10 −2 , the diffracted Bessel beam propagates relatively stably within the beam propagation range. If we calculate the asymmetric intensity difference ∆I xy = I x − I y = ∆I x − ∆I y , the intensity asymmetric property will be quantified. The intensity asymmetry does exist since the two curves in Fig. 4(d) are clearly separated. However, the maximum asymmetric intensity difference ∆I xy is only 0.89%. Especially, for applications of Bessel beams we are generally interested in a few lobes near the center. The asymmetric intensity difference ∆I xy is below 0.12% for the central five intensity lobes. Consequently, the intensity asymmetry can almost be neglected in the paraxial region.

Conclusion and discussions
In this paper, we present a new feasible way to flatten the diffracted axial intensity distribution of a finite-sized Bessel beam. Numerical results by the complete Rayleigh-Sommerfeld method have validated our strategy. It is found that the Bessel beam diffracted from a cardioid-like hole maintains a very good propagating invariance property along the axial direction. The boundary formula of the cardioid-like hole is given analytically. Compared with the previous SA technique, the proposed HA technique provides a more feasible choice in practical applications with higher accuracy, since the error only comes from the hole boundary.
In fact, the proposed HA technique is a universal method for suppressing the axial intensity oscillations, by transforming the transmittance distribution into spatial transmitted filling factors. Therefore, for other amplitude transmittance functions like the flattened Gaussian shape, by analogy we can obtain another closed transmitted hole on the input plane. Moreover, the proposed method may be extended to suppress axial intensity oscillations for other kinds of incident waves, such as plane waves or Gaussian waves. It is expected to have practical applications in many optical processing systems.