Double-slit interferometry as a lossy beam splitter

To represent diffraction as a transformation in a Hilbert space, we use solutions of the Helmholtz equation [20, 21, 29]. The Helmholtz equation is a self-adjoint linear partial differential equation. Therefore, its solutions or fields have a vector in a Hilbert space associated with them. Moreover, the projections of the 3D fields on 2D planes, say the xy plane also form a Hilbert space. It should be noted that there is a difference between 2D fields and the projections of 3D fields on 2D planes. The 2D projections are referred to as slices of the 3D fields.

Diffraction of light is understood by solving the Maxwell's equations [20, 21, 29], specifically the wave equation, with appropriate boundary conditions. Generally, the time-dependence of solutions (or fields) is considered harmonic, i.e. of the form e^iωt, where ω is the angular frequency. Consequently, the wave equation reduces to the time-independent Helmholtz equation [20, 21, 29]. For a source $\rho ({\boldsymbol{r}})$ in a volume ${ \mathcal V }$ enclosed by a surface $\partial { \mathcal V }$ on which the boundary condition is specified, the most general solution of the Helmholtz equation is

$$\begin{eqnarray}&&E({\boldsymbol{r}})={\iiint }_{{ \mathcal V }}{{\rm{d}}}^{3}{\boldsymbol{r}}^{\prime} \,G({\boldsymbol{r}},{\boldsymbol{r}}^{\prime} )\rho ({\boldsymbol{r}}^{\prime} )+{\iint }_{\partial { \mathcal V }}{{\rm{d}}}^{2}{\boldsymbol{r}}^{\prime} \hat{n}({\boldsymbol{r}}^{\prime} )\cdot \left(E({\boldsymbol{r}}^{\prime} ){\rm{\nabla }}^{\prime} G({\boldsymbol{r}},{\boldsymbol{r}}^{\prime} )-G({\boldsymbol{r}},{\boldsymbol{r}}^{\prime} ){\rm{\nabla }}^{\prime} E({\boldsymbol{r}}^{\prime} ,),)\right),\end{eqnarray} \tag{ 1 }$$

where $\hat{n}({\boldsymbol{r}}^{\prime} )$ is unit normal to the surface and $G({\boldsymbol{r}},{\boldsymbol{r}}^{\prime} )$ is the Green's function for the Helmholtz equation.

We apply the Fraunhofer approximation to equation (1) (see appendix A for details) to find solutions of the Helmholtz equation and project them onto the xy plane by fixing z. Using the surface term in equation (1) we find a slice-to-slice map (see appendix B). To get a matrix representation of the slice-to-slice map, we represent each slice as a column vector in a suitable basis on that slice. One set of orthonormal vectors that span the Hilbert space on a slice, can be found by finding the eigensolutions of the homogeneous Helmholtz equation, with the appropriate boundary conditions [21], and using them as modes. For example, the eigensolutions are standing sinusoidal waves if the boundary condition is, say reflective. However, these modes are not localized and thus unsuitable for finite number of detectors with a given size. Therefore, we choose a basis of 2D functions with compact support, that spans the Hilbert space on a slice.

We use Haar wavelets and scaling function [30, 31] to construct orthonormal basis for a slice (projections of fields on a plane, see section 2.1). Compact support of the wavelets makes them suitable modes for detectors that have finite size. The orthogonality of the wavelets ensures that there is no overlap between measurements by two detectors. For a detector with square-shaped window, the 2D Haar wavelets [32] are chosen.

Wavelets are square-integrable functions with compact support over a finite interval. The simplest example is the Haar wavelet [30–32], which is defined by its wavelet function ψ (or mother function) and a scaling function ϕ (or father function)

$$\begin{eqnarray}&&\psi (x):= \sqcap \left(2\left(x-\displaystyle \frac{1}{4}\right)\right)-\sqcap \left(2\left(x-\displaystyle \frac{3}{4}\right)\right),\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\phi (x):= \sqcap \left(x-1/2\right),\end{eqnarray} \tag{ 3 }$$

respectively, where $\sqcap$ is the box function. These functions are dilated and translated to create other Haar wavelets and scaling functions

$$\begin{eqnarray}&&{\psi }_{m,n}(x):= \displaystyle \frac{1}{\sqrt{{2}^{-m}}}\psi \left(\displaystyle \frac{x}{{2}^{-m}}-n\right),\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{\phi }_{j,k}(x):= \displaystyle \frac{1}{\sqrt{{2}^{-j}}}\phi \left(\displaystyle \frac{x}{{2}^{-j}}-k\right),\end{eqnarray} \tag{ 5 }$$

respectively, where m and j are dilation parameters whereas n and k are translation parameters, and all take integer values. Examples of dilated and translated wavelet and scaling functions are shown in figures 1 and 2 respectively.

Zoom In Zoom Out Reset image size

The Haar wavelets are orthogonal with respect to the inner product

$$\begin{eqnarray}&&\left\langle {\psi }_{m,n},{\psi }_{m^{\prime} ,n^{\prime} }\right\rangle := {\int }_{-\infty }^{\infty }{\rm{d}}x\,{\psi }_{m,n}^{* }(x){\psi }_{m^{\prime} ,n^{\prime} }(x)={\delta }_{m,m^{\prime} }{\delta }_{n,n^{\prime} },\end{eqnarray} \tag{ 6 }$$

and the Haar scaling functions are also orthogonal, i.e.

$$\begin{eqnarray}&&\left\langle {\phi }_{{j}_{0},k},{\phi }_{{j}_{0},k^{\prime} }\right\rangle ={\delta }_{k,k^{\prime} },\end{eqnarray} \tag{ 7 }$$

at a particular scale, say, j = j₀ (see figure 2). In signal processing, any signal f(x) can be decomposed into Haar wavelets and scaling function of a particular scale j₀ [30–32] as

$$\begin{eqnarray}&&f(x)=\displaystyle \sum _{k=-\infty }^{\infty }{a}_{k}{\phi }_{{j}_{0},k}(x)+\displaystyle \sum _{m\geqslant {j}_{0}}\displaystyle \sum _{n=-\infty }^{\infty }{b}_{m,n}{\psi }_{m,n}(x),\end{eqnarray} \tag{ 8 }$$

where the scaling function is equivalent to a low-pass filter and the wavelet functions are equivalent to band-pass filters.

Zoom In Zoom Out Reset image size

By constructing a basis using Haar wavelets on each slice, we represent the projected field as a column vector in that basis, and express the slice-to-slice map as a transfer matrix between two slices. To validate our formalism we use this approach to show that a double-slit system is effectively a lossy beam splitter and verify it by studying the correlation of the outputs of the double-slit setup and compare it with that of a cubic beam splitter.

Wavelets are used as a set of orthogonal functions in a similar manner as a Fourier basis [32]. The advantage of Haar wavelets is that, unlike the Fourier basis functions which are localized only in frequency but are non-local in space, Haar wavelets are localized in both frequency and space. This makes Haar wavelets suitable postselected modes for detectors with rectangular-shaped windows.

A beam-splitter is a ubiquitous two-input-two-output component in interferometry. In optics, a beam splitter is commonly in the form of a glass cube, half-silvered mirror or fiber based, which have two input and two output modes corresponding to each of their ports. In a 50:50 cubic beam splitter, for example, the modes are the k-vectors corresponding to the plane wave entering each of its ports, forming a basis to represent the inputs and outputs as 2D column vectors in a Hilbert space. In such a representation, the beam splitter transformation has a 2D transfer matrix representation [5, 6, 33], denoted here by

$$\begin{eqnarray}{U}_{\mathrm{BS}}=\displaystyle \frac{1}{\sqrt{2}}\left(\begin{array}{cc}1 & {\rm{i}}\\ {\rm{i}} & 1\end{array}\right),\end{eqnarray} \tag{ 9 }$$

where each row corresponds to the superposition of the two input modes to form the outputs, and complex elements of the matrix denote the phase-shift introduced in each input.

