Demonstration of an optimised focal field with long focal depth and high transmission obtained with the Extended Nijboer-Zernike theory

Abstract: In several optical systems, a specific Point Spread Function (PSF) needs to be generated. This can be achieved by shaping the complex field at the pupil. The Extended Nijboer-Zernike (ENZ) theory relates complex Zernike modes on the pupil directly to functions in the focal region. In this paper, we introduce a method to engineer a PSF using the ENZ theory. In particular, we present an optimization algorithm to design an extended depth of focus with high lateral resolution, while keeping the transmission of light high (over 60%). We also have demonstrated three outcomes of the algorithm using a Spatial Light Modulator (SLM).


Introduction
Engineering an elongated focal spot by manipulating the incident wavefront has been done for a long time. It has been known that by putting an annulus in the lens pupil the depth of focus is increased [1]. In fact, it turned out that if the annulus is made extremely narrow, a diffractionfree Bessel beam is approximated [2]. Another method which has been proposed to extend the depth of focus is by dividing the lens pupil in rings and modulate the phase and amplitude of each ring [3]. These and other methods are discussed in more detail in [4][5][6][7][8][9][10].
In this article, a new method of designing pupil function is proposed. It relies on a result of the Extended Nijboer-Zernike theory, as has been suggested by [11]. By exploiting a one-toone correspondence between the Zernike polynomials which compose the pupil function and functions which compose the focal field [12], the number of degrees of freedom reduces to the number of polynomials with which the pupil is constructed, while the computation of the focal field is done by a quick matrix multiplication, rather than a time consuming diffraction integral. In particular, we present in this work an algorithm to find a pupil function which gives rise to an extended depth of focus for the case of a 0.4 numerical aperture, with low loss (less than 40%) of light intensity.
The relevance of such a focal field with an extended depth of focus is apparent from the many applications of Bessel beams, which range from optical tweezers to certain forms of microscopy to barcode scanners. The solution which is presented here however, has the advantage that light passes through the pupil in more areas than just the outer ring, which is a standard method to create a Bessel beam, and thereby achieving higher transmittance.
Finally, three outcomes of the algrotihm have been obtained experimentally using a phaseonly Spatial Light Modulator (SLM) to shape the pupil field [13].

Theory
In this section we briefly discuss theoretical results required for the optimization algorithm that is later presented in this paper. First we give the definition of the complex Zernike polynomials, in which any pupil function can be decomposed. Then, we explain a result from ENZ-theory [12,14], which relates the pupil Zernike modes Z m n to basic functions V m n in the field in the focal region.

Complex Zernike polynomials
Any complex function defined on the unit disk can be expanded in terms of complex Zernike polynomials where the coefficients β m n can be complex and where the complex Zernike polynomials are defined as where This expansion determines both the phase and amplitude of the pupil function. The Zernike polynomials Z m n are a complete set of orthogonal functions on the unit disk [12].

Result from ENZ-theory
In [12] it is shown that there exists a one-to-one correspondence between the complex Zernike polynomials Z m n (ρ, θ ) which compose the pupil field E 0 (ρ, θ ) and functions V |m| n (r, f ) which compose the electric field in the focal region E(r, φ , f ): where (ρ, θ ) are polar coordinates, (r, φ , f ) cylindrical coordinates, and β m n complex coefficients. An expression for the V m n -functions is given in [15] and is included in the Appendix. The importance of this result lies in the fact that now the intensity distribution can be computed directly using the Zernike coefficients β m n and the precalculated V m n rather than using the input field E 0 (ρ, θ ) and the diffraction integral [16].

The optimization algorithm
In order to find good solutions efficiently, we rely on Eq. (4). This equation states that once one has found β such that is close to a desired pattern, the corresponding pupil function follows immediately from Note that choosing m = 0 implies that we only consider circularly symmetric pupil functions. Since the V m n 's only need to be precalculated once, the intensity I can computed this way significantly faster than using the diffraction integral.
In the algorithm, only the intensity distributions along the lateral axis (x-axis) and along the optical axis (z-axis) are considered. More specifically, let (0, 0) be the Gaussian focal point.
Then I x = I(x, z = 0) is the intensity along the lateral axis in the focal plane and I z = I(x = 0, z) is the intensity along the optical axis. Using a least squares method, the β m n coefficients are found such that I x and I z approximate two target functions I x,target and I z,target as shown in Fig. 1. After that, the transmittance of the pupil function is optimized using the Nelder-Mead algorithm. The algorithm is explained in more detail in the Appendix.
In Fig. 2 the cross sections of the phase and amplitude of the obtained pupil function are shown, accompanied by the intensity profiles in the focal region along the z-axis and x-axis. The focal spot for this pupil is longer in the axial direction than in the aberration-free case, meaning that the depth of focus is increased. Since this pupil involves a quite structured amplitude illumination, we have made some approximations in such a way that the pupil could more easily be implemented in an experimental setup. This is treated in detail in the following section.   Table 1 in the row C2.

