Double-slit interferometry as a lossy beam splitter

TL;DR

This paper models double-slit diffraction as a lossy beam splitter within a quantum framework, enabling potential applications in high-dimensional quantum information processing and offering a robust, miniaturizable alternative to traditional interferometry.

Contribution

It introduces a novel framework for describing diffraction-based interferometry as a post-selected unitary transformation suitable for quantum information processing.

Findings

Diffraction can be represented as an infinite-dimensional transformation.

Post-selection truncates the transformation to a finite-dimensional lossy beam splitter.

Classical diffraction optics can be extended to quantum information applications.

Abstract

We cast diffraction-based interferometry in the framework of post-selected unitary description towards enabling it as a platform for quantum information processing. We express slit-diffraction as an infinite-dimensional transformation and truncate it to a finite-dimensional transfer matrix by post-selecting modes. Using such a framework with classical fields, a customized double-slit setup is effectively a lossy beam splitter in a post-selected sense. Diffraction optics provides a robust alternative to conventional multi-beam interferometry with scope for miniaturization, and also has applications in matter wave interferometry. In this work, the classical treatment sets the stage for quantization of fields and implementing higher dimensional quantum information processing like that done with other platforms such as orbital angular momentum.

Tables

Table 1. 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 Eqs. ( 18 ) and ( 19 ). Therefore the choice of modes is justified.

Integrated intensity intercepted	Integrated projections on chosen modes	Ratio
$\underset{D_1}{\iint }dxdy{\left|{\overline{E}}_1^{\left(1\right)}(x,y;z_2)\right|}^2=0.633115$	$\underset{D_1}{\iint }dxdy{\left|⟨e_1\left(z_2\right),{\overline{E}}_1^{\left(1\right)}\left(z_2\right)⟩e_1(x,y;z_2)\right|}^2=0.633087$	$99.9956\%$
$\underset{D_2}{\iint }dxdy{\left|{\overline{E}}_2^{\left(1\right)}(x,y;z_2)\right|}^2=0.61200$	$\underset{D_2}{\iint }dxdy{\left|⟨e_2\left(z_2\right),{\overline{E}}_2^{\left(1\right)}\left(z_2\right)⟩e_2(x,y;z_2)\right|}^2=0.611971$	$99.9952\%$
$\underset{D_1}{\iint }dxdy{\left|{\overline{E}}_1^{\left(2\right)}(x,y;z_2)\right|}^2=0.61200$	$\underset{D_1}{\iint }dxdy{\left|⟨e_1\left(z_2\right),{\overline{E}}_1^{\left(2\right)}\left(z_2\right)⟩e_1(x,y;z_2)\right|}^2=0.611971$	$99.9952\%$
$\underset{D_2}{\iint }dxdy{\left|{\overline{E}}_2^{\left(2\right)}(x,y;z_2)\right|}^2=0.633115$	$\underset{D_2}{\iint }dxdy{\left|⟨e_2\left(z_2\right),{\overline{E}}_2^{\left(2\right)}\left(z_2\right)⟩e_2(x,y;z_2)\right|}^2=0.633087$	$99.9956\%$

Equations

$$E(\mathbf{r})=\underset{\mathcal{V}}{\iiint }{\mathrm{d}}^3{\mathbf{r}}^{\prime }\text{leavevmode}G(\mathbf{r},{\mathbf{r}}^{\prime })\rho ({\mathbf{r}}^{\prime })+\underset{\partial \mathcal{V}}{\iint }{\mathrm{d}}^2{\mathbf{r}}^{\prime }\text{leavevmode}\hat{n}({\mathbf{r}}^{\prime })\cdot \left(E({\mathbf{r}}^{\prime }){\nabla }^{\prime }G(\mathbf{r},{\mathbf{r}}^{\prime })-G(\mathbf{r},{\mathbf{r}}^{\prime }){\nabla }^{\prime }E({\mathbf{r}}^{\prime })\right),$$

$$\psi (x):=$$

$$\varphi (x):=$$

$${\psi }_{m,n}(x):=$$

$${\varphi }_{j,k}(x):=$$

