Inverse-designed metastructures that solve equations

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 g(u)=∫abK(u,v)g(v)dv. An input signal I_in(u) is then introduced into the integral equation (in the format of the Fredholm integral equation of the second kind) by exploiting a set of coupling elements along the feedback waveguides (Fig. 1A). In general, various types of linear and nonlinear integral equations are used across different scientific fields. The Fredholm integral equation of the second kind represents a general case of linear integral equations used in several fields of science and technology, such as in antenna theory (magnetic field integral equation) and in the perturbation theory in quantum mechanics. With this class of examples, we want to provide a roadmap for how to build an equation-solver possibly for other types of integral equations, through properly modifying this approach. For example, one may explore how the solution to a nonlinear Fredholm integral equation of the second kind may be obtained by judiciously introducing nonlinearity in the metamaterial platforms.

• Download high-res image
• Open in new tab
• Download Powerpoint

The general relation governing the system at the steady state then followsg(u)=Iin(u)+∫abK(u,v)g(v)dv(1)where g(u) is the (unknown) solution of this integral equation for the arbitrary input function I_in(u) and the known operator K(u, v). K(u, v) can be nonseperable and complex valued and is in general translationally variant; it can not be expressed in the form of K(u – v).

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 ∫abK(u,v)g(v)dv at N points across the [a, b] interval (36). In this form, the system implements the equivalent N × N matrix equation of the integral equation in Eq. 1. This approach provides four immediate advantages: First, this platform is suitable for direct integration into an optical processor by using currently available silicon photonics technology. Second, the overall footprint of the structure is in the order of few wavelengths, which is orders of magnitude smaller than the conventional optical signal processing schemes based on Fourier transform (23, 31–33). Third, the feedback mechanism and the excitation are fully controlled; both the phase and amplitude of the monochromatic signals are accessible and controllable, circumventing the limitations associated with fiber-based configurations (34) and/or free-space phase retrieval. In general, it is not easy to preserve the phase information of the signal along a fiber-based or free-space setup because the overall size of the system is much larger than an integrated setup. Fourth, because the device entirely operates in the spatial domain, more general translationally variant kernels are straightforwardly accessible and implementable, a feature that is not fully present in Fourier-based designs.

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 e^iω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)K₁(u, v) = 0.25icos(u)exp[–2iuv – 0.5(v – 1)²] + 0.25exp(–u²v²)(2)over [a, b] = [–2, 2]. Elements of this kernel’s scattering matrix, S^K, are calculated based on the distribution of the kernel K₁(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 elementsA=(1.4+0.25i0.25−0.35i−0.35+0.15i−0.3+0.15i−0.5+0.15i00.6−0.2+0.1i−0.65+0.2i0.35−0.3i0.3−0.15i0.3+0.55i1.25+0.2i−0.35−0.25i0.3i0.2+0.5i0.2−0.2i0.25+0.4i1.2−0.1i0.45−0.2i−0.5+0.2i−0.15−0.2i0.55+0.25i−0.30.65+0.1i)(3)In order to compute the inverse of matrix A, A^–1, the equation system A · X = I_in is to be solved for an arbitrary input vector I_in. Following Eq. 1, the relevant kernel operator is thus associated with the transfer matrix (I – A), in which I is the identity matrix of order N = 5. We considered a realistic and physically accessible material design to numerically implement this system. Subsystems in Fig. 1A are each separately designed to operate at microwave frequency with λ₀ ≈ 7 cm (36). We considered conventional rectangular waveguides to carry information to and across the system (via complex amplitude of the out-of-plane component of the electric field in the fundamental TE₁₀ mode of the waveguide). The electrical lengths of all waveguides connecting the kernel to couplers are multiples of the TE₁₀ guide wavelength, ensuring zero phase difference between the two ends of each waveguide. We adapted an inverse design optimization algorithm (41) to design an inhomogeneous permittivity distribution inside a predefined, compact region to yield the transfer properties of the kernel (36). We used the objective-first optimization technique, which has proven advantageous for the design of compact photonic structures such as demultiplexers and couplers (41). In this technique, the transfer function of the kernel is maintained as a constraint, and with each iteration, the relative permittivity becomes successively closer to that which achieves our goal. The results of this process are shown in Fig. 2A, where the metamaterial kernel is designed for the transfer matrix (I – A). The dimensions of the kernel were chosen to be large enough for there to be sufficient degrees of freedom in the optimization process yet still on the order of multiple wavelengths to develop the design in a reasonable amount of time. These dimensions are not unique and may be chosen differently. Mathematically, it is shown to be possible to achieve any arbitrary linear kernel by using a set of input and output waveguides (42); thus, the optimization is expected to converge for general kernels. The physical specifications of the structure are not directly related to the equation parameters u and v, which merely determine the sampling points in the integral equation. For example, the physical positions of the waveguides connecting to the kernel (as in Fig. 2A) can be arbitrarily chosen, even though they represent specific values of u.

• Download high-res image
• Open in new tab
• Download Powerpoint

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 I_in,1 and phase ϕ_in,1, whereas the other four input couplers are not excited. This represents an analog input signal of [I_in,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 isK3(u,v)=0.06{(4−4i)J0(uv)+3exp[i2π5(u+v)]−(1−2i)}(4)defined over [a, b] = [–2, 2]. We provide two distinct arrangements of the feedback system to solve the integral equation associated with K₃(u, v). In general, in our technique the solution of the integral equation is attained through creation of the proper recursive path for the propagation of the electromagnetic waves. The structure in Fig. 1A illustrates the suitable setup to realize an external feedback mechanism. Alternatively, feedback process can be executed internally and in a reflective setup, as conceptually shown in Fig. 3A. In this configuration, the metastructure is designed—by using the inverse design optimization technique (41)—to modify the reflected wave based on the kernel operator. A set of N coupling elements (one in each waveguide) was then used to apply the input signal and probe the response of the system in the steady state. This configuration provides an alternative arrangement to physically implement Eq. 1, and the solution of the integral equation is accessible through the coupling elements (36).

• Download high-res image
• Open in new tab
• Download Powerpoint

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, S^K = (S^K)^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 TE₁₀ mode at λ₀ = 6.85 cm, f₀ = 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, ∫abK(u,v)g(v)dv, is defined as the complex amplitudes of waves reflecting from the kernel (Fig. 3D, red arrows). The frequency of operation of the physical system, size of the metamaterial kernel, and the materials used (the polystyrene) are not special choices dictated by the characteristics of the equation; u and v are dimensionless. In other words, our choices merely facilitated the design and experimental procedure, and the same equation can be implemented and solved at, for example, optical or terahertz frequencies. Each of the five waveguides were terminated with a shorted waveguide stub so that the waves emanating from the kernel were fed back into the kernel with identical phase (they have “zero” electrical length), serving the same role as the feedback loops in the two previous examples. Two coaxial probes were placed on each stub, one to introduce the input signal I_in and the other to sample the solution g(u). In total, 90% of the power entering the coupling elements cycles back into the kernel.

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 I_in = [0, 0, I_in,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.