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 [Formula: see text] Feynman paths, quantum neuron implementations and solutions of nonlinear Schrödinger equation.

The scalability of quantum resources including qubits and quantum gates, improved error correction capabilities and practical problem solving ability are the most important challenges for modern quantum computing (QC). Recent quantum supremacy experiments (QS) of Google as a success milestone for the human history of computing emphasize the importance of these properties in their set-up with 53 qubits and 20 cycles reaching Hilbert space size of Feynman paths making it significantly difficult to classically calculate their result[1,2]. QS experiments target to show the computational capability of QCs such that feasible computations obtained with QCs require significant resources to perform with classical computers[3]. In the global and highly competitive race including technology giants, a wide variety of but quite complex hardware architectures are used. For example, Google, IBM and Rigetti Computing use superconducting circuits while Microsoft using topological anions generated by frozen nanowires, Ion-Q using ion traps at room temperature and D-Wave using quantum annealing technology[4]. On the other hand, both the challenges of scalability and practical problem solving capability continue to exist such that QCs that can show significant advantages in practical problem solving compared to conventional computers require much more resources such as thousands of logical qubits and hundreds of thousands physical qubits[5]. Therefore, building QC system architectures which are more tolerant to noise and decoherence combined with capabilities of error correction, practical problem solving, low hardware complexity and resource scalability is significantly important for near-term advantages of QC.

Linear optic (LO) based QCs (LOQCs) have a special place among the existing QC architectures due to easy manipulation of photons, unique features of photons not interacting with the environment in terms of decoherence, working at room temperature and maturity in classical optics for centuries[7-9]. For example, boson sampling[10] as a candidate for LOQC based QS promises experimental implementation of QS in the future while the recent experimental achievements in Ref.[11] improve the complexity of the solved problems gradually reaching Hilbert space size of . However, the existing LOQC candidates for QS have still the fundamental challenges of scalability of the resources, e.g., the requirement of single photon sources and receivers, and practical problem solving capability, e.g., matrix permanents in boson sampling. Furthermore, multi-photon entanglement resources and quantum circuits are challenging to create due to the difficulty in the interaction of photons with each other. Multi-plane diffraction (MPD) based QC system denoted by quantum path computing (QPC) as shown in Fig. 1 is recently proposed in Ref.[6] as one of the simple LOQC architectures targeting to realize scalability and practical problem solving capabilities. Sampling from the interference of exponentially increasing number of propagation paths and the freely entangled projections at different time instants[12] are utilized for QC purposes, for the first time, by exploiting coherent and classical sources such as standard laser sources. Exponentially increasing number of Feynman paths with respect to the given amount of slits and diffraction planes makes the classical calculation of the interference output significantly hard[6]. The unique form of temporal correlation freely available among the exponentially increasing number of Feynman paths in the MPD set-up is denoted as quantum path entanglement (QPE) in Ref.[12] as a novel resource to exploit for QC based on the coherence and superposition of the classical light source[13,14]. In fact, a quantum mechanical propagator for photons with the form of the classical Fresnel diffraction integral is verified in Ref.[15] such that the classical intensity of the field is proportional to the probability density of photon detection for the position observable transversal to the propagation and in the limit for large number of quanta. Therefore, FO based set-up exploits classical light source and its intensity measurement while exploiting quantum temporal correlations in a unique MPD design. QPC promises a significant alternative to cope with the fundamental challenges of scalability of multi-photon entanglement resources and the complex requirements of single photon sources and detection mechanisms observed in conventional linear optical QC systems[9,16-18]. It has important all-in-one advantages combining the utilization of the classical sources, i.e., either fermion or boson, hardware simplicity based on diffraction slits and detection with conventional photon counting intensity detectors without requiring simultaneous detection of multiple photons in multiple registers[10].

Figure 1

QPC set-up composed of diffraction planes, FSP between the planes and a single sensor plane on which exponentially increasing number of propagation paths interfere[6].

