**Characterization and Optimal Design of Silicon-Rich Nitride Nonlinear Waveguides for 2** μ**m Wavelength Band**

#### **Zhihua Tu 1, Daru Chen 2, Hao Hu 3, Shiming Gao 1,4 and Xiaowei Guan 3,\***


Received: 9 October 2020; Accepted: 12 November 2020; Published: 15 November 2020

**Abstract:** Optical communication using the 2 μm wavelength band is attracting growing attention for the sake of mitigating the information 'capacity crunch' on the way, where on-chip nonlinear waveguides can play vital roles. Here, silicon-rich nitride (SRN) ridge waveguides with different widths and rib heights are fabricated and measured. Linear characterizations show a loss of ~2 dB/cm of the SRN ridge waveguides and four-wave mixing (FWM) experiments with a continuous wave (CW) pump reveal a nonlinear refractive index of ~1.13 <sup>×</sup> <sup>10</sup>−<sup>18</sup> <sup>m</sup>2/W of the SRN material around the wavelength 1950 nm. With the extracted parameters, dimensions of the SRN ridge waveguides are optimally designed for improved nonlinear performances for the 2 μm band, i.e., a maximal nonlinear figure of merit (i.e., the ratio of nonlinearity to loss) of 0.0804 W−<sup>1</sup> or a super-broad FWM bandwidth of 518 nm. Our results and design method open up new possibilities for achieving high-performance on-chip nonlinear waveguides for long-wavelength optical communications.

**Keywords:** four-wave mixing; nonlinear figure of merit; silicon-rich nitride; ridge waveguide; conversion bandwidth

#### **1. Introduction**

The capacity of the current optical communication systems operating at the mainstream wavelengths, i.e., around 1310 nm and 1550 nm is approaching a shortage due to the explosive growth in information being transmitted. In order to exploit new frequency resources for telecommunications, transmitting signals at mid-infrared (MIR) wavelengths has been put on the agenda [1–6], among which the 2 μm wavelength band (1900 to 2100 nm) attracts special attention as signals in this band can be significantly amplified by thulium-doped fiber amplifiers (TDFAs) [7] and glass fiber in this band is still transparent [8]. Indeed, besides the amplifier, many other functional devices necessary for the 2-μm-band optical communication system have been also commercially available or laboratory developed including the laser [9,10], modulator [11], low-loss hollow-core fiber [12], wavelength-division multiplexer [13,14], detector [15,16], etc. Meanwhile, nonlinear devices like the wavelength converter are also indispensable and, in order to leverage the manufacturing scalability for such devices, it is preferred to use the on-chip nonlinear waveguides with the fabrication processes compatible with the mature silicon-based complementary metal-oxide-semiconductor (CMOS) technology.

As to the CMOS-compatible nonlinear waveguides, silicon (Si) waveguides and stoichiometric silicon nitride (Si3N4) waveguides are currently the workhorses at the telecommunication wavelengths [17–21]. However, Si still suffers from the two-photon absorption (TPA) below the wavelength 2200 nm [22], despite that linear Si waveguides can possess a low loss in the 2 μm band [23] and the Si nonlinear waveguides have been used for parametric conversion around 2 μm [24]. Indeed, key nonlinear devices like the parametric amplifier [25] or the frequency comb generator [26] have been demonstrated using the Si nonlinear waveguides in the MIR wavelengths but just over the wavelength 2200 nm. For Si3N4, its low linear refractive index (RI) inevitably gives a large footprint of any Si3N4-based photonic device and the low nonlinear RI usually produces a poor energy efficiency of devices based on the Si3N4 nonlinear waveguides. Moreover, in spite of demonstrations of high-performance Si3N4 waveguides fabricated by using some advanced processes [27], strain issues are still present in Si3N4 film, especially in the thick Si3N4 film, which are yet necessary for achieving compact devices at long wavelengths like 2 μm and for waveguides with proper dispersion. Besides, titanium dioxide (TiO2) waveguides [4] and silicon germanium (SiGe) waveguides [28] were also demonstrated to transmit optical signals at 2 μm, but, up to now, no reports can be found for them being used for nonlinear processes in this band.

Recently, silicon-rich nitride (SRN) has emerged as a promising candidate for the on-chip nonlinear optics [29–39] since SRN possesses larger linear and nonlinear RIs from a Si excess compared to the stoichiometric Si3N4 and, more importantly, the SRN film can have a substantially reduced stress [40]. While the SRN waveguides have been applied for many nonlinear processes like the octave-spanning supercontinuum generation (SCG) [31] and the frequency comb generation (FCG) [32], they were almost operating at the near infrared wavelengths. A SRN waveguide was very recently used to transmit optical signals at 2 μm and the wavelength conversion was also demonstrated with a pulse pump in the 2 μm band, but the waveguide was quite thin (300 nm) and the dispersion was not optimized [39]. Here, we design and fabricate SRN ridge waveguides and perform linear and nonlinear characterizations of them, which show a propagation loss of ~2 dB/cm and a moderate nonlinearity of ~4.62 W−<sup>1</sup> m−<sup>1</sup> of the fabricated SRN waveguides, corresponding to a nonlinear RI of ~1.13 <sup>×</sup> 10−<sup>18</sup> m2/W of the SRN material at the wavelength 1.95 μm. Then the SRN nonlinear ridge waveguides are optimized by adjusting the width and the rib height to achieve the maximal nonlinear figure of merit (FOM) and to achieve a broad four-wave mixing (FWM) bandwidth of 518 nm.

#### **2. Structure and Fabrication**

