1 Introduction↩︎

During fabrication errors are introduced on the surfaces of lenses. It is common to distinguish low, mid and high spatial frequency errors. These are best characterized by considering the spatial frequencies of the errors that are induced in the exit pupil of the imaging system. Many authors consider errors on a particular surface to be low frequency errors when the corresponding errors in the pupil can be accurately represented by Zernike polynomials up to the 37th order [1]. Low frequency errors can give considerable broadening of the point spread function (PSF), causing deterioration of the image quality. When a surface contains high spatial frequency errors of the order of the wavelength, usually a considerable amount of the light is scattered under large angles leading to scattering out of the optical system. This causes a reduction of intensity in the image, however the sharpness of the image is mainly retained. Surface errors intermediate to the low and high spatial frequencies are called mid-spatial frequency (MSF) errors. The corresponding errors in the exit pupil are so high that very high order Zernike polynomials would be needed to capture them. However, these expansions in high order Zernike polynomials typically lead to instabilities and will require special care, as numerical precision will become an issue unless calculated with recursion formulas [2], [3]. Also, increasing the number of Zernike polynomials will give diminishing returns on the quality of the fit [4]. Furthermore, it was discussed that for high-order Zernike polynomials a finer grid may be required than for the measurement data, thus necessitating interpolation of the measurement data before these high-order polynomials can be used. Otherwise, the orthogonality of the Zernike polynomials is not maintained. Without orthogonality, adding terms can change the coefficients of the previous terms, potentially changing the quality of the fit. Due to these challenges, the quantification of the influence of MSF errors on the PSF of an imaging system is a relatively difficult problem.

Several methods to quantify the influence of MSF errors on image quality have been proposed. One of the methods is using perturbation methods as described in [5]. In this method, nominal rays are used, and small path length differences are added or subtracted from rays depending on the surface defect. That is to say that the ray paths are not affected by the structure, but only the phase information is changed. However, it was shown that this method only works for small propagation distances [6]. An alternative method is to implement MSF errors in ray tracing software and also change the intersections, resulting in change of ray direction and path. This introduces other challenges, such as the above-mentioned fitting of surfaces containing MSF errors. Instead of increasing Zernike polynomials for surface fits, it is also possible to use specially designed functions to better represent the surface errors introduced by machining [7], [8]. In this way, a better surface representation can be obtained without significantly increasing the required terms in the surface fit. However, this method does not account for any diffractive effects which can occur for the MSF structures with small feature sizes. Another way of analysing the impact of MSF in incoherent imaging systems is done through analysis of the OTF [9]–[13]. In the current paper, the influence of MSF surface errors on coherent imaging is studied.

To study the influence of MSF errors of a certain surface on the PSF, we follow the method proposed in [6], to first image the surface with MSF errors by the components to the right of it in the optical system. Then, similar to the work [14], we propagate from the exit pupil to the image plane. However, note that we do not apply the MSF in the exit pupil, but instead take them into account in the plane of the geometrical image of the surface under the components to the right of it. Thus, the image of the surface in image space is illuminated by the converging spherical wave from the exit pupil and then this field is perturbed by the MSF errors in the image. This field is further propagated using Fresnel diffraction integrals, and the PSF is evaluated in the image plane.

The aim of the paper is to demonstrate the feasibility of the model by comparing with experiments. In Section 2 we show that for a quite general class of MSF structures, the PSF can be computed analytically. However, when the structures are more complicated, evaluation of the analytical formula becomes too time-consuming and numerical solutions are more appropriate. In Section 3 the experimental setup is described and in Section 4 the predictions of the model are compared with the experimental results.

2 Theory↩︎

Optical systems are normally designed using geometrical optics. In the case of, e.g., a lens with MSF structures on its surfaces we cannot use geometrical optics, because the MSF structures will diffract light. On the other hand, modelling with purely physical optics is also not possible due to computational constraints with regard to memory.

Hence, we try to take advantage of the geometrical optics in modelling the optical system as a whole, but introduce physical optics modelling when calculating the PSF. Therefore, we only introduce the MSF structures in the last step of the computation. A similar method was proposed in [6], but instead of our hybrid approach, only geometrical optics was used in that paper.

An optical imaging system consists of one or more optical components, such as lenses and mirrors. We can generalise such a system as ‘black box’ as shown in Fig. 1, and assume that the propagation within it is governed by geometrical optics [15], [16]. The black box is bounded by the entrance and exit pupil. Diffractive effects are incorporated when propagating from the exit pupil to the image plane. When the imaging system is diffraction limited and does not suffer from aberrations, the field in the exit pupil will be a perfect spherical wave, albeit bounded by the aperture size. In practice, the field exiting the exit pupil will not be a perfect sphere. The deviations from the spherical field are wavefront aberrations that reduce the imaging performance [15], [17], [18].

To take account of diffraction effects caused by the MSF errors, we use the approach introduced in [6], where the MSF surface is assumed to be perfectly imaged by the components to the right of it. However, our method of field propagation consists of Fresnel propagation integrals applied to perturbed fields, instead of propagating the phase errors using geometrical optics to the exit pupil. This way, we also can directly propagate the field to the image plane of the optical system and evaluate the PSF. Our method is illustrated in Fig. 2, where the image formation is disturbed by the MSF structures