Besides that, another QC architecture related to multi-slit structures is denoted as duality computer (DC) which exploits duality parallelism for performing different gate operations on the sub-wave functions through sub-waves corresponding to each slit[19-22] while utilized in various machine learning[23] and photonic chip applications[24]. The quantum wave divider (QWD) divides wave functions into sub-waves. In addition, sub-sub-waves can also be obtained in a multi-level QWD while quantum wave combiner (QWC) is utilized to combine the waves after performing operations on each path. It is firstly required as an open issue to explicitly model the theoretical computational complexity and QC capabilities of a multi-level version of DC with respect to the targeted QWD/QWC configurations before comparing with MPD design. For example, combined operations of uniform QWD/QWC in a DC leaving the state unchanged is not comparable with MPD modifying the state with diffractions as time evolves. Single level complexity discussions in the DC related literature include the requirement of an extra qudit[20,21] for simulating DC device with an ordinary quantum computer which is not comparable with the tensor product structure of MPD requiring multiple qudits. The division of each wave into a sub-wave and then into a sub-sub-wave could require an exponentially increasing number of slit resources making it challenging as discussed next for Ref.[25]. On the other hand, DC based architectures have a long history and maturity along two decades with many capabilities and applications including practical database search without requiring extra resources[19] verified with quantum circuit modeling in Ref.[20]. DC utilizes linear combination of unitaries (LCU) for QC while becoming one of the five major techniques for designing quantum algorithms[26], in multi-party secure computation[27] and HHL quantum algorithm[28].

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 input signal array with an output signal creating 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 inter-planar 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-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 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-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 all-optical 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 ( and period ) 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 diffraction planes with slits on each plane for and a single sensor plane indexed with N while the central position of a slit is given by for as shown in Fig. 1a. Each slit is assumed to apply a spatial filtering of , i.e., slit mask function, where determines the slit width. The wave function on jth plane is denoted with which is the wave form after diffraction through the previous planes, i.e., with the indices , 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., , while nth path includes the diffraction through a single slit on each plane with the corresponding wave function on jth plane. Assume that nth path passes through the slit indexed with on jth plane and we define the path vectors and where and is the transpose operation. The mapping between the path index n and slit index for the path is defined with the function where is predefined for each n. Furthermore, is the column vector of length k with all zeros and is the square matrix of all zeros with the size . Similarly, rectangular matrices are shown with . In the rest of the article, a parameter B depending only on but not on is denoted with on each jth plane including over the symbol. Therefore, if the slits are chosen with the same specific to each plane, then becomes independent of n and is converted to the notation .

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 with the following special form to utilize in solutions of important and classically hard number theoretical problems:where , is a sampling interval, , and . The complex valued matrix and the vector have the values depending on for and corresponding to the specific selection of slits in the path , inter-plane durations for the particle propagation, particle mass m (for electron based set-ups in Refs.[6,12]), beam width 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-51]. The second method utilizes MPD with the phase term in for period finding[52] and the solution of specific instances of SDA problems[53].

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 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., . 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 represents the projection operator through the slit indexed with , as 0 or 1 determines a compound set of trajectories, denotes tensor product operation and 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 and 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 , 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., , is described as the interference of four different wave function histories corresponding to each trajectory, i.e., for : corresponds to 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].

Figure 2

MPD based two-qubit state represented with four quantum histories of a single photon with the tensor product structure in time domain in analogy with the entangled state of two spatial qubits of two photons.

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.

Quantum path computing with Fourier optical systems

The set-up in Fig. 1 defined with FSP 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 to the output is represented as follows:In matrix notation, it is shown with the following unit-determinant matrix:where and the matrix for the composition of two consecutive systems represented by and is calculated by the multiplication . The kernel matrices and 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 between jth and th planes are given as follows[34]:where scales with , the kernel for is and is defined in the Methods section after discussing (47).

Figure 3

Photonic QPC architecture composed of classical light source, MPD set-up composed of diffraction planes with slits on jth plane, general LCT phase-space optics represented with the matrix elements between the planes indexed with j and , and a single sensor plane on which exponentially increasing number of propagation paths interfere. Each is implemented with sections of FSP for the lengths of and , and a thin lens of focal length between them.

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). 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 is implemented in phase-space optics by consecutive applications of FSP of length , then thin lens of focal length , and another FSP of length [35]. LCT matrix is calculated as follows:where and , 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 of 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 and 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 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 results in the final intensity distribution of MPD with the same form of while with different algorithms for calculating the iteration parameters as shown in Table 1 in the Methods section. The kernel has the same form with used for QPC modeling in Ref.[6] by replacing the electron mass m with the photon equivalent mass . 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.

Table 1

Iteration parameters for FPI modeling of MPD with the kernels and for Gaussian sources.

	Formula for K_HO based MPD (ωt_j,j+1/π∉Z)	Formula for KLCT(a,b,c,d) based MPD (b_j,j+1≠0)
