## Computational Metamaterials

Optical signal processing of light waves can represent certain mathematical functions and perform computational tasks on signals or images in an analog fashion. However, the complex systems of lenses and filters required are bulky. Metamaterials can perform similar optical processing operations but with materials that need only be a wavelength thick. **Silva et al.** (p. 160; see the Perspective by

**Sihvola**) present a simulation study that shows how an architecture based on such metamaterials can be designed to perform a suite of mathematical functions to create ultrathin optical signal and data processors.

## Abstract

We introduce the concept of metamaterial analog computing, based on suitably designed metamaterial blocks that can perform mathematical operations (such as spatial differentiation, integration, or convolution) on the profile of an impinging wave as it propagates through these blocks. Two approaches are presented to achieve such functionality: (i) subwavelength structured metascreens combined with graded-index waveguides and (ii) multilayered slabs designed to achieve a desired spatial Green’s function. Both techniques offer the possibility of miniaturized, potentially integrable, wave-based computing systems that are thinner than conventional lens-based optical signal and data processors by several orders of magnitude.

Metamaterials—artificially engineered structures designed to exhibit unconventional wave properties—have recently opened exciting venues for unprecedented control and manipulation of electromagnetic fields, as well as other forms of waves, such as acoustic, elastic, mechanical, and matter waves (*1*–*5*). The spatial and/or temporal engineering of their bulk parameters, such as the permittivity and permeability of an artificial material, enables shaping of local and global field distributions to achieve a desired functionality. Transformation optics (*6*, *7*), for instance, has shown that by locally tailoring the electric and magnetic response of a material, one may be able to achieve functionalities of interest, such as cloaking (*6*–*8*), field concentration or rerouting (*9*–*14*), and optical illusions (*15*).

Here, we introduce the concept of “computational metamaterials” that can be used to perform mathematical operations—such as spatial differentiation, integration, or convolution—as electromagnetic waves propagate through them. Traditional analog computers, or “calculating machines” in the form of mechanical, electronic, and hybrid analog computers, have been designed in the past to perform such mathematical computations (*16*–*18*), but these solutions have suffered from substantial limitations, including relatively large size and slow response. Mechanical deformation of arrays of bars and pivots to achieve a certain functional response (*19*) and construction of linkages to synthesize a portion of an algebraic curve (*20*) have also been explored. Our idea for metamaterial-based computing is conceptually sketched in Fig. 1, in which metamaterial blocks produce a wave profile at their output proportional to the desired operation performed on an arbitrary input wave profile by suitably manipulating the impinging wave as it propagates through them. Although a selected set of these functionalities may be achieved with well-known optical signal processing or Fourier optics using conventional lenses (*21*–*23*), our approaches provide the possibility of highly compact, potentially integrable architectures within much smaller volumes and, in some cases, over subwavelength thicknesses, ensuring controlled manipulation and processing of the incoming signal.

We introduce two approaches to computational metamaterials: (i) a metasurface (MS) approach, in which thin planar metamaterial blocks perform mathematical operations in the spatial Fourier domain and (ii) a Green’s function (GF) approach, in which multilayered metamaterial slabs realize a desired GF (i.e., spatial impulse response) at the output, allowing the synthesis of mathematical operations of interest. For simplicity, we focus on metamaterials that realize linear, shift-invariant operations on transverse spatial variables, assuming exp(–*i*ω*t*) monochromatic signals. Without loss of generality, we assume two-dimensional scenarios in which wave signals and the metamaterial parameters depend only on the longitudinal direction of propagation *x* and one transverse coordinate *y*.

We base our first (MS) approach on the fact that the linear convolution *g*(*y*) = ∫ *f*(*u*)*G*(*y* – *u*)*du* between an arbitrary field distribution *f*(*y*) and the GF *G*(*y*) associated with the operator of choice may be expressed in the spatial Fourier space as , where the tilde indicates Fourier transform. In this approach, our metamaterial system consists of three cascaded subblocks: (i) a Fourier transform subblock, (ii) a suitably tailored MS spatial filter that applies the operation in the Fourier domain, and (iii) an inverse Fourier transform subblock (Fig. 2A). For the first subblock, we consider a two-dimensional graded-index (GRIN) dielectric slab with permeability μ = μ_{0} and a parabolic variation of permittivity
(1)where ε_{c} is the permittivity at the central plane of the GRIN and *L*_{g} is the characteristic length (*24*). In Fig. 2A, GRIN(+) indicates a conventional GRIN, with (+) referring to permittivity and permeability being positive. In the paraxial approximation, the GRIN(+) slab operates as a Fourier transformer at the “focal” length *L*_{g} along the direction of propagation (*24*). For the inverse Fourier transform [subblock GRIN(–)], we use an ideal GRIN structure with negative parameters: permeability μ = –μ_{0} and permittivity ε = –ε(*y*). Following the concept of complementary materials (*25*), the GRIN(–) has the inverse functionality of the GRIN(+); that is, it acts as an inverse Fourier transformer.