In general, if the source of light does not emit in a single mode (say, a divergent beam), the vector representation of the inputs and outputs can be infinite-dimensional, yielding an infinite-dimensional transfer matrix of the beam splitter. In such a case, two suitable input and two output modes can be post-selected to reduce the infinite-dimensional transfer matrix to a post-selected 2 × 2 transfer matrix as in equation (9).

A consequence of such a transformation is that the outputs of the beam splitter are correlated, as discussed in the next subsection. The correlation of the outputs, as a function of a parameter that distinguishes the inputs, is a signature of a beam splitter. We use this signature to verify the claim that a double-slit setup is effectively a beam splitter.

In semi-classical theory of photo-detection [34, 35, 29], the probability of coincident photo-detections is proportional to the intensity–intensity cross-correlation of light in the post-selected modes, falling on the detectors. Such a correlation of the outputs of a 50:50 beam splitter plotted as a function of some distinguishability parameter shows a dip [36–38] for identical pulses (or photon states in the quantum regime [39]) at the input ports. A parameter, say, time-delay between the input pulses, distinguishes the otherwise identical input pulses. The correlation depends on the shape of the input pulses and the fluctuations in the light field. Specifically, if the fluctuation is uniform, the correlation shows a dip of 50% as the distinguishability parameter (like time-delay) approaches zero. A brief discussion of this concept is in appendix C and a detailed analysis of this phenomenon is presented in [40].

Combination of the above concepts have been used to cast diffraction optics in the framework of post-selected unitary description. Consequently, the mathematical framework of diffraction optics becomes at par with that of other types of interferometry. The classical treatment outlined in the coming sections sets the stage for quantization of fields enabling diffraction optics as an alternative platform for QIP.

As discussed in section 2.1, non-local eigenfunctions of the Helmholtz equation do not make suitable modes for detectors with finite-sized windows. Haar functions (section 2.2) on the other hand, have compact support over a given interval and therefore 2D Haar functions make suitable modes for the square-shaped detector windows. The Haar wavelets and the Haar scaling functions, however, form an overcomplete set of orthonormal functions [30, 31].

To remove the overcompleteness, we divide each slice into non-overlapping square patches, each with side-length equal w, as shown in figure 3. The square patches are labeled using two indices k and $k^{\prime}$ which take integer values.

Zoom In Zoom Out Reset image size

Each square patch supports a countably infinite set of Haar wavelets that fall entirely within the patch. Together with the Haar scaling function that covers the square patch, all the supported Haar wavelets form a basis for any function that has support over the patch. The first element of this basis is the Haar scaling function that covers the entire patch, i.e.

$$\begin{eqnarray}&&{g}_{1}(x,y;{z}_{1},{j}_{0},k,k^{\prime} ):= {\phi }_{{j}_{0},k}(x)\,{\phi }_{{j}_{0},k^{\prime} }(y-{j}_{0}),\end{eqnarray} \tag{ 10 }$$

where the dilation parameter of the scaling function, i.e. j₀ is set so that ${g}_{1}(x,y;{z}_{1},{j}_{0},k,k^{\prime} )$ covers the entire patch (see equation (5)), and z = z₁ is the plane on which the slice is considered. The other elements of the basis are all Haar wavelets with compact support over the square patch, i.e.

$$\begin{eqnarray}&&{g}_{{\imath }}(x,y;{z}_{1},{j}_{0},k,k^{\prime} ):= {\psi }_{m,n}(x){\psi }_{m^{\prime} n^{\prime} }(y-{j}_{0})\,\forall \,{\imath }\gt 1\,\in \,{{\mathbb{Z}}}^{+},\end{eqnarray} \tag{ 11 }$$

where ${{\mathbb{Z}}}^{+}$ is the set of positive integers. The subscript ${\imath }$ is a meta-index for m, n, $m^{\prime}$ and $n^{\prime}$, and

$$\begin{eqnarray}&&m\geqslant {j}_{0},\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{2}^{m-{j}_{0}}k\leqslant n\lt {2}^{m-{j}_{0}}(k+1),\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&m^{\prime} \geqslant {j}_{0},\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&{2}^{m^{\prime} -{j}_{0}^{{\prime} }}k^{\prime} \leqslant n^{\prime} \lt {2}^{m^{\prime} -{j}_{0}^{{\prime} }}(k^{\prime} +1),\end{eqnarray} \tag{ 15 }$$

where the ranges ensure that all the wavelets have compact support over the square patch chosen. If such Haar functions for all the square patches are combined, the slice can be resolved in terms of these functions using equation (8).

Let ${{\mathbb{H}}}_{k,k^{\prime} }({z}_{1})$ denote the space of functions on a slice at z = z₁, with compact support on the patch with index $(k,k^{\prime} )$ and ${{ \mathcal B }}_{k,k^{\prime} }({z}_{1})$ be the basis that spans it, which is given by equations (10) and (11) with constraints given by equations (12)–(15). The space of functions with support over the entire slice is

$$\begin{eqnarray}&&{\mathbb{H}}({z}_{1})=\mathop{\displaystyle \bigoplus }\limits_{k,k^{\prime} }{{\mathbb{H}}}_{k,k^{\prime} }({z}_{1})\end{eqnarray} \tag{ 16 }$$

with a complete basis

$$\begin{eqnarray}&&{ \mathcal B }({z}_{1})=\mathop{\displaystyle \bigoplus }\limits_{k,k^{\prime} }{{ \mathcal B }}_{k,k^{\prime} }({z}_{1})\end{eqnarray} \tag{ 17 }$$

that spans the entire slice.

We elaborate on the slice modes concept using a customized double-slit setup with two sources and two detectors as shown in figure 4, where $\hat{{\boldsymbol{y}}}$ extends into the plane of the paper. The slits are parallel to the xy plane and so are the sources and the detectors at different values of z. The width of the apertures and other distances are chosen such that far-field approximations can be applied to solutions of the Helmholtz equation. A perfectly absorbing barrier is added that runs along the z direction and separates a slit from the detector across the barrier. The purpose of the barrier is to prevent the high diffraction orders [20] from reaching the detectors and also to isolate one detector from another to avoid an overlap of fields between the two.

Zoom In Zoom Out Reset image size

The sources ${{\rm{S}}}_{1}$ and ${{\rm{S}}}_{2}$ are monochromatic point-like sources (practically a spherical source with diameter ∼λ) emanating linearly polarized light as spherical waves (${Y}_{0}^{0}(\theta ,\phi )$ spherical harmonic [21]) with wavelength λ. They are placed at ${{\boldsymbol{r}}}_{{{\rm{S}}}_{1}}=\left(d/2,0,-L\right)$ and ${{\boldsymbol{r}}}_{{{\rm{S}}}_{2}}=\left(-d/2,0,-L\right)$ respectively, where d = 20λ is the distance between the two sources and L = 800λ is the distance between the sources and the slit plane along the z direction. Without loss of generality, we choose λ = 1.

In the far-field regime [20], these sources can be approximated by Dirac-delta functions ${\delta }^{3}({\boldsymbol{r}}-{{\boldsymbol{r}}}_{{{\rm{S}}}_{1}})$ and ${\delta }^{3}({\boldsymbol{r}}-{{\boldsymbol{r}}}_{{{\rm{S}}}_{2}})$. We multiply the source-term with a factor of 10⁵ so that the simulation results do not suffer precision errors. As the sources are linearly polarized, the field from source S_i at points far from the slits, before diffraction, can be found by solving the Helmholtz equation for scalar fields and can be approximated by

$$\begin{eqnarray}&&{E}^{(i)}({\boldsymbol{r}})\approx G({\boldsymbol{r}},{{\boldsymbol{r}}}_{{{\rm{S}}}_{i}}),\end{eqnarray} \tag{ 18 }$$

where the use of the scalar equation is justified because the polarizations of field from both the sources are collinear. Note that the approximation in equation (18) is valid only because the slit plane is far enough from the plane at z₁, so that the surface-effects are negligible.

