††thanks: These authors contributed equally††thanks: These authors contributed
equally††thanks: Corresponding author<EMAIL_ADDRESS>
# Inverse-designed Metastructures together with Reconfigurable Couplers to
Compute Forward Scattering

Vahid Nikkhah  University of Pennsylvania, Department of Electrical and
Systems Engineering, Philadelphia, PA 19104 U.S.A.  Dimitrios C. Tzarouchis
University of Pennsylvania, Department of Electrical and Systems Engineering,
Philadelphia, PA 19104 U.S.A.  Ahmad Hoorfar  Villanova University, Department
of Electrical and Computer Engineering, Villanova, PA 19085 U.S.A.  Nader
Engheta  University of Pennsylvania, Department of Electrical and Systems
Engineering, Philadelphia, PA 19104 U.S.A.

###### Abstract

Wave-based analog computing in the forms of inverse-designed metastructures
and the meshes of Mach–Zehnder interferometers (MZI) have recently received
considerable attention due to their capability in emulating linear operators,
performing vector-matrix multiplication, inverting matrices, and solving
integral and differential equations, via electromagnetic wave interaction and
manipulation in such structures. Here, we combine these two platforms to
propose a wave-based metadevice that can compute scattered fields in
electromagnetic forward scattering problems. The proposed device consists of
two sub-systems: a set of reconfigurable couplers with a proper feedback
system and an inverse-designed inhomogeneous material block. The first sub-
system computes the magnitude and phase of the dipole polarization induced in
the scatterers when illuminated with a given incident wave (matrix inversion).
The second sub-system computes the magnitude and phase of the scattered fields
at given detection points (vector-matrix multiplication). We discuss the
functionality of this metadevice, and through several examples, we
theoretically evaluate its performance by comparing the simulation results of
this device with full-wave numerical simulations and numerically evaluated
matrix inversion. We also highlight that since the first section is
reconfigurable, the proposed device can be used for different permittivity
distributions of the scatterer and different incident excitations without
changing the inverse-designed section. Our proposed device may provide a
versatile platform for rapid computation in various scattering scenarios.

##  I Introduction

Optical analog computing shows great promise for next-generation computing
platforms. Owing to some of the features of light-matter interaction, optical
metastructures gives way for low-power, and high-speed analog computations [1,
2, 3, 4, 5]. The artificially-engineered electromagnetic structures known as
metamaterials [6] combined with inverse-design techniques enable engineering
of wave propagation through material structures down to subwavelength scales
with unprecedented efficiency [7, 8]. These structures are capable of doing
computation with near the speed of light for complex computational problems
such as solving integral equations [9]. Moreover, the wave propagation through
material structures exhibits parallelism in terms of encoding and processing
information [10]. Considering the capabilities of optical analog computing,
here we propose the idea of wave-based metadevices for performing specialized
computational tasks, particularly for certain physical problems such as the
electromagnetic forward-scattering problems, which can be computationally very
intensive.

In this work, we encode the information about the quantities of interest,
i.e., the complex-valued electric field and induced polarization in our
scattering problem, into the complex amplitudes of the waveguide modes. Then
we process them by a metadevice, schematically shown in Fig. (1), that
consists of two sub-systems: (a) a set of reconfigurable couplers, e.g.,
Mach–Zehnder interferometers (MZIs), and (b) an inverse-designed
metastructure. In the first part, a network of reconfigurable couplers is
designed, which when endowed with a proper feedback loop, can compute the
inverse of a given matrix. This part is designed to compute the magnitude and
phase of the dipole polarization induced inside the target scatterers by the
local electric field, for a given incident field. As discussed below, since
these couplers are tunable the permittivity of implemented scatterers and the
incident fields can be modified at will. The second part is an inverse-
designed metastructure that can perform a matrix-vector product with waves for
computing the scattered fields, generated by the polarization induced in the
target, at the designated detection points using the Green’s function as the
propagator. This part of our metadevice is formed using an inhomogeneous
distribution of dielectric materials, whose optimized distribution is obtained
using the method of inverse-design. Finally, we validate the performance of
the proposed design through simulations of a series of forward scattering
examples.

The article is organized as follows: In section II, a brief overview of the
mathematical formalism for solving forward scattering problems using the
method of discrete dipole approximation (DDA) is given. In section III, the
details of the metadevice design are presented, i.e., the tunable couplers and
the inverse-designed stages. We then test, via simulation, the proposed device
against two particular examples of a forward scattering problem in section IV,
i.e., (a) an example of a scattering problem with various incident fields and
various permittivity profiles and (b) an example of a through-the-wall radar
imaging problem. The summary and the concluding remarks are given in section
V.

##  II Formulation Overview

We start with the standard integral representation of the electromagnetic
scattering for a two-dimensional (2D) problem:

$\mathbf{E}(\mathbf{r})=\mathbf{E}_{\text{inc}}(\mathbf{r})+\int_{S}\mathbf{G}(\mathbf{r},\mathbf{r}^{\prime})\chi(\mathbf{r}^{\prime})\mathbf{E}(\mathbf{r}^{\prime})d^{2}\mathbf{r}^{\prime}$
(1)

where $\mathbf{E}(\mathbf{r})$ and $\mathbf{E}_{\text{inc}}(\mathbf{r})$ are,
respectively, the total and the incident electric fields at the observation
point $\mathbf{r}$, $\mathbf{G}(\mathbf{r},\mathbf{r}^{\prime})$ is the 2D
Green’s function of the Helmholtz’s equation and
$\chi(\mathbf{r}^{\prime})=\varepsilon_{0}(\varepsilon_{r}(\mathbf{r}^{\prime})-1)$
is the susceptibility (sometimes referred to as contrast or scattering
potential) of the object with relative permittivity
$\varepsilon_{r}(\mathbf{r}^{\prime})$ and enclosed by the cross-sectional
surface $S$ as shown in Fig. (1) (a). Here we have assumed that the
surrounding medium is air with permittivity $\varepsilon_{0}$. Indeed, Eq. (1)
is a Fredholm integral equation of the 2nd kind that is often referred to as
the Lippmann–Schwinger equation [11] and lies within the heart of many
numerical electromagnetic methods, such as the discrete dipole approximation
(DDA) [12, 13], the contrast source inversion method (CSI) [14, 15], the
coupled dipole method [16], and the Method of Moments (MoM) [17].

We need to express Eq. (1) in terms of its equivalent discrete version using
matrix/vector quantities. For this purpose we closely follow the DDA approach,
a general method for solving electromagnetic scattering problems of dielectric
targets of arbitrary shapes and compositions [18]. The domain $S$ that
contains the 2D dielectric object is discretized in small scattering cells
$i=1,2,3,\dots N$, where each 2D cell exhibits an induced dipole moment
$\mathbf{p}_{i}$ that is generally proportional to the geometrical shape (here
surface) and the associated permittivity at the $i$-th cell (assuming
piecewise homogeneous scatterers) and the local electric field (Fig. (1) (a)).
The evaluation of the unknown induced dipole moments is the first goal of our
approach.

