## Abstract

Due to analytical and numerical difficulties, the propagation of optical fields in any state of spatial coherence is traditionally computed under severe approximations. The paraxial approach in the Fresnel–Fraunhofer domain is one of the most widely used. These approximations provide a rough knowledge of the actual light behavior as it propagates, which is not enough for supporting applications, such as light propagation under a high numerical aperture (NA). In this paper, a non-approximated model for the propagation of optical fields in any state of spatial coherence is presented. The method is applicable in very practical cases, as high-NA propagations, because of its simplicity of implementation. This approach allows for studying unaware behaviors of light as it propagates. The light behavior close to the diffracting transmittances can also be analyzed with the aid of the proposed tool.

© 2011 OSA

## 1. Introduction

The exact description of the propagation of a stationary random optical field in any state of spatial coherence was formalized by Zernike in an equation that involves the mutual intensity of the field on two surfaces in space [1,2]. That formula has no restrictions in the shapes and sizes of the surfaces, neither the distances between points on them nor the involved inclination angles. By accounting each spectral component of the optical field, Zernike’s formula can be replaced by an expression introduced by Wolf [2,3] that describes the field in terms of its cross-spectral densities on the considered surfaces. Wolf’s equation has the above mentioned attributes of Zernike’s formula as well, and it has become the most used tool for modeling stationary random optical fields in any state of spatial coherence.

However, its exact calculation is very involved, and therefore, some strategies have been applied in order to approach it to specific practical situations. The most known of them is the paraxial approach [2,3], which is very useful when the propagation occurs in the Fresnel–Fraunhofer domain, as in conventional image formation for instance. The stationary-phase approach [3,4] is also used in the synthesis of diffractive optical elements [4]. Nevertheless, non-paraxial modeling and calculations are more required nowadays for supporting the novel developments, such as in digital holography microscopy [5] and beam shaping [6] for micro-lithography or optical tweezers.

To this aim, an alternative modeling strategy is proposed in this work. It is based on the non-paraxial definition of the marginal power spectrum [7], obtained from Wolf’s exact formula for the propagation of the cross-spectral density. It allows modeling optical fields in any state of spatial coherence in terms of carrier and modulating 0-*π* rays that respectively propagate the radiant and the modulating energies of the field, emitted by point sources distributed on two different layers associated to the emitting surface, named the radiant and the virtual layers [8]. This strategy leads to a very efficient numerical tool [9] that allows for computing the light propagation in domains where the Fresnel–Fraunhofer approximation does not apply. Specifically, the exact calculation of the power spectrum distribution produced on the observation plane is analyzed. For the sake of simplicity and conciseness, the presented modeling method is illustrated with one-dimensional examples.

## 2. Exact calculation of the marginal power spectrum

Wolf’s formula for the exact propagation of the cross-spectral density of an optical field at the frequency $\nu $ from the aperture plane (AP) to the observation plane (OP), placed at the arbitrary distance *z* to each other, is [3]