Ψ_0(x_0)=exp(-x02/(2σ02)) / σ_0π & Ψ_1(x_1)=χ_0exp(A_0x12+ıB_0x12)
A_0	-m^02σ02/(2cos^2(α_0)m^02σ04+2λ02)	-2π^2σ02/(4π^2a012σ04+b012)
B_0	cos(α_0)m^_0λ02-sin^2(α_0)m^02σ04/(2λ_0(cos^2(α_0)m^02σ04+λ02))	πd_01/b_01-4π^3a_01σ04/(b_01(4π^2a012σ04+b012))
χ_0	π^-1/4m^_0σ_0/(cos(α_0)m^_0σ02+ıλ_0)	exp-ıπ/42πσ_0/(b_01-2ıπa_01σ02)
Ψ_2,n(x_2)=χ_0χ_1,nexpA~_1,nx22+ıB~_1,nx22+C_1,nx_2+ıD_1,nx_2
A~_1,n	β~1,n2m^1222A_0β~1,n2-1/(2ζ~_1,n)	2π^2β~1,n2(2A_0β~1,n2-1)/ζ~_1,n
B~_1,n	(2B_0β~1,n4m^122cos(2α_1)+cos(α_1)λ_1m^_12ϱ_1)/(2ζ~_1,n)-β~1,n4m^123cos(α_1)sin^2(α_1)/(2λ_1ζ~_1,n)	πd_12ζ~_1,n-4πβ~1,n4(πa_12+B_0b_12)/(b_12ζ~_1,n)
χ_1,n, C_1,n, D_1,n	ξ~_1,nexp(p~_1,1,nx1,n2), ζ~_1,c,nx_1,n, ζ~_1,d,nx_1,n
Ψ_j+1,n(x_j+1)=χ_0(∏k=1jχ_k,n)e^(A~_j,n+ıB~_j,n)xj+12+(C_j,n+ıD_j,n)x_j+1 for j∈[2,N-1]
p~_2,j,n,p~_3,j,n,p~_4,j,n,p~_5,j,n	β~j,n2p~_3,j,n/2,-λ_j/(ıς~_j,n),β~j,n2ζ~_j,c,n,-β~j,n2ζ~_j,d,n	β~j,n2p~_3,j,n/2,ς~_j,n/ζ~_j,n,β~j,n2ζ~_j,c,n,-β~j,n2ζ~_j,d,n
A~_j,n	β~j,n2m^j22A~_j-1,nβ~j,n2-1/(2ζ~_j,n)	2π^2β~j,n2(2A~_j-1,nβ~j,n2-1)/ζ~_j,n
B~_j,n	m^_j(2B~_j-1,nβ~j,n4cos(2α_j)m^_j+cos(α_j)λ_jϱ~_j,n)/(2ζ~_j,n)-(β~j,n4m^j3cos(α_j)sin^2(α_j))/(2λ_jζ~_j,n)	πd_j,j+1ζ~_j,n-4πβ~j,n4(πa_j,j+1+B~_j-1,nb_j,j+1)(b_j,j+1ζ~_j,n)
C_j,n, D_j,n	ζ~_j,c,nx_j,n+p~_4,j,nC_j-1,n+p~_5,j,nD_j-1,n, ζ~_j,d,nx_j,n-p~_5,j,nC_j-1,n+p~_4,j,nD_j-1,n
χ_j,n	ξ~_j,nexp(p~_1,j,nxj,n2)×exp(p~_2,j,n(C_j-1,n+ıD_j-1,n)^2)×exp(p~_3,j,n(C_j-1,n+ıD_j-1,n)x_j,n)
The following variables defined for j∈[1,N-1]
λ_j or λ~_j,n	ħt_j,j+1	b_j,j+1(A~_j-1,n+ıB~_j-1,n)+ıπa_j,j+1
p~_1,j,n	-(2λ_j(A~_j-1,n+ıB~_j-1,n)+ıcos(α_j)m^_j)/(2ıς~_j,n)	λ~_j,n(b_j,j+1-2β~j,n2λ~j,n∗)/ζ~_j,n
ς~_j,n, ξ~_j,n	β~j,n2(cos(α_j)m^_j+2λ_j(B~_j-1,n-ıA~_j-1,n))+ıλ_j,β~j,n2m^_j/ς~_j,n	b_j,j+1(b_j,j+1-2β~j,n2λ~j,n∗),2πβ~j,n2/(ı(b_j,j+1-2β~j,n2λ~_j,n))
ϱ~_j,n	4β~j,n4(A~j-1,n2+B~j-1,n2)-4A~_j-1,nβ~j,n2+1
ζ~_j,n	4B~_j-1,nβ~j,n4cos(α_j)λ_jm^_j+β~j,n4cos^2(α_j)m^j2+λj2ϱ~_j,n	bj,j+12ϱ~_j,n+4πa_j,j+1β~j,n4(πa_j,j+1+2B~_j-1,nb_j,j+1)
ζ~_j,c,n	β~j,n2m^_j(2B~_j-1,nλ_j+cos(α_j)m^_j)/ζ~_j,n	4πβ~j,n2(πa_j,j+1+B~_j-1,nb_j,j+1)/ζ~_j,n
ζ~_j,d,n	λ_jm^_j(2A~_j-1,nβ~j,n2-1)/ζ~_j,n	2πb_j,j+1(2A~_j-1,nβ~j,n2-1)/ζ~_j,n

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 dependent parts and summing the contributions from each path as follows:where , and the complex vector and the matrix 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 with for corresponding to integer multiples of FRFT order 2 since the result is (inserting into )[34]. This case can be simply realized by assuming that spatial filtering operations of the slits on jth and 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.

