## Abstract

It is reported an approximation to the scalar and vectorial Rayleigh-Sommerfeld diffraction integrals performed with a Fourier transform operation, as done in Fresnel and Fraunhoffer approximations, but suitable for paraxial, non-paraxial and off-axis regimes. High accuracy is obtained even for waves with low f numbers. The approximation becomes exact on-axis.

© 2017 Optical Society of America

## 1. Introduction

Diffraction calculus involves the use of Fresnel-Kirchoff or Rayleigh-Sommerfeld integrals in its scalar or vectorial forms [1]. Such calculus is generally performed numerically, which involves a double integral of a complex function with a fast oscillating argument [2]. In order to facilitate their calculus or to get some insight about diffraction properties, such integrals are approximated in the so called Fresnel or Fraunhoffer regimes, which are adequate for paraxial regions [3]. For non-paraxial optical systems other approximations have been developed like that due to Debye, well suited for optical systems with high Fresnel numbers [3,4].

So, apparently, there is not a reliable approximation to calculate the diffraction properties of high, intermediate and low numerical aperture optical systems. In this work we will propose an approximation that could fulfill such requirements. We propose to use tilted spheres to approximate the shifted spheres of the diffraction integral; in Fraunhoffer and Fresnel integrals planes and parabolas are used to perform such approximation. The beauty of our approximation is that we perform the calculations with a simple Fourier transform operation, as done with the previous two approximations, but without having the paraxial restrictions. Some immediate novel results, described in this work, obtained as a consequence of our proposal are: a non-paraxial formula for rotationally symmetric functions via a Hankel transform; off-axis diffraction analysis; non-paraxial vectorial formulas performed with Fourier transforms; and even the study of evanescent waves. Since our approximation becomes exact on-axis it is well behaved near the optical axis and could be a useful tool to study optical systems with high focus depth. An important feature is that our approximation is calculated with a Fast Fourier Transform algorithm and 2D images or even 3D volumes are easily obtained.

## 2. Alternative approximation to a diffraction integral

Let us consider the First Rayleigh-Sommerfeld formulation (R-S) in rectangular coordinates, as described by [3]. The complex amplitude field distribution $G({x}_{o},{y}_{o})$ of a function *g(x,y)* at a certain distance $z$ is given by

_{$g(x,y)=f(x,y){g}_{I}(x,y)$}, being ${g}_{I}(x,y)$the illuminating wave and $f(x,y)$ an aperture function; $({x}_{0},{y}_{0})$ are the image coordinates; $\Sigma $ is the diffraction aperture and $k=2\pi /\lambda $, the wavenumber, where $\lambda $ means the wavelength of the illuminating light (see Fig. 1). Note that $\mathrm{exp}(jkR)/R$are spheres with its center shifted to $({x}_{0},{y}_{0})$.

As it is known, instead of a direct evaluation of Eq. (1), which involves various numerical difficulties, usually the R in the argument of the exponential (phase function) is replaced by the first terms of the Taylor series

providing $|x-{x}_{0}|$and $|y-{y}_{0}|$be smaller than*z*(paraxial case) or

*z*be large enough compared with such differences (very low numerical aperture systems). The amplitude term 1/R in Eq. (1) is generally approximated by 1/z. Assuming also that $z>>\lambda $ we obtain the well-known Fresnel approximation (paraxial) of Eq. (1) given by

It can be noticed that the arguments in the exponential functions of Eq. (4) can be interpreted as tilting parabolas. When it is assumed that $z>>r$, such parabolas behaves like planes and we obtain the Fraunhoffer approximation.

In the following descriptions we will use the above explanations for comparison purposes.

Now, Taylor’s formula for a differentiable function of two-variables is used to expand $R$(Eq. (2)) up to second order to obtain

Since it seems to me that no further simplifications to Eq. (6) have been done before, I will show that retaining only the first two terms of Eq. (6) and replacing the respective${R}_{n}$in the denominator by a constant${A}_{z}$will be enough to simplify the calculus of Eq. (1) without sacrificing its accuracy significantly within a neighborhood of $({x}_{0},{y}_{0})$. Thus, our approximation would be simply given by

where${R}_{n}$is given in Eq. (7).Substituting Eq. (8) in Eq. (1) and setting $R\approx {R}_{n}$ (the first term of Eq. (8)) in the denominators we obtain our new approximation to the diffraction integral of Eq. (1) given by

