Quantum Path Computing and Communications with Fourier Optics

Multi-plane diffraction (MPD) systems with classical sources and conventional intensity detection are recently proposed for scalable quantum computing (QC) and communications (QComm) with time domain entanglement resources and by exploiting the energy efficient interference of exponentially increasing number of propagation paths. MPD provides unique advantages for the challenges of scalability of qubits and complex set-ups including single photon generation and detection mechanisms in state-of-the-art linear optics implementations. However, MPD based QC architectures denoted by quantum path computing (QPC) are theoretically modeled for only electron based system set-up with Gaussian sources while proposed classical communication architectures are defined for free space propagation without modeling for arbitrary Fourier optical set-ups being mathematically equivalent to linear canonical transforms (LCTs). In this article, MPD architectures are defined, theoretically modeled and numerically analyzed for Fourier optics with arbitrary LCTs between diffraction planes while utilizing both Gaussian and Hermite-Gaussian laser modes. Photonic MPD proposes QC and QComm based on the mature science of Fourier optics significantly developed since the last century with globally available resources for fast and wide spread development of the proposed design. It promises the simplest linear optical system with important QC applications while promising novel resources for classical and quantum communications.


I. INTRODUCTION
Multi-plane diffraction (MPD) based quantum computing (QC) and communication (QComm) systems are recently introduced in [1]- [3] as the simplest and scalable linear optical architectures with significantly important advantages combining the utilization of the classical sources, simplicity of the set-up composed of only planes with diffraction slits and detection with conventional intensity photodetectors. MPD proposes a promising alternative architecture coping with the fundamental challenges of scalability of multi-photon entanglement resources and the requirements of single photon sources and detection mechanisms in conventional linear optical QC [4]- [7]. Exponentially increasing number of propagation paths as a unique form of non-Gaussian transformation of the classical source provides QC, i.e., denoted as quantum path computing (QPC) in [1], and time domain entanglement resources for QComm by utilizing tensor product Hilbert spaces in time domain [2]. Besides that, diffractive and phase space optics are also getting attention with periodic single plane diffraction (SPD) structures for implementing quantum logic gates using quantum Talbot effect [8], to test D-dimensional (qudit) Bell inequality with free space entangled quantum carpets [9], and evaluating of entanglement over the entire transverse field distribution of the photons [10] while without any discussion regarding the MPD based advantages. MPD not only promises the solutions of important number theoretical problems such as partial sum of Riemann theta function and Diophantine approximation as a QC application [1] but also brings a novel form of quantum spatial modulation (QSM) exploiting the quantum properties of light including large Hilbert space and unique interference pattern in [3] for classical optical communications. However, MPD interference pattern is formulated for electron based set-up and Gaussian sources in [1], [2] while defined only for free space propagation (FSP) of light in [3] without generalizing to arbitrary set-ups of Fourier optics, i.e., first order centered optical or quadratic-phase systems including arbitrary sections of free space, thin lenses, graded index media and spatial filters [11]. Quadratic-phase systems are mathematically characterized as linear canonical transforms (LCTs) as discussed in Section III [12].
In this article, MPD set-up is defined, theoretically modeled and numerically analyzed for Fourier optical systems with arbitrary LCTs between diffraction planes. The quantum nature of Fourier optics is discussed based on the recent experimental [13] and theoretical [14] studies verifying the validity of Fresnel diffraction formulation for quantum optical propagation. LCTs are linear integral transforms including the Fresnel and fractional Fourier transform (FRFT), scaling, chirp multiplication and some other operations as special cases while being equivalent to spatial distribution of light in phase-space optics for quadratic-phase systems [11]. Classical monochromatic light sources of both Gaussian and Hermite-Gaussian (HG) beams are utilized compatible with the standard laser sources within the paraxial approximation [15]. LCT based design which provides more flexibility is numerically compared with FSP in terms of improvement on the detection efficiency, i.e., probability of the detection or intensity at the sensor plane, and the interference complexity defined with the magnitudes of the interfering paths and negative volume of Wigner distribution function [1], [16]. It is observed that photonic MPD with Gaussian sources results in unique mathematical forms of intensity distribution on the sensor plane in (18) which are promising to be exploited for the solutions of the partial sum of Riemann theta function [17]- [20], period finding [21] and Diophantine approximation [22] similar to the algorithms and methods in [1] but with much more design flexibility due to LCTs, diversity of the tools and maturity in Fourier optics. Furthermore, HG sources result in different forms in (22) and (24) while closely related to the standard Riemann theta form and requiring future studies to exploit for the solutions of numerical problems in various scientific disciplines.
Proposed theoretical modeling and system design with widely available optical components, e.g., thin lenses, free space and diffraction planes as a form of spatial filtering, provide a unique opportunity to exploit conventional Fourier optics for QC and QComm purposes. The large amount of theoretical and experimental maturity in Fourier optics since the last century is combined with MPD based system design to realize scalable and low complexity QC, QComm and classical communications systems with important capabilities and global resources for efficient implementation and development.
The remainder of the paper is organized as follows. In Section II, MPD based design and advantages for QC and QComm are briefly reviewed. In Section III, quantum nature of Fourier optics and modeling of light propagation are discussed. Then, in Section IV, theoretical modeling of MPD with arbitrary LCT set-ups is presented for Gaussian and HG sources. The negativity of Wigner distribution function and the path magnitudes are discussed in Section V while numerical analysis is provided in Section VI. Finally, in Sections VII and VIII, open issues and conclusions are given.  [2] and (c) a scenario of interference in time where the probability of the particle to be diffracted through the first and the second planes, i.e., P E (2) and P E (3), increases with the superposition of the paths or histories with indices 4 and 5 while decreasing through the third plane, i.e., P E (4) on the fourth (detector) plane, due to the destructive interference between the two paths [2]. (d) The constellation diagram and QSM block diagram where blocking varying combinations of slits results in the superposition of varying paths which are utilized as different QSM symbols [3].

