Introduction. Fourier optics is a theoretical and computational framework for understanding and simulating the propagation of light through free space and dielectric optical elements, such as lenses, zone plates, gratings, and apertures [1]. Its uses and applications are many and varied, underpinning the analysis of classical imaging systems—from microscopes to telescopes [1], the realization of super-resolution imaging via Fourier ptychography [2], the implementation of digital holography for phase-sensitive imaging and retrieval [3], the computationally tractable production of computer-generated holograms (CGH) for 3D display applications [4], and the design of optical computing devices for low-energy, high-throughput, and signal-integrated inferential computer vision systems [5]. In addition to providing a satisfactory model of electromagnetic wave propagation in many applications, Fourier optics provides highly computationally efficient formulations of diffraction by leveraging fast Fourier transforms (FFTs), making it ideal for training large physics-aware deep-learning models [6], incorporating into intensive gradient descent procedures [7], or performing real-time optimizations on conventional computer hardware [8]. In other cases, the Fourier optical framework demonstrates the potential to extract useful computational operations from light propagation itself, since in specific optical layouts (such as that produced by a single lens), the action of propagation through the optical system results, to a reasonable accuracy, in a Fourier relationship between input and output planes [9]. Alternatively, the ready availability of an optical Fourier transform can be exploited to perform Fourier domain filtering using, for example, spatial light modulators [10] or metasurfaces [11] for the hardware acceleration of computationally expensive convolutions.

In recent work in the field of photonic computing, some of the familiar models from Fourier optics in 3D free space have been adapted to the setting of slab (or planar) waveguides [12–20], and surface plasmon polaritons [21], providing diffractive processing resources in chip-sized form factors. In Refs. [15], an angular spectrum model is used to optimize the settings of a series of 1D metasurfaces inserted between diffractive regions, implementing the weights of consecutive layers of a neural network, by analogy with work in 3D free space [22]. In Refs. [14,17], a Fourier domain Fresnel approximation is used to similar effect. By contrast, scalar diffraction can be described directly in terms of Green’s functions, loosely representing the propagation of point sources or Huygens’ wavelets. This style of analysis can be used to show that a star coupler architecture functions as a Fourier transforming device [16,23], a tool that has been used to implement the convolution steps of neural networks [13,16], miniaturizing free-space designs that rely on lens optics [10].

While the downconversion of 3D formulas into 2D equivalents has had some practical success so far in the literature, there are certain gaps in the theoretical treatments presently available. We address these gaps in this Letter, providing theoretical parity with the key results of free-space optics. In all the derivations, care is taken around the manipulation of Green’s functions. In 3D, we are familiar with the field of a point source resembling a complex exponential. Not so in 2D—we begin with a Hankel function to represent the expanding point source in planar contexts, and it is only after asymptotic approximations are applied that we recover a complex exponential to high accuracy for sufficiently large argument [24].

The following results are obtained. First, for completeness, we present a simple proof of the (2D) Weber integral formula in a manner that generalizes to 3D (the Kirchhoff–Helmholtz integral formula; see Ref. [25] for typical approaches). Next we derive equivalents to the Rayleigh–Sommerfeld (RS) diffraction formulas, from which a Fresnel approximation is derived, along with an appropriate radiation condition (for corresponding results in 3D, see Refs. [1,24]). We then present a 2D angular spectrum formulation, show that it is equivalent to a direct formulation of diffraction, one of the RS formulations (for corresponding results in 3D, see Ref. [26]), and, further, that it yields a Fourier-domain Fresnel approximation. Finally, we observe the equivalence of the Fresnel approximations in Fourier and direct forms (for corresponding results in 3D, see Ref. [1]).

Results. Proof of Weber–Kirchhoff–Helmholtz integral formula. Suppose $G$ is a Green’s function satisfying an inhomogeneous Helmholtz wave equation,

We introduce an auxiliary surface $C_{\epsilon }$, a circle of radius $\epsilon$ around $\mathbf {r}$ that enables us to avoid applying the divergence theorem at a point of discontinuity. Radius $\epsilon$ can be made arbitrarily small so that $C_{\epsilon }$ fits entirely in $S$, as in Fig. 1.

 

Fig. 1. Diagram corresponding to the variables in the Weber integral formula and its modification to incorporate $C_{\epsilon }$.

Download Full Size | PDF

Applying the divergence theorem, then adding zero, we have