Expanding ${R}_{n}$ of Eq. (8) in Taylor series and using only its first terms in Eq. (9) we reach again, up to a constant factor, the Fresnel and Fraunhoffer approximations, Eq. (4). This way, our approximation can be regarded as a generalization of both.

Since we can define $u={x}_{0}/\lambda {A}_{z}$and $v={y}_{0}/\lambda {A}_{z}$, the ambiguity in choosing ${A}_{z}$ means that the diffracted fields, for a given z, are correctly calculated up to a scale factor. That scale could be unimportant when relative calculations are done or when the shape of the PSF is the goal. In general, we can estimate ${A}_{z}$ either numerically or analytically and their value will depend on the criteria used to calculate it. In order to achieve the calculations it must be noticed that a good choice of ${A}_{z}$ means that the linear function $x/{A}_{z}$ looks like the function $x/{({z}^{2}+{x}^{2}+{y}^{2})}^{1/2}$. With a least squares formulation such similarity can be optimized with the following expression

## 3. Error analysis in the approximation

To perform our error analysis, for simplicity we will contemplate${x}_{0}\ne 0$,${y}_{0}=0$and $\Sigma $ a circular aperture of diameter $D$_{.} Now, considering again the first two terms of Eq. (6), ${R}_{n}-(x{x}_{0}+y{y}_{0})/{R}_{n}$, and the Taylor series expansion of $x/{R}_{n}$we obtain

*x*we roughly estimate the maximum error between $x{x}_{0}/{R}_{n}$ and $x{x}_{0}/z$ like $o({x}_{0}/{f}_{num}^{3})$, being ${f}_{num}=z/D$. Thus, the higher the ${f}_{num}$or the lower ${x}_{0}$ (or ${y}_{0}$), the higher the exactitude of Eq. (8). An important issue is that errors depends only on${f}_{num}$, ${x}_{0}$ and ${y}_{0}$and not on the aperture D itself. This is an important result because it means that converging waves having equal ${f}_{num}$ and different Fresnel numbers $Fn=a/(2\lambda {f}_{num})$ must have the same accuracy.

Now we will compute integrals specified by Eqs. (1), (4) and (9) to compare the intensities${I}_{RS}({x}_{0},{y}_{0})$, ${I}_{F}({x}_{0},{y}_{0})$and ${I}_{N}({x}_{0},{y}_{0})$ given by the squared modulus of the respective amplitudes. To perform such comparison we have chosen the following parameters: a circular aperture$\Sigma $of diameter $D=5\text{\hspace{0.17em}}mm$, $\lambda =0.6328x{10}^{3}\text{\hspace{0.17em}}mm$, $z=0.6D$, $f(x,y)=1$, and ${g}_{I}(x,y)$defined by a converging spherical wave of radius equal to *z*,

In Fig. 2, dotted-continuous line corresponds to ${I}_{RS}$, continuous line to ${I}_{F}$ and dashed to ${I}_{N}$. We show various graphs of ${I}_{N}$ corresponding to values of ${A}_{z}$ = 0.85z, 1.0z, 1.08z and 1.2z, being 1.2z the most accurate (this value was computed with Eq. (11)). With this graphics it is clarified that ${A}_{z}$ is just a scale factor of the diffraction pattern. Defining our error of calculus as$E=\mathrm{max}|({I}_{RS}/\mathrm{max}{I}_{RS})-({I}_{N}/\mathrm{max}{I}_{N})|$, we found for the best curve of ${I}_{N}$ shown in Fig. 2 that$E=1.85\%$, which may be considered small, taking into account that a high numerical wave has been analyzed. In addition to this, we show in Fig. 3 a curve specifying that errors decreases when ${f}_{num}$increases, as predicted with the error analysis of Eq. (12). Errors of Fig. 3 can be reduced when the analysis range for the spatial coordinates ${x}_{0}$ is reduced; for example the previous error $E({f}_{num}=0.6)=1.85\%$ is reduced to $E=1\%$ when $\mathrm{max}(|{x}_{0}|)\le 1AiR$was considered. This result is expected since our approximation is truly exact on axis.

