1. Introduction

The study of propagation of light continues to play a central role in Optics and Photonics. The increasing complexity of modern optical systems poses new challenges to light/matter interaction modeling. In particular, sharp-edge diffraction represents a key problem to be tackled whenever light is left to pass through a linear optical system. Pupils, filters, diaphragms, lenses, unavoidably limit the transverse distribution of the incoming wavefield. As a consequence, edge diffraction effects at all relevant boundaries have to take into account in order for the field emerging from the exit pupil to be adequately characterized (in amplitude and phase).

Sharp-edge diffraction is as old as wave theory of light: in 1802, Thomas Young first suggested the idea that the rim of an illuminated aperture could act as a secondary light source [1,2]. Due to their physical appeal, Young’s ideas received a continuous, growing attention. But only after the pioneering works by Maggi [3] and by Rubinowicz [4], these ideas have definitely found a quantitative formulation in the form of the so-called BDW theory [5, Ch. 8], thought for spherical and/or plane wave diffraction by arbitrarily shaped sharp-edge apertures (or obstacles). The first attempt of extending BDW’s theory to deal with general impinging wavefields was done by Miyamoto and Wolf at the beginning of sixties [6,7]. Despite its formal elegance, the practical applicability of BDW/Miyamoto/Wolf’s theory has mostly been limited to Gaussian beams [8–12]. To develop a more manageable theory, paraxial approximation was then invoked from the beginning. In this way, a “genuinely paraxial” version of the original Young/Maggi/Rubinowicz theory was proposed in [13,14], based on some results published in [15,16]. This paraxial revisitation of BDW theory soon revealed its predictive potential [17,18], especially once placed within the context of the so-called Catastrophe Optics [19,20]. Later on, further generalizations aimed at dealing with sharp-edge diffraction under Gaussian and Bessel illuminations have also been proposed in [21,22] and in [23], respecively.

The basic issue of sharp-edge diffraction theory is the transformation of the two-dimensional (2D) Kirchhoff integral, which implements Huygens’ superposition principle, into a contour (i.e., 1D) integral over the aperture boundary. The most known mathematical tool to achieve such a conversion is Green’s theorem. Gordon [24] and, independently, Asvestas [25,26], Forbes and Asatryan later [27], emphasized the importance of Poincaré’s elegant construction of vector potentials [28] as a practical tool for achieving surface-to-line conversion of Kirchhoff’s integral. Now, the paraxial limit of Kirchhoff’s integral is Fresnel’s integral. It would then be natural to ask whether it is possible again to invoke paraxial approximation from the beginning, in order for a general sharp-edge paraxial diffraction theory to be developed. Up to my knowledge, this does not seem to have yet been proposed.

The idea is simple: to employ the Poincaré vector potential construction directly into Fresnel’s integral, in order to convert it into a contour integral over the edge. While, in principle, such a conversion turns out to be always possible, its practical applicability depends on the capacity of analytically solving certain 1D integrals. In the present paper it is proved that such conversion is possible for an important subclass of the Laguerre-Gauss beam family. This would be enough to explore a virtually infinite variety of different scenarios. To give a single example, the near-field produced by the sharp-edge diffraction of vortex beams by triangular apertures is explored. Similar scenarios have already been analyzed in the past, but limitedly to the far-field zone.

If the conversion to a single contour integral were not analytically achievable, the proposed approach unavoidably leads to a double integral representation of the diffracted wavefield. However, differently from Fresnel’s integral, whose domain coincides with the aperture, the new representation turns out to always be defined onto a rectangular domain, a fact that greatly simplify its numerical computation. To give evidence of this, an iconic example will be illustrated: the focal wavefield distribution produced by a collimated laser beam impinging onto a water droplet lens whose boundary is forced to assume an equilateral triangular shape. The simulations are aimed at reproducing some of the beautiful experimental results obtained forty years ago by Berry, Nye, and Wright in a seminal paper which became part of the Catastrophe Optics manifesto [29]. All historical considerations aside, this example constitutes an important test to check the practical applicability of the proposed approach in rather extreme situations, with Fresnel numbers of the order of thousands, and where the 2D integral can be numerically evaluated by employing standard Montecarlo integration packages.

2. Theoretical analysis

2.1 Preliminaries: the Fresnel integral