$$⟨{\psi }_{m,n},{\psi }_{m^{\prime },n^{\prime }}⟩:=\underset{-\infty }{\overset{\infty }{\int }}\mathrm{d}x\text{leavevmode}{\psi }_{m,n}^{\ast }(x){\psi }_{m^{\prime },n^{\prime }}(x)={\delta }_{m,m^{\prime }}{\delta }_{n,n^{\prime }},$$

$$⟨{\varphi }_{j_0,k},{\varphi }_{j_0,k^{\prime }}⟩={\delta }_{k,k^{\prime }},$$

$$f(x)=\sum _{k=-\infty }^{\infty }{a}_k{\varphi }_{j_0,k}(x)+\sum _{m\ge j_0}\sum _{n=-\infty }^{\infty }{b}_{m,n}{\psi }_{m,n}(x),$$

$$U_{\mathrm{BS}}=\frac{1}{\sqrt{2}}\left(\begin{array}{cc}1& \mathrm{i}\\ \mathrm{i}& 1\end{array}\right),$$

$$g_1(x,y;z_1,j_0,k,k^{\prime }):={\varphi }_{j_0,k}(x)\text{leavevmode}{\varphi }_{j_0,k^{\prime }}(y-j_0),$$

$$g_{ı}(x,y;z_1,j_0,k,k^{\prime }):={\psi }_{m,n}(x){\psi }_{m^{\prime }{n}^{\prime }}(y-j_0)\text{leavevmode}\forall \text{leavevmode}ı>1\text{leavevmode}\in \text{leavevmode}{\mathbb{Z}}^{+},$$

$$m\ge j_0,$$

$$2^{m-j_0}k\le$$

$$m^{\prime }\ge j_0,$$

$$2^{m^{\prime }-j_0^{\prime }}{k}^{\prime }\le$$

$$E^{(i)}(\mathbf{r})\approx$$

$$E_1^{(1)}(x,y;z_2)=\sum _{ı=1}^{\infty }{A}_{ı}(z_2)\text{leavevmode}{g}_{ı}(x,y;z_2,-2,2,0),$$

$$e_1(x,y;z_2):=$$

$$e_2(x,y;z_2):=$$

$${\mathbf{Y}}^{(i)}(z_2)=\left(\begin{array}{c}⟨e_1(z_2),E_1^{(i)}(z_2)⟩\\ ⟨e_2(z_2),E_2^{(i)}(z_2)⟩\end{array}\right),$$

$$⟨e_1(z_2),E_1^{(i)}(z_2)⟩=$$

$$e_1(x,y;z_1):=$$

$$e_2(x,y;z_1):=$$

$${\mathbf{X}}^{(i)}(z_1)=\left(\begin{array}{c}⟨e_1(z_1),E_1^{(i)}(z_1)⟩\\ ⟨e_2(z_1),E_2^{(i)}(z_1)⟩\end{array}\right),$$

$$T(z_2,z_1)\text{leavevmode}{\mathbf{X}}^{(i)}(z_1)={\mathbf{Y}}^{(i)}(z_2),$$

$$T(z_2,z_1)\left(\begin{array}{cc}{\mathbf{X}}^{(1)}(z_1)& {\mathbf{X}}^{(2)}(z_1)\end{array}\right)=$$

$$T(z_2,z_1)=\left(\begin{array}{cc}{\mathbf{Y}}^{(1)}(z_2)& {\mathbf{Y}}^{(2)}(z_2)\end{array}\right)\text{leavevmode}{\left(\begin{array}{cc}{\mathbf{X}}^{(1)}(z_1)& {\mathbf{X}}^{(2)}(z_1)\end{array}\right)}^{-1},$$

$${\Vert E_{+}(z_2,\theta ,\phi )\Vert }^2:=$$

$$=$$

$$+2\mathrm{Re}\left\{{⟨e_1(z_2),E_1^{(1)}(z_2)⟩}^{\ast }⟨e_1(z_2),E_1^{(2)}(z_2)⟩\right\}\cos \phi \cos \theta ,$$

$$C(\theta ;z_2):=\frac{\int \mathrm{d}\phi \text{leavevmode}p(\phi ){\Vert E_{+}(z_2,\theta ,\phi )\Vert }^2{\Vert E_{-}(z_2,\theta ,\phi )\Vert }^2}{\int \mathrm{d}\phi \text{leavevmode}p(\phi ){\Vert E_{+}(z_2,\theta ,\phi )\Vert }^2\int \mathrm{d}\phi \text{leavevmode}p(\phi ){\Vert E_{-}(z_2,\theta ,\phi )\Vert }^2}$$