With the purpose to show that effectively the exactitude of Eq. (9) depends only on ${f}_{num}$ and not on D, we fixed ${f}_{num}$ = 0.6, 1.0 and 2.0 and varied the corresponding aperture diameters D from 0.05mm to 1000mm, calculating the respective diffraction intensities and estimating the respective errors as done in curve of Fig. 3. Results are shown in Fig. 4 where it can be seen that, effectively, errors do not depend on the aperture diameter itself, being constant for the same${f}_{num}$. Logarithmic scale (natural logarithm) has been used in order to show the wide range of values of D considered. Notice that log(D) = 3 in Fig. 4 means Fresnel number $Fn\approx ({10}^{5})$. Additionally, we have obtained curves for ${f}_{num}$ = 100 and ${f}_{num}$ = 10000 (not shown) finding similar behavior: constant errors independently of D. The curve obtained for ${f}_{num}$ = 10000, with errors$E\approx 0.001\%$, comprises analysis when Fresnel number goes from 0.002 to 39.

From these results we conclude that our proposal can be used confidently in the analysis of low, medium or high numerical aperture optical systems with errors on the order of 2%.

## 4. Vectorial diffraction integrals

One application of our approximation is in the study of diffraction properties of polarized beams. Since the first kind vectorial Rayleigh-Sommerfeld formulas expressed in rectangular coordinates are given by [7]

*n*is the refractive index of the medium space; ${E}_{x}(x,y)$and ${E}_{y}(x,y)$denote the electric-field components along the x and y directions, respectively and ${E}_{z}(x,y,{x}_{0},{y}_{0})={E}_{x}(x,y)(x-{x}_{0})+{E}_{y}(x,y)(y-{y}_{0})$. Clearly,

*n = 1,*a circular aperture$\Sigma $of diameter $D=5\text{\hspace{0.17em}}mm$, $\lambda =0.6328x{10}^{3}\text{\hspace{0.17em}}mm$, $z=0.6D$, ${A}_{z}=1.2z$and a spherical convergent function${g}_{I}(x,y)$ given by Eq. (13) . The results of using our approximation with a fast Fourier algorithm are presented in Figs. 5(a-c) and using the integration routines of MatLab in Fig. 5d, where the comparison with exacts integral is shown. Good agreement has been obtained again, with

*E = 2.2%.*

## 5. Comparison with Debye approximation

Using the results of [9,10] we can write the vectorial diffraction integrals within the Debye approximation like

with

It is interesting to study the behavior of our proposal when defocus is introduced. Without going deeper into this matter we show Fig. 7, obtained with the previous data of Fig. 6 ($D=5\text{\hspace{0.17em}}mm$, $\lambda =0.6328x{10}^{3}\text{\hspace{0.17em}}mm$,${A}_{z}=1.2z$), not at the focus plane but a distance of $f/60=0.05\text{\hspace{0.17em}}mm$ from it towards the aperture, that is, at$z=0.59D$.

It can be appreciated that at this defocused position our results (dotted-continuous line) were more accurate than those obtained with the Debye approximation (dashed line). When choosing another value for ${A}_{z}$, i.e, ${A}_{z}=1.29z$, the curves of Fig. 8 were obtained; in that figure it can be noticed that our approximation (I_N) was improved. At this point it must be emphasized that different values for ${A}_{z}$means only different scale factors of the whole diffraction pattern.

Although we have derived a formula to calculate ${A}_{z}$ for a focused plane, Eq. (11), a general expression must be found yet when the diffraction patterns are imaged onto others plane positions.

## 6. Diffraction integrals for rotationally symmetric functions

When a circular aperture$\Sigma $of semi diameter $a$ is used and $g(x,y)=g(r=\sqrt{{x}^{2}+{y}^{2}})$ is defined as a rotationally symmetric function, it can be obtained straightforwardly a one dimensional integral form of Eq. (9), given by

_{${R}_{r}=\sqrt{{z}^{2}+{r}^{2}}$}, ${r}_{0}=\sqrt{{x}_{0}^{2}+{y}_{0}^{2}}$and ${J}_{0}$ a Bessel function of the first kind, zero order.

## 7. Off-axis diffraction analysis

We have demonstrated our approximation for small values of $({x}_{0},{y}_{0})$, which means calculations near the optical axis. However our proposal is not restricted to those values since we can define $({x}_{n0},{y}_{n0})=({x}_{n}+{x}_{0},{y}_{n}+{y}_{0})$with $({x}_{0},{y}_{0})$ small, as in previous analysis, and $({x}_{n},{y}_{n})$any constant point on the image plane, not necessarily near $(0,0)$. Then, it can be shown that our approximation can be described by

## 8. Diffraction integrals using planar illumination