Consider a scalar disturbance impinging onto a planar, opaque screen having an aperture $\mathcal {A}$ delimited by the boundary $\Gamma =\partial \mathcal {A}$. The screen is placed at the plane $z=0$ of a suitable cylindrical reference frame $(\boldsymbol {r};z)$. On denoting $\psi _0(\boldsymbol {r})$ the disturbance distribution at $z=0^{-}$, the field distribution $\psi$ at the observation point $P\equiv (\boldsymbol {r};z)$, with $z>0$ is given, within paraxial approximation and apart from an overall phase factor $\exp (\mathrm {i} k z)$, by the Fresnel integral

The whole paraxial scalar diffraction theory is based on Eq. (1), which can also be recast in terms of the new integration variable ${ \boldsymbol {R}}=\boldsymbol {\rho }-\boldsymbol {r}$ as follows:

 

Fig. 1. Describing the sharp-edge aperture from different reference frames.

Download Full Size | PDF

2.2 Reducing Fresnel’s integral to contour integrals: the Poincaré algorithm

Under very general hypotheses the Fresnel integral, written both in the forms (1) and (2), can be reduced, in principle, to a (1D) contour integral defined onto the boundary $\Gamma$. Hannay first noted that, for plane-wave illumination (i.e., $\psi _0=1$), the conversion of Fresnel’s integral in the form given by Eq. (2) can be achieved in a trivial way simply by expressing the integration variable ${ \boldsymbol {R}}$ through its polar coordinates (see Fig. 1), in such a way that [13,14,16]

For a typical incident disturbance $\psi _0$, the surface-to-line conversion could be achieved by interpreting the integrals into Eqs. (1) and (2) as fluxes of suitable divergence-free transverse vectorial fields, say ${ \boldsymbol {f}}(\boldsymbol {\rho })=f(\boldsymbol {\rho })\hat {\boldsymbol {k}}$ and ${ \boldsymbol {F}}({ \boldsymbol {R}})=F({ \boldsymbol {R}})\hat {\boldsymbol {k}}$, respectively, where

The same holds by formally letting $F\,\to \,f$, ${ \boldsymbol {R}}\,\to \,\boldsymbol {\rho }$, $A\,\to \,a$. This theorem follows from a more general theorem about vector potential representation of divergence-free vectorial fields in the three-dimensional Euclidean space. Readers are encouraged to go through [24–27] to appreciate the simplicity and the elegance of the proof, first due to Poincaré [24], which, as Forbes and Asatryan pointed out in their work, “deserves to be in all the handbooks and texts, but it is not” [27]. A notable exception is the beautiful textbook by Wilfred Kaplan [30].

Now, using $F({ \boldsymbol {R}})$ or $f(\boldsymbol {\rho })$ will allow different interpretation schemes for the paraxial field $\psi$ to be given. Before doing this, it is worth checking Eq. (3). To this end, it is sufficient to employ the potential vector given by Eq. (6), together with Stokes’ theorem, which gives at once

Then, on further introducing the infinitesimal angle $\mathrm {d}\varphi ={ \boldsymbol {R}}\,\times \,\mathrm {d}{ \boldsymbol {R}}\,\cdot \,\hat {\boldsymbol {k}}/R^{2}$ sketched in Fig. 2, Eq. (8) can be recast as

On substituting from Eq. (10) into Eq. (9), trivial algebra leads to Eq. (3).

The use of $\boldsymbol {\rho }$ instead of ${ \boldsymbol {R}}$ as integration variable, thus of ${ \boldsymbol {f}}$ in place of ${ \boldsymbol {F}}$ into Eq. (5), provides new and interesting results. To this aim, it is sufficient to recast the first of Eq. (4) through the identity $|\boldsymbol {\rho }-\boldsymbol {r}|^{2}=\rho ^{2}+r^{2}-2\boldsymbol {\rho }\cdot \boldsymbol {r}$, so that the following integral representation of $\psi$ is obtained:

 

Fig. 2. The geometrical meaning of $\mathrm {d}\varphi ={ \boldsymbol {R}}\,\times \,\mathrm {d}{ \boldsymbol {R}}\,\cdot \,\hat {\boldsymbol {k}}/R^{2}$ and of $\mathrm {d}\theta =\boldsymbol {\rho }\,\times \,\mathrm {d}\boldsymbol {\rho }\,\cdot \,\hat {\boldsymbol {k}}/\rho ^{2}$.