Two square-shaped slits ${{\rm{A}}}_{1}$ and ${{\rm{A}}}_{2}$, each with side-length w = 4λ are placed at z = 0 with the positions of their centers ${{\boldsymbol{r}}}_{{{\rm{A}}}_{1}}=\left(d/2,0,0\right)$ and ${{\boldsymbol{r}}}_{{{\rm{A}}}_{2}}=\left(-d/2,0,0\right)$ respectively, and therefore aligned with the respective sources. In the Fraunhofer regime, the diffracted field ${E}_{j}^{(i)}({\boldsymbol{r}})$ from source S_i through slit A_j, is calculated by simplifying the surface term in equation (1) by applying the appropriate approximations (see appendix A). Note that, due to the opaque barrier between the two detectors, detector port D_k is blocked from the field ${E}_{j}^{(i)}({\boldsymbol{r}})$ if $k\ne j$.

Another pair of square-shaped slits ${{\rm{D}}}_{1}$ and ${{\rm{D}}}_{2}$, aligned with ${{\rm{A}}}_{1}$ and ${{\rm{A}}}_{2}$ respectively, are placed at z = L, behind each of which is a square-law detector (that measures the integrated squared magnitude of fields) whose window is of the same shape and size as those of the slits. The detector is 100% efficient for light of wavelength λ. The detectors could have been placed without the second pair of slits which play a role only when two such double-slit setups are concatenated to construct interferometers.

Two imaginary planes are at z₁ = −0.9L and z₂ = L on which the input and output slices are considered respectively, for the double-slit. The input slice is not placed at −L where the sources are, as the solution of the Helmholtz equation diverges at the sources. The 2D wavelets are used as basis functions (as discussed in section 3.1) on each of these slices so that each slice can be represented as a column vector. By choosing the input and the output slices on either side of the slits, a transfer matrix mapping the input slice to the output slice can be calculated. However, the slice-to-slice map can also be done in a continuous manner where a propagator sequentially maps one slice to another very close to it, gradually moving forward in the z direction. Such a map is constructed using the surface term of the formal solution of the Helmholtz equation (see equation (1)), which is discussed in appendix B. However, for the purpose of showing that a double-slit is effectively a beam splitter, a direct transfer matrix between the input and output slices suffices.

The two detector-windows in the double-slit setup cover only two of the square patches. Consequently, they do not intercept the entire slice, but only a portion of it. Nevertheless, each detector window supports a countably infinite number of Haar functions. We justify the post-selection of two input and two output modes.

3.3.1. Output modes

As the width of the detector is w = 4, we set j₀ = −2 so that the Haar scaling function covers the entire patch. According to the positions of the detectors in figure 4, the square patch occupied by the detector at port D₁ is the one with indices $k=2,k^{\prime} =0$. Similarly, the patch covered by the detector at port D₂ is the one with $k=-3,k^{\prime} =0$. From the infinite set of Haar functions supported by the square patches, two have to be post-selected. There are two ways to achieve that.

One way is to design a detector that responds to the projection of light on a particular Haar wavelet or scaling function. Although possible in principle, making such a detector is practically challenging because of the jump discontinuities in the Haar wavelet functions. Another and more tractable approach is to construct the double-slit setup in such a way that most of the light intercepted by the detectors has projection on a single Haar wavelet or the scaling function, which becomes a detector mode. Consequently, even if the detector is multimode, the detection is in single mode. The latter is the case with the double-slit setup considered in this work.

To find such modes, the diffracted fields intercepted by the detector windows (ports D₁ and D₂) are resolved in terms of the corresponding Haar functions. For example, the field ${E}_{1}^{(1)}(x,y;{z}_{2})$ can be expanded as

$$\begin{eqnarray}&&{E}_{1}^{(1)}(x,y;{z}_{2})=\displaystyle \sum _{{\imath }=1}^{\infty }{A}_{{\imath }}({z}_{2})\,{g}_{{\imath }}(x,y;{z}_{2},-2,2,0),\end{eqnarray} \tag{ 19 }$$

where ${A}_{{\imath }}({z}_{2})$ are the projections of the field on the corresponding Haar function. For convenient visualization, the field at y = 0, i.e. ${E}_{1}^{(1)}(x,0;{z}_{2})$ is shown in figures 5 and 6.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Although it takes more than one basis function to capture the fine features of the field, most of the power of the light resides in just one mode, i.e. g₁(x, y; z₂, −2, 2, 0). The proof of this fact is in table 1, where the total power intercepted by port D₁ is compared with the total power in the projection on g₁(x, y; z₂, −2, 2, 0). The calculation shows that about 99.9956% of the total power is in the said projection. The field ${E}_{1}^{(2)}(x,y;{z}_{2})$ also has most of the power in this mode. The large projection of the field, intercepted by the detectors, on the Haar scaling function is a consequence of the far-field design of the setup. As the detector is small compared to its distance from the slits, the field intercepted by it varies very slowly over its window. As a result, the overlap of the field with the Haar scaling function, which is constant over the detector window, is dominant.

Table 1.  A comparison of the total integrated intensity detected by the detectors and the square of the magnitudes of the projections of the fields on the chosen modes. The ratios show that most of the light intercepted by the detectors are in the chosen modes as in equations (20) and (21). Therefore the choice of modes is justified.

Integrated intensity intercepted	Integrated projections on chosen modes	Ratio
${\iint }_{{{\rm{D}}}_{1}}{\rm{d}}x{\rm{d}}y\,{\left|{\bar{E}}_{1}^{(1)}(x,y;{z}_{2})\right|}^{2}=0.633\,115$	${\iint }_{{{\rm{D}}}_{1}}{\rm{d}}x{\rm{d}}y\,{\left|,\left\langle {e}_{1}({z}_{2}),{\bar{E}}_{1}^{(1)}({z}_{2})\right\rangle {e}_{1}(x,y;{z}_{2})\right|}^{2}=0.633\,087$	99.995 6%
${\iint }_{{{\rm{D}}}_{2}}{\rm{d}}x{\rm{d}}y\,{\left|{\bar{E}}_{2}^{(1)}(x,y;{z}_{2})\right|}^{2}=0.612\,00$	${\iint }_{{{\rm{D}}}_{2}}{\rm{d}}x{\rm{d}}y\,{\left|,\left\langle {e}_{2}({z}_{2}),{\bar{E}}_{2}^{(1)}({z}_{2})\right\rangle {e}_{2}(x,y;{z}_{2})\right|}^{2}=0.611\,971$	99.995 2%
${\iint }_{{{\rm{D}}}_{1}}{\rm{d}}x{\rm{d}}y\,{\left|{\bar{E}}_{1}^{(2)}(x,y;{z}_{2})\right|}^{2}=0.612\,00$	${\iint }_{{{\rm{D}}}_{1}}{\rm{d}}x{\rm{d}}y\,{\left|,\left\langle {e}_{1}({z}_{2}),{\bar{E}}_{1}^{(2)}({z}_{2})\right\rangle {e}_{1}(x,y;{z}_{2})\right|}^{2}=0.611\,971$	99.995 2%
${\iint }_{{{\rm{D}}}_{2}}{\rm{d}}x{\rm{d}}y\,{\left|{\bar{E}}_{2}^{(2)}(x,y;{z}_{2})\right|}^{2}=0.633\,115$	${\iint }_{{{\rm{D}}}_{2}}{\rm{d}}x{\rm{d}}y\,{\left|,\left\langle {e}_{2}({z}_{2}),{\bar{E}}_{2}^{(2)}({z}_{2})\right\rangle {e}_{2}(x,y;{z}_{2})\right|}^{2}=0.633\,087$	99.995 6%

Therefore, we ignore all the other Haar functions which have negligible contribution to the field. Similarly, for detector port D₂, the dominant contribution is from g₁(x, y; z₂, −2, −3,0). Henceforth, two post-selected output modes are

$$\begin{eqnarray}&&{e}_{1}(x,y;{z}_{2}):= {g}_{1}(x,y;{z}_{2},-2,2,0),\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&{e}_{2}(x,y;{z}_{2}):= {g}_{1}(x,y;{z}_{2},-2,-3,0),\end{eqnarray} \tag{ 21 }$$

which reduce the infinite-dimensional representation of the slice at z = z₂ to a 2D column vector