2.1 Description of the model↩︎

We first introduce some notations and make some general remarks about the Fresnel propagator. The implicit time dependence is throughout given by \(e^{-i\omega t}\), with \(\omega>0\). The positive \(z\)-direction is the direction in which the light propagates. Suppose that between plane \(z=z_1\) and \(z=z_2\) there is a homogenous medium (usually air) where the wavelength and the wave number of light are \(\lambda\) and \(k\), respectively. If \(U_{z_1}\) is the field in the plane \(z=z_1\), the field in plane \(z=z_2\) is given by

\[U_{z_2}(x,y) = \iint U_{z_1}(x',y')p_{z_1z_2}(x-x', y-y') \text{d}x'\text{d}y', \label{eq46fresnelz1z2}\tag{1}\] where \(p_{z_1z_2}(x,y)\) is the Fresnel propagator kernel from the plane \(z=z_1\) to the plane \(z=z_2\): \[\require{physics} \label{eq:propagator} p_{z_1z_2}(x,y) = \frac{\text{e}^{\text{i}k\qty(z_2-z_1)}}{\text{i}\lambda \qty(z_2-z_1)} \text{e}^{\text{i}\pi \frac{x^2+y^2}{\lambda (z_2-z_1)}},\tag{2}\]

Note that (1 ) is valid both the case \(z_2>z_1\) and \(z_2< z_1\), i.e. for both forward and backward propagation. To simplify notation, we let \({\cal P}_{z_1z_2}\) denote the Fresnel propagation operator from \(z=z_1\) to \(z=z_2\):

\[U_{z_2} = {\cal P}_{z_1z_2}(U_{z_1}). \label{eq46Fresneloperator}\tag{3}\] Let the Fourier transform \({\cal F}\) be defined by \[{\cal F}(f)(\xi,\eta) = \iint f(x,y) e^{-2\pi i (\xi x+ \eta y)} \, \text{d}x \text{d}y. \label{eq46defFFT}\tag{4}\] The convolution theorem of the Fourier transform implies \[{\cal F}(U_{z_2})(\xi,\eta) = {\cal F}(p_{z_1z_2})(\xi, \eta) {\cal F}(U_{z_1})(\xi,\eta). \label{eq46FPz1z2}\tag{5}\] Because \[{\cal F}(p_{z_1z_2})(\xi,\eta) = e^{i k (z_2-z_1) e^{-i\pi \lambda(z_2-z_1)(\xi^2+\eta^2)}}, \label{eq46Fpz1z2}\tag{6}\] it follows that when \(z=z_3\) is any third plane: \[{\cal F}(p_{z_1z_3})(\xi,\eta) = {\cal F}(p_{z_2z_3})(\xi,\eta) {\cal F}(p_{z_1z_2})(\xi,\eta). \label{eq46Fpz1z2z3}\tag{7}\] Then (5 ) and (7 ) imply \[{\cal P}_{z_1z_3}= {\cal P}_{z_2z_3}\circ{\cal P}_{z_1z_2}. \label{eq46Pcomposition}\tag{8}\] This means that in homogenous matter, the propagation of a field from plane \(z_1\) to plane \(z_3\) is the same as propagating first from plane \(z_1\) to \(z_2\) and then from \(z_2\) to \(z_3\).

The following formula for the propagation of the product of two fields \(U_{z_1}\) and \(V_{z_1}\) will be used below and is derived in Appendix 6: \[\begin{align} {\cal P}_{z_1z_2}(U_{z_1}V_{z_1})(x,y) = \frac{e^{\pi i \frac{x^2+y^2}{\lambda(z_2-z_1)}}}{\lambda^2 (z_2-z_1)^2} \iint {\cal P}_{z_1z_2}(U_{z_1})(x-x',y-y')\, e^{-i\pi\frac{(x-x')^2 + (y-y')^2}{\lambda(z_2-z_1)}} \nonumber \\ \cdot {\cal F}(V_{z_1})\left(\frac{x'}{\lambda (z_2-z_1)}, \frac{y'}{\lambda (z_2-z_1)} \right) \, \text{d}x' \text{d}y'. \label{eq46propprod} \end{align}\tag{9}\]

We now turn to the situation shown in Fig. 2 where the \(z\)-axis is the optical axis of the system. We limit ourselves to optical imaging systems that have sufficiently low numerical aperture (NA\(\leq 0.6\)) so that polarisation effects can be neglected and paraxial scalar diffraction theory can be used.

For a point source on the optical axis in object space, the field in the exit pupil \(z=z_1\) in Fig. 2 is given by:

\[U_{z_1}(x,y) = P(x,y)\text{e}^{-\text{i}k\frac{x^2+y^2}{2(z_3-z_1)}}, \label{eq46Uz1}\tag{10}\] where \(z=z_3\) is the Gaussian image plane, \(k=\lambda/(2\pi)\) is the wave number with \(\lambda\) the wavelength in image space and \(P(x,y)\) is the pupil function: \[P(x,y) = \begin{cases} e^{i \Phi(x,y)}, \quad \sqrt{x^2+y^2}\leq a, \\ 0, \quad \sqrt{x^2+y^2} > a, \end{cases} \label{eq46pupil}\tag{11}\] where \(a\) is the radius of the exit pupil. The phase function \(\Phi\) may contain low spatial frequency aberrations and can be expanded into Zernike polynomials. For a system without aberrations: \(\Phi(x,y)=0\) for \(\sqrt{x^2+y^2}\leq a\).