Figure 1: Proposed metadevice that emulates the forward scattering problem:
(a) the physical problem that consists of the incident excitation (green
arrows), the scattered field (red) and the scattered field obtained at a given
set of detection points. Computing the scattered fields requires the matrix
inversion for the determination of the induced dipole polarization vectors at
the scatters due to the incident excitation, and a matrix-vector operation for
the evaluation of the scattered field at the desired detection domain. (b) Our
proposed device consists of the reconfigurable direct-complex-matric (DCM)
architecture involving reconfigurable multipliers (for matrix inversion) and
the inverse-designed metastructure (for vector-matrix multiplication). The
first emulates the interaction between the induced dipoles of the scatterers
(following DDA), whereas the second evaluates the propagator.

Each cell acquires its dipole moment due to the local electric field, i.e.,

$\mathbf{p}_{i}=\alpha_{i}\mathbf{E}^{\text{loc}}_{i}$ (2)

where the $\mathbf{E}_{i}^{\text{loc}}$ is the local electric field at the
center of the cell $i$ and $\alpha_{i}$ is the polarizability that depends on
the shape and the material composition of each 2D cell. In turn, the local
field is generated by the incident field and the secondary fields generated
from all the rest of the dipoles (other discrete cells) such that:

$\mathbf{E}_{i}^{\text{loc}}=\mathbf{E}_{i}^{\text{inc}}+\sum_{k\neq
i}^{N}\mathbf{G}_{ik}\mathbf{p}_{k}$ (3)

where $\mathbf{E}^{\text{inc}}_{i}$ is the incident field at the target point
$i$, $\mathbf{p}_{k}$ is the polarization induced inside the target cell $k$
and $\mathbf{G}_{ik}$ is the 2D Green’s function for calculating the fields,
generated from $\mathbf{p}_{k}$, at the location of the target point $i$. In
our case we consider a two-dimensional (2D) problem with a transverse electric
(TE) excitation. Therefore the corresponding Green’s function reads

$\mathbf{G}_{ik}=\mathbf{G}(\mathbf{r}_{i}-\mathbf{r}_{k})=-j\frac{k_{0}^{2}}{4\pi\varepsilon_{0}}H^{(2)}_{0}(k_{0}|\mathbf{r}_{i}-\mathbf{r}_{k}|)$
(4)

where $H^{(2)}_{0}(k_{0}|\mathbf{r}_{i}-\mathbf{r}_{k}|)$ is Hankel function
of the second type and 0-th order and
$k_{0}=\omega_{0}\sqrt{\mu_{0}\varepsilon_{0}}$ is the free-space wavenumber
[19] (here we use the $e^{j\omega t}$ convention). By combining Eqs. (2) and
(3) we get the following expression:

$\alpha^{-1}_{i}\mathbf{p}_{i}=\mathbf{E}^{\text{inc}}_{i}+\sum_{i\neq
k}^{N}\mathbf{G}_{ik}\mathbf{p}_{k}$ (5)

for which the unknown quantity is the dipole polarization induced inside the
cells, i.e. $\mathbf{p}_{i},\quad i=1,..,N$ . Using simple algebra, the above
system can be arranged using the matrix formulation as follows

$p=\left(\text{diag}(\alpha^{-1})-G\right)^{-1}e_{\text{inc}}$ (6)

where the lowercase quantities
$p=[\mathbf{p}_{1},\mathbf{p}_{2},...,\mathbf{p}_{N}]^{T},\alpha=A_{\text{cell}}\varepsilon_{0}[\varepsilon_{1}-1,\varepsilon_{2}-1,...,\varepsilon_{N}-1]$
($A_{\text{cell}}$ is the cross-sectional area of a cell) and
$e_{\text{inc}}=[\mathbf{E}^{\text{inc}}_{1},\mathbf{E}^{\text{inc}}_{2},...,\mathbf{E}^{\text{inc}}_{N}]^{T}$
are $\mathbb{C}^{N\times 1}$vectors, $\text{diag}(\cdot)$ is the diagonal
matrix operator, and $G$ is a $\mathbb{C}^{N\times N}$ Toeplitz matrix with
zero diagonal entries (assuming a uniformly spaced discrete grid) [20].

Finally, the scattered field observed at $M$ specified discrete detection
points (in general $M\neq N$) is given by:

$e_{\text{sca}}=G^{\text{pr}}~{}p=G^{\text{pr}}\left(\text{diag}(\alpha^{-1})-G\right)^{-1}e_{\text{inc}}$
(7)

where $G^{pr}\in\mathbb{C}^{M\times N}$ is the “propagator” Green’s function
matrix connecting the induced dipole polarization vectors of subscatterers
(cells) and the detection points. Based on the above matrix representations of
the scattering problem we propose a specially designed wave-based
reconfigurable metadevice that is able to simultaneously solve Eq. (6) and
implement the matrix-vector operation of Eq. (7) for different excitations and
for different scattering scenarios.

##  III Metadevice design

In this section we present the main results of our study. As sketched in Fig.
(1) (b), our proposed metadevice consists of two parts: (A) the matrix
inversion part and (B) the vector-matrix multiplication part. For the first
part we utilize a recently introduced Direct Complex Matrix (DCM)
reconfigurable architecture [21] properly connected to a feedback loop,
capable of performing the desired matrix inversion in Eq. (6). For the second
part we propose an inverse-designed metastructure that is able to conduct the
matrix-vector product for emulating the propagator matrix, $G^{\text{pr}}$ in
Eq. (7).

###  III.1 Reconfigurable DCM network stage for inverting matrices

This stage, shown as the left panel of Fig. (1) (b), follows the architecture
introduced in [21]. In particular, a given $N\times N$ complex-valued matrix
can be implemented in this device for inversion that has three sections: (a)
the ingress, (b) the middle, and (c) the egress sections. The ingress section
consists of $N$ number of $1\to N$ signal splitters, each of which splits each
of the $N$ inputs into $N$ outputs. All such $N^{2}$ outputs of the ingress
section are then properly routed to the middle section that consists of
$N^{2}$ dedicated tunable components, each of which consists of a combination
of a phase shifter and an amplifier/attenuator, which we call a “multiplier”.
These $N^{2}$ multipliers directly represent $N^{2}$ elements of the matrix in
Eq. (6) to be inverted. Finally, the $N^{2}$ outputs of the middle sections
are then routed to the egress section, i.e., $N$ number of $N\to 1$
combiners111As an aside, it is worth pointing out that such combiners are not
power combiners. Instead, the output of each $N\to 1$ combiner is proportional
to the fraction of the algebraic sum of complex-valued amplitudes of $N$ input
signals, and consequently the output power is usually less than the sum of the
input power; As a result, the remaining power goes to some “loss channel”
e.g., either dissipates in the element, scattered away from the system,
reflected, or goes out of the system in unused channels. Consequently, we need
to have amplifiers to bring the output signal level up to the desired value
for the feedback system. As such the egress section is basically inverse of
the egress section. Moreover, there is a feedback loop connecting the outputs
to the inputs.