By assuming the center and difference coordinates $\left({\xi}_{A},{\xi}_{D}\right)$ on the AP and $\left({r}_{A},{r}_{D}\right)$ on the OP [7], the points ${Q}_{1,2}$ and ${P}_{1,2}$ are denoted as ${\xi}_{A}\pm {\xi}_{D}/2$ and ${r}_{A}\pm {r}_{D}/2$, respectively, with + corresponding to the suffix 1 and – to the suffix 2 in all the expressions. Consequently, ${s}_{\pm}=\left|z+{r}_{A}\pm {r}_{D}/2-\left({\xi}_{A}\pm {\xi}_{D}/2\right)\text{\hspace{0.17em}}\right|$ holds, with **z** being the vector of magnitude *z,* which is normal to both the AP and the OP. Furthermore, the cross-spectral density at the AP takes the form $W\left({Q}_{1},{Q}_{2};\nu \right)=W\left(+,-;\nu \right)=\sqrt{{S}_{0}(+)}\text{\hspace{0.17em}}t(+)\text{\hspace{0.17em}}\sqrt{{S}_{0}(-)}\text{\hspace{0.17em}}{t}^{*}(-)\text{\hspace{0.17em}}\mu \left(+,-\right)$ [3], with $\pm $ being a short notation for ${\xi}_{A}\pm {\xi}_{D}/2$, ${S}_{0}(\pm )$ being the illumination power spectrum, $t(\pm )=\left|\text{\hspace{0.17em}}t(\pm )\text{\hspace{0.17em}}\right|\text{\hspace{0.17em}}\mathrm{exp}[i\text{\hspace{0.17em}}\phi (\pm )]$ being the complex transmission and $\mu \left(+,-\right)=\left|\mu \left(+,-\right)\text{\hspace{0.17em}}\right|\text{\hspace{0.17em}}\mathrm{exp}\left[i\text{\hspace{0.17em}}\alpha \left(+,-\right)\right]$ being the complex degree of spatial coherence at the AP [2,3] related to the positions ${\xi}_{A}\pm {\xi}_{D}/2$ on the AP. After replacing these definitions in Eq. (1) and taking into account that $d{\sigma}_{1}\text{\hspace{0.17em}}d{\sigma}_{2}={d}^{2}{\xi}_{A}\text{\hspace{0.17em}}{d}^{2}{\xi}_{D}$ holds, it yields the cross-spectral density at the OP in center and difference coordinates $W\left({r}_{A}+{r}_{D}/2,{r}_{A}-{r}_{D}/2;\nu \right)$. Thus, the power spectrum at the OP is obtained by evaluating the cross-spectral density there for ${P}_{1}={P}_{2}$, i.e., by evaluating Eq. (1) in center and difference coordinates for ${r}_{D}=0$, that is

Equation (3) is the exact definition of the marginal power spectrum [7] for the propagation of optical fields in any state of spatial coherence from the AP to the OP, along $z\ge 0$. It has no restrictions in the sizes of the illuminated areas on both the AP and the OP, neither the distances between points on such planes or the involved inclination angles. Its integration region is the smallest area between the aperture and the structured spatial coherence support [10] centered at ${\xi}_{A}$. The geometrical parameters involved in Eq. (3) are illustrated in Fig. 1 . In the following, results obtained from Eq. (3) and from the conventional paraxial approached calculation in the Fresnel–Fraunhofer domain will be compared with each other. Furthermore, unaware behavior of the light propagation, even from the vicinities of the AP, will be analyzed with use of the non-approximated formulation.

## 3. Diffraction features shown by the exact calculation software

#### 3.1 The free-space diffraction envelope

The expression $t(+)\text{\hspace{0.17em}}{t}^{\ast}(-)=\delta \left({\xi}_{A}-a\right)\text{\hspace{0.17em}}\delta \left({\xi}_{D}\right)$, with $\delta (\u2022)$ being the Dirac’s delta, denotes an opaque mask with a pinhole placed at the position ${\xi}_{A}=a$ at the AP, which isolates a single radiant point source. Accordingly, Eqs. (3) and (2) yield

In contrast, the exact calculation denoted by Eq. (5) predicts the evolution of the power spectrum from a delta-like function for $z=0$ to a Lorentzian-like profile for *z* > 0. Both distributions coincide only in the paraxial region, i.e., within a relatively slight region around the optical axis. It is illustrated by the comparison shown in Fig. 2
for $\lambda =0.632\mu m$ and $z={10}^{4}\mu \text{\hspace{0.17em}}m$. Figure 3(a)
illustrates the evolution of the power spectrum at the OP within the interval $0\le z\le 0.1\mu \text{\hspace{0.17em}}m$. The Lorentzian-like profile of the power spectrum remains shape-invariant for the light propagation along arbitrary distances, but it spreads and its maximum decays follow the $1/{z}^{2}$ law [2]; thus, it diffracts because of its propagation in free space, which is quite different from the paraxial approach prediction (Fig. 2). Therefore, it is reasonable to call the Lorentizan-like profile a *free-space diffraction envelope*, because it depends only on the geometrical parameters of the propagation. Figure 3(b) shows the comparison between the free-space diffraction envelope and conventional Lorentzian and Gaussian distributions. The free-space diffraction envelope deviates about 1.5% from the Lorentzian function and about 5% from the Gaussian distribution.