The plane \(z=z_2\) is the image plane of the surface in the optical system with MSF errors. As mentioned above, it is assumed that this image is an in general (de)-magnified perfect image of the surface. Blurring caused by the imaging of the surface by the optical components to the right of it, are secondary effects. Note that the image plane need not be to the right of the exit pupil as in Fig. 2. If it is not, \(z_2<z_1\) and the following formula below are still valid. The field in the \(z=z_2\) plane is obtained by applying the Fresnel propagator to the pupil field (10 ): \[U_{z_2} = {\cal P}_{z_1z_2}(U_{z_1}), \label{eq46fresneloperatorz1z2}\tag{12}\] The field perturbed by the the MSF errors is: \[\tilde{U}_{z_2}(x,y) = \tau(x,y) U_{z_2}(x,y), \label{eq46Utildez2}\tag{13}\] where the transmission function \(\tau\) is related to the surface errors described by the function \(z=h(x,y)\) as: \[\tau(x,y) = \exp\left[ \frac{2\pi i \Delta n}{\lambda}h(x',y')\right], \label{eq46tau}\tag{14}\] with \(\Delta n\) the difference in refractive index between the material of the optical component with the surface with the MSF errors and the surrounding medium.

The PSF is obtained by propagating the field (13 ) to the Gaussian image plane \(z=z_3\): \[\begin{align} \widetilde{\text{PSF}} =\tilde{U}_{z_3} & = & {\cal P}_{z_2z_3}(\tilde{U}_{z_2}) \nonumber \\ & = & {\cal P}_{z_2z_3}( \tau U_{z_2}) \nonumber \\ & = & {\cal P}_{z_2z_3}\left( \tau {\cal P}_{z_1z_2}(U_{z_1})\right). \label{eq46Utildez3} \end{align}\tag{15}\] The tilde indicates that the field contains the effect of the MSF errors. \(U_{z_3}\) (without tilde) is the field of the PSF in the image plane without MSF errors, which is the Airy function if there are no aberrations.

We now rewrite this result in a more appealing form. It is seen that (15 ) is equal to the propagator operator \({\cal P}_{z_2z_3}\) applied to the product of \(\tau\) and \(U_{z_2}\). We use formula (9 ) with \(z_1\) and \(z_2\) replaced by \(z_2\) and \(z_3\), respectively and with \(U_{z_1}\) and \(V_{z_1}\) replaced by \(U_{z_2}\) and \(\tau\), respectively. This gives \[\begin{align} \widetilde{\text{PSF}}(x,y) = \frac{e^{\pi i \frac{x^2+y^2}{\lambda (z_3-z_2)}}}{\lambda^2(z_3-z_2)^2} \iint {\cal P}_{z_2 z_3}(U_{z_2})(x-x',y-y') e^{-\pi i \frac{(x-x')^2+(y-y')^2}{\lambda (z_3-z_2)}}\nonumber \\ {\cal F}(\tau)\left(\frac{x'}{\lambda (z_3-z_2)}, \frac{y'}{\lambda (z_3-z_2)} \right) \, \text{d}x' \text{d}y'. \label{eq46PSFtildeb} \end{align}\tag{16}\] Using (8 ), (10 ) and \[\begin{align} {\cal P}_{z_2z_3}(U_{z_2})(x,y) &={\cal P}_{z_2z_3} \circ {\cal P}_{z_1z_2}(U_{z_1})(x,y) = {\cal P}_{z_1,z_3}(U_{z_1})(x,y) \nonumber \\ &=\frac{e^{ik(z_3-z_1)}}{i \lambda (z_3-z_1)} e^{\pi i \frac{x^2+y^2}{\lambda (z_3-z_1)}}\iint P(x',y') e^{ - 2\pi i\frac{x x' + yy'}{\lambda(z_3-z_1)}}\text{d}x'\text{d}y' \nonumber \\ &=\frac{e^{ik(z_3-z_1)}}{i \lambda (z_3-z_1)} e^{i \pi \frac{x^2+y^2}{\lambda(z_3-z_1)}} {\cal F}(P)\left( \frac{x}{\lambda(z_3-z_1)}, \frac{y}{\lambda(z_3-z_1)}\right) \nonumber \\ &=\text{PSF}(x,y), \label{eq46PSF} \end{align}\tag{17}\] where \(\text{PSF}\) is the point spread function without MSF errors. Substituting (17 ) into (16 ) gives: \[\begin{align} \widetilde{\text{PSF}}(x,y) = \frac{e^{\pi i \frac{x^2+y^2}{\lambda (z_3-z_2)}}}{\lambda^2(z_3-z_2)^2} \iint \text{PSF}(x-x',y-y') e^{-\pi i \frac{(x-x')^2+(y-y')^2}{\lambda (z_3-z_2)}} \nonumber \\ {\cal F}(\tau)\left(\frac{x'}{\lambda (z_3-z_2)}, \frac{y'}{\lambda (z_3-z_2)} \right) \, \text{d}x' \text{d}y'. \label{eq46PSFtilde} \end{align}\tag{18}\] Hence, apart form the quadratic phase factors, the PSF in the presence of MSF errors is equal to the PSF without MSF errors convoluted with the Fourier transform of the transmission function \(\tau\), which represents the MSF errors, evaluated at spatial frequencies \(x/(\lambda (z_3-z_2))\) and \(y/(\lambda(z_3-z_2))\): \[{\cal F}(\tau)\left(\frac{x}{\lambda (z_3-z_2)}, \frac{y}{\lambda (z_3-z_2)} \right) , \label{eq46convol}\tag{19}\] where \(z=z_3\) is the plane where the PSF is evaluated and \(z=z_2\) is the plane with the perfect image of the surface with MSF errors due to the optical component tot eh right of that surface. The function (19 ) is broader when the distance \(z_3-z_2\) between the PSF plane and the plane of the image of the surface with MSF errors is larger and also when the MSF errors contain more higher frequencies. For the special case that there are no other errors than MSF errors, the PSF without MSF errors is the Airy disk: \[\text{PSF}(x,y)= \frac{e^{ik(z_3-z_1)}}{i \lambda (z_3-z_1)} e^{\frac{ik}{2(z_3-z_1)}(x^2+y^2)} a \frac{ J_1\left(2\pi a \sqrt{\frac{x^2}{\lambda^2(z_3-z_1)^2}+ \frac{y^2}{\lambda^2(z_3-z_1)^2}}\right)}{\sqrt{\frac{x^2}{\lambda^2(z_3-z_1)^2}+ \frac{y^2}{\lambda^2(z_3-z_1)^2}}}. \label{eq46Airy}\tag{20}\]