For the feedback loop section we used two identical couplers that split the
input and output signal with a coupling ratio $\beta_{1,2}$ (and the through
is $|\alpha_{1,2}|^{2}=1-|\beta_{1,2}|^{2}$, see Fig. (1) (b)) and the
supporting information. The output of the DCM stage estimates a value that is
proportional to the induced dipole polarization vector $p$. In order to
retrieve the actuall value of the polarization vector we need (a) to scale the
input vector $e_{inc}$, i.e., $e_{\text{inc}}/\beta_{1}$, and (b) scale the
output, i.e., $p_{\text{DCM}}/\beta_{2}$. Therefore, from Eq. (6) the total
scaling factor $\frac{1}{\beta_{1}\beta_{2}}$ should be applied to the DCM
output. Notice that in all examples presented below we assume couplers with
$-30$dB coupling ($|\alpha_{1}|^{2}=|\alpha_{2}|^{2}=0.999$ and
$|\beta_{1}|^{2}=|\beta_{2}|^{2}=0.001$) requiring a scaling factor of
$\approx 10^{3}$ to be applied at the output of this stage for the retrieval
of the induced dipole polarization vector. (See the supplementary information
for evaluation of dimensionless quantities in our analysis.)

Overall, such device can implement the desired matrix inversion, required in
Eq. (6) [21]. Since the $N$ diagonal terms of this matrix depend on the
relative permittivities of $N$ scattering cells, these diagonal modules of the
middle section can be reconfigured by changing the phase and magnitude of
their transfer functions (using the phase shifters and the
amplifiers/attenuators in them) in order to accommodate different
permittivities for these scatterers. In short, when the complex values of the
incident fields at the location of $N$ scatterers are inserted into the DCM
stage of our metadevice, we evaluate the induced polarization vectors for
different scattering scenarios represented by different material properties,
e.g., permittivities, of the individual scatterers and different incident
excitations.

###  III.2 Inverse-designed metastructures for emulating propagator matrix

In this section, we utilize the method of inverse design [22, 23, 24, 25, 26,
27, 28], which is a density-based topology optimization technique using a
gradient-descent approach implemented by adjoint sensitivity analysis in
COMSOL, to propose the wave-based inhomogeneous metastructures that emulate
the propagator matrix $G^{\text{pr}}$ (see also supporting information).
Figure 2(a) presents an example of such a structure that represents the
propagator matrix for a 2D scattering example, shown in Fig. (3) (a), with 5
scatterers and 7 detection points. This structure may consist of a metallic
box with perfect electric conducting (PEC) walls, with 5 single-mode
waveguides as the input ports and 7 single-mode waveguides as its output
ports. The inputs are excited with the complex amplitudes of the induced
dipole polarization vectors (which are the outputs of the reconfigurable DCM
network), and the complex-valued scattered fields at the detection points can
be obtained (after applying the scale factor) at the outputs of this
metastructure, according to $e_{\text{sca}}=G^{\text{pr}}p$ (see Fig. (2)
(b)).

The size of the design domain is $6\lambda_{0}\times 9\lambda_{0}\times
0.3\lambda_{0}$ ($\lambda_{0}$ is the free space wavelength corresponding to
the design frequency $f_{0}=3^{\text{GHz}}$). The size of the waveguides cross
section is $0.6\lambda_{0}\times 0.3\lambda_{0}$ supporting $\text{TE}^{10}$
propagating mode only. Since the $G^{\text{pr}}$ is not a unitary matrix, a
waveguide on top of the metastructure takes the excess energy to outside. This
waveguide has a wider width of $5\lambda_{0}$ to support multiple propagating
modes for exiting the excess energy. The optimized material distribution
inside the design region, for the example shown in Fig. (3) (a), is shown in
Fig. (2) (a). The dark grey regions are assumed to be polystyrene with
$\varepsilon_{r,\text{poly}}=2.53$ and the light grey regions are air. The
material distribution is only a function of $(x,y)$ and it is invariant along
the $z$ direction. Fig. (2) (b) shows the norm of the simulated electric field
distribution inside the optimized structure when the input port 2 is excited,
for the example shown in Fig. (3) (a).

Since the elements of $G^{\text{pr}}$ only depend on the relative distances
between locations of scattering cells and detection points, as long as the
locations of the cells and the detection points are unchanged the designed
metastructure can be used to compute the scattered fields for various forward
scattering problems with different relative permittivity of the scatterers
and/or various types of incident excitations.

In the next section, the simulation results for our proposed metadevice, which
is the combination of the DCM stage and inverse-designed metastructure, are
presented and evaluated for computing forward scattering in several examples.
We show that the metadevice is capable of computing the complex amplitudes of
the scattered fields at the detection points, for different incident
radiations and various relative permittivities of the target scatterers.

Figure 2: Inverse-designed metastructure: Inverse-designed structure designed
for $f_{0}=3^{\text{GHz}}$ comprising 5 input and 7 output ports for emulating
the propagator matrix, $G^{\text{pr}}$, of the configuration shown in Fig. (3)
(a). The size of the design region is
$(L_{x},L_{y})=(6\lambda_{0},9\lambda_{0})$ and $0.3\lambda_{0}$ along $z$
direction. The material distribution is a function of $(x,y)$ only. The input
and output ports are single mode $\text{TE}^{10}$ waveguides. A multimode
waveguide on top of the structure lets the excess power out (a) The optimized
distribution of polystyrene with $\varepsilon_{r}=2.53$ (dark grey regions)
and air with $\varepsilon_{r}=1$ (light grey regions) (b) The magnitude of
simulated electric field distribution inside the inverse-designed
metastructure when the input port 2 is excited.

##  IV Examples and Discussion