QPC with arbitrary linear canonical transforming optical systems

Gaussian Sources The resulting final intensity of MPD propagation for the case of () with Gaussian sources has the same form with while with different algorithms for calculating iteration parameters in Table 1 in the Methods section and replacing with . Therefore, all the derivations utilized for including the explicit forms of the wave function are applicable. We have not included with for simplicity. Two simple cases are scaling and chirp multiplication with resulting in [34]. These cases further improve the flexibility of the LCT system to realize the desired transformation on the wave function.

A simple example is presented with explicit expressions as follows for , , , for , fixed slit width parameters and for simplicity and :where the following parameters are derived by using Table 1 in the Methods section: where it is observed that the parameters are some rational complex polynomials of , and with the coefficients in terms of polynomial expressions of , and for with varying orders of (, , ): reaching (0, 2, 2) for ; (2, 2, 2) for ; (0, 0, 2) for ; (2, 0, 2) for ; (6, 4, 8) for , , and ; (8, 4, 8) for , (4, 2, 5) for ; (4, 2, 6) for ; (4, 4, 4) for , and .

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.

Table 2

Polynomial expressions in (18-21) for the case of and Gaussian source with ( for ).

pol_1	π(-2ıb_01q_7πβ22-2π(2π(a_01q_7-b_12q_12)β22+b_12b_23ıq_19)σ02+b_01b_12b_23q_11)
pol_2	2π(4β12π^2(a_01q_7-b_12q_12)β22-a_01b_01b122b_23+2b_12ıπq_17)σ02+b_01(4β12ıπ^2q_7β22+b_01b122b_23(-ı)-2b_12q_23π)
pol_3, pol_4	-2b_01b_12b_23π(2a_01πσ02+b_01ı), π(b_01(b_01b_12q_12-2ıβ12q_7π)-2π(2π(a_01q_7-b_12q_12)β12+a_01b_01b_12ıq_12)σ02)
pol_5	-4β22b012b123b_23π^2q_20(b122b014+4β12π^2q_13b012+4π^2σ02(2β12q_14b012+(4π^2q_28β12+a012b012b122)σ02))
pol_6	(4b012π^2q112β14+b014b122+4π^2σ02(2β12b012b122+(4π^2q192β14+a012b012b122)σ02))×((16β14π^4q72β24+b012b124b232+4b122π^2q_8)b012+4π^2σ02(2β12b012q_9b122+q_22σ02))
pol_7, pol_10	16β12β24b_01b_12π^4(4a_01π^2q_19σ04+b012q_11)q_29, ıβ12β22b_01b_12σ_0(b_01-2ıa_01πσ02)(b_01(b_01b_12-2ıβ12q_11π)-2πq_18σ02)
pol_8	-8β12β22b_01b122b_23π^3(4a_01π^2q_19σ04+b012q_11)(b122b014+4β12π^2q_13b012+4π^2σ02(2β12q_14b012+(4π^2q_28β12+a012b012b122)σ02))
pol_9, pol_12	-8β24b012b122π^3q_20q_29, 2β22b122π^2(-4β12q_13π^2b012-8β12π^2q_14σ02b012-4(4π^4q_28β12+a012b012b122π^2)σ04-b014b122)
pol_11	(2a_01πσ02+b_01ı)(2πq_15σ02+b_01(2πq_11β12+b_01b_12ı))×(b_01q_16-2π(4β12π^2(a_01q_7-b_12q_12)β22-a_01b_01b122b_23+2b_12ıπq_17)σ02)
pol_13	(16β14π^4q72β24+b012b124b232+4b122π^2q_8)b012+4π^2σ02(2β12b012q_9b122+q_22σ02)
pol_14	4π^3σ02(σ02(a012b012b124b232d_23+16π^4β14β24q_4(a_01q_7-b_12q_12)+4π^2b122q_1)+2β12b012b122(4π^2β22q_6+b122b232d_23))+πb012(16π^4β14β24q_5q_7+b012b124b232d_23+4π^2b122q_2)
The functionsq_jforj∈[1,30]utilized while defining the polynomials are defined as follows:q_1≡b122b232d_23β14-2a_01b_12b232q_11d_23β14+a012q_2,q_2≡a232b012b122d_23β24-a_23b012b_12q_27β24+b_23q_3,q_3≡b_23d_23q112β14+2β22b012b_23d_23β12+β24b012d_12q_26,q_4≡b_12(b_12-q_12d_23)-a_01b_12q_11+a_01d_23q_7,q_5≡-b_01b_23d_23+a_12b_01q_25+b_12d_01q_25,q_6≡(β12b232+β22q122)d_23-β22b_12q_12,q_7≡-b_01b_23+a_12b_01q_12+b_12d_01q_12,q_8≡a232b012b122β24+2a_23b012b_12b_23d_12β24+b232q_24,q_9≡4π^2(β12b232+β22q122)β22+b122b232,q_10≡(b_12q_12-a_01q_7)^2,q_11≡a_12b_01+b_12d_01,q_12≡a_23b_12+b_23d_12,q_13≡b122d012β12+2a_12b_01b_12d_01β12+(a122β12+β22)b012,q_14≡2β12π^2β22+b122,q_15≡a_01b_01b_12-2ıβ12q_19π,q_16≡-4ıβ12q_7π^2β22+b_01b122b_23ı+2b_12πq_23,q_17≡a_01q_23-β12b_12b_23,q_18≡2πq_19β12+a_01b_01b_12ı,q_19≡a_01a_12b_01-b_12+a_01b_12d_01,q_20≡b012+4π^2σ02(β12+a012σ02),q_21≡b122b232β14-2a_01b_12b232q_11β14+a012q_8,q_22≡16β14π^4q_10β24+a012b012b124b232+4b122π^2q_21,q_23≡b_23q_11β12+a_23β22b_01b_12+β22b_01b_23d_12,q_24≡a122b012β14+b122d012β14+2a_12b_01b_12d_01β14+2β22b012β12+β24b012d122,q_25≡b_23d_12d_23+b_12(a_23d_23-1),q_26≡b_23d_12d_23-b_12,q_27≡b_12-2b_23d_12d_23,q_28≡q_13a012-2β12b_12q_11a_01+β12b122,q_29≡4b012π^2q_11q_7β14+b014b122q_12+4π^2σ02(2β12b012q_12b122+q_30σ02)and finallyq_30≡4π^2q_19(a_01q_7-b_12q_12)β14+a012b012b122q_12.

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 (m) and (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 with (ns) as shown in Fig. 4a. The shift is modeled with the ratio . 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.

Figure 4

(a) The set-up with a single thin-lens of focal length 60 (mm), (b) Gaussian source with (m) and (nm), and (c) the distribution of the wave function on the first plane by shifting the lens inside the spatial interval of fixed total length with (ns) by varying the ratio of .

Hermite-Gaussian Sources If the source is chosen as the standard HG waveform of for () where is the lth order Hermite polynomial[34,65], then for nth path is obtained as follows by using the integral equality of HG functions in the Methods section:where the parameters , , , , and obtained in an iterative manner for are calculated with simple algebra for nth path and shown in Table 3 in the Methods section. Simple algebraic manipulations of (22) to extract dependent parts result in the following simplification:where for is defined in Table 3, the vectors and , and the matrix are defined in the Methods section. It is observed in (23) that each different path results in a different shift on Hermite polynomial determined with even for the uniform for each slit on jth plane. As a result, the final wave function on the sensor plane denoted with for the general case of non-uniform slit widths defined with for and is given by the following:where 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 and the matrix varying for each path making it much harder compared with the computation of conventional partial sum of RTF. Therefore, the complexity characterization of computing 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.

Table 3

Iteration parameters for FPI modeling of MPD with the kernel for Hermite-Gaussian sources ().

Symbol	Formula	Symbol	Formula	Symbol	Formula
Ψ_0(x_0)=(2^1/4/W_02^ll!)exp(-πx02/W02)H_l(2πx_0/W_0) and Ψ_1(x_1)=χ_01expu_01x12H_l(g_01x_1)
g_01	2πW02a012W04+b012	u_01	-πW02a012W04+b012+ıπd_01b_01-πW02a_01W02b_01a012W04+b012	χ_01	2^1/4W_02^ll!a_01W02-ıb_01a012W04+b012a_01W02-ıb_01a012W04+b012^l

If the uniform slit width case is chosen and the path independent variables are denoted with , , , , , , , then (25) is transformed into the following:It is further simplified as follows by using the useful identity in Ref.[66] where :where . The computational complexity of calculating 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 , , and matrices , and in (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 is defined as Fourier transform of spatial representation of wave function 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 . On the other hand, the probability of the particle to be detected on jth plane, i.e., to be diffracted through th plane, is computed as follows:In this article, the increasing interference complexity and non-classical nature of MPD based time-domain entanglement resources are assumed to be characterized by utilizing and by observing the magnitudes of the interfering paths defined as for each nth path. Therefore, more paths with large magnitudes and emphasize increasing interference and non-classicality. Characterizing the correlation between the distribution of path magnitudes and is an open issue since the behavior is highly set-up specific as observed in the Numerical Results section. On the other hand, path magnitudes throughout the whole plane may not reflect their localized characteristics such as effecting some sample locations more compared with the others. Therefore, it is an open issue to characterize the interference complexity in terms of the intensity distribution of the paths while the path magnitudes are taken as a reference for simplicity in this article. Next, the potential future applications of QPC architecture based on FO are presented.

Future applications

Applications for quantum supremacy and certified random number generation

The scalable structure of QPC with coherent light sources and low complexity FO promise the large scale implementation of QS experiments. Recently, some experiments are performed with 53 qubits with superconductor based architectures of Google as a milestone for the human history of computational capabilities[1,2]. 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 () to a total of and the number of significant paths decreases by 1/r () leading to a decrease in the total number of effective paths by approximating the final intensity distribution. Therefore, if we define the number of paths that can be realized for unit source energy with , then the following is obtained:If we measure the Hilbert space size created by the total number of paths by defining the virtual qubit number and assuming , , (assuming the minimum of with ), , then the following definition is obtained:The basic expression determining the size of the Hilbert space is denoted with the gain . 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 , even at very low gain rates, e.g., with , 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, , which can be designed flexibly in a multi-slit architecture, is adapted against the low gain. In addition, even with 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 . 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 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., , as observed in MPD simulation studies in Ref.[6]. Then, if the condition is satisfied, the number of effective paths will increase with the multiplication by through j consecutive planes. For example, effective number of paths in Ref.[6] is observed to be increasing with even by removing many effective paths (Fig. 7(b) in Ref.[6]). In other words, assuming allows to reach effective Feynman paths with 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.

Figure 5

(a) The special QPC design with planes diffracts the quantum wave function expanding with planes consisting of a linearly increasing number of slits. The positions of the planes (), the number of slits (k) in the first plane and the linear increment ratio m, i.e., the number of slits growing with , allow flexible design according to the width of the laser beam. (b) The virtual qubit Hilbert space size corresponding to the number of Feynman paths versus the number of planes for varying gain variable (G). QS experiments of Google have recently been performed with Feynman paths[1]. Even at significant loss rates of the QPC system, thanks to dozens of planes, the Hilbert space size reaches hundreds of virtual qubit levels and achieves very strong QS capability. It is promising as an alternative system design for next generation QS experiments.

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 different Feynman paths for 53 qubits and 20 cycles (Table XI in Ref.[1]) as shown in Fig. 5b. Compared to the , 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 inter-plane 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-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.

Neuromorphic applications with quantum neuron implementations

In a classical artificial neuron implementation, the output is represented as where f(.) is the nonlinear activation function, is the ith synaptic weight, is the ith input and b is the bias as shown in Fig. 6a. QPC based neuron has inherent nonlinearity with the form of between the input slit position vector and output (O) as shown in Fig. 6b based on the MPD formulation for Gaussian sources in (16). The quadratic weighting , linear weighting and bias b are fed into the nonlinear function for the LCT set-up where , , the bias , and is the measurement position.

Figure 6

(a) Classical artificial neuron implementation with the weight and input vectors of and , respectively, activation function f(.), bias b and the output O. (b) QPC based design of QNNs with quadratic weighting relationship of and nonlinear activation having quantum superposition and interference combining the inputs.

The example output through the slits for 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., , and :where is the slit mask function depending on the output slit for of the Nth plane and 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 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.

Solution of nonlinear Schrödinger equation

Finite-band solutions of NLSE in (1) are expressed with RTF as follows[60,61]:where Riemann spectrum is , Riemann period matrix [71] is calculated using and together with nonlinear spectrum data (Appendix to Section 24 in Ref.[51]) and the partial sum of RTF denoted as converging to for is defined as follows:where is a complex matrix, is a complex vector and . As shown in Ref.[6], if and , we select and also if the slit widths are kept constant for each plane, then , , , , and values become path independent, i.e., , , , , and while the superscript removed for simplicity, and (2) is transformed as follows for :where and the diagonal matrix is formed of the elements . In other words, we can achieve the absolute value of the partial sum of particular RTF by using the measurement result on the sensor plane. Although the calculation of RTF function is classically quite difficult[6,51], it has important applications in areas including geometry, arithmetic and number theory[50], nonlinear spectral theory for ocean and water sciences[51], cryptography and the solution modeling of NLSEs[60,61].

Measurements in the QPC system allow to obtain information about . The most important challenge for utilizing QPC in the solutions of NLSE is to determine the set equation parameters which can be implemented with specific QPC design. Theoretical modeling and extensive simulations are required to examine all the practical sample parameters and systems in the literature where NLSE solutions are achieved with RTF based solution.

Next, numerical simulations are achieved to analyze the effects of FO based components such as lenses and HG sources on the intensity distribution obtained with QPC. Simulation studies for large scale implementations of QPC for future applications are open issues.

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 m and m, respectively. The set-up is shown in Fig. 7 with planes. The source waveforms are shown in Fig. 8a,b, respectively. HG order is set to with highly oscillatory and negative initial of 1.076. The wavelength of the light is chosen in the red spectrum of nm while the low cost laser sources are commercially available in a wide spread manner.

Figure 7

Photonic QPC set-up with Gaussian and HG classical monochromatic light source with nm, two planes for diffraction with the number of slits and , respectively, and specific set-up of (a) LCT and (b) FSP design with the only difference of the existence of thin-lenses of focal length of 63 mm in the LCT system.

Figure 8

(a) Gaussian and (b) Hermite–Gaussian (order ) source waveforms with m and m, respectively. (c) Slit positions with and slits (d) and the widths for and .

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. The kernel based on HO in (47) giving FRFT as a special case of general LCT formulation is not numerically analyzed since LCT based system is a more general version while various combinations including FRFT systems are applicable with the formulation in Tables 1 and 3 in the Methods section. Therefore, two different set-ups with the kernels in (48) and in (6) for the inter-plane propagation are compared for the same design of set-up in terms of the properties of the slits and planes. Inter-plane distance vector is given by (cm). The distances of the first plane to the first lens and the second plane to the second lens are denoted with the vector where both the distances are set to 21 cm while the lenses of the focal length 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. and 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 m and the width of m. and 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 () and also the negative volume of Wigner function compared with FSP. The vectors of and composed of the values on the first, second and third (sensor) planes for Gaussian sources are denoted with and , respectively, for FSP while with and 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. 9a–c,e–g for FSP and LCT cases, respectively. It is observed that LCT focuses the light better on the slit locations while FSP reduces significantly. is much more improved compared with . The magnitudes of the interfering paths are shown in Fig. 9d,h, for and , 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 is also improved compared with . Observe that HG source has already negative Wigner volume of 1.076 which is much further improved by LCT set-up compared with FSP.

Figure 9

HG source with the order of and m is utilized where the resulting spatial domain waveforms on the planes with the indices (a) , (b) and (c) for FSP, and (e) , (f) and (g) for LCT. for (d) and (h) . Scaled Wigner distribution for FSP on the planes with (i) and (j) , and for LCT with (k) and (l) .

Gaussian sources

The waveforms in spatial domain for Gaussian sources are shown in Fig. 10a–c,e–g for FSP and LCT cases, respectively. and values are obtained where the magnitudes of the interfering paths are shown in Fig. 10d,h, for and , 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. and 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. 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 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.

Figure 10

Gaussian source with m is utilized where the resulting spatial domain waveforms on the planes with the indices (a) , (b) and (c) for FSP, and (e) , (f) and (g) for LCT. for (d) and (h) . Scaled Wigner distribution for FSP on the planes with (i) and (j) , and for LCT with (k) and (l) .

Open issues and discussion

There are some open issues to best exploit photonic QPC method based on MPD and FO. Mathematical formulation correlating specific set-up parameters to path magnitude distribution and negative volume of Wigner function is an open issue. Iterative formulation of the vectors , , and the matrices , and in (16), (25) and (27) are complicated as shown in the Methods section with complicated parameters and iterations in 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, 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 is represented as follows:where and 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 of the single slit. The functional form with partial sum of RTF should be calculated and summed for each combination of 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 is the wave amplitude at the point , 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, is the Green’s function vanishing on the diffraction surface for the first type of solution of Rayleigh-Sommerfeld formula, and for the monochromatic light source of wavelength as shown in Fig. 11[77]. Assuming that , the following approximation holds in rectangular coordinates:where the kernel for FSP is defined as follows: where . The kernel for Fresnel diffraction integral is obtained by further approximation of in the near-field for large z resulting in (45). This expression is the convolution integral conventionally used in phase-space optics for FSP[34].

Figure 11

Formulation set-up for first type of solution of Rayleigh–Sommerfeld diffraction through a slit [77].

Recently, scalar diffraction theory and Fresnel diffraction integral are discussed in Ref.[15] to be validly representing the evolution of light wave function modeled with the Hamiltonian of the quantized electromagnetic field as the Feynman’s path integral (FPI) solution of the quantum mechanical HO[78]. 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 on a plane is modeled to propagate into the amplitude on another plane (Eq. 16 in Ref.[15] 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, , for , is the propagation duration between the planes and is the defined equivalent mass of photon propagation. In addition, the approximated FSP kernel in (45) is simply converted to the following in 1D system:The kernel for massive particles with the mass m such as an electron is expressed as follows[6, 78]:In other words, the formulation based on phase-space optics for photon and electron propagation wave amplitudes have the similar form in (48) and (49) except an overall phase factor. The form in (48) is utilized in Ref.[33] for defining QSM while targeting only classical communications.

On the other hand, both the kernels and 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 phase-space 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 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 and are given as follows: has the same form with the propagation of light in quadratic graded-index media of having the refractive index distribution of where and are the medium parameters. The parameter matrix of the propagation through the quadratic graded-index medium of length is given by the following (Section 8.3.3 in Ref.[34]):where . There is a FRFT relation between scaled versions of the input and output with FRFT order as where and denotes the ath order FRFT of f(x). FRFT operation of order is represented with the parameter matrix of and . As a result, represents a FRFT relation between the input and output scaled with the parameter where the parameters are and while as a special case of LCTs.

Matrix formulation for HO/LCT system with Gaussian sources

The following formulation is valid for both HO and LCT based design with Gaussian sources where the corresponding iteration parameters are defined in Table 1. The elements in the vector are defined as follows:where for is given as follows:Here, the matrix multiplication denotes for any matrix for and , , and for are defined in Table 1. Assume that and define the operators creating block diagonal matrices by putting the vectors and the matrices for , respectively, (all the vectors or the matrices having the same dimensions) to the main diagonal and making zero the remaining elements. The matrix is more simplified as follows compared with the more complicated form achieved for electron based FSP in Ref.[6]:where the diagonal matrices are defined as follows: block and vector for are defined as follows: is a lower triangular block matrix defined as follows:Expanding in terms of real and imaginary parts is achieved by finding the real and imaginary parts of for and for , and and for since is a real matrix. This is easily achieved by using the explicit forms of and in Table 1. Some variables and constants used in Table 1 (not defined in the table) are the following: , denotes , , for , and as the conjugate of for real values of .

Polynomials for the example in (18-21) are presented in Table 2 for the simple case of and Gaussian source. It is possible by using the explicit modeling to make various gedanken experiments and to perform complexity theoretical calculations.

Matrix formulation for LCT system with Hermite–Gaussian sources

Iteration parameters utilized in (22–26) are presented in Table 3. Some variables and constants used in Table 3 (not defined previously or for Table 1) are the following: for , and . The formulation in (22) is obtained by using the integral equality for Hermite polynomials (Section 16.5 in Ref.[65]) in an iterative manner along the planes:The parameters and in utilized in (23) are defined as follows by using the iterations in Table 3:where the following are defined:where is defined for and while and are defined as follows:where and for , and . Finally, becomes equal to the following by using the iterations in Table 3:where is the following: