1. Introduction

Quantum computation [1] has significant advantages in speed over classical computation since the ability of parallel processing by exploiting features of quantum mechanics as superposition of quantum states and quantum entanglement. A typical example is the Deutsch-Jozsa algorithm (DJ algorithm), which is proposed by David Deutsch for the first quantum algorithm in 1985 [2] and continuously improved by Richard Jozsa [3], Cleve, Ekert, Macchiavello and Masca [4] et al. While the DJ algorithm is a proof-of-principle study designed purely to exhibit the great potential with speed-up of quantum computation over the classical algorithm, it lays the foundation for later practical quantum algorithms such as Shor’s algorithm [5] and Grover’s search algorithm [6]. Due to its importance in the field of quantum computation, leading hitherto to an enormous number of experimental realizations in different physical systems such as QED cavities [7,8], Nuclear Magnetic Resonance [9,10], quantum dots [11], trapped ions [12] and quantum optical systems [13,14]. However, in order to put quantum computing into practice, it requires a number of quantum bits (qubits) in a superposition state which is difficult to prepare and coherently control. Although it is highly desirable to realize such significant quantum algorithms, the extreme requirements for performing the DJ algorithm in quantum regime obstruct it far from practical applications.

Since the inherent wave nature shared by both classical optics and quantum mechanics, such as superposition principle and interference effect, it is possible to simulate certain quantum algorithms with classical light. For instance, by encoding quantum bits into different degrees of freedom for the electromagnetic field (frequency, polarization, orbital angular momentum, space and time bins), many quantum computations can be efficiently simulated in optics, such as DJ algorithm [15–18], Grover’s search algorithm [19–22] and quantum walks [23,24].

Recently, taking advantage of the unprecedented flexibility of metamaterial, Silva et al [25,26]. proposed the concept of “computational metamaterials”. Metamaterial refers to a kind of artificial microstructure not found in nature, which provides additional degree of freedom to arbitrarily manipulate the propagation of light beams with locally redefined effective medium [27–33]. The computational metamaterials can perform various mathematical operations, for example, spatial differentiation, integration and convolution, etc. We note that, many kinds of artificial structures have been designed [34–39], such as the grating [40–43], plasmonics [44], photonic crystal slabs [45], all-dielectric metasurface [46–49] and a single optical interface with spin-Hall [50] or Brewster effect [51]. This inspired us that it is also feasible to realize a certain quantum algorithm with the help of metamaterials.

In this paper, we propose and experimentally demonstrate a metamaterial-based system which can simulate DJ algorithm. Three different metamaterial blocks with either constant or balanced phase modulation in the oracle subblock are used to test the validity of the scheme. In this case, whether the marked function being constant or balanced, we can directly determine it by measuring the output electromagnetic field with only one time. Furthermore, three-dimensional polyjet printing technology [52] is introduced to guarantee the processing finenesses of the metamaterials with micron class resolution. Theoretically, the general protocol of this strategy can be applied in multiple spectra including microwave, terahertz and optical region. Such a metamaterial-based DJ algorithm simulator may lead to remarkable achievements in integrable and miniaturized wave-based signal processors.

2. Scheme to simulate Deutsch-Jozsa algorithm with metamaterials

We consider a binary function $f \colon \{0,1\}^{n} \rightarrow \{0,1\}$ with $N=2^{n}$ arguments of $y$, where $n$ is an integer to give the number of input bits. Provided that either $f(y)=0$ or $f(y)=1$ runs for all arguments, the function f can be regarded as a constant function. Otherwise, if the value of f is 0 for N/2 input arguments and 1 for the other N/2 arguments, f is called balance function. The DJ algorithm is to determine what kind of the function that f is. For a classical Turing machine, in the worst scenario, it required to access and compute the function values at least N/2+1 times to solve this problem, more than half capacity of the database. Consequently, when N is extremely large, the computation will spend vast resources and time. In another word, classical algorithms are usually “ergodic”. Nonetheless, the DJ algorithm requires only one step to certainly determine the function type by using quantum superposition principle and interference effect, providing an exponential speedup as compared to classical algorithms.

The general strategy for simulating the DJ algorithm with metamaterials is graphically shown in Fig. 1. Firstly, to simulate the DJ algorithm with classical optics, a possible method is required for preparing classical analogue for the n qubits input state. In our protocol, a system of n-qubit in superposition states $\Psi _{1}$ can be simulated with a classical electric field “E(y)”.

 

Fig. 1. Metamaterial-based scheme to simulate DJ algorithm. The whole metamaterial model (width of w) consists of two functional subblocks. The oracle subblock (length of $d_{o}$) with the gradient color encodes the detecting function f(y) into the input states $\Psi _{in}$ by assigning a phase shift on each spatial position along y-aixs. The Fourier subblock (length of l) with blue color acts as Fourier transformer. The output profile indicates the function type.

Download Full Size | PDF