As the first example, depicted in Fig. (3) (a), we consider the 2D scattering
from 5 infinitely long cylindrical dielectric rods (parallel with the $z$
axis) that are separated by $\lambda_{0}/2$, placed along a linear grid on the
$y$ axis. The diameter of each rod is chosen to be deeply subwavelength,
$d=\lambda_{0}/50$, to model each rod with a polarizability of
$\alpha_{i}=A_{\text{rod}}\varepsilon_{0}(\varepsilon_{i}-1)$ where
$A_{\text{rod}}=\pi d^{2}/4$ is the cross-sectional area of the rod and
$\varepsilon_{i}$ is the relative permittvity of the rod $i,\quad
i=1,2,...,5$. In this example, the relative permittivities of the rods are
chosen to be $\varepsilon_{r}=\\{2,3,4,3,2\\}$ as shown in the figure. An
infinitely long z-oriented electric line source, located $\lambda_{0}/2$ away
from the center of the array, is generating the incident electromagnetic field
with an out-of-plane E-field polarization. The scattered fields generated from
these 5 rods are to be evaluated at 7 designated detection points positioned
on a quarter section of a circle with radius of $2\lambda_{0}$ as shown in
Fig. (3) (a). As we discuss below our proposed metadevice can evaluate the
scattered fields from this structure, evaluated at the detection points. The
DCM stage of the metadevice is designed to allow the inversion of the matrix
containing the information about the relative permittivities and locations of
these 5 scatterers (see Eq. (6)), and the inverse-designed metastructure
section is designed to emulate the propagator part of Eq. (7) for this
scattering problem. Fig. (3) (c) shows this inverse-designed metastructure.
The DCM stage is simulated using the AWR Microwave Office® and the wave
interaction in the inverse-designed metastructure is simulated using the
COMSOL Multiphysics® [29] (see the supporting information). As depicted in
Fig. 1 (b), the complex-valued incident electric fields at the locations of 5
rods are considered as the inputs to the DCM stage, and the complex-valued
induced dipole polarization amplitudes excited from the DCM stage are the
inputs to the inverse-designed metastructures. Fig. (3) (b) depicts the
magnitude of the scattered electric field, obtained from full-wave simulations
using the COMSOL Multiphysics® [29]. The solid curve indicates this magnitude
on the quarter circle of radius of $2\lambda_{0}$ where the detection points
are located. Moreover, we show (as blue triangle) the theoretical results for
such scattered field amplitude as obtained from the analytical expression (Eq.
(7)) using Matlab ®, and the simulation results from the outputs of our
proposed metadevice (shown as red circles). We note the good agreement among
all three approaches (for the phase information and the comparison, see the
supporting information). Fig. (3) (c), shows the electric field distribution
inside the inverse-designed metastructure and the relative mode amplitudes at
the output ports (the results of which are shown by red circles in Fig.3 (b)).

Figure 3: Example 1 and its simulation results (a) Geometry of the forward
scattering problem under investigation as the first example. 5 two-dimensional
(2) D dielectric rods, with subwavelength cross sections, are illuminated by
an infinitely long line electric current source with the out-of-plane electric
field. The scattered fields are obtained at 7 designated detection points on
the circular path highlighted by green (b) The scattered field distribution
computed by full-wave simulation using the COMSOL Multiphysics®. The scattered
field data on the detection path obtained by the full-wave simulation (black
solid curve), from the DDA calculation in Matlab® (blue triangles), and from
the simulated outputs of our proposed metadevice (with AWR Microwave Office,
red circles) are compared on a polar plot (c) Inverse-designed metastructure
for emulating the propagator Green function. The inputs to this structure are
excited with the complex amplitudes proportional to the p vector obtained by
the DCM architecture. The mode amplitudes at the outputs are proportional to
the scattered field complex amplitudes at the detection points

Since our metadevice is designed based on Green’s function formulation in Eq.
(7), the incident field could be arbitrarily chosen. Therefore, the same
metadevice is capable of computing the scattered fields for all sorts of
incident radiation, without a need to be re-designed. For demonstrating this
capability, as the second set of examples here, we displace the line source
and evaluate the new scattered fields. Fig. (4) (a) & (b) show the results for
the line source displaced $0.2\lambda_{0}$ and $0.4\lambda_{0}$ downward
relative to the center of the array, respectively. The panels next to the
scattered field results show the electric field inside the inverse-designed
metastructure and the mode amplitudes at the output ports. The good agreement
among the scattered field data confirms the ability of our metadevice to
evaluate scattering fields for this scattering problem for an arbitrary
incident radiation. (For the phase information, see the supporting
information.)

As mentioned earlier, our DCM stage is reconfigurable, and consequently it can
be adjusted for different values of relative permittivities of the rods as
long as the rods stay at their given locations. This latter condition, when
added to the fact that the detection points are also fixed, guarantees that
the inverse-designed metastructure can be re-used without any need for
redesigning. This allows our metadevice to compute the scattered fields at the
detection points for arbitrary relative permittivities of rods and arbitrary
incident radiation. This capability could be useful in optimization problems
where one needs to evaluate the EM scattering for numerous different
distributions of the relative permittivity of the object, for example in the
ground-penetrating radar (GPR) and subsurface scenarios involving inverse-
profiling of composite targets [30, 31]. As the 3rd set of examples, for the
same geometrical configuration of Fig. (3) (a), we consider
$\varepsilon_{r}=\\{4,3,2,3,4\\}$ as a new set of dielectric constants for the
rods. Next, we modify some of the reconfigurable diagonal elements in the
middle stage of our DCM architecture to implement the new matrix to be
inverted in order to provide us with a new set of the polarization induced
inside the objects according to their given relative permittivities. Fig. (4)
(c) shows the scattered field results for the new set of relative
permittivities of the objects being illuminated by a line source
$\lambda_{0}/2$ away from the center of the array. Fig. (4) (d) shows the
corresponding results for anther set of relative permittivities of
$\varepsilon_{r}=\\{-4+j0.02,\ -2+j0.01,\ j0.01,\ -2+j0.01,\ -4+j0.02\\}$.
Again good agreement is noted among the scattered field magnitude results for
both sets of relative permittivities, hence, verifying the capability of our
metadevice to compute the scattered fields in forward-scattering problems with
various relative permittivity profiles of the target scatterer. (for the phase
information and the comparison, see the supporting information.)

Figure 4: Examples 2 (first row) and 3 (second row): Different incident
excitation and different sets of relative permittivities of the rods. Similar
to the example of Fig. (3), here we have the simulations results for the
scattered fields (blue triangle Matlab® results, red circle AWR Microwave
Office® simulation results, black line full wave simulations with COMSOL
Multiphysics®) for different incident radiations produced by a line source
that is displaced by (a) $0.2\lambda_{0}$ and (b) $0.4\lambda_{0}$ downward
relative to the center of the array. The second row shows the corresponding
results for the modified relative permittivities for (c)
$\varepsilon_{r}=\\{4,3,2,3,4\\}$ and for (d) $\varepsilon_{r}=\\{-4+j0.02,\
-2+j0.01,\ j0.01,\ -2+j0.01,\ -4+j0.02\\}$.