Figure 1a shows the schematic structure of the present SRN ridge waveguide with the total height, rib height and width denoted as *H*, *h*, and *w*, respectively. A thick SRN layer was needed for tightly confining light in the SRN core and obtaining a proper dispersion profile. Here, we used a low-pressure chemical vapor deposition (LPCVD) machine to deposit an 810 nm SRN film on 2.5 μm thermal silicon dioxide (SiO2). The SRN film was first deposited by a reaction between dichlorsilane (SiH2Cl2 at 80 sccm) and ammonia (NH3 at 20 sccm) at 830 ◦C and a pressure of 120 mTorr, and then annealed for three hours at 1150 ◦C. After annealing, the Si excess was measured to be ~18.7% and the SRN film was condensed to 800 nm. A thick resist (CSAR, ~1.1 μm) was then spun and patterned by the electron beam lithography. The patterns were then transferred to the SRN film by using an inductive coupled plasma machine at –10 ◦C with flows of the etching gases SF6/Ar/CF4/CH4 being 4/10/4/2 sccm. SRN ridge waveguides with various widths and rib heights were fabricated with a fixed *H* of 800 nm. After residual resist stripping, the ridge waveguides were finalized with air as the upper cladding and the dimensions were measured, as shown by an exemplified scanning electron microscopy (SEM) image in Figure 1b. RIs at near infrared wavelengths of the deposited and annealed SRN film were measured by an ellipsometry method and fitted with the Sellmeier formula to extract the RIs around 2 μm wavelengths. Here, the RI was calculated to be 2.123 at 1950 nm. Figure 1c shows the simulated mode profile of the fundamental transverse-electric (TE0) mode for a SRN ridge waveguide with

*w* = 1.32 μm and *h* = 640 nm at the wavelength 1950 nm. The power confinement ratio in the SRN core was calculated to be 89.2%, indicating strong confinement of the SRN ridge waveguide.

**Figure 1.** The fabricated silicon-rich nitride (SRN) ridge waveguide. (**a**) Schematic diagram of the structure. (**b**) Scanning electron microscopy image. (**c**) Simulated fundamental transverse-electric (TE0) mode profile at 1950 nm of the waveguide with a width of 1.32 μm and a rib height of 640 nm.

#### **3. Linear and Nonlinear Characterizations**

#### *3.1. Linear Characterizations*

Linear transmission spectra of the fabricated SRN ridge waveguides with lengths of 5.5 mm, 13 mm and 21 mm were measured from 1860 nm to 1980 nm. Then the cut-back method was used to extract the propagation losses and the coupling losses to a tapered lensed fiber for waveguides with different *w* and *h*. In order to reduce the accidental error, six identical waveguides were fabricated and measured for each waveguide dimension configuration. To serve as an example, Figure 2a,b exhibits the propagation and coupling losses with respect to the wavelength for the fabricated SRN ridge waveguides with *h* = 380 nm and *w* = 950 nm, and *h* = 640 nm and *w* = 1320 nm, respectively. Due to the Fabry–Perot interference effect between the two facets of a waveguide, the nonperfect uniformity of the six identical waveguides and the measurement inaccuracy, there were statistical errors for the extracted losses, as shown in the figures. Nevertheless, the propagation losses of the fabricated SRN ridge waveguides were around 2 dB/cm in the investigated wavelength range. Meanwhile, the coupling losses were around 4 dB/facet. Results on the other waveguide dimensions are summarized in Table 1 in the end of this section.

**Figure 2.** The measured propagation loss and coupling loss with respect to the wavelength for the fabricated silicon-rich nitride ridge waveguides with (**a**) *h* = 380 nm and *w* = 950 nm and (**b**) *h* = 640 nm and *w* = 1320.

The propagation loss α of the fabricated SRN ridge waveguide includes linear absorption loss α<sup>1</sup> and scattering loss α2. α<sup>1</sup> can be derived from the imaginary part of the effective RI, *Ne*ff*,r* of the SRN ridge waveguide by considering the absorption of the SRN material and expressed as [41]

$$a\_1 = \frac{2\omega}{\varepsilon} \text{Im}(N\_{eff,r}) \tag{1}$$

where ω is the optical angular frequency and *c* is the speed of light in vacuum. The scattering loss α<sup>2</sup> depends on the sidewall roughness, half the waveguide width *w*/2, RI steps from the SRN core to the surroundings, and scaling factor of a ridge waveguide to a channel waveguide (i.e., *h* = *H*). The sidewall roughness is characterized by the correlation length *L*<sup>c</sup> and the mean squared error σ<sup>2</sup> deviated from the sidewall surface. Thus, the scattering loss α<sup>2</sup> can be expressed as Equation (2) [42]

$$a\_2 = 4.343 \frac{\sigma^2}{\sqrt{2} k\_0 (w/2)^4 \mathcal{N}\_{eff,r}} \text{g} \cdot f \cdot s \tag{2}$$

where the unit of α<sup>2</sup> is dB/m, *k*<sup>0</sup> is the wave vector in vacuum, and function *g* is completely determined by the dimensions of the SRN channel waveguide and *f* is relevant to *Lc* and RI steps. Details for calculating *g* and *f* can be found in Payne and Lacey's work for analyzing a planar waveguide's scattering loss using the exponential autocorrelation function [43]. The scaling factor *s* can be calculated as [42]

$$s = \frac{\delta \mathcal{N}\_{eff,r} / \delta d}{\delta \mathcal{N}\_{eff,r} / \delta d} \tag{3}$$

where *Ne*ff*,c* is the effective RI of the SRN channel waveguide and *d* is half the waveguide width *w*/2.

A larger *h* and a smaller *w* mean more SRN being etched and more light fields overlapping the rough sidewalls and thus a larger scattering loss. In contrast, for the SRN ridge waveguide with a small *h* and a large *w*, almost all the light is confined inside the SRN core and, thus, α is dominantly determined by the absorption loss α1. Here, we reasonably assume α<sup>1</sup> to be 1.7 dB/cm considering it is the minimal propagation loss value we have extracted for the waveguides with different dimensions and from the ridge waveguide with a very shallow etch, i.e., *h* = 230 nm and *w* = 840 nm. Thus, with Equation (2), we can use the measured propagation losses of the SRN ridge waveguides with different dimensions to fit *Lc* and σ2. The fitted values of (5, 11.3) nm for (σ, *Lc*) were found to be able to provide nice fittings of the calculated losses using Equation (2) to the measured propagation losses, as shown in Figure 3, and meanwhile consistent with the values (σ = 5 nm, *Lc* = 45 nm) reported in other literature where the sidewall roughness of a fabricated SRN waveguide was characterized [30]. Thus, it is reasonable to use Equation (2) to predict propagation losses of the SRN ridge waveguide with any *w* and *h*.