$$\begin{eqnarray}&&{{\boldsymbol{Y}}}^{(i)}({z}_{2})=\left(\begin{array}{c}\left\langle {e}_{1}({z}_{2}),{E}_{1}^{(i)}({z}_{2})\right\rangle \\ \left\langle {e}_{2}({z}_{2}),{E}_{2}^{(i)}({z}_{2})\right\rangle \end{array}\right),\end{eqnarray} \tag{ 22 }$$

where ${E}_{j}^{(i)}(x,y;{z}_{2})$ is the field from source S_i diffracted by slit A_j and therefore intercepted by the detector port D_j, and

$$\begin{eqnarray}&&\left\langle {e}_{1}({z}_{2}),{E}_{1}^{(i)}({z}_{2})\right\rangle ={\int }_{-\infty }^{\infty }{\rm{d}}x{\int }_{-\infty }^{\infty }{\rm{d}}y\,{e}_{1}^{* }(x,y;{z}_{2})\,{E}_{1}^{(i)}(x,y;{z}_{2})\end{eqnarray} \tag{ 23 }$$

and similarly for the second entry in the column vector.

3.3.2. Input modes

Each source in the double-slit setup (figure 4) emits light in a particular mode, for example spherical mode (approximately). The projection of the source modes on the slice at z = z₁ are projected on the Haar functions on that plane. The input ports chosen are square patches on the plane at z₁ centered at the same positions on the plane as the output ports are placed on plane at z₂. Moreover, the post-selected input modes are similar to those chosen for the output, i.e.

$$\begin{eqnarray}&&{e}_{1}(x,y;{z}_{1}):= {g}_{1}(x,y;{z}_{1},-2,2,0),\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}&&{e}_{2}(x,y;{z}_{1}):= {g}_{1}(x,y;{z}_{1},-2,-3,0),\end{eqnarray} \tag{ 25 }$$

with z₁ denoting that the modes are for a slice on plane z = z₁. Note that the input and output modes are distinguished using the parameter that denotes the plane on which the slice is, i.e. z. Therefore, the post-selected 2D vector representation of the input, i.e. slice at z = z₁

$$\begin{eqnarray}&&{{\boldsymbol{X}}}^{(i)}({z}_{1})=\left(\begin{array}{c}\left\langle {e}_{1}({z}_{1}),{E}_{1}^{(i)}({z}_{1})\right\rangle \\ \left\langle {e}_{2}({z}_{1}),{E}_{2}^{(i)}({z}_{1})\right\rangle \end{array}\right),\end{eqnarray} \tag{ 26 }$$

where the superscript denotes that source S_i is turned on.

The transfer matrix T(z₂, z₁) must map the input vector ${{\boldsymbol{X}}}^{(i)}({z}_{1})$ to the output vector ${{\boldsymbol{Y}}}^{(i)}({z}_{2})$, i.e.

$$\begin{eqnarray}&&T({z}_{2},{z}_{1})\,{{\boldsymbol{X}}}^{(i)}({z}_{1})={{\boldsymbol{Y}}}^{(i)}({z}_{2}),\end{eqnarray} \tag{ 27 }$$

for i ∈ {1, 2}, which gives a set of four simultaneous equations for the four elements of T(z₂, z₁). To solve the equations, the double-slit setup is characterized numerically, by calculating ${{\boldsymbol{X}}}^{(i)}({z}_{1})$ and ${{\boldsymbol{Y}}}^{(i)}({z}_{2})$ for each source turned on at a time. The transformation equation for both sources are combined into one matrix equation

$$\begin{eqnarray}&&T({z}_{2},{z}_{1})\left({{\boldsymbol{X}}}^{(1)}({z}_{1})\ \ {{\boldsymbol{X}}}^{(2)}({z}_{1},),)\right)=\left({{\boldsymbol{Y}}}^{(1)}({z}_{2})\ \ {{\boldsymbol{Y}}}^{(2)}({z}_{2},),)\right),\end{eqnarray} \tag{ 28 }$$

where $\left({{\boldsymbol{X}}}^{(1)}({z}_{1})\ \ {{\boldsymbol{X}}}^{(2)}({z}_{1},),)\right)$ is a 2 × 2 matrix with ${{\boldsymbol{X}}}^{(1)}({z}_{1})$ and ${{\boldsymbol{X}}}^{(2)}({z}_{1})$ as columns, and similar for the right-hand side. Inverting equation (28) yields

$$\begin{eqnarray}&&T({z}_{2},{z}_{1})=\left({{\boldsymbol{Y}}}^{(1)}({z}_{2})\ \ {{\boldsymbol{Y}}}^{(2)}({z}_{2}\right)\,{\left({{\boldsymbol{X}}}^{(1)}({z}_{1}){{\boldsymbol{X}}}^{(2)}({z}_{1},),)\right)}^{-1},\end{eqnarray} \tag{ 29 }$$

provided that $\left({{\boldsymbol{X}}}^{(1)}({z}_{1})\ \ {{\boldsymbol{X}}}^{(2)}({z}_{1},),)\right)$ is invertible. In general, $\left({{\boldsymbol{X}}}^{(1)}({z}_{1})\ \ {{\boldsymbol{X}}}^{(2)}({z}_{1},),)\right)$ is a symmetric matrix because of the symmetry in the setup, and the diagonal elements are slightly different from the off-diagonal elements as the projections of field from one source is not equal on both the post-selected modes. Such a matrix is always invertible.

However, the effective transfer matrix is not unitary because diffraction is intrinsically a lossy process in which most of the light incident on the slits are blocked by the opaque areas. Moreover, the wavelets, chosen as bases of the Hilbert space of each slice, are not eigenfunctions of the Helmholtz equation. Therefore, there is cross-talk between different modes as one moves from one slice to another. To reveal the underlying unitary transformation, a polar decomposition of the transfer matrix is done. Such a decomposition factorizes the non-unitary transfer matrix into a unitary matrix and a Hermitian matrix.

To check the efficacy of this beam splitter, we study the cross-correlation of the outputs. As the solutions of the Helmholtz equation are time-independent, the cross-correlation is modified such that the distinguishing parameter between the inputs is the angle of polarization [41] instead of the time-delay. Further, we concatenate two such double-slit beam splitters to construct a MZI.

3.5.1. Cross-correlation of the post-selected output fields

Let the field from source S₂ have a phase φ with respect to that from source S₁. Also, we rotate the polarization of source S₂ so that θ is the angle between the directions of polarizations of the fields from the two sources. When both the sources are turned on, port D₁ intercepts the vector superposition of fields from sources S₁ and S₂ through slit A₁. The integrated intensity in the post-selected mode on port D₁, i.e. e₁(x, y; z₂) is the projection of the vector sum of the fields intercepted by the port, i.e.

$$\begin{eqnarray}\begin{array}{rcl}{\parallel {E}_{+}({z}_{2},\theta ,\varphi )\parallel }^{2} & := & {\iint }_{{z}_{2}}{\rm{d}}x{\rm{d}}y\,{\left|{E}_{+}(x,y,{z}_{2},\theta ,\varphi )\right|}^{2}\\ & = & {\left|,\left\langle {e}_{1}({z}_{2}),{E}_{1}^{(1)}({z}_{2})\right\rangle ,\right|}^{2}+{\left|,\left\langle {e}_{1}({z}_{2}),{E}_{1}^{(2)}({z}_{2})\right\rangle ,\right|}^{2}\\ & & +2\mathrm{Re}\left\{{\left\langle {e}_{1}({z}_{2}),{E}_{1}^{(1)}({z}_{2})\right\rangle }^{* }\left\langle {e}_{1}({z}_{2}),{E}_{1}^{(2)}({z}_{2})\right\rangle \right\}\cos \varphi \cos \theta ,\end{array}\end{eqnarray} \tag{ 30 }$$

where θ is the distinguishability parameter between the two sources. Similarly, the integrated intensity in the post-selected mode on port D₂ is calculated by projecting E₋ (x, y; z₂, θ, φ) on to e₂(x, y; z₂).

The cross-correlation between the two outputs as a function of θ in equation (77) is modified for time-independent fields as