II. QUANTUM PATH COMPUTING AND COMMUNICATIONS WITH CLASSICAL SOURCES
MPD set-up is introduced and the applications for QC, QComm and classical communications proposed in [1]- [3] are briefly reviewed in the following while emphasizing the QPC capability.
The set-up is composed of N−1 diffraction planes with K j slits on each plane for j ∈ [1, N−1] and a single sensor plane indexed with N while the center of the slits is given by X j,i for i ∈ [1, K j ] as shown in Fig. 1(a). Each slit is assumed to apply a spatial filtering of exp − (x j − X j,i ) 2 / (2 β 2 j,i ) where β j,i determines the slit width. The wave function on jth plane is denoted with Ψ j (x j ) which is the wave form after diffraction through the previous planes, i.e., with the indices k ∈ [1, j − 1], while before diffraction through the slits on jth plane. There is an exponentially increasing number of propagation paths through the slits until to the final sensor plane, i.e., N p ≡ N−1 j=1 K j , while nth path includes the diffraction through a single slit on each plane with the corresponding wave function Ψ j,n (x j ) on jth plane. Assume that nth path passes through the slit indexed with s n, j on jth plane and define the path vectors ì x T N−1,n ≡ X 1,s n,1 X 2,s n,2 . . . X N−1,s n, N−1 and ì The mapping between the path index n and slit index s n, j for the path is defined with the function n = f s2n (ì s) where s n, j is predefined for each n. Furthermore, ì 0 k is the column vector of length k with all zeros and 0 k is the square matrix of all zeros with the size k × k. Similarly, rectangular matrices are shown with 0 k,l . In the rest of the article, a parameter B depending only on ì β N−1,n ≡ [β 1,s n,1 ... β N−1,s n, N−1 ] but not on ì x N−1,n is denoted with B j,n on each jth plane including . over the symbol. Therefore, if the slits are chosen with the same β j,n = β j specific to each plane, then B j,n becomes independent of n and is converted to the notation B j .
MPD set-up shown in Fig. 1(a) is utilized for QC denoted by QPC in [1] with exponentially increasing number of interfering paths as quantum resources. Consistent and entangled histories approaches are utilized in [2] to show Leggett-Garg inequality violation as the temporal analogue of Bell inequality requiring further studies to exploit the proposed time domain entanglement resources for practical QComm. An example including seven different histories of the propagation composed of diffraction through slits and measurements on the planes is shown in Fig. 1(b). A scenario showing the interference between different path histories is simulated in [2] as shown in two paths allows an increase in the probability of the particle to be diffracted through the first two planes, i.e., the probability of P E ( j) to be detected on the planes with the indices j = 2 and j = 3, compared with the fifth path alone while the destructive interference reduces the probability to be diffracted through the third plane as a counter-intuitive observation. Furthermore, a novel form of QSM is defined in [3] where opening and closing slits result in varying quantum superposition waveform message symbols at the transmitter which can be detected at the receiver as shown in Fig. 1(c). Interference of exponentially increasing number of paths is exploited which is not possible in classical domain. As a result, the capability to theoretically characterize MPD with quantum Fourier optics provides future applications of QC, QComm, quantum information theory and classical communications by exploiting energy efficient combination of optical elements.
MPD set-up measurement on the sensor plane in [1], [2] generates a black-box (BB) function with the following special form to utilize in solutions of important and classically hard number theoretical problems: where k ∈ Z, T s ∈ R + is a sampling interval, A ì s ∈ R − , B ì s ∈ R + and Υ ì s ∈ C. The complex valued matrix H ì s ≡ H R,ì s + ı H I,ì s and the vector ì h ì s ≡ ì c ì s + ı ì d ì s have the values depending on β j,i for j ∈ [1, N − 1] and i ∈ [1, K j ] corresponding to the specific selection of slits in the path ì s, inter-plane durations for the particle propagation, particle mass m (for electron based set-ups in [1], [2]), beam width σ 0 of the Gaussian source wave packet and Planck's constant . The calculation of (1) in an efficient manner is significantly hard while two different methods utilizing (1) for practical problems are introduced. Solution for the partial sum of Riemann theta or multidimensional theta function is the first application with importance in number theory and geometry [17]- [20]. The second method utilizes MPD with the phase ì d T ì s ì x N−1,ì s in exp ì h T ì s ì x N−1,ì s k T s for period finding [21] and the solution of specific instances of SDA problems [22].
QPC based on Fourier optics promises expanding the set of problems both with LCT based general system design and also the sources including HG beams. Furthermore, a discussion is included to utilize non-Gaussian slits with the proposed mathematical modeling in Section VII.
Next, the propagation through Fourier optical systems based on Fresnel diffraction is modeled emphasizing the quantum nature of Fresnel diffraction and Fourier optics.

