Quantum interference enables constant-time quantum information processing

We found that quantum interference allows processing large sets of data faster and more accurately than computer algorithms.


This PDF file includes:
Kravchuk polynomials are related to Kravchuk matrices that are generated by a binomial expressions of the form (1 + x) N −j (1 − x) j , with positive integer N and j = 0, . . . , N . Such an expression expands to a set of N +1 polynomials a 0,j +a 1,j x+a 2,j x 2 +. . .+a i,j x i +. . .+a N,j x N Their coefficients can be used to construct a square (N + 1) × (N + 1) matrix K (N ) with an entry in ith row and jth column given by K (N ) i,j = a i,j . Thus As an example, let us examine the expression(1+x) N −j (1−x) j for N = 3 and j = 0, . . . , 3.
Kravchuk polynomials are defined by the rows of K (N ) . Their domain is x = 0, . . . , N .
They fulfill the following orthogonality relation where n and m are integers and δ m,n is the Kronecker delta equal 1 where n = m and 0 otherwise.
In case of an unsymmetric binomial expression, the Kravchuk matrix can be generated like where 0 ≤ p ≤ 1 and the Kravchuk polynomials fulfill the following orthogonality relation Kravchuk polynomials belong to the family of special functions and can also be expressed in terms of a Gauss hypergeometric function 2 F 1 [a, b; c; z] (see Section 4) n (x, N ) are polynomials of degree n and that the parameter n and variable x can be exchanged in the following way Kravchuk polynomials can be regarded as a discrete and finite counterpart of Hermite polyno-

Kravchuk functions
n are Kravchuk polynomials that are normalized and centered at the maximum of a binomial distribution N  Fig. S1. ravchuk polynomials . and the variables n and i can be used interchangeably

Symmetric K
Interestingly, Kravchuk functions are solutions of finite oscillator wave equation where H (N ) is a finite-difference operator identified as a discrete Hamiltonian of this system. In the limit of N → ∞ Kravchuk functions tend to the harmonic oscillator wave functions ψ n (x) (Hermite-Gauss polynomials) shows Kravchuk polynomials and Kravchuk functions for N = 5.

Kravchuk transform
The Kravchuk transform (KT) is defined by means of Kravchuk functions. KT converts an input sequence (x 0 , x 1 , . . . , x S ) into a new string (X 0 , X 1 , . . . , X S ) in the following way (cf. the main text) where α is the fractionality of the transform and p = sin 2 (πα/4). Effectively, it decomposes the input string in the basis of Kravchuk functions. The KT can also be seen as multiplication of the input vector by a rescaled Kravchuk matrix K (N ) with additional phase terms.
The parameter α modifies the Kravchuk transform (Eq. S16) in the following way 1. for α = 1, it is a forward transform, 2. for α = 2, it is an inverse (backward) transform, allowing one to compute the original sequence (x l ) from (X k ) using the same algorithm or a physical system, 3. for α = 3, the transform negates the input, X k = −x k , 4. for α = 0 or for α = 4, the transform is an identity operation, X k = x k .
Similar properties are observed for the (integral) Fourier transform.
In addition, the Fourier transform and the KT can be seen as circular rotations of the data in the time-frequency space by an angle θ = πα/2. This rotation is also well-defined for 0 < α < 1, leading to an α-fractional transforms. The KT is additive with respect to α, i.e.

Discrete Fourier transform
The discrete Fourier transform (DFT) as well as the Fast Fourier Transform (FFT) algorithm decompose input data string in the basis of plane waves Since plane waves e −i2π kl S+1 are periodic and their domain consist of negative and positive integers and zero, the discrete Fourier transform defined in Eq. S17 correctly decomposes only samples of periodic data. Since the Fourier transform of a discrete periodic signal with period of length S + 1 results also in a discrete periodic output of the same length, it is enough to process one period of the input string, (x 0 , . . . , x S ), and store one period of the output sequence (X 0 , . . . X S ). However, one should remember that effectively the input sequence is seen by this algorithm as infinite . . . , x 0 , x 1 , . . . , x S , x 0 , x 1 , . . . , x S , x 0 , x 1 , . . . , x S , . . . Basis states for a 16-point discrete FT Plots (A)-( P ) depict 16 basis states exp{−i2πk/(S + 1)} for S = 15 and k = 0, . . . , S. Blue circles denote real, while red triangles -imaginary components of the states. KT versus DFT. The left column depicts exemplary input data (x l ) for l = 0, . . . , 10, while the middle and right columns -the results of computation of |X k | 2 using the DFT and the Kravchuk transform, respectively. and so is the resulting sequence . . . , X 0 , X 1 , . . . , X S , X 0 , X 1 , . . . , X S , X 0 , X 1 , . . . , X S , . . .
As a result of this, two shifted input sequences, e.g.

Image processing
In computer image processing, image moments are vectors or matrices which describe interesting properties of the source data. They can be obtained e.g. by decomposition of pixel intensities in the basis of orthogonal functions. Kravchuk polynomials have been found to be extremely useful for this purpose, as they are well defined on a finite domain of raster images and the computed moments carry out a lot of information about characteristic image features. In practice, computing Kravchuk moments is equivalent to performing the Kravchuk transform of the input image. This method has been already shown to be useful in optical character recognition (15), autonomous reading of the sign language, hand signature discrimination, writer identification, automatic face and gesture recognition (32) as well as facial expression and gait analysis. Importantly, this approach allows one to extract both local and global features of the input data and to recognize images corrupted with noise or tilt and possessing change in e.g. facial expression (32).
The α-fractional KT has been shown to be very useful in data watermarking schemes (watermark insertion and detection) (33). It allows to produce results which are invariant to image rotation, scaling and translation. For these reasons, the KT can be used in e.g. detection of copy-move forgery of images.

Search engines
The vectors obtained by the decomposition of the data in the basis of the Kravchuk functions have been also successfully tested in search engines. Such descriptors are especially able to capture sharp changes in the source data. By using only low-order Kravchuk moments one is able to perform e.g. search for similarly-looking 2-and 3-dimensional shapes (34), detect types of objects seen in a radar signal, perform automatic classification of videos, still images (including medical images) and sound files. Kravchuk transform has been also proposed as a method of efficient lossy compression of images, similarly to the discrete cosine transform used e.g. in JPEG file format.

Medical image recognition
In biology, Kravchuk image moments have been tested for determination of drugs in human plasma microscopic images and prediction of phosphorylation sites in cells. This approach has been thoroughly compared with the state-of-the-art methods towards various medical imaging applications (35). It was chosen as the best performing for breast mammography images, where it allowed to identify benign and malign masses with 90% accuracy, compared to 81% offered by the other techniques. This scheme was proven to outperform also other algorithms in analysis of computer tomography and ultrasound scans towards recognition of liver and prostate tumors.

Medical image reconstruction
KT has been already tested as a potential replacement of the FFT in generation of diagnostic ultrasound and magnetic resonance images (MRI). In these systems, the data regarding patients' bodies are collected in the frequency domain (k-space) and next, are processed by the inverse Fourier transform to obtain a real-space image. Tests were performed with MRI data coming from open repositories of brain and knee examinations. The images were reconstructed with the Kravchuk, Zernike, Pseudo-Zernike, Fourier-Merlin, Legendre and Chebyshev kernels (15, 35).
It has been pointed out that only the Kravchuk and Chebyshev transforms are discrete and allow to operate in the original Cartesian image coordinates. The Kravchuk-based method presented the best behavior in most test cases, giving the smallest reconstruction error and the highest peak signal-to-noise ratio as the moment order increased thus, is best suited for processing of high resolution data.
To assess the advantage of the KT for the MRI diagnostics, we performed similar steps to the ones presented in (35). Fig. S5 shows a comparative numerical study of the KT and FFT for a "pirate" test image and a brain scan from the OASIS database. Both source figures have resolution of 512 × 512 pixels. The source image of the brain is made of k-space raw data from an MRI system, while the "pirate" was originally prepared in the real space and next converted to the k-space using the NumPy numerical library to keep both data sets in the same form. A third data set consisted of the original brain scan but truncated to 256×256 values in the k-space by removing (zeroing) the higher-frequency components.
Subsequently, a 1% white Gaussian noise was added in parallel to all three data sets (still in the k-space). No further distortion was applied to any of the figures. Finally, the k-space images were transformed to the real space (reconstructed) with the FFT and KT.
In our test, the FFT produced artifacts, while some details (which could be tumor cells) were missing. This is best captured by the structural similarity index (SSIM). For the "pirate" it was 0.73 (FFT) and 0.92 (KT), and for the 512 × 512 brain it was 0.84 and 0.98, respectively.
In case of the 256 × 256 brain figure, SSIM achieved was 0.72 for the FFT and 0.91 for the KT. The mean square error (MSE) was ten times smaller and the peak signal-to-noise ratio (PSNR) was 10 dB larger for the KT than for the FFT. The FFT led to degradation of the usable resolution from 1-2 mm per voxel to over 5 mm. Our findings confirmed the results of the previous research (35).
Example of FFT and KT image processing Two 512 × 512-pixel test images, a "pirate" (A) and a brain scan, the latter in a form of a raw k-space data from the OASIS database (D), were used. First, the "pirate" image has been transformed into the k-space with the FFT algorithm to keep both inputs in the same form. Additionally, a third set of data was created by truncating the brain raw data to 256 × 256 values by removing higher-frequency components (G). Next, all the images in the k-space were supplemented with a 1% additive white Gaussian noise, and reconstructed with corresponding inverse transforms to model the operation of an MRI analysis.

Theory: multi-photon Hong-Ou-Mandel interference
We will now analyze in a detailed manner a generalized multi-photon HOM effect. As explained in the main text, we consider two interfering modes a and b on a beam splitter device with a tunable reflectivity r (defined as the probability of reflection of a single photon). As the input states we take photon number (Fock) states

The Schwinger representation
One may represent the su(2) Lie algebra in terms of the annihilation and creation operators of the harmonic oscillator -the Schwinger representation. For a single spin two independent oscillators a and b are required. The spin operators are then constructed in the following way S 0 is the Casimir operator S 0 (S 0 + 1) = S 2 x + S 2 y + S 2 z . The spin components fulfill the standard su(2) commutation relations

Beam splitter
Interference of two independent modes a and b on a beam splitter is governed by the following H 0 is the free quantum oscillator energy and H BS -the beam splitter interaction (19). ϕ is the phase difference between the reflected and transmitted fields behind the beam splitter. H 0 commutes with H BS .
Using the Schwinger representation, we express H in terms of the spin operators S 0 , S x , S y , S z The Hamiltonian generates the evolution operator The evolution in the Heisenberg picture allows to establish a linear relation between the input (a, b) and the output (a r , a t ) annihilation operators The relation takes the following matrix form where U BS U † BS = 1 and U † BS = U BS −1 hold true. U 0 amounts to a global phase.
We now substitute sin θ 2 = √ r and cos θ 2 = √ 1 − r to relate the evolution directly to the beam splitter reflectivity This brings us to the following relation between the input and output creation operators, to be used in the next section

