Solving equations with waves
Signal processing of light waves can be used to represent certain mathematical functions and to perform computational tasks on signals or images in an analog fashion. Such processing typically requires complex systems of bulk optical elements such as lenses, filters, and mirrors. Mohammadi Estakhri et al. demonstrate that specially designed nanophotonic structures can take input waveforms encoded as complex mathematical functions, manipulate them, and provide an output that is the integral of the functions. The results, demonstrated for microwaves, provide a route to develop chip-based analog optical computers and computing elements.
Science, this issue p. 1333
Abstract
Metastructures hold the potential to bring a new twist to the field of spatial-domain optical analog computing: migrating from free-space and bulky systems into conceptually wavelength-sized elements. We introduce a metamaterial platform capable of solving integral equations using monochromatic electromagnetic fields. For an arbitrary wave as the input function to an equation associated with a prescribed integral operator, the solution of such an equation is generated as a complex-valued output electromagnetic field. Our approach is experimentally demonstrated at microwave frequencies through solving a generic integral equation and using a set of waveguides as the input and output to the designed metastructures. By exploiting subwavelength-scale light-matter interactions in a metamaterial platform, our wave-based, material-based analog computer may provide a route to achieve chip-scale, fast, and integrable computing elements.
Equations are the ubiquitous means to express the fundamental characteristics of a system, and solving them is to unravel and predict the behavior of the system. Most scientific and technological phenomena are described through systems of differential and/or integral equations (1, 2), and today, numerous analytical and numerical methods are available to compute and obtain solutions of a given equation (3–5). Computers are among the most popular tools used to simulate, model, and numerically solve systems of equations. However, as computer technology reaches its physical limitations (6), there is an ongoing quest to discover other platforms and disruptive approaches that may afford improved performances, particularly for specialized tasks such as object classification and edge detection (7–9).
Optical analog computing has recently gained renewed attention as an alternative paradigm to contribute to computing technology. Temporal data processors enable high-speed pulse manipulation (10–15), which is relevant for solving differential equations in real time (16, 17). Ultrafast optical systems provide a surprisingly rich platform to model and probe the dynamics of stochastic and complex nonlinear systems (18, 19). Progress in programmable microwave/photonic processors provides the opportunity to implement arbitrary linear transformations and operators, such as spatial differentiators (20–22). Moreover, analog computing in the spatial domain offers computational parallelism, enabling huge amounts of data to be processed simultaneously (23). Recent developments in this field point to the benefits of structured materials—metamaterials (24)—to achieve basic mathematical operations such as spatial differentiation and integration over wavelength-comparable dimensions (25–27). The spatial properties of the artificial medium or surface are suitably engineered to enable computation occurring upon interaction over substantially smaller volumes as compared with conventional Fourier-based devices (23). In extreme scenarios, the simplest configuration, such as a Bragg reflector or metal-dielectric interface, has also been shown to hold some limited form of computational capabilities (9, 28).
We introduce a platform for analog computing in which specially designed metastructures can solve linear integral equations with waves. Information is carried as complex-valued electromagnetic fields, and the solution is attained in the steady state through propagation in a recursive path inside the designed medium. Historically, analog mechanical devices suited to solve differential equations were presented more than a century ago (29, 30). More recently, coherent optical feedback systems (31–33) and fiber-optic networks (34) have emerged as optical computing machines capable of solving integral and differential equations and performing matrix inversion. These proposals, however, face challenges in terms of large size and incompatibility with integrated devices, as well as maintaining phase information over large fiber-optic networks. As will be shown below, our platform does not face these challenges.
A conceptual representation of our idea for solving integral equations with wave propagation in metastructures is shown in Fig. 1A. An integral equation is characteristically an inverse problem, and a feedback mechanism—internal or external—can be used to generate the solution to such an equation in the form of a Neumann series (35). Let us assume that the operator, implemented as a metamaterial block in Fig. 1A, is designed to transfer any arbitrary monochromatic wave g(u) based on a given integral operator with kernel K(u, v). By routing the output of the operator back to the input side, the system can be forced to obey the source-free equation
(A) A sketch of a closed-loop network, consisting of a suitably designed kernel operator (such as a metamaterial block that performs the desired integral operator on the arbitrary input wave), feedback mechanism to establish recursive wave propagation inside the network, and in/out coupling elements to excite and probe the waves. The N waveguides are used to create the feedback network, in this case an external feedback mechanism. Arrows indicate the direction of the wave flow across the structure. The scattering matrices of the blocks can be found in (36). Green and red curves indicate distribution of g(u) and
The general relation governing the system at the steady state then follows
Using the metamaterial kernel, the system is visualized in Fig. 1A in a planar, discretized configuration. We used N feedback waveguides that sample the input of the operator g(u) and its output
We now discuss three examples of equation solving by use of our proposed platform: (i) numerical block diagram analysis of a general integral equation, (ii) design and simulation of a physical structure to calculate inverse of a N × N matrix, and (iii) design and experimental validation of a metamaterial structure that can solve integral equations. Without loss of generality, we focused on one-dimensional linear integral equations (as in Eq. 1). The eiωt time dependence is assumed throughout.
Our first example is based on the ideal block diagram configuration depicted in Fig. 1A. The entire system consists of three distinct subsystems: (i) the prescribed kernel, (ii) feedback waveguides, and (iii) coupling elements. Subsystems can be designed and assembled depending on the frequency of operation and the technological platform of choice (as seen in examples two and three). To decouple the limitations associated with the specific physical implementation of the device from those related to the system-level design, for the first example we considered ideal blocks described with known scattering matrices (37). We chose a general kernel, Eq. 2, that is asymmetric, nonseparable, translationally variant, and complex-valued (Fig. 1B)K1(u, v) = 0.25icos(u)exp[–2iuv – 0.5(v – 1)2] + 0.25exp(–u2v2)(2)over [a, b] = [–2, 2]. Elements of this kernel’s scattering matrix, SK, are calculated based on the distribution of the kernel K1(u, v) over the chosen spatial domain (36). Our implementation also entails reciprocity; the scattering matrix associated with any chosen kernel corresponds to a reciprocal system that facilitates the physical implementation of the solver. The integral equation in Eq. 1 may originate from a physical problem itself (for example, a thermodynamic, biological, or mechanical problem) or may be a purely mathematical one. Our representation of the equation in a material platform (with the complex amplitude of the electromagnetic waves corresponding to the mathematical quantities described by Eq. 1) treats u and v as dimensionless variables, independent from the original concept. Discussions on the characteristics of the scattering matrix, including reciprocity and passivity, can be found in (36).
Once the scattering matrices are known, for any arbitrary input signal we used numerical, system-level simulations (38) to calculate the solution of the integral equation associated with the kernel shown in Eq. 2—fields generated on the waveguides at the points that the wave enters the kernel (Fig. 1A). The solution was then extracted by using ideal directional couplers embedded along the feedback loops, each with a coupling coefficient of 10%. Ideally, probing couplers are expected to impose minimal perturbation to the flow of the wave inside the loop, and the coupling coefficient is desired to be as small as possible. We confirmed , however, that the process is stable, and the solution is highly resilient to such probing leakages (36). Shown in Fig. 1C is an arbitrary input function, which is proportional to the distribution of the complex-valued fields applied to the couplers’ input ports. In this example, N = 20 waveguides are used to sample the waves across the u ∈ [–2, 2] range. The extracted output signal is shown in Fig. 1D (real and imaginary parts; solid lines), corresponding to g(u) in Eq. 1. The solution faithfully follows the same trajectory as the expected exact solution of the corresponding N linear equations, calculated with conventional numerical tools (39) (Fig. 1D, dashed lines). The input signal, which represents the input function of the Fredholm integral equation with the kernel given in Eq. 2, was arbitrarily chosen for illustration of the ability of our system to generate the solution at the steady state. In general, and to assess the performance of the solver, the structure may be excited at all N ports (one by one), as seen in the following two examples.
To demonstrate the performance of this technique in a physical system, for the second example we considered realization of matrix inversion with a designed metastructure. In general, calculating the inverse of a known, nondegenerate, N × N matrix is equivalent to simultaneously solving N linear equations. The process, therefore, is also implementable with our platform because it solves the equivalent system of equations associated with the given linear integral equation in Eq. 1. While it is an interesting mathematical problem by itself, efficient calculation of the inverse of matrices is also of great importance for many applications across physics, mathematics, and computer science. In general, our system can handle any linear problem if the eigenvalues of the scattering matrix associated with the kernel are smaller than unity (40). This condition ensures the convergence of the fields across the system at the steady state (32, 40) and may be straightforwardly satisfied by normalizing the kernel with an appropriate scaling factor related to its eigenvalues [a detailed discussion is provided in (36)]. We considered an arbitrary 5 × 5 asymmetric matrix with complex-valued elements
(A) Inhomogeneous metamaterial kernel designed for the original matrix (I – A, where A is given in Eq. 3). Distribution of the relative permittivity ε(x, y) inside the two-dimensional region enclosed with absorbing layers [perfectly matched layer (PML)] on top and bottom with the thickness of one-wavelength at the initial design frequency, and perfect electric conducting walls on all other sides (dashed lines). Five equidistance waveguide ports are placed on each side of the block. Arrows indicate the direction of wave flow in the structure. A thin air layer is placed between the metamaterial and the absorbing layer to minimize local reflections. (B) Time snapshot of the simulation results for the out-of-plane component of the electric field in the closed-loop network, when excited at port 1, with the other four ports unexcited: Iin = [Iin,1cos(ωt + ϕin,1), 0, 0, 0, 0]. First column of the inverse matrix is built up at the position of vertical black arrows as the complex-valued amplitudes of the fields. The kernel in (A) is highlighted in the orange box. Input and output ports of the couplers are sequentially numbered from 1 to 5 and 6 to 10, respectively. (C) Comparison between the simulated results (blue circles) and the exact theoretical values (gray circles) for the elements of the inverse matrix in the polar complex plane, showing a good agreement.
Illustrated in Fig. 2B are our numerical simulation results for the distribution of the out-of-plane component of the electric field (snapshot in time) in the assembled system, which includes our designed metastructure and the five feedback waveguides and directional couplers. The structure is excited at the first waveguide (through the corresponding coupler) with a monochromatic signal with electric field of amplitude Iin,1 and phase ϕin,1, whereas the other four input couplers are not excited. This represents an analog input signal of [Iin,1cos(ωt + ϕin,1), 0, 0, 0, 0] at ports 1 through 5 marked on Fig. 2B. The complex values of the field distribution constructed on the waveguides (at the position of small vertical arrows in Fig. 2B) represents the first column of A–1 (extracted via ports 6 through 10 of the couplers shown in Fig. 2B). The process may be repeated with excitation on all five ports one by one to obtain all the columns of this inverse matrix. We compared our simulation results with those obtained through a conventional numerical tool (39), as shown in Fig. 2C, which reveals good agreement. Details on the design of couplers, obtaining the solution, and calculating errors can be found in (36).
In the second example, with the operation on the complex-valued field (rather than intensity), we successfully obtained the inverse of the arbitrary complex-valued matrix. The computation time after the implementation of the setup (designing metamaterial kernel and assembling the feedback structure) is linearly dependent on the size of the matrix because we needed to excite the structure on each of the N waveguides to calculate one column of the inverse matrix. The specific metamaterial platform exploited here—inhomogeneous distribution of permittivity ε(x, y)—is not the only possible approach to physically realize the kernel (Fig. 2A). Alternatively, we could use techniques such as self-configuring and programmable optical networks (20–22, 42) and power-flow engineering (43) for the implementation of the same subsystem. Each form of implementation of kernel provides advantages that may suit different applications. Our designs here are fixed (for a specific kernel, but for arbitrary inputs) yet provide a small footprint. Alternatively, reconfigurable networks can provide similar functionalities, yet over larger footprints (20–22, 42). Self-configuring systems can be also used to find the singular value decomposition of a physical system (kernel in our approach), which may be used to calculate the inverse of the transfer function of the physical system (44). By embedding the feedback mechanism inside the system, we directly generated the inverse matrix.
In the third example, we validated the performance of our equation-solving metastructures by designing, numerically simulating, and experimentally demonstrating a proof-of-concept structure at microwave frequencies. Although the design approach is applicable to any wavelength, we chose microwaves to facilitate the fabrication and measurement processes. The kernel of choice is
(A) Conceptual sketch of a closed-loop network in the reflective set-up, consisting of the suitably designed metamaterial kernel and in/out coupling elements to excite and probe the waves. The N waveguides are used to connect the kernel and couplers, forming an internal feedback mechanism. The scattering matrices of the blocks can be found in (36). Green and red lines indicate distribution of g(u) (ingoing wave) and
The advantage of the reflective configuration is threefold: First, in comparison with the transmissive configuration with external feedback lines, the kernel’s footprint is approximately shrinkable by half. Second, the feedback paths are equilength in a planar structure, yielding a better control of their effective electrical size, which is particularly useful for the experimental validation. Third, there are no possibilities for scattering and propagation in the reverse direction, and every available mode in the system is directly used. In general, however, the reflective configuration requires a symmetric kernel so that K(u, v) = K(v, u), or equivalently, SK = (SK)tr. This restriction may be deduced from the reciprocity that exists in conventional linear structures. In the case of asymmetric kernel operators, the reflective configuration may require nonreciprocal physical implementation (which is in general a disadvantage), and therefore, the transmissive configuration with external feedbacks (as in Fig. 1A) would instead be more preferred for asymmetric kernels.
For our third example, and to explore the reflective arrangement, we intentionally chose a symmetric kernel (yet otherwise a general one), Eq. 4. Additionally, the same kernel was also numerically studied by using the transmissive configuration with an external feedback setup, the results for which are reported and compared with this experiment in (36).
A schematic rendering and photographs of the fabricated structure are shown in Fig. 3, B and C. Here, and for our experiment, the inverse design optimization algorithm constraints are properly modified to achieve a binarized metastructure; only two materials are present inside the structure: air and low-loss polystyrene, indicated with light and dark orange, respectively, in Fig. 3D. Five sampling points are used across the u ∈ [–2, 2] range in our kernel designed based on Eq. 4, each represented by the monochromatic wave propagating in a metallic waveguide operating at the fundamental TE10 mode at λ0 = 6.85 cm, f0 = 4.38 GHz.
The input of the kernel, g(u), which is the unknown solution of the integral equation, is defined as the complex amplitudes of the waves incident on the five waveguides (Fig. 3D, green arrows) and the output of the integral operator,
To characterize the performance of the system, first, the scattering properties of the implemented kernel (before mounting it into the feedback loop) were estimated through full-wave simulation and compared with the desired distribution (Fig. 3E). Next, the entire system was assembled and excited with input signals on all five input ports, one by one (simulation results for excitation at port 3 are shown in Fig. 3F). The input signal may, in general, be an arbitrary complex distribution of monochromatic waves at the five input ports. To fully characterize the system, we used five orthogonal input signals (each case consists of exciting one of the five input ports and leaving the other four ports unexcited). Any arbitrary input signal may be attained through linear combination of these “orthogonal” bases. The measured complex-valued fields on the coaxial probes (which are directly proportional to the solution of the equations) were then used to extract the solution (Fig. 3G). Although small variations compared with the ideal theoretical results were observed (Fig. 3G, dashed lines), we attained excellent agreement between theoretical and experimental results. The experimental results are almost indistinguishable from the full-wave simulations, indicating the possibility of removing such small differences through de-embedding the effect of couplers (36). These slight differences are attributed to (i) probing the reflected fields with coupling elements; (ii) variations in the length of connecting waveguides, resulting in nonzero electrical length; (iii) unintended differences in the coupling coefficients; and (iv) minor differences between the implemented and desired kernel, highlighted in Fig. 3E. A detailed discussion on the effect of noise due to variations in different subsystems and the tolerance of the solution is presented in (36).
The structure was numerically simulated by using commercially available CST Studio Suite full-wave simulation software (38). The out-of-plane component of the electric field (snapshot in time) for the given input signal Iin = [0, 0, Iin,3cos(ωt + ϕin,3), 0, 0] is shown in Fig. 3F. The solution was the portion of the wave propagating toward the kernel from the coaxial ports (Fig. 3D, green arrows), extracted via coupling elements [more results can be found in (36)].
Our inverse-designed metamaterial platform provides a powerful tool and methodology for solving the general Fredholm integral equation of the second kind in an analog manner and at the hardware level. Similar configurations may potentially be explored to estimate eigenvalues of an integral operator (a matrix in the discretized form) by adding a controllable amplifier/attenuator and phase-shifter in the feedback path. Adding nonlinearity may also provide another degree of freedom to tackle nonlinear equations. The designs we explored are orders of magnitude smaller than Fourier-based systems, with potentially very low energy consumption. Our numerical analysis (36) indicates that the steady-state response may be achievable in less than 300 cycles of the monochromatic wave used, predicting nanosecond and picosecond estimated times for obtaining the solutions at microwave and optical frequencies, respectively. Implementing the technique on a silicon photonic platforms may lead to chip-scale, ultrafast, integrable, and reconfigurable analog computing elements.
Supplementary Materials
www.sciencemag.org/content/363/6433/1333/suppl/DC1
Materials and Methods
Supplementary Text
Figs. S1 to S23
Movies S1 and S2
This is an article distributed under the terms of the Science Journals Default License.
References and Notes
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵Materials and methods are available as supplementary materials.
- ↵
- ↵CST Design Studio; We used the schematic solver (www.cst.com/products/cstds) to model ideal blocks on the basis of their scattering matrixes. We used the full-wave solver (www.cst.com/products/csts2) for modeling the metamaterial structures.
- ↵Matlab (www.mathworks.com/products/matlab.html) and Mathematica (www.wolfram.com/mathematica) numerical solvers.
- ↵
- ↵
- ↵
- ↵
- ↵