III. QUANTUM FOURIER OPTICS
In scalar diffraction theory, the first Rayleigh-Sommerfeld formula of the Huygens-Fresnel principle for the propagation of light on planar surfaces is described as follows by using the Green's theorem [23]: where U I (P 1 ) is the wave amplitude at the point P 1 , U is the distribution on the planar screen where diffraction occurs, Σ denotes the integration over the slit including its multiplicative 2 [23]. Assuming that r 01 ≫ λ, the following approximation holds in rectangular coordinates: where the kernel K FS (ì r 1 , ì r 0 ) for FSP is defined as follows: where The kernel for Fresnel diffraction integral is obtained by further approximation of r 01 in the near-field for large z resulting in (6). This expression is the convolution integral conventionally used in phase-space optics for FSP [11].
Recently, scalar diffraction theory and Fresnel diffraction integral are discussed in [13] to be validly representing the evolution of light wave function modeled with the Hamiltonian of the quantized electromagnetic field H = (p 2 + ω 2q2 ) / 2 as the Feynman's path integral (FPI) solution of the quantum mechanical harmonic oscillator (HO) [24]. Fresnel diffraction nature of the propagation is verified with experimental photon counting studies for single photons. The wave function amplitude of light field in one dimension (1D) on a plane Ψ(x 0 ) is modeled to propagate into the amplitude Ψ(x 1 ) on another plane (Eq. 16 in [13] transformed into a simpler form) with the following formulation: where the kernel based on HO is the following: where c is the velocity of light, ω ≡ 2 π c / λ, ω t n π for n ∈ Z, t 01 is the propagation duration between the planes and m λ ≡ k / c is the defined equivalent mass of photon propagation. In addition, the approximated FSP kernel in (6) is simply converted to the following in 1D: The kernel for massive particles with the mass m such as an electron is expressed as follows [1], [24]: In other words, the formulation based on phase-space optics for photon and electron propagation wave amplitudes have the similar form in (9) and (10) except an overall phase factor. The form in (9) is utilized in [3] for defining QSM while targeting only classical communications.
On the other hand, both the kernels K HO (x 1 , x 0 ) and K FS (x 1 , x 0 ) are special cases of LCTs defined for quadratic-phase optics [11]. As a class of linear integral transforms, they include as special cases the Fresnel and fractional Fourier transform (FRFT), simple scaling, chirp multiplication and some other operations. Spatial distribution of light in phase-space optics for the class denoted by quadratic-phase systems is mathematically equivalent to LCTs (Chapters 3 and 8 in [11]). These optical systems include arbitrary combinations of the sections of free space in the Fresnel approximation, thin lenses and sections of quadratic graded-index media.
In [13], FRFT nature of the kernel K HO (x 1 , x 0 ) is shown both theoretically and experimentally while emphasizing the applicability of all the properties of Fourier analysis to quantum optics. In this article, propagation of the wave function is extended to the general case of LCTs providing flexibility to utilize arbitrary optical set-ups by enlarging the functional structures and number theoretical problems exploited in QPC. Furthermore, a better control is obtained for the energy flow of the light through the slits. The kernel of 1D quadratic-phase system or LCT is represented as follows: In matrix notation, it is shown with the following unit-determinant matrix: where a d − b c = 1 and the matrix for the composition of two consecutive systems represented by M 1 and M 2 is calculated by the multiplication M 2 M 1 . The matrices for K HO (x 1 , x 0 ) and K FS (x 1 , x 0 ) are given as follows: M HO has the same form with the propagation of light in quadratic graded-index media of having where n 0 and χ are the medium parameters. The parameter matrix of the propagation through the quadratic graded-index medium of length d gri is given by the following [11]: where α = d gri / χ. There is a FRFT relation between scaled versions of the inputf (x) and output and f a (x) denotes the ath order FRFT of f (x). FRFT operation of order α is represented with the parameter matrix of a = d = cos(α) and b = sin(α). As a result, M HO represents a FRFT relation between the input and output scaled with the parameter λ χ where the parameters are α = ω t and λ χ ≡ 2 π t 01 / m λ while as a special case of LCTs. Next, MPD modeling is proposed for Fresnel diffraction and arbitrary LCT based optical systems by utilizing the proposed kernels.

