Paraxial Sharp-Edge Diffraction: A General Computational Approach

A general reformulation of classical sharp-edge diffraction theory is proposed within paraxial approximation. The, not so much known, Poincar\'e vector potential construction is employed directly inside Fresnel's 2D integral in order for it to be converted into a single 1D contour integral over the aperture boundary. Differently from the recently developed paraxial revisitation of BDW's theory, such approach can be applied to arbitrary wavefield distributions impinging onto arbitrarily shaped sharp-edge planar apertures. A couple of interesting examples of application of the proposed method is presented.


I. 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, sharpedge 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 the BDW/Miyamoto/Wolf theory turned out to not to be as easy as it could have been expected, even for apparently simple incoming disturbances, like for instance Gaussian beams. 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 [8,9], based on some results published in [10,11]. This paraxial revisitation of BDW theory soon revealed its predictive potential [12,13], especially once placed within the context of the so-called Catastrophe Optics [14,15]. Later on, further generalizations aimed at dealing with sharp-edge diffraction under Gaussian and Bessel illuminations have also been proposed in [16,17] and in [18], respecively.
The basic issue of sharp-edge diffraction theory is the transformation of the two-dimensional (2D) Kirchhoff in-tegral, 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 [19] and, indipendently, Asvestas [20,21], Forbes and Asatryan later [22], emphasized the importance of Poincaré's elegant construction of vector potentials [23] 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 sharpedge 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 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 equilater 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 [24]. All historical considerations aside, arXiv:2204.05685v1 [physics.optics] 12 Apr 2022 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.

