Optical Realization of Wave-Based Analog Computing with Metamaterials

: Recently, the study of analog optical computing raised renewed interest due to its natural advantages of parallel, high speed and low energy consumption over conventional digital counter-part, particularly in applications of big data and high-throughput image processing. The emergence of metamaterials or metasurfaces in the last decades offered unprecedented opportunities to arbitrarily manipulate the light waves within subwavelength scale. Metamaterials and metasurfaces with freely controlled optical properties have accelerated the progress of wave-based analog computing and are emerging as a practical, easy-integration platform for optical analog computing. In this review, the recent progress of metamaterial-based spatial analog optical computing is brieﬂy reviewed. We ﬁrst survey the implementation of classical mathematical operations followed by two fundamental approaches (metasurface approach and Green’s function approach). Then, we discuss recent developments based on different physical mechanisms and the classical optical simulating of quantum algorithms are investigated, which may lead to a new way for high-efﬁciency signal processing by exploiting quantum behaviors. The challenges and future opportunities in the booming research ﬁeld are discussed. with improved robustness with incoherent light. Kwon et designed a nonlocal metasurface by engineering nonlocality in momentum space that can perform both even- and odd-order differenti-ations. This study provides a high-quality, efﬁcient and polarization-insensitive image processing platform for 2D edge detection. We believe that as it continuously evolves, the metamaterial-based analogous optical computing will be extensively applied in signal and image processing.


Introduction
Exploring novel approaches to improve the computational capacity and efficiency is a goal that humans have been continuously pursuing. Early computers are constructed mechanically [1] or electronically [2][3][4] on the principle of analog, aiming to perform mathematical operations. Despite the impressive success achieved in the fields of weather prediction, aerospace and nuclear industry, these machines faced significant obstacles of slow response and large size [5]. In the 20th century, digital computing emerged with a rapid development of semiconductor technology and large-scale integrated circuits. It began to gradually substitute for their conventional analogue counterparts on the strength of easy programmability, high speed and flexibility. However, in solving specialized computational tasks, such as imaging-processing and edge detection, digital computers are often inefficient and hindered by high-power consumptions. As Moore's law is approaching its physical limitations, the long-abandoned analog approach as an alternative paradigm has attracted renewed attention for its potential abilities to overcome these shortcomings [6][7][8].
In addition to classical mathematical operations, researchers have also extended the concept of analog computing to quantum algorithms that promise an exponential speedup that is far beyond classical ones in solving problems of large integer factorization (Shor's algorithm) [34] and combinatorial optimization (Grover's algorithm) [35,36]. Some fundamental properties of quantum computing, including superposition principle and interference phenomena, are the essence of wave nature, which are not exclusive to quantum mechanical but are common to classical waves. By encoding quantum bit (qubits) into different degrees of freedom for the electromagnetic field (e.g., frequency, polarization, orbital angular momentum, space and time bins), many quantum computations can be efficiently simulated in optics .
In this review, we will focus on recent developments in both the physics and applications of the spatial analog optical computing. The paper is organized in four parts as described below. In Section 2 we overview the basic concept of computational metamaterials and the general principles for design and implementation for classical mathematical operations including integration, differentiation and integration equation solution. In Section 3 we introduce a type of wave-based signal processors that can mimic quantum mechanism with classical optical waves. These studies show that the quantum algorithms like Grover's search algorithm and Deutsch-Jozsa algorithm can be simulated by cascading metamaterial functional blocks. In the last section, we conclude with an overview of computational metamaterials based on different mechanisms, the main challenges and the opportunities in the field for future research.