The transverse positions “y” are used to label the different states (input arguments) [53,54].The Gaussian beam profile on the position “y” can be understood as the analog of the probability amplitude of the equivalent quantum state [19,55]. In this configuration, the maximum capacity of quantum states depends on the width of the incidence beam “D”, however, due to the diffraction effect, the spatial resolution of the phase modulation can also restrict the actual usage of the quantum states. It is worthy to note that the speed up effect of DJ algorithm is originated from wave interference instead of entanglement, the later is not necessary for the efficiency of quantum DJ algorithm [56].

Then, the detecting function f(y) is encoded into the input states by assigning a phase modulation on each spatial position “y” along the input transversal direction. The refractive index of the first subblock is designed to achieve 0 and $\pi$ phase distribution depending on the value of the function f. When the incident wave illuminates on the first subblock (the area with gradient color in Fig. 1), the phase of the transimitted wave will be modulated with the factor being $k_{0}n(y)d_{o}$. Here, $k_{0}=2\pi /\lambda _{0}$ is the vaccum wave vector, $d_o$ is the thickness of the first subblock and n(y) is the effective refractive index at position “y”. Equivalent to the original quantum state $\Psi _{1}$, the electromagnitic field, which has been marked by the first subblock, can be transformed into $\Psi _{2}$:

From optical point of view, the fourier transform is the key operator for parallel speedup based on the essence of the wave interference. Here, we omit the first Hadamard transform which is used in the original DJ algorithm to yield the superposition states, the reason is that the incidence beam profile already stands for the $2^{n}$ superposition transverse modes of a single photon. We use the spatial freedom of optical field which can be regarded as cbits (c for classical) instead of qubits. Since the operators such as phase shifter and Fourier transformer acting on cbits will not leading to decoherence which is a universal essential problem in true quantum computation, the concept of fidelity may not apply to this scheme. For clarity, the differences of general protocol for performing DJ algorithm in quantum and simulating in classical realm are listed in Table 1.

Table 1. Comparison of the quantum realization and the classical simulation for the DJ algorithm.

View Table

Based on the discussion above stated, the designed metamaterials shall be comprised of two cascaded functional subblocks. The first part is an oracle subblock $U_{o}$, which can mark the function f(y) to the incident electric field by imprinting a spatially-dependent phase factor $exp(-i\Delta \varphi )=1$ or -1 $(\Delta \varphi =0,\pi )$ on the beam. In this way, we can modulate the symmetry of the function by changing the phase distribution. The second part is a Fourier transform subblock, a graded-index (GRIN) dieletric slab with a parabolic variation of permittivity is introduced to realize the Fourier transform operator [25].

To facilitate the simulation, we set permeability $\mu _{o(g)}=\mu _{0}$ for the whole system. In this case, the effective permittivity of $U_{o}$ with $\Delta \varphi$ being 0 and $\pi$ can be expressed as:

3. Numerical results

In order to qualitatively demonstrate the feasibility of our strategy to simulate the DJ algorithm, we propose an elaborately designed metamaterial with the effective permittivity varying continuously along the y-axis according to the Eq. (3) and Eq. (4). For the constant function, the oracle subblock should provide a same phase shift along the y-axis. Here, we choose the identical phase factor to be $\pi$ in the following simulations (shown in Fig. 2(a)). As for the balance function, in priciple, the oracle subblock should be filled up equally with 0 and $\pi$ phase modulation. Here, we restrict the oracle subblock to form an anti-symmetrical phase distribution, in this way, it can ensure the wavefront of $\Psi _{2}$ is balanced (with both equal phase and amplitude distribution) under the Gaussian distribution incident wave. To verify our hypothesis, two kinds of balance functions, which correspond to simple and complex cases, are designed.

 

Fig. 2. Theoretical simulations of metamaterial-based DJ algorithm analogy with continuously varied permittivities. The incident waves $E_{in}(y)$ with Gaussian profile propagates throughout the metamaterial where three different detecting functions are applied: (a) constant, (d) simple balance and (g) complex balance. The oracle subblock (with thickness $d_{o}$) is comprised of two permittivity areas which provide 0 (yellow) and $\pi$ (gray) phase shift on transmitted waves. The Fourier transform subblock is represented by graduated green. The simulation results of snapshots of the field intensity in the xy-plane and output electric fields for three different detecting functions are illustrated in (b),(e),(h) and (c),(f),(i), respectively.

Download Full Size | PDF