$$\frac{1}{2}\left(\begin{array}{cc}1& \mathrm{i}\\ \mathrm{i}& 1\end{array}\right)\left(\begin{array}{cc}1& \mathrm{i}\\ \mathrm{i}& 1\end{array}\right)=\left(\begin{array}{cc}0& i\\ i& 0\end{array}\right),$$

$$e_1(x,y;z_1):=$$

$$e_2(x,y;z_1):=$$

$$e_3(x,y;z_1):=$$

$$T(z_2,z_1)\left(\begin{array}{ccc}{\mathbf{X}}^{(1)}(z_1)& {\mathbf{X}}^{(2)}(z_1)& {\mathbf{X}}^{(3)}(z_1)\end{array}\right)=$$

Full text

*Email: simanraj@rri.res.in

$\ddagger$Email: sandersb@ucalgary.ca

$\dagger$Email: usinha@rri.res.in

1 Introduction

Linear optical quantum information processing (QIP) [1, 2] has a mathematical representation in the form of finite-dimensional unitary transfer matrices operating on a Hilbert space of vectors that represent qubits/qudits [3, 4]. The qubits are usually encoded in the polarization degree of freedom of a single photon, and optical components like beam splitters [5, 6, 7] and phase-shifters are used to implement the unitary transformations on them. For higher-dimensional QIP, systems such as orbital angular momentum [7, 8, 9, 10, 11] of photons are used. We map diffraction optics over a finite-dimensional unitary representation and connect it to qubit/qudit processing.

The novel interpretation of slit-diffraction that we present here sets the stage for extending the scope of application of diffraction interferometry to modern problems like higher-dimensional information processing. Such a formalism provides an alternative to the implementation of higher-dimensional QIP using the orbital angular momentum (OAM) of light [7]. Slit-diffraction based optical interferometers can be used to construct qudits encoded in spatial modes [12, 13], with robustness, unlike in the case of OAM based qudits which have practical limitations in state-preparation and state-readability [7]. Another potential advantage of the multi-slit-diffraction-based interferometer is scalability of table-top experiments. Moreover, a finite-dimensional unitary description of diffraction also has applications in the field of matter-wave interferometry [14, 15, 16, 17, 18, 19].

We deal with diffraction of classical fields and show a formalism in which slit-diffraction is represented as a finite-dimensional unitary transfer matrix [20] (in the post-selected sense). We project the three-dimensional solutions of the Helmholtz equation [20, 21] on two-dimensional imaginary planes and call these projections slices. The propagation and diffraction of the fields is expressed as a slice-to-slice map as one goes from one slice to another from the sources to the detectors through the slits. By choosing an appropriate basis for the slices, we get an infinite-dimensional transfer matrix representation of such a map. The transfer matrix is reduced to an effective finite-dimensional matrix by post-selecting a finite number of basis elements on the slices as post-selected modes. We show that such a truncated matrix is in general not unitary because of the losses in diffraction, and that the underlying unitary transfer matrix can be revealed by performing a polar decomposition [22] on the effective transfer matrix separating it into unitary and lossy (Hermitian) components.

Using the post-selected unitary transfer matrix formalism of diffraction, we show that a customized double-slit setup is effectively a lossy beam splitter in the classical regime. A cubic beam splitter is a two-input-two-output optical device that has a $2\times 2$ unitary transfer matrix that transforms the fields entering its input ports to the fields exiting its output ports [5, 6]. This 4-port device, along with a phase-shifter which is a 2-port device that imparts a phase, serves as the building blocks of any $N$-channel interferometer [23, 24, 25, 26, 27]. The novelty and importance of our work lies in connecting one of the most elegant and fundamental experiments in scientific history, i.e., double-slit-diffraction with other types of interferometries which are used to solve some of the most important problems in modern physics, like QIP.

To verify the beam splitter like behaviour of the double-slit setup, we compare the correlation of the classical outputs with that of the cubic beam splitter. Moreover, by concatenation of two such double-slit based beam splitters and using a phase-shifter, we construct an effective Mach-Zehnder interferometer [28]. The two-dimensional transfer matrix representation of double-slit-diffraction validates the formalism and allows us to extend to higher-dimensional system and find a transfer matrix representation for the same. Here we show such an application by finding the transfer matrix for a triple-slit system, demonstrating the way to extend the formalism from two slits to a higher number of slits.