Download Full Size | PDF

Eqs. (11) and (12) contitute the main result of the present paper. Before continuing, it is worth recalling some important aspects of the present formulation. First of all, it allows severe numerical problems characteristic of the classical BDW theory to be naturally bypassed. In fact, as pointed out by Forbes and Asatryan [27], within the BDW theory both the geometrical and the BDW components of the total diffracted field are discontinuous along the geometrical shadow, and this necessarily implies the integrand of the integral over the rim to be singular, in order to compensate the discontinuity in the geometric field. But the reduction of Fresnel’s integral provided by Eqs. (11) and (12) does not include any discontinuous geometric field, thus allowing singularities to be avoided. If the initial field $\psi _0(\boldsymbol {r})$ does allow the analytical exact evaluation of the integral (12), then the diffracted field $\psi$ will be expressed by a single contour integral. If not, the propagated field $\psi$ will unavoidably be expressed through a double integral, whose numerical evaluation should expect to be easier than that of the original Fresnel’s integrals (1) and (2). However, if $\Gamma$ were circular (unit radius), Eqs. (11) and (12) would be equivalent to compute Fresnel’s integral (1) via polar coordinates $(\rho,\theta )$, being $\tau \equiv \rho$ . The computational novelty of Eqs. (11) and (12) can then be unveiled by exploring sharp-edge diffraction from noncircular apertures. In fact, differently from Fresnel, whose integration domain coincides with the aperture $\mathcal {A}$, Eqs. (11) and (12) imply that, regardless the aperture shape, the new double integral representation of $\psi$ will be always defined onto the Cartesian product between the $\tau$ integration interval $[0,1]$ and the (finite) integration interval related to the parametrization of the boundary $\Gamma$. As we shall see, this allows standard Montecarlo integration techniques on hypercubes to be efficiently employable.

3. Numerical results

3.1 Preliminaries

A couple of examples will now be illustrated which, for simplicity, involve the same aperture, namely the equilateral triangle $ABC$ shown in Fig. 3. The characteristic length will then be identified by the radius $\overline {OA}$ of the circumscribed circle. Accordingly, in the following it will be set $\overline {AB}=\sqrt 3$. Each side of the boundary $\Gamma$ will be parametrized according to the most natural choice. So, the parametrization of the side $AB$ reads (see Fig. 3)

 

Fig. 3. The equilateral triangular diffracting sharp-edge aperture.

Download Full Size | PDF

3.2 Paraxial diffraction of Gaussian vortex beams

The first example deals with the effect of a sharp-edge aperture on the propagation of light beams carrying on vortices of a given topological charge, say $m\ge 1$. A simple analytical model for the impinging beam is

Now, it can be proved that the whole family of integrals:

On using Eq. (15) into Eq. (11) and on recalling the above described triangle parametrization, the diffracted wavefield $\psi (\boldsymbol {r};U)$ can be written (apart from unessential amplitude and overall phase factors) as

The results of our numerical experiment are shown in Fig. 4. A collimated LG beam carrying on vortex with unitary topological charge impinges onto the triangular shape of Fig. 3. The spot-size of the incident beam will be set to $w_0=\pi \,a$, i.e., larger enough than the aperture size to simulate a plane wave carrying on a unitary charge vortex, as suggested in [32]. In Fig. 4, 2D maps of the modulus (a) and the phase (b) of the diffracted field $\psi (\boldsymbol {r};U)$, numerically evaluated via Eq. (18), are plotted for $U=5,\,10,\,20,\,50,\,100,\,200$ (note that the $U$ scale is logarithmic).

From the figure it is possible in particular to appreciate the evolution of the topological complexity of the 2D field distribution from the near ($U=200$, bottom) to the far ($U=5$, top) zone. The possibility of modeling near-field diffraction of vortex beams in such a simple way should be positively acknowledged as a further powerful investigation tool to explore light’s orbital angular momentum. In fact, since the pioneering work [33], experimental and theoretical investigations have mainly been focused on Fraunhofer diffraction, a considerably easier scenario to be numerically simulated with respect to Fresnel [32,34,35].

 