For the simple case (shown in Fig. 2(d)), the oracle subblock processes anti-symmetrically distributed permittivity with respect to the x-axis and is able to form a $\pi$ phase difference between $y\geq 0$ and $y<0$ regions. The two phase elements indicate 2 quantum states (1 qubit), in one way, it can be recognized simply as the simulation of the Deutsch algorithm. When the incident wave is propagating through the first subblock, the phase distribution of the wavefront is in a form of anti-symmetric balance function about x-aixs. To further verifying our scheme can be used to simulate the DJ algorithm, we divide the oracle subblock into 6 parts (about 2 qubits) where regions of $-86\,mm<y<-46\,mm$, $2\,mm<y<42\,mm$, $90\,mm<y<100\,mm$ process $\pi$ phase variation and $-100\,mm<y<-90\,mm$, $-42\, mm<y<-2\,mm$, $46\,mm<y<86\,mm$ process 0 phase variation. As shown in Fig. 2(g), since the total areas of 0 and $\pi$ phase modulation are the same, the wavefront is still marked by a balance function. The electric fields with phase modulated by the oracle subblock will pass through the GRIN subblock, which can perform the Fourier transformation. Then, the type of the function can be directly distinguished by a single measurement of the output wave profile $E_{out}(y)$, where the zero (maximum) intensity at the focal plane accompanies with the balanced (constant) function marked by the oracle subblock.

We perform finite element method based 2D numerical simulations by using COMSOL Multiphysics. Perfectly matched layers technique is applied at truncated boundary area to simulate the open boundary condition. Two ports are placed at the front and rear end of the metamaterial, one is acting as the excitation source with Gaussian beam profile $E_{in}(y)=exp(-y^{2}/40 mm^{2})$ propagating along the +x direction, and the other serves as a receiver. The geometrical parameters of the metamaterial are $d_{o}=28$ mm (the length of oracle subblock), $l=228$ mm (the length of GRIN subblock) and $w=204$ mm (the width). Figure 2(b), Fig. 2(e) and Fig. 2(h) present the snapshots of the electric field intensity throughout the metamaterial with the detecting function being constant, simple balance and complex balance, respectively. The corresponding field distribution at the output plane are shown in Fig. 2(c), Fig. 2(f) and Fig. 2(i). It can be clearly seen from Fig. 2(c) that the incident wave is primarily focused on the output port. The maximum intensity at $y=0$ position means the encoded function y is constant. On the other hand, the center intensity is zero, as shown in Fig. 2(f) and Fig. 2(i). This indicates that the oracle subblock carries balance function. Since the 0 and $\pi$ phase elements correspond with different effective permittivities, it makes the oracle subblock processes spatially varying impedance which results in different transmissions. Therefore, the intensity of the output electric field shown in Fig. 2(f) is not symmetric.

In addition, we analyze the 3-qubit system for balance function with 8, 10, 12 phase regions of the oracle subblock. The top parts of Fig. 3 show the three configurations for the oracle operator. From the simulated field intensity of Fig. 3(a), Fig. 3(b) and Fig. 3(c), it is found when the oracle subblock is divided into more phase elements, the interval between two peaks of the output electric field will increase and clutters will arise (see Appendix A for details). The separation of two peaks can be interpreted as the result of the Fourier transform. The oracle subblock loads a sinusoidal-like oscillatory signal (with equally distributed $e^{-i0}(+1)$ and $e^{-i\pi }(-1)$ phase factors) to the input Gaussian waveform, the more amount of phase elements correspond to the higher oscillation frequency. When we perform the Fourier transformation (through the GRIN subblock) on the oracle-modulated wave-front, the interval of the peaks will become wider as the oscillation frequency being increased. The number of simulated quantum states $N$ can be further increased by improving the spatial resolution of the phase modulation to make the oracle subblock contains more phase pixels, or by adopting two-dimensional structures such as dielectric metasurfaces to utilize the whole information of the beam.

 

Fig. 3. The simulation results for 8 (a), 10 (b) and 12 (c) phase elements of the oracle subblock. Top: the phase distributions of three balance functions. the yellow and gray areas correspond to phase shift of 0 or $\pi$, respectively. Middle: the field intensity in the metamaterial. Bottom: the field intensity of the output plane.

Download Full Size | PDF

Here, as a demonstration, a practical metamaterial based on the effective medium theory is designed by drilling sub-wavelength hole-arrays with different radii on a single piece of dielectric layer. The benefit of adopting this method is that the phase difference can be changed continuously and ease of integration and fabrication. The unit cell of the designed metamaterial is shown in Fig. 4(a), an air hole with diameter being R in the center is drilled on a dielectric host medium of relative permittivity being $\varepsilon _{sub}=3$. The period and height for the unit cell are chosen to be $a=4$ mm and $h=30$ mm, respectively. R is the diameter of the hole and can be adjusted depending on the effective permittivities. According to the effecitve medium theory [57], the operating wavelength should be large enough with respect to the size of the unit cell to avoid the bragg scattering. Here, the wavelength of the incident wave is $\lambda _{0}=15$ mm. The effective permittivities of the unit cell can be expressed as:

 

Fig. 4. Perspective view of unit cell and the full-wave simulation of our designed metamaterial. (a) Diagram of a unit cell with $a=4$ mm, $h=30$ mm and R being adjustable. (b) Simulated transmittance $T$ (red line) and phase $\varphi$ (blue dots) responses of the designed meta-atom for the oracle subblock (the length is 7a) with variation of diameter R from 0 to 4 mm. The 0($2\pi$) and $\pi$ phase responses correspond to $R=2.06$ mm, $T=0.81$ and $R=3.5$ mm, $T=0.81$, respectively. The top view of the metamaterials with the actual structural design (c),(e),(g) and the simulation results of the field intensity (d),(f),(g) for constant, simple balance and complex balance functions.

