1. Introduction

The optical response of devices at micro and nano-scales is governed by the local phase and amplitude fields [1,2]. In particular, phase is central to a complete understanding of plasmonic-based structures: light-matter interaction and the near field response of such devices can be greatly enhanced [3,4]. Indeed, for optical antenna and antenna-coupled detectors, a wide range of applications in sensing or non-linear optics, have arisen [5–8]. Specifically, characterization of antenna-coupled detectors, both at IR and visible frequencies, depends on the spatial response mapping of these devices and also on the phase distribution of the incident beam [9]. To characterize such devices a probe beam, (a high quality, almost Gaussian focused laser beam) is scanned across the photosensitive region. Thus, the device under test is located under the beam waist plane, where the wavefront is assumed to be represented by a plane wave. The measured output signal is, in general, the convolution of the sensor’s spatial response and the beam spatial distribution [10]. Furthermore, it has been shown that their characterization is the largest source of uncertainty in this deconvolution process [11]. Therefore, a reliable characterization, both in amplitude and phase, is required.

However, measuring both phase and amplitude of the probe beam at focus is a challenging task. Light sources for these characterization processes and the use of focusing optics yields a beam waist of a few wavelengths in size. Therefore, at IR wavelengths, conventional image-based laser beam characterization is no longer applicable: the available focal plane array detectors have a pixel size comparable to the beam size, if not greater. Consequently low spatial resolution can introduce strong uncertainties in the intensity distribution measurement. On the other hand, phase measurement demands, in principle, an interferometric arrangement to perform a phase measurement at the plane where the detector is going to be probed, which is frequently at focus, as mentioned above. Hence, the same problem arises, and lack of resolution prevents an accurate measurement of the beam wavefront. A well established method to determine phase distributions from irradiance maps is based in the Gerchberg-Saxton method that relates the irradiances at the plane of interest and at the plane where the Fourier transform appears [12]. This strategy has been analyzed, refined and adapted even to include the Fresnel diffraction maps or eliminate the irradiance measurement at the plane of interest using a set of lenses [13–15]. An interesting extension of the diffractive imaging techniques has been exploited to find the phase distribution at a given plane of interest for coherent and incoherent x-ray sources [16, 17]. However, these methods require a propagation of the beam from the plane of interest to another acquisition plane where the irradiance is measured by a CCD or similar array detector. When this propagation is not possible, or an array detector is not available, we need to work with the accessible irradiance distributions around the plane of interest, as it happens in this paper.

Tomographic techniques based on the knife-edge technique have been applied to the determination of the intensity distribution of focused beams [18, 19]. This method avoids the use of image forming devices which are not suitable under low spatial resolution conditions and allows reconstruction of irradiance planes with sub-microm resolution [20, 21]. In this paper we evaluate the goodness of this method. The transport of intensity equation (TIE) was also used to derive the wavefront from these irradiance measurements. As a result, a complete description of the propagating field in the focal region of the lens is obtained, and an analysis of the reliability of the numerical calculation is explained based on the experimental procedure.

This contribution is sectioned as follows: In section 2 the theoretical foundations with the main approximations that will be employed are presented. Section 3 is devoted to the description of the experimental set up. We used a multi angle knife-edge technique to obtain a series of intensity measurements along the path of a laser beam. Section 4 evaluates the goodness of this method and its reliability. Then a tomographic process is performed: the intensity profiles obtained are processed via the Radon transform to form 2D maps of irradiance distribution at planes perpendicular to the beam axis. With these planes of irradiance, measured at known distances, the TIE can derive the phase of the propagating beam. However, before these images are fed into the TIE algorithm, a number of considerations must be taken into account. Our main interest has been to understand and evaluate the different uncertainties of the measurement and processing method. The analysis discussed in section 5 leads to a specific set of polynomials to be employed in the Zernike expansion. Assuming these considerations the results of the wavefronts retrieved are shown in section 6 and the wavefront reconstruction error is evaluated. Finally, a summary of the main conclusions is presented in section 7.