$$\begin{eqnarray}&&C(\theta ;{z}_{2}):= \displaystyle \frac{\int {\rm{d}}\varphi \,p(\varphi ){\parallel {E}_{+}({z}_{2},\theta ,\varphi )\parallel }^{2}{\parallel {E}_{-}({z}_{2},\theta ,\varphi )\parallel }^{2}}{\int {\rm{d}}\varphi \,p(\varphi ){\parallel {E}_{+}({z}_{2},\theta ,\varphi )\parallel }^{2}\int {\rm{d}}\varphi \,p(\varphi ){\parallel {E}_{-}({z}_{2},\theta ,\varphi )\parallel }^{2}}\end{eqnarray} \tag{ 31 }$$

which is the intensity–intensity correlation of the two output modes of the slice at z = z₂.

3.5.2. Effective MZI

The double-slit beam splitters, discussed in this work, can be concatenated to construct more sophisticated interferometers. As an example, figure 7 shows the schematic of an effective MZI [28] made by concatenating two such double-slit modules.

Zoom In Zoom Out Reset image size

The detectors behind ports D₁ and D₂ as in figure 4 are removed, and these ports now serve as inputs to the second double-slit module. The fields from these ports get diffracted by slits ${{\rm{A}}}_{1}^{{\prime} }$ and ${{\rm{A}}}_{2}^{{\prime} }$ and reach the output ports ${{\rm{D}}}_{1}^{{\prime} }$ and ${{\rm{D}}}_{2}^{{\prime} }$, of the second double-slit module behind each of which is a detector.

A phase shifter (see appendix D on how the phase-shifter is implemented numerically) causes an interference pattern at the output ports ${{\rm{D}}}_{1}^{{\prime} }$ and ${{\rm{D}}}_{2}^{{\prime} }$ as the phase, say α is changed in one of the arms of the MZI. The interfernce pattern obtained is used as a signature to verify the double-slit based MZI.

The MZI is essentially a concatenation of two beam splitters. If both the beam splitters are identical 50:50 splitters with transfer matrix in equation (9) and the phase in one arm, say α is set to zero (the arm lengths are considered equal), the transfer matrix for the MZI is [28]

$$\begin{eqnarray}\displaystyle \frac{1}{2}\left(\begin{array}{cc}1 & {\rm{i}}\\ {\rm{i}} & 1\end{array}\right)\left(\begin{array}{cc}1 & {\rm{i}}\\ {\rm{i}} & 1\end{array}\right)=\left(\begin{array}{cc}0 & i\\ i & 0\end{array}\right),\end{eqnarray} \tag{ 32 }$$

and therefore the transfer matrix for the double-slit MZI should be close to this. We use a similar approach as that used for the double-slit setup, to find the transfer matrix for the effective MZI, with the output slice at z₄ (as shown in figure 7).

A beam splitter is a two-input-two-output device which, as discussed above, can be effectively constructed using double-slit diffraction. However, one of the key potential uses of slit-diffraction and the framework outlined in this work is extension to higher dimensions. As an example, figure 8 shows a triple-slit setup in which a third source S₃ is placed at (−3d/2, 0, −L), a slit A₃ centered at (−3d/2, 0, 0) and a detector D₃ centered at (−3d/2, 0, L). Similar to the double-slit case, the transfer matrix approach can be applied to this system.

Zoom In Zoom Out Reset image size

For the triple-slit setup, three Haar scaling functions are chosen as post-selected input modes and three for the post-selected output modes. According to the positions of the detectors (and slits) these modes are

$$\begin{eqnarray}&&{e}_{1}(x,y;{z}_{1}):= {g}_{1}(x,y;{z}_{1},-2,2,0),\end{eqnarray} \tag{ 33 }$$

$$\begin{eqnarray}&&{e}_{2}(x,y;{z}_{1}):= {g}_{1}(x,y;{z}_{1},-2,-3,0),\end{eqnarray} \tag{ 34 }$$

$$\begin{eqnarray}&&{e}_{3}(x,y;{z}_{1}):= {g}_{1}(x,y;{z}_{1},-2,-8,0),\end{eqnarray} \tag{ 35 }$$

and similarly for slice at z₂.

Like in the case of two slits, the post-selected input and output will have a 3D column representation similar to those in equations (22) and (26). The equation for the effective transfer matrix is

$$\begin{eqnarray}&&T({z}_{2},{z}_{1})\left({{\boldsymbol{X}}}^{(1)}({z}_{1})\ \ {{\boldsymbol{X}}}^{(2)}({z}_{1})\ \ {{\boldsymbol{X}}}^{(3)}({z}_{1},),)\right)=\left({{\boldsymbol{Y}}}^{(1)}({z}_{2})\ \ {{\boldsymbol{Y}}}^{(2)}({z}_{2})\ \ {{\boldsymbol{Y}}}^{(3)}({z}_{2},),)\right),\end{eqnarray} \tag{ 36 }$$

where the superscripts denote the source that is turned on.

As discussed in section 3.4, the effective 2 × 2 transfer matrix is calculated by characterizing the double-slit setup by turning one source on at a time. For each source the input and output slices have a post-selected vector representation. For the double-slit setup under consideration, the input and output vectors when source S₁ is on are

$$\begin{eqnarray}&&{{\boldsymbol{X}}}^{(1)}({z}_{1})\approx \left(\begin{array}{c}-394.761-41.473{\rm{i}}\\ 10.284-13.398{\rm{i}}\end{array}\right),\end{eqnarray} \tag{ 37 }$$

$$\begin{eqnarray}&&{{\boldsymbol{Y}}}^{(1)}({z}_{2})\approx \left(\begin{array}{c}0.008-0.796{\rm{i}}\\ 0.782\,+\,0.008{\rm{i}}\end{array}\right),\end{eqnarray} \tag{ 38 }$$

with respect to the post-selected modes (see section 3.3). Similarly, the post-selected vector representations of the input and output slice when source S₂ is on are

$$\begin{eqnarray}&&{{\boldsymbol{X}}}^{(2)}({z}_{1})\approx \left(\begin{array}{c}10.284-13.398{\rm{i}}\\ -394.761-41.473{\rm{i}}\end{array}\right),\end{eqnarray} \tag{ 39 }$$

$$\begin{eqnarray}&&{{\boldsymbol{Y}}}^{(2)}({z}_{2})\approx \left(\begin{array}{c}0.782\,+\,0.008{\rm{i}}\\ 0.008-0.796{\rm{i}}\end{array}\right).\end{eqnarray} \tag{ 40 }$$

We use the above results in equation (29) to calculate the effective transfer matrix.

With respect to the post-selected input and output vectors, the effective transfer matrix is

$$\begin{eqnarray}T({z}_{2},{z}_{1})\approx {{\rm{e}}}^{0.476\pi {\rm{i}}}\left(\begin{array}{cc}2.07 & 1.90{{\rm{e}}}^{0.486\pi {\rm{i}}}\\ 1.90{{\rm{e}}}^{0.486\pi {\rm{i}}} & 2.07\end{array}\right)\times {10}^{-3},\end{eqnarray} \tag{ 41 }$$

which is a symmetric matrix as expected from the symmetry of the double-slit setup (figure 4). Note that apart from a factor of about ${{\rm{e}}}^{0.476\pi {\rm{i}}}\times 2\sqrt{2}\times {10}^{-3}$, the matrix T(z₂, z₁) is approximately (but not exactly) a 50:50 beam splitter matrix. The deviation from the ideal 50:50 beam splitter is due to the non-unitary nature of the transfer matrix. To reveal the exact unitary transformation the polar decomposition of the transfer matrix is performed.

As expected, the transfer matrix T(z₂, z₁) is not unitary as can be seen from

$$\begin{eqnarray}T({z}_{2},{z}_{1}){T}^{\dagger }({z}_{2},{z}_{1})\approx \left(\begin{array}{cc}7.90 & 0.35\\ 0.35 & 7.90\end{array}\right)\times {10}^{-6}\ne \left(\begin{array}{cc}1 & 0\\ 0 & 1\end{array}\right),\end{eqnarray} \tag{ 42 }$$

because of reasons discussed in section 3.4. To reveal the underlying unitary transformation a polar decomposition of T(z₂, z₁) is carried out which yields

$$\begin{eqnarray}&&T({z}_{2},{z}_{1})=U({z}_{2},{z}_{1})P({z}_{2},{z}_{1}),\end{eqnarray} \tag{ 43 }$$

where the result for the double-slit system considered in this paper, i.e. the polar decompostion of the transfer matrix in equation (41) is

$$\begin{eqnarray}U({z}_{2},{z}_{1})\approx {{\rm{e}}}^{0.47\pi {\rm{i}}}\times \displaystyle \frac{1}{\sqrt{2}}\left(\begin{array}{cc}1.04 & 0.95{\rm{i}}\\ 0.95{\rm{i}} & 1.04\end{array}\right),\end{eqnarray} \tag{ 44 }$$

$$\begin{eqnarray}P({z}_{2},{z}_{1})\approx \left(\begin{array}{cc}2.81 & 0.06\\ 0.06 & 2.81\end{array}\right)\times {10}^{-3},\end{eqnarray} \tag{ 45 }$$

where U is the transformation of a 54:46 beam splitter upto a global phase (which is irrelevant as the detectors are square-law type). The Hermitian component P(z₂, z₁) captures the non-unitarity of the transfer matrix. Its diagonal elements show the fraction of the input that is detected by the detectors after post-selection. The off-diagonal terms show cross-talk between the two modes.

Therefore, the double-slit setup in figure 4 is effectively a lossy beam splitter with respect to the post-selected input and output modes. To verify this result, cross-correlation of the post-selected outputs is calculated and the result is compared with what is expected from an ideal beam splitter (see appendix C).

Here we show the result of the numerically calculated intensity–intensity cross-correlation of the post-selected outputs using equation (31) with

$$\begin{eqnarray}&&p(\varphi )=\displaystyle \frac{1}{2\pi },\end{eqnarray} \tag{ 46 }$$

i.e. the relative phase between the two sources is uniformly random. Figure 9 shows the values of the correlation as a function of the distinguishability parameter, which in this case is the relative polarization angle θ between the two sources. The function that fits the result is

$$\begin{eqnarray}&&{C}_{50}(\theta )=0.75-0.25\,\cos 2\theta ,\end{eqnarray} \tag{ 47 }$$

the visibility of which is 0.5.

Zoom In Zoom Out Reset image size

On the other hand, for a 50:50 cubic beam splitter (with transfer matrix as in equation (9)), if the relative phase φ between the two sources is distributed uniformly over the interval [0, 2π), the cross-correlation function shows a visibility of 0.5 as shown in figure 10 (see appendix C for more details on the correlation of the outputs of a cubic beam splitter). On comparing the variation of the correlation of the outputs of the double-slit setup with that of the cubic beam splitter, we confirm the beam splitter like behavior of the double-slit setup.

Zoom In Zoom Out Reset image size

For completeness, appendix E discusses the 100% dip in the correlation by using a suitable probability distribution of phase, as suggested in [40].

Here we show the interference pattern at the output of the double-slit based MZI as shown in figure 7. Similar to the case of one double-slit setup, adopting the Fraunhofer approximation and calculating the integrated intensities at the two detectors for different values of α yields an interference pattern shown in figure 11. Such an interference pattern is a signature of an MZI.

Zoom In Zoom Out Reset image size

The curves that fit the resultant integrated intensities at ports ${{\rm{D}}}_{1}^{{\prime} }$ and ${{\rm{D}}}_{2}^{{\prime} }$ are approximately

$$\begin{eqnarray}&&{I}_{1}(\alpha )=2.41\times {10}^{-8}\left(1-\sin \left(x-2.45\times {10}^{-4}\right)\right),\end{eqnarray} \tag{ 48 }$$

$$\begin{eqnarray}&&{I}_{2}(\alpha )=2.41\times {10}^{-8}\left(1+\sin \left(x+2.45\times {10}^{-4}\right)\right),\end{eqnarray} \tag{ 49 }$$

respectively. The visibility V_MZI of both the curves is

$$\begin{eqnarray}&&V=99.94 \% ,\end{eqnarray} \tag{ 50 }$$

which means that the effective MZI using double-slit modules closely emulates an MZI with cubic beam splitters. Therefore, the transfer matrix formalism applied to the setup in figure 7 should yield the transformation in equation (32).

Using the method outlined in this paper, the transfer matrix for the double-slit Mach–Zehnder in figure 7 is

$$\begin{eqnarray}{T}_{\mathrm{MZ}}\approx {{\rm{e}}}^{0.08\pi {\rm{i}}}\left(\begin{array}{cc}0.255 & 3.903\,{{\rm{e}}}^{0.39\pi {\rm{i}}}\\ 3.903\,{{\rm{e}}}^{0.39\pi {\rm{i}}} & 0.255\end{array}\right)\times {10}^{-7},\end{eqnarray} \tag{ 51 }$$

which after polar decomposition yields

$$\begin{eqnarray}{U}_{\mathrm{MZ}}\approx \left(\begin{array}{cc}0.061 & 0.998{\rm{i}}\\ 0.998{\rm{i}} & 0.061\end{array}\right)\times {{\rm{e}}}^{-0.03\pi {\rm{i}}},\end{eqnarray} \tag{ 52 }$$

$$\begin{eqnarray}{P}_{\mathrm{MZ}}\approx \left(\begin{array}{cc}3.910 & 0.088\\ 0.088 & 3.910\end{array}\right)\times {10}^{-7}.\end{eqnarray} \tag{ 53 }$$

Up to a global phase the transfer matrix method successfully reveals the underlying unitary operator for the double-slit Mach–Zehnder which can be checked by comparing equations (32) and (52).

We have verified our formalism using a double-slit setup and by comparing it to a well-known optical device, the cubic beam splitter. Here, we apply the transfer matrix formalism to get the post-selected transfer matrix of a triple-slit setup, demonstrating the extensibiltiy of the framework to higher-dimensional systems.

The effective 3 × 3 transfer matrix for the triple-slit setup in figure 8, with respect to the post-selected modes discussed in section 3.6 is

$$\begin{eqnarray}{T}_{3}({z}_{2},{z}_{1})\approx {{\rm{e}}}^{-{\rm{i}}0.16\pi }\,\left(\begin{array}{ccc}0.636 & 0.629{\rm{i}} & 0.594\\ 0.626{\rm{i}} & 0.637 & 0.626{\rm{i}}\\ 0.594 & 0.629{\rm{i}} & 0.636\end{array}\right)\times {10}^{-3},\end{eqnarray} \tag{ 54 }$$

where subscript 3 denotes that the transfer matrix is for a triple-slit setup. The polar decomposition of the 3 × 3 transfer matrix results in

$$\begin{eqnarray}{U}_{3}({z}_{2},{z}_{1})\approx \displaystyle \frac{{{\rm{e}}}^{-{\rm{i}}0.15\pi }}{\sqrt{3}}\,\left(\begin{array}{ccc}1.493 & 0.844{\rm{i}} & -0.242\\ 0.844{\rm{i}} & 1.255 & 0.844{\rm{i}}\\ -0.242 & 0.844{\rm{i}} & 1.493\end{array}\right),\end{eqnarray} \tag{ 55 }$$

$$\begin{eqnarray}{P}_{3}({z}_{2},{z}_{1})\approx \left(\begin{array}{ccc}0.772 & 0.146{\rm{i}} & 0.730\\ -0.146{\rm{i}} & 1.075 & -0.146{\rm{i}}\\ 0.730 & 0.146{\rm{i}} & 0.772\end{array}\right)\times {10}^{-3},\end{eqnarray} \tag{ 56 }$$

which reveals the underlying unitary transformation along with the losses captured by the Hermitian component.

By further increasing the number of slits in a similar fashion, even higher-dimensional transfer matrices can be realized using slit-diffraction. Therefore, unlike with four-port devices like beam splitters, which have to be concatenated to implement higher-dimensional transformations, a single N-slit setup can be used for an N-dimensional transformation. This way of implementing higher-dimensional transfer matrices should prove to be easier than that using OAM of light, because of the practical limits on obtaining high OAM states [7].