Creating binary pupil functions
In order to simplify the experimental implementation, the pupil functions found by the algorithm (for example the one shown in Fig. 2) are made binary. Consider the β m n shown in Table  1. According to theory these coefficients should give pupil functions which give rise to an extended depth of focus and a lateral resolution below the diffraction limit. To construct the binary pupil functions we first defineẼ  where Z m n are the complex Zernike polynomials as defined in Eq. (2). The pupil function for 'Phase' (named so because only the phase is modulated) is obtained using Note that this is well-defined, because using only Zernike modes with m = 0 means thatẼ 0 (ρ) is real. The pupil functions for 'Complex 1' and 'Complex 2' (named so because both amplitude and phase are modulated) are obtained using where it turns out by trial-and-error that t = 0.57 is a good choice for 'Complex 1' and t = 0.55 for 'Complex 2'. The resulting pupil functions which will be assigned to the SLM are shown in Figs. 3(a)-3(c). Note that in the regions where the amplitude should be 0, a phase ramp is added so that light hitting that region of the SLM is tilted away. The corresponding theoretical intensity distributions in the focal region produced by these pupil functions are shown in Figs.
, the pupil function that gives the longest focal depth, 61% of the light is transmitted.

Experiment
In the experiment we used a Holoeye PLUTO liquid crystal Spatial Light Modulator (SLM) to assign the pattern to the laser beam. The specifications of the SLM are given in [17]. A schematic of the experimental setup is shown in Fig. 4. A He-Ne laser with a wavelength of 633nm is used as the source. Its light is coupled to a single-mode optical fiber after which it emerges as a point source (by approximation). The light is then collimated using a lens with a focal length of 80cm and cut off by an aperture. The collimated light passes through the beamsplitter, hits the SLM perpendicularly, and is reflected by the same beamsplitter. The plane of the SLM is conjugated with the plane of the entrance pupil of a 0.4 NA microscope objective using two lenses with a focal length of 30cm. The light in the focal region of the 0.4 NA microscope objective is then imaged onto a CCD camera using a 0.9 NA microscope objective and a tube lens. The 0.4 NA microscope objective can be translated back and forth using a piezo-stage, which allows us to scan different planes of the focal field. An interferometer is used to determine the exact position of the piezo-stage [18]. Several important details of the experimental setup are summarized in Table 2.