2. Theoretical background

The wavefront of a propagating beam can be recovered using the TIE [22, 23]. Consider the inset of Fig. 1. In a paraxial light beam propagating in a given z-direction the TIE associates the irradiance images at two different planes (I_P+1 and I_P−1), located symmetrically with respect to the plane of interest, with the phase front φ at I_P.

 

Fig. 1 Experimental setup diagram. The inset shows the working region where the beam path is cut by the knife.

Download Full Size | PDF

To solve this equation we have used the method derived by Gureyev and Nugent [24]. A fundamental aspect of this method is that both the data and the solution are expanded into a series of N Zernike polynomials. For instance, the phase is retrieved as decomposition in Zernike terms.

These polynomials are widely recognized as a good estimation of the aberration content of the phase front [25, 26]. The TIE is thus reduced to a set of N algebraic equations,

M is a N × N matrix that depends on the irradiance at the plane P_F where the phase is to be recovered (see inset of Fig. 1), which elements M_ij are defined as [24]

δ_zI corresponds to the axial derivative of intensity. This can be obtained by applying the central finite difference approximation:

That is, subtracting the irradiance of two planes and dividing by the distance separating them, (this last factor 1/2Δz was moved to the constant N_f). This is acceptable as long as Δz is small so that the ray deflections are not too large. That is indeed the case within the focal region of a laser beam, and we will check this condition in our experimental data.

δ_zI is then expanded in a series of N Zernike polynomials, and applying Eq. (2) the N Zernike coefficients of the phase can be retrieved. The reliability of this process depends strongly on the experimental estimation of δ_zI [27,28]. However, the number of polynomials used in a Zernike expansion is typically limited to the secondary order aberrations, around 30 polynomials, which is enough to describe wave aberration functions of optical systems with circular pupils. Nevertheless, we consider that this is not sufficient for our application. We believe there is room for a criterion, and we will based it in our experimental conditions.

3. Experimental set- up

Figure 1 shows the experimental set up. A power-stabilized CO₂ laser (LASY-5 Access Laser) is used as a light source. It has an output continuous power of 5 W maximum at a 10.6 μm wavelength, with a high quality gaussian beam (M² ∼ 1). The beam is directed to a suitable IR lens (ASPH-ZC-50-50 from ISP Optics) with an f = 50.4 mm of focal length. The lens was slightly tilted (∼ 6°) to introduce some aberrations. The knife edge is located at the vicinity of the waist of the beam, and by means of DC-motor stages (PI M-415.DG) it can be moved in any x, y and z-direction. Intensity variations, due to the obstructing knife, are observed in a thermopile based detector (Thorlabs S310C).

Prior to the tomographic process the focal plane P_F, or region of minimum width, has to be located. Using the knife-edge technique, the beam width was measured at a number of different z-locations. A typical signal curve obtained when the knife is moved through the beam is shown in Fig. 2, which corresponds to a measured intensity accumulated profile. An approximated profile of the beam is reconstructed by fitting a Gaussian curve to the derivative of the measured intensity profile. The actual beam irradiance and evolution will be obtained using the tomographic method. The calculation of the parameters of this preliminary measurements provides a value of the minimum beam radius of ω_min = 38 μm. If the beam was rotationally symmetric we could identify this value as the size of the beam waist. However, in our case, the beam could show asymmetry and this parameter is only considered to locate the region where the beam is focused. This is done by obtaining an axial distance Δz, related to the usual definition of the Rayleigh range $\mathrm{\Delta }z=\left(\pi {\omega }_{\mathit{min}}^2\right)/4\lambda$, and we will use this value as the distance between measurement planes.

 

Fig. 2 Normalized intensity measurement (dashed line), and deconvolved signal computed with the PSF of the detector (dotted line). Difference between both curves is plotted below.

Download Full Size | PDF

