Theory of quantum path computing with Fourier optics and future applications for quantum supremacy, neural networks and nonlinear Schrödinger equations

The scalability, error correction and practical problem solving are important challenges for quantum computing (QC) as more emphasized by quantum supremacy (QS) experiments. Quantum path computing (QPC), recently introduced for linear optic based QCs as an unconventional design, targets to obtain scalability and practical problem solving. It samples the intensity from the interference of exponentially increasing number of propagation paths obtained in multi-plane diffraction (MPD) of classical particle sources. QPC exploits MPD based quantum temporal correlations of the paths and freely entangled projections at different time instants, for the first time, with the classical light source and intensity measurement while not requiring photon interactions or single photon sources and receivers. In this article, photonic QPC is defined, theoretically modeled and numerically analyzed for arbitrary Fourier optical or quadratic phase set-ups while utilizing both Gaussian and Hermite-Gaussian source laser modes. Problem solving capabilities already including partial sum of Riemann theta functions are extended. Important future applications, implementation challenges and open issues such as universal computation and quantum circuit implementations determining the scope of QC capabilities are discussed. The applications include QS experiments reaching more than \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2^{100}$$\end{document}2100 Feynman paths, quantum neuron implementations and solutions of nonlinear Schrödinger equation.

Scientific RepoRtS | (2020) 10:10968 | https://doi.org/10.1038/s41598-020-67364-0 www.nature.com/scientificreports/ Furthermore, simple optical setups exploiting wave-particle duality and interferometers have the cost of exponential complexity of resources either in time, space or energy domains to achieve QC advantages as discussed in Ref. 6 . For example, wave particle (WP) computer 29 exploits full optical interconnections of an N × N input signal array with an N × N output signal creating N 4 channels with a tensor product. WP computer also utilizes a filter array between the input and output to increase the number of connections in an additive manner with respect to the connections in each inter-planar region. Such architectures, including also Ref. 25 for a slit based modeling, provide advantages of parallelism compared with classical models without exploiting temporal correlations of quantum histories and their tensor product structure 12 . They utilize the tensor product only for a single interplanar propagation, i.e., a single measurement plane directly detecting propagation from the input array. The rich set of two and three dimensional alternative optical interconnection architectures and opto-electronic computing are discussed in detail in Ref. 30 by also including multi-stage interconnection topology. Analog Fourier optics (FO) and its digital equivalent, i.e., digital FO architectures composed of smart pixel arrays of two-dimensional electronic processing units connected with optical interconnections, exploit speed and parallelism advantages of the optical design 31 . Furthermore, programmable directed logic networks are discussed in Ref. 32 by emphasizing the energy efficiency of optical architectures.
On the other hand, QPC formulation is performed for electron based set-up in Ref. 6 while theoretical studies modeling QPE in Ref. 12 and classical optical communications in Ref. 33 formulate free space propagation (FSP) of light. They do not generalize to arbitrary set-ups of FO, i.e., first order centered optical or quadratic-phase systems including arbitrary sections of free space, thin lenses, graded index media and spatial filters 34 and mathematically characterized as linear canonical transforms (LCTs) 35 . 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 34 . Besides that, previous MPD studies utilize Gaussian sources without extending to Hermite-Gaussian (HG) beams compatible with the standard laser sources within the paraxial approximation 36 . Photonic QPC formulation is not available while important applications of QPC other than the partial sum of Riemann theta function (RTF) and period finding presented in Ref. 6 are not discussed and theoretically analyzed yet.
Diffractive and phase space optics are also getting attention in quantum technologies with periodic single plane diffraction for implementing quantum logic gates using quantum Talbot effect 37 , for testing D-dimensional (qudit) Bell inequality with free space entangled quantum carpets 38 and for the evaluation of entanglement over the entire transverse field distribution of the photons 39 while without any discussion regarding the MPD based advantages. Proposed theoretical modeling and system design of photonic QPC 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 FO for QC. The large amount of theoretical and experimental maturity in FO since the last century is combined with MPD based system design to realize scalable and low complexity QC systems with important capabilities and global resources for efficient implementation and development.
In this article, QPC set-up is defined, theoretically modeled and numerically analyzed for FO with arbitrary LCTs between diffraction planes. QPC system exploiting diffraction in an unconventional manner maintains photonic advantages including decoherence and noise while avoids the need to interact with multiple photons by eliminating many problems encountered in multi-photon entanglement and circuit implementations. Furthermore, the quantum nature of FO is discussed based on the experimental 15,40 and theoretical 41 studies verifying the validity of Fresnel diffraction formulation for quantum optical propagation. Classical monochromatic light sources of both Gaussian and Hermite-Gaussian (HG) beams are utilized compatible with the standard laser sources within the paraxial approximation 36 . LCT based design which provides more flexibility is numerically compared with FSP in terms of improvement on the detection efficiency and the interference complexity defined with the magnitudes of the interfering paths and negative volume of Wigner distribution function 6,42 .
Important future applications of photonic QPC are, for the first time, introduced and theoretically modeled in an introductory and brief manner. These include the feasibility of QS experiments compared with alternative technologies, adapting certified random number generation protocols for the photonic QPC architecture [43][44][45][46] , quantum neural network (QNN) implementations and making the solutions of nonlinear Schrödinger equation (NLSE) easier. The detailed modeling and utilization of photonic QPC for these applications are presented as open issues.
The potential of QS experiments with photonic QPC is presented in this article to reach more than 2 100 Feynman paths in a scalable set-up with several tens of diffraction planes while requiring experimental implementations for better modeling and verifying the scalability for large scale QPC set-ups. A feasible method is proposed to exponentially increase the number of Feynman paths with the cost of linearly increasing number of planes and slits allowing to obtain significantly large Hilbert space. However, it is an open issue to verify QS capability both complexity theoretically and experimentally based on the promising results in Ref. 6 and the modeling in this article such as by performing analogous modeling and experiments in Refs. 10,11 and 47 achieved for Boson sampling. Moreover, QPC with Gaussian sources results in unique mathematical forms of wave functions on the sensor plane in (16) to be exploited for the solutions of the partial sum of RTF [48][49][50][51] , period finding 52 and Diophantine approximation 53 similar to the algorithms and methods in Ref. 6 but with much more design flexibility due to LCTs, diversity of the tools and maturity in the science of FO. HG sources result in different forms in (25) and (27) while closely related to the standard RTF form and requiring future studies to exploit for the solutions of numerical problems in various scientific disciplines. On the other hand, open issues and challenges for FO based QPC design are discussed to determine the scope of the proposed design for QC purposes, e.g., universal quantum computation capability, implementations of quantum circuit gates and basic search algorithms such as Grover search.
Neural networks (NNs) exploiting the quantum advantages, i.e., QNNs, improve the capabilities of classical NNs with quantum interference and superposition for deep learning applications 54 in various disciplines 55 . On the other hand, linear and unitary framework of quantum mechanics results in the challenges of implementing non-linear and dissipative dynamics of classical neural networks 55 . The state-of-the-art neuron implementations utilize various methods to introduce non-linearity including quantum measurements 55,56 . The quantum interference among the exponentially increasing number of paths and the entanglement denoted as QPE in Ref. 12 are promising for designing and practically implementing novel design of QNNs. QPC set-up has inherently nonlinear formulation with respect to slit positions to encode the input and it operates on the quantum superposition of the inputs. Besides that, implementations of diffractive NNs utilizing single-layer 57,58 and alloptical multi-layer diffractive architectures 59 do not exploit interference among the paths or quantum domain advantages. Photonic QPC succeeds to combine the implementations of QC and QNNs with the same hardware design of MPD as a uniquely valuable unconventional hardware architecture.
NLSE solution is very important in the analysis and performance measurement of fiber optic cables 60,61 . It is also necessary for the solution of nonlinear Fourier transform (NLFT) which is a transformation that finds a wide range of applications with increasing importance 60 . NLSE and NLFT play a similar role for nonlinear and integration equations compared with the role of the FT in linear systems. NLFT is a transformation system for expressing the signal in the time plane by using nonlinear periodic waves or solitons 60 . It is also referred to as scattering transform. NLSE is expressed as follows 60,61 : where q(x, t) is the solution wave function that provides the periodic boundary condition ( q(x + l, t) = q(x, t) and period l > 0 ) and κ is some variable. In this article, the speed up in NLSE solution is conjectured by exploiting RTF summations in QPC.
The remainder of the paper is organized as follows. We firstly define and theoretically model photonic QPC and its extension for FO followed by the discussion of the performance based on Wigner distribution function. Then, future applications including QS, quantum neuron implementation and solution of NLSE are introduced, theoretically modeled and the challenges are discussed. Numerical analysis for photonic MPD is provided and then open issues for realizing photonic QPC are presented.

