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
Diffraction calculus involves the use of Fresnel-Kirchoff or Rayleigh-Sommerfeld integrals in its scalar or vectorial forms . Such calculus is generally performed numerically, which involves a double integral of a complex function with a fast oscillating argument . 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 . 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 . The complex amplitude field distribution of a function g(x,y) at a certain distance is given byFig. 1). Note that are spheres with its center shifted to .
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 seriesEq. (1) is generally approximated by 1/z. Assuming also that 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 , 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 (Eq. (2)) up to second order to obtainEquation (6) has been used in  to obtain analytical expressions for a focused Gaussian beam, and later in , with the first two terms, to study radially polarized beams. Unfortunately, in those forms the resulting expressions are difficult to deal with, particularly due to the existence of in the denominator.
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 respectivein the denominator by a constantwill be enough to simplify the calculus of Eq. (1) without sacrificing its accuracy significantly within a neighborhood of . Thus, our approximation would be simply given byEq. (7).Equation (9) means that the diffracted field is just the Fourier transform of the original function, which is multiplied by the complex function . Note that the phase means tilted spheres. It can be appreciated that no quadratic phase term appears in the resulting diffracted field as it happens in Eq. (4). It can be noticed also that our approximation becomes exact on-axis since .
Expanding 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 and , the ambiguity in choosing 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 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 means that the linear function looks like the function . With a least squares formulation such similarity can be optimized with the following expressionEq. (10) can be computed numerically. When circular apertures are chosen, settingin Eq. (9) provides errors lower than 2% for , being . In following sections the performance of Eq. (11) will be shown.
3. Error analysis in the approximation
To perform our error analysis, for simplicity we will contemplate,and a circular aperture of diameter . Now, considering again the first two terms of Eq. (6), , and the Taylor series expansion of we obtainEq. (8). Taking into account that is the maximum value that can reach x we roughly estimate the maximum error between and like , being . Thus, the higher the or the lower (or ), the higher the exactitude of Eq. (8). An important issue is that errors depends only on, and and not on the aperture D itself. This is an important result because it means that converging waves having equal and different Fresnel numbers must have the same accuracy.
Now we will compute integrals specified by Eqs. (1), (4) and (9) to compare the intensities, and given by the squared modulus of the respective amplitudes. To perform such comparison we have chosen the following parameters: a circular apertureof diameter , , , , and defined by a converging spherical wave of radius equal to z,Fig. 2, where it can be seen that among those PSFs, Fresnel is the most different, as expected.
In Fig. 2, dotted-continuous line corresponds to , continuous line to and dashed to . We show various graphs of corresponding to values of = 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 is just a scale factor of the diffraction pattern. Defining our error of calculus as, we found for the best curve of shown in Fig. 2 that, 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 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 is reduced; for example the previous error is reduced to when 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 and not on D, we fixed = 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. 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 . Additionally, we have obtained curves for = 100 and = 10000 (not shown) finding similar behavior: constant errors independently of D. The curve obtained for = 10000, with errors, 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 Eq. (1), Eqs. (14) and (16) can be approximated with the following substitution8]. Then, we used Eq. (18) to calculate Eqs. (14), and (16) together with the parameters given previously like n = 1, a circular apertureof diameter , , , and a spherical convergent function 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
withFig. 5, as described in section 4, we evaluated Eq. (19) and show in Fig. 6, with dashed lines, the corresponding calculated intensities . The other two curves correspond to those shown in Fig. 5d. It can be appreciated that both approximation, at focus, behaves similarly even that we are using a converging wave with a NA = 0.64. The error obtained with the Debye approximation was which is a little lower than the obtained with our method .
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 (, ,), not at the focus plane but a distance of from it towards the aperture, that is, at.
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 , i.e, , 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 means only different scale factors of the whole diffraction pattern.
Although we have derived a formula to calculate 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 apertureof semi diameter is used and is defined as a rotationally symmetric function, it can be obtained straightforwardly a one dimensional integral form of Eq. (9), given by
, and a Bessel function of the first kind, zero order.
7. Off-axis diffraction analysis
We have demonstrated our approximation for small values of , which means calculations near the optical axis. However our proposal is not restricted to those values since we can define with small, as in previous analysis, and any constant point on the image plane, not necessarily near . 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), , a unit planar wave parallel to the optical axis, and, a constant aperture function. Thus, what we have is the Fourier transform of . Using the following result [11,12]Eq. (7), it is readily shown from Eq. (22) that
From Eq. (24) it can be seen that for , evanescent waves are obtained since and are all real valued and positive. So, Eq. (22) may be helpful to describe the behavior of evanescent waves for general values of .
It can be noticed from Eq. (9) (and those of Eqs. (14-18)) that doing the usual assumption the second term of the evanescent function 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)Eq. (25) with a = 5mm, , and a circular aperture, an intensity error as little as 1.5% is introduced.
Neglecting also the evanescent term in Eq. (21), setting and using the proper converging spherical wave, it turns out to beEq. (26) is important since their zeros positions are strongly related with the definition of the resolution of optical systems.
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.
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).
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]