Measurement method
First, we need to realize that not all the light that hits the SLM gets modulated by it due to interpixel space and a 'wrong' polarization. Thus in order to separate the modulated from the unmodulated beam an additional phase ramp is added on top of the images shown in Fig. 3 [19].
To view different focal planes, a lens phase is added to the image assigned to the SLM. The strength of this lens phase can be expressed as the Zernike defocus coefficient α 0 2 (note that  Table 1 which are assigned to the SLM. In the black regions the phase is 0. In the grey regions the phase is π. In the regions with a phase ramp the light is tilted away, which corresponds to the amplitude being modulated to 0.    However, if we want to compare measurement results to the theoretical predictions, we first need to find out how the Zernike defocus coefficient α 0 2 relates to the Rayleigh unit R E , i.e. which value for α 0 2 corresponds to a shift of 1 R E . This has been done by scanning through focus without a phase mask, and observing for what value for α 0 2 the intensity has dropped to 1/e times the maximum intensity. This value of α 0 2 corresponds to 1 R E . In Fig. 5 the result of this measurement is shown.

Results and discussion
Once the calibrations are done, we obtained the focal field intensity distributions for the three pupil functions (as in Table 1) and compared them with the theoretical predictions as shown in Fig. 6. In Fig. 7 the lateral width (along the x-axis) of the focal spot for the aberration-free case and for the pupil function 'Complex 2' are compared. Although the measurement results clearly indicate an extended depth of focus of a length similar to the prediction, there are still some  discrepancies. To find out the possible causes, we compared the cases where the through focus scans were made by adding a lens phase, or by moving the 0.4 NA microscope objective. The comparison is shown in Fig. 8. It turns out that these measurements give somewhat different results, but it still does not account for the mismatch with the theory. A possible explanation may be that because the 0.9 NA objective used for imaging is infinity corrected, so whenever a plane other than the focal plane is imaged, aberrations may come into play. This is particularly relevant in our case since the focal depth is around 18 R E . Nevertheless, the experimental results have demonstrated that pupil masks indeed create elongated focal spots and that the lateral resolution at the focal plane is below the diffraction limit.

Conclusion
In conclusion, we used the Extended Nijboer-Zernike theory to obtain a focused field with elongated focal length (up to 18 Rayleigh distances) with diffraction-limited spot size, while keeping the transmittance over 60%. Since this theory is based on functions that can be precalculated, the optimisation parameters are reduced to a limited number of Zernike coefficients that compose the pupil, allowing fast calculation of the focal field. We have found not only a way to create elongated focal spots, but also we have shown more generally that ENZ theory has the potential to be used for pupil engineering, as had been suggested by [11]. We also demonstrated experimentally three pupil functions generated by an SLM that produce elongated focal spots.
⋆ We only consider real Zernike coefficients β . This implies that I(x, y, z) is symmetric in z.
⋆ The cross-sections I(x, z) = I(x, 0) and I(x, z) = I(0, z) are a sufficiently good representation of I(x, z) in the entire positive xz-plane, hence we only consider the cross-sections.
⋆ The transverse coordinates in the focal region are normalized by the Airy radius R A = 0.61λ NA , and z is normalized by the Rayleigh unit R E = nλ NA 2 . R A denotes the lateral width of the focal spot of a focused plane wave, while R E denotes the length of the spot along the optical axis for this field.
We are now able to define our decision variables: • Define x 0 to be the location of the local minimum of I(x, 0) closest to the optical axis.
• Define z 0 to be the location of the point with smallest value where I(0, z) < 0.9.
• Define T to be the transmittance of the pupil: Note that the maximum T is obtained when |E 0 (ρ)| = 1, in which case T = 1.
Using these definitions, the problem can be formulated as follows: Find β such that:

B.2. Finding solutions
In order to find good solutions efficiently, we rely on Eq. (4). This equation states that once one has found β such that provides a good result, the correspondig pupil function follows immediately using

B.2.1. Precalculations
We start by defining grid points on x and z where we want to calculate the focal field intensity, and the points ρ where we want to find the pupil function We then proceed by calculating the V 2k and Z 0 2k for those points If we now define the matrices X, Z and R as then I(x, 0), I(0, z) and |E 0 | 2 are calculated and normalized according to Note that these matrices need to be calculated only once, which makes this method of finding the focal field distribution much more efficient than calculating a diffraction integral repeatedly.

B.2.2. The algorithm
We start the algorithm by defining the target functions I target (x, 0) and I target (0, z) which roughly resemble our desired intensity distribution as in Fig. 1. This is shown in the flowchart in Fig. 9. Then we find β which closely match these targets using the built-in Matlab algorithm lsqnonlin.
The merit function f is defined as: If we feed lsqnonlin a random starting β 0 , it will apply the trust-region-reflective algorithm using f and returns a β 1 . We keep creating random starting β 0 until the output is satisfactory.
Having found an initial solution we proceed by improving x 0 , z 0 , T and D using another Matlab function, fminsearch. Given a starting point, it applies the Nelder-Mead algorithm, a heuristic search method, to find minima. Of course, while improving T , we do not want x 0 or z 0 to worsen. Therefore we use the following scheme: • Given a starting solution β 1 , we calculate x * 0 and z * 0 . • We define a function MT ( β ) which, given a certain β , calculates x 0 , z 0 , and T , and returns So basically x * 0 and z * 0 are constraints for x 0 and z 0 . • We apply the Nelder-Mead algorithm to the function MT ( β ), using the starting solution β 1 . It then returns a new β , say β 2 , with improved T .
The algorithm is summarized in Fig. 9.