Computational Metamaterials
In 2014, Silva et al. [102] proposed a concept of "computational metamaterials" by locally tailoring the electromagnetic parameters of the metamaterials, the elaborately designed meta-structures can reshape the spatial profile of the input signal to perform mathematical operations including spatial differentiation, integration and convolution. The basic idea of computational metamaterials is schematically illustrated in Figure 1a as an example. The stacked multi-layered structure is functional as a first-order differentiator, for arbitrary input wavefronts f 1 (y) and f 2 (y), the corresponding output profiles are proportional to df 1 (y)/dy and df 2 (y)/dy, respectively. Compared with conventional analog signal processors and Fourier optics system, metamaterial-based approach provides more flexible mechanism for manipulation, more integrable volume and subwavelength thickness. In general, there are two major protocols to realize metamaterial-based analog computing: metasurface (MS) approach and Green's function (GF) approach. Similar to the classical 4f systems, the MS approach is based on spatial Fourier transformation, with a single-layered metasurface or multi-layered meta-transmit(reflect)-array to realize the desired transfer function, sandwiched between two subblocks performing the Fourier and inverse Fourier transform (see Figure 1b). In the GF approach, by optimizing the transmitted (reflected) response of the metamaterial slabs it can also act as certain mathematical operations without involving additional Fourier lenses. The former has the advantage of implementation simplicity and the latter has a more compact size. In the following sections, we will briefly review the design principle and recent progress pioneered by these two approaches.

Metasurface Approach
Initially, let us consider a pure mathematical transformation of convolution operation in Fourier domain. For an given input function f (x, y), the corresponding output function g(x, y) is the result of a desired mathematical operator which is indicated by the Green's function h(x, y), x and y denote the spatial coordinates on the two-dimensional (2D) plane. By applying the linear convolution operation, g(x, y) can be described by Based on convolution theorem, the Fourier transform of g(x, y) is equivalent to the multiplication of functions f and h in the Fourier domain: where (F −1 {·})F {·} represents the (inverse) Fourier transform, abbreviated (I)FT. k x and k y are frequency variables in the spatial Fourier space, H(k x , k y ) is spatial FT of the transfer function h(x, y). It is actually easy to implement Equation (2) in the practical systems by simulating the input (output) function with a classical electric field E in (x, y) (E out (x, y)), in this case, the spatial variables (x, y) can play the role of (k x , k y ). In wave-based computing system, Equation (2) can be rewritten as: According to Equation (3), the optical computing system should comprise of three cascade metamaterial functional subblocks which are designed to accomplish FT, H(x, y) and IFT sequentially. In order to realize the FT subblock, an analog of traditional lenses is employed since the converging lenses can achieve 2D Fourier transform in the focal plane. In [102], a graded-index (GRIN) dielectric slab with a parabolic variation of permittivity is introduced to realize FT operator [103,104]. Besides of GRIN medium, meta-lens can also be applied to perform FT operator [105][106][107][108][109][110]. It worth nothing that, although we can define a GRIN(−) subblock to perform IFT in which the effective constitutive parameters should satisfy ε(µ) GRI N(−) (x, y) = −ε(µ) GRI N(+) (x, y), it is not feasible for natural materials or practical for metamaterials to realize simultaneously negative permeability and permittivity. Therefore, as an alternative method, we can use additional FT subblock to approximate output profile according to the relation: In this way, the electric field distribution is proportionate to the mirror image of the desired output function. Based on above discussion, the fundamental components of MS approach are represented in Figure 1b. Next, we will go to the details of how to design transfer functions when implementing different mathematical operations, such as differentiation, integration and convolution.
Firstly, for simplicity, we discuss the cases of one-dimensional nth derivation operators. Based on derivative property of the Fourier transform, d n f (y)/dy n = F −1 (ik y ) n F [ f (y)] , according to Equation (3), the expression ik y n is the required transfer function, where H(y) ∝ ik y n or H(y) ∝ (iy) n . It means that the amplitude and phase of impinging optical wave will be modulated by e ik∆ ∝ (iy) n when propagating through the transfer function subblock (suppose the thickness of transfer function subblock is ∆), where k = (2π/λ 0 ) ε ∆ (y)µ ∆ (y), with free space wavelength λ 0 , ε ∆ and µ ∆ are relative permittivity and permeability of transfer function subblock. For the cases of integration and convolution operations, the transfer functions become to H(y) ∝ (iy) −1 and H(y) ∝ sin c(y), respectively. In addition, considering the limitation of lateral dimension, the transfer function H(y) has to be normalized to ensure the energy conservation (without gain materials) and keep the reflection (transmission) coefficient blow unity. Now, the question is, how to find a proper spatially-variant meta-structure that fulfills the required electromagnetic response. Here, we divide recent achievements of the MS approach into two categories according to the construction type.