2.2 Some special cases↩︎

Suppose that the surface error corresponds to a sinusoidal surface perturbation in one direction, say the \(x\)-direction, with amplitude \(h_0\) and spatial frequency \(\kappa\): \[h(x,y) = h_0 \cos(2\pi \kappa x). \label{eq46surfharm}\tag{21}\] The perturbation has been chosen symmetric around \(x=0\). The corresponding transmission function (14 ) is then \[\tau(x)= e^{ i k \Delta n h_0 \cos(2\pi \kappa x)} = e^{i \bar{h}_0 \cos(2\pi \kappa x)}, \label{eq46specialtau}\tag{22}\] where the dimensionless \(\bar{h}_0\) is given by \[\bar{h}_0 = k \Delta n h_0. \label{eq46barh0}\tag{23}\] With this special case of MSF surface, the influence of a spatial frequency in the phase error on the PSF can be studied. As we will show, an analytic formula for the modified PSF can be derived. Using the Jacobi-Anger expansion \[\tau(x) = \sum_{m=-\infty}^\infty i^m J_{m}( \bar{h}_0) e^{2\pi m \kappa x}.\] Hence, \[{\cal F}(\tau)(\xi,\eta) = \sum_{m=-\infty}^\infty i^m J_{m}(\bar{h}_0) \delta(\xi-m \kappa)\delta(\eta). \label{eq46Ftauspecial}\tag{24}\] Substitution in (18 ) gives \[\widetilde{\text{PSF}}(x,y) = \sum_{m=-\infty}^\infty i^m J_m(\bar{h}_0) \,\text{PSF}(x-\lambda(z_3-z_2)m\kappa, \, y) \, e^{-\pi i \lambda(z_3-z_2) m^2 \kappa^2} e^{2\pi i m \kappa x} \label{eq46PSFharmonic}\tag{25}\] It is seen that the PSF, in the case of a single frequency MSF phase errors, is a linear superposition of PSFs without MSF errors translated over the distances \(\lambda(z_3-z_2)m\kappa\), with \(m\) an integer. The relative contribution of the term corresponding to a given value of \(m\) depends on the value of \(J_m(\bar{h}_0)\).