Results
Quantum path computing with optical multi-plane diffraction and coherent light sources. MPD set-up introduced in Ref. 6 as shown in Fig. 1 is extended to optical implementations for QC by using coherent laser sources and conventional photodetectors. 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 central position of a slit is given by X j,i for i ∈ [1, K j ] as shown in Fig. 1a. Each slit is assumed to apply a spatial filtering of G(X j,i , β j,i , x j ) ≡ exp − (x j − X j,i ) 2 / (2 β 2 j,i ) , i.e., slit mask function, 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 we 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 � x T N−1,� s ≡ X 1,s 1 X 2,s 2 . . . X N−1,s N−1 where � s ≡ s 1 s 2 . . . s N−1 and (.) T is the transpose operation. 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,s n,j = β 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 is utilized for QC denoted by QPC in Ref. 6 by sampling the interference of exponentially increasing number of interfering paths. The capability to theoretically characterize QPC with quantum FO provides future applications for both QC and quantum information theory by exploiting energy efficient combination of optical elements. Intensity sampling on the sensor plane (I[k]) for the MPD set-up in Refs. 6,12 generates a black-box (BB) function f BB [k] 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 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 Refs. 6,12 ), beam width σ 0 of the Gaussian source wave packet and Planck's constant . The calculation of (2) in an efficient manner is significantly hard while two different methods utilizing (2) for practical problems are introduced. Solution for the partial sum of RTF or multi-dimensional theta function is the first application with importance in number theory and geometry [48][49][50][51] .
The second method utilizes MPD with the phase term 52 and the solution of specific instances of SDA problems 53 .
Scientific RepoRtS | (2020) 10:10968 | https://doi.org/10.1038/s41598-020-67364-0 www.nature.com/scientificreports/ The basic unit of QC systems, i.e., the qubit, is defined on a two-state system where discretized degrees of freedom (DoF) of photons including path, transverse-spatial modes and time/frequency bins are exploited to create high-dimensional entanglement 62 . For example, multi-slit structures are already utilized to define spatial qudits by projecting the wave function into the transverse position and momentum Hilbert spaces through slits and characterizing their properties using their propagation in free space 63,64 . The qubit states in Ref. 64 are expressed in the basis |l�, |r� representing the photons passing through either the left or the right slit. However, entangled multiple photons, e.g., photon pair A and B, are conventionally generated through spontaneous parametric down-conversion (SPDC) to realize multi-photon entangled state, e.g., 1 / √ 2 |l A �|r B � + |r A �|l B � . The fundamental difference of MPD based qudits from multi-photon slit based entangled spatial qudits is the utilization of the tensor product structure for each single photon in time domain rather than spatially among multiple photons obtained through SPDC 6 . The projection events through the slits of consecutive planes are freely entangled at two different time instants denoted as QPE with the detailed modeling presented in Ref. 12 based on consistent histories and entangled histories frameworks. The presented free entanglement in time domain provides an important advantage exploiting directly the classical light sources and not requiring the difficult coupling of multiple photons. The concept of free entanglement is introduced for boson sampling exploiting boson statistics of a number of indistinguishable bosons while they still require multiple photons, and generation and detection mechanisms for single photons 10 .
History state in MPD is composed of diffraction events as follows 12 : where P j,s n,j represents the projection operator through the slit indexed with s n,j , π n as 0 or 1 determines a compound set of trajectories, ⊙ denotes tensor product operation and [ρ 0 ] denotes the initial state. The analogy of MPD based multiple qubits with the general two qubit state of two photons is represented as shown in Fig. 2 in the basis of |U� and |L� for the upper and lower slits, respectively. The general state for the projection through two diffraction planes indexed with A and B is represented as follows: where the amplitudes are denoted by a ij , and i and j denote the projection through upper or lower slits. The A and B in the MPD set-up denote the indices of planes for the projection of a single photon at different time instants rather than the indices of two photons as in the entangled state of two spatial qubits of two photons. There are four different projection history states where the wave function whose intensity to be measured on the final detection plane, i.e., � 3 (x) , is described as the interference of four different wave function histories corresponding to each trajectory, i.e., � 3,j (x) for j ∈ [1,4]: corresponds to a uu |U A U B � and the other components are defined as shown in Fig. 2. Each component depends in a complex manner on the slit geometries as modeled by the RTF. QPC applications of MPD based high dimensional entangled states do not include any measurement by closing or opening slits but a final measurement on the detector plane obtaining the complicated interference pattern of exponentially many Feynman paths 6 . QPC based on FO promises expanding the set of solvable 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 the Open Issues and Discussion section. Propagation through Fourier optical systems based on Fresnel diffraction is modeled emphasizing the quantum nature of Fresnel diffraction and FO in the Methods section. Next, MPD modeling is proposed for Fresnel diffraction and arbitrary LCT based optical systems by utilizing the proposed kernels.
(3) n π n P N−1,s n,N−1 ⊙ P N−2,s n,N−2 ⊙ . . . ⊙ P 1,s n,1 ⊙ [ρ 0 ] and electron based theoretical formulation is extended to optical systems of LCT as shown in Fig. 3. The kernel of one dimensional (1D) quadratic-phase system or LCT converting the input signal f (x 0 ) to the output 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 kernel matrices K FS (x 1 , x 0 ) and K HO (x 1 , x 0 ) denoting FSP kernel in phase-space optics 34 and the kernel based on quantum mechanical harmonic oscillator (HO) modeling of the evolution of light wave function 15 (in analogy with FRFT implementation), respectively, are defined in the Methods section while discussing quantum FO. Some simple examples of LCT matrices for propagation of length L j,j+1 between jth and (j + 1) th planes are given as follows 34 : where M S scales with � j+1 (x j+1 ) ≡ (1 / √ a j,j+1 ) � j (x j+1 / a j,j+1 ) , the kernel for M L is e −ı π x 2 1 / ( f ) and m ≡ 2 π / ( c) is defined in the Methods section after discussing (47).
The varying forms of wave functions on the measurement plane extending (2) are modeled which are promising to be utilized in QC applications. It is presented next such that obtained forms are similar to (2) while having higher flexibility of system design. The theoretical modeling of BB functions for quantum HO based or FRFT based light propagation modeling with Gaussian sources is presented next with the wave function in (16). www.nature.com/scientificreports/ Gaussian source case is also extended to arbitrary LCTs. Similarly, the wave functions for arbitrary LCTs with HG sources are presented in (25) and (27) next. 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 35 . 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 34 . FRFT with scaling is a special case of LCT as discussed in the Methods section. Therefore, FSP, FRFT and arbitrary LCT based QPC set-ups are implemented with the universal configuration in Fig. 3.
QPC with Fresnel diffraction and FRFT by using Gaussian sources. Firstly, two special cases of LCTs are considered, i.e., FSP of light and propagation modeled with FRFTs denoting graded-index media propagation as the solution of the quantum HO in Ref. 15 . Furthermore, we assume that the source distribution has a Gaussian form π while HG waveforms as eigenfunctions of FRFTs 34 are considered for the general case of LCTs in the next section. It is assumed that the optical system between the planes results in the kernels K FS (x 1 , x 0 ) and K HO (x 1 , x 0 ) defined in (45) and (47) based on Fresnel diffraction integral for free space and quantum HO solution 15 , respectively. The definition and the derivation of HO based kernel with the following kernel matrix for the propagation duration of t 01 are detailed in the Methods section while we are discussing the quantum mechanical modeling of FO, i.e., denoting with quantum FO: The important observation is that 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 algorithms for calculating the iteration parameters as shown in Table 1 in the Methods section. The kernel K FS (x 1 , x 0 ) has the same form with K m,FS (x 1 , x 0 ) used for QPC modeling in Ref. 6 by replacing the electron mass m with the photon equivalent mass m . Therefore, the same formulations are utilized for the cases of FS and HO solutions while modeling the sampled wave function on the sensor plane with iterations and the resulting structure of problem solving capabilities.
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: ξ j,n , and the complex vector � h N−1,n and the matrix H HO/G N−1,n are defined in the Methods section for the HO case with simplified formulation compared with the case for electron based FSP set-up in Ref. 6 . The corresponding iteration parameters are given in Table 1 in the Methods section.
We have not included the effects of special forms of K HO (x j+1 , x j ) with ω t j,j+1 = k π for k ∈ Z corresponding to integer multiples of FRFT order 2 since the result is 34 . 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. ( b j,j+1 � = 0 ) with Gaussian sources has the same form with K HO while with different algorithms for calculating iteration parameters in Table 1 in the Methods section and replacing H HO/G N−1,n with H LCT/G N−1,n . Therefore, all the derivations utilized for K HO including the explicit forms of the wave function are applicable. We have not included K (a j,j+1 , b j,j+1 , c j,j+1 , d j,j+1 ) LCT with b j,j+1 = 0 for simplicity. Two simple cases are scaling and chirp multiplication with a j,j+1 = d j, 34 . These cases further improve the flexibility of the LCT system to realize the desired transformation on the wave function.
The highly complicated expressions for the polynomials are explicitly shown in Table 2 in the Methods section such that they are obtained by using the iteration method in Table 1 in the same section. It is possible by using the explicit expressions directly to perform various gedanken experiments and computational complexity analysis with any number of slits and LCTs.
A simple numerical example is presented as shown in Fig. 4. The scaling property of thin lens is utilized in the Numerical Results section to improve the intensity of the diffraction on the final detection plane. For example, a simple Gaussian source beamwidth of σ 0 = 20 ( µ m) and = 650 (nm) shown in Fig. 4b is scaled by shifting the position of the lens of focal length 60 (mm) inside the interval of fixed length of L 01 = L a,01 + L b,01 = t 01 × c with t 01 = 1 (ns) as shown in Fig. 4a. The shift is modeled with the ratio r L = L a,01 / L 01 . It is observed in Fig. 4c that the intensity of the wave function can be focused with respect to the positions of the slits on the first plane.
Hermite-Gaussian Sources If the source is chosen as the standard HG waveform of 34,65 , then � N,n (x N ) for nth path is obtained as follows by using the integral equality of HG functions in the Methods section: 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 3 in the Methods section. Simple algebraic manipulations of (22) to extract � x N−1,n dependent parts result in the following simplification:    (23) that each different path results in a different shift on Hermite polynomial determined with � η T N−1,n � x N−1,n even for the uniform β j for each slit on jth plane. As a result, the final wave function on the sensor plane denoted with � HG N (x N ) for the general case of non-uniform slit widths defined with β j,n for j ∈ [1, N − 1] and n ∈ [0, N p − 1] is given by the following: where ϒ HG N,n ≡ χ 01 N−1 j=1 χ a,j,j+1,n e − ı π (N−2) / 4 and with the similarity to the form in (16) for the HO solution except multiplicative Hermite polynomial for each nth path. The complexity of calculating the Gaussian form in (16) is classically hard as thoroughly discussed in Ref. 6 which requires to compute a special form of partial sum of RTF while the complex vector � h N−1,n and the matrix H HO/G N−1,n varying for each path making it much harder compared with the computation of conventional partial sum of RTF. Therefore, the complexity characterization of computing � HG N (x N ) is an open issue while it is expected to be significantly hard since each summation term depends on path index n with varying vector and matrix parameters while also including a product term of Hermite polynomial for each path making it harder.
If the uniform slit width case is chosen and the path independent variables are denoted with ϒ HG N,n = ϒ HG N , (25) 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 Ref. 66 where is similarly expected to be significantly hard since the mathematical form is more complicated compared with the partial sum of RTF.
The set-up parameters including the slits, lenses and inter-plane distances are required to be tuned in order to obtain the desired vectors  (16), (25) and (27) for the targeted number theoretical problems. Next, Wigner distribution is defined where its negative volume is regarded as an indicator of non-classicality.
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 6,42 : The negative volume of Wigner function defined in Ref. 42 and utilized in Ref. 6 to describe the increasing non-classicality or time-domain entanglement resources in Ref. 12 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: In this article, we propose that a similar experiment could be formulated in the QPC set-up as the problem of finding the distribution of light intensity on the photodetector array plane by randomly generated slit positions and widths analogical to random circuit sampling 45,67 . Figure 5a shows the system architecture that is similar to the randomly generated quantum gates. The aim is to perform a complexity analysis of a randomly generated QPC architecture. The total number of Feynman paths is expressed as follows for L diffraction planes: Suppose that energy decreases by 1/s ( s > 1 ) to a total of 1/s L and the number of significant paths decreases by 1/r ( r > 1 ) leading to a decrease in the total number of effective paths by 1/r L approximating the final intensity distribution. Therefore, if we define the number of paths that can be realized for unit source energy with N path , then the following is obtained: www.nature.com/scientificreports/ If we measure the Hilbert space size created by the total number of paths by defining the virtual qubit number and assuming s ≡ 2 s * , r ≡ 2 r * , m ≡ 2 m * 2 2 (assuming the minimum of m = 8 with m * = 1 ), k ≡ 2 k * , then the following definition is obtained: The basic expression determining the size of the Hilbert space is denoted with the gain G = m * + k * − s * − r * . It provides the cumulative effect of increasing number of paths due to the linearly increasing number of the slits with the coefficient m and the initial number of the slits k combined with the decreasing number of paths due to the inter-plane attenuation coefficients s and r for the effective number of the paths. As shown in Fig. 5b for L = N − 1 , even at very low gain rates, e.g., G = −5 with m * = 1 , that is, where the spreading energy drops very quickly and the number of significant paths is too low, the number of virtual qubits reaches hundreds. Furthermore, m * + k * , which can be designed flexibly in a multi-slit architecture, is adapted against the low gain.
In addition, even with N = 10 planes, Hilbert space size of approximately reaching hundreds of virtual qubits is obtained for the case of high gain G. For example, assume the worst case situation such that diffracted photon forms a large amplitude path by diffracting through a locally limited number of slits on the next plane denoted by the parameter r . In other words, the slit locations distant apart on consecutive planes will not form a large amplitude path and the number of effective paths increases as the multiplication by r j after diffracting through j consecutive planes (by assuming there is enough number of closely spaced slits on the next plane). This can be adjusted by increasing the inter-plane distance such that each diffracted beam will expand to a larger area on the consecutive plane. Moreover, assume that the slits are placed close enough to keep the probability of the diffraction through the next plane roughly constant, i.e., 1 / s , as observed in MPD simulation studies in Ref. 6 . Then, if the condition s < r is satisfied, the number of effective paths will increase with the multiplication by ( r / s) j through j consecutive planes. For example, effective number of paths in Ref. 6 is observed to be increasing with ( r / s) 3 > 70 even by removing many effective paths ( Fig. 7(b) in Ref. 6 ). In other words, assuming r / s ≈ 2 2 allows to reach 2 100 effective Feynman paths with N ≈ 51 planes. Therefore, without requiring extensive simulations, it is clear that the number of paths increases exponentially with specially adjusted set-up parameters of inter-plane distance, slit widths and distributions. Each path will form a unique contribution to the overall intensity. There is no apparent way of calculating the exact final intensity other than identifying and summing the contribution of each path. Then, increasing the complexity of MPD set-up makes it harder to calculate the contribution of each path until reaching to the QS scale. It is possible to compare roughly with Google QS experiment where the computational complexity requires the calculation of 4 31 × 2 4 = 2 66 different Feynman paths for 53 qubits and 20 cycles (Table XI in Ref. 1 ) as shown in Fig. 5b. Compared to the q path = 66 , the proposed QPC system architecture suggests to perform future QS experiments with a simple system structure.
The open issues include the rigorous characterization of the computational complexity of sampling from MPD exploring the relations among the paths in terms of magnitude and distribution. Furthermore, the interplane gain G is required to be both theoretically modeled and experimentally measured for random and large diffraction architectures. Another open issue is to analyze the modeling of the sampling problem of QPC with universal quantum circuits and to determine the computational complexity class, e.g., the relations with BQP 10 and complexity theoretical fundamentals of QS experiments 47 . Experimental implementation requires slit design and manufacturing, sensitive photon detection due to the attenuation after large number of planes and spatially coherent light sources covering all the paths reaching to the detector plane. The number of slits, i.e., determined by the parameters k, m and N, is limited by the capability to realize significant number of small width slits, e.g., in micrometer scale, on an appropriate planar surface such as by patterning metallic slits on glass substrate 68 . The beam width and inter-planar distances should be adapted for spatial coherence of the light diffracting through planes 12 . However, the linear modeling of the architecture in Fig. 5a and the expansion of light beam through propagation allow to realize a feasible architecture in future experimental implementations.
On the other hand, achieving QS experiments allows to adapt certified random number generation protocols for the QPC architecture [43][44][45][46] . Although there are recent high speed, e.g., on the orders of several Gbit/s, random number generation protocols working in a local manner and exploiting the sampling of interference based intensity fluctuations of laser pulses such as Refs. 69 and 70 , the idea of randomness extraction from QS experiments proposed by Aaronson 43 in a way allowing to download from a remote and trusted public source is new. The user interacts with a remote QC and makes it to generate random bits without any trust to the QC itself. Similar to Aaronson's protocol, it is possible to firstly collect random numbers from a trusted computer. Then, using these numbers, the widths and planar distributions of the slits are determined to have a random diffraction set-up by assuming that the mechanical structure of the device can be modified remotely. Then, intensity distribution in the photodetector array is measured. Therefore, both the random structure of MPD set-up and interference of the exponential number of paths result in a very difficult measurement output intensity to efficiently calculate with classical computers. However, it is an open issue whether it is possible to utilize the proposed remote QC device based on QPC similar to the Aaronson's protocol which realizes sampling from the n-qubit output of the quantum circuit and performs Heavy Output Generation (HOG) tests. On the other hand, QPC does not allow to sample the probability of a single path but the interference of exponentially many number of paths. Therefore, it is a challenge and open issue to adapt the interference sampling operation in QPC for a similar complexity theoretical proof of randomness generation.
The example output through the slits Y j for j ∈ [1,2] in Fig. 6b depends on interfering quantum superposition of input combinations as follows while assuming path independent forms of the variables for simplicity, i.e., where G(y, β N , x N ) is the slit mask function depending on the output slit y = Y i for i ∈ [1,2] of the Nth plane and β N is the fixed slit width with path-independent assumption for simplicity. The parameters are possible to depend on each path with variable slit masks. Exponentially large number of synaptic chains (paths) through slit inputs, their quantum interference and simplicity to sample the intensity output |O| 2 provide significant opportunities to exploit for quantum advantages. The challenges include designing the quantum neuron based on the slit positions as inputs while changing the weight in a controllable manner. Besides that, extensive simulation studies are required to practically observe the quantum advantages for various problems. The positions of the slits are required to be modified dynamically with special designs. In addition, large scale QNN implementations both in simulations and experiments are required to observe the performances in various problems, e.g., pattern recognition or machine learning for very large problem sizes. (1) are expressed with RTF as follows 60,61 : 71 is calculated using k and ω together with nonlinear spectrum data (Appendix to Section 24 in Ref. 51 ) and the partial sum of RTF denoted as M converging to for M → ∞ is defined as follows:

Solution of nonlinear Schrödinger equation. Finite-band solutions of NLSE in
where Ŵ is a complex matrix, y is a complex vector and � a T ≡ [a 1 a 2 . . . a N−1 ] . As shown in Ref. 6 , if j ∈ [1, N − 1] and a j ∈ S M ≡ [−M, M], we select X j,i ∈ S M x j and also if the slit widths are kept constant for each plane, then . . .

numerical results
MPD set-up with two diffraction planes and single sensor plane is numerically analyzed for both Gaussian and HG sources with beam width and waist sizes of σ 0 = 20 µ m and W 0 = 200 µ m, respectively. The set-up is shown in Fig. 7 with N = 3 planes. The source waveforms are shown in Fig. 8a,b, respectively. HG order is set to l = 10 with highly oscillatory and negative initial V 0 of 1.076. The wavelength of the light is chosen in the red spectrum of = 650 nm while the low cost laser sources are commercially available in a wide spread manner.
Two different set-ups composed of LCT and FSP systems as shown in Fig. 7a,b, respectively, are compared where the LCT system includes thin lenses between the diffraction planes while not included in the FSP system.  ] mm focus the light intensity to more compact areas on the consecutive planes compared with FSP. It is assumed that the propagation between the second and third planes includes only FSP without any thin lens to simplify the set-up. K 1 = 11 and K 2 = 27 slits are used on the first and second planes, respectively. The slit positions and widths on the first plane, as shown in Fig. 8c,d, respectively, are adapted to the maximum intensity locations of HG source propagation on the first plane in LCT system while the ones on the second plane are chosen uniformly with the separation of 40 µ m and the width of β 2,j ≡ β 2 = 8 µ m. K a 01 , b 01 , c 01 , d 01 LCT and K a 12 , b 12 , c 12 , d 12 LCT are calculated by using (13). FSP has less control over the propagation of light compared with LCT based FO. FSP spreads the light without any tuning to the slit positions by reducing the probability of the photon to reach to the consecutive planes after diffraction. Therefore, in numerical analysis, LCT is shown to improve the probability of photon detection on the sensor plane ( P E ) and also the negative volume of Wigner function compared with FSP. The vectors of P E (j) and V j composed of the values on the first, second and third (sensor) planes for Gaussian sources are denoted with P G,FSP E and V G,FSP , respectively, for FSP while with P G,LCT E and V G,LCT for LCT. The cases with HG sources are denoted with the superscript of HG. It is an open issue 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 photon reaching to the sensor plane.
Hermite-Gaussian sources. The waveforms on the three planes in spatial domain are shown in  Fig. 9d,h, for j = 2 and j = 3 , respectively, while Wigner distributions on the second and third planes scaled with are shown in Fig. 9i,j for FSP and, (k) and (l) for LCT. It is observed that LCT provides significantly larger path magnitudes while  Fig. 10d,h, for j = 2 and j = 3 , respectively. Similar to the HG sources, LCT improves diffraction probabilities significantly compared with FSP. Wigner distributions scaled with are shown in Fig. 10i-l having different characteristics compared with HG case in Fig. 9i-l. It is similarly observed that LCT provides significantly larger path magnitudes. � V G,LCT = [0 0.842 1.426] and � V G,FSP = [0 1.21 0.93] are obtained with increasing interference complexity through diffraction on consecutive planes in LCT case while starting with purely classical Gaussian source of zero negative Wigner volume. V 2 of 0.842 for LCT is smaller than 1.21 for FSP on the second plane. This is due to the both the specific set-up parameters and more diverse distribution of the path magnitudes in LCT after diffraction from the first plane as shown in Fig. 10d. It becomes more difficult on the third plane to correlate the distribution of the path magnitudes shown in Fig. 10h with V N = V 3 shown in Fig. 10j,l. In other words, complexity behaves differently compared with the transmission probability while requiring simultaneous maximization depending on specific set-up as an open issue as discussed in the Results section while presenting Wigner distribution analysis.