Reflective Configurations of MS Approach
Plasmonic metasurfaces [83,[111][112][113][114][115][116][117][118][119] promise the abilities to tailor the local phase and amplitude by exploiting the strong light-matter interactions due to the plasmonic resonances. In 2015, an inhomogeneous plasmonic metasurface approach, which consists of three-layer sandwich structure was introduced to optical analog computing for the first time by Anders Pors et al. [120]. A periodic arrangement of gold nanobrick arrays is on the top layer, an optically thick gold ground and subwavelength dielectric spacer is used to excite gap-surface plasmon (see Figure 2a). By properly designing the position and size of each nanobrick, this metal-insulator-metal structure can perform integration at visible wavelength with high efficiency. Similar to this work, in 2017, Chen et al. [121] used a dendritic structure instead of nanobrick to implement first-order differential operation. By varying the geometrical shapes, the desired reflection coefficients could be achieved (see Figure 2c). A different approach for reflectarray configuration is dielectric metasurface, which is composed of high refractive index dielectric nanoparticles supporting both magnetic and electric dipole-like resonances based on Mie theory. This approach can significantly reduce the ohmic losses caused by metallic structures in the optical spectrum. Motivate by this, Ata Chizari et al. [122] propose a dielectric meta-reflect-array where silicon nanobricks are placed on the silica spacer and silver substrate to manipulate am-plitude and phase of reflected cross-polarized field by varying the lengths of nanobricks (see Figure 2b).

Transmittive Configurations of MS Approach
Compared with reflective configurations, the optical computing systems with transmission mode are easily applied to optical devices, however most of them are facing the challenges of low transmission efficiency. In 2015, Amin Khavasi et al. [32] proposed a concept of "metalines" in which three symmetrically stacked graphene-based building blocks were utilized to perform differentiation and integration in the transmitted way (see Figure 3a). By manipulating plasmons surface wave via surface conductivity variation of the graphene, the amplitudes and phases of the transmitted wave could be controlled completely. While it was a proof-of-principle study, it opened up an avenue of the graphene-based structure which are designed to implement the mathematical operators with advantages of dynamic tunable and high-compact characteristics. Then, in 2016, Zhang et al. [123] introduced a more practical and flexible computing metamaterial system based on effective medium theory by drilling sub-wavelength hole-arrays with different radiuses, the desired effective permittivity distribution can be obtained to realized functional subblocks (see Figure 3b). To further improve the parallel processing abilities, a multi-way analog computing device combined with additional transformation optical subblock is investigated [124]. Two identical signals on the same input port propagate along opposite directions and perform first-and second-order differentiators simultaneously. Very recently, a one-dimensional high-contrast transmitarray (HCTA) metasurface was experimentally demonstrated by Zhang et al. [100] to realize Fourier transform and spatial differentiation with feature size of 140 nm. The on-chip meta-system provides low insertion loss and broad operating bands by tailoring widths and lengths of the rectangular slot arrays (see Figure 3c,d). Inspired by radial Hilbert transform filter, Huo et al. [125] proposed a spin-dependent dielectric metasurface, which can perform a spiral phase filtering operation in which a spiral phase profile is imprinted to the input wave. The tricky part in this proposal is to utilize the π phase difference in opposite azimuth leading to a destructive interference when processing convolution transformation. The authors demonstrate that spiral phase filtering operation is equivalent to the 2D spatial differentiation of incident light field (see Figure 3e).

Green's Function Approach
In the original GF approach proposed by Silva et al. [102], the transmission or reflection coefficient of multi-layered stacked slabs with transversely homogenous and longitudi-nally inhomogeneous are engineered to approximate the desired transfer function without getting into the Fourier domain. Since additional Fourier lenses are not required, the GF approach has advantages over the MS approach with a more compact size and an easier fabrication process. The basic idea of the GF approach is to find the proper transmittance coefficient T(k y ) of the multi-layered slab to fit the desired GF kernel G(k y ) within 0 < k y < k 0 , in which k y is the transverse wavevector, and k 0 is the free space wavevector. A fast synthesis algorithm is adopted to minimize the difference between T(k y ) and G(k y ), which can be expressed as: where M is the number of layers, w r and w i are weight coefficients for real and imaginary parts. The optimized transmission spectrum is in good agreement with the desired GF which correspond to 2th spatial derivative provided that M = 10, w r = 2, w i = 1, and relative permeability of each layer keep at unity.

Diffraction Gratings
Diffraction gratings have been widely used in processing optical signals with integration or differentiation operations, most of them are based on Fano resonance or guided-mode resonance in which the reflection or transmission coefficient can approximate the transfer function in the vicinity of the resonance. In [131], Victor A. Soifer et al. experimentally demonstrate that a guided-mode resonant grating with period TiO 2 slab waveguide on a quartz substrate can perform spatial differentiation with an obliquely incident beam. As shown in Figure 4a, the transverse profile of the incident wave P inc (x, y) is formed by zeroth order diffraction to the transmitted profile of p tr (x, y). Figure 4b is the scanning electron microscopy image of the designed diffraction grating. A similar strategy is to adopt high-contrast grating, as illustrated in Figure 4c [132]. By observing the amplitude and phase (indicated by red dot and blue line in Figure 4d) of the spatial spectral transfer function evolving with different incident angle, it is found that the transfer function has a phase variation of π around a certain angle θ 0 . Within slight deviation of θ 0 , the transfer function satisfies first-order differentiation. In a recent work, Alexios Parthenopoulos et al. [134] propose a novel dielectric subwavelength grating by employing high quality suspended Si 3 N 4 films, where the first-and second-order spatial differentiation can be implemented at oblique and normal incidence, respectively. The even guided mode is excited with normal incident wavefront, and change dramatically to the odd guided mode as increasing the incidence angle. In addition to these dielectric gratings, metallic gratings are also investigated and experimentally demonstrated to realize analogue spatial differentiators. It is well-known that when a metallic grating is illuminated by an incident wave, at certain incident angle, the surface plasmon resonance will be excited. Yang et al. [133] found that by adjusting the geometric parameters of the subwavelength gratings, in the vicinity of surface plasmon resonance, the transfer function can be tailored to a linear function which satisfies the requirements of first-order differentiation.

Plasmonic Structure
By exploiting the critical coupling condition of the surface plasmon polaritons (SPPs), many plasmonic analog computing devices have been realized. In 2017, Zhu et al. [135] introduced a classical Kretschmann prism configuration to realize the plasmonic spatial differentiation where the oblique incident wave (transverse profile of S in (k x )) with TMpolarized are used to excite the surface plasmon between glass substrate and silver film, k x is the transverse component of the wavevector (see Figure 5a). In this case, the amplitude of reflected wave (transverse profile of S out (k x )) is the result of the interference between the direct reflection and the leakage of SPPs. The spatial spectral transfer function around k x = 0 can be expressed as: where ϕ is the phase change related to the process of the direct reflection at the glass-metal interface. A = (α spp − α 1 )/ cos θ 0 and B = (α spp + α 1 )/ cos θ 0 , where α 1 and α spp are the radiative leakage rate of the SPP and the intrinsic material loss rate, respectively. When the critical coupling condition is satisfied (α 1 = α spp ), Equation (6) can be simplified and approximated as H(k x ) ∝ (e iϕ /B)ik x , which corresponds to the first-order differentiation. However, the SPPs is sensitive to the defect of the interface. As such, to further improve its robustness, Zhang et al. [137] proposed a unidirectional SPPs-based spatial differentiator which can avoid backscattering by applying an external static magnetic field (shown in Figure 5b). The magnetic field breaks the time-reversal symmetry and makes the system have a nonreciprocal property, which is critical to the unidirectional SPPs leaky mode. As shown in the bottom of Figure 5b, when applying external magnetic field, the SPPs surface mode changes from bidirectional to unidirectional. To verify the practicability when facing real-time image processing, the impact of plasmonic spatial differentiation under time-variant optical signals is analyzed by Zhang et al. [136]. By adopting similar setup described in [135], the reflected field profile with central frequency ω 0 in the spatiotemporal coordinate can be expressed as: where c s and c t are two constants where c s = e iφ cos θ 0 2a 1 phase change of direct reflection without excitation of SPP, θ 0 , a 1 , v g , v gls are the incident angle, the leakage rate, the group velocity of the SPP, and the speed of light in the prism at central frequency ω 0 , respectively. The second term of Equation (7) indicates that the output field profile should consider the change rate of input time-modulated signals. As a conclusion, the authors estimate that the processing speed of plasmonic differentiator is up to the maximum of 10 13 frame/s which is restricted by the time of establishing the SPP leaky radiation.

Two-Dimensional Dielectric Metasurfaces
The above discussions are mainly focused on 1D mathematical operations. However, anisotropic edge detection is not sufficient for real-time, high-throughput 2D image processing since multiple measurements are required. In this part, metasurfaces with symmetrically distributed dielectric nano-resonators are briefly reviewed to illustrate their capacity to reveal all information from the boundaries of objects. The basic physical mechanism of these 2D dielectric metasurfaces originates from the modal evolution where a bound state resonance mode with normal incident plane wave will change to a leaky waveguide mode at oblique incidence. The transmitted amplitudes encoded with information of spatial variation are sensitive to the transmission when increasing the incidence angle [151].
A type configuration is illustrated in Figure 6a [143], silicon nitride patch-arrays are embedded in a homogeneous medium of silicon dioxide. By tuning the height of patch-arrays, the spatial bandwidth and resolution of the edge detection can be modulated. The transmission coefficient of TE incident wave can be expressed as T TE (ϕ, k r ) ≈ α(ϕ)k 2 r , where ϕ is the azimuthal angle, α(ϕ) is the gain of the system, and k r = k 2 x + k 2 y . From Figure 6b, it can be seen that the transmission coefficient on the xoy plane formed the symmetric parabolic distribution due to the geometric symmetric, thereby it can be used to perform the second-order derivative with the quadratic approximation (see Figure 6c). In [142], a similar strategy is adopted to realize a 2D Laplace operator which can be combined with traditional imaging systems with a numerical aperture (NA) up to 0.32. However, since the quasi-guided leaky modes can only couple with p-polarized incident waves due to the modal symmetry, this scheme can only operate for one polarization. To overcome this limitation, Wan et al. [141] introduce a spatial differentiator can be performed for arbitrary polarized (unpolarized) wave by tailoring the spatial dispersion of electric dipole resonance supported by silicon nanodisks (see Figure 6d). Figure 6e shows that for both sor p-polarized incident waves, the transfer function have similar parabolic shape which indicate the functionality of spatial differentiation for arbitrary polarization. The corresponding experiments verify the performance of this scheme (see Figure 6f).

Nonlocal Metasurface
Kwon et al. [144] introduce a new mechanism with increasing the nonlocal response of metasurfaces, which is generally considered to be undesirable and detrimental, by sinusoidally modulating the permittivity of each split-ring resonator (SRR) (Figure 7a). Within the modulation, the transmission formed a Fano response where a sharp variation of transmission coefficient emerges on the resonance frequency accompanied with the incident angle (see Figure 7b). By adding a horizontally misplaced metallic wire-arrays, the requirements of breaking both vertical and horizontal mirror symmetry are fulfilled to perform first-order derivative operation (see Figure 7c). Furthermore, 2D edge detection is demonstrated with combining the 1D computing metasurface and its 90 • rotational symmetric structure (see Figure 7d-f).

Random Medium
Another counterintuitive approach is by employing random medium rather than deliberately designed structures with certain scattering properties to perform wave-based analog computing. In the study of Hougne et al. [146], it consists of two subsystems: a chaotic cavity as the random medium and a reflect-array metasurface playing the role of wavefront shaping device. As shown in Figure 8a, a plane impinging wave with input vector X is modulated by a wave front shaping device and subsequently propagating through a random medium to obtain the desired output wave front. For controlling spatial degree of freedom, the incident wave front is divided into four segments, marked by A to D (see Figure 8b). Here, an indoor cavity of irregular geometry is used as the random environment characterized with the Green's function, which plays the key role of analog computing by enabling each segment of wave front to contribute to each output point. The reconfigurable metasurface contains 88 units whose reflection coefficient can be tuned from −1 to 1. Figure 8c shows the experimental configuration in which microwave absorbers are installed to limit the reverberation of the cavity.

Inverse Design
Estakhri et al. [147] implement the discrete numerical method of solving integral equation in an optical way by using a recursive approach with a metamaterial block, feedback waveguides and coupling elements (see Figure 9a). The Fredholm integral equation of the second kind is expressed as: is the function to be solved, K(u,v) and I in (u) are integral kernel and input signal. By sampling the complex values of electric fields on the input and output plane (with N points), the corresponding N × N matrix equation can be established. To obtain the function g(u) is equivalent to calculate the inverse N × N matrix (A), which can be realized by an inhomogeneous permittivity distributed metamaterial block. In the study, the authors adopt the objective-first optimization technique to approach the desired relative permittivity in each iteration constrained by criterion of transfer function I-A, where I is the identity matrix (the simulation results is shown in Figure 9b). The experiment of reflective configuration verifies that it is feasible to apply the inverse-designed metamaterial platform in solve integral equation (see Figure 9c,d).

Brewster Effect
As mentioned above, the transversely homogenous multilayer slabs with even symmetry in the spatial Fourier domain cannot be applied to perform first-order derivative or integration operators since the odd symmetry is needed for these operations. In the study of Youssefi et al. [148], a rotated configuration is theoretically demonstrated to overcome this constraint by breaking the refection symmetry. The schematic of the rotated structure is shown in Figure 10a, when the oblique incident wave illuminating on the boundary of two dielectric medium (in this study, they are air and a dielectric substrate with different refractive indices of 1 and 2.1) at the Brewster angle, the GF around k y = 0 can be approximated with a linear function which is used to implement first-order derivative (see Figure 10b). The proposed simple structure enables effective implementation of the mathematical operations, however, there are some drawbacks in this configuration. The approximation requires the reflection spectrum of the interface to become equal to zero, which means only the signals around the Brewster angle can perform the derivation, leading to a relatively narrow spatial spectrum and small reflected energy.

Goos-Hänchen Effect
Motivated by the high sensitivity of Goos-Hänchen (GH) effect for the interface of total internal reflection, Xu et al. [149] propose a GH-based approach where the kernel transfer function is acting on the air-glass interface to perform first-order derivative. After the total internal reflection, the reflected angular spectrum contains two linearly polarized components which are converted by a quarter-wave plate into two opposite circular-polarized components. The destructive interference of the reflected waves leads to the functionality of edge detection. The experimental setup is illustrated in Figure 11a, the phase distribution and spatial spectral transfer function are shown in Figure 11b. the results of edge detection for different images are shown in Figure 11c. Figure 11. (a) The experimental setup for measurement of the spatial spectral transfer function of the first order spatial differentiator based on GH effect; (b) top: the phase distribution for the spatial spectral transfer function. Bottom: the comparison of theoretical and experimental results for the spatial spectral transfer function for k x = 0; (c) output image of edge detection with different rotating angles. Adapted with permission from [149], AIP Publishing, 2020.

Spin Hall Effect
In the study of Zhu et al. [150], the authors demonstrate that realizing spatial differentiation is a natural effect of spin Hall effect (SHE) for the reflected or refracted wave at any planar interface. In this proposal, spin-dependent transverse shifts of −δ and +δ corresponding to parallel and antiparallel spin states, the destructive interference of the opposite shifts by selecting an orthogonally polarized reflected wave lead to the final output wave function, which can be expressed as: |ϕ out = i 2 dy[ϕ in (y + δ) − ϕ in (y − δ)]|y , where |ϕ in and |ϕ out denote the input and output wave functions. The limit of |ϕ out is approximately proportional to the first-order spatial differentiation when δ is much smaller than the initial wavefunction profile. As an example, Figure 12a shows that an obliquely incident paraxial beam illuminates on the interface of two isotropic media, two orthogonal polarizers are installed between the incident and reflected waves. The transfer function around k y = 0 is depicted in Figure 12b. the experimental results for 1D edge detection are shown in Figure 12c.

Quantum Computing with Metamaterials
Quantum computation [38], by employing the principles of quantum mechanics, such as superposition and entanglement, provides the basis for quantum algorithms that enable tremendous improvements over classical computing technology for certain complex tasks. For example, the large integer factorization problems are considered practically impossible for classical algorithms as the numbers get larger than 2048 bits. Therefore, it is essential for modern RSA encryption because the codes are virtually unbreakable. However, this cryptosystem is facing a daunting challenge by Shor's algorithm [34], which promises the abilities of factoring integers in polynomial time. Another famous example is Grover's search algorithm [35,36] in which searching an unsorted database only takes O √ N time compared with classical algorithm of O(N), where N is the number of entries in the database. Although the quadratic speedup is slightly less impressive than other quantum algorithms with exponential acceleration capacities, the time saved is significant when N is considerably large. In fact, the extraordinary power of quantum computation comes from the smallest information carrier named quantum bits (or "qubits"). The quantum superposition phenomenon allows a single qubit to hold 0 and 1 states at the same time, as multiple qubits is coherently controlled, they process inherent parallel computing abilities that are far beyond the fastest classical computers. However, despite these advantages over classical computers, up to now, to realize large-scale universal quantum computer is extreme difficult. Until 2016, a proof-of-principle demonstration of Shor's algorithm can only be achieved to factor the number 15 [152]. Apparently, quantum computation still has a long way to go before it becomes practical. The biggest challenges that we faced were preparing, control and measurement of qubits, actions on qubits should be very carefully handled that even the slightest interaction with surrounding environments may leading to decoherence and destroy the quantum information. Among enormous numbers of experimental realizations of quantum algorithms in different physical systems (including nuclear magnetic resonance, quantum dots, trapped ions or QED cavities), the optical implementation is a promising way since the photons are more robust to external perturbation and thus corresponds to long decoherence times. In addition, the quantum optical system do not need extremely cold temperatures to function. The central idea of quantum optical computing is that quantum information can be encoded by different degrees of freedom of photons (e.g., polarization, orbital angular momentum, spatial and temporal modes), by utilizing the common properties shared by both classical optics and quantum mechanics, such as superposition and interference, it is possible to simulate certain quantum behaviors with classical light [37,39,40,42,43,45,[48][49][50][51][52][53][55][56][57][58][59]62,67,153].
As a new platform for arbitrarily manipulating wavefront of light, metamaterials or metasurfaces can also be applied in simulating quantum algorithms with classical wavebased optics. In 2017, Zhang et al. [154] proposed a metamaterial-based quantum algorithm analog to perform Grover's search algorithm. The designed metamaterial consisted of four cascade subblocks (an oracle subblock U m , two Fourier transform subblocks F, a phase plate subblock I − 2|0 0|) corresponding to the operators of oracle, Walsh-Hadamard transform and inversion-about-average (IAA) operation in Grover's search algorithm (see Figure 13a). In this scheme, the quantum states are directly mapped onto classical optical fields. For instance, the incident electric field amplitude "E(y)" is the analogue of quantum probability amplitude, transversal coordinate "y" is used to label the item of the database and the maximum number of the database is depending on the full width at halfmaximum of the beam "D" (listed in Table 1). When light occurs on the metamaterial, each functional subblock it passed is equivalent to performing a quantum operation; therefore, each roundtrip represents one iteration of the search algorithm. After multiple iterations in the metamaterial, the marked item is found by measuring the field distribution on the output plane (see Figure 13b). Recently, Cheng et al. [155] verified that Deutsch-Jozsa (DJ) algorithm could be simulated by using a similar strategy (see Figure 13c,d).

Quantum Classical
Items in the database |i Probability amplitude of the equivalent quantum state "y" "E(y)" The maximum number of the database N "D" In Table 1, we compare the general protocol of performing Grover's search algorithm in quantum and classical realm to clarify the differences between the two strategies. As the spatial freedom of the optical field (cbit) is adopted to simulate the qubit, the classical wavebased architecture shows a simple and enlightening way to the field of signal processing. It should be mentioned that despite that it has been proved that entanglement is not necessary for the efficiency of some quantum algorithms [156], a lack of entanglement will essentially limit the physical resources for these classical analogies, leading to an exponential growth of the width of the beam as the numbers of qubits increase. In addition, the size of the database is also limited by the spatial resolution, diffraction effect and paraxial approximation of the optical system. With the rapid development of the field of artificial electromagnetism, we can expect more novel functionalities on simulating quantum behaviors, such as three-dimensional all-optical search and data classification.

Conclusions and Outlook
In this review, we have presented some of the most exciting developments in the field of analogous optical computing taking advantage of metamaterials, offering an integrated and miniaturized platform for image processing and edge detection to overcome the limitations of digital electronic computers. With continuous exploitation, tremendous novel approaches have been proposed based on different physical mechanisms, including Brewster effect, surface plasmonic, photonic spin Hall effect, local/excitation mode coupling to implement mathematical operations. Despite the promising unprecedent advances of meta-structures over traditional optical system for analog computing, to date, it is still extremely challenging to further improve the performance which is limited by spatial resolution, efficiency, operating bandwidth, strong polarization/angle/symmetry dependence. In general, computational metamaterials can be divided into two major categories: plasmonics and dielectrics. Leveraging highly confined surface plasmonic modes, plasmonic spatial differentiators process advantages of simplicity and ultracompact size. However, they suffer from some fundamental limitations, for instance, the critical coupling condition limits the incident angles, the presence of higher-order Taylor coefficients lead to a narrow spatial bandwidth and further effect the resolution of the computational system. In addition, the conversion efficiency is restricted by intrinsic material loss. To improve the efficiency and achieve high-resolution edge detection, all-dielectric metasurface is emerging as an alternative platform, which is a promising candidate for further broadening of the operational spatial bandwidth and reducing the absorption losses. It should be mentioned that dielectric metasurfaces also have some drawbacks-some of them only work for one polarization while the energy of other polarized waves is completely wasted, which may lead to relatively low signal-to-noise ratio. Moreover, determining how to balance the features of numerical aperture, efficiency and spatial resolution requires an elaborate design [157].
Although some limitations exist in this emerging field, we are pleased to see that many achievements have been proposed and demonstrated to overcome these drawbacks. Wang et al. [158] proposed a photonic crystal differentiator with improved robustness with incoherent light. Kwon et al. [145] designed a nonlocal metasurface by engineering nonlocality in momentum space that can perform both even-and odd-order differentiations. This study provides a high-quality, efficient and polarization-insensitive image processing platform for 2D edge detection. We believe that as it continuously evolves, the metamaterial-based analogous optical computing will be extensively applied in signal and image processing.

Conflicts of Interest:
The authors declare no conflict of interest.