As the 4th example, we demonstrate the capability of our metadevice for
performing forward-scattering computations for through-the-wall imaging
applications. The through-the-wall radar imaging (TWRI) is a technique for
making images of targets behind visually opaque dielectric walls using
microwaves. The scattered fields from targets propagate back to the receiver
for post-processing [32, 33, 34, 35, 36]. One of the challenges of the TWRI is
the computation time (particularly when many iterations in the forward
scattering are needed in such computation), which is a determining factor in
the applicability of an imaging method for real-time applications [37, 32].
Here we show how an extension of our proposed metadevice can be utilized to do
the forward-scattering computation for the TWRI. Since our device computes
with waves (and thus may exhibit high wave speed), the required time for image
reconstruction may be reduced considerably, especially if one is also
interested in estimating the permittivity profiles of the targets behind or
within the wall. As an aside, it is worth noting that in inverse profiling
(e.g., GPR and subsurface applications), where one is concerned with finding
both location and shape of an object as well as its material composition, the
forward scattering problem becomes time consuming, even when the target
location is fixed [38]. In such scenarios, our proposed meta-device may become
even more useful in reducing computation time.

A typical EM model of through-the-wall forward scattering for a single-layer
dielectric wall is shown in Fig. (5) (a). A 2D infinitely long line source
with an out-of-plane electric field illuminates the scene and the detection
points are positioned close to the source in a linear configuration, for
receiving the scattered fields from the 2D target scatterers behind the wall.
In order to use the far-field approximation of the wall’s Green function the
line source, scatterers, and the detection points are located several
wavelengths away from the wall. The far-field Green function for a single-
layer dielectric wall is as follows:

$\mathbf{G}^{\text{w}}(\mathbf{r}-\mathbf{r}_{s})\approx
T_{W}\mathbf{G}^{\text{pr}}(\mathbf{r}-\mathbf{r}_{s})=-jT_{W}\frac{k_{0}^{2}}{4\pi\varepsilon_{0}}H^{(2)}_{0}(k_{0}|\mathbf{r}-\mathbf{r}_{s}|)$
(8)

where $\mathbf{G}^{\text{w}}$ and $\mathbf{G}^{\text{pr}}$ are the 2D Green
function in the presence of the wall and the 2D Green function in free-space,
respectively. $T_{\text{w}}$ is the wall’s transmission coefficient for a
normal incident plane wave and $\mathbf{r}_{s}$ is the position of the line-
source. The transmitted fields through the wall illuminate the objects behind
it and the scattered fields propagate back through the wall to reach the
detection points. Therefore, using the notation we used earlier, we may write
the following formulation for calculating the scattered fields at the
designated detection points:

$e_{\text{sca}}=T_{\text{w}}G^{\text{pr}}~{}p=T_{\text{w}}G^{\text{pr}}\left(\text{diag}(\alpha^{-1})-G\right)^{-1}T_{\text{w}}e_{\text{inc}}$
(9)

Compared to the Eq. (7) for a scattering problem in a free-space background,
the only change we need to make for the present case with the wall is to
include the transmission coefficient twice, one for the incident wave going
through the wall to reach the scatterers and the other is for the scattered
fields from the objects to reach the detection points. We note that the Green
function in Eq. (8)-(9) may be appropriately modified to redesign and tailor
the metadevice for other forward scattering problems of interest such as those
in intra-wall imaging [32] on in GPR subsurface profiling [39].

In Fig. (5) (a), we consider a case of through-the-wall forward scattering.
The thickness and relative permittivity of the wall are assumed to be
$0.75\lambda_{0}$ and $\varepsilon_{r}=9$, respectively. Five 2D dielectric
rods with diameter of $d=\lambda_{0}/50$ are positioned as a linear array
$10\lambda_{0}$ away from the wall. The distance between adjacent cylinders is
$\lambda_{0}/2$ and the relative permittivities of the rods are
$\varepsilon_{r}=\\{2,\ 5,\ 3,\ 6,\ 1.5\\}$. Seven detection points with
spacing of $\lambda_{0}/2$ are on the other side of the wall close to the
source. The source is positioned $\lambda_{0}$ away from the detection points.
Fig. (5) (d) depicts the results of the full-wave simulation for the scattered
field distribution using COMSOL Multiphysics ® [29]. Fig. (5) (b) presents the
inverse-designed metasrtcuture for emulating the far-field Green function with
the wall present in the configuration shown in Fig. (5) (a). Fig. (5) (c)
shows the magnitude of scattered fields at the detection points using the
COMSOL Multiphysics ® full-wave simulation (black solid curve), the numerical
results of Eq. (9) implemented in Matlab® (blue triangles), and the simulation
results for our metadevice (red circles). Close agreement is observed,
verifying the capability of our proposed metadevice for performing forward-
scattering computations for through-the-wall imaging. (For phase information
and the comparison, see the supplementary information.)

Figure 5: Through-the-wall-radar imaging (TWRI) application: (a) The geometry
of a through-the-wall forward scattering problem. 5 dielectric rods with
spacing of $\lambda_{0}/2$ and relative permittivities of
$\varepsilon_{r}=\\{2,5,3,6,1.5\\}$ from top to bottom, respectively, are
hidden behind a dielectric wall with $\varepsilon_{r}=9$ and thickness
$d=0.75\lambda_{0}$. An infinitely line electric line source illuminates the
scene and the scattered fields are detected at 7 designated detection points
with spacing of $\lambda_{0}/2$. The line source, detection points, and the
scatterers are several wavelengths ($\sim 10\lambda_{0}$) away from the wall
in order to use the far-field Green’s function approximation in the presence
of the wall. (b) the inverse-designed metastructure for emulating such Green’s
function, $\mathbf{G}^{\text{w}}$, for this problem (a). The input ports are
excited with the complex amplitudes proportional to the induced diepole
polarizations computed from the DCM stage. (c) The scattered field data at the
detection points are evaluated with full-wave simulations in COMSOL
Multiphysics® (black solid curve), in Eq. (9) in Matlab® (blue triangles), and
in the outputs of the metastructure simulated in AWR Microwave Office® (red
circles) (d) The magnitude of the scattered electric field distribution
simulated for the EM model shown in panel(a).

##  V Summary and Conclusion

In this work, we proposed a wave-based metadevice for computing
electromagnetic forward scattering problems. The scattering object on which
the integral equation is defined is discretized into small cells and the
induced dipole polarizations inside the cells can be computed by inverting a
matrix. We proposed a mesh of tunable couplers based on the DCM architecture
for computing the matrix and unity feedback on top of it for inverting it.
Then the computed polarizations are fed into an inverse-designed metastructure
for determining the scattered fields at the designated detection points using
the propagator Green’s function. The proposed metadevice is capable of
computing the scattered fields for various scattering problems with different
relative permittivity profiles of the object and various different incident
excitations. Although our proposed platform is customized for solving
electromagnetic scattering problems, we can use the general concept behind
this architecture to design other wave-based metastructure for solving other
physical problems. One of the interesting potentials of this wave-based
computing platform is the possibility for enhanced computation speed offered
by the high speed of electromagnetic waves. Therefore, in applications in
which numerical computation of a specific physical problem needs to be
performed many times, our architecture may decrease the wait time for
computations significantly.