**Figure 3.** The measured and calculated propagation losses of the fabricated silicon-rich nitride ridge waveguides with different dimensions.

#### *3.2. Nonlinear Characterizations*

With the experimental setup shown in Figure 4, FWM experiments were implemented to characterize the nonlinear parameter γ of the fabricated SRN ridge waveguides with different dimensions. The pump signal was provided by a commercial continuous-wave (CW) laser (AdValue Photonics AP-CW1) with a fixed wavelength of 1950.1 nm. The probe signal was generated by a homemade CW laser having a thulium doped fiber as the gain media pumped at 793 nm. The wavelength of our homemade CW laser was set at 1953 nm. High power isolators (HP-ISOs) were used to protect the CW lasers and the pump and signal lights were combined by a 9:1 coupler. The combination of a fiber polarization beam splitter (PBS) and a fiber polarization controller (PC) was implemented to guarantee both the pump light and the signal light were TE-polarized when injected into the waveguides by a tapered lensed fiber. The output light was collected by another lensed fiber and the FWM spectra were recorded by an optical spectrum analyzer (OSA, Yokogawa AQ6375B).

**Figure 4.** The experimental setup for four-wave mixing (FWM) experiments in the fabricated silicon-rich nitride (SRN) ridge waveguides. CW: continuous-wave; HP-ISO: high power isolator; PBS: polarization beam splitter; PC: polarization controller; VOA: variable optical attenuator; OSA: optical spectrum analyzer.

The idler signal with ω*<sup>I</sup>* = 2ω*<sup>P</sup>* − ω*<sup>S</sup>* can be generated in the SRN ridge waveguide via the degenerate FWM, where ω*I*, ω*P*, ω*<sup>S</sup>* are the frequencies of the idler, pump and signal lights, respectively. Since the deposited SRN possessing a moderate excess of Si, the TPA and free-carrier absorption (FCA) can be neglected in the 2 μm band. Considering the linear loss, self-phase modulation, cross-phase modulation and FWM, the interaction relationships among these three lights can be described by the following full-vectorial coupling equations [41]

$$\frac{\partial A\_P}{\partial z} = -\frac{1}{2}apAp + j\gamma p \left( |Ap|^2 + 2|A\_S|^2 + 2|A\iota|^2 \right)Ap + 2j\gamma p A\_S A\_I A\_P^\* \exp\left( j\Delta\beta z\right) \tag{4}$$

$$\frac{\partial A\_S}{\partial z} = -\frac{1}{2} a\_S A\_S + j\gamma\_S \left( |A\_S|^2 + 2|A\nu|^2 + 2|A\gamma|^2 \right) A\_S + j\gamma\_S A\_P^2 A\_I^\* \exp\left( -j\Delta\beta z \right) \tag{5}$$

$$\frac{\partial A\_I}{\partial z} = -\frac{1}{2} a\_I A\_S + j\gamma\_I \left( |A\_I|^2 + 2|A\_P|^2 + 2|A\_S|^2 \right) A\_I + j\gamma\_I A\_P^2 A\_S^\* \exp(-j\Delta\beta z) \tag{6}$$

where *Am* (*m* = *P, S, I*) is the field amplitude of the pump, signal, or idler, α*<sup>m</sup>* is the linear loss, Δβ is the linear phase mismatch, and γ*<sup>m</sup>* is the nonlinear parameter, respectively. When the interacting waves are all in the same wavelength band, γ*<sup>m</sup>* can be expressed as Equation (7) [44]

$$\gamma\_{\rm mf} = \frac{\omega\_{\rm mf}}{c A\_{eff}(\omega\_{\rm mf})} \overline{n}\_2(\omega\_{\rm mf}) \tag{7}$$

where *Ae*ff is the effective mode area and *n*<sup>2</sup> is the effective nonlinear RI. They can be calculated as [45]

$$A\_{eff}(\boldsymbol{\omega}\_m) = \frac{\left| \iint \text{Re} \left[ \overrightarrow{\boldsymbol{c}} \left( \mathbf{x}, \boldsymbol{y}, \boldsymbol{\omega}\_m \right) \times \overrightarrow{\boldsymbol{h}}^\* \left( \mathbf{x}, \boldsymbol{y}, \boldsymbol{\omega}\_m \right) \right] \cdot \boldsymbol{\varepsilon}\_z \mathbf{dx} dy \right|^2}{\iint \left\{ \text{Re} \left[ \overrightarrow{\boldsymbol{c}} \left( \mathbf{x}, \boldsymbol{y}, \boldsymbol{\omega}\_m \right) \times \overrightarrow{\boldsymbol{h}}^\* \left( \mathbf{x}, \boldsymbol{y}, \boldsymbol{\omega}\_m \right) \right] \cdot \boldsymbol{\varepsilon}\_z \right\}^2 d\mathbf{x} dy} \tag{8}$$

$$\overline{m}\_{2}(\omega\_{m}) = \frac{\varepsilon\_{0}}{\mu\_{0}} \frac{\iint n\_{2}(\mathbf{x}, \mathbf{y}, \omega\_{m}) n^{2}(\mathbf{x}, \mathbf{y}, \omega\_{m}) \Big| \overline{\hat{\varepsilon}}(\mathbf{x}, \mathbf{y}, \omega\_{m}) \Big|^{2} \Big| \overline{\hat{\varepsilon}}^{\*}(\mathbf{x}, \mathbf{y}, \omega\_{m}) \Big|^{2} d\mathbf{x} dy}{\iint \left\{ \text{Re} \left[ \overline{\hat{\varepsilon}}(\mathbf{x}, \mathbf{y}, \omega\_{m}) \times \overline{\hat{h}}^{\*}(\mathbf{x}, \mathbf{y}, \omega\_{m}) \right] \cdot \hat{\varepsilon}\_{z} \right\}^{2} d\mathbf{x} dy} \tag{9}$$