Photon number amplitude
Let the input states in modes a and b be the Fock states |l and |S − l , respectively. Then, Let us substitute m + n = k to change the summation variables. Then |m + n, S − m − n = |k, S − k and the ranges of k and m are as follows The probability amplitude of detecting k and S − k photons behind the beam splitter provided that l and S − l were injected into it is The inner sum over m in Eq. S54 is a hypergeometric series. In order to simplify it, the identities from Section 4 are used. The four cases below (A-D) correspond to different summation ranges.
For simplicity, let us assume that l ≤ S − l, i.e. l ≤ S 2 .
the empty set.
Case C: min{l, k} = k and max{0, k Case D: min{l, k} = l and max{0, k + l − S} = k + l − S. This implies l ≤ S − l ≤ k. To compute the sum, the following substitution is used: Summarizing, the inner sum in Eq. S54 equals S k 2 F 1 −l, −k; −S; 1 r under the assumption that l ≤ S 2 . The probability amplitude can be rewritten into the following form The photon number statistics behind the beam splitter is given by the probability p S (k, l) = |A S (k, l)| 2

Kravchuk transform
The α-fractional Kravchuk transform of an input sequence x n = f (ξ n ), where n = 0, 1, . . . , N and ξ n = (n − N/2), is defined as follows (4) (cf. Eq. 5.2) We used the following relations (36) k (p) n (n , N ) = (−1) n N n p n 2 F 1 −n, −n ; −N ; 1 as well as the fact that the Kravchuk functions are orthonormal Now we turn A S (k, l) shown in Eq. S87 to the form of Eq. S92 where r = sin 2 θ 2 .
In specific, if we take ϕ = − π 2 and rearrange terms

Quantum Kravchuk transform on a beam splitter
Let us send a superposition S l=0 x l |l, S − l to a BS. The superposition amplitudes encode the sequence (x 1 , . . . , x S ) to be transformed. We will compute the probabilities of detecting |k and |S − k photons behind the BS It is clear now that multi-photon interference on a beam splitter followed by photon-counting detection implements α = 2θ π -fractional QKT of the input probability amplitudes where |X k | 2 are experimentally determined photon number statistics for k = 0, . . . , S.

Gauss hypergeometric function
Definition. The Gauss hypergeometric function is a special function defined with the following hypergeometric series where a, b and c are parameters, z is an argument and (x) k is the Pochhammer symbol In general, all 2 F 1 arguments and the parameter may be complex, a, b, c, z ∈ C however, within this note the arguments are always integer, a, b, c ∈ Z and the parameter is real, z ∈ R.
Properties. The Pochhammer symbol can be expressed as a division of factorials The form of Eq. S108 implies that the arguments a and b can be swapped In case of a negative a or b, the infinite sum in Eq. S108 is truncated because (x) k = 0 if x is a negative integer and k > −x. Let us assume that a < 0 and b ≥ 0 ∨ b < a. Then, let m = −a Moreover, for the same assumptions as in case of Eq. S113, the following transformation can be used to change z to 1 − z [NIST Digital Library of Mathematical Functions, 15.8.7] Identities analogous to Eqs. S113 and S114 are also valid for negative b and a ≥ 0 ∨ a < b, due to Eq. S112.
Finally, the following Pfaff's hypergeometric transformation is valid for any a, b, c and z

Characterization of the setup
In order to estimate transmission losses, we performed Klyshko efficiency measurements on the setup. In a Klyshko measurement with one SPDC source and binary detectors, one counts single events C A , C B from either output channel and coincidence clicks C AB between both channels and defines the Klyshko efficiencies η A and η B and vice versa. For low pump powers, these Klyshko efficiencies show a linear pump power dependency, and their intercept is a measure at zero pump power of total transmission efficiency (including both propagation and detection losses) of the associated spatial mode (37).
We pumped each of our SPDC sources, one at a time, with the variable beam-splitter in position 50 : 50, at successively lower power values. The resulting four-mode correlated photon statistics were then transformed into binary "photon(s)/no-photon" datasets to emulate standard binary detectors such as avalanche photo-diodes, and we determined the total efficiencies of the heralding modes to be η 1 = 50.3% and η 4 = 48.5%. The beam-splitter modes, carrying each a 3 dB loss from the splitter itself and an additional 1 dB due to splitter insertion loss and fiber-to-fiber coupling loss, exhibit a total efficiency of η 2 = 21.6% and η 3 = 20.6%. Taking into account the additional optical elements in the splitter modes, the efficiencies are consistent.
We account for the transmission losses of approximately 50% ≈ 3 dB with 1 dB initial fiber in-coupling loss due to spatial mode mismatch, 0.25 dB from imperfect detectors, and the rest from three FC/PC fiber-to-fiber couplers per mode as well as bending losses in the transmission fibers between the experimental setup and the detectors. 1 + V HOM is consistent with this result. From this, we can infer an effective Schmidt mode number of K = 1 g (2) −1 = 1.16, (38) i.e. both of our SPDC sources are close to being singlemode.
The TES detectors used in the experiment were thoroughly characterized with quantum tomography methods (20). Their quantum efficiency is above 90%.
6 Analysis of the experimental data 6.1 HOM visibilities The second-order visibility exceeding the classical value of 50% certifies quantum nature of the HOM interference and thus, the fractional QKT. The visibility is computed with the following formula (39) v (2) = n max − n min n max + n min (S117) where n max and n min are the maximal and minimal number of events registered by the TES detectors for the given photon number S.
The obtained values are gathered in Tab. S1. For S = 5 it was always greater than 50%.

Computation of probability distributions and estimation of errors
Experimental demonstration of two-mode multi-photon HOM interference requires collecting photon-number statistics, which are then compared with theoretical probability distributions.
The statistics result from multiple measurements performed with the setup depicted in Fig. 1B in the main text. The heralding modes (A & D) inform about the input state fed into the variable BS and together with the output modes are measured by highly efficient photon counting TES detectors. Thus, each measurement results in a 4-tuple consisting of the number of photons registered by TES 1−4 , denoted as (n 1 , n 2 , n 3 , n 4 ) and corresponding to photon-number states in modes A-D (29). In a single run, the SPDC source produces input Fock states consisting of up to approximately 10 photons with probability governed by the pump power (see the Materials and Methods section in the main text). The detectors register all possible values of n i ∈ [0, 10], Table S1. Second-order interferometric visibilities in HOM interference. i = 1, . . . , 4. The automation software stores this data in a database and assigns the number of events to each possible tuple. During a single 400-second run, approx. 10 9 data points are collected.
In order to obtain a photon-number statistics for a given r = sin 2 θ 2 and input Fock state a post-processing is required. The database is searched for a given pair (n 1 , n 4 ) which determines the two-mode Fock state at the BS input. Then, only records fulfilling the condition n 1 + n 4 = n 2 + n 3 are selected as they may correspond to the case of no losses in all paths. For the given (n 1 , n 4 ) the individual probabilities are computed as p S (k, n 1 + n 4 − k) = N (n 1 , k, n 1 + n 4 − k, n 4 ) S(n 1 , n 4 ) where N (n 1 , n 2 , n 3 , n 4 ) denotes the number of events of registering the given 4-tuple, S(n 1 , n 4 ) = n 1 +n 4 m=0 N (n 1 , m, n 1 + n 4 − m, n 4 ) is the total number of contributing data points and k as well as n 1 + n 4 − k are the photon numbers registered at the BS outputs. The full probability distribution consists of n 1 + n 4 + 1 values for k ranging from 0 to n 1 + n 4 .
For the TES detectors, due to the overlap between the outcomes associated with neighboring photon numbers, an |n state results in a value of n ± 1, where n is registered with probability over 0.9 and the probabilities of n − 1 and n + 1 are below 0.1 with p(n − 1) p(n + 1). Therefore, the absolute error of a single measurement ∆n = ±1. As the computation of probability is based on S(n 1 , n 4 ) data points, the measurement uncertainty equals ∆p = |∆n| S(n 1 , n 4 ) ≈ 1 S(n 1 , n 4 ) The data post-processing and error estimation was done with a Python script, which prepared input files for the Asymptote plotting software. The probability distributions for an ideal system were computed with Eq. S88. Factorials and binomial coefficients were approximated with the standard lgamma(n) function.

Realistic theoretical model
Actual experimental results (Fig. 3 in the main text) were compared with an enhanced realistic theoretical model which allowed to assess the imperfections of the system. The model includes the following parameters: average photon numbers at the outputs of both SPDCs, strength of the fiber coupling, losses in heralded and interfering modes as well as efficiencies of individual TES detectors.