## Acknowledgements

This work is supported in part by the US Air Force Office of Scientific
Research (AFOSR) Multidisciplinary University Research Initiative (MURI) grant
number FA9550-17-1-0002, and in part by the US National Science Foundation
(NSF) MRSEC program under award No. DMR-1720530.

## Notes

The authors declare no competing financial interests. N.E. is a strategic
scientific advisor/consultant to Meta Materials, Inc. The numerical data,
codes, and other data that support the findings of this study are available
from the authors upon reasonable request.

## References

  * [1] D. A. B. Miller, “The role of optics in computing,” Nature Photonics, vol. 4, no. 7, p. 406–406, 2010.
  * [2] D. A. B. Miller, “Attojoule optoelectronics for low-energy information processing and communications,” Journal of Lightwave Technology, vol. 35, no. 3, p. 346–396, 2017.
  * [3] H. J. Caulfield and S. Dolev, “Why future supercomputing requires optics,” Nature Photonics, vol. 4, no. 5, p. 261–263, 2010.
  * [4] W. Bogaerts, D. Pérez, J. Capmany, D. A. B. Miller, J. Poon, D. Englund, F. Morichetti, and A. Melloni, “Programmable photonic circuits,” Nature, vol. 586, no. 7828, p. 207–216, 2020.
  * [5] F. Zangeneh-Nejad, D. L. Sounas, A. Alù, and R. Fleury, “Analogue computing with metamaterials,” Nature Reviews Materials, vol. 6, no. 3, pp. 207–225, 2021.
  * [6] N. Engheta and R. W. Ziolkowski, Metamaterials: physics and engineering explorations.  Wiley-Interscience, 2006.
  * [7] A. Y. Piggott, J. Lu, K. G. Lagoudakis, J. Petykiewicz, T. M. Babinec, and J. Vucković, “Inverse design and demonstration of a compact and broadband on-chip wavelength demultiplexer,” Nature Photonics, vol. 9, no. 6, pp. 374–377, 2015.
  * [8] T. W. Hughes, M. Minkov, I. A. D. Williamson, and S. Fan, “Adjoint method and inverse design for nonlinear nanophotonic devices,” ACS Photonics, vol. 5, no. 12, p. 4781–4787, 2018.
  * [9] N. M. Estakhri, B. Edwards, and N. Engheta, “Inverse-designed metastructures that solve equations,” Science, vol. 363, no. 6433, pp. 1333–1338, 2019\.
  * [10] M. Camacho, B. Edwards, and N. Engheta, “A single inverse-designed photonic structure that performs parallel computing,” Nature Communications, vol. 12, no. 1, 2021.
  * [11] D. Colton and R. Kress, Inverse Acoustic and Electromagnetic Scattering Theory.  Springer International Publishing, 2019.
  * [12] E. M. Purcell and C. R. Pennypacker, “Scattering and Absorption of Light by Nonspherical Dielectric Grains,” The Astrophysical Journal, vol. 186, p. 705, dec 1973.
  * [13] B. T. Draine and P. J. Flatau, “Discrete-Dipole Approximation For Scattering Calculations,” Journal of the Optical Society of America A, vol. 11, no. 4, p. 1491, 1994.
  * [14] Y. M. Wang and W. C. Chew, “An iterative solution of the two-dimensional electromagnetic inverse scattering problem,” International Journal of Imaging Systems and Technology, vol. 1, pp. 100–108, jun 1989.
  * [15] P. M. van den Berg and R. E. Kleinman, “A contrast source inversion method,” Inverse Problems, vol. 13, no. 6, pp. 1607–1620, 1997.
  * [16] A. Lakhtakia, “Macroscopic theory of the coupled dipole approximation method,” Optics Communications, vol. 79, no. 1-2, pp. 1–5, 1990.
  * [17] F. M. Kahnert, Numerical methods in electromagnetic scattering theory, vol. 79-80.  2003\.
  * [18] M. A. Yurkin and A. G. Hoekstra, “The discrete dipole approximation: An overview and recent developments,” Journal of Quantitative Spectroscopy and Radiative Transfer, vol. 106, no. 1-3, pp. 558–589, 2007.
  * [19] C. A. Balanis, “Advanced engineering electromagnetics,” vol. Hoboken: J. Wiley & Sons, 2012.
  * [20] S. P. Groth, A. G. Polimeridis, and J. K. White, “Accelerating the discrete dipole approximation via circulant preconditioning,” Journal of Quantitative Spectroscopy and Radiative Transfer, vol. 240, p. 106689, 2020.
  * [21] D. Tzarouchis, M. J. Mencagli, B. Edwards, and N. Engheta, “Mathematical Operations and Equation Solving with Reconfigurable Metadevices,” arXiv: 2206.02549, pp. 1–19, May 2022.
  * [22] S. Molesky, Z. Lin, A. Y. Piggott, W. Jin, J. Vucković, and A. W. Rodriguez, “Inverse design in nanophotonics,” Nature Photonics, vol. 12, no. 11, p. 659–670, 2018.
  * [23] M. P. Bendsøe and O. Sigmund, Topology optimization: theory, methods, and applications.  Springer, 2019.
  * [24] A. Y. Piggott, J. Lu, K. G. Lagoudakis, J. Petykiewicz, T. M. Babinec, and J. Vučković, “Inverse design and demonstration of a compact and broadband on-chip wavelength demultiplexer,” Nature Photonics, vol. 9, no. 6, p. 374–377, 2015.
  * [25] A. Y. Piggott, J. Petykiewicz, L. Su, and J. Vučković, “Fabrication-constrained nanophotonic inverse design,” Scientific Reports, vol. 7, no. 1, 2017.
  * [26] O. Sigmund and D. N. Chigrin, “Topology optimization in nano-photonics,” AIP Conference Proceedings, 2009.
  * [27] J. Lu, S. Boyd, and J. Vučković, “Inverse design of a three-dimensional nanophotonic resonator,” Optics Express, vol. 19, no. 11, p. 10563, 2011\.
  * [28] C. M. Lalau-Keraly, S. Bhargava, O. D. Miller, and E. Yablonovitch, “Adjoint shape optimization applied to electromagnetic design,” Optics Express, vol. 21, no. 18, p. 21693, 2013.
  * [29] “Comsol multiphysics® v. 5.6. www.comsol.com. comsol ab, stockholm, sweden.”
  * [30] W. Zhang and A. Hoorfar, “MIMO ground penetrating radar imaging through multilayered subsurface using total variation minimization,” IEEE Transactions on Geoscience and Remote Sensing, vol. 57, no. 4, pp. 2107–2115, 2018.
  * [31] R. Autieri, M. D’Urso, T. Isernia, and V. Pascazio, “Inverse profiling via an effective linearized scattering model and mrf regularization,” IEEE Geoscience and Remote Sensing Letters, vol. 8, no. 6, pp. 1021–1025, 2011.
  * [32] W. Zhang, A. Hoorfar, and C. Thajudeen, “Two- and three-dimensional fast intrawall imaging with microwave tomographic algorithm,” International Journal of Antennas and Propagation, vol. 2018, p. 1–12, 2018.
  * [33] M. G. Amin and F. Ahmad, “Compressive sensing for through-the-wall radar imaging,” Journal of Electronic Imaging, vol. 22, no. 3, p. 030901, 2013\.
  * [34] K. M. Yemelyanov, N. Engheta, A. Hoorfar, and J. A. Mcvay, “Adaptive polarization contrast techniques for through-wall microwave imaging applications,” IEEE Transactions on Geoscience and Remote Sensing, vol. 47, no. 5, p. 1362–1374, 2009.
  * [35] E. J. Baranoski, “Through wall imaging: Historical perspective and future directions,” 2008 IEEE International Conference on Acoustics, Speech and Signal Processing, 2008.
  * [36] W. Zhang and A. Hoorfar, “Three-dimensional synthetic aperture radar imaging through multilayered walls,” IEEE Transactions on Antennas and Propagation, vol. 62, no. 1, p. 459–462, 2014.
  * [37] W. Zhang and A. Hoorfar, “Three-dimensional real-time through-the-wall radar imaging with diffraction tomographic algorithm,” IEEE Transactions on Geoscience and Remote Sensing, vol. 51, no. 7, p. 4155–4163, 2013.
  * [38] M. F. Imani, J. N. Gollub, O. Yurduseven, A. V. Diebold, M. Boyarsky, T. Fromenteze, L. Pulido-Mancera, T. Sleasman, and D. R. Smith, “Review of metasurface antennas for computational microwave imaging,” IEEE Transactions on Antennas and Propagation, vol. 68, no. 3, p. 1860–1875, 2020\.
  * [39] M. Hajebi and A. Hoorfar, “A multiscale optimization technique for large-scale subsurface profiling,” IEEE Geoscience and Remote Sensing Letters, vol. 18, no. 10, p. 1706–1710, 2021.