Fig. 4. 2D maps of the modulus (a) and the phase (b) of the diffracted field $\psi (\boldsymbol {r};U)$, evaluated via Eq. (18), are plotted for $U=5,\,10,\,20,\,50,\,100,\,200$.

Download Full Size | PDF

3.3 Natural focusing of light by water droplets

The second, nearly iconic, example of application of Eqs. (11) and (12) will now be illustrated. The following quotation well describes the experimental situation [29]:

The experiment was aimed at producing what is known to be an elliptic umbilic diffraction catastrophe. Diffraction catastrophes are mathematical bricks with which, at the end of seventies, John Nye and Michael Berry founded the so-called Catastrophe Optics: a new, modern theoretical framework aimed at studying the so-called “natural focusing” of light [19,20]. CO’s description of focused wavefields is built up starting from light skeletons of bright caustics which are decorated, at the wavelength scale, by characteristic diffraction patterns organized according to a precise hierarchy [19]. The elliptic umbilic is one of them.

Some of the experimental results presented in [29] will now be reproduced thanks to a straightforward implementation of the general method here developed. To this aim, the radius of the circumscribed circle to the triangle $a=\ell /\sqrt 3$ is again introduced, while the refractive index of the droplet will be set to $n \simeq 4/3$. In [29] a simple and effective physical model of the droplet profile, based on the combined action of surface tension and gravity, was developed (see also [36]). In particular, on denoting $h(\boldsymbol {r})$ the height of the droplet upper surface below the plane $z = 0$, we have

A contour level map of the normalized profile $h/H$ is sketched in Fig. 5 to give readers an idea of the 3D profile, as well as the $2\pi /3$ symmetry. Once the droplet is illuminated by the laser beam, the diffracted field can be obtained (up to unessential overall phase factors) by evaluating the Fresnel integral into Eq. (1) with $\psi _0$ given by

 

Fig. 5. Contour level map of the droplet normalized profile $h(\boldsymbol {r})/H$ given in Eq. (19).

Download Full Size | PDF

The nonregular shape of $\mathcal {A}$, together with the nonsmall values of $U$ (of the order of thousands), would make the evaluation of the 2D Fresnel integral (Eq. (1)) a challenging numerical task. In [29] such technical difficulties were partially circumvented by suitably extending the integration domain $\mathcal {A}$ to cover the whole Euclidean plane $\mathbb {R}^{2}$, thus neglecting the edge wave contribution. In other words, $\psi (\boldsymbol {r};U)$ was approximated through an elliptic umbilic diffraction catastrophe, whose evaluation can be achieved for instance by using suitable asymptotic techniques. In the following, Eqs. (11) and (12) will be employed to reproduce some of the experimental results shown in Fig. 2 of [29], without any approximations but the paraxial one. To help the comparison, Fresnel’s number $U$ will be recast as follows:

Method $\to$ "AdaptiveMonteCarlo"

"MaxPoints" $\to 10^{5}$

In Fig. 6, 2D maps of the transverse intensity distribution $|\psi (\boldsymbol {r};U)|^{2}$ are shown for some values of the dimensionless absicssa $\zeta$. They are aimed at reproducing the experimental pictures reported in Fig. 2 of Ref. [29], precisely Figs. 2(a), 2(b), 2(c), 2(d), 2(g), 2(h), and 2(i), corresponding to $\zeta =0$ (a), 1 (b), 2 (c), 3 (d), 4 (g), 4.9 (h), and 5.81 (i), respectively. Each slice of Fig. 6 contains a $200\,\times \,200$ matrix of intensity values. In Fig. 7 some blowups of the upper slice of Fig. 6, which corresponds to $\zeta =5.81$, are also shown. The figure appears to be slightly noisy, due to Montecarlo. Nevertheless, the visual agreement between Figs. 6 and 7 with the experimental results presented in Fig. 2 of [29] is excellent.

 

Fig. 6. 2D maps of the optical intensity distribution of the field focused by the droplet lens of Fig. 5, evaluated via Eq. (18), are plotted for $\zeta =0,\,1,\,2,\,3,\,4,\,4.9,\,5.81$.

Download Full Size | PDF

 

Fig. 7. Blowups of the upper slice of Fig. 6, corresponding to $\zeta =5.81$.

Download Full Size | PDF

4. Conclusions

Fresnel’s diffraction theory represents a milestone of classical optics for more than two centuries. However, the numerical evaluation of 2D diffraction integrals remains a challenging task, due to the highly oscillatory behavior of the integrands and the shape of the integration domains. For plane and/or spherical waves, a reduction of Fresnel’s integral to 1D, singular-free contour phase integrals is always guaranteed as far as sharp-edge diffraction is concerned. Such a formulation, which has been developed during the last decade in the light of catastrophe optics, revealed to be an unorthodox, very interesting point of view from which diffraction phenomena can be explored.

In the present paper a further step toward a general paraxial sharp-edge diffraction theory dealing with, in principle, arbitrary wavefields impinging onto arbitrarily shaped planar apertures, has been proposed. The algorithmic engine is the beautiful mathematical theorem which allows a bonafide vector potential representation of a given solenoidal vector field to be straightforwardly obtained. Several authors in the past attempted to emphasize the importance of such theorem in order to give classical diffraction theory a new, modern dress. Unfortunately, the intrinsic mathematical complexity involved into Kirchhoff diffraction theory has attenuated the interest of researchers on this fascinating subject in the past four decades. What has been shown here is that Poincaré’s vector potential construction algorithm, when is employed in a mathematically more friendly context, and paraxial diffraction theory certainly is, leads to a “minimal” formulation which could be potentially very attractive in computational optics. In particular, by using Poincaré’s algorithm, the Fresnel integral has been converted into a contour integral over the aperture rim. In this way, it has been proved that sharp-edge diffraction of a whole subclass of Laguerre-Gauss beams carrying on vortices of arbitrarily high topological charges can numerically be dealt only with 1D integrals.

But also when the desired “2D$\to$1D” conversion cannot be analytically achieved, the present approach is always able to convert the Fresnel integral, defined onto arbitrarily shaped 2D domains, into a (2D) integral defined onto a simple square, which is particularly suitable for Montecarlo integration. This has been pointed out in the second numerical experiment, carried out on an iconic example of natural focusing of light by water droplets. An excellent agreement with the experimental results obtained by Berry, Nye, and Wright was achieved. All numerical simulations presented in the paper have deliberately been carried out with a “low profile” approach. The computing machine employed to generate all figures was a commercial laptop equipped with a 3 GHz Intel Core i7 processor and 16 GB RAM. Moreover, all numerical integrations (both 1D and 2D) presented have been performed via standard native Mathematica routines, which are nowadays available also (and especially) to nonspecialists. This was done also to convince our readers about the feasibility, the implementation easiness, and the numerical effectiveness of the proposed method also when “extreme” scenarios are dealt with. We do hope that researchers interested in computational optics could find here useful material for developing new, more effective numerical methods.

Future studies are in progress. Among these, the application to cascaded diffraction in optical systems is one of the most relevant, as witnessed by the recent literature [37–39]. When diffraction by a sequence of sharp-edge apertures has to be tackled, the approach here developed is expected to be highly promising. In particular, the emerging wavefield would naturally be represented in terms of multiple integrals defined onto hypercubes, a perfect scenario for effectively employing Montecarlo-based computational techniques.

Light scattering from “large” tridimensional objects is another topic which would be worth exploring with the above computational tools. In particular, the basic features of the scattered wavefields could, in principle, be grasped by replacing the scatterers by “equivalent” planar apertures. Raman and Krishnan first experimentally explored this topic [40], which has continued to receive attention, also due to its important astronomical implications [41–44].

Finally, the proposed approach could also be helpful in numerically exploring the role played by sharp-edge diffraction in the study of the fractal nature of the light, pioneered in [45] and presently a topic of central interest in theoretical and applied optics (see for instance the review in [46]).

“There is pleasure in recognizing old things from a new viewpoint,” Richard Feynman loved to say. Such a pleasure is sometimes also accompanied by the discover of unexpected, useful, new results. During the last decade, the change of perspective offered by Catastrophe Optics in describing paraxial sharp-edge diffraction provided new light on old, nearly forgotten experiments, as well as new, unorthodox interpretation schemes of diffraction problems by now considered obsolete. We hope what has here been presented could be helpful in continuing to explore new, still unveiled aspects of classical diffraction theory.