It is worth noting that the free-space diffraction envelope spreads within a specific angle, which can be estimated by determining the encircled energy curves [2] of the power spectrum at different planes. Figure 3(c) shows the profiles of the encircled energy curves for diffraction of light $\left(\lambda =0.632\mu \text{\hspace{0.17em}}m\right)$ emitted by a single radiant point source on the OP at $z=1,\text{\hspace{0.17em}}10,\text{\hspace{0.17em}}20,\text{\hspace{0.17em}}30\text{\hspace{0.17em}}\mu \text{\hspace{0.17em}}m$, respectively. Then, the aperture angle of the cone of light emitted by the point source can be determined by specifying the size of the region of the power spectrum at a given plane that concentrates at least 95% of the total energy. Denoting the radius of this region as *R*, which is directly provided by the encircled energy curve for the considered *z*, the aperture angle of the cone is $\vartheta =\mathrm{arctan}\left(R/z\right)$. It was proved that about 95% of the total energy spread within an semi-angle of 75.9° for $1\text{\hspace{0.17em}}\mu \text{\hspace{0.17em}}m\le z\le {10}^{4}\mu \text{\hspace{0.17em}}m$.

The results above lead to the conclusion that the free-space diffraction envelope determines a finite observation window for the power spectrum emitted by a single radiant point source at any OP, which has a very physical meaning with an important effect on diffraction of light through finite apertures, as discussed later. The paraxial approach is not able to predict such results.

#### 3.2 Spatial frequency modulation by interference

The term spatial frequency is well-known and applied in paraxial optics [11] and can be physically evidenced by means of Young’s experiment with radiant point sources. In this case, the power spectrum at the OP is given by

Equation (6) points out the unique spatial frequency of the paraxial field, which is corresponding to the frequency of the interference fringes determined by the well-defined parameter $f=k\text{\hspace{0.17em}}\left|b\right|/z$ in the argument of the cosine modulation. The remainder terms of the cosine argument determine the shift of the modulation with respect to the origin of the ${r}_{A}$-coordinate at the OP.

For the non-paraxial model, the Young’s mask is represented by $t(+)\text{\hspace{0.17em}}\text{\hspace{0.17em}}{t}^{*}(-)=$ $\delta \left({\xi}_{D}\right)\text{\hspace{0.17em}}\left[\delta \left({\xi}_{A}-a+b/2\right)\text{\hspace{0.17em}}+\text{\hspace{0.17em}}\delta \left({\xi}_{A}-a-b/2\right)\right]+\delta \left({\xi}_{A}-a\right)\text{\hspace{0.17em}}\left[\delta \left({\xi}_{D}+b\right)\text{\hspace{0.17em}}+\text{\hspace{0.17em}}\delta \left({\xi}_{D}-b\right)\right]$, where **a** is the midpoint between the two pinholes and **b** is their separation vector. Consequently, Eqs. (3) and (2) yield

It is worth noting that this cone of light is different from that used for characterizing the spatial frequency in Fourier optics [11]. The vertex and basis of the light cone in the analysis above are at the AP and the OP, respectively, while they are at the OP and the AP, respectively, for the light cone used in Fourier optics. Thus, the angle of the last cone grows as the separation of the sources increases, and consequently the spatial frequency grows as well [11]. Nevertheless, this growth of the spatial frequency is also manifested in the increase of the spatial frequency chirping in the exactly calculated light cone.

Figure 4 illustrates this statement a for fully spatially coherent Young’s experiment in the paraxial region for $z={10}^{4}\mu \text{\hspace{0.17em}}m$, with $\lambda =0.632\mu \text{\hspace{0.17em}}m$. Each row shows (a) the exactly calculated and (b) the paraxial approached marginal power spectra [7] and (c) the comparison between the respective power spectra. The top row is corresponding to the low (quasi-null) spatial frequency, i.e., the source separation if much smaller than the wavelength $\left(\text{\hspace{0.17em}}\left|\text{\hspace{0.17em}}b\right|=0.1\mu \text{\hspace{0.17em}}m\right)$ and there is no cosine modulation in the paraxial region at the OP. The source separation is larger than the wavelength for the middle and the bottom rows, but the spatial frequency is higher for the last one, so that the source separation and the cosine modulation in the paraxial region at the OP is larger.

Young’s experiment was modeled with two identical point sources separated $5\mu \text{\hspace{0.17em}}m$ to each other, which emit an optical field of $\lambda =0.632\mu \text{\hspace{0.17em}}m$. After assuming a complex degree of spatial coherence with magnitude 0.3 and phase *π*, Fig. 5(a)
shows the exact calculation of the power spectrum propagation along the first $100\mu \text{\hspace{0.17em}}m$ on the *z* axis, and Fig. 5(b) sketches the exactly calculated profile of the power spectrum at the OP in the Fresnel–Fraunhofer domain $\left(z={10}^{4}\mu \text{\hspace{0.17em}}m\right)$. Bright and dark fringes denote the regions of constructive and destructive interference, respectively; their visibility and the fact that a dark fringe is aligned with the *z* axis are determined by the complex degree of spatial coherence. The decay of the fringe brightness and the increase in the fringe width from the cone center to the cone edge is due to the free-space diffraction envelope and to the spatial frequency chirping, respectively.

A comparison between the exactly calculated and the paraxial approached power spectra in the Fresnel–Fraunhofer domain is shown in Fig. 5(c) for fully spatially coherent sources separated $50\mu \text{\hspace{0.17em}}m$ to each other. It allows estimating the size of the paraxial region in which the spatial frequency remains fixed as conventionally defined.

#### 3.3 Diffraction with discrete sets of point sources

An important feature of the proposed model is its compatibility with the field representation in terms of two separated layers of point sources, named the radiant and the virtual layers, respectively [8]. The former one is responsible for the emission of the radiant energy of the field, independently of its spatial coherence state, while the latter one is responsible for the emission of modulating energy in close relationship with the spatial coherence state of the field. This representation allows developing very efficient algorithms of calculation and simulation of Fresnel–Fraunhofer interference and diffraction of optical fields in any state of spatial coherence [9]. Now, it is used in order to solve the exact calculation of interference and diffraction efficiently.

In this context, the marginal power spectrum can be expressed as $S\left({\xi}_{A},{r}_{A};\nu \right)=$ ${S}_{rad}\left({\xi}_{A},{r}_{A};\nu \right)+{S}_{virt}\left({\xi}_{A},{r}_{A};\nu \right)$, and following Eq. (2), the power spectrum at the OP takes the form $S\left({r}_{A};\nu \right)={S}_{rad}\left({r}_{A};\nu \right)+{S}_{virt}\left({r}_{A};\nu \right)$ [7,8]. Therefore,

*N*point sources with pitch

**b**, allocated in the radiant layer. The power spectrum that this array provides at the OP takes the form

In order to determine the distribution of point sources on the virtual layer, the concept of *classes of radiator pairs* [12] is used. Any member of a specific class provides the characteristic Young cosine-like modulation of the class, which is a component of the modulating energy emitted by the virtual point source placed at the midpoint of such member. Thus,

*m*th class of radiator pairs, i.e., the class with separation vector $m\text{\hspace{0.17em}}b$ with $1\le m\le N-1$; then, the contributed power spectrum at the OP provided by this class is

The parameter $\beta $ in Eq. (12) takes only the values $\beta =1/2$ for pure virtual point sources, placed at the midpoints between consecutive radiant point sources, and $\beta =0$ for virtual components of dual point sources, i.e., virtual point sources that coincide in position with the radiant point sources [8].

The values of the indices *P* and *Q* depend on the number of radiator pairs in each class, and their sequences are illustrated in Table 1
for 15 classes constituted by 16 radiant point sources, i.e., $1\le m\le 15$ with separation vectors $m\text{\hspace{0.17em}}b$, that is 29 virtual point sources turned on when all the radiant point sources are correlated in some extent.

