Efficient light propagation algorithm using quantum computers

Quantum algorithms can potentially overcome the boundary of computationally hard problems. One of the cornerstones in modern optics is the beam propagation algorithm, facilitating the calculation of how waves with a particular dispersion relation propagate in time and space. This algorithm solves the wave propagation equation by Fourier transformation, multiplication with a transfer function, and subsequent back transformation. This transfer function is determined from the respective dispersion relation, which can often be expanded as a polynomial. In the case of paraxial wave propagation in free space or picosecond pulse propagation, this expansion can be truncated after the quadratic term. The classical solution to the wave propagation requires O(NlogN) computation steps, where N is the number of points into which the wave function is discretized. Here, we show that the propagation can be performed as a quantum algorithm with OlogN2 single-controlled phase gates, indicating exponentially reduced computational complexity. We herein demonstrate this quantum beam propagation method (QBPM) and perform such propagation in both one- and two-dimensional systems for the double-slit experiment and Gaussian beam propagation. We highlight the importance of the selection of suitable observables to retain the quantum advantage in the face of the statistical nature of the quantum measurement process, which leads to sampling errors that do not exist in classical solutions.

Quantum algorithms in particular have been expedited to overcome the boundary of the problems unsolvable by classical computation [10,15,16] and also speed up classical algorithms [17][18][19].Recent accomplishments in wave physics have e.g.demonstrated quantum advantage for solving the Poisson equation [20], computing the electromagnetic scattering cross sections [21], and simulating non-classical interference phenomena [22].Some of these have demonstrated exponential speedup over the bestknown classical algorithms, highlighting the efficiency with which quantum algorithms can treat interference problems.
Among the most prominent quantum algorithms are those which leverage the quantum Fourier transform (QFT) [23,24] to efficiently convert between real space and Fourier space representations.These include quantum phase estimation (QPE) [25], and Shor's algorithm [26] among many others.From a physical point of view, Fourier transformations are a particularly powerful tool to solve linear differential equations, such as those encountered in wave propagation.Hence, solving wave propagation problems in quantum computers is an active area of research, as can be seen from recent works, consisting of wave equations under Dirichlet and Neumann boundary conditions [22] and the singleparticle Schrödinger equation with various potential barriers [27,28].Nevertheless, solving the wave propagation for modern optics remains unrealized so far.Unlike prior studies of wave propagation by quantum computers, we herein focus on solving wave propagation from the Helmholtz equation by implementing a Fourier transformation technique inspired by the classical beam propagation method (BPM).Classically, the BPM solves the Helmholtz equation by intermediate transfer of the wavefunction, which is discretized into N points in each dimension, into the Fourier domain using the Fast Fourier Transformation (FFT), the multiplication by a specific transfer function, and subsequent back transformation into real space.By implementing this procedure in quantum computers, we can exploit QFT and an efficient quadratic phase operator, introduced in this work, to solve the Helmholtz equa-Figure 1.Schematic overview of classical and quantum propagation algorithm.In both cases, to solve the Helmholtz equation, we convert the initial wavefunction from real space to Fourier space via either the classical or quantum Fourier transform, then apply the transfer function, and finally take the inverse Fourier transformation to obtain the propagated wave in real space.The operator P j is indexed by the qubit j and P jl is indexed by qubits j and l with the specific phase shift value.
tion in paraxial approximation, achieving computational complexity speed up compared to the classical BPM.
In addition, in this work, the difference between classical and quantum algorithms is demonstrated and compared.It is shown that only n = log 2 N qubits are required to store the complete wavefunction, and we propose a method to implement the transfer functions whose the dispersion relation has been expanded up to the order of p-th with O((log N ) p ) operations efficiently using a gate-based quantum computation algorithm.This implies the generalization of our QBPM being able to apply for an arbitrary light propagation.For demonstration purposes, we used Qiskit [29] to implement the QBPM and simulated it in a noiseless environment, i.e. on a quantum computer simulator.We showcase paraxial wave propagation (p = 2), in both one and twodimensional environments.That is, we simulated the double slit experiment for one-dimensional propagation and Gaussian beam propagation in two lateral dimensions.

II. CLASSICAL WAVE PROPAGATION
The propagation of electromagnetic, monochromatic waves in free space or in an isotropic media is described by the Helmholtz equation, where U (r, ω) is an arbitrary field and k(ω) is a wave vector, which is based on some dispersion relation ω = ω(k).A typical problem in wave propagation is finding the solution of the wave field U (r ⊥ ; z) given a specific wave frequency ω and a known initial value U 0 (r ⊥ ) = U 0 (r ⊥ ; z = 0) defined along the transverse coordinates r ⊥ , which are, r ⊥ = xe x in case of one transverse dimension and r ⊥ = xe x + ye y in case of two.
One common approach for solving this problem is to first take the Fourier transform of the initial field with respect to the transverse coordinates [30] as depicted in Fig. 1.In case of a problem with one transverse dimension, this is given as: where α is the transverse wave number in x direction.
To attain the propagated field in the frequency domain, we need to take the transfer function into account.If the transverse wave numbers are much smaller than the longitudinal wave numbers, i.e. α 2 ≪ k 2 , the transfer function can be expressed by the quadratic polynomial approximation of the dispersion relation, i.e.
The propagated field in the frequency domain is then given by where the first exponent is just a fixed phase and can be ignored for many problems.Likewise, if we consider the dispersion of a temporal pulse, we can expand the dispersion relation into a Taylor series and truncate after the quadratic term, if the pulse is not too short (typically sub-picosecond) and the linear term can be removed by the introduction of a co-moving frame of reference.Finally, we take the inverse Fourier transform to obtain the propagated field in space U (x; z).The generalization of this to two or more transverse dimensions is straightforward.
For a numerical solution, we discretize the wavefunction U (x) on N equidistant points in space is the discretized wavefunction.N = 2 n is usually selected as a power of two, consistent with the number of qubits n in the quantum algorithm.The Fourier transformation with the FFT algorithm requires O(N log(N )) steps and yields the discretized Fourier amplitudes Ũ (0) γ at the frequencies α γ = γ∆α, where ∆α = 2π/(N ∆x).
The solution then requires the multiplication of the transfer function for each discrete frequency H γ = H(α γ ; z, k), which requires O(N ) steps, yielding the Fourier amplitudes Ũγ (z) and a subsequent application of the inverse FFT to calculate U γ (z).The overall cost of applying the algorithm is thus O(N log(N )).
It should be noted that the algorithm does not find a perfect solution to the Helmholtz equation, because its discretization in real and Fourier spaces imply periodic boundary conditions in both spaces, yielding a non-vanishing discretization error.If an analytic solution I Ref (x γ ; z), or any other reference solution for a particular initial field is known, we can measure the impact of the error on a wavefunction's intensity using e.g. a rootmean-square-error ϵ(z): where is the absolute value square of the analytic solution and is the absolute value square of the numerically solved discrete solution.
We shall later see that the QBPM introduces random errors in the measurement process and hence ϵ(z) is a statistically distributed quantity in its own right.Its distribution can be sampled over many repetitions N sim and described by its mean error µ = ⟨ϵ⟩ and the standard error σ = ( ϵ 2 − ⟨ϵ⟩ 2 ) 1/2 .

III. QUANTUM ALGORITHM FOR PROPAGATION
To formulate the quantum algorithm for the wave propagation discussed above, we map the classical BPM steps onto quantum gates.All procedures are illustrated in Fig. 1.The step-by-step derivation for our proposed algorithm is as follows.

Selection of qubits: we create a wavefunction
|ψ⟩ from n = log N qubits, the state of which can be described by an arbitrary superposition of its computational basis states |γ⟩ with γ ∈ [−N/2, ..., 0, 1, ..., N/2 − 1], which can be represented as a binary number γ = −a n−1 2 n−1 + n−2 j=0 a j 2 j composed of its binary digits a j ∈ {0, 1} using the two's complement to indicate negative values.Mapping the N point wavefunction onto the n qubits already entails exponential enhancement of memory efficiency.
2. Initialization: We represent the initial wavefunction γ |γ⟩.An efficient and general solution to this problem is a subject of active research and not discussed in this work.However, for various special cases, some specific solutions have been shown [31].For the frequently encountered case of Gaussian beams, we can, e.g., make use of the Kitaev-Webb algorithm [32].

Fourier domain:
We perform a standard QFT on |ψ 0 ⟩ transforming the state into ψ0 = γ Ũ (0) γ |γ⟩.Note that γ now indicates the wave number index.As shown in Fig. 2, H is the Hadamard gate, followed by the phase gate R j , which leads to phase shift by the following expression; R j (ϕ) = 1 0 0 exp(2πi/2 j ) .
4. Second-order propagation: In the Fourier space, we need to multiply each Ũ (0) γ term with the transfer . This represents a phase shift, which only depends on the fixed value of ϕ and the square of the state index γ.An efficient implementation can be done by rewriting γ in its binary form.This leads to (See details in Supplementary Material S1): Note that j and l are qubit indices and that either a j = 1 or a j = 0 for all j.This allows for an efficient implementation of the transfer function operator, as we note the summands of the first sum have the value 4 j under the condition that qubit j is in the |1⟩ state and zero otherwise.Likewise, all summands of the middle sum have the value 2 j+l+1 under the condition that the qubit-pair j and l are both in their |1⟩ state and zero otherwise.The last term is for two's complement and of identical structure.The multiplication of the transfer operator Û can thus be decomposed into a series of phase gates and single-controlled phase gates, such that ψ(z) =  , where the j index indicates that it acts on the j-th qubit.Likewise, the controlled phase gates are defined using P lj (ϕ) = 1 0 0 exp(iϕ) , where the control and target qubits are respectively the l-th and j-th qubits.Note that the phases that each phase gate or controlled phase gate apply correspond to the coefficients in the two's-complement decomposition in Eq. 6.The matrix representation of these gates can be seen in Supplementary Material S2.
5. The calculation of the final state |ψ(z)⟩ is done by subsequent application of the inverse QFT (IQFT).
By measuring the final state, we collapse the field into a computational basis state with a probability of collapse being equal to the absolute value square of the wavefunction at this point in space, hence its intensity: This is equivalent to measuring a single photon in the final wavefield on a camera with N pixels.The complete intensity distribution can be sampled from a series of N s "single photon clicks", termed "shots".If the intensity distribution is not required, but some other measurable is sought-after, then a different approach for measurement may be better suited.
Note, however, that any kind of sampling induces a sampling error, with -as opposed to the discretization error -is statistical in nature.The impact of both types of error can hence be separated from an ensemble average by calculating the statistical features of the RMSE with a fixed number of shots N s over many simulations N sim .If the sample based estimator is unbiased then a mean error µ is indicative of the discretization error and the standard error σ is indicative of sampling noise.
Both QFT as well as the transfer function multiplication operation require O(n 2 ) operations, thus the entire calculations scales with O((log N ) 2 ).As such, the QBPM is exponentially faster than the classical BPM algorithm.
Although all above BPM and QBPM are derived based on the second-order approximation, it is crucial to note that our algorithm can also be applied for higher order polynomials α p , as are required to e.g.model pulse propagation in dispersive media or to go beyond the paraxial approximation for large-angle wave propagation.A general formula can be found in the Supplementary Material S1 and it can be seen that the speedup is still maintained, as the scaling for these operations is O(n p ).

IV. ONE-DIMENSIONAL PROPAGATION
We now turn to implementing the above algorithm and presenting the double-slit experiment as a test case to validate the method's functionality and to investigate the relative impact of discretization error µ, which is a generic feature of the discretized solver, vs. the sampling noise σ, which is a feature of the measurement procedure inherent in quantum computers.
The initial state was chosen as: Where the separation of the slits d is 0.5 mm, the width of each slit w = 0.1 mm, and the wavelength of the monochromatic light λ = 532 nm.The number of qubits is 15 qubits.The intensity pattern at the output is a function of z and the analytical solution is [33]: where tan θ = x/z and I 0 is the maximum intensity in the center.

A. Simulated interference pattern and Error analysis
We now compare our quantum algorithm to the classical analytic solution for the double-slit experiment.As depicted in Fig. 3, we can see the result of the quantum simulation after different z propagation distances compared to the analytical solution.The result agrees well with the expected interference pattern.
Both discretization error µ and sampling error σ depend on the propagation distance z and the number of shots N s .In Fig. 4(a), both are plotted as a function of the propagation distance.The result indicates that the discretization error increases linearly with the propagation distance, and that it is larger than the sampling error except for z = 0.This is an important finding, as it naturally gives a rise to a selection rule for the number of shots N s,max , beyond which a reduction of the shot noise will not improve overall precision, beyond the systematic discretization error.
For the given simulation, such a number can be estimated from the data shown in Fig. 4(b).While at z = 0 the total error can be reduced with the inverse of the square root of N s to an arbitrarily low value, a minimal error is evident for all z > 0. For all investigated distances this likely occurs in our simulation at N s,max ≈ 10 6 .As a complete reconstruction of the absolute value of the wavefunction by means of sampling does not appear to be a suitable mode of operation for a quantum simulator anyway, a more systematic investigation of this boundary is not undertaken.

V. TWO-DIMENSIONAL PROPAGATION
In addition to the simulation of one-dimensional propagation as demonstrated in Sec.IV, we herein extend our propagation algorithm to two-dimensional propagation and select a Gaussian beam propagation in free space as a case study.

A. Computational model
The wave function is initialized by where A 0 is the normalization constant and w 0 is the beam waist.As such, here, α and β are discretized by 2πγ/∆xN and 2πζ/∆yN , accordingly.In this example, we define the width at focal plane (w 0 ) at 0.05 m, wavelength (λ) at 532 nm, the number of qubits in each direction set to 5 qubits (10 qubits total), and the Gaussian is centered at the domain's center, such that x 0 = y 0 = 0.For two-dimensional propagation, we can straightforwardly add one more coordinate into the circuit as the algorithm is separable along x and y directions.This also applies to non-separable input states, which, however, have not been demonstrated in this work.This implies that all steps remain the same as those of the one-dimensional case, hence the quantum algorithm will be duplicated for the other axis as depicted in Fig. 5.It should be noted that the dimensional independence is a feature of the dimensional independence of diffraction from the quadratic polynomial approximation and not of the quantum algorithm.
The Gaussian beam is propagated with the propagating distance (z) as illustrated in Fig. 6.It is obvious that when z r is larger, the Gaussian beam broadens, as expected.

B. Observable and Error analysis
As opposed to the 1D-case, we here aim to highlight the impact of a well-selected observable on the number of shots N s required for a good estimation of some observable.As a test case, we estimate the beam waist w(z) as shown in Eq. 12 as a function of the propagation distance, (12) where x i and y j run from 0 to the position domain ranged N .Where at each position (x i , y j ) we have measured the number of times the circuit collapses into this final state P ij = N ij /N s .We compare it with the value obtained from classical discretized solution as shown below, (13) where |U ij | is the absolute value of the wavefunction.As opposed to the last section we here reference against the discrete approximation; hence the difference ϵ = w Q − w ref is exclusively caused by the sampling error.Again,  we measure the standard error σ w = ( ϵ 2 − ⟨ϵ⟩ 2 ) 1/2 .Both values depend on the number of shots N s and the propagation distance z and have been determined from N sim =100 independent simulations.
We now analyze the size of the standard error as a function of the propagation length and the number of shots.This data is displayed in Fig. 7.The data indicates that only N s = 100 already allows for a very good sampling precision of the beam waist radius.If a better approximation is required then the number of shots can be increased, and we observe an increase in precision of 1/ √ N s .Larger propagation distances inherit the more considerable error, consistent with the 1D case.This shows that the fact, that quantum computers require sampling, does not systematically ill-affect the efficiency of the algorithm; if a sufficiently well-behaved measureable is required, which we argue is a common use case, a very moderate number of shots is sufficient for a good estimation of the measurable.
Nonetheless, as the propagation in the 2D case is more sophisticated than that in 1D, we also benchmarked our QBPM among four different cases, consisting of noisy quantum simulator, noiseless quantum simulator, classical discretized algorithm, and classical analytic solution to guarantee the accuracy of our algorithm.All of these comparisons agree well as can be seen in Supplementary Material S3.

VI. CONCLUSION
We have developed a quantum beam propagation algorithm based on the paraxial approximation.We show that our quantum algorithm can efficiently simulate beam propagation in both one and two dimensions by implementing only phase gates and controlled phase gates with an exponential speedup over the classical solution, based on the fast Fourier transformation.This allows for a low number of qubits and operations used for the propagation when compared with the classical counterpart.We also carry out a rigorous error analysis to give more insights into the role of the sampling error, which is a new source of errors, when quantum computers are used.We analyze two scenarios: the sampling of the absolute value squared of the wavefunction and find that the sampling error falls below the discretization error after only a few thousand shots N s .In a second scenario, we only analyze a scalar measurable, in this case, the waist of a broadening gaussian beam.Here we find that a very good  accuracy can be achieved with only a few hundred shots and hence the quantum advantage can be maintained for a broad range of scenarios.
As the 2D Helmholtz equation can have factorized and non-factorized initial states, the demonstration for this initialization problem can fulfill this important gap.In principle, this will allow one to simulate a system bounded by Dirichlet or Neumann boundary conditions where its initial states are not separable.
We further argue that our work does not merely add one more algorithm to the toolbox of quantum simulation.Instead, it bridges the gap between Fourier optics and quantum computation and lays the foundation for quantum computers as a tool to efficiently simulate more complex problems in photonics, such as Boson sampling or spontaneous conversion in photonic nanostructures.

Figure 2 .
Figure 2. One-dimensional QBPM logic circuit for the double-slit experiment.

Figure 3 .
Figure 3. Interference pattern sampled from Ns = 100, 000 shots and 15 qubits as a result of simulating the double-slit experiment using the QBPM in the 1D case compared to the analytical solution with different propagation distances.

Figure 4 .
Figure 4. Mean error µ (points) and standard error σ (bars) are accumulated from Nsim = 100 simulations.(a) As a function of propagation distance z where the number of shots Ns=100,000.(b) As a function of Ns at different propagation distances.

Figure 5 .
Figure 5. Two-dimensional QBPM logic circuit for Gaussian beam propagation in free space where the upper part of the circuit is for propagation in x direction and the lower part is for y direction.

Figure 7 .
Figure 7. Sampling error (measured by the standard error σw) of the measurable wQ, i.e. the radius of the gaussian wave packet, determined over Nsim = 100 as a function of the propagation distance zr and a function of the number of shots Ns.Even for longer propagation lengths, wQ is a very good approximation of wc