overlapping structure is the same as for the standard repetition code. Note
the inclusion of effective X-basis measurement and Z-basis classical feedback
at boundaries being stepped away from; these correct the distance-1 X
observable, and are unnecessary when considering the repetition code as a
classical code. By combining various cycles, the repetition code state can
step in either direction by one physical qubit per cycle.

###  B.2 CXSWAP Repetition Code

We can similarly use freedom in the shape of detecting regions to provide
cycle circuits that use CXSWAP gates instead of CX gates. This substitution
also allows both standard circuits that preserve measure and data qubit roles,
and stepping cycles where they are exchanged. This is illustrated in Figure
B.2.

Figure B.2:  Repetition code circuits with CXSWAP gates. The bit-flip
repetition code can also be implemented using CXSWAP gates. a) Four possible
cycle circuits. The Standard cycles preserve qubit roles, while the stepping
cycles exchange qubit roles. b) All detecting regions for a distance-3 step
code. In order, the cycles are [Standard Down, Standard Down, Step up,
Standard Up]. The combination of these cycles produces a variety of shapes of
detecting regions, but their size and overlapping structure is the same as for
the standard repetition code. Note the inclusion of effective X-basis
measurement and Z-basis classical feedback at boundaries being stepped away
from; these correct the distance-1 X observable, and are unnecessary when
considering the repetition code as a classical code.

##  Appendix C Sliding Surface Codes for Logical Compilation

In Section 5, we focused on explaining the walking circuit and describing its
benefits in simplifying hardware. The capability of walking logical patches
also provides new possible freedoms at a higher level, presenting a new
primitive operation at the level of lattice surgery which can permit some
logical operations to be cheaper. Here, we detail one such application, using
the sliding behaviour to perform a lateral shift on a densely packed register
of logical qubit patches.

In algorithms compiled down to lattice surgery primitives, shifting a dense
linear register of logical qubit patches is a common and desirable operation
[GE21]. A typical use case is when processing register bits iteratively,
treating the register as a queue to be consumed; here, being able to shift the
register enables the compilation to operate on the least-significant bit and
then shift the register, rather than needing to operate on each bit in its
original location.

However, register shifting is not a very natural operation in lattice surgery,
as the shifting of each logical patch must be done serially, as shown in
Figure C.1a. For a linear register of $N$ logical patches with cost distance
$d$ attempting to shift $S$ patch distances along its length, the naive total
cost of this strategy in space-time volume is $N(N+S)=N^{2}+NS$ in units of
$d^{3}$. This operation has a cost quadratic in the size of the register, and
is particularly expensive for shifting a large register over a short distance.

Comparatively, we can use the sliding circuit for $2d$ cycles to shift all the
logical patches in parallel by a distance $d$, as shown in Figure C.1b. This
requires the existence of an additional qubit row parallel to the direction of
shifting for the walking circuit to execute, but does not effect the overall
volume consumed by the register as it walks, especially if surrounding patches
are wiggling in-phase with the sliding. This can greatly reduces the cost of
some logical compilations, such as using a large register as a FIFO queue.

Figure C.1:  Cost comparison for register shifting using lattice surgery vs
using sliding. Left: A 2D projection of the naive strategy using lattice
surgery expansions and contractions to shift a logical register of $N=6$
qubits by a distance $S=1$, consuming $(N^{2}+NS)$ space-time blocks each
$d\times d\times d$ in volume. Right: the same operation achieved using the
sliding behaviour on all patches making up the register, consuming only $2NS$
space-time blocks. For shifts $S<N$, this strategy lowers the overall space-
time cost.

##  Appendix D Numerical Benchmarking and Noise Model

For numerical benchmarking of our constructions, we use Stim [Gid21]. In the
main text, we compile the relevant circuit for superconducting hardware using
CZ gates, and use the SI1000 noise model introduced in [GNM22]. This noise
model is inspired by the hardware errors experienced by superconducting
transmon qubit arrays. Here, we augment this error model with an ISWAP gate
with the same leading fidelity. This elides any different in fidelity between
performing a CZ or ISWAP gate, whereas we expect the error models in hardware
to be noticeably different.

Other constructions using CNOT or entangling measurements are benchmarked
using the simpler UniformDepolarizing noise model, which does not have ratios
of errors adjusted to minic real devices. The definitions of the noise
operations in all three models are given in Table D.1, and the models
themselves are detailed in Table D.2.

The code used to run the numerical simulations is included in our data and
code repository [MBG23]. We used the Sinter interface and took up to $1\times
10^{8}$ shots for each circuit implementation, completing early if we
collected $1\times 10^{3}$ logical errors for a specific circuit case. For
decoding, we used an internal decoder that performs correlated-error-aware
minimum weight perfect matching [Goo+22]. We combined the decoded results for
memory experiments in the X and Z basis individually into statistics for a
combined XZ basis by projecting the likelihood that neither X nor Z bases
experience a logical error. The benchmarking for all circuits collectively
took $65345126$ core-seconds in total, which is $2.07$ core-years. It executed
in a little under a 8 days on a 96 core machine.