The metastructure screen is designed using two different approaches (insets in Fig. 2A): (i) a “conceptual” thin MS with volumetric parameters ε_{ms}(*y*)/ε_{0} = μ_{ms}(*y*)/μ_{0} (left inset) and (ii) a realistic three-layered thin metamaterial based on an appropriate alternation of two materials—in this case, aluminum-doped zinc oxide (AZO) (*26*) and silicon (Si)—with different levels of loss (right inset). Our “meta-transmit-array” (MTA) design is inspired by (*27*) and is extended to the present case (*28*) so as to control both phase and magnitude of the transmitted wave according to the desired local transfer function while minimizing reflections. We apply these concepts to realize some fundamental mathematical operators.

First-order derivative: *g*(*y*) ∝ *df*(*y*)/*dy*. In this case, and we first consider a transversely inhomogeneous MS with thickness Δ < λ_{0} (where λ_{0} is the free-space wavelength), in which the transverse coordinate *y* plays the role of *k _{y}*. The relative permittivity and permeability of such a MS (laterally limited to |

*y*| ≤

*W/*2), required to achieve the normalized transfer function proportional to

*iy/*(

*W/*2), are (

*28*) (2)Once this subblock is sandwiched between a GRIN(+) and a GRIN(–), each one with length

*L*

_{g}, we obtain a first-differentiator metamaterial system. Figure 2B shows our numerical simulation for the input function

*f*(

*y*) =

*ay*exp(–

*y*

^{2}/

*b*) with

*a*= 0.7 μm

^{–1}= and

*b*= 10 μm

^{2}= , used to test all the operators in Fig. 2. In this set of simulations, we assume λ

_{0}= 3 μm,

*W*= 29.75 μm ≈ 10λ

_{0},

*L*

_{g}= 35 μm ≈ 12λ

_{0}, ε

_{c}= 2.01ε

_{0}, and Δ = 1 μm = λ

_{0}/3. Although this design principle can be applied to any wavelength, here we have chosen to operate at infrared frequencies, which are of interest for various applications; this allows us to use realistic CMOS (complementary metal-oxide semiconductor)–compatible materials. The input function, depicted in perspective at the left of Fig. 2B, is the

*z*-component of the electric field distribution at the entrance of this metamaterial system. We also show the distribution (snapshot in time) of this field component as the wave propagates through the entire system. At the output plane (2

*L*

_{g}+ Δ), we obtain the electric field distribution

*g*(

*y*) shown also in Fig. 2C, which is indeed closely proportional to the first spatial derivative of the input function (green line).

To get closer to a practical realization, we also consider the more realistic MTA design and replace the ideal GRIN(–) subblock with a GRIN(+), exploiting the relation F{F[*g*(*y*)]} ∝ *g*(–*y*), with F(·) denoting Fourier transform, thereby expecting the output function to be proportional to a “mirror image” of the desired result. Simulation results for this scenario (Fig. 2D) are clearly proportional to the analytical result.

Second-order derivative: *g*(*y*) ∝ *d*^{2}*f*(*y*)/*dy*^{2}. A metamaterial that performs second spatial differentiation could be designed as a cascade of two first-differentiator blocks. However, to demonstrate the generality of our concept, we design a MS that directly performs the second differentiation. Thanks to our “modular” scheme, we simply replace the MS subblock by another one with the desired transfer function [*iy/*(*W*/2)]^{2}, which requires volumetric permittivity and permeability functions
(3)for |*y*| ≤ *W/*2. Figure 2E presents our numerical simulation for this case, confirming excellent functionality.

Integration: *g*(*y*) ∝ ∫ *f*(*y*)*dy*. For this operation, the MS subblock needs to implement the transfer function (*iy/d*)^{–1}, where *d* = 1 μm is an arbitrary normalizing length. To avoid gain requirements for transmission coefficients with magnitude larger than unity for |*y*| < *d*, we truncate the required transfer function at |*y*| = *d* and assume its magnitude to be unity within this range. Therefore, the MS has volumetric constitutive parameters
(4)for |*y*| > *d*, and
(5)for |*y*| < *d*. Figure 2F presents the simulation results for the integral operation performed by this structure compared with the expected result, again demonstrating good functionality.

Convolution: *g*(*y*) = ∫ *f*(*u*)*G*(*y* – *u*)*du*. We also design a MS that performs convolution with a rectangular kernel of width *W*_{k} = 16 μm = 5.33λ_{0}. As before, our MS realizes a transfer function corresponding to the Fourier transform of the kernel
(6)where
(7)is a GRIN scale factor defined in (*24*). For such a transfer function, the required constitutive parameters are
(8)The results of our simulation are shown in Fig. 2G, demonstrating good agreement with the analytical results.

In our second (GF) approach, we design a multilayered metamaterial slab, transversely homogeneous but longitudinally inhomogeneous, that realizes an output field distribution *g*(*y*) to an input function *f*(*y*), consistent with the GF *G*(*y*) associated with a desired operator of choice (Fig. 3A). In this approach, we avoid the need of going into the Fourier domain, hence avoiding the GRIN subblocks that perform Fourier and inverse Fourier transforms. Because this GF slab is transversely symmetric, it can in principle realize an arbitrary GF with even symmetry, such as the second spatial derivative *g*(*y*) ∝ *d*^{2}*f*(*y*)/*dy*^{2}, which can be written in terms of an even kernel as *g*(*y*) ∝ ∫ *f*(*u*)δ′′(*y* – *u*)*du* [where the required GF is proportional to the second spatial derivative of the Dirac delta function δ′′(*y*)]. Our goal is then to tailor the transmission coefficient for impinging plane waves as a function of the transverse wavenumber *k _{y}* to be . In (

*29*), we addressed this problem within the framework of nonlocal transformation optics, assuming infinite media. Here, we consider a finite longitudinal thickness divided into

*N*parallel layers with subwavelength thicknesses

*d*,

_{i}*i*= 1, 2, …,

*N*(Fig. 3A). We developed a fast synthesis approach to obtain the permittivity, permeability, and thickness of each layer required to tailor the overall plane-wave transmission coefficient to match for all impinging angles (

*28*). It is also possible to obtain the same result using nonmagnetic layered metamaterials by suitably increasing the number of layers, arguably simplifying the design and bringing it closer to practical realization. The 10-layered design (

*28*) in Fig. 3B is a GF slab with nonmagnetic (μ = μ

_{0}) constitutive parameters and an overall thickness comparable to λ

_{0}.

Our numerical results confirm that the designed metamaterial provides a GF kernel that resembles, given the unavoidable limited resolution, the desired δ′′(*y*) at the output plane when the input is an approximate delta function. In Fig. 3, B and C, we apply a set of quadratic polynomial functions and the Austin city skyline borders, respectively, as inputs for our second-spatial-derivative GF slab, showing the evolution of the *z*-component of the magnetic field distribution along the slab, and demonstrating that the designed structure indeed approximately performs the second-derivative operation as the wave propagates through it. Edge-detection spatial analog filters [e.g., Laplace filters (*30*)] and image-processing devices may largely benefit from this system. We also design a five-layered GF slab (with overall thickness of 0.37λ_{0}) that performs the convolution operation of an input function with a one-wavelength rectangular spatial kernel (*28*). Figure 3D shows the simulation results using another rectangular function as input.

Our proposed computational metastructures are smaller than conventional lens-based optical signal-processing systems by several orders of magnitude. Moreover, with no analog-to-digital conversion or other systematic delays, here mathematical operations get processed as the electromagnetic signals propagate through the highly miniaturized and compact computational metamaterials. Such designs may potentially lead to direct, ultrafast, wave-based analog computation, equation solving, and signal processing at the hardware level.

## Supplementary Materials

www.sciencemag.org/content/343/6167/160/suppl/DC1

Materials and Methods

Supplementary Text

Figs. S1 to S23

Tables S1 and S2