Analytical Model of Point Spread Function under Defocused Degradation in Diffraction-Limited Systems: Confluent Hypergeometric Function

Abstract

In recent years, optical systems near the diffraction limit have been widely used in high-end applications. Evidently, an analytical solution of the point spread function (PSF) will help to enhance both understanding and dealing with the imaging process. This paper analyzes the Fresnel diffraction of diffraction-limited optical systems in defocused conditions. For this work, an analytical solution of the defocused PSF was obtained using the series expansion of the confluent hypergeometric functions. The analytical expression of the defocused optical transfer function is also presented herein for comparison with the PSF. Additionally, some characteristic parameters for the PSF are provided, such as the equivalent bandwidth and the Strehl ratio. Comparing the PSF obtained using the fast Fourier transform algorithm of an optical system with known, detailed parameters to the analytical solution derived in this paper using only the typical parameters, the root mean square errors of the two methods were found to be less than 3% in the weak and medium defocus range. The attractive advantages of the universal model, which is independent of design details, objective types, and applications, are discussed.

1. Introduction

The point spread function (PSF) is a term used to describe the intensity distribution of the image that is received at the detector plane when the light emitted by a point source passes through an optical system [1,2,3]. This function can be utilized to evaluate the performance of an optical system. The image captured by the optical system can be seen as the convolution of the actual object and the PSF. Thus, the PSF is often used as a physical reference in image processing and computer vision, particularly for image restoration, image denoising, and target detection [4,5,6]. Chen et al. [7,8] employed ray tracing to compute the PSF and utilized the spatial frequency response measurements to optimize the perturbation parameters of the optical system. They successfully developed a virtual image camera that produced imaging simulation results closely matching real-world photography. Furthermore, by integrating a deep learning algorithm, they were able to mitigate optical aberrations arising from processing errors. Researchers from the Chinese Academy of Sciences, Tsinghua University, and the Howard Hughes Medical Institute fused physical priors, such as the PSF of the optical system, with a neural network architecture design to achieve high-performance and high-fidelity microscopic image denoising and super-resolution reconstruction [9]. Additionally, Chan et al. from Carnegie Mellon University [10] have put forth a groundbreaking approach to compressive high-speed imaging. By time-encoding the camera’s PSF, their method provides better optical efficiency and can be applied across a broader spectrum of scene categories. Yoav Shechtman et al. [11] proposed a framework for pupil-plane modulation in 3D imaging applications that require precise localization, based on the PSF of the system. This framework aims to extract maximum physical information about the position of a single nanoscale object from an image. The research group of Zhang Yunhai at the Suzhou Institute of Biomedical Engineering Technology, Chinese Academy of Sciences, has proposed a 2π-DH-PSF using the Fresnel zone approach that can rotate 2π radians, showcasing its superior performance in the 3D localization imaging of nanoparticles [12].

Diffraction-limited systems are those in which the diffraction dominates the image, and the geometric aberrations are of the same order compared to the Airy disk. With high-level design algorithms and the advancement of precision optical and mechanical machining capabilities, an increasing number of diffraction-limited optical systems are being designed and utilized in high-end technology fields such as aerospace, medicine, and astronomy [13,14,15]. Due to the more rigorous requirements regarding imaging quality within diffraction-limited optical systems, physical intelligence-driven image post-processing algorithms are frequently utilized to enhance image quality, in addition to focusing. The implementation of intelligent image-processing techniques, grounded in physical priors, hinges on the precise determination of the PSF for each defocusing position [16]. Therefore, there is an immediate need to investigate a high-precision, low-complexity solution for the PSF.

The most widely used method for calculating the PSF is via the Fourier transform of the wavefront using the fast Fourier transform algorithm, necessitating the modeling of the optical system to obtain the optical path difference through ray tracing [7,17,18,19,20]. While this method is commonly employed for PSF calculations in optical systems, its specific implementation remains challenging. Firstly, the ray tracing method necessitates precise knowledge of the optical system’s detailed parameters, such as material properties, curvature, thickness, spacing, and other intricate lens details, which may be difficult to acquire for various reasons. Additionally, complex optical systems require significant computational resources for the process of ray tracing, which leads to a slow calculation process; analytical solutions can alleviate these problems to a certain extent. Various research papers have explored the analysis and analytical calculation of the PSF. Lommel first provided an analytical expression for the PSF in the case of aberration-free and out-of-focus positioning in 1885 [21,22]. In 1942, Nijboer expanded the aberration functions using Zernike polynomials, which allowed for the analytical evaluation of the PSF of small aberration optical systems [23]. Based on the research of Nijboer, Janssen et al. formulated a general expression using the power-Bessel series, solving the PSF for optical systems with large aberrations [24,25]. In 2019, Miks derived analytical formulas that enable the approximate calculation of the PSF for optical systems affected by defocus and spherical aberration [26].