To produce the teraquop footprint, we fit a line to the XZ logical error rate
versus code distance to extrapolate to a logical error rate of $1\times
10^{-12}$ to find the code distance and therefore footprint required, which we
plot as the teraquop footprint in the main text. Our plots generally include
maximum likelihood highlighting to indicate confidence in our numerical
estimates. These correspond to the region of hypotheses up to 1000 times less
likely than the maximum likelihood hypothesis. The procedure we used also
included in our data and code repository.

Noisy Gate | Definition
---|---
$\text{AnyClifford}_{2}(p)$ | Any two-qubit Clifford gate, followed by a two-qubit depolarizing channel of strength $p$.
$\text{AnyClifford}_{1}(p)$ | Any one-qubit Clifford gate, followed by a one-qubit depolarizing channel of strength $p$.
$\text{R}_{Z}(p)$ | Initialize the qubit as $\ket{0}$, followed by a bitflip channel of strength $p$.
$\text{R}_{X}(p)$ | Initialize the qubit as $\ket{+}$, followed by a phaseflip channel of strength $p$.
$M_{Z}(p,q)$ | Measure the qubit in the $Z$-basis, followed by a one-qubit depolarizing channel of strength $p$,
| and flip the value of the classical measurement result with probability $q$.
$M_{X}(p,q)$ | Measure the qubit in the $X$-basis, followed by a one-qubit depolarizing channel of strength $p$,
| and flip the value of the classical measurement result with probability $q$.
$M_{PP}(p,q)$ | Measure a Pauli product $PP$ on a pair of qubits,
| followed by a two-qubit depolarizing channel of strength $p$,
| and flip the classically reported measurement value with probability $q$.
$\text{Idle}(p)$ | If the qubit is not used in this time step, apply a one-qubit depolarizing channel of strength $p$.
$\text{ResonatorIdle}(p)$ | If the qubit is not measured or reset in a time step during which other qubits are
| being measured or reset, apply a one-qubit depolarizing channel of strength
$p$.

Table D.1:  Modified from [GNM22]. Noise channels and the rules used to apply them. Noisy rules stack with each other - for example, Idle($p$) and ResonatorIdle($p$) can both apply depolarizing channels in the same time step.  Name |  | Uniform
---
Depolarizing
| Superconducting
---
Inspired (SI1000)
Noisy Gateset |  | $\text{CX}(p)$
---
$\text{CXSWAP}(p)$
$\text{AnyClifford}_{1}(p)$
$\text{R}_{Z/X}(p)$
$M_{Z/X}(p,p)$
$M_{PP}(p,p)$
$\text{Idle}(p)$
| $\text{CZ}(p)$
---
$\text{ISWAP}(p)$
$\text{AnyClifford}_{1}(p/10)$
$\text{Init}_{Z}(2p)$
$M_{Z}(p,5p)$
$M_{ZZ}(p,5p)$
$\text{Idle}(p/10)$
$\text{ResonatorIdle}(2p)$
Table D.2:  Modified from [GNM22]. Details of the error models used in this
paper. Circuits compiled for CX or CXSWAP gates are benchmarked using
UniformDepolarizing. Circuits compiled for CZ or ISWAP gates are benchmarked
using SI1000. The superconducting-inspired acronym refers to an expected cycle
time of about $1000$ nanoseconds for the standard surface code circuit cycle
[Goo+21, Goo+22]. See Table D.1 for definitions of the noisy gates.

##  Appendix E Further Benchmarking

In this section, we present benchmarking for each of the constructions listed
in Table 1.

We provide teraquop footprints for each circuit in groups sharing an error
model. Figure E.1 includes all planar circuits using CZ or ISWAP gates and
benchmarked using the SI1000 error model. Figure E.2 includes all planar
circuits using CX or CXSWAP gates and benchmarked using the
UniformDepolarizing error model. Figure E.3 includes all toric circuits using
CX gates, benchmarked using the UniformDepolarizing error model. Finally,
Figure E.4 includes all circuits using entangling measurements. As detailed in
Appendix D, our noise models both treat entangling measurements as having the
same depolarizing error as two qubit gates, so circuits using entangling
measurements perform significantly better than circuits using only two qubits
gates.

In the ancillary files available with the paper, and in the data repository
[MBG23], we additionally provide set of supplementary figures. For each
circuit, we provide a visualization of the circuit schedule, a plot of the
logical error versus physical error rate and the logical error rate versus
code distance. The circuit schedules were automatically generated using Stim’s
diagram functionality, and can be reproduced from the circuits also included
in the data repository. The raw data for the benchmarking plots, and the code
that produces them, is also available in our code repository.

Figure E.1:  Teraquop footprints for planar circuits using the SI1000 error
model.  The number of physical qubits required for a single code patch to
achieve a logical error rate of $1\times 10^{-12}$ over a $d\times d\times d$
space-time block. Each curve averages memory experiments in both Z and X
basis, and uses the SI1000 noise model as described in Appendix D. Circuit
details for each curve are included in Table 1.

Figure E.2:  Teraquop footprints for planar circuits using the
UniformDepolarizing error model.  The number of physical qubits required for a
single code patch to achieve a logical error rate of $1\times 10^{-12}$ over a
$d\times d\times d$ space-time block. Each curve averages memory experiments
in both Z and X basis, and uses the UniformDepolarizing noise model as
described in Appendix D. Circuit details for each curve are included in Table
1.

Figure E.3:  Teraquop footprints for circuits benchmarked with toric boundary
conditions.  The number of physical qubits required for a single code patch to
achieve a logical error rate of $1\times 10^{-12}$ over a $d\times d\times d$
space-time block. Each curve averages memory experiments in both Z and X
basis, and uses the UniformDepolarizing noise model as described in Appendix
D. Circuit details for each curve are included in Table 1.

Figure E.4:  Teraquop footprints for circuits benchmarked with hybrid
entangling operations.  The number of physical qubits required for a single
code patch to achieve a logical error rate of $1\times 10^{-12}$ over a
$d\times d\times d$ space-time block. Each curve averages memory experiments
in both Z and X basis. 3-CZ_MZZ uses the SI1000 noise model, 3-CX_MZZ_MXX and
TORIC-3-CX_MZZ_MXX use the UniformDepolarizing noise model, as described in
Appendix D. Circuit details for each curve are included in Table 1. Note that
TORIC-3-CX_MZZ_MXX has no data point for $p=1\times 10^{-4}$; No errors were
sampled for distances above $d=4$, so extrapolation to the teraquop regime was
not possible.