The complete solution of the Helmholtz equation (equation (57)) in a volume ${ \mathcal V }$ enclosed by a surface $\partial { \mathcal V }$ is

$$\begin{eqnarray}&&E({\boldsymbol{r}})={\int }_{{ \mathcal V }}{{\rm{d}}}^{3}{\boldsymbol{r}}^{\prime} G({\boldsymbol{r}},{\boldsymbol{r}}^{\prime} )\,{\delta }^{3}({\boldsymbol{r}}^{\prime} -{{\boldsymbol{r}}}_{{{\rm{S}}}_{1}})+{\int }_{\partial { \mathcal V }}{{\rm{d}}}^{2}{\boldsymbol{r}}^{\prime} \hat{{\boldsymbol{n}}}({\boldsymbol{r}}^{\prime} )\cdot \left(E({\boldsymbol{r}}^{\prime} ){\rm{\nabla }}^{\prime} G({\boldsymbol{r}},{\boldsymbol{r}}^{\prime} )-G({\boldsymbol{r}},{\boldsymbol{r}}^{\prime} ){\rm{\nabla }}^{\prime} E({\boldsymbol{r}}^{\prime} ,),)\right).\end{eqnarray} \tag{ 58 }$$

Consider the surface $\partial { \mathcal V }$ to be at infinity such that there is no contribution from the surface term in equation (58). As we have assumed that the opaque portions of the slit-plane are perfectly absorbing, the field within the area of the slit is

$$\begin{eqnarray}&&E({\boldsymbol{r}})=G({\boldsymbol{r}},{{\boldsymbol{r}}}_{f})\end{eqnarray} \tag{ 59 }$$

and zero outside the area of the slit. Here

$$\begin{eqnarray}&&G({\boldsymbol{r}},{\boldsymbol{r}}^{\prime} )=-\displaystyle \frac{1}{4\pi }\,\displaystyle \frac{{{\rm{e}}}^{{\rm{i}}k\parallel {\boldsymbol{r}}-{\boldsymbol{r}}^{\prime} \parallel }}{\parallel {\boldsymbol{r}}-{\boldsymbol{r}}^{\prime} \parallel }\end{eqnarray} \tag{ 60 }$$

is the Green's function of the Helmholtz equation in three dimension. To find the diffracted field from the slit, consider a semi-infinite volume with the slit-plane as one part of the surface $\partial { \mathcal V }$ and the other parts at infinity as discussed in appendix B. This volume does not contain any source, but the surface on the slit-plane gives contribution from the slits. Therefore, we use the surface term of equation (58) to calculate the diffracted field as

$$\begin{eqnarray}&&{\int }_{z={{\boldsymbol{R}}}_{s}\cdot \hat{{\boldsymbol{z}}}}{{\rm{d}}}^{2}{\boldsymbol{r}}^{\prime} \hat{{\boldsymbol{z}}}\cdot \left(G({\boldsymbol{r}}^{\prime} ,{{\boldsymbol{r}}}_{f}){\rm{\nabla }}^{\prime} G({\boldsymbol{r}},{\boldsymbol{r}}^{\prime} )-G({\boldsymbol{r}},{\boldsymbol{r}}^{\prime} ){\rm{\nabla }}^{\prime} G({\boldsymbol{r}}^{\prime} ,{{\boldsymbol{r}}}_{f},),)\right),\end{eqnarray} \tag{ 61 }$$

where the only contribution is from the slits-plane. To further simplify the expression for the diffracted field, note that

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{2}G({{\boldsymbol{r}}}_{2},{{\boldsymbol{r}}}_{1})=-\displaystyle \frac{1}{4\pi }\,\displaystyle \frac{{\rm{i}}k{{\rm{e}}}^{{\rm{i}}k| {{\boldsymbol{r}}}_{2}-{{\boldsymbol{r}}}_{1}| }}{| {{\boldsymbol{r}}}_{2}-{{\boldsymbol{r}}}_{1}| }-\displaystyle \frac{{{\rm{e}}}^{{\rm{i}}k| {{\boldsymbol{r}}}_{2}-{{\boldsymbol{r}}}_{1}| }}{| {{\boldsymbol{r}}}_{2}-{{\boldsymbol{r}}}_{1}{| }^{2}}\,{{\rm{\nabla }}}_{2}| {{\boldsymbol{r}}}_{2}-{{\boldsymbol{r}}}_{1}| =-\displaystyle \frac{1}{4\pi }\,\displaystyle \frac{{{\rm{e}}}^{{\rm{i}}k| {{\boldsymbol{r}}}_{2}-{{\boldsymbol{r}}}_{1}| }}{| {{\boldsymbol{r}}}_{2}-{{\boldsymbol{r}}}_{1}| }\left({\rm{i}}k-\displaystyle \frac{1}{| {{\boldsymbol{r}}}_{2}-{{\boldsymbol{r}}}_{1}| }\right)\displaystyle \frac{{{\boldsymbol{r}}}_{2}-{{\boldsymbol{r}}}_{1}}{| {{\boldsymbol{r}}}_{2}-{{\boldsymbol{r}}}_{1}| },\end{eqnarray} \tag{ 62 }$$

which means that ${{\rm{\nabla }}}_{1}G({{\boldsymbol{r}}}_{2},{{\boldsymbol{r}}}_{1})=-{{\rm{\nabla }}}_{2}G({{\boldsymbol{r}}}_{2},{{\boldsymbol{r}}}_{1})$ and from equation (60) $G({{\boldsymbol{r}}}_{2},{{\boldsymbol{r}}}_{1})=G({{\boldsymbol{r}}}_{1},{{\boldsymbol{r}}}_{2})$. The far-field is applied by assuming that $| {\rm{i}}k| \gg \tfrac{1}{| {{\boldsymbol{r}}}_{2}-{{\boldsymbol{r}}}_{1}| }$. Then

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{2}G({{\boldsymbol{r}}}_{2},{{\boldsymbol{r}}}_{1})\approx {\rm{i}}k\,G({{\boldsymbol{r}}}_{2},{{\boldsymbol{r}}}_{1})\,\displaystyle \frac{{{\boldsymbol{r}}}_{2}-{{\boldsymbol{r}}}_{1}}{| {{\boldsymbol{r}}}_{2}-{{\boldsymbol{r}}}_{1}| }.\end{eqnarray} \tag{ 63 }$$

With this approximation, equation (61) simplifies to

$$\begin{eqnarray}&&E({\boldsymbol{r}})=-{\rm{i}}k{\int }_{{ \mathcal S }}{{\rm{d}}}^{2}{\boldsymbol{r}}^{\prime} G({\boldsymbol{r}}^{\prime} ,{{\boldsymbol{r}}}_{f})G({\boldsymbol{r}},{\boldsymbol{r}}^{\prime} )\left(\displaystyle \frac{z-z^{\prime} }{| {\boldsymbol{r}}-{\boldsymbol{r}}^{\prime} | }+\displaystyle \frac{z^{\prime} -{z}_{f}}{| {\boldsymbol{r}}^{\prime} -{{\boldsymbol{r}}}_{f}| }\right).\end{eqnarray} \tag{ 64 }$$

The expression for the diffracted field simplifies further when the small-slit approximation is considered. Let ${{\boldsymbol{R}}}_{s}$ be the center of the slit with a width very small compared with its distance from the point of detection. A position within the area of the slit can be written as ${\boldsymbol{r}}^{\prime} ={{\boldsymbol{R}}}_{s}+{{\boldsymbol{\Delta }}}_{s}$. Therefore,

$$\begin{eqnarray}&&G({\boldsymbol{r}}^{\prime} ,{{\boldsymbol{r}}}_{f})=-\displaystyle \frac{1}{4\pi }\,\displaystyle \frac{{{\rm{e}}}^{{\rm{i}}k| {{\boldsymbol{R}}}_{s}-{{\boldsymbol{r}}}_{f}+{{\boldsymbol{\Delta }}}_{s}| }}{| {{\boldsymbol{R}}}_{s}-{{\boldsymbol{r}}}_{f}+{{\boldsymbol{\Delta }}}_{s}| },\end{eqnarray} \tag{ 65 }$$