For a tomographic reconstruction the measurements must be performed at different angles. Thus, the knife is also attached to a rotary stage (New Focus 8410M). This allows a change in its orientation, around z–axis, so that the beam can be cut from any angle θ. We have performed 9 cuts along different orientations with Δθ = 20° separation. Clearly the more cuts the more reliable the tomographic reconstruction will result. We will discuss the influence of the number of cuts on the angular resolution in section 5. A set of normalized measurements at a function of angle are illustrated in Fig. 3(a). These are 9 intensity profiles, arranged in columns (representation also known as sinogram). Once differentiated, these profiles will form a set of projections of the beam intensity. Because the projections at 0° and 180° are mirror images of each other, the measurements are carried out to the last angle increment before 180°. The sinogram clearly shows a strong asymmetry in the beam shape produced by the misaligned lens.

 

Fig. 3 (a) Normalized intensity profiles obtained at different angles (arranged in columns), also called sinogram. (b) Beam irradiance map generated from this sinogram via the inverse Radon transform. Astigmatism is clearly seen.

Download Full Size | PDF

Since the detector is a slow thermal sensor (t_r ≈ 2.3 s rise time) the knife edge must move accordingly. The speed of the knife (v) was thus set to 3.4 μm s⁻¹ and recording each intensity profile took around 4 mins (knife-edge going back and forth). The power meter console (Thor-labs PM100D), allowed a sampling frequency of ∼ 13 Hz, and was connected via USB to a computer where data was processed using MATLAB.

4. Measurement constraints and limitations

The first source of inaccuracy analyzed here is related to the bandwidth of the thermal detector. If the speed of the knife was too high, the detector may not follow: it acts as a low pass filter. Analogously, v sets the smallest feature (Δρ) that can be resolved. We have chosen a conservative figure Δρ = t_r · v ≈ 7.8 μm as the smallest distance interval that will cause a change in the measurement system.

This resolution can be improved by reducing speed at a cost of increasing measuring time. For high v the bandwidth of the detector sets the limit. We have derived the point spread function (PSF) of the detector by measuring the response to a step impulse. Then, by applying a deconvolution procedure to the PSF and a measured intensity profile, we obtained a restored signal. For these calculations we used the Lucy-Richardson algorithm in an iterative process that converges to the maximum likelihood solution [10]. The restored, deconvoluted signal and the measured signal are plotted in Fig. 2. The difference between these two plots is also represented. This difference shows a relative error of ∼ 1%, and corresponds with an uncertainty in the beam width of about 0.5 μm, well below Δρ.

In order to obtain the intensity profile of the beam, a derivative process is applied to the accumulated intensity measured profiles. This process will unavoidable add high frequency noise. We will address this issue in the next section.

The intensity profiles obtained at different angles form a set of projections. They are used to reconstruct a 2D map via the inverse Radon transform [29]. For the tomographic image reconstruction we performed an unfiltered back-projection operation. The ramp filter used in [20] increases linearly with frequency, thus making it susceptible to high frequency noise. Filtered back-projection algorithms also use a square window to suppress high frequencies (making the response more of a low-pass filter). That achieves noise reduction at the expense of spatial resolution along the radial axis [30]. In this paper we keep all the spectral information and will deal with high frequency noise, whichever its origin, in the next section.

Figure 3(b) shows a 2D irradiance map reconstructed from the data shown in the sinogram of Fig. 3(a). It corresponds to the beam irradiance at a location close to the focal plane in a square region of 400 μm² with a step grid size of ∼ 1.4 μm. The elliptic shape of the beam is due to astigmatism. Some artifacts from the transform, as rays emanating from the center, can be observed.

5. Zernike expansion considerations

An important property of the Zernike polynomials is that they can be written as a product of a radial and an azimuthal part [31]