RS solutions. If we can orchestrate a Green’s function where $G(\mathbf {r} - \mathbf {r}')$ vanishes for $\mathbf {r}$ on the boundary $\partial S$, we are left with a 2D equivalent of the first RS diffraction formula:

Equivalently, the second RS diffraction formula is formed by causing $\partial G / \partial n$ to vanish on the boundary $\partial S$:

Constructing a Green’s function with the appropriate vanishing conditions is simple enough if $S$ is a half-space, say $z\ge 0$. Then from $\mathbf {r} = (x,z)$, one can construct its reflection in $\partial S$, $\boldsymbol {\tilde {r}}=(x,-z)$. It is easy to show that $G^{-}_{2D}$ and $G^{+}_{2D}$ satisfy the conditions for the first and second RS formulas respectively if defined as follows:

We have glossed over an important detail, however. The Helmholtz equation is elliptic and hence requires Dirichlet (or Neumann) conditions over a closed boundary. We attain the half-space $z\ge 0$ as the limit as $R\to \infty$ of a parameterized region $S_{R} = S \cap C_{R}$, where $C_{R}$ is a circle of radius $R$ centered at $\mathbf {r}$. The curved part of the boundary is given by $\partial C_{R} \cap S$, while the straight part is given by $\partial S \cap C_{R}$. One such $S_{R}$, with relevant boundary portions labeled, is shown in Fig. 2. By imposing the following condition (a “radiation condition”) on $U$ on the boundary portion $\partial C_{R} \cap S$ as $R\to \infty$, we ensure that the contribution of the curved boundary to the integral tends to zero and we are free to make use of the limiting region $S$ with straight boundary at $z=0$:

To see that the radiation condition is sufficient to nullify the contribution of the curved boundary to the Weber integral formula, we label this contribution as $\mathcal {I}$:

Equivalence with angular spectrum formulation. Assuming $U$ satisfies the radiation condition, we can make use of the first RS formula in the case where $S$ is the half-space $z \ge 0$. Let us suppose, as is often the case in applications, we wish to characterize $U$ at values of $\mathbf {r}$ in a plane with some fixed value $z$, propagating the wave forward from $\partial S$. Modifying notation to reflect the fact that $\partial S$ is constant in $z'$, we obtain

 

Fig. 2. Construction of $S_R$ to accommodate the radiation condition.

Download Full Size | PDF

An argument appears in Ref. [27] (pp. 112–115, with reference to pp. 105–106 and 103) to compute the Fourier transform of $G_{2D}$ as defined in Eq. (16). But the argument depends only on the Green’s function defining equation and the decay of the function and its derivatives at infinity. $G_{2D}^{-}$ satisfies both of these requirements and so we can carry over the results (the derivation is explicitly adapted in Supplement 1). We conclude that

Fresnel approximation. To form an approximation to the integral in Eq. (18), we first observe that the normal derivative of $G_{2D}^{-}$ matches that of $G_{2D}$ on the boundary $\partial S$, but for a factor of $2$:

If, further, we make a paraxial assumption, that $|x-x'|$ is small compared with $|z-z'|$, we arrive at a Fresnel approximation for 2D diffraction:

Once again, we have a convolution in the $x$ dimension. Separating out the propagation kernel, we obtain

The Fresnel approximation can equally be expressed in the spatial frequency domain by way of the angular spectrum formulation in Eq. (19), whereupon a truncated binomial series representation of $(1-k^{2}/\alpha ^{2})$ yields

Discussion. This Letter has provided some of the missing analytical pieces in the use of time-independent modeling of diffraction in planar contexts. It closes the gap between prior works that may seem to take disparate approaches, but may now be seen to be equivalent or closely related. For example, it is seen that the angular spectrum formulation used in layered diffractive photonic neural network modeling and optimization [15] can be thought of as a Fourier-domain solution to a Helmholtz equation. The Fourier-domain Fresnel modeling approach adopted in other layered architecture optimizations is an approximate form of this same solution [14,17]; that this formulation is an approximate solution, increasingly accurate for larger distances, may account, in part, for the fact that the modeling was found to require longer diffraction distances to ensure device accuracy (order of $250~\mu$m with a Fresnel model [17] versus distances as short as $20~\mu$m in the angular spectrum regime [15]). In contrast to Fourier-domain Fresnel approximations, derivations of the Fourier transforming properties of star couplers [13,16] and other curved wavefront designs [21] make use of a direct Fresnel model of diffraction, relying on the Weber integral formula, as we do.

The assumption that the wave function in question is time-separable as $u(x,z, t) = U(x,z)\exp (i\omega t)$, however, while a reasonable assumption in systems of low operating frequencies, fails to reflect the physical reality of 2D wave propagation, which is diffusive in nature, with wavefronts spreading as they propagate [24,28]. Indeed, for systems with high operating frequencies (in the range of tens or hundreds of gigahertz for standard slab waveguide configurations), a time-based error emerges in addition to other terms dependent on the spatial variables [29].

In future work, we wish to examine further distinctions between time-aware and time-independent analytical methods. The radiation condition derived in this work is formally necessary to find solutions for open regions, since the Helmholtz equation is an elliptic partial differential equation, which requires boundary conditions on a closed boundary [24]. The radiation condition can be thought of, however, as accounting for the fact that the wave phenomena in question are temporal, and, due to finite propagation speeds, will only have propagated over a finite expanse of space if measured at a given instant in time. If the full (time-dependent) wave equation, a hyperbolic partial differential equation, is employed, the requirement of a closed boundary is no longer required, and the artificial introduction of a radiation condition is no longer needed, yielding fewer conditions on the wave function $u$. This suggests that a time-aware approach is a more physically meaningful methodology, and, indeed, we note that these considerations apply in the analysis of 3D propagation just as they do in planar contexts.