Analytical Fresnel imaging models for photon sieves.

Photon sieves are a fairly new class of diffractive lenses that open unprecedented possibilities for high resolution imaging and spectroscopy, especially at short wavelengths such as UV and x-rays. In this paper, we model and analyze the image formation process of photon sieves using Fourier optics. We derive closed-form Fresnel imaging models that relate an input object to the image formed by a photon sieve system, both for coherent and incoherent illumination. These analytical models also provide a closed-form expression for the point-spread function of the system for both in-focus and out-of-focus cases. All the formulas are expressed in terms of Fourier transforms and convolutions, which enable easy interpretation as well as fast computation. The derived analytical models provide a unified framework to effectively develop new imaging modalities enabled by diffractive lenses and analyze their imaging capabilities for different design configurations, prior to physical production. To illustrate their utility and versatility, the derived formulas are applied to several important special cases such as photon sieves with circular holes and pixelated diffractive lenses generated by SLM-type devices. The analytical image formation models presented in this paper provide a generalizable and powerful means for effective analysis and simulation of any imaging system with a diffractive lens, including Fresnel zone plates, Fresnel phase plates, and other modified Fresnel lenses and mask-like patterns such as coded apertures.


Introduction
Photon sieves are modified Fresnel zone plates in which a large number of holes are distributed over the transparent Fresnel zones. Similar to Fresnel zone plates, photon sieves offer diffractionlimited resolution, but with several additional advantages [1]. These include tighter focus for a given smallest structure size (by enlarging the outer diameter) and suppression of secondary maxima and higher diffraction orders (by varying the density and the diameter of the holes) [1,2]. Furthermore, photon sieves consist of a single connected piece, and are robust to manufacturing errors, therefore enable low-cost and simpler fabrication [3][4][5][6].
Photon sieves are particularly useful at short wavelengths such as UV and x-rays. At these wavelengths, refractive lenses do not exist due to the strong absorption of available materials, and reflective mirrors with near diffraction-limited resolution are costly to manufacture. In fact, figure errors and surface roughness cause the resolution of reflective optics to be considerably worse than the diffraction limit [3,7,8]. As a result, lightweight photon sieves open new possibilities for high resolution microscopy, spectroscopy, lithography, and ultra-large telescopes [1,3,[9][10][11][12][13][14].
As the advent of photon sieve imaging systems with new configurations expands, development of analytical and computationally efficient imaging models becomes crucial for the design and analysis of such reconfigurable systems, as well as for digital image reconstruction. In this paper, we present analytical imaging models for photon sieves using Fourier optics [39]. To accomplish this, we consider the imaging system shown in Fig. 1 where a photon sieve resides between an extended object and an image plane. We derive closed-form Fresnel imaging models that provide the relationship between the input object and the image formed at the measurement plane, both with coherent and incoherent illumination. The derivations also provide a closedform expression for the point-spread function (PSF) of the overall photon-sieve system for both in-focus and out-of-focus cases.
The presented imaging formulas and point-spread functions are all expressed in terms of Fourier transforms and convolutions (instead of complicated integrals), which provide an easier understanding and insight into the imaging process. Moreover, this enables the fast and accurate simulation of the overall imaging system using fast Fourier transform (FFT) [40], which removes the need to use computationally intensive softwares, numerical integration, or diffraction computations from multiple planes [41,42]. These analytical models and point-spread functions are crucial for effective development, design, and analysis of the new imaging modalities enabled by photon sieves, prior to their physical production. In particular, the implications of various design choices, such as different hole shapes [43], pixelization [28], and apodization [30], on imaging attributes can be readily analyzed and simulated. To illustrate their utility and versatility, the derived formulas are applied to several important special cases such as photon sieves with circular holes and pixelated diffractive lenses generated by SLM-type devices. These imaging models can also be used as an accurate model of the imaging process for image deconvolution and reconstruction tasks involved in computational imaging systems with diffractive lenses [37,38].
The focusing properties of photon sieves have been earlier analyzed through Fresnel-Kirchhoff diffraction integrals [1], Fresnel integrals [44,45] and Rayleigh-Sommerfeld diffraction [46]. While the calculations in [1] were limited to point sources only, the approximate treatment of Fresnel diffraction [44], numerical Fresnel propagation [45], and Rayleigh-Sommerfeld diffraction [46] were computed for the half of the imaging system considered here (i.e. only from the photon sieve plane to the image plane), hence requiring the knowledge of the complex-valued field at the photon sieve plane. A similar approach has also been developed to compute the diffracted fields of rings instead of individual pinholes for faster simulation [47].
In this paper, by considering the overall photon-sieve system in Fig. 1, we derive the complete analytical image formation model applicable to extended input objects. Our models are based on Fresnel (near-field) diffraction, which is a part of scalar diffraction theory, and hence do not require far-field (or Fraunhofer) conditions to hold. Moreover, while the approximate Fresnel treatment in [44] are only for photon sieves with circular holes, our framework is generalizable in that it can be applied to imaging systems with any photon sieve configuration. In fact, the imaging models presented in this paper are powerful in that they can be used to effectively analyze and simulate any imaging system with a diffractive lens, including modified Fresnel lenses and any other mask-like patterns such as coded apertures.
The paper is organized as follows. In Section 2, analytical Fresnel imaging models and PSFs are presented for both coherent and incoherent illumination. Preliminary versions of some of these models were presented in [48]. In Section 3, the derived formulas are applied to several important special cases with commonly used design configurations. To illustrate their use, in Section 4, numerical simulations are presented for the analysis of the imaging properties of the photon sieve system under various design and observing scenarios.

Imaging setting
As shown in Fig. 1, a photon sieve resides between an extended object and an image plane. The distances from the photon sieve plane to the object and image planes are denoted by d s and d i , respectively. Suppose the light originated from the source (object) plane is a spatially coherent, space-varying, and monochromatic wave with wavelength λ. (Spatially incoherent case will be treated later in Section 2.3.) Then the complex baseband representation of the input field can be expressed as a linear combination of plane waves propagating in a continuum of directions [49], given by where A(α, β) represents the angular spectrum of the input field u(x, y), and is related to its with the two-dimensional (2D) Fourier transform defined as Any photon sieve can be fully characterized by its amplitude transmittance (aperture) function. The transmittance function, t(x, y), is the ratio of the transmitted field amplitude to the incident field amplitude at each point (x, y) in the photon sieve plane. The aperture t(x, y) of a photon sieve (or in general a diffractive lens) is often binary (taking values 1 or 0). However, here we allow t(x, y) to be an arbitrary complex-valued function to keep the derived models general enough and applicable to any diffractive imaging element or lens, including Fresnel phase plates, or any other mask-like pattern that has both amplitude and phase modulation.
For instance, for the conventional photon sieve configuration with circular-shaped holes, the transmittance function is given by where t n (x, y) is the transmittance of the nth pinhole with diameter d n and central location (x n , y n ), and N denotes the total number of holes. This transmittance t n (x, y) takes value 1 inside the nth circular pinhole, which is above expressed in terms of the circle function [49] circ(x, y) = Note that different design configurations such as different shapes of holes, apodization, and pixelization can be easily included to the aperture function by using a proper mathematical form for t(x, y). This will be illustrated in Section 3. A photon sieve can form images either with coherent or incoherent illumination, which will be successively discussed below.

Fresnel imaging formula for coherent sources
Theorem 1 Within the Fresnel approximation, the image r(x, y) formed at a distance d i after a photon sieve is related to the complex input wavefront u(x, y) located at a distance d s behind the photon sieve (as shown in Fig. 1) through where * denotes 2D convolution,ũ(x, y) is an inverted, scaled (magnified) and chirp modulated version of the input wavefront u(x, y) with a position-dependent phase shift, given bỹ and h d i ,λ (x, y) represents the coherent point-spread function (PSF) of the photon sieve at distance d i and wavelength λ: Here which is simply a magnified and chirp modulated version of the aperture function. As a result, the formed image can be interpreted as a scaled and inverted version of the input signal u(x, y) with a position-dependent phase shift, filtered in the frequency domain by the chirp modulated, scaled aperture function H d i ,λ ( f x , f y ). Naturally, the magnification of the image relative to the object is d i /d s , and the image is upside down relative to the object's orientation. Also the intensity i(x, y) of the image is given by Second, when d s ≫ d i , the coherent PSF of the photon sieve reduces to with the corresponding coherent transfer function being that is, they only change with λd i . Hence, the PSF is the same for different values of λ and d i as long as λd i remains the same. Third, the PSF is circularly symmetric whenever the aperture function is circularly symmetric. This is because (i) a circularly symmetric aperture function makes the transfer function circularly symmetric, and (ii) this in turn makes the PSF circularly symmetric since the Fourier transform of a circularly symmetric function is also circularly symmetric [49].
Fourth, the presented imaging formula is general enough to accommodate for any type of aperture function for a photon sieve or more generally for a diffractive imaging element. Moreover, the imaging relation expressed in terms of Fourier transforms and convolutions enables an easier understanding of the imaging process, as well as fast and accurate simulation of the imaging system.
In particular, the given imaging formula suggests two different ways for the numerical computation of the coherent PSF. A direct computation of Eq. (8) requires first computing the Fourier transform of the aperture function and then convolving the result with a complex exponential (i.e. chirp). Alternatively, based on Eq. (9), one can first compute the coherent transfer function by multiplying the scaled aperture function with a complex exponential, and then compute the inverse Fourier transform of the transfer function to obtain the PSF. Because the convolution in the direct approach is a time-consuming operation, the approach based on Eq. (9) is the clear choice for faster computation.
Hence, on a computer, the samples of the PSF can be computed by first sampling the transfer function in Eq. (9) and then computing its inverse FFT. That is, the PSF of the overall imaging system can be computed through a single FFT computation, without resorting to computationally intensive softwares, numerical integration, or diffraction computations from multiple planes [41,42]. This fast computation method works accurately provided that the separation between the samples of the transfer function is sufficiently small (to satisfy the Nyquist sampling criterion) and the total number of samples is sufficiently large (to contain the significant amount of total energy of the PSF in the computation) [40,50].
Proof. The field v(x, y) just before the photon sieve is related to the input field u(x, y) by Here the input field is convolved with the PSF of the free space propagation under Fresnel approximation [39,49] (where the constant phase term is dropped). If the open form of this convolution is used together with the 2D Fourier transform expression in Eq. (3), an equivalent form can be obtained as follows [39,Chap.4]: Hence the relation between v(x, y) and u(x, y) can be expressed as a chirp (quadratic phase exponential) multiplication, followed by the evaluation of the Fourier transform at the frequency { f x = x/(λd s ), f y = y/(λd s )}, followed by a second chirp multiplication. Secondly, the field amplitude s(x, y) just after the aperture of the photon sieve is related to v(x, y) as where t(x, y) is the aperture function.
Lastly, the complex amplitude, r(x, y), in the image plane is related to s(x, y) through another free space propagation: By using the above three relations together with the properties of the Fourier transform, the image amplitude r(x, y) can be obtained as follows: where T( f x , f y ),ũ(x, y), and ∆ are as defined before.
This result, derived from first principles, could also be obtained by invoking the lens law. The lens law [49] states that for a lens with transmittance t(x, y) = e −iπ x 2 +y 2 λ f l q(x, y) where q(x, y) describes the finite aperture and f l is the focal length, the image formed at the distance d i behind the lens is, within the Fresnel approximation, given by where in the above expression to arrive at the same result: Note that the presented models are valid for diffractive elements operating at the Fresnel (near-field) regime of scalar diffraction theory. As well-known [39], the scalar theory (such as Rayleigh-Sommerfeld diffraction) is accurate when both (i) the dimensions of the smallest diffracting structure, and (ii) the distance from the aperture to the observation plane are large compared to a wavelength. Fresnel approximation is an additional approximation on top of that made for scalar theory, and there are also sufficient conditions for its validity [39,49]. These conditions check whether the distance from the diffractive lens to the image plane, d i , as well as to the source plane, d s , are large enough. On the other hand, these sufficient conditions are known to be overly stringent, and not necessary for the accuracy of the Fresnel approximation. Based on different analyses [39,51,52], it has been observed that the accuracy of the widely used Fresnel approximation is extremely good to distances that are very close to the diffracting aperture (or equivalently, to field angles that are too large). Hence, we can conclude that the models in this paper are accurate provided that the hole diameters are large compared to a wavelength, and the image and source planes are not very close to the diffractive lens.
These two requirements for the validity of the models limit the maximum numerical aperture (NA) of the diffractive lens and hence its resolution. That is, there exist practical cases for which the required conditions are not satisfied. As an example, consider a binary Fresnel zone plate with outer diameter D, first-order focus f , and outer zone width w. In this case, f = Dw/λ is satisfied [7], and NA can be approximated as D/2 f = λ/2w. Since the validity of scalar diffraction theory requires w ≫ λ, it is required that f ≫ D and N A ≪ 0.5 for this type of Fresnel lens. In general, when the size of the smallest diffracting structure approaches the wavelength and polarization effects become significant, the scalar diffraction theory and hence the models in this paper are inadequate. In these cases, the vector diffraction theory or the finite-difference time-domain based numerical computations are needed [53][54][55][56][57]. Moreover, for cases where Fresnel approximation fails only, a more richer scalar diffraction theory such as Rayleigh-Sommerfeld diffraction can be used [46]. Nevertheless, even when the presented analytical models are not valid, they will still be useful to understand the trends in the imaging behavior that may be difficult to achieve with numerical computations.

Theorem 2 For incoherent imaging within the Fresnel approximation, the average intensity
whereũ i (x, y) is the inverted and scaled object intensity given bỹ and g d i ,λ (x, y) denotes the incoherent PSF of the photon sieve at distance d i and wavelength λ: The angular bracket ⟨.⟩ here denotes a statistical (ensemble) average (with the assumption on ergodicity, which can be interpreted as time average). Hence for incoherent illumination, the image intensity is given by the convolution of the ideal image intensityũ i with the impulse response g d i ,λ . That is, the intensities are convolved, rather than the complex amplitudes as in the coherent case. Naturally, the incoherent PSF of the photon sieve is just the squared magnitude of its coherent PSF given in Eq. (8): Note that the intensity recorded in a stationary system is actually the sample time average given by i(x, y) = 1 T ∫ T 0 |r(x, y, t)| 2 dt, whereas here we are dealing with the ensemble average i(x, y) = ⟨|r(x, y, t)| 2 ⟩. When the ergodicity assumption holds, this ensemble (statistical) average is equivalent to the observed time average at the detector.
Proof. A spatially incoherent wave emanating from the object u(x, y) can be described by a random complex amplitude that varies with spatial position and time [49] as follows: Here ρ(x, y, t) represents a spatially white random process satisfying ⟨ρ(x, where δ(.) denotes the Dirac delta function. The spatial autocorrelation function of the input wave w(x, y, t) can then be written in terms of the object intensity u i (x, y) as that is, there is no correlation between the values at different spatial positions. Using this model for spatial incoherence and the derived imaging model for the coherent case (Eq. (6)), an expression can be obtained for the image intensity as given by the statistical average of the squared image magnitude: where Lastly, we note that the case of imaging partially coherent scenes could be treated similarly if the two-dimensional spatial correlation function is known for the input. In this case, the derivation above should be repeated by replacing the Dirac delta function in Eq. (25) with the known spatial correlation function.

Some examples of diffractive lenses
Here we illustrate the generality of the presented formulas by applying them to some important special cases with commonly used design configurations. These also illustrate that analytical PSFs can be obtained when analytical expressions exist for the transmittance function.

Conventional photon sieve with circular holes
The conventional photon sieve (PS) configuration consists of circle-shaped holes, so its transmittance can be expressed as sum of circle functions as in Eq. (4). By inserting this form for t(x, y) in Eqs. (8) and (9), the coherent PSF and transfer function of this type of photon sieve can be respectively obtained as follows: Here ∆ = 1/d i + 1/d s , the total number of holes is N, and the nth pinhole's diameter is denoted by d n and central location by (x n , y n ). The two-dimensional jinc function appears in the PSF expression because it is the Fourier transform of the circle function, and defined as [49] jinc(x, y) = jinc where J 1 (.) is the first-order Bessel function of the first kind. For this particular photon sieve with circular holes, the Fresnel PSF has been given in [44] in an integral form for the case d s ≫ d i . It can be shown that the PSF given here in the convolution form is the general form of this. For this, the convolution in Eq. (27) should be written in open form; after some calculations and assuming d s ≫ d i , it can be simplified as Here J 0 (.) is the 0th order Bessel function of the first kind and R 2 = (x − x n ) 2 + (y − y n ) 2 . Other than the complex exponential (chirp) term in front, this expression is same as the one in [44] when the illumination there is taken as a plane wave of normal incidence.

Pixelated structures
Another important special case is a pixelated diffractive lens, which is often realized by a spatial light modulator (SLM) or a digital micromirror device (DMD) [28]. Because a pixelated diffractive lens also acts as a diffraction grating, the PSF will have contributions from many diffraction orders. Analytical expressions can be obtained for the PSF and the transfer function after properly formulating the pixelated aperture function. Without loss of generality, we consider rectangular shaped pixels. Let t(x, y) be the continuous transmittance function and t p (x, y) be its pixelated version with the pixel pitch in x and y directions given by d x and d y , and the pixel dimensions given by a x and a y , respectively. Then the mathematical relation between t p (x, y) and t(x, y) can be expressed as where the 2D comb function is the infinite train of impulses given by and the 2D rectangle function is defined as The Fourier transform of this aperture function can be obtained as where the 2D sinc function is the Fourier transform of the rectangle function and defined as sinc(x, y) = sinc(x)sinc(y), with sinc(t) = sin(πt)/πt. By inserting this pixelated aperture function t p (x, y) and its Fourier transform in Eqs. (8) and (9), the coherent PSF and transfer function can be respectively obtained as follows: (37) Here ∆ = 1/d i + 1/d s as before. As expected, the PSF is composed of infinitely many terms, each associated with a different diffraction order. The (m, n) diffraction order is centered at (mλd i /d x , nλd i /d y ), and each order contains modulation with a 2D sinc envelope. The PSF due to (0, 0) diffraction order mainly differs from the non-pixelated case (given in Eq. (8)) by this sinc envelope, and the separation between different diffraction orders increases as the pitch size decreases, as expected.

Binary Fresnel zone plate
Photon sieves are modified Fresnel zone plates in which a large number of holes are distributed over the transparent Fresnel zones. Because of this connection, many studies about the photon sieve and its design criteria are closely related to that of a Fresnel zone plate [2,3,28,58,59]. Here we consider a binary Fresnel zone plate and provide closed-form expressions for its PSF and transfer function. To the best of our knowledge, these expressions have not appeared elsewhere before. Since the binary Fresnel zone plate can be viewed as an infinite series of thin lenses (as discussed below), its PSF and transfer function have simpler analytical forms, and, as will be illustrated in Section 4, these simpler formulas closely approximate the PSF/transfer function of a conventional photon sieve, as one would expect. Hence an equivalent binary Fresnel zone plate can also be used to approximately simulate conventional photon sieve systems.
To obtain a closed-form expression for the PSF/transfer function of a binary Fresnel zone plate, its transmittance function needs to be properly formulated. As an example, we consider a binary Fresnel zone plate (FZP) with transparent odd zones, and whose first order focus is f and outer diameter is D. The aperture function of this FZP can be mathematically expressed as follows [39,49]: where the sign function sgn(x) takes value 1 if x is positive and −1 otherwise. By using the Fourier series expansion of a periodic square wave [60], this aperture function can be rewritten as where f m = f /m. Note that each odd term in the above summation corresponds to the transmittance function of a thin lens with focal length f m . That is, the FZP can be viewed as an infinite series of thin lenses, each associated with a different diffraction order m. Hence, FZP generates infinitely many diffraction orders, and each odd order m comes to focus at a focal distance f m [7].
As a result, the PSF of the FZP should consist of infinite terms, each due to a different diffraction order and representing whether the mth diffraction order comes to focus or not at the distance d i . The analytical expression for the coherent PSF can be obtained by either using the known PSF formulas for a thin lens with circular aperture [49,Chap.4] or by evaluating Eq. (8) for the aperture function in Eq. (39). This results in the following PSF: where each contributing term to the PSF due to a different diffraction order is Here ∆ = 1/d i + 1/d s as before, and ϵ m = 1/d i + 1/d s − 1/ f m indicates whether the mth lens (i.e. mth order diffraction) is in focus (when ϵ m = 0) or out of focus (when ϵ m 0) at the distance d i . Hence, as expected, the PSF of the binary FZP is composed of infinitely many terms, each associated with a different diffraction order and depending on whether the corresponding order is in focus. In fact, h m (x, y) is the coherent PSF of the mth contributing lens in the aperture function, and its form depends on whether the mth diffraction order is in focus or not. In the focused case ( ϵ m = 0), the intensity |h m (x, y)| 2 is simply the well-known airy pattern [39,49]. For example, if a monochromatic collimated beam is focused at the first order focus f (i.e. d i = f , d s = ∞, ϵ 1 = 0), the airy pattern is observed when all the contributing terms in Eq. (40) other than the first order diffracted field h 1 (x, y) have negligible energy. Similarly, the coherent transfer function is with The diffraction efficiencies for different orders also follow as given in [7]: Hence 50% of the incident energy is absorbed by the zone plate through the non-transparent zones (associated with the vanishing even orders). 25% is passed without being focused (m = 0), 10% appears on the first order focus f 1 , (another 10% goes to the divergent first order, i.e. m = −1), 1% appears on the third order focus f 3 , etc. If this lens model will be used as an approximate model for the conventional photon sieve, it should be taken into account that white zones of the Fresnel zone plate are replaced with holes. This in turn affects the amount of transmitted light through the diffractive imaging element, and hence the diffraction efficiencies given above. If γ is the fraction of the area of the open holes (in the photon sieve) to the area of the white zones (in the Fresnel zone plate), then the diffraction efficiency should approximately scale by γ 2 . However, note that if the levels of the side-lobes are different for the PSFs of the PS and the equivalent FZP, then the peak intensity can also differ from this predicted value.

Binary diffractive lenses with reversed transparency
When working with binary diffractive lenses, it may be of interest to exchange the transparent and opaque regions. A well-known example is an FZP with transparent even zones, which can be obtained from an FZP with transparent odd zones, described in Eq. (38), by switching the transparent and opaque regions. If two diffractive lenses are related to each other in this manner, the PSF/transfer function of one of them can be directly obtained from the PSF/transfer function of the other.
To explicitly show this, let us consider two diffractive lenses with aperture functions t 1 (x, y) and t 2 (x, y) such that their sum is a circle function: where D denotes the outer diameter of the diffractive lenses. Then their PSFs must sum to the PSF of a circular aperture. One can prove this version of Babinet's principle [49] as follows.
The PSF of a circular aperture can be obtained by inserting the Fourier transform of its aperture function in Eq. (8): Let T 1 ( f x , f y ) and T 2 ( f x , f y ) denote the Fourier transforms of the aperture functions t 1 (x, y) and t 2 (x, y), and h 1 (x, y) and h 2 (x, y) denote their respective PSFs. Using Eq. (45) in replace of the circular aperture function, one can also obtain the following: Hence, as shown, the PSFs sum to the PSF of the circular aperture given in Eq. (46). Similarly, the transfer functions of such diffractive lenses sum to the transfer function of a circular aperture; that is, (48) where the first equality is obtained using Eq. (9).
To illustrate the usefulness of these relationships, we consider an FZP with transparent even zones, which is obtained from an FZP with transparent odd zones described in Eq. (38) by switching the transparent zones with opaque ones. Each term contributing the coherent PSF and transfer function of this type of FZP can be obtained from the PSF and transfer function terms given in Eqs. (41) and (43), respectively, by using the relations in Eqs. (47) and (48): Note that when the above formulas are compared with the ones in Eqs. (41) and (43), it is clear that, due to the contribution from the 0th order term, the overall PSF (OTF) of an FZP does not remain exactly the same (although similar) after reversing the transparent and opaque regions.

Numerical results
For demonstration of the usefulness and generality of the presented imaging models, numerical simulations are performed to illustrate the imaging properties under various different design configurations and observing scenarios. For this, we consider and compare four different designs of diffractive lenses. The first one is a binary FZP design with transparent even zones. Based on this, two sample PS designs, one with circular holes and another with rectangular holes, are generated. A pixelated version of the PS design with circular holes is also included in the comparison. For all designs, the wavelength is chosen as 630nm, the outer diameter as 10mm, and the width of the outer zone as 50µm. The resulting FZP has a focal length of 0.7937m and 25 transparent zones. For the PS designs, the fraction of the transparent area due to holes (relative to the zone area) is chosen as 0.6, and the hole size is selected to be 1.53 times the underlying zone width. This resulted in 2662 holes, distributed uniformly within each zone (i.e. separated by the same angle) and with the smallest hole size being 76.5µm. For the PS design with square holes, the width of the square holes is taken to be equal to the diameter of the circular holes in the first PS design. In the pixelated version of the first PS design (with circular holes), the pixel pitch, d x , and the effective pixel dimension, a x , are respectively chosen as 50µm and 40µm, and same in the y direction, resulting in a fill factor of 64%. Fig. 2 shows the resulting aperture functions for all designs.
For comparison, two different observing scenarios are considered. In both cases, the photon sieve is located at a distance d s = 15m from the source. In the first scenario, the image plane at a distance d i = 0.838m after the photon sieve is considered, which corresponds to the plane of (first-order) focus. In the second scenario, the plane that is three depth of focus (DOF) away from the focal plane is considered, resulting in a defocus case. More specifically, DOF [3] is equal to 0.7937mm for the above designs, and hence in the second scenario the distance of the image plane to the photon sieve is chosen as d i = 0.8618m.
For all designs and scenarios, first the samples of the coherent transfer function is computed using the derived formula in Eq. (9) and the respective aperture function, then the coherent PSF is obtained from the sampled transfer function using a single FFT. As mentioned earlier, for accurate computation of the PSF, the separation between the samples of the transfer function should be sufficiently small (to satisfy the Nyquist sampling criterion) and the total number of samples should be sufficiently large (to contain the significant amount of total energy of the PSF in the computation) [40]. That is, if The reason for the decrease of δ f in the pixelated case is the increase in the extent of the PSF due to the additional diffraction orders arising from pixelization. This dictates a higher Nyquist rate to sample the transfer function. Once the coherent PSFs are computed as explained, the incoherent PSFs are also obtained from their squared magnitude (see Eq. (23)). The magnitude of the normalized coherent PSFs for the first (focus) and second (defocus) scenarios are shown in Figs. 3 and 4, respectively. The incoherent PSFs for the first (focus) scenario is also shown in Fig. 5. As seen in these figures, the normalized PSFs of different designs appear to be very similar in their first few lobes. The major difference is observed in the peak intensity values due to the different amount of light transmission. In fact, the peak intensity for the circular PS design is lower than that of the FZP by a factor of 6, and the PS design with square holes is slightly lower than this. For the pixelated one, due to the additional diffraction orders arising from pixelization, the peak intensity is lower than that of the circular PS by a factor of 3. To better see the relative efficiency of the imaging, insets of line-scan across the PSF center are also shown in Figs. 3, 4 and 5. These PSF slices in the insets are normalized to the peak PSF value of the zone-plate. Although not shown here, in the computed PSFs of the pixelated PS, the additional diffraction orders are also observed at the multiples of λd i /d x = 0.0109m, as expected.
To see the differences in the normalized PSFs more clearly, the central slices of the PSFs are also shown in Fig. 6  case, the distance from the peak to the first zero crossing (i.e. Rayleigh resolution) is nearly 1.22∆ = 61µm for all designs, as expected, since all designs have the same outer diameter (i.e. same numerical aperture) [2,3,46,47]. On the other hand, in the three DOF away case, the location of the first zero crossing increases to 450µm. As mentioned before, the PSFs appear to be almost identical in the first few lobes for all designs. The only exception is that the PS with rectangular holes has slightly higher secondary side-lobes than others in the focused case, which is also expected. However, although the PSFs are very similar in the first few lobes, higher side-lobes are visible at the tails for all FZP-like designs compared to the traditional FZP. More specifically, the side-lobes of pixelated PS are the highest, which is followed by rectangular and circular PS.
To prevent possible confusion, we emphasize that the same outer diameter and zone plate geometry are used for all designs in order to perform a fair comparison. As a result, the PS designs considered here do not show the well-known advantages over the zone plates such as sharper focusing and suppression of sidelobes. As mentioned in many earlier works [2,46,47], simply replacing the open rings of a zone plate with the pinholes does not enable these advantages. Our results clearly illustrate this. Nevertheless, using pinholes, one can extend the diameter of a photon sieve beyond the diameter of a zone plate with the same smallest structure, and sharper focusing can be achieved due to this larger aperture size (i.e. increased numerical aperture) [1,2,46,47]. Moreover, using pinholes also allow to smoothly vary the density of pinholes over the zones, and suppression in sidelobes can be achieved when such a smoothly varying aperture (transmission) function is used [1,2,46]. As illustrated in Fig. 2, the PS designs considered here neither employ a larger aperture size nor a smoothly varying aperture function; as a result, the considered PS designs provide neither sharper focusing nor suppression of sidelobes compared to the zone plate.
As stated before, the Fresnel models used in the simulations are accurate provided that the pinhole diameters are large compared to a wavelength, and the image and source planes are not very close to the diffractive lens. For the latter requirement, there are well-known sufficient conditions [39,49,52] that can be checked. For example, the sufficient condition given in [49] for the validity of the Fresnel approximation is as follows: Here, the diffractive aperture is assumed to lie in the (ξ, η) plane, the source/image is in the (x, y) plane, and z is the distance between these two planes. In our results, the width of the observation region is 4mm, the diameter of the diffractive element is 10mm, and the wavelength is 630nm. Inserting these values into the above equation gives the condition that z > 0.16m. Hence the condition is satisfied for both distances, d i and d s , in the focused and defocused cases. Since this is a sufficient condition, the Fresnel approximation and our models derived based on this approximation yield accurate results here. However, it should also be noted that this condition is usually overly stringent, and not necessary for the accuracy of the Fresnel approximation [39,51,52]. Therefore, the results may as well be satisfactory for much shorter distances. Lastly, as discussed in Section 3.3.3, the PSF of the FZP is composed of infinitely many terms, each associated with a different diffraction order. As a final remark, we note that the single PSF term in Eq. (41) corresponding to the 1st diffraction order (i.e. m = 1) provides a good approximation not only to the overall PSF of the FZP at focus but also to the PSF of other FZP-like diffractive lenses. Hence the single lens model (with m = 1) provides a good approximation for these diffractive lenses, while it requires less computations for the simulation of the system. Hence, this model provides an approximate, but a simpler analysis.

Conclusion and discussions
In this paper, we presented closed-form Fresnel imaging formulas that relate the input field distribution to its image formed by a diffractive element such as a photon sieve. Both spatially coherent and incoherent cases are considered, and analytical expressions for the associated PSFs are also provided. All the formulas are given in terms of Fourier transforms and convolutions, which enable easy interpretation as well as fast computation. These formulas are derived based on the scalar diffraction theory and Fresnel approximation, and hence they are accurate provided that the hole diameters are large compared to a wavelength, and the image and source planes are not very close to the diffractive element.
The presented imaging models are powerful in that they can be used to effectively analyze and simulate any diffractive imaging system operating at the Fresnel regime. The system can contain an on-axis or off-axis photon sieve, Fresnel zone plate, Fresnel phase plate, or other mask-like pattern such as a coded aperture. To illustrate their utility and generality, the derived formulas are applied to several important special cases with common design considerations including photon sieves with circular holes, pixelated diffractive lenses generated by SLM type devices, and binary diffractive lenses with exchanged transparent and opaque regions.
The presented analytical models provide a unified framework to effectively develop new high-resolution imaging modalities enabled by diffractive elements, as well as to analyze their imaging capabilities under various design and observing scenarios. To illustrate this, numerical simulations are performed that show the focusing performance of different FZP-like diffractive lenses that are based on the same binary FZP design. It is observed that all of these diffractive elements (including photon sieves with circular or square holes, and a pixelated version) possess similar PSFs for both focused and defocused cases, with the major difference observed at peak intensity values and at the tails. Hence, the ease of computing the derived PSF formulas enables testing how changing the design parameters affects the resolution and the imaging performance (prior to constructing the physical system), as well as attaining desirable designs that sharpen the PSF. Moreover, using the computed PSFs and the derived input-output relationships involving convolutions, the output of the imaging system can be easily simulated for a given object with coherent or incoherent illumination, using FFTs. These imaging models can also be used as an accurate model of the imaging process for image deconvolution and reconstruction tasks involved in novel computational imaging systems with diffractive lenses [37,38].