2 Background

The transfer matrix representation of diffraction presented in this paper uses concepts from classical optics (Helmholtz equation), signal processing (wavelets) as well as linear algebra. A brief discussion of these concepts and their relevance in this work is presented in this section.

2.1 The Helmholtz equation and Hilbert space

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 three-dimensional (3D) fields on two-dimensional (2D) planes, say the $xy$ plane also form a Hilbert space. It should be noted that there is a difference between two-dimensional 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 $\mathrm{e}^{\mathrm{i}\omega t}$, where $\omega$ is the angular frequency. Consequently, the wave equation reduces to the time-independent Helmholtz equation [20, 21, 29]. For a source $\rho(\bm{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,

[TABLE]

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

We apply the Fraunhofer approximation to Eq. (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 Eq. (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 two-dimensional functions with compact support, that spans the Hilbert space on a slice.

2.2 Haar wavelets

We use Haar wavelets and scaling function [30, 31] to construct orthonormal basis for a slice (projections of fields on a plane, see §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 two-dimensional 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, 31, 32], which is defined by its wavelet function $\psi$ (or mother function) and a scaling function $\phi$ (or father function)

[TABLE]

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

[TABLE]

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 Figs. 2 and 2 respectively.

The Haar wavelets are orthogonal with respect to the inner product

[TABLE]

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

[TABLE]

at a particular scale, say, $j=j_{0}$ (see Fig. 2). In signal processing, any signal $f(x)$ can be decomposed into Haar wavelets and scaling function of a particular scale $j_{0}$ [30, 31, 32] as

[TABLE]

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

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.

2.3 Beam splitter and its transfer matrix

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 fibre 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 two-dimensional column vectors in a Hilbert space. In such a representation, the beam splitter transformation has a two-dimensional transfer matrix representation [5, 6, 33], denoted here by

[TABLE]

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\times 2$ transfer matrix as in Eq. (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.

2.4 Correlated outputs of 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, 37, 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.

3 Approach and method

Here we discuss the approach towards the transfer matrix formalism of diffraction using a double-slit setup as an example, and then extend its application to find the transfer matrix of a triple-slit setup. We discuss the slice modes using Haar functions and the column vector representation of the slices. Then we truncate the dimensionality by choosing certain Haar functions as post-selected in input/output modes, and finding an effective $2\times 2$ transfer matrix for double-slit-diffraction, showing that it behaves like a lossy beam splitter. We verify the efficacy of the double-slit beam splitter by studying the correlation of the outputs and, by making an MZI by concatenating two double-slit beam splitters. Finally, we use the transfer matrix formalism to find the transfer matrix of a triple-slit system demonstrating its application to higher-dimensional systems.

3.1 The slice modes

As discussed in §2.1, non-local eigenfunctions of the Helmholtz equation do not make suitable modes for detectors with finite-sized windows. Haar functions (§2.2) on the other hand, have compact support over a given interval and therefore two-dimensional 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 [31, 30].

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

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.,

[TABLE]

where the dilation parameter of the scaling function, i.e., $j_{0}$ is set so that $g_{1}(x,y;z_{1},j_{0},k,k^{\prime})$ covers the entire patch (see Eq. (5)), and $z=z_{1}$ 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.,

[TABLE]

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

[TABLE]

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 Eq. (8).

3.2 The double-slit setup

We elaborate on the slice modes concept using a customized double-slit setup with two sources and two detectors as shown in Fig. 4, where $\hat{\bm{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.

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

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

[TABLE]

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 Eq. (16) is valid only because the slit plane is far enough from the plane at $z_{1}$, so that the surface-effects are negligible.

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

Another pair of square-shaped slits $\mathrm{D_{1}}$ and $\mathrm{D_{2}}$, aligned with $\mathrm{A_{1}}$ and $\mathrm{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 $\lambda$. 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_{1}=-0.9L$ and $z_{2}=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 §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 Eq. (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.

3.3 The post-selected slice modes

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_{0}=-2$ so that the Haar scaling function covers the entire patch. According to the positions of the detectors in Fig. 4, the square patch occupied by the detector at port D1 is the one with indices $k=2,\leavevmode\nobreak\ k^{\prime}=0$. Similarly, the patch covered by the detector at port D2 is the one with $k=-3,\leavevmode\nobreak\ 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 D1 and D2) are resolved in terms of the corresponding Haar functions. For example, the field $E^{(1)}_{1}(x,y;z_{2})$ can be expanded as

[TABLE]

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 Figs. 6 and 6.

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_{1}(x,y;z_{2},-2,2,0)$. The proof of this fact is in Table 1, where the total power intercepted by port D1 is compared with the total power in the projection on $g_{1}(x,y;z_{2},-2,2,0)$. The calculation shows that about $99.9956\%$ of the total power is in the said projection. The field $E^{(2)}_{1}(x,y;z_{2})$ also has most of the power in this mode. Therefore, we ignore all the other Haar functions which have negligible contribution to the field. Similarly, for detector port D2, the dominant contribution is from $g_{1}(x,y;z_{2},-2,-3,0)$. Henceforth, two post-selected output modes are

[TABLE]

which reduce the infinite-dimensional representation of the slice at $z=z_{2}$ to a two-dimensional column vector

[TABLE]

where $E^{(i)}_{j}(x,y;z_{2})$ is the field from source Si diffracted by slit Aj and therefore intercepted by the detector port Dj, and

[TABLE]

and similarly for the second entry in the column vector.

3.3.2 Input modes

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

[TABLE]

with $z_{1}$ denoting that the modes are for a slice on plane $z=z_{1}$. 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 two-dimensional vector representation of the input, i.e., slice at $z=z_{1}$

[TABLE]

where the superscript denotes that source Si is turned on.

3.4 The effective $2\times 2$ transfer matrix

The transfer matrix $T(z_{2},z_{1})$ must map the input vector $\bm{X}^{(i)}(z_{1})$ to the output vector $\bm{Y}^{(i)}(z_{2})$, i.e.,

[TABLE]

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

[TABLE]

where $\left(\begin{matrix}\bm{X}^{(1)}(z_{1})&\bm{X}^{(2)}(z_{1})\end{matrix}\right)$ is a $2\times 2$ matrix with $\bm{X}^{(1)}(z_{1})$ and $\bm{X}^{(2)}(z_{1})$ as columns, and similar for the right-hand side. Inverting Eq. (26) yields

[TABLE]

provided that $\left(\begin{matrix}\bm{X}^{(1)}(z_{1})&\bm{X}^{(2)}(z_{1})\end{matrix}\right)$ is invertible. In general, $\left(\begin{matrix}\bm{X}^{(1)}(z_{1})&\bm{X}^{(2)}(z_{1})\end{matrix}\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.

3.5 Verifying the double-slit beam splitter

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 Mach-Zehnder interferometer.

3.5.1 Cross-correlation of the post-selected output fields

Let the field from source $\mathrm{S}_{2}$ have a phase $\varphi$ with respect to that from source $\mathrm{S}_{1}$. Also, we rotate the polarization of source S2 so that $\theta$ is the angle between the directions of polarizations of the fields from the two sources. When both the sources are turned on, port $\mathrm{D}_{1}$ intercepts the vector superposition of fields from sources $\mathrm{S}_{1}$ and $\mathrm{S}_{2}$ through slit $\mathrm{A}_{1}$. The integrated intensity in the post-selected mode on port $\mathrm{D}_{1}$, i.e., $e_{1}(x,y;z_{2})$ is the projection of the vector sum of the fields intercepted by the port, i.e.,

[TABLE]

where $\theta$ is the distinguishability parameter between the two sources. Similarly, the integrated intensity in the post-selected mode on port $\mathrm{D}_{2}$ is calculated by projecting $E_{-}(x,y;z_{2},\theta,\varphi)$ on to $e_{2}(x,y;z_{2})$.

The cross-correlation between the two outputs as a function of $\theta$ in Eq. (78) is modified for time-independent fields as

[TABLE]

which is the intensity-intensity correlation of the two output modes of the slice at $z=z_{2}$.

3.5.2 Effective Mach-Zehnder interferometer

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

The detectors behind ports $\mathrm{D}_{1}$ and $\mathrm{D}_{2}$ as in Fig. 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 $\mathrm{A}^{\prime}_{1}$ and $\mathrm{A}^{\prime}_{2}$ and reach the output ports $\mathrm{D}^{\prime}_{1}$ and $\mathrm{D}^{\prime}_{2}$, 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 $\mathrm{D}^{\prime}_{1}$ and $\mathrm{D}^{\prime}_{2}$ as the phase, say $\alpha$ 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 Eq. (9) and the phase in one arm, say $\alpha$ is set to zero (the arm lengths are considered equal), the transfer matrix for the MZI is [28]

[TABLE]

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_{4}$ (as shown in Fig. 7).

3.6 Extending to three dimensions

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, Fig. 8 shows a triple-slit setup in which a third source S3 is placed at $(-3d/2,0,-L)$, a slit A3 centered at $(-3d/2,0,0)$ and a detector D3 centered at $(-3d/2,0,L)$. Similar to the double-slit case, the transfer matrix approach can be applied to this system.

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

[TABLE]

and similarly for slice at $z_{2}$.

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

[TABLE]

where the superscripts denote the source that is turned on.

4 Results

Numerical calculations of the solutions of the Helmholtz equation (with Fraunhofer approximations, see appendix A) gives the resultant transfer matrices for the double and triple-slit setups. We also show the results of the variation of the correlation of the post-selected outputs of the double-slit setup. Further, we show the numerically calculated interference pattern at the two outputs of the MZI made by concatenating two double-slit beam splitters, and also find its transfer matrix.

4.1 Effective transfer matrix for the double-slit setup

As discussed in §3.4, the effective $2\times 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 S1 is on are

[TABLE]

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

[TABLE]

We use the above results in Eq. (27) to calculate the effective transfer matrix.

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

[TABLE]

which is a symmetric matrix as expected from the symmetry of the double-slit setup (Fig. 4). Note that apart from a factor of about $\mathrm{e}^{0.476\pi\mathrm{i}}\times 2\sqrt{2}\times 10^{-3}$, the matrix $T(z_{2},z_{1})$ 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.

4.2 The underlying unitary transformation

As expected, the transfer matrix $T(z_{2},z_{1})$ is not unitary as can be seen from

[TABLE]

because of reasons discussed in §3.4. To reveal the underlying unitary transformation a polar decomposition of $T(z_{2},z_{1})$ is carried out which yields,

[TABLE]

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

[TABLE]

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_{2},z_{1})$ 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 Fig. 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).

4.3 Cross-correlation of the post-selected outputs

Here we show the result of the numerically calculated intensity-intensity cross-correlation of the post-selected outputs using Eq. (29) with

[TABLE]

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

[TABLE]

the visibility of which is $0.5$.

On the other hand, for a 50:50 cubic beam splitter (with transfer matrix as in Eq. (9)), if the relative phase $\varphi$ between the two sources is distributed uniformly over the interval $[0,2\pi)$, the cross-correlation function shows a visibility of $0.5$ as shown in Fig. 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 behaviour of the double-slit setup.

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

4.4 The effective MZI

Here we show the interference pattern at the output of the double-slit based MZI as shown in Fig. 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 $\alpha$ yields an interference pattern shown in Fig. 11. Such an interference pattern is a signature of an MZI.

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

[TABLE]

respectively. The visibility $V_{\mathrm{MZI}}$ of both the curves is

[TABLE]

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 Fig. 7 should yield the transformation in Eq. (30).

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

[TABLE]

which after polar decomposition yields

[TABLE]

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 Eqs. (30) and (50).

4.5 Transfer matrix for the triple-slit setup

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\times 3$ transfer matrix for the triple-slit setup in Fig. 8, with respect to the post-selected modes discussed in §3.6 is

[TABLE]

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

[TABLE]

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 orbital angular momentum of light, because of the practical limits on obtaining high OAM states [7].

5 Conclusion

A post-selected unitary representation of slit-diffraction is achieved by projecting the solutions of the Helmholtz equation on two-dimensional plane and finding a transfer matrix that maps one slice to another. The Haar wavelets and scaling functions are used as orthonormal modes that span the slice on each plane. From the infinite set of modes, two input and two output modes are post-selected depending on the area and position of the detectors. The non-unitary transfer matrix is polar decomposed to reveal the underlying unitary transformation and a Hermitian component that captures the losses. Using this approach, a double-slit setup, with appropriate modification, is effectively a 54:46 beam splitter.

The beam splitter behaviour is verified by calculating the intensity-intensity cross-correlation of the outputs and getting a Hong-Ou-Mandel like variation. Two such double-slit beam splitters are concatenated to construct a Mach-Zehnder interferometer showing that sophisticated interferometers can be constructed using slit based diffraction. Similar to a double-slit setup, a higher number of slits can be used to construct a multi-input-multi-output devices, an example of which is shown by finding the post-selected transfer matrix for a triple-slit setup. The future work involves quantizing the fields and making a quantum version of slit-diffraction-based interferometers, which can be used for implementing QIP protocols.

SS would like to thank A. Sinha for valuable discussions and suggestions; S.N. Sahoo and A. Singh for discussions and proof-reading of the manuscript. SS also thanks the Canadian Queen Elizabeth II Diamond Jubilee Scholarships program (QES) for funding his visit to the University of Calgary during the course of this project. BCS appreciates the VAJRA fellowship support from SERB, Govt. of India.

\appendixpage

Appendix A Fraunhofer approximation to solutions of Helmholtz equation

Here we derive the solutions of the Helmholtz equation when only source S1 is switched on in the double-slit setup in Fig. 4, i.e.,

[TABLE]

with the boundary conditions specified in appendix B. Here we have ignored the $10^{5}$ factor that is the amplitude of the source, as it will only re-scale the diffracted field by that factor. The first subsection considers the far-field approximation and the second subsection considers the slit width to be very small compared to the distance between the center of the slit and point at which the field is calculated.

A.1 Far-field approximations applied to the propagator

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

[TABLE]

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

[TABLE]

and zero outside the area of the slit. Here

[TABLE]

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 Eq. (59) to calculate the diffracted field as

[TABLE]

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

[TABLE]

which means that $\nabla_{1}G(\bm{r}_{2},\bm{r}_{1})=-\nabla_{2}G(\bm{r}_{2},\bm{r}_{1})$ and from Eq. (61) $G(\bm{r}_{2},\bm{r}_{1})=G(\bm{r}_{1},\bm{r}_{2})$. The far-field is applied by assuming that $|\mathrm{i}k|>>\frac{1}{|\bm{r}_{2}-\bm{r}_{1}|}$. Then

[TABLE]

With this approximation, Eq. (62) simplifies to

[TABLE]

A.2 Small slit approximation

The expression for the diffracted field simplifies further when the small-slit approximation is considered. Let $\bm{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 $\bm{r}^{\prime}=\bm{R}_{s}+\bm{\Delta}_{s}$. Therefore,

[TABLE]

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

[TABLE]

Furthermore, for far-field and small slit

[TABLE]

which results in

[TABLE]

Similarly,

[TABLE]

and using this approximation in Eq. (65) we get

[TABLE]

For a rectangular slit $w_{x}$ wide along $x$ and $w_{y}$ along $y$, the integral is a $\mathrm{sinc}$ function and the final expression for the field of the slice is

[TABLE]

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

Appendix B Slice-to-slice map

The surface term of the formal solution of the Helmholtz equation is used to find the slice-to-slice map as follows. Consider the double-slit setup shown in Fig. 4 with only source S1 switched on. The Helmholtz equation for the field within this boundary is Eq. (58). As the slits are far from the source, the surface term of the formal solution of the Helmholtz equation in Eq. (1) can be dropped and the field on the plane $z=z_{1}$ is approximately $G(\bm{r},\bm{r}_{\mathrm{S}_{1}})$ (discussed in the main text in Eq. (16)) where $\bm{r}\cdot z=z_{1}$, i.e., the field is projected on this plane and hence is a slice of the 3D solution.

The field at another plane, say $z=z_{2}$ can be calculated from the slice at $z_{1}$. For this, consider a semi-infinite volume enclosed by the planes $z=z_{1}$, $z\rightarrow\infty$, $x\rightarrow-\infty$, $x\rightarrow\infty$, $y\rightarrow-\infty$ and $y\rightarrow\infty$. This volume includes the double-slit, the source is now excluded, leading to a homogeneous Helmholtz equation within the volume, albeit with complications due to the presence of the slits, whose opaque parts will have some dielectric constant other than one. Because of this the Green’s function, say $\tilde{G}(\bm{r},\bm{r}^{\prime})$ is no longer the free-space Green’s function, near the slit plane.

However, this approach yields the slice-to-slice map directly, as the solution at any point within the volume will have contribution only from the slice at $z_{1}$ because the other surfaces are at infinity. Therefore if one defines a propagator

[TABLE]

the field within the volume can be calculated from the slice at $z_{1}$ as

[TABLE]

where the integration is over the plane $z=z_{1}$. Finally, the slice at $z_{2}$ is the projection of the field on the plane $z=z_{2}$, i.e., $\left.E(\bm{r})\right|_{z_{2}}$.

The basis constructed using Haar scaling functions and the wavelet functions form a discrete orthonormal basis for the slices. This facilitates a matrix representation of the propagator that maps one slice to another. For example, the matrix representation of the free slice-to-slice propagator $\mathcal{P}(\bm{r}_{\perp},\bm{r}^{\prime}_{\perp};z,z^{\prime})$ in the new basis is

[TABLE]

where $g_{i}(\bm{r};z)$ is a basis function on slice at $z$ and $g_{j}(\bm{r}^{\prime};z^{\prime})$ is that on slice at $z^{\prime}$. Note that the basis is infinite-dimensional and therefore so is the matrix representation of the free propagator.

Appendix C The cross-correlation of outputs in 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 the outputs, a normalized version of which is

[TABLE]

where $E_{+}$ and $E_{-}$ are the output pulses when the input pulses have a time delay $\tau$ between them and a relative phase $\varphi$. The phase $\varphi$ fluctuates with a probability distribution $p(\varphi)$ and, $T_{\text{on}}$ and $T_{\text{off}}$ are detector on and off times respectively. The delay $\tau$ plays the role of a distinguishability parameter between the two input pulses. The cross-correlation is a measure of the fourth-order interference between the two outputs. For a 50:50 beam splitter, $C(\tau)$ shows a variation dependent on the shape of the pulse. If the probability distribution $p(\varphi)$ is uniform over the interval $\left[0,2\pi\right)$, i.e.,

[TABLE]

the curve shows a visibility of $0.5$. A detailed analysis of this cross-correlation with classical pulses is discussed in [40].

If the distinguishability parameter is the angle of polarization $\theta$ between the two input pulses instead of time-delay $\tau$, the cross-correlation of the intensities can be redefined as

[TABLE]

where the time-delay between the input pulses is zero. For a 50:50 beam splitter, $C(\theta)$ shows a sinusoidal variation as shown in Fig. 10, for uniformly randomized phase as in Eq. (77).

The equation that fits the data in the plot of Fig. 10 is

[TABLE]

which has a visibility of $0.5$.

As $C(\tau)$ and $C(\theta)$ are dependent on the probability distribution of the phase $\varphi$, the visibility of the curves can exceed $0.5$, with an appropriate choice of the probability distribution. For some distribution, the visibility can reach $1$, classically [40]. The variation of the correlation as a function of the distinguishability parameter is used as a signature of a beam splitter, which the double-slit setup, as is discussed in this work, also exhibits.

Appendix D The implementation of the phase shifter in MZI

The phase shifter in the double-slit MZI is modelled as a medium of thickness $t$ and with refractive index $n$. Within the medium the Green’s function (and hence the propagator) will change to

[TABLE]

where the refractive index of the medium causes a change in the propagation constant resulting in bending of light and a change in the phase. The thickness of the medium is small enough that the effect can be approximated by an extra phase

[TABLE]

imparted to the field and the net effect is captured by simply multiplying the output at port D2 by $\mathrm{e}^{\mathrm{i}\alpha}$.

Appendix E 100% dip in correlation of the outputs of the double-slit beam splitter

The visibility of the correlation is dependent on the probability distribution $p(\varphi)$. In particular, if

[TABLE]

the correlation function shows a visibility of $1.0$, as shown in Fig. 13. The function that fits the result is

[TABLE]

which shows that although the fields are classical, a 100% dip or a Hong-Ou-Mandel like dip is achieved if the probability distribution of the relative phase between the inputs are chosen carefully. A study of such an effect is discussed in [40].