It seems that our approximation could be used to analyze diffraction apertures illuminated with planar waves. To show this, let us define in Eq. (22), ${g}_{I}(x,y)=1.0$, a unit planar wave parallel to the optical axis, and$f(x,y)=1.0$, a constant aperture function. Thus, what we have is the Fourier transform of $(jk-1/{R}_{n0})(1/{R}_{n0}^{2})\mathrm{exp}(jk{R}_{n0})$. Using the following result [11,12]

From Eq. (24) it can be seen that for ${({x}_{0}^{2}+{y}_{0}^{2})}^{1/2}\ge {B}_{z}$, evanescent waves are obtained since $k,\text{\hspace{0.17em}}\text{\hspace{0.17em}}z$ and${B}_{z}$ are all real valued and positive. So, Eq. (22) may be helpful to describe the behavior of evanescent waves for general values of $g(x,y)$.

## 9. Discussions

It can be noticed from Eq. (9) (and those of Eqs. (14-18)) that doing the usual assumption$z>>\lambda $ the second term of the evanescent function $(jk-1/{R}_{n})$can be neglected. If in addition to this we suppose that the illuminating light corresponds to a converging spherical wave behaves like Eq. (13) we obtain from Eq. (9)

*a = 5mm*, ${f}_{num}=1$, $f(x,y)=1$ and a circular aperture, an intensity error as little as 1.5% is introduced.

Neglecting also the evanescent term in Eq. (21), setting $f(r)=1$ and using the proper converging spherical wave, it turns out to be

It is convenient to mention that, according with the explanation given in section 7, we can calculate the diffraction patterns in a neighborhood of each out-off axis point, and by a stitching process we could obtain such diffraction, with the desired exactitude, in bigger areas beyond those neighborhoods at the cost of performing various Fourier transforms.

## 10. Conclusions

We have proved that an alternative approximation can be used to analyze paraxial and non-paraxial systems without restriction on the Fresnel number. The approximation simplifies scalar and vectorial diffraction integrals to be calculated with a Fourier transform operation. An expression for off-axis analysis is obtained as well. Consequences of our approximation like one-dimensional integral for rotationally symmetric functions or evanescent waves have been considered in the analysis. Although we have focused our approximation into RS a similar analysis can be carried out with Fresnel-Kirchhoff integrals since they share similar expressions.

## References and links

**1. **M. Born and E. Wolf, *Principles of Optics*, 7th ed. (Cambridge University Press, 1999).

**2. **J. J. Stamnes, *Waves in Focal Regions* (IOP Publishing Limited, 1986).

**3. **J. W. Goodman, *Introduction to Fourier Optics*, 2nd ed. (McGraw-Hill, 1996).

**4. **C. J. R. Sheppard, “Limitations of the paraxial Debye approximation,” Opt. Lett. **38**(7), 1074–1076 (2013). [CrossRef] [PubMed]

**5. **A. Ciattoni, B. Crosignani, and P. Di Porto, “Vectorial analytical description of propagation of a highly nonparaxial beam,” Opt. Commun. **202**(1-3), 17–20 (2002). [CrossRef]

**6. **D. Deng, “Nonparaxial propagation of radially polarized light beams,” J. Opt. Soc. Am. B **23**(6), 1228–1234 (2006). [CrossRef]

**7. **H. Ye, C.-W. Qiu, K. Huang, J. Teng, B. Luk’yanchuk, and S. P. Yeo, “Creation of a longitudinally polarized subwavelength hotspot with an utra-thin planar lens: vectorial Rayleigh-Sommerfeld method,” Laser Phys. Lett. **10**(6), 1–8 (2013). [CrossRef]

**8. **X. Hao, C. Kuang, T. Wang, and X. Liu, “Effects of polarization on the de-excitation dark focal spot in STED microscopy,” J. Opt. **12**(11), 115707 (2010). [CrossRef]

**9. **B. Richards and E. Wolf, “Electromagnetic diffraction in optical systems II. Structure of the image field in an aplanatic system,” Proc. R. Soc. Lond. **253**(1274), 358–379 (1959). [CrossRef]

**10. **L. E. Helseth, “Focusing of atoms with strongly confined light potentials,” Opt. Commun. **212**(4-6), 343–352 (2002). [CrossRef]

**11. **E. Lalor, “Conditions for the validity of the angular spectrum of plane waves,” J. Opt. Soc. Am. **58**(9), 1235–1237 (1968). [CrossRef]

**12. **W. H. Southwell, “Validity of the Fresnel approximation in the near field,” J. Opt. Soc. Am. **71**(1), 7–14 (1981). [CrossRef]