Double-slit interferometry as a lossy beam splitter

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.


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 finitedimensional 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 twoinput-two-output optical device that has a 2 × 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.

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.

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 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 ρ(r) in a volume V enclosed by a surface ∂V on which the boundary condition is specified, the most general solution of the Helmholtz equation is, wheren(r ) is unit normal to the surface and G(r, r ) 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.

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 ψ (or mother function) and a scaling function φ (or father function) respectively, where is the box function. These functions are dilated and translated to create other Haar wavelets and scaling functions, and the Haar scaling functions are also orthogonal, i.e., 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 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.

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

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

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

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 which take integer values. Figure 3: An example of how a plane at some z can be segmented into non-overlapping square patches indexed by two integers k and k . The width of each patch is equal to the width of the available detector. To know the entire slice at z, all the patches must be considered. But usually, there are a finite number of detectors so that only a few patches can be covered. In that case, only those Haar wavelets and scaling functions are considered which have a compact support on the considered patches.
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., 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 ) 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., where Z + is the set of positive integers. The subscript ı is a meta-index for m, n, m and n , and 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).

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ŷ 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. Figure 4: Schematic of the double-slit setup considered in this paper. Two point sources S 1 and S 2 emanate monochromatic linearly polarized light with harmonic time-dependence. The imaginary plane at z = z 1 (represented by a dotted line) is for the input slice. Two slits A 1 and A 2 are placed at z = 0 where each slit is aligned center-to-center with one of the sources. A second pair of slits D 1 and D 2 are placed at z = z 2 where each port is aligned center-to-center with one of the source. Behind each of these slits is a square-faced detector which measures the integrated intensity of the light falling on it. The plane z = z 2 is also for projecting the output slice. A perfectly absorbing barrier runs between z = 0 and z = z 2 that prevents the field from slit A 1 (A 2 ) from reaching port D 2 (D 1 ). The dashed arrows represent the ray approximations of the fields from the sources to the detectors.
The sources S 1 and S 2 are monochromatic point-like sources (practically a spherical source with diameter ∼ λ) emanating linearly polarized light as spherical waves (Y 0 0 (θ, φ) spherical harmonic [21]) with wavelength λ. They are placed at r S1 = (d/2, 0, −L) and r S2 = (−d/2, 0, −L) 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 δ 3 (r − r S1 ) and δ 3 (r − r S2 ). 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 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 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 A 1 and A 2 , each with side-length w = 4λ are placed at z = 0 with the positions of their centers r A1 = (d/2, 0, 0) and r A2 = (−d/2, 0, 0) respectively, and therefore aligned with the respective sources. In the Fraunhofer regime, the diffracted field E (i) j (r) from source S i through slit A j , 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 D k is blocked from the field E Another pair of square-shaped slits D 1 and D 2 , aligned with A 1 and 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 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.

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.

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 D 1 is the one with indices k = 2, k = 0. Similarly, the patch covered by the detector at port D 2 is the one with k = −3, k = 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 1 and D 2 ) are resolved in terms of the corresponding Haar functions. For example, the field E (1) 1 (x, y; z 2 ) can be expanded as where A ı (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 (x, 0; z 2 ) is shown in Figs. 5 and 6.
1 (x, 0, z 2 ) using the wavelets supported over the patch k = 2, k = 0, using Eq. (17). The values of dilation parameter m for the Haar wavelets (as in Eq. (4)) is taken from −2 to 2 so that the Haar wavelets are visible. A finer reconstruction can be done by taking m upto higher values. For m upto 6, the reconstruction is almost perfect.
1 (x, 0, z 2 ) using the wavelets supported over the patch k = 2, k = 0, using Eq. (17). The values of dilation parameter m for the Haar wavelets (as in Eq. (4)) is taken from −2 to 2 so that the Haar wavelets are visible. A finer reconstruction can be done by taking m upto higher values. For m upto 6, the reconstruction is almost perfect.
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 D 1 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 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 D 2 , the dominant contribution is from g 1 (x, y; z 2 , −2, −3, 0). Henceforth, two post-selected output modes are which reduce the infinite-dimensional representation of the slice at z = z 2 to a two-dimensional column vector where E (i) j (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 and similarly for the second entry in the column vector.

Integrated intensity intercepted Integrated projections on chosen modes Ratio
2 (z 2 ) e 2 (x, y; z 2 ) 2 = 0.633087 99.9956% 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.

Input modes
Each source in the double-slit setup (Fig. 4) 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 where the superscript denotes that source S i is turned on.

The effective 2 × 2 transfer matrix
The transfer matrix T (z 2 , z 1 ) must map the input vector X (i) (z 1 ) to the output vector Y (i) (z 2 ), i.e., for i ∈ {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 X (i) (z 1 ) and Y (i) (z 2 ) for each source turned on at a time. The transformation equation for both sources are combined into one matrix equation where X (1) (z 1 ) X (2) (z 1 ) is a 2 × 2 matrix with X (1) (z 1 ) and X (2) (z 1 ) as columns, and similar for the right-hand side. Inverting Eq. (26) yields provided that X (1) (z 1 ) X (2) (z 1 ) is invertible. In general, X (1) (z 1 ) X (2) (z 1 ) 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.

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.

Cross-correlation of the post-selected output fields
Let the field from source S 2 have a phase ϕ with respect to that from source S 1 . Also, we rotate the polarization of source S 2 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 1 intercepts the vector superposition of fields from sources S 1 and S 2 through slit A 1 . The integrated intensity in the post-selected mode on port 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., 1 (z 2 ) cos ϕ cos θ, where θ is the distinguishability parameter between the two sources. Similarly, the integrated intensity in the post-selected mode on port D 2 is calculated by projecting E − (x, y; z 2 , θ, ϕ) on to e 2 (x, y; z 2 ). The cross-correlation between the two outputs as a function of θ in Eq. (75) is modified for timeindependent fields as which is the intensity-intensity correlation of the two output modes of the slice at z = z 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.
Phase shifter y z = 0 Figure 7: Using two double-slit setups a Mach-Zehnder interferometer is constructed by concatenating them such that both are aligned center-to-center and parallel to each other. The output ports D 1 and D 2 of the first double-slit beam splitter serve as the inputs for the second one. The two detectors are placed behind the output ports D 1 and D 2 of the second double-slit setup. A phase shifter is placed at the output port D 2 which changes the phase of the field from that port by α.
The detectors behind ports D 1 and 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 A 1 and A 2 and reach the output ports D 1 and D 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 D 1 and D 2 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 Eq. (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] 1 2 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).

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 S 3 is placed at (−3d/2, 0, −L), a slit A 3 centered at (−3d/2, 0, 0) and a detector D 3 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 e 1 (x, y; z 1 ) :=g 1 (x, y; z 1 , −2, 2, 0), 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 where the superscripts denote the source that is turned on. Figure 8: Schematic of a triple-slit setup constructed in a similar way as the double-slit setup in Fig. 4, by adding a source, slit and detector to the latter.

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.

Effective transfer matrix for the double-slit setup
As discussed in §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 1 is on are with respect to the post-selected modes (see §3. 3). Similarly, the post-selected vector representations of the input and output slice when source S 2 is on are 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 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 e 0.476πi × 2 √ 2 × 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.

The underlying unitary transformation
As expected, the transfer matrix T (z 2 , z 1 ) is not unitary as can be seen from 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, 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 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 postselected 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).

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 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 the visibility of which is 0.5.  Figure 10: The cross-correlation function plotted as a function of the relative polarization angle θ between the input pulses, for a cubic beam splitter. The correlation is minimum when both the sources are indistinguishable, i.e., θ = 0, and maximum when they are completely distinguishable, i.e., θ = π/2. The detailed calculations are presented in appendix C On the other hand, for a 50:50 cubic beam splitter (with transfer matrix as in Eq. (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 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].

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 α yields an interference pattern shown in Fig. 11. Such an interference pattern is a signature of an MZI. Output at D 1 Output at D 2 Figure 11: Interference pattern at the output of the MZI made of concatenated double-slit beam splitters.
The curves that fit the resultant integrated intensities at ports D 1 and D 2 are approximately respectively. The visibility V MZI of both the curves is 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

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

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.

Appendices A Fraunhofer approximation to solutions of Helmholtz equation
Here we derive the solutions of the Helmholtz equation when only source S 1 is switched on in the double-slit setup in Fig. 4, i.e., 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. (55)) in a volume V enclosed by a surface ∂V is Consider the surface ∂V to be at infinity such that there is no contribution from the surface term in Eq. (56).
As we have assumed that the opaque portions of the slit-plane are perfectly absorbing, the field within the area of the slit is E(r) = G(r, r f ).
and zero outside the area of the slit. Here 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 ∂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. (56) to calculate the diffracted field as where the only contribution is from the slits-plane. To further simplify the expression for the diffracted field, note that which means that ∇ 1 G(r 2 , r 1 ) = −∇ 2 G(r 2 , r 1 ) and from Eq. (58) G(r 2 , r 1 ) = G(r 1 , r 2 ). The far-field is applied by assuming that |ik| >> 1 |r2−r1| . Then With this approximation, Eq. (59) simplifies to

A.2 Small slit approximation
The expression for the diffracted field simplifies further when the small-slit approximation is considered. Let 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 r = R s + ∆ s . Therefore, where the small slit approximation is applied as |∆ s |/|R s − r f | is approximately G(r, r S1 ) (discussed in the main text in Eq. (16)) where r · 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 semiinfinite volume enclosed by the planes z = z 1 , z → ∞, x → −∞, x → ∞, y → −∞ and y → ∞. 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, sayG(r, r ) 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 the field within the volume can be calculated from the slice at z 1 as 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., E(r)| z2 . 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 P(r ⊥ , r ⊥ ; z, z ) in the new basis is where g i (r; z) is a basis function on slice at z and g j (r ; z ) is that on slice at z . Note that the basis is infinite-dimensional and therefore so is the matrix representation of the free propagator.

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 where E + and E − are the output pulses when the input pulses have a time delay τ between them and a relative phase ϕ. The phase ϕ fluctuates with a probability distribution p(ϕ) and, T on and T off are detector on and off times respectively. The delay τ 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(τ ) shows a variation dependent on the shape of the pulse. If the probability distribution p(ϕ) is uniform over the interval [0, 2π), i.e., 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 θ between the two input pulses instead of time-delay τ , the cross-correlation of the intensities can be redefined as C(θ) := dϕ p (ϕ) dt |E + (t; θ, ϕ)| 2 dt |E − (t ; θ, ϕ)| 2 dϕ p (ϕ) dt |E + (t; θ, ϕ)| 2 dϕ p (ϕ ) dt |E − (t; θ, ϕ )| 2 , where the time-delay between the input pulses is zero. For a 50:50 beam splitter, C(θ) shows a sinusoidal variation as shown in Fig. 10, for uniformly randomized phase as in Eq. (74).  Figure 12: The cross-correlation function plotted as a function of the relative polarization angle θ, which is the distinguishability parameter. The correlation is minimum when both the sources are indistinguishable, i.e., θ = 0, and maximum when they are completely distinguishable, i.e., θ = π/2. The plot has been generated using Eq. (75) for a regular cubic 50:50 beam splitter with two identical input pulses having zero delay between them. When the distinguishability parameter is θ the shape of the pulses does not affect the correlation.
The equation that fits the data in the plot of Fig. 10 is C fit (θ) = 0.75 − 0.25 cos 2θ, which has a visibility of 0.5. As C(τ ) and C(θ) are dependent on the probability distribution of the phase ϕ, 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.

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 G(r, r ; n) = − 1 4π 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 imparted to the field and the net effect is captured by simply multiplying the output at port D 2 by e iα . The visibility of the correlation is dependent on the probability distribution p(ϕ). In particular, if the correlation function shows a visibility of 1.0, as shown in Fig. 13. The function that fits the result is C 100 (θ) = 0.5 − 0.5 cos 2θ, 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].