where *n*(*x*, *y*, ω*m*) and *n*2(*x*, *y*, ω*m*) are the linear and nonlinear RIs of the material at position (*x, y*) at the frequency ω*m*, respectively. <sup>→</sup> *<sup>e</sup>* (*x*, *<sup>y</sup>*, <sup>ω</sup>*m*) and <sup>→</sup> *h* (*x*, *y*, ω*m*) are the electric and magnetic field distributions on the waveguide transverse plane. *e*ˆ*<sup>z</sup>* is the unit vector along the propagation direction. ε<sup>0</sup> and μ<sup>0</sup> are the dielectric constant and the permeability constant, respectively. Finally, the conversion efficiency η (in the unit of dB) is defined as the ratio of the output idler power to the output signal power, that is

$$\eta(\text{dB}) = -10 \lg \frac{\left| A\_I(L) \right|^2}{\left| A\_S(L) \right|^2} \tag{10}$$

where *L* is the physical length of the waveguide.

We have measured the output FWM spectra of 21-mm-long SRN ridge waveguides with different waveguide dimensions under various coupled pump powers and extracted the conversion efficiencies (CEs), when the signal wavelength and the incident signal power were fixed at 1953 nm and 10 mW, respectively. For example, Figure 5a,b shows the measured FWM spectra of the waveguide with *h* = 380 nm and *w* = 950 nm under a coupled pump power of 61.5 mW, and *h* = 640 nm and *w* = 1320 under a coupled pump power of 52.4 mW, respectively. CEs of -53.2 dB and -51.1 dB can be extracted from the spectra. It should be noted that there was already some idler light generated in the incident fiber. We have normalized the output idler powers to that coming from the fiber and found little difference on the CEs between the normalizations before and after. Figure 5c,d shows the measured and normalized CEs with respect to the coupled pump power for the waveguides with *h* = 380 nm and *w* = 950 nm, and *h* = 640 nm and *w* = 1320 nm, respectively. By solving the equations from (4) to (10), we can fit the measured CEs versus the coupled pump power and the nonlinear parameter γ can be extracted to be 2.79 W−<sup>1</sup> m−<sup>1</sup> for the waveguide with *h* = 380 nm and *w* = 950 nm and 4.62 W−<sup>1</sup> m−<sup>1</sup> for the waveguide with *h* = 640 nm and *w* = 1320, respectively. With the extracted γ value, the nonlinear index *n*<sup>2</sup> of the SRN material can also be calculated by using Equation (7). Here, we assume the whole nonlinear effect was contributed by the SRN material since the silica surroundings have a nonlinear RI orders lower than that of the SRN material [45].

We have summarized the linear and nonlinear properties of the fabricated SRN ridge waveguides with different waveguide dimensions in Table 1. While there were variations for these measured or extracted values, the fabricated waveguides generally exhibited a linear loss of ~ 2 dB/cm and a nonlinear index of ~1.13 <sup>×</sup> 10−<sup>18</sup> m2/W. These values are indeed consistent with that from literatures where the SRN material and waveguide were measured at 1550 nm [30].

**Figure 5.** Measured output four-wave mixing spectra (**a,b**) and the wavelength conversion efficiencies with respect to the coupled pump powers (**c,d**) of the fabricated silicon-rich nitride ridge waveguides with a length of 21 mm. (**a**) and (**c**) are for the waveguide with *h* =380 nm and *w* = 950 nm. (**b**) and (**d**) are for the waveguide with *h* = 640 nm and *w* = 1320 nm.

**Table 1.** Summary of the linear and nonlinear properties of the fabricated silicon-rich nitride ridge waveguides with various rib heights (*h*) and widths (*w*) at 1950 nm.


#### **4. Optimal Design for the Maximal Nonlinear Energy E**ffi**ciency**

For waveguides without nonlinear losses, we can use a figure of merit (FOM) to evaluate the nonlinear energy efficiency of the waveguide, this is the ratio of the waveguide nonlinear parameter to the linear loss, expressed as [46]

$$\text{FOM} = \frac{\mathcal{V}}{a} \tag{11}$$

For the proposed SRN ridge waveguide, both γ and α are dependent on the width and rib height. While a deeper etching and a moderately smaller width can give a stronger confinement and hence a larger nonlinearity, this also leaves more light overlapping the rough sidewalls and therefore gives a larger linear propagation loss. Thus, the FOM can be used to fairly evaluate the overall nonlinear efficiency of the SRN ridge waveguide with different dimensions. We have calculated the nonlinear

parameter, linear loss and FOM of the SRN ridge waveguides for many combinations of *w* and *h* by knowing the linear and nonlinear RIs of SRN and the fabrication quality (i.e., *L*<sup>c</sup> and σ) and shown them in Figure 6a–c, respectively. As expectations, γ and α generally have the similar varying trend as the waveguide width or etch depth changes. This finally yielded a maximal value of 0.0804 W−<sup>1</sup> of the FOM at the rib height *h* = 700 nm and width *w* = 1100 nm. Here, γ an α were calculated to be 5.50 W−<sup>1</sup> m−<sup>1</sup> and 2.97 dB/cm, respectively.

**Figure 6.** The calculated waveguide nonlinear parameter γ (**a**), linear loss α (**b**) and figure of merit (FOM) (**c**) of a silicon-rich nitride ridge waveguide with respect to the width *w* and rib height *h*. Here, the total height *H* is fixed at 800 nm.