Figure 6
illustrates the modeling of Fraunhofer interference and diffraction $\left(z={10}^{4}\mu \text{\hspace{0.17em}}m\right)$ with seven fully spatially coherent radiant point sources of total length *L* and pitch $\left|\text{\hspace{0.17em}}b\right|$. This array contains 6 classes of radiator pairs that turn on 11 virtual point sources. The power spectra at the end of the top and the middle rows show the expected distribution of 5 secondary maxima between consecutive main maxima. The effects of the free-space diffraction envelope and the spatial frequency chirping are apparent in the decay of the peak powers and the peak spreading toward the patterns’ edges. The increase in the high spatial frequency components as the array pitch increases yields the shifting of the main maxima of the patterns toward the pattern edge. Consequently, the maxima corresponding to positions beyond the cutoff of the free-space diffraction envelope are dropped out, as shown by such power spectra profiles. In other words, the free space behaves as a Lorentzian low-pass spatial filter.

It is worth noting that both the pitch and total length of the array are greater that the wavelength in such panels. If the pitch becomes shorter but the total length remains longer than the wavelength, then only the central main maximum and some secondary maxima appear in the power spectrum pattern at the OP. It is corresponding to the diffraction pattern of a uniform and spatially coherent plane wave produced by a slit of width equal to the total length of the radiant source array. Indeed, Fig. 7(a) shows a comparison between the exactly calculated power spectrum pattern and the squared circular sinus obtained by paraxial approach in the Fresnel–Fraunhofer domain under the same parameters. The fit of those results is significantly improved by increasing the radiant point source density in the array, as shown in Fig. 7(b). In this case, 16 radiant point sources are distributed along the array of $3\mu \text{\hspace{0.17em}}m$ length, so that the pitch is reduced to $0.2\mu \text{\hspace{0.17em}}m$. This array contains 15 classes of radiator pairs that turn on 29 point sources on the virtual layer.

This result points out that the experimental diffraction patterns can be modeled as produced by a discrete set of relative few high correlated radiant point sources, and pitch smaller than the wavelength, instead of a continuum wavefront. Furthermore, if the array length becomes smaller than the wavelength, then the diffraction pattern will be determined by the free-space diffraction envelope in the same way as the light emitted by a single radiant point source, as shown in Fig. 7(c). This feature is crucial in order to develop efficient algorithms for numerical calculations and simulations.

#### 3.4 Modulating energy at the AP

Equation (2) refers to an important feature of the proposed model, because it is applicable to any distance *z* between the AP and the OP, including $z=0$ which is corresponding to the situation in which the OP coincides with the AP. It means that it allows determining both the power distribution and the correlation properties of the optical field at the AP without appealing to the Fourier transform methods, conventionally used in the paraxial approach. Specifically, it allows estimating the distribution of the modulating energy at the AP, which reveals the correlation state of the radiant point sources at this plane.

Figure 8 panel (a) (top figure) shows the exactly calculated marginal power spectrum $\left(\lambda =0.632\mu \text{\hspace{0.17em}}m\right)$ at the mask plane of a Young’s experiment with two fully spatially coherent radiant point sources, separated $5\mu \text{\hspace{0.17em}}m$ to each other. The corresponding power spectrum is shown in the bottom figure. It has the expected distribution for the two radiant point sources at the mask plane. Figure 8 panels (b) and (c) show the profiles of the modulating energy associated to the virtual point source placed at the midpoint between the radiant point sources, which is the center of the structured spatial coherence support that includes the pair of radiant point sources. It means that this energy spreads over an extended region of the mask plane although it is addressed at the position of the virtual point source. The main peaks of the modulating energy select the correlated pair of radiant point sources, in such a way that its height and sign indicate the magnitude and phase of their correlation respectively. Indeed, panel (b) (top figure) shows the profile of the modulating energy when the phase difference of the emissions of the radiant point sources is $2n\pi $, with $n=0,\text{\hspace{0.17em}}1,\text{\hspace{0.17em}}\cdots $, and therefore both peaks take positive values. Such phase difference is $\left(2n+1\right)\pi $ in the bottom profile, where both peaks take negative values. In panel (c) (top figure) the phase difference is $\left(2n+1/2\right)\pi $, and then the value of the peak on the left is positive, while that of the peak on the right is negative; in the bottom figure the distribution is reversed because the phase difference is $\left(2n+3/2\right)\pi $. It is worth noting that this feature cannot be predicted by paraxial approached calculations.

## 4. Conclusion