Download Full Size | PDF

The top view of our designed discrete models based on air holes in all-dielectric metamaterials are shown in Fig. 4(c), Fig. 4(e) and Fig.  4(g), which correspond to the case of constant, simple balance and complex balance functions, respectively. Full-wave simulations are performed to further assess and prove the reliability of the proposed method. The field distributions in xoy-plane clearly show that there are excellent agreements between the actual structural design and the ideal medium, as shown in Fig. 4(d), Fig. 4(f) and Fig. 4(h).

4. Experimental results

To further validate our design, experimental measurements are carried out. Figure 5(a) shows the experimental setup with a linearly polarized standard-gain horn antenna (working in Ku band) being placed 50 mm away from the sample. In this case, the incident wave with near Gaussian profile can be realized. The antenna gain is 24.3 dB at 19.93 GHz and the polarization direction is perpendicular to the xoy-plane. To ensure the polarization orientation perpendicular to the metamaterial, top and bottom surfaces of the samples are covered with copper foil. The near field distribution of transmitted wave is measured using a small monopole antenna which is positioned 15 mm away from the output plane of the sample. The monopole antenna is installed on a two-dimensional scanner driven by computer program to improve the sampling accuracy. Since the antenna will continuously scan the field intensity during the test, it cannot be tightly attached to the sample. We use an Agilent PNA-X 5242A network analyzer to generate and receive the signals. An absorbing plate with an aperture, which has the same size with the cross-section of the sample, is used to prevent unwanted electromagnetic scattering being received by the antenna.

 

Fig. 5. Experimental verification of metamaterial-based DJ algorithm analogy. (a) The experimental setup with a 2D near field scanning system. (b), (c), (d) The photographs of 3D-printed metamaterial samples. The zoom-in photographs show the details of the holes of diameter $R=3.5$ mm for constant fucntion (b) and $R=3.5,2.91,2.06$ mm which corresponding phase variation from 0 to $\pi$ for balance functions (c), (d). The zoomed red area indicate the GRIN subblock with a group of graduated apertures distributed along y-aixs. (e),(f),(g) Simulated (blue dots) and measured (red line) results of the electric field intensity which is normalized to a range from 0 to 1.

Download Full Size | PDF

The proposed three metamaterial-based samples are fabricated by a 3D printer (Objet 1000 Plus) with “High Quality” option in which a print resolution of up to $300\times 300\times 1600$ dpi is achieved and the total consuming time is 19.47 hours. $\textrm {VeroClear}^{TM}$ RGD810 is chosen as modeling material with 30 mm-thick where the dielectric constant and loss tangent being 3.0 and 0.01, respectively, at the operating wavelength $\lambda _{0}=15$ mm. $\textrm {Support}^{TM}$ SUP705 is used as the supporting material to connect the platform and samples. Each sample processes the same size with $256\times 204\times 30 \ mm^{3}$. The photographs of the 3D-printed metamaterial samples with the detecting function being constant, simple balance and complex balance are shown in Fig. 5(b)-Fig. 5(d), respectively. The main difference among the three samples is the variation of diameter for air holes in the oracle subblock. The amplifying images of Fig. 5 show the details of the subwavelength structures. It is clearly shown that the same air holes are distributed within the whole oracle subblock for the constants function. Differently, there are two sizes of hole-arrays locating symmetrically with respect to the x-axis for the balance function. Using the fabricated samples, we measured the transmitted wave for constant, simple balance and complex balance functions, as shown in Fig. 5(e) -Fig. 5(g), respectively. The output profiles at $y=0$ reach a maximum value indicate the detecting function is constant, on the contrary, a near-zero intensities indicate the balance functions. It worth nothing that, the background material used in the metamaterial model is of relative low dielectric constant, such that the whole system is highly transparent without a multiple of internal reflection. Therefore, the steady-state field distribution can be regarded as the result of one round propagation of the incident wave throughout the system. In order to make experimental results be comparable with simulated solutions, the field intensities are normalized to range from 0 to 1. In this case, it is clearly shown that a good consistency between the measured (red line) and simulated (blue dots) spatial distribution of the output electric field plane is achieved. The little deviation may result from the fabricating precision, diffraction effects, oblique incidence and minor variation for the incident wave profile compared to the ideal condition. When the incidence wave is oblique, the position of maximum (zero) intensity will slightly deviate from the central point at y=0, but it has little impact on the profile of the field intensity. Furthermore, due to the design of metamaterial structure is based on two-dimensional effective medium approximation, therefore, The deviation of designed radius causing by the fabrication variation may affect the performance of Oracle and GRIN subblocks. For Oracle subblock, the changes of diameter will make the phases of transmitted wave not equal to 0 ($\pi$) and it will change the focal length for GRIN subblock, both of which will impact the accuracy of the final results. Based on the numerical and experimental results, it can be concluded that the proposed strategy with all-dielectric metamaterials can simulate the DJ algorithm at microwave frequency.