This paper presents the first derivation of an analytical solution for the defocused PSF in diffraction-limited systems. The solution is expressed as a confluent hypergeometric function with rapid converging series expression. The findings exhibit strong consistency across various defocusing scenarios. Further exploration of this analytical solution can yield insights into the image degradation of diverse design configurations within diffraction-limited systems under defocusing conditions, thereby establishing a unified model for high-definition imaging systems in proximity to the focal point.

2. Fresnel Diffraction and Fourier–Bessel Transform of Converging Spherical Waves

The fundamental function of an incoherent optical imaging system is to transform the spherical wave that emerges from a point source on the object plane through the entrance pupil so that it converges into a focused wave and is directed onto the image plane via the exit pupil. In the case of a diffraction-limited system without any geometric aberrations, these converging waves are spherical. However, if any geometric aberration exists, the wavefront of the converging wave deviates from being perfectly spherical; we will only discuss those cases where such aberrations are absent.

Assuming that a spherical wave undergoes diffraction through a circular exit pupil (coordinated in the ξ-η plane), it converges at point P (0, 0, z₀), as illustrated in Figure 1. When observed in the x-y coordinate system on the observation plane located at a distance δ beyond point P, it results in defocusing. Here, δ represents the amount of defocus, typically where δ $\ll$ z₀.

The propagation effect of a light wave from the exit pupil to the observation plane is Fresnel diffraction. Considering a point denoted as Q (ξ, η) on the exit pupil and a point Q′ (x, y, z) on the observation plane, according to the Fresnel diffraction theory, the field distribution ψ (x, y, z) is given by the following integral:

$$\psi \left(x,y,z\right)=\frac{e^{ikz}}{i\lambda z}{\displaystyle \underset{\Sigma }{\iint }p\left(\xi ,\eta \right)}\exp \left\{i\frac{\pi }{\lambda z}\left[{\left(x-\xi \right)}^2+{\left(y-\eta \right)}^2\right]\right\}d\xi d\eta ,$$

(1)

where p (ξ, η) is the pupil function, k is 2π/λ, and λ is the wavelength.

Working from Figure 1, let r be the distance between point Q and point P on the ideal image plane axis. This yields the following approximate expression:

$$r={\left[z_0^2+\left({\xi }^2+{\eta }^2\right)\right]}^{1/2}=z_0{\left[1+\frac{{\xi }^2+{\eta }^2}{z_0^2}\right]}^{1/2}\approx z_0+\frac{{\xi }^2+{\eta }^2}{2z_0}.$$

(2)

The pupil function is the complex amplitude distribution of converging spherical waves on the pupil surface (see Figure 2), given by:

$$p\left(\xi ,\eta \right)=\{\begin{array}{cc}\frac{e^{-ikr}}{r}\approx \frac{1}{z_0}\exp \left[-i\frac{\pi }{\lambda z_0}\left({\xi }^2+{\eta }^2\right)\right]\hfill & \sqrt{{\xi }^2+{\eta }^2}<a,\hfill \\ 0\hfill & \sqrt{{\xi }^2+{\eta }^2}\ge a,\hfill \end{array}$$

(3)

where $a$ is the radius of the exit pupil.

Substituting the complex amplitude distribution $p\left(\xi ,\eta \right)$ of the light wave on the exit pupil plane into the Fresnel diffraction formula, as shown in Equation (1), leads to complex amplitude distribution on the observation plane (defocused image plane) x-y:

$$\begin{array}{cc}\hfill \psi \left(x,y,z\right)=& \frac{e^{ikz}}{i\lambda z{z}_0}\exp \left[i\frac{\pi }{\lambda }\left(\frac{x^2+y^2}{z}\right)\right]\times \hfill \\ & {\displaystyle \underset{\Sigma }{\iint }\exp \left[-i\frac{\pi \delta }{\lambda z_0\left(z_0+\delta \right)}\left({\xi }^2+{\eta }^2\right)\right]\exp \left[-i\frac{2\pi }{\lambda z}\left(x\xi +y\eta \right)\right]}d\xi d\eta .\hfill \end{array}$$

(4)

According to the rotational symmetry of the system, polar coordinates are used and normalized:

$$\begin{array}{l}\xi =ar\cos \theta ,\eta =ar\sin \theta ,\frac{x}{\lambda z}=\frac{\rho }{a}\cos \varphi ,\frac{y}{\lambda z}=\frac{\rho }{a}\sin \varphi ,\\ C=\frac{e^{ikz}}{i\lambda z{z}_0}\exp \left[i\frac{\pi }{\lambda z}\left(x^2+y^2\right)\right],\kappa =\pi A^2\left(\delta /\lambda \right),A=a/z_0,\end{array}$$

(5)

where A is the numerical aperture, and r and ρ are the normalized dimensionless radius in the spatial domain (exit pupil coordinates) and frequency domain (observation plane coordinates), respectively. Thus, Equation (4) can now be expressed as:

$$\begin{array}{cc}\hfill \psi \left(x,y,z\right)& =\psi \left(\rho \right)\hfill \\ & =C{\displaystyle {\int }_0^{2\pi }\left[{\displaystyle {\int }_0^1\exp \left(i\frac{\kappa }{\left(1+\delta /z_0\right)}{r}^2\right)rdr}\right]}\exp \left[-i2\pi r\rho \cos \left(\theta -\varphi \right)\right]d\theta \hfill \\ & \propto {\displaystyle {\int }_0^{2\pi }\left[{\displaystyle {\int }_0^1\exp \left(i\kappa r^2\right)rdr}\right]}\exp \left(-i2\pi r\rho \cos \theta \right)d\theta \hfill \\ & =\mathfrak{B}\left\{{\displaystyle {\int }_0^1\exp \left(i\kappa r^2\right)rdr}\right\}.\hfill \end{array}$$

(6)

Because the δ/z₀ $\ll$ 1, δ/z₀ is omitted in the first-order approximation, and so is the constant, C. Equation (6) indicates that due to the rotation symmetry of the system, the complex amplitude distribution $\psi \left(\rho \right)$ on the defocusing plane is a Fourier–Bessel transform, denoted by symbol B of the integral ${\displaystyle {\int }_0^1\exp \left(i\pi \kappa r^2\right)rdr}$, with the parameter κ under the first-order approximation.

The integral expression of the zero-order Bessel function is:

$$J_0\left(z\right)=\frac{1}{2\pi }{\displaystyle {\int }_0^{2\pi }{e}^{iz\cos \theta }d\theta },$$

(7)

Substituting it into Equation (6) by setting $z=2\pi r\rho$, one obtains:

$$\begin{array}{cc}\hfill \psi \left(A,\delta ;\rho \right)& =2\pi {\displaystyle {\int }_0^1\exp \left[i\kappa r^2\right]J_0\left(2\pi r\rho \right)rdr}\hfill \\ & =2\pi {\displaystyle {\int }_0^1\exp \left[i\pi A^2\left(\delta /\lambda \right)r^2\right]J_0\left(2\pi r\rho \right)rdr}.\hfill \end{array}$$

(8)

3. Analytical Expression of the PSF

3.1. Zero-Order Approximate (Focal Plane) PSF—Airy Disk

When δ = 0 (with no defocus), the intensity distribution (PSF) corresponding to Equation (8) is the Airy disk:

$$h\left(\rho \right)={\left|\psi \left(\rho \right)\right|}^2={\left[2\pi {\displaystyle {\int }_0^1{J}_0\left(2\pi r\rho \right)rdr}\right]}^2={\left[\frac{J_1\left(2\pi \rho \right)}{\rho }\right]}^2,$$

(9)

Let μ_1,1 be the first zero point of J₁ (x), then:

$$\rho =\frac{{\mu }_{1,1}}{2\pi }=\frac{3.832}{2\pi }=0.61.$$

(10)

It can be established that the radius of the first dark ring of the Airy disk is:

$$x=\rho \frac{\lambda z_0}{a}=0.61\frac{\lambda }{A}.$$

(11)

3.2. First-Order Approximation (Defocus) PSF—Confluent Hypergeometric Function

When $\delta \ne 0$, one can obtain from Equation (8):

$$\begin{array}{cc}\hfill h\left(A,\delta ;\rho \right)& ={\left|2\pi a^{2}C{\displaystyle {\int }_0^1\exp \left(i\kappa r^2\right)J_0\left(2\pi r\rho \right)rdr}\right|}^2\hfill \\ & \propto {\left(2\pi \right)}^2\left\{{\left[{\displaystyle {\int }_0^1\cos \left(\kappa r^2\right)J_0\left(2\pi r\rho \right)rdr}\right]}^2+{\left[{\displaystyle {\int }_0^1\sin \left(\kappa r^2\right)J_0\left(2\pi r\rho \right)rdr}\right]}^2\right\}.\hfill \end{array}$$

(12)

The series expression of the Bessel function of order zero is as follows:

$$J_0\left(x\right)={\displaystyle \sum _{n=0}^{\infty }\frac{{\left(-1\right)}^n}{{\left(n!\right)}^2}}{\left(\frac{x}{2}\right)}^{2n}.$$

(13)

Substituting Equation (13) into Equation (12), setting $u=r^2$ and omitting the unimportant constants before the integral, yields:

$$\begin{array}{cc}\hfill h\left(A,\delta ;\rho \right)& ={\left[{\displaystyle \sum _{n=0}^{\infty }\frac{{\left(-1\right)}^n{\left(\pi \rho \right)}^{2n}}{{\left(n!\right)}^2}\left[{\displaystyle {\int }_0^1\cos \left(\kappa u\right)u^{n}du}\right]}\right]}^2+{\left[{\displaystyle \sum _{n=0}^{\infty }\frac{{\left(-1\right)}^n{\left(\pi \rho \right)}^{2n}}{{\left(n!\right)}^2}\left[{\displaystyle {\int }_0^1\sin \left(\kappa u\right)u^{n}du}\right]}\right]}^2\hfill \\ & ={\left[{\displaystyle \sum _{n=0}^{\infty }\frac{{\left(-1\right)}^n{\left(\pi \rho \right)}^{2n}}{{\left(n!\right)}^2}{I}_C^{\langle n\rangle }}\right]}^2+{\left[{\displaystyle \sum _{n=0}^{\infty }\frac{{\left(-1\right)}^n{\left(\pi \rho \right)}^{2n}}{{\left(n!\right)}^2}{I}_S^{\langle n\rangle }}\right]}^2,\hfill \end{array}$$

(14)

where $I_C^{\langle n\rangle },I_S^{\langle n\rangle }$ can be written as:

$$\{\begin{array}{c}{I}_C^{\langle n\rangle }={\displaystyle {\int }_0^1\cos \left(\kappa u\right)u^{n}du}=\frac{1}{2\left(n+1\right)}\left[F\left(n+1,n+2;i\kappa \right)+F\left(n+1,n+2;-i\kappa \right)\right],\\ I_S^{\langle n\rangle }={\displaystyle {\int }_0^1\sin \left(\kappa u\right)u^{n}du}=\frac{-i}{2\left(n+1\right)}\left[F\left(n+1,n+2;i\kappa \right)-F\left(n+1,n+2;-i\kappa \right)\right],\end{array}$$

(15)

and where $F\left(\alpha ,\gamma ;z\right)$ is the confluent hyper-geometric function [27,28], defined as follows:

$$\begin{array}{cc}\hfill F\left(\alpha ,\gamma ;z\right)& ={\displaystyle \sum _{k=0}^{\infty }\frac{{\left(\alpha \right)}_k}{k!{\left(\gamma \right)}_k}{z}^k}\hfill \\ & =1+\frac{\alpha }{\gamma }\frac{z}{1!}+\frac{\alpha \left(\alpha +1\right)}{\gamma \left(\gamma +1\right)}\frac{z^2}{2!}+\frac{\alpha \left(\alpha +1\right)\left(\alpha +2\right)}{\gamma \left(\gamma +1\right)\left(\gamma +2\right)}\frac{z^3}{3!}+\cdots .\hfill \end{array}$$

(16)

Thus, the analytical solution of the defocus PSF of the diffraction-limited imaging system can be presented as:

$$\begin{array}{l}h\left(A,\delta ;\rho \right)={\left\{{\displaystyle \sum _{n=0}^{\infty }\frac{{\left(-1\right)}^n{\left(\pi \rho \right)}^{2n}}{{\left[2{\left(n+1\right)}^{2}n!\right]}^2}}\left[F\left(n+1,n+2;i\kappa \right)+F\left(n+1,n+2;-i\kappa \right)\right]\right\}}^2+\\ \hspace{1em}\hspace{1em}\hspace{1em}{\left\{{\displaystyle \sum _{n=0}^{\infty }\frac{{\left(-1\right)}^n{\left(\pi \rho \right)}^{2n}}{{\left[2{\left(n+1\right)}^{2}n!\right]}^2}}\left[F\left(n+1,n+2;i\kappa \right)-F\left(n+1,n+2;-i\kappa \right)\right]\right\}}^2.\end{array}$$

(17)

For $\alpha =n+1$ and $\gamma =n+2$, these functions can be expressed as series that converge very rapidly:

$$\{\begin{array}{c}F\left(n+1,n+2;i\kappa \right)+F\left(n+1,n+2;-i\kappa \right)=2\left(n+1\right){\displaystyle \sum _{k=0}^{\infty }\frac{{\left(-1\right)}^k}{\left(n+2k+1\right)}\frac{{\kappa }^{2k}}{\left(2k\right)!}},\\ F\left(n+1,n+2;i\kappa \right)-F\left(n+1,n+2;-i\kappa \right)=-2i\left(n+1\right){\displaystyle \sum _{k=1}^{\infty }\frac{{\left(-1\right)}^k}{\left(n+2k\right)}\frac{{\kappa }^{2k-1}}{\left(2k-1\right)!}}.\end{array}$$

(18)

Substituting Equation (18) into Equation (17) yields a series expression for the PSF:

$$h\left(A,\delta ;\rho \right)={\left[{\displaystyle \sum _{n=0}^{\infty }{\left(-1\right)}^n\frac{{\left(\pi \rho \right)}^{2n}}{{\left(n!\right)}^2}}{\displaystyle \sum _{k=0}^{\infty }\frac{{\left(-1\right)}^k{\kappa }^{2k}}{\left(n+2k+1\right)\left(2k\right)!}}\right]}^2+{\left[{\displaystyle \sum _{n=0}^{\infty }{\left(-1\right)}^n\frac{{\left(\pi \rho \right)}^{2n}}{{\left(n!\right)}^2}}{\displaystyle \sum _{k=1}^{\infty }\frac{{\left(-1\right)}^k{\kappa }^{2k-1}}{\left(n+2k\right)\left(2k-1\right)!}}\right]}^2.$$

(19)

When A = 0.242 and δ/λ = 0,6,12,18,24,30, the PSF is shown in Figure 2. It can be seen that by increasing the defocus, the brightness of the central point of the PSF decreases, and more energy diffuses to the sidelobes, resulting in increased spreading.

In evaluating the convergence speed of the series equation, Equation (19), for a system with typical parameter values, seven sampling values on the PSF curve are taken into account, as depicted in Figure 3. The abscissa represents the truncation term number N. All terms from n = 0 to n = N are retained and are included in the summation, while terms with n > N are omitted. The ordinate represents the corresponding values of the PSF curves with labels. The various curves in the graph represent the convergence status of the functions at the corresponding positions in the small graph. Evidently, all the curves converge rapidly and approach constants when N > 12, indicating that the confluent hypergeometric function series converges smoothly and quickly.

3.3. Optical Transfer Functions for Diffraction-Limited Systems with Defocus

The optical transfer function (OTF) is a crucial, complex-valued function that describes the response of an imaging system in relation to spatial frequency; it is a crucial parameter in analyzing and evaluating the performance of imaging systems. The OTF of the diffraction limit system is expressed as:

$$\begin{array}{cc}\hfill H\left(\nu \right)=& \frac{1}{\pi \kappa \nu }\cos \left(4\kappa {\nu }^2\right)\{J_1\left(4\kappa \nu \right){\cos }^{-1}\left(\nu \right)-\hfill \\ & {\displaystyle \sum _{n=0}^{\infty }\frac{{\left(-1\right)}^n}{2n}\left[J_{2n-1}\left(4\kappa \nu \right)-J_{2n+1}\left(4\kappa \nu \right)\right]\sin \left[2n{\cos }^{-1}\left(\nu \right)\right]}+\hfill \\ & \frac{1}{2\kappa \nu }{\displaystyle \sum _{n=0}^{\infty }{\left(-1\right)}^{n}2n{J}_{2n}\left(4\kappa \nu \right)\sin \left[2n{\cos }^{-1}\left(\nu \right)\right]}\}+\hfill \\ & \frac{1}{\pi \kappa \nu }\sin \left(4\kappa {\nu }^2\right)\{{\displaystyle \sum _{n=0}^{\infty }\frac{{\left(-1\right)}^n}{2n-1}\left[J_{2n-2}\left(4\kappa \nu \right)-J_{2n}\left(4\kappa \nu \right)\right]\sin \left[\left(2n-1\right){\cos }^{-1}\left(\nu \right)\right]}-\hfill \\ & \frac{1}{2\kappa \nu }{\displaystyle \sum _{n=0}^{\infty }{\left(-1\right)}^n\left(2n-1\right)J_{2n-1}\left(4\kappa \nu \right)\sin \left[\left(2n-1\right){\cos }^{-1}\left(\nu \right)\right]}\},\hfill \end{array}$$

(20)

where $J_n\left(z\right)$ is the Bessel function of the first kind of the order n. Each subfigure in Figure 2 corresponds to an OTF curve in Figure 4. The OTF gradually decreases at each frequency as the degree of defocus increases.

4. The Characteristic Parameters of the Solution

4.1. Equivalent Bandwidth and Equivalent Linewidth

Both the OTF and the PSF are complex functions, for which it is convenient to introduce a single parameter to characterize the image quality for comparison. Goodman proposed equivalent bandwidth and equivalent linewidth as evaluation indexes [17]. We can assume that the main components of the signal (the image and its spatial spectrum) in both the spatial and frequency domains are concentrated only in a finite region near the center. To measure the spreading of the function g(x) and its bandwidth G(ν), the relation between the function g(x) and G(ν) is given by Fourier transform:

$$\{\begin{array}{c}G\left(\nu \right)={\displaystyle {\int }_{-\infty }^{+\infty }g\left(x\right)\exp \left(-2i\pi \nu x\right)dx},\\ g\left(x\right)={\displaystyle {\int }_{-\infty }^{+\infty }G\left(\nu \right)}\exp \left(2i\pi \nu x\right)d\nu .\end{array}$$

(21)

It can be established from the above two equations that:

$$G\left(0\right)={\displaystyle {\int }_{-\infty }^{+\infty }g\left(x\right)dx},\hspace{1em}g\left(0\right)={\displaystyle {\int }_{-\infty }^{+\infty }G\left(\nu \right)}d\nu .$$

(22)

Equation (22) has a clear geometric meaning: if the height of a rectangle is equal to G(0) and the area is the same as the area under the curve G(ν), its width is called “equivalent bandwidth”; see Figure 5a. The reciprocal of the slope of the object image’s light intensity curve at the object’s edge is called the “equivalent line width”; see Figure 5b. The equivalent bandwidth and equivalent linewidth can be expressed as follows:

$$\tilde{\omega }=\frac{{\displaystyle {\int }_{-\infty }^{\infty }G\left(\nu \right)d\nu }}{G\left(0\right)},\hspace{1em}\tilde{s}=\frac{{\displaystyle {\int }_{-\infty }^{\infty }g\left(x\right)dx}}{g\left(0\right)}.$$

(23)

The product of equivalent linewidth and equivalent bandwidth can be obtained from Equations (22) and (23):

$$\tilde{s}\tilde{\omega }\approx 1.$$

(24)

The above equation represents the mutual restriction between the signal spreading in the spatial domain and that in the frequency domain, which is also called “uncertainty”.

Due to the symmetry, both equations $H(\alpha ,\nu )$ and $h(\alpha ,\rho )$ have rotationally symmetrical surfaces. One may make a cylinder where the top surface is at a tangent to the surface and where the bottom surface coincides with the coordinate plane, such that the volume of the cylinder is equal to the volume under the surface (Figure 6). It is physically reasonable to take the diameter 2R of the cylinder as either the equivalent bandwidth of the OTF or the equivalent linewidth of the PSF.

When the surface expression is G(r), then the volume of the cylinder is:

$$\pi R^{2}h={\displaystyle {\int }_0^{2\pi }d\phi }{\displaystyle {\int }_0^{\infty }G\left(r\right)rdr}=2\pi {\displaystyle {\int }_0^{\infty }G\left(r\right)rdr}.$$

(25)

From Equation (25), the following result is obtained:

$$R=\sqrt{2{\displaystyle {\int }_0^{\infty }G\left(r\right)rdr}/h}=\sqrt{{\displaystyle {\int }_{-\infty }^{\infty }G\left(r\right)rdr}/h},$$

(26)

where h is the maximum value of the surface; that is the center value.

When defocus occurs, the equivalent bandwidth will become smaller, while the equivalent linewidth will become larger, indicating increased spreading. By substituting the expression for the OTF (Equation (20)) and the expression for the PSF (Equation (19)) into Equation (26), a series of equivalent bandwidth and equivalent linewidth values are obtained, as depicted in

Table 1. Equivalent linewidth $\tilde{s}$ and equivalent bandwidth $\tilde{\omega }$ values at defocus.

$\mathit{\delta }/\mathit{\lambda }$	$\tilde{\mathit{s}}$ (μm)	$\tilde{\mathit{\omega }}$ (mm−1)	$\tilde{\mathit{s}}\tilde{\mathit{\omega }}$
0	2.27	440	0.998
6	2.39	418	0.999
12	2.80	356	0.997
18	3.74	265	0.991

.

The last row of the above table shows that, corresponding to a given amount of defocus, the product of equivalent linewidth and equivalent bandwidth is approximately equal to 1.

4.2. Strehl Ratio and Resolution

The Strehl ratio (SR) is a measure of the relative peak intensity of an optical system’s PSF, compared to that of an ideal, aberration-free system. In Equation (19), the SR can be calculated by setting $n=0$, which requires considering only the term $\rho =0$ in the following equation:

$$SR={h\left(A,\delta ;\rho \right)|}_{\rho =0}={\left[{\displaystyle \sum _{k=0}^{\infty }\frac{{\left(-1\right)}^k{\kappa }^{2k}}{\left(2k+1\right)!}}\right]}^2+{\left[{\displaystyle \sum _{k=1}^{\infty }\frac{{\left(-1\right)}^k{\kappa }^{2k-1}}{\left(2k\right)!}}\right]}^2.$$

(27)

This formula can be simplified to a trigonometric function:

$$SR\left(A,\delta \right)={\left[\frac{2\sin \left(\kappa /2\right)}{\kappa }\right]}^2.$$

(28)

In Figure 7, the red curve represents the confluent hypergeometric function, and the blue is calculated by Equation (28); their values are highly consistent. The values of each point are shown in

Table 2. Defocus amount versus the SR.

$\mathit{\delta }/\mathit{\lambda }$	0	3	6	9	12	15	18	21	24	27	30
SR	1	0.975	0.902	0.791	0.654	0.506	0.362	0.235	0.133	0.061	0.018

.

The relationship between the SR and resolution will be discussed in the following content. Resolution typically denotes the maximum spatial frequency that can be distinguished. Within industry settings, the resolution of the design output when the OTF falls between 0.15 and 0.2 is commonly employed as a reference value.

Figure 8 shows a series of defocus OTF curves, obtained with an optical system numerical aperture of 0.242 and a main wavelength of 550 nm, the defocus range being 0~21 times the wavelength. The intersection points between these curves and the line where OTF = 0.2 are regarded as the resolution in each defocus state.

Figure 9a shows the curve of spatial resolution as a function of the SR. The abscissa of each discrete point represents the intersection point between the curve of different defocus states (see Figure 8) and the line OTF = 0.2, which is the resolution. The ordinate is the SR value at the corresponding defocus position. The ordinate in Figure 9b is the logarithm of the SR; the results show that the resolution is more closely linearly related to ln(SR).

5. Comparison of the Analytical Solution with the Ray-Tracing Results of the Real System

To verify the accuracy of our analytical results, we implemented ray tracing with the optical design software, Zemax (version 18.4.1), to generate both the FFT PSF and FFT MTF. The FFT used here refers to the fast Fourier transform technique, which is the method employed by Zemax software for calculating the PSF and MTF of optical systems. The MTF is the absolute value of the OTF, known as the modulus of the OTF. The calculation method for the FFT PSF is as follows. A set of grid rays emitted from the source point propagate to the exit pupil; the complex amplitude of the corresponding point on the wavefront grid at the exit pupil is calculated using the amplitude and path difference. After appropriate scaling and Fourier transformation, the square of the resultant values yields the real value of the PSF. The FFT MTF essentially represents the Fourier transform of the FFT PSF. In subsequent content, these will be referred to simply as PSF and MTF.

The 2D system layout of the three systems used for comparison is illustrated in Figure 10. Light rays from the left side of the system converge on the image plane on the right side.

Table 3. Main parameters of the optical systems.

System	Wavelength	F/#	NA′	f′	Diameter (Entrance Pupil)	RMS/Waves
1	588 nm	2.00	0.2425	20	10.0	0.0039
2	550 nm	7.14	0.0698	500	70.0	0.0124
3	550 nm	2.54	0.1961	100	40.0	0.0005

lists the main optical parameters of systems 1–3. For systems 1–3, which have respective A values of 0.2425, 0.0698, and 0.1961, step values of δ/λ were chosen as 5, 50, and 5 to accurately reflect the variations in MTF and PSF with defocus.

We compared the results obtained by the Zemax software’s ray tracing with the PSF and MTF calculated by the analytical solution. The defocused MTF and PSF curves of the three systems are shown in Figure 11, Figure 12 and Figure 13, respectively, where the solid line represents the result calculated by the analytical solution and the dashed line represents the result obtained by ray tracing.

Table 4. Comparison of the PSF and the MTF for different systems.

	System 1 (RMSE)		System 2 (RMSE)		System 3 (RMSE)
MTF	δ = 0	0.25%	δ = 0	0.33%	δ = 0	0.08%
	δ = 5λ	0.39%	δ = 50λ	0.23%	δ = 5λ	0.22%
	δ = 10λ	0.92%	δ = 100λ	0.27%	δ = 10λ	0.60%
	δ = 15λ	1.48%	δ = 150λ	0.63%	δ = 15λ	1.10%
	δ = 20λ	1.72%	δ = 200λ	0.94%	δ = 20λ	1.57%
	δ = 25λ	1.54%	δ = 250λ	1.05%	δ = 25λ	1.89%
	δ = 30λ	1.34%	δ = 300λ	0.97%	δ = 30λ	2.01%
PSF	δ = 0	1.22%	δ = 0	1.23%	δ = 0	1.23%
	δ = 5λ	1.04%	δ = 50λ	1.22%	δ = 5λ	1.16%
	δ = 10λ	0.76%	δ = 100λ	1.24%	δ = 10λ	1.00%
	δ = 15λ	0.95%	δ = 150λ	1.31%	δ = 15λ	0.79%
	δ = 20λ	1.17%	δ = 200λ	1.35%	δ = 20λ	0.60%
	δ = 25λ	0.99%	δ = 250λ	1.25%	δ = 25λ	0.50%
	δ = 30λ	0.56%	δ = 300λ	0.97%	δ = 30λ	0.48%

presents the root mean square error (RMSE) values calculated between the MTF and PSF, obtained through two different methods at different defocus levels.

6. Discussion

Many imaging systems are precision-designed and fabricated to achieve imaging quality that is close to the diffraction limit, in order to deliver exceptional performance. These systems include cameras equipped in aerospace vehicles, high-end plan-apochromatic microscope objectives used in medicine and life sciences, object identification systems, and imaging systems employed in semiconductor applications. For high-end imaging systems, the most commonly encountered imaging error is defocusing rather than inherent aberrations. Defocusing may be caused by variations in object distance, fabrication and assembly tolerances, environmental temperature fluctuations, and other factors.

By applying the above confluent hypergeometric functions, this paper offers an analytic solution for the Fresnel diffraction of a diffraction-limited system. The universal model describes the behaviors of systems under defocusing. The analytic model’s most compelling feature is its complete independence from intricate design details, objective types, and various application scenarios. It is anticipated that such an analytic formula will facilitate the analysis of image quality deterioration, the establishment of suitable tolerances for design and fabrication, and the prediction of defocusing effects.

7. Conclusions

In this paper, the analytic formula of the defocused PSF of the diffraction-limited system utilizing confluent hypergeometric function is first derived, and then the characteristic parameters, including equivalent bandwidth, equivalent linewidth, and the Strehl ratio, are derived. The analytical expression of the Strehl ratio of defocusing systems can be reduced to a simple trigonometric function. Examples of three different types of diffraction-limited systems are analyzed and compared with the analytic solution. The PSF of each example is very close to the analytic solution, and the S values of defocus are consistent with those calculated from the above analytical expressions.

It follows that the defocused PSF of the diffraction-limited system, which includes the confluent hypergeometric function, is indeed a typical expression of a diffraction-limited system. This analytical formula can give the defocus characteristics and behavior of all diffraction-limited systems with known or unknown structural parameters. This analytical solution describes the properties of all diffraction-limit systems, including defocus performance, resolution, etc.