where the small slit approximation is applied as $| {{\boldsymbol{\Delta }}}_{{\boldsymbol{s}}}| /| {{\boldsymbol{R}}}_{s}-{{\boldsymbol{r}}}_{f}| \ll 1$ so that

$$\begin{eqnarray}\begin{array}{rcl}| {{\boldsymbol{R}}}_{s}-{{\boldsymbol{r}}}_{f}+{{\boldsymbol{\Delta }}}_{s}| & = & \sqrt{| {{\boldsymbol{R}}}_{s}-{{\boldsymbol{r}}}_{f}{| }^{2}+| {{\boldsymbol{\Delta }}}_{s}{| }^{2}+2{{\boldsymbol{\Delta }}}_{s}\cdot ({{\boldsymbol{R}}}_{s}-{{\boldsymbol{r}}}_{f})}=| {{\boldsymbol{R}}}_{s}-{{\boldsymbol{r}}}_{f}| \sqrt{1+\displaystyle \frac{| {{\boldsymbol{\Delta }}}_{s}{| }^{2}}{| {{\boldsymbol{R}}}_{s}-{{\boldsymbol{r}}}_{f}{| }^{2}}+2{{\boldsymbol{\Delta }}}_{s}\cdot \displaystyle \frac{({{\boldsymbol{R}}}_{s}-{{\boldsymbol{r}}}_{f})}{| {{\boldsymbol{R}}}_{s}-{{\boldsymbol{r}}}_{f}{| }^{2}}}\\ & \approx & | {{\boldsymbol{R}}}_{s}-{{\boldsymbol{r}}}_{f}| \sqrt{1+2{{\boldsymbol{\Delta }}}_{s}\cdot \displaystyle \frac{({{\boldsymbol{R}}}_{s}-{{\boldsymbol{r}}}_{f})}{| {{\boldsymbol{R}}}_{s}-{{\boldsymbol{r}}}_{f}{| }^{2}}}\approx | {{\boldsymbol{R}}}_{s}-{{\boldsymbol{r}}}_{f}| \left(1+{{\boldsymbol{\Delta }}}_{s}\cdot \displaystyle \frac{({{\boldsymbol{R}}}_{s}-{{\boldsymbol{r}}}_{f})}{| {{\boldsymbol{R}}}_{s}-{{\boldsymbol{r}}}_{f}{| }^{2}}\right).\end{array}\end{eqnarray} \tag{ 66 }$$

Furthermore, for far-field and small slit

$$\begin{eqnarray}&&\displaystyle \frac{1}{| {\boldsymbol{r}}^{\prime} -{{\boldsymbol{r}}}_{f}| }\approx \displaystyle \frac{1}{| {{\boldsymbol{R}}}_{s}-{{\boldsymbol{r}}}_{f}| }\end{eqnarray} \tag{ 67 }$$

which results in

$$\begin{eqnarray}&&G({\boldsymbol{r}}^{\prime} ,{{\boldsymbol{r}}}_{f})=-\displaystyle \frac{1}{4\pi }\,\displaystyle \frac{{{\rm{e}}}^{{\rm{i}}k| {{\boldsymbol{R}}}_{s}-{{\boldsymbol{r}}}_{f}| }}{| {{\boldsymbol{R}}}_{s}-{{\boldsymbol{r}}}_{f}| }{{\rm{e}}}^{{\rm{i}}k{{\boldsymbol{\Delta }}}_{s}\cdot \tfrac{{{\boldsymbol{R}}}_{s}-{{\boldsymbol{r}}}_{f}}{| {{\boldsymbol{R}}}_{s}-{{\boldsymbol{r}}}_{f}| }}=G({{\boldsymbol{R}}}_{s},{{\boldsymbol{r}}}_{f})\,{{\rm{e}}}^{{\rm{i}}k{{\boldsymbol{\Delta }}}_{s}\cdot \tfrac{{{\boldsymbol{R}}}_{s}-{{\boldsymbol{r}}}_{f}}{| {{\boldsymbol{R}}}_{s}-{{\boldsymbol{r}}}_{f}| }}.\end{eqnarray} \tag{ 68 }$$

Similarly,

$$\begin{eqnarray}&&G({\boldsymbol{r}},{\boldsymbol{r}}^{\prime} )=G({\boldsymbol{r}},{{\boldsymbol{R}}}_{s})\,{{\rm{e}}}^{-{\rm{i}}k{{\boldsymbol{\Delta }}}_{s}\cdot \tfrac{{\boldsymbol{r}}-{{\boldsymbol{R}}}_{s}}{| {\boldsymbol{r}}-{{\boldsymbol{R}}}_{s}| }},\end{eqnarray} \tag{ 69 }$$

and using this approximation in equation (64) we get

$$\begin{eqnarray}&&E({\boldsymbol{r}})=-{\rm{i}}k{\int }_{{ \mathcal S }}{{\rm{d}}}^{2}{{\boldsymbol{\Delta }}}_{s}\,G({\boldsymbol{r}},{{\boldsymbol{R}}}_{s})\,G({{\boldsymbol{R}}}_{s},{{\boldsymbol{r}}}_{f})\,{{\rm{e}}}^{-{\rm{i}}k{{\boldsymbol{\Delta }}}_{s}\cdot \left(\tfrac{{\boldsymbol{r}}-{{\boldsymbol{R}}}_{s}}{| {\boldsymbol{r}}-{{\boldsymbol{R}}}_{s}| }-\displaystyle \frac{{{\boldsymbol{R}}}_{s}-{{\boldsymbol{r}}}_{f}}{| {{\boldsymbol{R}}}_{s}-{{\boldsymbol{r}}}_{f}| }\right)}\,\times \left(\displaystyle \frac{z-{z}_{s}}{| {\boldsymbol{r}}-{{\boldsymbol{R}}}_{s}| }+\displaystyle \frac{{z}_{s}-{z}_{f}}{| {{\boldsymbol{R}}}_{s}-{{\boldsymbol{r}}}_{f}| }\right).\end{eqnarray} \tag{ 70 }$$

For a rectangular slit w_x wide along x and w_y along y, the integral is a sinc function and the final expression for the field of the slice is

$$\begin{eqnarray}\begin{array}{rcl}E({\boldsymbol{r}}) & = & -{\rm{i}}{{kw}}_{x}{w}_{y}\,\left(\displaystyle \frac{z-{z}_{s}}{| {\boldsymbol{r}}-{{\boldsymbol{R}}}_{s}| }+\displaystyle \frac{{z}_{s}-{z}_{f}}{| {{\boldsymbol{R}}}_{s}-{{\boldsymbol{r}}}_{f}| }\right)\,G({\boldsymbol{r}},{{\boldsymbol{R}}}_{s})\,G({{\boldsymbol{R}}}_{s},{{\boldsymbol{r}}}_{f})\times \mathrm{sinc}\,\left(\displaystyle \frac{{{kw}}_{x}}{2}\left(\displaystyle \frac{{\boldsymbol{r}}-{{\boldsymbol{R}}}_{s}}{| {\boldsymbol{r}}-{{\boldsymbol{R}}}_{s}| }-\displaystyle \frac{{{\boldsymbol{R}}}_{s}-{{\boldsymbol{r}}}_{f}}{| {{\boldsymbol{R}}}_{s}-{{\boldsymbol{r}}}_{f}| }\right)\cdot \hat{{\boldsymbol{x}}}\right)\\ & & \times \mathrm{sinc}\,\left(\displaystyle \frac{{{kw}}_{y}}{2}\left(\displaystyle \frac{{\boldsymbol{r}}-{{\boldsymbol{R}}}_{s}}{| {\boldsymbol{r}}-{{\boldsymbol{R}}}_{s}| }-\displaystyle \frac{{{\boldsymbol{R}}}_{s}-{{\boldsymbol{r}}}_{f}}{| {{\boldsymbol{R}}}_{s}-{{\boldsymbol{r}}}_{f}| }\right)\cdot \hat{{\boldsymbol{y}}}\right),\end{array}\end{eqnarray} \tag{ 71 }$$

which with appropriate position vectors of the sources and slit in the double-slit setup yields the diffracted field in that system.