5. Conclusion

In summary, we have proposed a wave-based scheme for simulating DJ algorithm by using all-dielectric metamaterials. The numerical simulations and experiment results demonstrate the viability of our protocol to simulate DJ algorithm properly by a single measurement. Our strategy further expands the metamaterial application fields of analog computation in classical optics. In addition, the metamaterial samples are fabricated through three-dimensional polyjet printing technology, which has advantages in fabricating complicated 3D subwavelength structures. Our design may have various applications involving all-optical analog computing devices, wave-based signal processing and communication system. Future research will be worthy to realize programmable or reconfigurable metamaterial system which are able to simulate quantum algorithms, such as Generalized DJ algorithm [58] and Grover’s algorithm.

Appendix

A. The influence of diffraction and impedance mismatch on the simulation of DJ algorithm

In our metamaterial-based scheme, due to the 0 and $\pi$ phase pixels possess different effective refractive indexes, the corresponding impedances are different. The different degrees of the impedance mismatch exist not only between the oracle subblock and the air, but also between the oracle subblock and the GRIN subblock. Therefore, scattering and diffraction in the course of wave propagation can not be ignored. When the oracle subblock is divided into more phase pixels (see Fig. 3), more regions of space will exist impedance mismatch, thus the intensity distribution near the origin (y=0) in the focal plane is not ideal zero but with clutters. The maximum deviation occurs in the situation of 12 phase pixels (Fig. 3(f)). The impedance mismatch problem can be solved, if the permeability is properly tailored to improve impedance characteristics. In addition, to reduce the influence of the optical diffraction, we choose at least 3 unit cells as a phase pixel (0 or $\pi$) and use a transition unit cell ($\pi /2$) between different phase pixels.

B. The influence of aberrations on the simulation of DJ algorithm

For optical system with large aperture size, the result can not be simply understand as the zero mode of the Fourier transform, the deformations introduced by high order aberrations must be considered, especially for the cases of symmetric balance function. In the main text, because the excitation source has a Gaussian beam profile $E_{in}(y)=exp(-y^{2}/40 mm^{2})$, we restrict our discussion on the situations of oracle subblock with anti-symmetrical phase distribution to ensure the wavefront is balanced both in phases and amplitudes. For symmetrical phase distribution, owe to the maximal amplitude of the wavefront (around $x=0$ area) have the same phases, the intensity distribution of the output port will be similar to the results of the constant function. Since the power of the incidence beam concentrate within the beam waist (FWHM of the incident intensity $\left | E(y) \right |^{2}$ equals 7.45 mm), it can reduce the high order aberrations impact on the final results. In order to further study the influence of the aberration, we analyze the intensity of the output plane under different width of the incidence beam “D”. Figure  6(b)-Fig. 6(d) plot the simulation results when increasing “D” from 100 mm (FWHM of 11.8 mm) to 1000 mm (FWHM of 37.2 mm). It is shown that the clutters will arise as the beam width increases, but the amplitudes at the focal point remain near-zero which indicate the input functions are balanced.

 

Fig. 6. The simulation results for different incidence wave. (a), (b) and (c), (d) are the simple and complex balance function models excited by the incidence beam of$E(y)=exp(-y^{2}/100 mm^{2})$ and $E(y)=exp(-y^{2}/1000 mm^{2})$, respectively. The output field intensities are normalized to range from 0 to 1.

Download Full Size | PDF

C. Evaluating the performance of the simulation of DJ algorithm

We define a parameter $\eta =\sum _{y}E_{focal}^2/\sum _{y}E_{total}^2$ to quantify the performance of the metamaterial, where $\sum _{y}E_{focal}^2$ is the intensity around the focal spot (two periods of the unit cell) and $\sum _{y}E_{total}^2$ is the total intensity of the output plane. We calculate $\eta$ with various numbers of the phase region marked by the oracle subblock, corresponding to the different function. It is found that $\eta$ is $5.97\%$, $4.16\%$, $4.20\%$, $1.91\%$, $5.61\%$, $5.52\%$ with 2, 4, 6, 8, 10, 12 phase regions and $61.27\%$ for constant function. In this work, we choose at least 3 unit cells as a phase region (0 or $\pi$) and use a transition unit cell ($\pi$/2) between different phase regions, the maximum number of phase regions is limited to 12. If the width of the metamaterial further enlarges, the number of phase regions cannot grow infinitely since the system is no longer satisfy the requirement of paraxial approximation and the deformations introduced by high order aberrations will affect the accuracy of the metamaterial system. It is note that the simple balance function (2 phase region) processes a large $\eta$. As shown in Fig. 2(f), the 3dB bandwidth of adjacent peaks may occupy the area around focal spot and cause a relatively large $\eta$. With increasing the number of phase regions, the impact of impedance mismatch and diffraction will heightened and can also lead to large $\eta$. In this case, it may infer that the function is balanced when $\eta <6\%$ and is constant with $\eta >60\%$ for this metamaterial model.