An efficient strategy for the exact calculation of the propagation integral of optical fields in any state of spatial coherence was shown. It is based on the exactly numerically calculated marginal power spectrum and does not have approximations nor restrictions in the sizes of the illuminated areas at both the AP and the OP, the distances between such planes, the apertures of the involved light cones, and the inclination factors. This modeling strategy reveals diffraction features that cannot be predicted or modeling by the conventional paraxial approach, such as the Lorentzian-like free-space diffraction envelope, the spatial frequency modulation by interference (Young’s fringe frequency chirping), the diffraction with discrete sets of point sources, and the modulating energy at the AP. Such features point out that (i) the free space behaves as a Lorentzian low-pass spatial filter, (ii) the set of point sources on the radiant layer at the AP must be discrete, because pure virtual point sources are necessary in order to accurately perform the mentioned exact calculations (moreover, discrete sets of few radiant point sources are able to reproduce the diffraction patterns of continuous wavefronts), and (iii) the calculation for $z=0$ is performed by evaluation of the distance *z* instead of the application of the conventional Fourier transform methods. The modulating energy in this case points out interesting details of the source correlations. The simplicity and the efficiency of the numerical algorithm are remarkable. It makes this strategy especially useful for modeling practical applications beyond the conventional paraxial approach in the Fresnel–Fraunhofer domain, such as partially coherent imaging, digital holography microscopy, and beam shaping.

## Acknowledgments

This work was partially supported by the Patrimonio Autónomo Fondo Nacional de Financiamiento para la Ciencia, la Tecnología y la Innovación, Francisco José de Caldas, Colciencias Grant number 111852128322, and by the Universidad Nacional de Colombia, Vicerrectoría de Investigación grants numbers 12932 and 12934. The authors also acknowledge the support of DIME (Dirección de Investigación Medellín UNAL) and DINAIN (Dirección Nacional de Investigación, UNAL).

## References and links

**1. **F. Zernike, “The concept of Degree of Coherence and its application to optical problems,” Physica **5**(8), 785–795 (1938). [CrossRef]

**2. **M. Born and E. Wolf, *Principles of Optics*, 6th. ed. (Pergamon Press, 1993), Chap. 10.

**3. **L. Mandel and E. Wolf, *Optical Coherence and Quantum Optics* (Cambridge University Press, 1995), Chap. 4.

**4. **Z. Jaroszewicz, *Axicons: Design and Propagation Properties* (Research and development treatises, SPIE Polish chapter) (SPIE, 1997), Vol. 5.

**5. **D. C. Alvarez-Palacio and J. Garcia-Sucerquia, “Lensless microscopy technique for static and dynamic colloidal systems,” J. Colloid Interface Sci. **349**(2), 637–640 (2010). [CrossRef] [PubMed]

**6. **K. Jahn and N. Bokor, “Intensity control of the focal spot by vectorial beam shaping,” Opt. Commun. **283**(24), 4859–4865 (2010). [CrossRef]

**7. **R. Castañeda, “The optics of spatial coherence wavelets, in *Advances in Imaging and Electron Physics*, P. Hawkes, ed. (Academic Press, 2010), Vol. 164.

**8. **R. Castañeda, G. Cañas-Cardona, and J. Garcia-Sucerquia, “Radiant, virtual, and dual sources of optical fields in any state of spatial coherence,” J. Opt. Soc. Am. A **27**(6), 1322–1330 (2010). [CrossRef] [PubMed]

**9. **R. Castañeda, H. Muñoz-Ossa, and J. Garcia-Sucerquia, “Efficient numerical calculation of interference and diffraction of optical fields in any state of spatial coherence in the phase-space representation,” Appl. Opt. **49**(31), 6063–6071 (2010). [CrossRef]

**10. **R. Castañeda, H. Muñoz, and G. Cañas-Cardona, “The structured spatial coherence support,” J. Mod. Opt. , **58**(11) (2011), doi: . [CrossRef]

**11. **K. Iizuka, *Engineering Optics* (Springer Verlag, 1985).

**12. **R. Castañeda and J. García, “Classes of source pairs in interference and diffraction,” Opt. Commun. **226**(1-6), 45–55 (2003). [CrossRef]