We have also calculated the nonlinear performances of the FWM wavelength conversion for the designed SRN ridge waveguide with the maximal FOM (*h* = 700 nm, *w* = 1100 nm). Figure 7a shows the CEs as the waveguide length increases, when the input pump wavelength/power and input signal wavelength/power were set to be 1950 nm/500 mW and 1951 nm/1 mW, respectively. Here, CE is defined as the ratio of the output idler power to the input signal power. The nonlinear effect for generating new idler photons dominated the power of the idler for a short waveguide while for a long waveguide, the linear loss took charge. Thus, there was an optimal length for maximal CE and, here in Figure 7a, a length of 14.6 mm was obtained for a maximal CE of –36.2 dB. It is worth noting that, to obtain a higher CE, one can further increase the pump power or use a resonating structure to enhance the nonlinear interactions between light and the SRN ridge waveguide. Figure 7b shows the calculated wavelength dependence of the CE for the optimized waveguide with a length of 14.6 mm at a fixed pump wavelength of 1950 nm and an adjusted signal wavelength. The 3 dB conversion bandwidth was found to be 94 nm.

**Figure 7.** The calculated four-wave mixing conversion efficiency for the designed silicon-rich nitride ridge waveguide (*h* = 700 nm, *w* = 1100 nm) with a maximal figure of merit (FOM). (**a**) Dependence on the waveguide length. Here, difference of the pump and signal wavelengths is 1 nm. (**b**) Dependence on the signal wavelength. Here, the waveguide length is 14.6 mm and the pump wavelength is 1950 nm.

#### **5. Optimal Design for a Superbroad FWM Conversion Bandwidth**

Although the optimally designed waveguide above exhibits a maximal nonlinear efficiency, the bandwidth is still very limited due to a dispersion which is not small enough (164.3 ps/nm/km) at 1950 nm as shown by the calculated wavelength dependence of the dispersion in Figure 8a (green solid line). Figure 8a also shows the calculated wavelength-dependent dispersion curves for the SRN ridge waveguides with some other dimensions. The dispersion was found to be close to 0 (−6.0 ps/nm/km) at 1950 nm for the waveguide with *h* = 625 nm and *w* = 1050 nm. Such a low dispersion is expected to give a broad nonlinear bandwidth. For such an SRN ridge waveguide, the waveguide nonlinear parameter and the linear loss were calculated to be 5.24 W−1m−<sup>1</sup> and 2.88 dB/cm, respectively, thus producing a FOM of 0.0790 W<sup>−</sup>1. This FOM is only a little compromised compared with the maximal one (0.0804 W<sup>−</sup>1). Furthermore, the calculated FWM CE was calculated to exhibit a maximal value of –36.22 dB for a waveguide length of 15 mm, when the input pump wavelength/power and input signal wavelength/power were set to be 1950 nm/500 mW and 1951 nm/1 mW, respectively. Figure 8b shows the calculated wavelength dependence of the FWM CE for the designed SRN ridge waveguide and the 3 dB conversion bandwidth was found to be as broad as 518 nm thanks to the small waveguide dispersion around 2 μm. Besides, one can also find that it may be not a good idea to etch the SRN layer through, i.e., *h* = 800 nm, see the solid brown curve in Figure 8a, for a low dispersion in the 2 μm band, which is nevertheless mostly implemented in the 1550 nm band.

**Figure 8.** (**a**) Calculated wavelength dependence of the dispersion for the silicon-rich nitride ridge waveguides with various waveguide dimensions. (**b**) Calculated four-wave mixing conversion efficiency with respect to the signal wavelength for a SRN ridge waveguide (*h* = 625 nm, *w* = 1050 nm) with the zero-dispersion wavelength at 1950 nm.

#### **6. Discussion and Conclusions**

We have proposed and fabricated the silicon-rich nitride ridge waveguides and characterized their linear and nonlinear performances at 2 μm wavelengths. SRN ridge waveguides with different rib heights and widths were fabricated, exhibiting a linear loss around 2 dB/cm. Four-wave mixing experiments with a CW pump at 1950 nm were performed in the fabricated SRN ridge waveguides and revealed a waveguide nonlinear parameter of ~3-6 W<sup>−</sup>1m−<sup>1</sup> around 2 μm. With the measured and extracted parameters which characterize the fabrication quality and the material property, optimal design of the SRN ridge waveguides was carried out and a maximal FOM of 0.0804 W−<sup>1</sup> was found among ridge waveguides with various dimension configurations at 1950 nm. Meanwhile, the ridge waveguide could also be designed to achieve the FWM conversion with a superbroad bandwidth (518 nm) and little compromise of the nonlinear FOM. First, these results show that a stripe silicon nitride waveguide may not be the best choice under some fabrication quality and material properties when moving the operating wavelengths from the conventional telecommunication band to longer wavelengths like the 2 μm band. Second, although the nonlinear conversion efficiency (–51.1 dB) achieved in our fabricated SRN ridge waveguide was much smaller than that of a silicon waveguide (–10 dB [24]) and an SRN waveguide (–18 dB [39]) pumped with pulses, it is indeed the first time, to the best of our knowledge, the nonlinear properties of a silicon nitride waveguide in the 2 μm spectral window have been revealed using a CW pump. Last, the proposed design method for optimal nonlinear performances, in terms of either FOM or the bandwidth, is expected to open new avenues for achieving better on-chip nonlinear waveguides for applications involving longer wavelengths like the 2 μm band.

**Author Contributions:** Conceptualization, design and fabrication, Z.T. and X.G.; Measurement, Z.T., and D.C.; Data analysis, Z.T., S.G. and X.G.; Writing—original draft preparation, Z.T.; writing—review and editing, D.C., H.H., S.G. and X.G.; All authors have read and agreed to the published version of the manuscript.

**Funding:** National Natural Science Foundation of China (Grant No. 61875172 and 61475138), Zhejiang Provincial Natural Science Foundation of China (Grant No. LD19F050001 and LY19F050014), Det Frie Forskningsråd (Danish Council for Independent Research) (DFF-7107-00242) and Villum Fonden (Villum Foundation) (023112, 00023316 and 15401).

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

#### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Review* **Optical Realization of Wave-Based Analog Computing with Metamaterials**

**Kaiyang Cheng 1, Yuancheng Fan 2,\*, Weixuan Zhang 3, Yubin Gong 1,4, Shen Fei 1,\* and Hongqiang Li 5,\***


**Abstract:** 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 counterpart, 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 briefly reviewed. We first 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-efficiency signal processing by exploiting quantum behaviors. The challenges and future opportunities in the booming research field are discussed.

**Keywords:** analog optical computing; metamaterials; metasurfaces; quantum algorithm; edge detection

#### **1. 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–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–8].

**Citation:** Cheng, K.; Fan, Y.; Zhang, W.; Gong, Y.; Fei, S.; Li, H. Optical Realization of Wave-Based Analog Computing with Metamaterials. *Appl. Sci.* **2021**, *11*, 141. https://dx.doi.org/ 10.3390/app11010141

Received: 26 November 2020 Accepted: 23 December 2020 Published: 25 December 2020

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license ( https://creativecommons.org/ licenses/by/4.0/).

Analog optical computing offers unique advantages of real-time, power-efficient, highthroughput imaging processing abilities originating from the wave-based nature [7,9–11]. Compared with standard digital processes, all-optical analog approaches do not refer to photoelectric converters but directly manipulate optical signals both temporally and spatially [12]. In the temporal domain, high-speed pulse waveform modulation enables numerous applications, including analog computing [13–19], differential equations solving [19–23], optical memory [24,25], photonic neural networks [26,27] and complex nonlinear system simulation [28–30]. In the spatial domain, the input functions are indicated by the spatial wavefronts which will be mathematically transformed with pre-designed optical systems. Therefore, this computing platform has intrinsic parallel characteristic, showing great potential for accelerating the processing of megalo-capacity datasets and images [31–33].

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 [37–68].

However, traditional analog optical computing requires bulky optical components, such as diffractive lenses and spatial filters, which are inconvenient for miniaturization and integration of modern ultracompact optics. Empowered by recent development of nanofabrication technologies, metamaterial [69–79] or its two-dimensional counterpart, metasurface [80–99], is able to tailor subwavelength building blocks on the scales of micro- or nanometers, providing unprecedented flexibility to arbitrarily control the electromagnetic waves that are unattainable in the nature. Owning to their powerful wave manipulation abilities and deeply subwavelength characteristics, those meta-structures can significantly reduce the complexity and shrink the size of computing systems, making it possible to realize chip-level all-optical information processing systems [100,101].

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.

#### **2. 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.

**Figure 1.** (**a**) The conceptual sketch of computational metamaterials; schematic of the general protocol for (**b**) MS approach; and (**c**) GF approach; adapted with permission from [102], AAAS, 2014.

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 4*f* 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.

#### *2.1. 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

$$\mathcal{S}(\mathbf{x}, \mathbf{y}) = \iint f(\mathbf{x}', \mathbf{y}') h(\mathbf{x} - \mathbf{x}', \mathbf{y} - \mathbf{y}') d\mathbf{x}' d\mathbf{y}' \tag{1}$$

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:

$$\log(\mathbf{x}, \mathbf{y}) = \mathcal{F}^{-1}\left\{ H(k\_x, k\_y) \mathcal{F}[f(\mathbf{x}, \mathbf{y})] \right\} \tag{2}$$

where (<sup>F</sup> <sup>−</sup>1{·})F{·} represents the (inverse) Fourier transform, abbreviated (*I*)*FT*. *kx* and *ky* are frequency variables in the spatial Fourier space, *H*(*kx*, *ky*) 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 *Ein*(*x*, *y*) (*Eout*(*x*, *y*)), in this case, the spatial variables (*x*, *y*) can play the role of (*kx*, *ky*). In wave-based computing system, Equation (2) can be rewritten as:

$$E\_{out}(\mathbf{x}, y) = \mathcal{F}^{-1}\left\{H(\mathbf{x}, y)\mathcal{F}[E\_{in}(\mathbf{x}, y)]\right\} \tag{3}$$

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–110]. It worth nothing that, although we can define a GRIN(−) subblock to perform *IFT* in which the effective constitutive parameters should satisfy *<sup>ε</sup>*(*μ*)*GRIN*(−)(*x*, *<sup>y</sup>*) = <sup>−</sup>*ε*(*μ*)*GRIN*(+)(*x*, *<sup>y</sup>*), 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:

$$\mathcal{F}\{\mathcal{F}\{E\_{out}(\mathbf{x},\mathbf{y})\}\} \propto E\_{out}(-\mathbf{x},-\mathbf{y})\tag{4}$$

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 *n*th derivation operators. Based on derivative property of the Fourier transform, *<sup>d</sup><sup>n</sup> <sup>f</sup>*(*y*)/*dy<sup>n</sup>* = F <sup>−</sup><sup>1</sup> (*iky*) *<sup>n</sup>*F[ *<sup>f</sup>*(*y*)] , according to Equation (3), the expression *iky <sup>n</sup>* is the required transfer function, where *H*(*y*) ∝ *iky <sup>n</sup>* or *H*(*y*) ∝ (*iy*) *<sup>n</sup>*. It means that the amplitude and phase of impinging optical wave will be modulated by *eik*<sup>Δ</sup> ∝ (*iy*) *<sup>n</sup>* 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 *<sup>λ</sup>*0, *<sup>ε</sup>*<sup>Δ</sup> and *<sup>μ</sup>*<sup>Δ</sup> 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*) <sup>−</sup><sup>1</sup> 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.

#### 2.1.1. Reflective Configurations of MS Approach

Plasmonic metasurfaces [83,111–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 amplitude and phase of reflected cross-polarized field by varying the lengths of nanobricks (see Figure 2b).

**Figure 2.** (**a**) The unit cell of plasmonic metasurfaces consist of gold nanobrick atop a dielectric spacer and thick metal ground; adapted with permission from [120], ACS Publications, 2015. (**b**) Schematic representation of dielectric meta-reflect-array for the first-order derivation; adapted with permission from [122] © The Optical Society. (**c**) Top: eleven silver dendritic structure units with different parameters; left: the unit cell of plasmonic metasurfaces with dendritic structure instead of nanobrick; right: comparison of the theoretical and simulated results of the position-dependent reflection coefficient. Adapted with permission from [121] © The Optical Society.

#### 2.1.2. 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).

**Figure 3.** (**a**) The sketch of graphene-based metaline; adapted with permission from [32] © The Optical Society. (**b**) Dielectric metamaterial with drilling hole arrays to solving second-order differential equation; adapted with permission from [123], IOP Publishing, 2016. (**c**) Schematic of on-chip 1D high-contrast transmitarray, the rectangular slot arrays are etched on the silicon-on-insulator (SOI) substrate with the period of 500 nm. (**d**) Comparison of the input (top) and output (bottom) intensity profile of an on-chip differentiator; adapted with permission from [100], Springer Nature, 2019. (**e**) Left: schematic of the concept for spin-dependent function control; middle: schematic of the dielectric metasurface spatial filter; right: photograph of the fabricated metasurface. Scale bar: 500 μm. Insets: scanning electron micrographs showing the top and oblique view of TiO2 nanopillar array. Scale bar: 1 μm. Adapted with permission from [125], ACS Publications, 2020.

#### *2.2. 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*(*ky*) of the multi-layered slab to fit the desired GF kernel *G*(*ky*) within 0 < *ky* < *k*0, in which *ky* is the transverse wavevector, and *k*<sup>0</sup> is the free space wavevector. A fast synthesis algorithm is adopted to minimize the difference between *T*(*ky*) and *G*(*ky*), which can be expressed as:

$$errr = \sum\_{i=1}^{M} \left( w\_{\tau} \left( \text{Re} \left[ \tilde{G}(k\_{y,i}) - \tilde{T}(k\_{y,i}) \right] \right)^{2} + w\_{i} \left( \text{Im} \left[ \tilde{G}(k\_{y,i}) - \tilde{T}(k\_{y,i}) \right] \right)^{2} \right) \tag{5}$$

where *M* is the number of layers, *wr* and *wi* 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, *wr* = 2, *wi* = 1, and relative permeability of each layer keep at unity.

In recent years, extensive studies with resonant or non-resonant structures were proposed with the GF approach, such as multi-layered slabs [126,127], photonic crystal [128], diffraction gratings [15,31,129–134], surface plasmon-based devices [130,135–137], dielectric metasurfaces [138–143], nonlocal metasurfaces [144,145], random medium [146] and inverse-designed meta-structures [147]. Moreover, various physical mechanisms including Brewster effect [148], Goos-Hänchen effect [149] and spin hall effect [150] can also be applied to perform optical analog computing. Therefore, the GF approach can be generalized to a fundamental protocol for designing the particular transfer function (defined as *<sup>H</sup>*(*kx*, *ky*) = *<sup>E</sup>out <sup>y</sup>* (*kx*, *ky*)/*Ein <sup>y</sup>* (*kx*, *ky*)) which should be approximately proportional to the desired mathematical operations.

#### 2.2.1. 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 TiO2 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 *Pinc*(*x*, *y*) is formed by zeroth order diffraction to the transmitted profile of *ptr*(*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 Si3N4 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.

**Figure 4.** (**a**) Schematic and (**b**) photograph of guided-mode resonant diffraction grating for firstorder differentiation; adapted with permission from [131] © The Optical Society. (**c**) The schematic and (**d**) spatial spectral transfer function of the optical spatial differentiator based on subwavelength high-contrast gratings. Adapted with permission from [132], AIP Publishing, 2018. (**e**) Schematic of the first-order differentiator based on subwavelength metallic gratings. (**f**) The experimental results for the incident and transmitted images (left) and the normalized field profiles at the position of white dashed line. Adapted with permission from [133] © The Optical Society.

#### 2.2.2. 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 *Sin*(*kx*)) with TMpolarized are used to excite the surface plasmon between glass substrate and silver film, *kx* is the transverse component of the wavevector (see Figure 5a). In this case, the amplitude of reflected wave (transverse profile of *Sout(kx*)) is the result of the interference between the direct reflection and the leakage of SPPs. The spatial spectral transfer function around *kx* = 0 can be expressed as:

$$H(k\_x) = \varepsilon^{i\varphi} \frac{ik\_x + A}{ik\_x + B} \tag{6}$$

where *ϕ* is the phase change related to the process of the direct reflection at the glass–metal interface. *A* = (*αspp* − *α*1)/ cos *θ*<sup>0</sup> and *B* = (*αspp* + *α*1)/ cos *θ*0, where *α*<sup>1</sup> and *αspp* are the radiative leakage rate of the SPP and the intrinsic material loss rate, respectively. When the critical coupling condition is satisfied (*α*<sup>1</sup> = *αspp*), Equation (6) can be simplified and approximated as *H*(*kx*) ∝ (*eiϕ*/*B*)*ikx*, which corresponds to the first-order differentiation.

**Figure 5.** (**a**) Top: schematic of the surface-plasmon-based scheme to perform spatial differentiation with the Kretschmann configuration. Bottom: photograph of the three samples with different thicknesses; adapted with permission from [135], Springer Nature, 2017. (**b**) Top: schematic of the backscattering-immune first-order spatial differentiator based on nonreciprocal plasmonic platform. Bottom: the simulated results of real and imaginary parts of the eigenvector with different external magnetic fields. Adapted with permission from [137], APS, 2019.

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 *ω*<sup>0</sup> in the spatiotemporal coordinate can be expressed as:

$$s\_{\rm out}(\mathbf{x}, t) = c\_s \frac{d\mathbf{s}\_{\rm in}(\mathbf{x}, t)}{d\mathbf{x}} + c\_t \frac{d\mathbf{s}\_{\rm in}(\mathbf{x}, t)}{dt} \tag{7}$$

where *cs* and *ct* are two constants where *cs* = *ei<sup>φ</sup>* cos *<sup>θ</sup>*<sup>0</sup> <sup>2</sup>*a*<sup>1</sup> and *ct* <sup>=</sup> *<sup>e</sup>iφ*(2*a*1) −1 <sup>1</sup> *vg* <sup>−</sup> sin *<sup>θ</sup>*<sup>0</sup> *vgls* , *φ* is phase change of direct reflection without excitation of SPP, *θ*0, *a*1, *vg*, *vgls* 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 1013 frame/s which is restricted by the time of establishing the SPP leaky radiation.

#### 2.2.3. 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 *TTE*(*ϕ*, *kr*) ≈ *<sup>α</sup>*(*ϕ*)*k*<sup>2</sup> *r*, where *ϕ* is the azimuthal angle, *α*(*ϕ*) is the gain of the system, and *kr* = *k*2 *<sup>x</sup>* + *k*<sup>2</sup> *<sup>y</sup>*. 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 *s*- or *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).

**Figure 6.** (**a**) Schematic of the square patch-arrays metasurface for 2D edge detection; (**b**) 2D spatial transfer function spectrum; (**c**) edge detection of square-shaped input beams; adapted with permission from [143], IEEE, 2018. (**d**) Schematic of the basic configuration of dielectric metasurface for 2D spatial differentiation; (**e**) the transfer functions for or both linearly polarized light fields; (**f**) the output image of 2D edge detection for *x* polarization light waves. Adapted with permission from [141] © The Optical Society.

#### **3. Other Emerging Approaches**

#### *3.1. 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).

**Figure 7.** (**a**) Schematic of sinusoidally modulated of split-ring resonators; (**b**) transmission spectrum with different incident angle; (**c**) schematic of asymmetric metasurface with a horizontally misplaced array of metallic wires for the first-derivative operation; (**d**) rotated metasurface for 2D secondderivative operation. (**e**) Input image and (**f**) output image of 2D edge detection for linear polarized wave. Adapted with permission from [144], APS, 2018.

#### *3.2. 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.

**Figure 8.** (**a**) Schematic the random medium scheme for wave-based analog computation; (**b**) tamed contribution of segment *D* to observation points; (**c**) experimental setup in a metallic cavity of irregular geometry. Adapted with permission from [146], APS, 2018.

#### *3.3. 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: *<sup>g</sup>*(*u*) = *Iin*(*u*) + *<sup>b</sup> <sup>a</sup> K*(*u*, *v*)*g*(*v*)*dv*, where *g*(*u*) is the function to be solved, *K*(*u*,*v*) and *Iin*(*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).

**Figure 9.** (**a**) The conceptual sketch of the inverse-designed metamaterial platform in solve integral equation; (**b**) time snapshot of the simulation results for the transmissive configuration with external feedback network; (**c**) the reflective configuration of the internal feedback metamaterial system. (**d**) Left: the photograph of the constructed metamaterial system. Right: time snapshot of the simulation results for the reflective configuration with internal feedback network. Adapted with permission from [147], AAAS, 2019.

#### *3.4. 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 *ky* = 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.

**Figure 10.** (**a**) The schematic of the rotated structure for realizing an even or odd Green's function; (**b**) the exact GF distribution and its approximation around *ky* = 0 based on the Taylor series. Adapted with permission from [148] © The Optical Society.

#### *3.5. 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 *kx* = 0; (**c**) output image of edge detection with different rotating angles. Adapted with permission from [149], AIP Publishing, 2020.

#### *3.6. 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: <sup>|</sup>*ϕout* <sup>=</sup> *<sup>i</sup>* 2 *dy*[*ϕin*(*<sup>y</sup>* <sup>+</sup> *<sup>δ</sup>*) <sup>−</sup> *<sup>ϕ</sup>in*(*<sup>y</sup>* <sup>−</sup> *<sup>δ</sup>*)]|*<sup>y</sup>*, 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 *ky* = 0 is depicted in Figure 12b. the experimental results for 1D edge detection are shown in Figure 12c.

**Figure 12.** (**a**) Schematic of spatial differentiation from the SHE of light on an optical planar interface between two isotropic materials; (**b**) the comparison of theoretical and experimental results for the spatial spectral transfer function for *kx* = 0; (**c**) the results of edge detection with different target images stored in *Einx* and *Einy*, respectively. Adapted with permission from [150], APS, 2019.

#### *3.7. 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* √*<sup>N</sup>* 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–53,55–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 *Um*, 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).

**Figure 13.** (**a**) Schematic of the general protocol of simulating quantum search algorithm with metamaterials. (**b**) Experimental results of the time traces (top) and the field intensity (bottom) of the output with different iterations. Adapted with permission from [154], Wiley, 2017. (**c**) Simulation results of metamaterial-based DJ algorithm analogy with constant, simple balance and complex balance functions, respectively. (**d**) The photograph (top) and experimental results of three samples for simulating DJ algorithm. Adapted with permission from [155] © The Optical Society.



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.

#### **4. 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.

**Author Contributions:** Y.F. conceived the review; K.C., Y.F. wrote the manuscript; W.Z., Y.G., S.F., H.L. revised the manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Natural Science Foundation of China (No. 11774057, 12074314, 11674266), This research was funded by Funds for Equipment Advance Research during the 13th Five-Year Plan under Grant 61400030305.

**Data Availability Statement:** Data sharing not applicable.

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

#### **References**


MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Applied Sciences* Editorial Office E-mail: applsci@mdpi.com www.mdpi.com/journal/applsci

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18

www.mdpi.com

ISBN 978-3-0365-3039-0