Funding

National Natural Science Foundation of China (11204218, 11674248, 11674266, 11774057, 11874285, 61505164); Hong Kong Scholars Program (XJ2017006).

Acknowledgments

Kaiyang Cheng and Weixuan Zhang contributed equally to this work. The authors are grateful to Prof. Lijun Jin for his valuable advice. This work is supported by the Central National Natural Science Foundation of China (NSFC) (Nos. 11774057). Z. Wei would like to acknowledge financial support from the NSFC (Nos. 11874285, 11674248, 11204218). Y. Fan would like to acknowledge financial support from the NSFC (Nos. 11674266, 61505164) and the Hong Kong Scholars Program (XJ2017006).

Disclosures

The authors declare no conflicts of interest.

References

1. M. A. Nielsen and I. L. Chuang, Quantum Computation and Quantum Information (Cambridge University Press, Cambridge and New York, 2000).

2. D. Deutsch, “Quantum theory, the church-turing principle and the universal quantum computer,” Proc. R. Soc. Lond. A 400, 97–117 (1985). [CrossRef]  

3. D. Deutsch and R. Jozsa, “Rapid solution of problems by quantum computation,” Proc. R. Soc. Lond.439, 553–558 (1992).

4. R. Cleve, A. Ekert, C. Macchiavello, and M. Mosca, “Quantum algorithms revisited,” Proc. R. Soc. Lond. A 454(1969), 339–354 (1998). [CrossRef]  

5. P. W. Shor, “Algorithms for quantum computation: discrete logarithms and factoring,” Proc. 35th Annu. Symp. on Foundations Comput. Sci. (Ed. S. Goldwasser,) Los Alamitos, CA (1994).

6. L. K. Grover, “Quantum computers can search arbitrarily large databases by a single query,” Phys. Rev. Lett. 79(23), 4709–4712 (1997). [CrossRef]  

7. H.-F. Wang, S. Zhang, and Y.-F. Zhao, “Implementing deutsch-jozsa algorithm with superconducting quantum interference devices in cavity qed,” Int. J. Theor. Phys. 48(8), 2384–2389 (2009). [CrossRef]  

8. S.-B. Zheng, “Scheme for implementing the deutsch-jozsa algorithm in cavity qed,” Phys. Rev. A 70(3), 034301 (2004). [CrossRef]  

9. J. Fitzsimons, L. Xiao, S. C. Benjamin, and J. A. Jones, “Quantum information processing with delocalized qubits under global control,” Phys. Rev. Lett. 99(3), 030501 (2007). [CrossRef]  

10. O. Mangold, A. Heidebrecht, and M. Mehring, “Nmr tomography of the three-qubit deutsch-jozsa algorithm,” Phys. Rev. A 70(4), 042307 (2004). [CrossRef]  

11. T. F. Watson, S. G. J. Philips, E. Kawakami, D. R. Ward, P. Scarlino, M. Veldhorst, D. E. Savage, M. G. Lagally, M. Friesen, S. N. Coppersmith, M. A. Eriksson, and L. M. K. Vandersypen, “A programmable two-qubit quantum processor in silicon,” Nature 555(7698), 633–637 (2018). [CrossRef]  

12. S. Debnath, N. M. Linke, C. Figgatt, K. A. Landsman, K. Wright, and C. Monroe, “Demonstration of a small programmable quantum computer with atomic qubits,” Nature 536(7614), 63–66 (2016). [CrossRef]  

13. M. Hor-Meyll, D. S. Tasca, S. P. Walborn, P. H. S. Ribeiro, M. M. Santos, and E. I. Duzzioni, “Deterministic quantum computation with one photonic qubit,” Phys. Rev. A 92(1), 012337 (2015). [CrossRef]  

14. S. M. Lee, H. S. Park, J. Cho, Y. Kang, J. Y. Lee, H. Kim, D.-H. Lee, and S.-K. Choi, “Experimental realization of a four-photon seven-qubit graph state for one-way quantum computation,” Opt. Express 20(7), 6915–6926 (2012). [CrossRef]  

15. B. Perez-Garcia, M. Mclaren, S. K. Goyal, R. I. Hernandez-Aranda, A. Forbes, and T. Konrad, “Quantum computation with classical light: Implementation of the deutsch-jozsa algorithm,” Phys. Lett. A 380(22-23), 1925–1931 (2016). [CrossRef]  

16. B. Perezgarcia, J. Francis, M. Mclaren, S. K. Goyal, R. I. Hernandezaranda, A. Forbes, and T. Konrad, “Implementation of deutsch and deutsch-jozsa algorithms with classical light,” in SPIE Optical Engineering + Applications,.