IV. MULTI-PLANE DIFFRACTION WITH QUANTUM FOURIER OPTICAL SYSTEMS
The set-up in Fig. 1(a) defined with FSP and electron based theoretical formulation is extended to optical systems of LCT as shown in Fig. 3. An arbitrary LCT with the matrix parameters {a j, j+1 , b j, j+1 , c j, j+1 , d j, j+1 } is implemented in phase-space optics by consecutive applications of FSP of length L a, j, j+1 , then thin lens of focal length f j, j+1 , and another FSP of length L b, j, j+1 [12]. LCT matrix M LCT is calculated as follows: where τ ⋆ a, j ≡ L a, j, j+1 / c and τ ⋆ b, j ≡ L b, j, j+1 / c, and the middle matrix is for the effect of thin lens [11]. FRFT with scaling is a special case of LCT as discussed in Section III. Therefore, FSP, FRFT and arbitrary LCT based QPC set-ups are implemented with the universal configuration in Fig. 3.

A. Fresnel Diffraction and FRFT with Gaussian Sources
It is assumed that the optical system between the planes results in the kernels defined in (6) and (8) based on Fresnel diffraction integral for free space and FRFT based propagation denoting graded-index media propagation or HO solution, respectively. It is assumed that the source distribution has a Gaussian form of Ψ 0 ( waveforms as eigenfunctions of FRFTs [11] are considered for the general case of LCTs in Section IV-B2. The kernel K FS (x 1 , x 0 ) has the same form with K e,FS (x 1 , x 0 ) used for QPC modeling in [1] by replacing the electron mass m with the photon equivalent mass m λ . Therefore, the same formulations are utilized for both the final wave function and problem solving capabilities. In addition, iterative integration with K HO (x 1 , x 0 ) results in the final intensity distribution of MPD with the same form of K FS (x 1 , x 0 ) while with different iteration parameters as shown in Table   II in Appendix D.
The wave function for nth path on the sensor plane for the general case of non-uniform slit widths is given by the following by using the iterative formulation: It is further simplified by extraction of ì x N−1,n dependent parts and summing the contributions from each path as follows: where Υ G N,n ≡ χ 0 ξ j,n , and the complex vector ì h N−1,n and the matrix H HO/G N−1,n are defined in Appendix A while having simpler structures compared with the formulation for electron based FSP set-up in [1]. The corresponding iteration parameters are given in Table II in Appendix D.
We have not included the effects of special forms of K HO (x j+1 , x j ) with ω t j, j+1 = k π for [11]. This case can be simply realized by assuming that spatial filtering operations of the slits on jth and ( j +1)th planes are combined on a single plane by also noting that whether the wave function is reversed or not. For example, multiple inter-plane propagation intervals can result in multiple reversals with the overall effect of the identity and combined spatial filtering of Gaussian slits.  [11]. These cases further improve the flexibility of the LCT system to realize the desired transformation on the wave function.
2) Hermite-Gaussian Sources: If the source is chosen as the standard HG waveform of 0) where H l (x) ≡ (−1) l e x 2 d l e −x 2 / dx l is the lth order Hermite polynomial [11], [25], then Ψ N,n (x N ) for nth path is obtained as follows by using the integral equality of HG functions in Appendix B: where the parameters χ 01 , χ j, j+1,n , u j, j+1,n , v j, j+1,n , g j, j+1,n and h j, j+1,n obtained in an iterative manner for j ∈ [1, N − 1] are calculated with simple algebra for nth path and shown in Table III in Appendix D. Simple algebraic manipulations of (19) to extract ì x N−1,n dependent parts result in the following simplification: where χ a, j, j+1,n for j ∈ [1, N − 1] is defined in Table III (22) is transformed into the following: It is further simplified as follows by using the useful identity H l (x + y) = (H + 2 y) l in [26] where H k ≡ H k (x): where is similarly expected to be significantly hard since the mathematical form is more complicated compared with the partial sum of Riemann theta function.
The set-up parameters including the slits, lenses and inter-plane distances are required to be tuned in order to obtain the desired vectors  (18), (22) and (24) for the target number theoretical problems. Next, Wigner distribution is defined where its negative volume is regarded as an indicator of non-classicality.