As we mention before, Eq. (2) is based on an expansion in Zernike polynomials. Since the polynomials are slow variant functions, a small amount of them may not provide the spatial spectral components necessary to reconstruct the wavefront we wish to represent. On the other hand increasing the number of terms can simply propagate high frequency noise generated in different processing stages, such as the artifacts generated by the inverse Radon transform. Thus, it is apparent that a compromise must be made in order to preserve the information that we pass to the Eq. (2) without introducing unnecessary terms. We describe two criteria that sets limits to the number of Zernike polynomials relevant to the expansion or, in other words, upper bounds for n and m.

5.1. Radial limit

The first constraint concerns the radial part of the Zernike polynomial. An important property of these polynomials is that the normalization has been chosen so that, for all permisible values of n and m, $R_n^m(1)=1$. We will therefore express the radial coordinate normalized by the aperture radius of our measurements R. Then, since the normalized minimum resolvable spatial feature is Δρ/R the associated angular spatial frequency is

Consequently we will neglect Zernike polynomials which spectral components fall out of the range [0, k_max]. The radial part of the Fourier transform of the Zernike polynomials (𝔕_n(k)) is a function that depends on the nth-order Bessel function of the first kind [32]. This spectral density have the shape of a bandpass spatial filter which central frequency increases with n. Therefore high order Zernike polynomials will contribute to carry high frequency noise. Our criterion is based on including only the terms necessary to describe the data.

At this point we define a ratio that provides a normalized power density of |𝔕_n(k)| in the mentioned range as a function of the radial degree:

For a given n, β is a measure of the amount of power spectral density that lies within the range, see Fig. 4. Consider a particular k_max, β decreases with increasing n, because the central frequency of the Fourier transform lies further. Alternatively, as k_max increases, the spectral density of more polynomials lay within the range and therefore the power ratio β increases. Four cases are plotted to illustrate this behavior.

 

Fig. 4 Evolution of β as a function of the radial polynomial degree n (evaluated within the range 0 < k < k_max). For the experimental conditions described in this paper k_max = 25.6. The β parameter is above 99% until n = 21.

Download Full Size | PDF

We will discard polynomials for which the power ratio β falls 1% (a −0.45 dB change). To some extent, this threshold is taken arbitrarily: due to the high spatial resolution, k_max is large, and the curves fall slowly. Since our measurement system is slow no information is carried by high spatial frequencies. Therefore including high order polynomials in the expansion becomes less and less relevant. More importantly, if we employ too few polynomials we are filtering the signal: smoothing of the maps and loss of contrast will result.

Zernike polynomials are typically used to expand the wavefront aberration of optical systems. Consequently, low order zernike polynomials are often sufficient to describe aberrations such as tilt, defocus, astigmatism or coma. However, in other applications, where more complicated surfaces are to be described, such as the topography of mechanical samples [33], a larger set of polynomials may be required. In our case the characteristics of the irradiance maps naturally requires higher order polynomials. According to Eq. (6) our experimental conditions lead to k_max = 25.6. In this case β is above the predefined threshold (β ≥ 0.99), for n_max = 21 (see Fig. 4). Therefore polynomials to the radial order 21 will be employed in the expansion. That accounts for the first 253 polynomials in the Noll notation [25].

5.2. Azimuthal limit

The second constraint is related to the axial symmetry of the Zernike polynomials. Eq. (5) shows that the angular part depends on the harmonic functions sin and cos. Therefore the polynomials have an angular period of 2π/|m| radians and this satisfies the requirements that rotating the coordinate system by an angle 2π/|m| does not change the form of the polynomial. Consequently high |m| order polynomials will show strong angular variations.

The tomographic reconstruction already implies a collection of knife-edge measurements angularly spaced (Δθ). This number of measurements is directly related to the angular variations that will be detected, and it therefore sets the maximum angular frequency of the reconstructed maps. If we perform m′ cuts in a tomographic knife-edge based process we will have measurements angularly spaced Δθ = π/m′. On the other hand, in a direct application of the Nyquist criterion, for a proper signal reconstruction it is required at least two (angular) samples (2Δθ) per period 2π/|m|. Or, in the limit