17. N. Johansson and J.-Å. Larsson, “Efficient classical simulation of the deutsch-jozsa and simon’s algorithms,” Quantum Inf. Process. 16(9), 233 (2017). [CrossRef]  

18. H.-R. Wei and J.-Z. Liu, “Deterministic implementations of single-photon multi-qubit deutsch-jozsa algorithms with linear optics,” Ann. Phys. 377, 38–47 (2017). [CrossRef]  

19. G. Puentes, C. L. Mela, S. Ledesma, C. Iemmi, J. P. Paz, and M. Saraceno, “Optical simulation of quantum algorithms using programmable liquid-crystal displays,” Phys. Rev. A 69(4), 042319 (2004). [CrossRef]  

20. T. W. Hijmans, T. N. Huussen, and R. J. C. Spreeuw, “Time and frequency domain solutions in an optical analogue of grover’s search algorithm,” J. Opt. Soc. Am. B 24(2), 214–220 (2007). [CrossRef]  

21. W. Zhang, K. Cheng, C. Wu, Y. Wang, H. Li, and X. Zhang, “Implementing quantum search algorithm with metamaterials,” Adv. Mater. 3, 1703986 (2018). [CrossRef]  

22. B. Perez-Garcia, R. I. Hernandez-Aranda, A. Forbes, and T. Konrad, “The first iteration of grover’s algorithm using classical light with orbital angular momentum,” J. Mod. Opt. 65(16), 1942–1948 (2018). [CrossRef]  

23. J. Yi, S. N. Burokur, G.-P. Piau, and A. de Lustrac, “Coherent beam control with an all-dielectric transformation optics based lens,” Sci. Rep. 6(1), 18819 (2016). [CrossRef]  

24. S. K. Goyal, F. S. Roux, A. Forbes, and T. Konrad, “Implementation of multidimensional quantum walks using linear optics and classical light,” Phys. Rev. A 92(4), 040302 (2015). [CrossRef]  

25. A. Silva, F. Monticone, G. Castaldi, V. Galdi, A. Alu, and N. Engheta, “Performing mathematical operations with metamaterials,” Science 343(6167), 160–163 (2014). [CrossRef]  

26. A. Sihvola, “Enabling optical analog computing with metamaterials,” Science 343(6167), 144–145 (2014). [CrossRef]  

27. X. Luo, “Subwavelength artificial structures: Opening a new era for engineering optics,” Adv. Mater. 31(4), 1804680 (2019). [CrossRef]  

28. Y. Liu and X. Zhang, “Metamaterials: a new frontier of science and technology,” Chem. Soc. Rev. 40(5), 2494–2507 (2011). [CrossRef]  

29. C. M. Soukoulis and M. Wegener, “Past achievements and future challenges in the development of three-dimensional photonic metamaterials,” Nat. Photonics 5(9), 523–530 (2011). [CrossRef]  

30. Y. Fan, N.-H. Shen, F. Zhang, Z. Wei, H. Li, Q. Zhao, Q. Fu, P. Zhang, T. Koschny, and C. M. Soukoulis, “Electrically tunable goos-hänchen effect with graphene in the terahertz regime,” Adv. Opt. Mater. 4(11), 1824–1828 (2016). [CrossRef]  

31. Y. Fan, F. Zhang, N.-H. Shen, Q. Fu, Z. Wei, H. Li, and C. M. Soukoulis, “Achieving a high-q response in metamaterials by manipulating the toroidal excitations,” Phys. Rev. A 97(3), 033816 (2018). [CrossRef]  

32. S. Chen, Z. Li, W. Liu, H. Cheng, and J. Tian, “From single-dimensional to multidimensional manipulation of optical waves with metasurfaces,” Adv. Mater. 31(16), 1802458 (2019). [CrossRef]  

33. Y. Fan, N.-H. Shen, F. Zhang, Q. Zhao, H. Wu, Q. Fu, Z. Wei, H. Li, and C. M. Soukoulis, “Graphene plasmonics: A platform for 2d optics,” Adv. Opt. Mater. 7(3), 1800537 (2019). [CrossRef]  

34. Y. Hwang and T. J. Davis, “Optical metasurfaces for subwavelength difference operations,” Appl. Phys. Lett. 109(18), 181101 (2016). [CrossRef]  

35. W. Zhang, C. Qu, and X. Zhang, “Solving constant-coefficient differential equations with dielectric metamaterials,” J. Opt. 18(7), 075102 (2016). [CrossRef]  

36. W. Wu, W. Jiang, J. Yang, S. Gong, and Y. Ma, “Multilayered analog optical differentiating device: performance analysis on structural parameters,” Opt. Lett. 42(24), 5270–5273 (2017). [CrossRef]  

37. P. del Hougne and G. Lerosey, “Leveraging chaos for wave-based analog computation: Demonstration with indoor wireless communication signals,” Phys. Rev. X 8(4), 041037 (2018). [CrossRef]  