Let us now consider a surface with errors given by an even but otherwise arbitrary function of \(x\): \(z=h(x)\). We then have \[\begin{align} h(x)& = &\int_0^\infty a(\xi) \cos(2\pi \xi x) \,\text{d}\xi \approx \sum_{j=0}^{N-1} a(\xi_j) \cos(2\pi j \Delta \xi x) \Delta \xi \nonumber \\ & = & \sum_{j=0}^{N-1} h_j \cos(2\pi j \kappa x), \label{eq46FFTh} \end{align}\tag{26}\] where \(h_j=a(\xi_j) \Delta \xi\) and \(\kappa=\Delta \xi\). The corresponding transmission function is \[\begin{align} \tau(x) & = & \prod_{j=0}^{N-1} e^{i \bar{h}_{j} \cos(2\pi j \kappa x)} \nonumber \\ & = & \prod_{j=0}^{N-1}\sum_{m=-\infty}^\infty i^m J_m(\bar{h}_j) e^{2\pi m j \kappa x} \nonumber \\ & = & \sum_{m_0=-\infty}^{\infty} \ldots \sum_{m_{N-1}=-\infty}^\infty \left[ \prod_{j=0}^{N-1} i^{m_j} J_{m_j}(\bar{h}_j)\right] e^{2\pi i \vec{m}\cdot \vec{\kappa} x}, \label{eq46tauxpanded} \end{align}\tag{27}\] where \(\bar{h}_j=k \Delta n h_j\), \(\vec{m}=(m_0,\ldots,m_{N-1})\) and \(\vec{\kappa}=(0,\kappa,\ldots, (N-1)\kappa)\). Hence, \[{\cal F}(\tau)(\xi,\eta) = \sum_{m_0=-\infty}^{\infty} \ldots \sum_{m_{N-1}=-\infty}^\infty \left[ \prod_{j=0}^{N-1} i^{m_j} J_{m_j}(\bar{h}_j)\right] \delta(\xi-\vec{m}\cdot\vec{\kappa})\delta(\eta). \label{eq46Ftaugneral}\tag{28}\] Substitution into (18 ) implies \[\begin{align} \widetilde{PSF}(x,y) = \sum_{m_0=-\infty}^{\infty} \ldots \sum_{m_{N-1}=-\infty}^\infty \left[ \prod_{j=0}^{N-1} i^{m_j} J_{m_j}(\bar{h}_j)\right] \text{PSF}(x-\lambda(z_3-z_2)\vec{m}\cdot\vec{\kappa}, y)& \\ e^{-\pi i \lambda (z_3-z_2)(\vec{m}\cdot\vec{\kappa})^2} e^{-2\pi \vec{m}\cdot\vec{\kappa} x}.& \end{align} \label{eq46PSFgeneral}\tag{29}\] This formula generalizes (25 ). When N is large, it is however not practical because its evaluation requires too many operations. Instead, the numerical method explained in the next section can be used.

2.3 Numerical computation of (18 )↩︎

Since the integral in (18 ) is a convolution of two functions, namely \[\text{PSF}(x,y) e^{-\pi i \frac{x^2+y^2}{\lambda(z_3-z_2)}}, \label{eq46function1}\tag{30}\] and \[{\cal F}(\tau)\left(\frac{x}{\lambda(z_3-z_2)}, \frac{y}{\lambda(z_3-z_2)}\right), \label{eq46function2}\tag{31}\] the integral is most conveniently computed by multiplying the Fourier transforms of these functions of \((x,y)\), then taking the inverse Fourier transform of the product and finally multiplying the result by the quadratic phase factor in front of the integral in (18 ). The Fourier transform of (30 ) should be computed numerically. The Fourier transform of (31 ) is \[\lambda^2(z_3-z_2)^2 \, \tau\left(\lambda(z_3-z_2)\xi, \lambda(z_-z_2)\eta\right). \label{eq46FFTfunction2}\tag{32}\]

If the grid of the image plane of the MSF and the grid of the image plane of the optical system are mismatched one can resort to computation of the convolution without using Fourier transforms or use interpolation before the Fourier transform.

3 Experimental Method↩︎

The model to evaluate the PSF in the presence of MSF errors as described in Section 2 has been validated with an experimental setup. The MSF errors induce perturbations of the phase of the field. A spatial light modulator (SLM) is used to create these phase perturbations and mimic the MSF errors.

3.1 Setup↩︎

A schematic view of the setup is shown in Fig. 3. The light source is a 10 mW polarised HeNe laser emitting at 633 nm. The light intensity is modulated using two polarisers. The second polariser is fixed to ensure that the polarisation direction is as required by the ‘Holoeye Pluto 2 VIS-096-D’ SLM. This SLM is a mirror with a layer of pixelated liquid crystals. Depending on the applied voltage, the orientation of the crystal axis of a pixel changes the orientation of the optical axis and hence the refractive index experienced by the reflected light. This way, in each pixel, the optical path length is adjusted to correspond to the desired phase change. The SLM has a resolution of (1080×1920) square pixels with a side length of 8 \(\mum\).

A pinhole of 20 \(\mu\)m is imaged onto the camera using a Thorlabs best form lens with focal distance \(f = 50\) mm. Using a beam splitter after the imaging lens, it is possible to place a reflective SLM in the converging spherical wave originating from the lens. Compared to a transmission SLM, the pixels of reflective SLM devices are often smaller, which is preferred in our setup. The aperture of the lens is limited so that we are as close to diffraction limited performance as possible. This is discussed in Appendix 7. The distance between the lens and the SLM is 18.5 cm, and the distance between the SLM and the camera is 32.7 cm.

3.2 Measurements procedure↩︎

The measurements have first been prepared by eliminating inaccuracies as much as possible. For example, the best form lens was chosen to minimise the chance to add extra MSF structures to the system [19], [20]. The phase response of the SLM has been calibrated. Also, the surface shape has been measured because SLM displays are typically not optically flat. The surface of the SLM was characterised by a procedure inspired by the one presented in [21]. The method makes use of a Twyman-Green interferometer to determine the height of the SLM surface.