Therefore |m|_max = m′, i.e. the maximum azimuthal index that should be employed in the expansion corresponds with the number of knife-edge cuts in our measurements.

Another way of seeing this relation follows from the back-projection process. The inverse Radon transform is based on the Fourier slice theorem [30]. This states that the complete back-projected image, resulting from a set of parallel-beam projections at Δθ angular variations, is computed from the discrete 1D Fourier transform of each projection. In the 2D spatial frequency domain that implies that these 1D transforms are arranged in lines at Δθ incremental angles. The final image is obtained by integrating all the inverse 1D transforms. By considering polynomials with angular frequencies higher than the ones present in the 2D transform we may be introducing signals that were not there to begin with.

By merging together the radial and azimuthal criteria we can select a set of Zernike polynomials that belongs to the accesible expansion in a reliable manner. Those polynomials having a larger value of their radial or azimuthal order fall beyond the measurement capabilities of the experimental set-up.

We show an example of the consequences of applying both conditions. Consider a set of measurements that lead to n_max = 4 and |m|_max = 2. If we arrange the Zernike polynomials in a triangle, as illustrated in Fig. 5, the vertical position sets its n order and the horizontal its m order. The radial constrain sets limits to the height of this triangle and the azimuthal condition establishes the maximum width. Only the functions within these limits, drawn in a gray box, should be called upon in the expansion. If we consider the radial constraint we would use a set of polynomials with n ≤ 4, a total number of N = 15. However taking into account the azimuthal limit we would have to remove four of them: $Z_3^{\pm 3}$ and $Z_4^{\pm 4}$ and therefore use a new set of N′ = 11 polynomials.

 

Fig. 5 Representation of Zernikes arranged according to their radial and azimuthal parameters ( $Z_n^m$). As an example, if experimental conditions are defined with n_max = 4 and |m|_max = 2, only the polynomials within the gray area should be employed.

Download Full Size | PDF

Analogously, in our case, we would have to remove from our list of polynomials with n ≤ 21 (N = 253) the ones with |m| > 9. Leaving the series expansion with N′ = 169 terms.

6. Phase retrieval results and discussion

Using the tomographic process explained in the previous sections, 7 irradiance maps were obtained. They were measured at planes where the beam reaches its minimum extent. The distance between planes has been chosen to properly sample the focal region, Δz = 108 μm as explained in section 3. Fig. 6(a) shows the contour irradiance maps with the propagation z-axis illustrated vertically. An astigmatic Gaussian beam clearly evolves from plane 1 to 7. Actually, from the obtained measurement we may conclude that the effect of the inclined lens is to produce a beam propagating as an orthogonal astigmatic beam having their orthogonal waists located at different planes [34]. The axis of the orthogonal beam are rotated about 45° with respect to the x and y-axis. One of the orthogonal beam waist is located around plane 2 and the other is around plane 6. This can be deduced from the evolution of the width, and more importantly, from the evolution of the phase front, that shows an almost cylindrical wavefront at planes 2 and 6, being the axis of the cylinders oriented along orthogonal directions. In between these two planes the wavefront has a saddle shape, as it should be expected for this type of beams, showing an almost circular beam in between the orthogonal waists. Clearly the astigmatism introduced by the misaligned lens dominates over other aberrations that may be present in the beam under study.

 

Fig. 6 (a) Irradiance contour maps reconstructed using the tomographic technique. The planes are separated Δz = 108 μm. (b) Five phase maps can be retrieved from seven irradiance maps.

Download Full Size | PDF

To check the consistency of experimentally measured irradiance data we have evaluated the energy conservation law by integrating δ_zI within the aperture (of radius R = 200 μm). Loss of intensity in a region arises through energy flow across the boundary of the region. For a phase solution to exist this value must be zero. δ_zI is computed from non-consecutive planes (I_P+1 and I_P−1), where P goes from 2 to 6. Integration of δ_zI yields a value below 10⁻⁶ for any P and it has a minimum, below 10⁻⁸, at P = 4. We will consider that negligible.