open issues and discussion
There are some open issues to best exploit photonic QPC method based on MPD and FO.  Tables 1 and 3. Therefore, adapting the physical set-up parameters to the desired form of partial sum of RTF for the target number theoretical problem, and characterizing the path distributions and the negative volume of Wigner function explicitly are important open issues.
In the proposed formulation, the qubits are obtained through the tensor product structure of projections at different time instants on the contrary with the spatial encoding and entanglement of multiple photons. However, www.nature.com/scientificreports/ the histories of photon trajectories are not formulated for realizing conventional quantum gate implementations. Implementations of the quantum circuit gates are required to obtain universal QC architectures. Therefore, proposed QC formulation is limited to utilization of the interference of exponentially increasing number of Feynman paths based on the superposition and coherence properties of light source. Implementations of quantum circuits and fundamental search algorithms such as Grover search are future works to clarify the potential future scope of QPC based computing architectures in terms of universal QC capabilities. On the other hand, FO based QPC implementation has two main advantages resembling Boson sampling advantages in a different context 10 : (a) not utilizing multiple photons as qubits getting rid of the coupling disadvantages while exploiting single photon trajectories and (b) utilizing the free entanglement, for the first time, of the classical light obtained through freely available temporal correlations among the projections at different time instants. Boson sampling compared with QPC utilizes still multiple indistinguishable photons (but not as qubits) while requires single photon generation and detection to exploit the free entanglement among indistinguishable photons through multi-mode interferometer with the regarding boson statistics. Realizing perfectly Gaussian slits compared with the conventional rectangular apertures is an important open issue for matching the experimental results with the proposed theoretical model. However, any slit structure can be represented as a composition of Gaussian slits by using the method defined in Ref. 72 and applied successfully in optical diffraction theory and experiments 73,74 . The one dimensional slit mask 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 (40), then the summations in (16), (25) and (27) should be made for each β i of the single slit. The functional form with partial sum of RTF should be calculated and summed for each combination of β i through all the slits. Therefore, non-Gaussian slits can possibly realize the solutions of much harder computational complexity problems as an open issue.
There are some factors effecting the degree of 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 Refs. 6,41,75 . All these considerations potentially lead to unavoidable errors requiring quantum error correction studies adapted to QPC architectures 76 . Another open issue is related to the utilization of the measurements on all the sensor planes for computational purposes not only the final sensor plane since they include diffraction through previous planes. Theoretical models are required to exploit the sensor measurement results.