V. WIGNER DISTRIBUTION, NEGATIVITY AND PATH MAGNITUDES
The momentum domain wave function Ψ p, j (p j ) is defined as Fourier transform of spatial representation of wave function Ψ j (x j ) on jth plane as follows: The distribution of energy through space-momentum phase-space is described by Wigner distribution function defined as follows [1], [16]: The negative volume of Wigner function defined in [16] and utilized in [1] to describe the increasing non-classicality or time-domain entanglement resources in [2] is described as V j ≡ ∫ ∫ |W j (x j , p j )| dx j dp j − 1 / 2. On the other hand, the probability of the particle to be detected on jth plane, i.e., to be diffracted through ( j − 1)th plane, is computed as follows:   Table I while  to adapt LCT parameters with respect to any given set-up including inter-plane distances, slit locations and widths in a way to maximize the interference and the probability of the particle reaching to the sensor plane. Observe that HG source has already negative Wigner volume of 1.076 which is much further improved by LCT set-up compared with FSP.  (18), (22) and (24)  However, any slit structure can be represented as a composition of Gaussian slits by using the method defined in [27] and applied successfully in optical diffraction theory and experiments [28], [29]. The one dimensional slit function G(x) is represented as follows: where a i and β i are found with optimization based on the experimental measurement results while increasing K provides more accurate results. If the perfect Gaussian slits are replaced with the superposition in (28), then the summations in (18), (22) and (24)  There are some factors reducing the compatibility between the theory and practice. These include imperfection in optical set-up, e.g., finite size lens effects, planar thickness, characterization of slit functions, sources and detector efficiency. The theoretical model should be extended including all the set-up parameters having diverging effects on the final intensity distribution.
Similarly, the effects of exotic paths, i.e., trajectories between the slits on the same plane, should be included in the mathematical model as thoroughly discussed in [1], [14], [30].
Another open issue is related to utilize measurements on all the sensor planes for computational purposes not only the final sensor plane since they include diffraction through previous planes. with the corresponding parameters defined in Table II in Appendix D. The elements in the vector ì h N−1,n = ì c N−1,n + ı ì d N−1,n are defined as follows: where ì v k, j,n for k ∈ [0, j − 1] is given as follows: where the matrix multiplication k i=1 U i denotes U 1 U 2 . . . U k for any matrix U i for i ∈ [1, k] and p 4, j,n , p 5, j,n , ζ j,c,n and ζ j,d,n for j ∈ [1, N − 1] are defined in Table II in Appendix D. Assume that diag{ ì y 1 , . . . , ì y K } and diag{y 1 , . . ., y K } define the operators creating block diagonal matrices by putting the vectors ì y j and the matrices y j for j ∈ [1, K], respectively, (all the vectors or the matrices having the same dimensions) to the main diagonal and making zero the remaining elements. The matrix H HO/G N−1,n is more simplified as follows compared with the more complicated form achieved for electron based FSP in [1]: where the diagonal matrices are defined as follows: where 2 × 2 block K j,b,n and 1 × 2 vector ì k T j,c,n are defined as follows for j ∈ [2, N − 1]: ì k T j,c,n = p 3, j,n 1 ı (36) V N−1,n is a lower triangular block matrix defined as follows: Expanding H HO/G N−1,n in terms of real and imaginary parts is achieved by finding the real and imaginary parts of p 1, j,n for j ∈ [1, N − 1] and p 3, j,n for j ∈ [2, N − 1], and K j,b,n and ì k j,c,n for j ∈ [2, N − 1] since V N−1,n is a real matrix. This is easily achieved by using the explicit forms of p 1, j,n and p 3, j,n in Table II in Appendix D.