After characterising the SLM, the next step is to achieve an image of the pinhole on the detector with as few aberrations as possible by aligning the components and introducing an aperture. Components such as the beam splitter, SLM and camera are aligned before adding the lens. Finally, the lens is added and the position is fine-tuned with stages until we cannot observe low spatial aberrations. Important is that we have a large imaging distance to make sure that the Airy disk is large enough and can be measured by the camera. Finally, imaging the PSF is done by turning the first of two polarisers and adjusting the exposure time of the camera so that clear images of the PSF can be taken.

4 Results↩︎

The model is validated by realizing specific MSF errors with the SLM. In the first subsection, a sinusoidal phase in one direction is realized, for which the PSF can be computed analytically with formula (25 ). In the subsequent subsection, a more complicated phase structure is realized.

4.1 Sinusoidal phase structure↩︎

For the sinusoidal MSF structure (22 ), and in the absence of low frequency aberrations, the resulting PSF (25 ) is a linear combination of translated Airy disks. These spots can be identified as diffracted orders of the periodic phase structure on the SLM. The \(mth\) spot is translated over the distance \(\lambda m\kappa(z_3-z_2)\), where \(\kappa\) is the spatial frequency of the sinusoidal structure, and is weighted by the factor \(J_m(\bar{h})\) where \(\bar{h}\) is the maximum of the phase modulation. Each spot contributes to the final perturbed PSF

Each Bessel function has known zeros depending on the value of the argument, e.g., the function \(J_0(h)\) for \(h=2.404\) is zero. This means that, according to the model, we can selectively eliminate orders from the image plane. The simulated structure is not related to a specific MSF structure typically found on optical components, but instead is a structure used to verify if the model works as we expect it to. In Fig. 4 four simulation results are shown using the phase perturbation in (22 ). In the simulations, the spatial frequency is kept constant: \(\kappa=6.25\) \(\text{mm}^{-1}\). The images show the intensity distribution in the normalised image plane, with a gamma correction to improve visibility. In (a) we do not apply a phase perturbation, and in (b) we take a phase perturbation of \(h=\pi/4\). In (c) and (d) we make use of the known zeros of the Bessel function to eliminate the 0th and 1st order by setting \(h\) to \(2.4048\) and \(3.8317\) respectfully.