Methods
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 77 : 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 effects on the wave amplitude, G − ≡ exp(ı k r 01 ) / r 01 − exp(ı k r 01 ) / r 01 is the Green's function vanishing on the diffraction surface for the first type of solution of Rayleigh-Sommerfeld formula, r 01 ≡ |� r 0 − � r 1 | and k ≡ 2 π / for the monochromatic light source of wavelength as shown in Fig. 11 77 . 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:  for Gaussian sources.
Scientific RepoRtS | (2020) 10:10968 | https://doi.org/10.1038/s41598-020-67364-0 www.nature.com/scientificreports/ 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 34 . As a class of linear integral transforms, they include as special cases the Fresnel transform and FRFT, simple scaling, chirp multiplication and some other operations. Spatial distribution of light in phasespace optics for the class denoted by quadratic-phase systems is mathematically equivalent to LCTs (Chapters 3 and 8 in Ref. 34 ). 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 Ref. 15 , 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 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 the refractive index distribution of n 2 (x) = n 2 0 (1 − (x / χ) 2 ) 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 (Section 8.3.3 in Ref. 34 ): where α = d gri / χ . There is a FRFT relation between scaled versions of the input f (x) and output ĝ(x) with FRFT order α as ĝ(x) = e −ı d gri / (2 χ) −1 /4 χf (x χ ) 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.  for Hermite-Gaussian sources ( b j,j+1 � = 0).