Then Eq. (2) was employed to retrieve the wavefront associated to this intensity field. We used the subset of N′ Zernike polynomials calculated previously and the parameter N_f = 110. From each group of three consecutive irradiance planes (I_P+1, I_P, I_P−1), a phase map could be generated. Hence a total number of five phase maps were recovered and they are plotted in Fig. 6(b). For instance, to retrieve the phase map at plane 2, the irradiance maps from planes 1, 2 and 3 were employed. Astigmatism clearly governs its propagation and change of curvature from positive to negative can also be seen as the wavefront crosses plane 4. Fig. 7(a) shows the phase map recovered at this plane.

 

Fig. 7 (a) Phase map reconstructed at the focal plane using N′ Zernike coefficients. (b) Difference between the phase map reconstructed with N′ and one with N (magnitude in radians).

Download Full Size | PDF

To evaluate the impact of a particular Zernike expansion in the retrieval algorithm we have compared the wavefront reconstructed with a N′ expansion set of polynomials, plotted in Fig. 7(a), and the one produced with the N set (at the same plane). Fig. 7(b) shows the difference between both wavefronts. Clearly the extra number of polynomials has an impact in the phase retrieval. The highest differences are in the order of 10%, and they are particularly located at the edges of the aperture. This supports the idea that, as we explained in section 5, the extra Zernike polynomials with |m| > |m|_max, (polynomials with strong angular changes), introduce more uncertainty in the final recovered wavefront.

Finally, it should be emphasized that expanding this approach to high-numerical apertures systems is possible. Indeed, vectorial effects could be describe within this experimental framework [35].

7. Conclusions

Characterization of light beams, both in amplitude an phase, is required to test nanophotonic devices such as optical antennas. This process becomes an issue when the width of the beam is of the same order than the pixel size of an imaging system, which happens differently depending on the available detection technology at the corresponding spectral range. Moreover, a reliable recovery of the wavefront typically reveals new problems when setting interferometric measurements. As a consequence, alternative methods have been developed.

In this contribution we have addressed measurements in the IR range and a tomographic method has been used to reconstruct the irradiance map. The tomography is obtained from multi-angle knife-edge scans across the beam. Then, the phase has been retrieved using a solution of the TIE, which is based on a truncated expansion of both the irradiance and the phase information in a basis of Zernike polynomials.

The results show how it can be recovered both phase and modulus of a laser beam at the focal region of a lens. We have illustrated the method using a lens that has been tilted to produce an orthogonal astigmatic beam. The orthogonal beam waists appear at different planes along the propagation axis. The obtained phase and irradiance maps are in accordance with this situation and a complete characterization of the electric field at this region was presented.

Furthermore, we have shown that there is a subset of polynomials that reduces numerical errors in the wavefront retrieval process. Thus, we have presented two criteria based on the radial and azimuthal properties of the Zernike polynomials, and related them with the experimental conditions: Along the radial direction, the limitations are linked to the rise time of the detector and the speed of the knife used in the experiment, showing a trade-off between the time necessary to obtain the measurement and its resolution. These parameters set a spatial frequency range that defines limits to the n-order of the polynomials employed. Alternatively, the azimuthal resolution is directly related with the number of knife-edge scans that are considered in the measurements. Both conditions define a set of polynomials to be used in a reliable expansion. The reconstructed phase front recovered should be free from the influence of artifacts or other high frequency noise that may come from the measurement process, and that could be propagated by an oversized set of polynomials.

Although we have worked this approach out under paraxial conditions, further extension to high numerical aperture systems are feasible, hence putting forward an experimental method to ease the characterization of antenna-coupled detectors and, in general, nanophotonic devices.