The simulation results from Fig. 4 were generated using parameters measured in the measurement setup, such as imaging distance, aperture size and camera pixel size, in order to get the same PSF in the measurement and simulation. The measurement shown in Fig. 4 (a') shows the PSF that we were able to obtain in the absence of any applied structures while compensating for the surface height variations of the SLM. In Fig. 4 (b'), we applied the same phase perturbation on the SLM as we had previously simulated with our model. In Fig. 4 (c') and (d') we were able to reproduce the removal of the 0th and 1st order. In 4 (c') we see that the 0th order does not completely disappear. This is likely because the SLM does not modulate perfectly.

For a clear comparison, we consider the horizontal cross-sections through the images of both simulated and measured PSF in Fig. 5. The data shown in Fig. 5 is normalised to the maximum of the measurement with \(h=0\), i.e. for a SLM without phase modulation, but no gamma correction was performed as was done in Fig. 4. The cross-sections show that the simulated and measured intensities match well. We can see that for the measurements where we eliminate the orders, we have an increase in the difference between the simulation and measurement. These differences are, however, small. One of the reasons that there is a difference can be the alignment of the structure of the SLM. In the simulations the structure is always perfectly centred but on the SLM this is not possible.

4.2 More general MSF structures↩︎

With the numerical solution, more general MSF structures can be simulated. A one-dimensional phase function consisting of multiple sine functions added is modelled. This way we work towards modelling real MSF structures in a system where we expect a spectrum of frequencies that make up the MSF.

We can write the phase function of the perturbation as a series. Each term contains a single frequency and a relative amplitude \(h_m\). The structure itself is then finally scaled to again the value \(h\) \[\label{eq:multiple95sin95perturbation} \Phi(x') = h\sum_{m=1}^3 h_m\sin(2 \pi \kappa_m x'),\tag{33}\] We insert the following transfer function into the general expression (18 ): \[\tau_{msf} = \text{e}^{\text{i}\Phi(x')}.\] The result obtained using the distortion of the field is shown in Fig. 6. We used three spatial frequencies, namely: \(6.25, 3.57\) and \(1.39\) \(\text{mm}^{-1}\), and the relative amplitudes of each contribution to the phase perturbation (\(h_m\)) are \(1/2, 1/6\) and \(2/6\) respectively, and the values of \(h\) are the same as in the results of Fig. 4. Just as with the results in Fig. 4 the graphs are all normalised to their respective maximum amplitude.

Compared to the results for a single frequency, we observe larger differences between the measurements and the model. By adding multiple sines with different frequencies, the periodicity of the projected structure changes. A shift of the illumination has less impact when a structure of a single frequency is used. When the three sines are added together, we lose periodicity over the SLM display, and any misalignment will alter the results to a higher degree.

5 Discussion and conclusions↩︎

In this work, we presented a model to simulate the effects of MSF structures on the surfaces of lenses with diffraction integrals. The general mathematical model can calculate the field on the image plane grid. If the convolution is performed without Fourier transformations, the grid in the plane with the MSF structures can be independently defined from the grid in the image plane. This gives more flexibility in calculating the PSF but as a result the computational complexity is higher. It is however possible to parallelise the computation.

The mathematical model has been verified in an experimental setup, where we observed agreement between the measurements and simulations. However, we should take into consideration that the model has been tested for a low NA. Increasing the NA in the setup is not easy because for higher NA the required resolution of the detector can become too large. In future works, this could be circumvented by working with higher NA setups and phase retrieval. The field can be measured in image space in a more convenient plane and then be propagated to the image plane of the optical system.

In the setup, we make use of an SLM to quickly test different structures. The advantages are that we can program the SLM with any custom structures, but the limitations are the limited resolution of structures that can be inserted in the model. For further experiments, more optical components can be inserted in the setup and the location of MSF structures can be changed to investigate the effect that it has on the imaging quality of the setup.

Some of the errors that we observe when testing the numerical model can most likely be attributed to alignment errors. Small variations in measured intensity compared to the simulations are expected to originate from misalignment of the beam on the SLM. The reason we observe this mostly for the experiments with multiple sinuses added is because the period of the combined frequencies is larger than the SLM. Therefore, the effects are much more pronounced when working with a phase perturbation that is a function of multiple structures compared to a single sinus.

The results show that the model can be used to estimate the effects of MSF in simple optical systems. However, the model needs to be tested in more complicated optical imaging systems with higher NA, and in the presence of multiple optical surfaces.

6 Derivation of Eq. (9 )↩︎

To simplify the notation we write \(U\) and \(V\) instead of \(U_{z_1}\) and \(V_{z_1}\), respectively. Furthermore, we will use for conciseness: \(\Delta z = z_2-z_1\). We recall that \[{\cal P}_{z_1z_2}(UV) = (U V) * p_{z_1 z_2}, \label{eq46defPprod}\tag{34}\] where * denotes the convolution. Hence \[{\cal F}[{\cal P}_{z_1z_2}(UV)] = {\cal F}(UV) {\cal F}(p_{z_1z_2})= [{\cal F}(U) * {\cal F}(V) ]\, {\cal F}(p_{z_1z_2}). \label{eq46FPprod}\tag{35}\] Now \[\begin{align} {\cal F}(p_{z_1z_2})(\xi,\eta) &= &e^{ik\Delta z} e^{-i\pi \lambda \Delta z (\xi^2+\eta^2)}\nonumber \\ & = & e^{ik\Delta z} e^{-i\pi \lambda \Delta z [(\xi-\xi')^2+(\eta-\eta')^2]} e^{i\pi \lambda \Delta z (\xi'^2+\eta'^2)} e^{-2\pi i \lambda\delta z(\xi\xi'+\eta\eta')} \nonumber \\ & = & e^{ik\Delta z} {\cal F}(p_{z_1 z_2})(\xi-\xi',\eta-\eta') {\cal F}(p_{z_2 z_1})(\xi',\eta')e^{-2\pi i \lambda \Delta z(\xi\xi'+\eta\eta')} \text{d}\xi'\text{d}\eta'.\nonumber \\ \label{eq46splitting} \end{align}\tag{36}\] Substitution of (36 ) into (35 ) gives \[\begin{align} {\cal F}[{\cal P}_{z_1z_2}(UV)](\xi,\eta) & = \iint {\cal F}(U)(\xi-\xi',\eta-\eta') {\cal F}(V)(\xi',\eta') \text{d}\xi'\text{d}\eta'\,\, \nonumber \\& = e^{ik\Delta z} \iint {\cal F}(U)(\xi-\xi',\eta-\eta') {\cal F}(p_{z_1 z_2})(\xi-\xi',\eta-\eta')\nonumber \\ & {\cal F}(V)(\xi',\eta') {\cal F}(p_{z_2 z_1})(\xi',\eta') \nonumber e^{-2\pi i \lambda \Delta z(\xi\xi'+\eta\eta')} \text{d}\xi'\text{d}\eta'. \nonumber \\ \label{eq46FPprodb} \end{align}\tag{37}\] We take the inverse Fourier transform with respect to the variable \((\xi,\eta)\) for fixed \((\xi',\eta')\) of the integrand of (37 ). The result of this inverse Fourier transform is a function of \((x,y)\): \[\begin{align} &{\cal F}^{-1}[ {\cal F}(U)(\xi-\xi',\eta-\eta') {\cal F}(p_{z_1 z_2})(\xi-\xi',\eta-\eta') e^{-2\pi i \lambda \Delta z(\xi\xi'+\eta\eta')}](x,y) \nonumber \\ & = \iint {\cal F}(U)(\xi-\xi',\eta-\eta') {\cal F}(p_{z_1 z_2})(\xi-\xi',\eta-\eta') e^{-2\pi i \lambda \Delta z(\xi\xi'+\eta\eta')} e^{2\pi i (x \xi + y \eta) } \text{d}\xi \text{d}\eta \nonumber \\ & = \iint {\cal F}(U)(\xi-\xi',\eta-\eta') {\cal F}(p_{z_1 z_2})(\xi-\xi',\eta-\eta') e^{2\pi i[(x-\lambda \Delta z\xi')\xi + (y-\lambda \Delta z \eta')\eta]}\text{d}\xi \text{d}\eta \nonumber \\ & =\iint {\cal F}[{\cal P}_{z_1 z_2}(U)](\xi-\xi',\eta-\eta') e^{2\pi i[(x-\lambda \Delta z\xi')\xi + (y-\lambda \Delta z \eta')\eta]}\text{d}\xi \text{d}\eta \nonumber \\ & = \iint {\cal F}[{\cal P}_{z_1 z_2}(U)](\tilde{\xi},\tilde{\eta}) e^{2\pi i[(x-\lambda \Delta z\xi')\tilde{\xi} + (y-\lambda \Delta z \eta')\tilde{\eta}]}\text{d}\tilde{\xi} \text{d}\tilde{\eta} \,\, e^{2\pi i[(x-\lambda \Delta z\xi')\xi' + (y-\lambda \Delta z \eta')\eta']} \nonumber \\ & = {\cal P}_{z_1 z_2}(U)( x-\lambda \Delta z \xi', y-\lambda \Delta z \eta') e^{-2\pi i \lambda \Delta z (\xi'^2+\eta'^2)} e^{2\pi i (x \xi'+ y \eta')}. \label{eq46Finvintegrand} \end{align}\tag{38}\] Hence, \[\begin{align} {\cal P}_{z_1 z_2}(UV)(x,y) & = {\cal F}^{-1} \circ {\cal F}[{\cal P}_{z_1 z_2}(UV)](x,y)&\nonumber \\ &\begin{aligned} = e^{ik \Delta z} \iint {\cal P}_{z_1 z_2}(U)( x-\lambda \Delta z \xi', y-\lambda \Delta z \eta') e^{-2\pi i \lambda \Delta z (\xi'^2+\eta'^2)} e^{2\pi i (x \xi'+ y \eta')}&\nonumber \\ {\cal F}(V)(\xi',\eta') {\cal F}(p_{z_2z_1})(\xi',\eta') \text{d}\xi' \text{d}\eta' &\nonumber \\ \end{aligned} \\ & \begin{align}= e^{ik \Delta z} \iint {\cal P}_{z_1 z_2}(U)( x-\lambda \Delta z \xi', y-\lambda \Delta z \eta') {\cal F}(V)(\xi',\eta') e^{-2\pi i \lambda \Delta z(\xi'^2+\eta'^2)} &\nonumber \\ e^{-i k \Delta z} e^{\pi i \lambda \Delta z (\xi'^2+\eta'^2)} e^{2\pi i(x \xi'+y \eta')} \text{d}\xi' \text{d}\eta' &\nonumber \\ \end{align} \\ & \begin{align} = \iint {\cal P}_{z_1 z_2}(U)( x-\lambda \Delta z \xi', y-\lambda \Delta z \eta') {\cal F}(V)(\xi',\eta') e^{-\pi i \lambda \Delta z (\xi'^2+\eta'^2)} \\ e^{2\pi i(x \xi'+y \eta')} \text{d}\xi' \text{d}\eta' \\ \end{align} \tag{39}\\ & \begin{aligned}= \frac{1}{\lambda^2 \Delta z^2}\iint {\cal P}_{z_1 z_2}(U)( x-x', y-y') {\cal F}(V)\left(\frac{x'}{\lambda \Delta z},\frac{y'}{\lambda \Delta z}\right) e^{-i\pi \frac{x'^2+y'^2}{\lambda \Delta z}} &\\ e^{2\pi i \left(\frac{x}{\lambda \Delta z} x'+\frac{y}{\lambda\Delta z}y'\right)} \text{d}x'\text{d}y' & \\ \end{aligned} \nonumber \\ & \begin{align} = \frac{e^{\pi i \frac{x^2+y^2}{\lambda \Delta z}}}{\lambda^2\Delta z^2}\iint {\cal P}_{z_1 z_2}(U)( x-x', y-y') e^{-\pi i \frac{(x-x'^2+(y-y')^2}{\lambda \Delta z}} & \\{\cal F}(V)\left(\frac{x'}{\lambda \Delta z},\frac{y'}{\lambda \Delta z}\right) \text{d}x'\text{d}y' & \\ \end{align} \tag{40} \end{align}\] We obtain two notable expressions which are the same as (9 ). The expression (40 ) is used to calculate the PSF as described in Section 2.3. If the grid for the MSF structure is mismatched with the grid in the image plane, it is possible to use the expression (39 ) or use interpolation between the grids.

7 Exit pupil↩︎

Limiting the aperture of the lens to ensure close to diffraction-limited performance was done by placing an iris almost against the lens surface. The distance was measured to be approximately 3.54 mm. The diameter of the iris is 4 mm, but the effective exit pupil is larger because the Airy disk imaged on the camera was smaller than would be the case for an exit pupil of 4 mm diameter. With ray optics, it is straightforward to determine that the exit pupil is located between the first surface of the lens and the pinhole. The exit pupil is effectively 9 mm in diameter. This is in agreement with the size of the Airy disk imaged on the camera.

References↩︎