II. THEORETICAL ANALYSIS
Consider a scalar disturbance impinging onto a planar, opaque screen having an aperture A delimited by the boundary Γ = ∂A. The screen is placed at the plane z = 0 of a suitable cylindrical reference frame (r; z). On denoting ψ 0 (r) the disturbance distribution at z = 0 − , the field distribution ψ at the observation point P ≡ (r; z), with z > 0 is given, within paraxial approximation and apart from an overall phase factor exp(ikz), by the Fresnel integral (1) where, in place of z, the Fresnel number U = ka 2 /z has been introduced. The symbol a denotes a characteristic length of the aperture size (for instance, if A were circular, a would coincide with its radius). In this way, only dimensionless quantities will be involved in the following.
The whole paraxial scalar diffraction theory is based on Eq. (1), which can also be recast in terms of the new integration variable R = ρ + r as follows: and which corresponds to observe the aperture A from a reference frame centred at the observation point P . In Fig. 1 the different viewpoints are shown. The z-axis will be supposed to be the mean propagation direction of the incident beam ψ 0 (r). Under very general hypotheses Fresnel's integral, written both in the forms (1) and (2), can be reduced, in principle, to a (1D) contour integral defined onto the boundary Γ. Hannay first noted that, for plane-wave illumination (i.e., ψ 0 = 1), the conversion of the Fresnel integral into Eq. (2) can be achieved in a trivial way simply by expressing the integration variable R through its polar coordinates (see Fig. 1), in such a way that [8,9,11] For a typical incident disturbance ψ 0 , the surface-toline conversion could be achieved by interpreting the integrals into Eqs. (1) and (2) as fluxes of suitable divergencefree transverse vectorial fields, say f (ρ) = f (ρ)k and and wherek denotes the unit vector of the z-axis. In this way Eqs. (1) and (2) can formally be recast as follows: where, in order to alleviate notation complexity, the explicit dependence of ψ on U and r will be tacitly assumed henceforth. The following mathematical theorem will then play a major role in the subsequent analysis: Consider the transverse vectorial field F (R) = F (R)k. Then, under very general hypotheses about the scalar field F (R), it is and The same holds by formally letting F → f , R → ρ, A → a. This theorem follows from a more general theorem about vector potential representation of divergencefree vectorial fields in the three-dimensional Euclidean space. Readers are encouraged to go through [19][20][21][22] to appreciate the simplicity and the elegance of the proof, first due to Poincaré [19], which, as Forbes and Asatryan pointed out in their work, "deserves to be in all the handbooks and texts, but it is not" [22]. A notable exception is the beautiful textbook by Wilfred Kaplan [25]. Now, using F (R) or f (ρ) will allow different interpretation schemes for the paraxial field ψ 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 dϕ = R × dR ·k/R 2 sketched in Fig. 2, Eq. (8) can be recast as while, on setting ψ 0 = 1 into the second row of Eq. (4), Eq. (7) gives On substituting from Eq. (10) into Eq. (9), trivial algebra leads to Eq. (3).
The use of ρ instead of R as integration variable, thus of f in place of 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 |ρ − r| 2 = ρ 2 + r 2 − 2ρ · r, so that the following integral representation of ψ is obtained: where now (12) Before continuing, it is worth recalling an important aspect of the present formulation. If the initial field ψ 0 (r) does allow the analytical exact evaluation of the integral (12), then the diffracted field ψ will be expressed by a single contour integral. If not, the propagated field ψ 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 Γ were circular (unit radius), Eqs. (11) and (12) would be equivalent to compute Fresnel's integral (1) via polar coordinates (ρ, θ), being τ ≡ ρ . 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 A, Eqs. (11) and (12) imply that, regardless the aperture shape, the new double integral representation of ψ will be always defined onto the Cartesian product between the τ integration interval [0, 1] and the (finite) integration interval related to the parametrization of the boundary Γ. As we shall see, this allows standard Montecarlo integration techniques on hypercubes to be efficiently employable.
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 OA of the circumscribed circle. Accordingly, in the following it will be set AB = √ 3. Each side of the boundary Γ will be parametrized according to the most natural choice. So, the parametrization of the side AB reads (see Fig. 3) and similarly for the other two sides. In this way, it is trivial to prove that ρ × dρ ·k = √ 3/2 dt inside Eq. (11), what considerably simplifies the evaluation of the contour integral Γ .
The first example deals with the effect of a sharp-edge aperture on the propagation of light beams carrying on vortices of given topological charge, say m ≥ 1. A simple analytical model for the impinging beam is Now, it can be proved that the whole family of integrals: can be evaluated exactly, for any complex L and real M , through the following notable recursive rule: which is one of the most relevant results of the present paper. Equation (17) provides an important generalization of the results recently found in [16,17] to sharp-edge diffraction of Gaussian beams [26].
On using Eq. (15) into Eq. (11) and on recalling the above described triangle parametrization, the diffracted wavefield ψ(r; U ) can be written (apart from unessential amplitude and overall phase factors) as where ρ j (t) = (ξ j (t), η j (t)) defines the parametrization of the jth triangle side.
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 = π a, i.e., larger enough than the aperture size to simulate a plane wave carrying on a unitary charge vortex, as suggested in [27]. In Fig. 4, 2D maps of the modulus (a) and the phase (b) of the diffracted field ψ(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 [28], experimental and theoretical investigations have mainly been focused on Fraunhofer diffraction, a considerably easier scenario to be numerically simulated with respect to Fresnel [27,29,30].
The second, nearly iconic, example of application of Eqs. (11) and (12) will now be illustrated. The following quotation well describes the experimental situation [24]: A hole whose shape approximated an equilateral triangle of side =2.6 mm was cut in adhesive tape stuck on to the horizontal surface of a glass microscope slide. A water droplet was allows to fall on to the slide, where it formed a thin lens. This lens was illuminated from below with a parallel beam of laser light (wavelength λ= 633 nm) broadened so as to fill the aperture of the lens. After refraction the focused light formed an elliptic umbilic diffraction catastrophe a few centimetres above the lens.
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 [14,15]. 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 [14]. The elliptic umbilic is one of them.
Some of the experimental results presented in [24] 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 = / √ 3 is again introduced, while the refractive index of the droplet will be set to n 4/3. In [24] a simple and effective physical model of the droplet profile, based on the combined action of surface tension and gravity, was developed (see also [31]). In particular, on denoting h(r) the height of the droplet upper surface below the plane z = 0, we have  A 3D visualization of the normalized profile h/H is sketched in Fig. 5. 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 ψ 0 given by ψ 0 (r) = exp[−iα (3 x 2 + 3 y 2 − 2 x 3 + 6 xy 2 )] , (20) where now α = k (n − 1) H 330.868 . . . The values of U will be chosen according to the prescriptions described in [24]. In particular, the droplet lens focus is located at a distance from the aperture plane equal to 2 /(18H(n − 1)) 11.2 mm, corresponding to U = 6α 2 × 10 3 . The subsequent simulations will then be carried out in the neighborhood of such a value.
The nonregular shape of A, together with the nonsmall values of U (of the order of thousands), would make the evaluation of the 2D Fresnel integral (1) a challenging numerical task. In [24] such technical difficulties were partially circumvented by suitably extending the integration domain A to cover the whole Euclidean plane R 2 , thus neglecting the edge wave contribution. In other words, ψ(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 [24], without any approximations but the paraxial one. To help the comparison, Fresnel's number U will be recast as follows: where the symbol ζ denotes a dimensionless, normalized abscissa whose origin is located at the above defined fo-cal plane. In particular, the definition (21) has been chosen in order for ζ to coincide with the parameter employed in [24] to individuate the position of the observation plane during their experiment. On substituting from Eq. (20) into Eq. (12) and then into Eq. (11), and on taking Eq. (21) into account, the diffracted wavefield ψ(r; U ) can formally be written (apart from unessential amplitude and overall phase factors) as follows:  Fig. 6, which corresponds to ζ = 5.81, are also shown. The figure appears to be slightly noisy, due to Montecarlo.

III. CONCLUSIONS
Fresnel's diffraction theory represents a milestone of classical optics since 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 sharpedge 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. By using Poincaré vector potential construction, Fresnel's integral has been converted into a contour integral over the aperture rim. In this way, it has been proved that sharpedge 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. When the analytical conversion to a single contour integral is no longer possible, a new double integral representation of the diffracted wavefield, suitable for Montecarlo integration, has been derived. It was tested on an iconic example of natural focusing of light by water droplets, with extremely promising results.
All numerical simulations presented in the paper have deliberately been carried out with a "low profile" ap-proach. 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. This should convince our readers about the feasibility, the implementation easiness, and the numerical effectiveness of the proposed method also when "extreme" scenarios are dealt with.
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 [32][33][34]. When diffraction by a sequence of sharpedge 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 Montecarlobased 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 [35], which has continued to receive attention, also due to its important astronomical implications [36][37][38][39].
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 [40] and presently a topic of central interest in theoretical and applied optics (see for instance the review in [41]).
"There is pleasure in recognizing old things from a new viewpoint," Richard Feynman loved to say. 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, unexpected, and still unveiled aspects of classical diffraction theory.