38. Y. Wu, Z. Zhuang, L. Deng, Y. Liu, Q. Xue, and Z. Ghassemlooy, “Arbitrary multi-way parallel mathematical operations based on planar discrete metamaterials,” Plasmonics 13(2), 599–607 (2018). [CrossRef]  

39. H. Kwon, D. Sounas, A. Cordaro, A. Polman, and A. Alù, “Nonlocal metasurfaces for optical signal processing,” Phys. Rev. Lett. 121(17), 173004 (2018). [CrossRef]  

40. D. A. Bykov, L. L. Doskolovich, E. A. Bezus, and V. A. Soifer, “Optical computation of the laplace operator using phase-shifted bragg grating,” Opt. Express 22(21), 25084–25092 (2014). [CrossRef]  

41. L. L. Doskolovich, D. A. Bykov, E. A. Bezus, and V. A. Soifer, “Spatial differentiation of optical beams using phase-shifted bragg grating,” Opt. Lett. 39(5), 1278–1281 (2014). [CrossRef]  

42. Y. Fang, Y. Lou, and Z. Ruan, “On-grating graphene surface plasmons enabling spatial differentiation in the terahertz region,” Opt. Lett. 42(19), 3840–3843 (2017). [CrossRef]  

43. Y. Fang and Z. Ruan, “Optical spatial differentiator for a synthetic three-dimensional optical field,” Opt. Lett. 43(23), 5893–5896 (2018). [CrossRef]  

44. Y. Hwang, T. J. Davis, J. Lin, and X.-C. Yuan, “Plasmonic circuit for second-order spatial differentiation at the subwavelength scale,” Opt. Express 26(6), 7368–7375 (2018). [CrossRef]  

45. C. Guo, M. Xiao, M. Minkov, Y. Shi, and S. Fan, “Photonic crystal slab laplace operator for image differentiation,” Optica 5(3), 251–256 (2018). [CrossRef]  

46. S. Abdollahramezani, A. Chizari, A. E. Dorche, M. V. Jamali, and J. A. Salehi, “Dielectric metasurfaces solve differential and integro-differential equations,” Opt. Lett. 42(7), 1197–1200 (2017). [CrossRef]  

47. A. Chizari, S. Abdollahramezani, M. V. Jamali, and J. A. Salehi, “Analog optical computing based on a dielectric meta-reflect array,” Opt. Lett. 41(15), 3451–3454 (2016). [CrossRef]  

48. A. Cordaro, H. Kwon, D. Sounas, A. F. Koenderink, A. Alu, and A. Polman, “High-index dielectric metasurfaces performing mathematical operations,” Nano Lett. 19(12), 8418–8423 (2019). [CrossRef]  

49. Z. Wang, T. Li, A. Soman, D. Mao, T. Kananen, and T. Gu, “On-chip wavefront shaping with dielectric metasurface,” Nat. Commun. 10(1), 3547 (2019). [CrossRef]  

50. T. Zhu, Y. Lou, Y. Zhou, J. Zhang, J. Huang, Y. Li, H. Luo, S. Wen, S. Zhu, Q. Gong, M. Qiu, and Z. Ruan, “Generalized spatial differentiation from the spin hall effect of light and its application in image processing of edge detection,” Phys. Rev. Appl. 11(3), 034043 (2019). [CrossRef]  

51. A. Youssefi, F. Zangeneh-Nejad, S. Abdollahramezani, and A. Khavasi, “Analog computing by brewster effect,” Opt. Lett. 41(15), 3467–3470 (2016). [CrossRef]  

52. S.-D. Zhao, A.-L. Chen, Y.-S. Wang, and C. Zhang, “Continuously tunable acoustic metasurface for transmitted wavefront modulation,” Phys. Rev. Appl. 10(5), 054066 (2018). [CrossRef]  

53. R. J. C. Spreeuw, “A classical analogy of entanglement,” Found. Phys. 28(3), 361–374 (1998). [CrossRef]  

54. R. J. C. Spreeuw, “Classical wave-optics analogy of quantum-information processing,” Phys. Rev. A 63(6), 062302 (2001). [CrossRef]  

55. N. Bhattacharya, H. B. van Linden van den Heuvell, and R. J. C. Spreeuw, “Implementation of quantum search algorithm using classical fourier optics,” Phys. Rev. Lett. 88(13), 137901 (2002). [CrossRef]  

56. B. P. Lanyon, M. Barbieri, M. P. Almeida, and A. G. White, “Experimental quantum computing without entanglement,” Phys. Rev. Lett. 101(20), 200501 (2008). [CrossRef]  

57. P. Halevi, A. A. Krokhin, and J. Arriaga, “Photonic crystal optics and homogenization of 2d periodic composites,” Phys. Rev. Lett. 82(4), 719–722 (1999). [CrossRef]  

58. D. Qiu and S. Zheng, “Generalized deutsch-jozsa problem and the optimal quantum algorithm,” Phys. Rev. A 97(6), 062331 (2018). [CrossRef]