# Inverse design of plasma metamaterial devices with realistic elements

Jesse A. Rodríguez<EMAIL_ADDRESS>Department of Mechanical Engineering,
Stanford University, Stanford, CA 94305  Mark A. Cappelli  Department of
Mechanical Engineering, Stanford University, Stanford, CA 94305

###### Abstract

In an expansion of a previous study [1], we apply inverse design methods to
produce two-dimensional plasma metamaterial devices with realistic plasma
elements which incorporate quartz envelopes, collisionality (loss), non-
uniform density profiles, and resistance to experimental error/perturbation.
Backpropagated finite difference frequency domain simulations are used to
design waveguides and demultiplexers operating under the transverse magnetic
polarization. Optimal devices with realistic elements are compared to previous
devices with idealized elements, and several parameter initialization schemes
for the optimization algorithm are explored. Demultiplexing and waveguiding
are demonstrated for microwave-regime devices composed of plasma elements with
reasonable space-averaged plasma frequencies $\sim 10$ GHz and a collision
frequency $\sim 1$ GHz, allowing for future in-situ training and experimental
realization of these designs.

##  I Introduction

Inverse design methodology is an algorithmic technique by which optimal
solutions are iteratively approximated from an a priori set of performance
metrics. When applied to electromagnetically-active systems [2, 3, 4, 5, 6, 7,
8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19], devices are built by
designating a design region wherein various numerical methods are used to
solve Maxwell’s equations for a given distribution of sources and material
domain. The constitutive properties (e.g., relative permittivity or
permeability) of the design region are then mapped to trainable parameters
which are modified to mold the propagation of source fields through the
designated region in a manner that fulfills the user’s performance criteria.
Typically, the performance criteria are encoded as an objective function that
needs to be maximized or minimized. In many problems, the objective functions
involve an inner product between the realized and the desired electromagnetic
fields in regions of interest. Given Maxwell’s equations as a constraint,
automatic differentiation allows the calculation of numerical gradients
relating the chosen objective function to the parameters that encode the
permittivity structure of the training region. The algorithm then iteratively
adjusts these parameters to maximize the overlap between the realized and
desired modes of propagation [2, 3, 4]. The result is a novel device geometry
specialized to the user’s performance criteria.

Unfortunately, the set of physically realizable devices for inverse design
methods is a small subset of the configuration space, leading to intensive
investigation of devices which are amenable to these constraints such as
photonic crystal-style devices [5, 8, 2, 6, 7] and metasurfaces [9, 10, 11,
12]. In general, these limitations require that certain restrictions be
applied during the training process. For example, utilizing a continuous range
of permittivities throughout the simulation domain greatly improves the
ability of the algorithm to extremize the objective function, but continuity
is not easily realizable with conventional manufacturing techniques and
materials. To resolve this issue, a map from the parameterized design region
to the set of actually manufacturable device configurations is used. These
maps encode restrictions on the range and spatial distribution of the domain
permittivity. A common technique in the literature is the use of nonlinear
projection to design less versatile, but more easily manufacturable, binarized
photonic devices [13]. Such a device can be printed or created via lithography
as a single-use design that has the requisite performance limitations of such
a constrained parameterization scheme.

Rather than restricting the range of permittivities, we focus instead on
limiting the spatial configuration by parameterizing the design region as a
plasma metamaterial (PMM). Our PMM is a periodic array of cylindrical gaseous
plasma elements (rods) which can be experimentally realized via the use of
discharge lamps or laser-generated plasmas. PMMs are ideal candidates for
inversely designed devices since the plasma density, and hence the
permittivity, of each element within the PMM can be dynamically and precisely
tuned through a continuous range of values. This yields an infinite
configuration space for training purposes [20] and allows a single device to
serve multiple functions. Past work with inverse design electromagnetic
devices achieved reconfigurability via methods such as refractive index
changes [9] in materials like liquid crystals [11]. In our case, the relative
permittivities of elements within the PMM device can be dynamically tuned by
varying the gas discharge current or laser power. In addition, when the PMM is
composed of elements that are small compared to the operating wavelength and
the source is polarized properly, we can access localized surface plasmon
(LSP) modes along the boundary between a positive permittivity background and
a negative permittivity plasma element. Such surface modes can yield more
complex and efficient device behavior [21], allowing for a particularly high
degree of reconfigurability, but for devices which seek to maximize
transmission, they can cause serious degradation of performance [1]. Prior
experimental studies of plasma photonic crystals (PPCs) have already
highlighted the richness of this geometry in waveguide [1, 22] and bandgap
[21, 23, 24] devices. The rich and varied physics of PMMs is well-matched to
inverse design methods which promise to allow more holistic and efficient
exploration of the configuration space.

The design processes for binarized and PMM geometries are contrasted in Fig.
1. Binarized geometries (Fig. 1a) are challenging to optimize because of the
requisite discontinuities. The jumps in permittivity from one material to
another are difficult to handle in gradient-based methods that require the
permittivity to vary continuously. Binarization can instead be approximated by
any number of techniques, including gaussian blurring and highly nonlinear
projections [13]. In these configurations, the trainable parameters $\rho$ (a
vector of dimension $n_{x}n_{y}$) maintain a one-to-one correspondence with
the individual pixels of the training region. In contrast, the number of
training parameters in our PMM (Fig. 1b) is determined by the dimensions of
the plasma rod array, where there is a one-to-one correspondence between each
plasma element and the $n_{x}n_{y}$ elements of $\rho$ (Fig. 1c). In the end,
both parameterization schemes result in a relative permittivity matrix
$\varepsilon_{r}$ (of size $N_{x}\times N_{y}$, dependent on the simulation
resolution) containing permittivity values at each pixel in the simulation
domain.

Figure 1: Flow chart describing the algorithmic design of either (a) binarized
material with penalized material use [25] or (b) array of plasma columns
either in the TM mode ($\mathbf{E_{x}}$) or the the TE mode
($\mathbf{E_{z}}$). Though the parameterization stages are different, both
configurations can be given identical sources and objective criteria. $J$
represents the modal source for the FDFD simulation. Examples of the training
parameter vector $\rho$, permittivity matrix $\varepsilon_{r}$ and simulated E
fields for a simple PMM device are provided in (c).

##  II Methods

Fig. 2 provides a diagram of our PMM device in a directional waveguide
configuration. A demultiplexer configuration, described later in the paper,
exhibits the exact same array of tunable plasma elements but with different
entrance/exit waveguides. The device is simply a 10$\times$10 array of tunable
plasma discharge tubes suspended in air and spaced according to the limits of
existing experimental facilites [22, 23]. A detailed description of the plasma
elements is found below.

Modal sources are introduced at the input waveguide(s) and allowed to scatter
through the training region before collecting at the desired output waveguide.
In this work, the fields are propagated through the domain via finite
difference frequency domain (FDFD) simulations computed with Ceviche, an
autograd-compliant electromagnetic simulation tool [3] that allows for
calculation of the gradients of optimization objectives with respect to the
input parameters that encode the permittivity domain. Note that the use of
FDFD implies that all computed devices represent the steady state solution
achieved after some characteristic time. The simulation domain was discretized
using a resolution of 75 pixels per lattice constant $a$ due to CPU memory and
runtime limitations, resulting in a domain size of $N_{x}=N_{y}=1500$ pixels
in each case. The ’optimal’ devices were then run at higher resolutions up to
150 pixels per lattice constant on a high-performance computing cluster to
ensure that the device functionality was preserved; ensuring that the training
resolution was high enough. A perfectly-matched boundary layer (PML) $2a$ in
width was applied along the domain boundaries.

The polarization of the input source has a strong effect in devices of this
nature, either $E_{z}$ ($E$ out of the page), which we call the TM case, or
$E_{y}$ ($H$ out of the page), which we call the TE case. Previous work with
PMMs indicates that both the TM and TE responses are tunable, with the latter
benefiting (and suffering, depending on the desired device functionality) from
the presence of LSP modes [26, 21], while the former makes more direct use of
dispersive and refractive effects [21, 22, 23]. In our prior work [1], we see
that for devices which seek to preserve an input signal like the waveguides
and demultiplexers we present below, the TE polarization leads to poor
performance, so we choose to focus on the TM mode alone in this study.

The simulated $\mathbf{E}$ fields are masked to compute the field intensity
and mode overlap integrals along planes of interest within the problem
geometry. These integrals are used to calculate the objective
$\mathcal{L}(\rho)$. Our simulation tool, Ceviche, is then used to compute
numerical gradients of the objective with respect to the training parameters
via forward-mode differentiation [3]. Ceviche uses the Adam optimization
algorithm [27] (gradient ascent, in essence) to iteratively adjust $\rho$ and
thereby maximize $\mathcal{L}$. Optimization was conducted with learning rates
ranging from $0.001-0.01$, and the default Adam hyperparameters
$\beta_{1}=0.9$, and $\beta_{2}=0.999$.

Figure 2: Schematic diagram for the PMM waveguide device consisting of a
$10\times 10$ array of plasma discharge bulbs suspended in air with metal
waveguides functioning as the entrance and exit(s). The bulb size is chosen to
be commensurate with current experimental facilities. The ’probe’ slices
represent where the mode overlap and intensity integrals that constitute the
objective function are calculated. $\varepsilon_{bg}$ is the background
permittivity. $\tilde{\omega}$ is the nondimensionalized operating frequency
$\omega/(2\pi c/a)$.

In summary, our optimization problem is

$\displaystyle\max_{\rho}\quad\mathcal{L}(\mathbf{E})$ $\displaystyle
given\quad\nabla\times\frac{1}{\mu_{0}}\nabla\times\mathbf{E}-\omega^{2}\varepsilon(\omega,\rho)\mathbf{E}=-i\omega
J.$

where $\mathcal{L}$ is the objective composed of a set of mode overlap
integrals ($L_{2}$ inner product of simulated field with desired propagation
mode), $\rho$ is an $n_{x}n_{y}$-dimensional vector that contains the
parameters that set the electron density of each of the PMM elements,
$\mu_{0}$ is the vacuum permeability, E is the electric field, $J$ is a
current density used to define a fundamental modal source at the input
waveguide, $\omega$ is the field frequency, and $\varepsilon(\omega,\rho)$ is
the spatially-dependent permittivity that is encoded by $\rho$. In practice,
this permittivity distribution is achieved by varying the plasma density
through control of the discharge current in each of the PMM elements according
to the Drude model, which is where the dependence on $\omega$ arises. The non-
dimensionalized plasma frequencies in $\rho$ are mapped to element
permittivities via the Drude model;

$\displaystyle\varepsilon=1-\frac{\omega_{p}^{2}}{\omega^{2}+i\omega\gamma}$

where $\gamma$ is the collision/damping rate,
$\omega_{p}^{2}=\frac{n_{e}e^{2}}{\varepsilon_{0}m_{e}}$ is the plasma
frequency squared, $n_{e}$ is the electron density, $e$ is the electron
charge, $m_{e}$ is the electron mass, and $\varepsilon_{0}$ is the free-space
permittivity.

The devices that follow have been explored previously in ref. [1], but mostly
as a proof of concept with very simplistic plasma elements that were
collisionless, uniform-density rods with no quartz envelope. In this study, we
attempt to model all major non-ideal factors in order to assess whether or not
such devices will be readily experimentally realized. Of course, real gas
discharge plasmas will be affected by a range of non-ideal factors in addition
to the damping rate included in the Drude model above. Such damping can be
brought to very low levels in practice by reducing the discharge pressure, but
this can come at the expense of plasma density (and correspondingly, plasma
frequency), which affects the range of values of the plasma dielectric
constant that we can access. Furthermore, the plasma density profile within
the gaseous elements is not uniform and dependent on the methods used to
generate the plasma. For example, within a discharge tube at nominal operation
conditions, the plasma has a core of relatively uniform density plasma and
then a sheath region near the bulb where the density drops to negligible
values at the inner bulb surface. In the past, our group has approximated this
with a uniform plasma column of radius $r_{i}/\sqrt{2}$ where $r_{i}$ is the
inner radius of the bulb’s quartz envelope. This leads to decent agreement
with experiment in photonic crystal configurations [22, 23, 26, 21], but may
not be appropriate for simulating devices with EM wave propagation structures
as complicated as are found in our inversely designed PMMs. To best account
for the non-uniform density profile, we choose as our model a 6th-order
polynomial as a function of radius, where the trainable parameters $rho$ map
to the space-averaged plasma frequency of each of the rods. This profile has
been shown to lead to the closest match between simulated and observed
transmission spectra for PPC devices [28] (when compared to other polynomial
profiles), and closely matches the qualitative behavior we expect from a wall-
bounded plasma where the density profile is mostly driven by diffusion. Future
work will further refine the density profile model via either kinetic or
hydrodynamic models. In addition, we include the $\varepsilon=3.8$ quartz bulb
around the plasma elements that was omitted in ref. [1]. Figure 3 illustrates
the difference between the simplistic plasma elements in ref. [1] and our more
realistic ones.

Figure 3: Illustration comparing our model (labeled ’new’) and the idealized
model (labeled ’old’) from ref. [1] for the plasma discharge tubes. The
idealized model does not include the quartz bulb.

In addition to all of this, while an accurate and precise control of element
permittivity by varying discharge current is possible in theory, in practice a
number of factors (such as thermal variation, mechanical/electrical error,
etc.) can lead to small unknown perturbations in discharge current $\sim 1$ \-
$5$%, which leads to density perturbations of about the same order of
magnitude relative to the nominal density. Thus, if one wishes to translate
these designs to a physical device, such perturbations must be accounted for.
This can be done in one of three ways: 1) random noise can be programmed into
the learning algorithm itself at greater computational cost as training would
take longer, 2) the device can be trained in-situ which would automatically
take experimental noise within the apparatus into account, or 3) one can
optimize the device without taking any noise into account and then do a
sensitivity test on the optimized parameters to determine how resistant the
resultant designs are to perturbations in plasma density. In this study, we
will be carrying out the third option.

The manner in which we conduct the sensitivity test is as follows: Since the
elements of the device would be individually controlled in practice, we will
add mean-zero independently distributed Gaussian noise to the density of each
element with a standard deviation $\sigma=p\sqrt{n_{e}}$ where $p$ represents
a ’perturbation factor’ that sets the variance of the random noise. If
$p=0.05$, we can expect variation in plasma density amongst the columns on the
order of $\sim 5$% of the plasma density in each individual column. This is
the same as reassigning the plasma density of each element according to

$\displaystyle
n_{e,perturbed}=\left|n_{e,opt}+\mathcal{N}(0,p^{2})n_{e,opt}\right|,$

where the absolute value is applied so we don’t end up with non-physical
negative densities when the Gaussian values are spuriously large and negative.
For each device, we run 10 different random perturbations with the same value
for $p$ and increase $p$ in increments of 0.01 until we see a significant
failure mode in at least 3 of the 10 runs. For a failure mode to be deemed
’significant’, we require a qualitatively obvious degradation of the objective
functionality where we either have much less signal in the correct output
waveguide or a large amount of loss to the surroundings/incorrect exit
waveguides, or both. The dominant and/or typical failure mode amongst the 3 or
more dysfunctional perturbed devices will be presented below alongside the
original design to illustrate how robust each device is and how the device
functionality will deteriorate if the plasma density is perturbed to a high
enough degree.

##  III Results and Discussion

###  III.1 Directional waveguide

Electromagnetic waveguides are fairly common in inverse design devices either
as the primary function or as important building blocks to more complicated
devices [6, 19]. Presented in Fig. 4 are field simulations and initial/final
domain permittivities of straight waveguide devices designed for the TM
polarization. Operating at $\omega=0.25\times(2\pi c/a)$, or a non-
dimensionalized source frequency $\tilde{\omega}=\omega/(2\pi c/a)=0.25$, this
represents the simplest type of optimization problem for our PMM
configuration. Because of the simplicity of this geometry, we limit the plasma
frequency to $f_{p}=7$ GHz, a modest operating condition that could be easily
reached and run at long time scales without damaging the discharge bulbs. In
this and all other simulations, the relative permittivity of the input and
output waveguide walls was set to $\varepsilon=-1000$ to serve as a lossless
metal. The waveguide objective is simply the $L_{2}$ inner product of the
simulated field with an $m=1$ propagation mode (in order to preserve the input
mode) at the location of the probe in the desired exit waveguide minus the
integrated field intensity at the incorrect exit to discourage leakage into
the wrong output:

$\mathcal{L}_{wvg}=\int E\cdot E_{m=1}^{*}dl_{\text{desired
exit}}-\int|E|^{2}dl_{\text{incorrect exit}}.$

One could likely improve performance even further by discouraging all leakage;
i.e. subtracting the field intensity integrated over all space other than
within the correct output, but for the devices in this study the improvement
is marginal.

In the prior idealized study [1], the initial parameters were set to be
uniform among all the plasma elements since the initial conditions were found
to be unimportant; i.e. the optimal parameters were very similar between runs
with different intial parameters. Upon adding all the non-ideal factors
described above, the optimization objective function became significantly more
complicated with many more local minima than before, and so the initial
conditions became very closely linked with the device performance once the
algorithm converged. In fig. 4 we can see that when initializing the domain
with all the plasmas turned off (fig. 4(b)), the functionality of the device
suffers, whereas if the device is initialized with the optimal parameters from
ref. [1] (fig. 4(a)), the final device performs much better even though the
parameters change significantly during training. The trained domain for the
naively intialized device is also significantly less explainable; that is,
given a plot of the relative permittivity, one would not be able to easily
guess where the source wave would be directed, while for the other case is it
easier to see that the wave would propagate directly through the device.

Figure 4: (a) Initial and trained relative permittivity domains and FDFD
simulated field magnitude $|\mathbf{E_{z}}|$ ($\hat{\mathbf{z}}$ out of the
page) showing the resultant performance for the straight waveguide initialized
with the optimal parameters from [1], where the maximum plasma frequency among
the plasma elements is about $\tilde{\omega}_{p}=0.35$ which translates to
$f_{p}=7$ GHz. (b) Initial and trained relative permittivity domains and field
magnitude for the device initialized with all plasma elements deactivated,
where the maximum plasma frequency among the plasma elements is likewise about
$\tilde{\omega}_{p}=0.35$.

Next, we defined an objective intended to produce a waveguide with a 90-degree
bend. The algorithm was run almost identically to the straight waveguide case,
the only difference being that the objective overlap integral was computed at
the bottom exit waveguide. The results for this case are found in Fig. 5.
Again we see that, while both cases result in decent performance, the device
intialized with the optimal parameters from ref. [1] performs better with
significantly less loss to the surroundings.

Figure 5: (a) Initial and trained relative permittivity domains and FDFD
simulated field magnitude $|\mathbf{E_{z}}|$ ($\hat{\mathbf{z}}$ out of the
page) showing the resultant performance for the bent waveguide initialized
with the optimal parameters from [1], where the maximum plasma frequency among
the plasma elements is about $\tilde{\omega}_{p}=0.35$ which translates to
$f_{p}=7$ GHz. (b) Initial and trained relative permittivity domains and field
magnitude for the device initialized with all plasma elements deactivated,
where the maximum plasma frequency among the plasma elements is likewise about
$\tilde{\omega}_{p}=0.35$.

The devices presented in Figs. 4 and 5 appear to make expert use of the
underlying electromagnetics. Inverse design allows a nuanced search of the
configuration space, yielding more effective designs that preserve the input
mode. We emphasize here that the final training results, even those of the
well-initialized devices, are not intuitive configurations. A human-reasoned
design for this device may more closely resemble that of Wang et al. in ref.
[22], where the elements along the desired path of propagation are made
equivalent to the background permittivity (i.e., the plasma is turned off) and
the remainder of the rods are activated at the lowest permittivity (highest
plasma frequency) possible to approximate a metallic waveguide. Instead,
similar to the previous study, the algorithm settled on a configuration where
a high contrast in permittivity of neighboring plasma columns is only present
near the edges of the device to prevent leaking. The remaining elements in the
devices are given a seemingly random distribution of plasma frequencies,
keeping in line with typical results of inverse design schemes where strange
structures often arise [13, 16]. In the end, the performance achieved for the
directional waveguide configurations with all the non-ideal factors present is
very similar to that achieved in ref. [1], with only slightly more loss to the
surroundings, and of course, to the plasma elements themselves.

In fig. 6 we see the results of the sensitivity analysis for the best straight
and bent waveguide configurations. In both cases, but particularly in the
straight waveguide case, we see that the optimal design is very robust to
random perturbations of the device parameters as $p$ needed to be increased to
$\geq 0.39$ to see an appreciable loss in functionality. In the straight
waveguide case, it appears that we have a combination of loss within the
columns and reflection into the input waveguide causing a decrease in the
relative field intensity in the output waveguide, where in the bent waveguide
case, the device tends to fail by simply leaking much more of the input signal
to the surroundings. Keeping experimental error within $\sim 39\%$ is well
within current capabilities and so we can assume that this particular non-
ideal factor will not hinder experimental realization of these designs.

Figure 6: Result for the waveguide stress test, showing a significant increase
in leakage into thesurroundings when the perturbation parameter $p$ is
increased to 0.45 in for the straight waveguide and 0.39 for the bent
waveguide.

The results depicted in Figs. 4 and 5 are encouraging. Unlike other inversely
designed devices which are static, the PMM geometry allows a single device to
achieve both the straight and bent propagation modes because the permittivity
of the device elements can be tuned electronically. Quickly switching between
the straight and bent waveguide behavior, which is easily enabled by the
controllable plasma current, opens up the possibility of transistor-like
optical switches that are experimentally feasible.

###  III.2 Demultiplexer

Next, we present a frequency demultiplexer designed to distinguish between
frequencies of $\tilde{\omega}_{1}=0.25$ and $\tilde{\omega}_{2}=0.27$.
Inverse design methods have been used before to optimize wavelength
demultiplexers in optical photonic structures [13, 10], and the demultiplexer
is a key component in many optical computing schemes [29]. The optimization
problem here is fundamentally different from the waveguiding case for two
reasons: (1) a single device must now adhere to an objective based on two
simulations, one for each operating frequency (2) the permittivity
distribution is slightly different for the two frequencies because of the
Drude dispersion relation. With this considered, the demultiplexer objective
is:

$\mathcal{L}_{mp}=\left(\int E_{\omega_{1}}\cdot
E_{m=1}^{*}dl_{\omega_{1}\text{ exit}}\right)\left(\int E_{\omega_{2}}\cdot
E_{m=1}^{*}dl_{\omega_{2}\text{
exit}}\right)-\left(\int|E_{\omega_{1}}|^{2}dl_{\omega_{2}\text{
exit}}\right)\left(\int|E_{\omega_{2}}|^{2}dl_{\omega_{1}\text{
exit}}\right),$

where $E_{\omega_{1}}$ is the simulated field for the $\tilde{\omega}=0.25$
source and $E_{\omega_{2}}$ is the simulated field for the
$\tilde{\omega}=0.27$ source. Since this objective is significantly more
complex than the waveguide case, we allow the plasma frequency of the rods to
increase up to $f_{p}=10$ GHz, an operating condition which can be maintained
for time scales $\leq 10$ seconds at a time without damaging the discharge
bulbs. If operated in short ($\sim 1$ $\mu$s) pulses with a duty cycle that
limits the proportion of time that the bulbs are active, this higher plasma
density will not damage the discharges.

As an additional result of the significant jump in complexity from the
waveguide case, it was not sufficient to simply take the optimal parameters
from the idealized study [1] and take those as our initial parameters.
Instead, we borrowed from this idea and trained the device parameters in
several steps, starting with all the plasmas off and training where the only
non-ideal factor is the presence of the bulbs, which has an objective function
with fewer local minima, allowing the algorithm to converge quickly on a
configuration with strong performance. Then, once an optimal design was
obtained for this simpler case, we introduced the non-uniform density profile
and trained using the optimal parameters from the previous run. Finally, we
added loss to the rods, trained once again and arrived at a final trained
domain not very dissimilar from the previous step with no collisions. The
parameter evolution along with the device performance is shown below in figure
7. While this design still exhibits leakage into the incorrect ports, the
field intensity at the correct output far outstrips the incorrect signal, and
the device performance is very similar to that achieved in ref. [1].

Figure 7: Evolution of the relative permittivity of demultiplexer simulation
domain throughout the training stages as the plasma model complexity is
increased, along with TM FDFD simulations showing final device performance at
each operating frequency. The maximum plasma frequency among the elements in
this case is about $\tilde{\omega}_{p}=0.6$ which translates to $f_{p}=10$
GHz.

To conclude our analysis of the demultiplexer, we ran our stress test
procedure on the final design and found that in order to see an appreciable
level of loss to the surroundings/incorrect exit waveguide, the perturbation
parameter $p$ had to be increased all the way to 0.34, which is very close to
the values found in the waveguide sensitivity test above. The optimized and
perturbed domain along with the degraded performance can be found below in
figure 8. This indicates that the plasma density distribution that we arrived
at in this case is very robust to experimental noise/error.

Figure 8: Result for the demultiplexer stress test, showing a significant
increase in leakage into the incorrect output waveguide at
$\tilde{\omega}=0.25$ when the perturbation parameter $p$ is increased to
0.34.

The results shown here are very compelling. In all cases, the field intensity
at the correct output waveguide is substantially higher than at the incorrect
output, and the designs are extremely robust to random perturbation of the
parameters. In addition, the insight gained from the optimization procedure
for the demultiplexer objective in particular paves the way for both
experimental realization of these and more complex devices, such as the logic
gates presented in ref. [1]. If one wants to achieve a complicated device
function, we have shown here that a good way of doing so is to optimize in
stages while adding one non-ideal factor at a time, and then perhaps perform
in-situ training of the device using the optimal parameters of the most
realistic optimized device as the initial parameters. In the end, these
results suggest that the outlook for experimental realization is very bright.

##  IV Conclusion

There are several opportunities for future work based on these results. The
waveguide and demultiplexer devices designed here can be readily implemented
and compared against simulated functionality using existing 2D PPC devices.
Those parameters which are here taken as static domain parameters
$(r/a,\varepsilon_{rod})$ can themselves be used as training parameters as
well. The lattice constant in particular has been noted to have a large impact
on the trade-off between transmission to the output waveguide and control of
propagation in the design region. Though in this study we chose the lattice
constant to be commensurate with existing experimental equipment, more
tightly-packed rods would likely lead to better performance, and this can be
accomplished by using a different bulb design. Consider also the prospect of a
hexagonal or a 3D PMM arrangement. The non-dimensionalized operation frequency
($\omega$ in terms of $c/a$) may also prove to be an important optimization
parameter as the dispersive physics which are the operational basis for these
devices are highly frequency-dependent. In this study, we simply chose the
lowest frequency possible to reduce the plasma density requirements. Of
particular interest in the future are the effects of magnetization on the
landscape of possible PMM devices. When used in conjunction with circularly
polarized light, gaseous plasmas are capable of achieving $\varepsilon_{r}>1$,
opening the door to even more complex logical operations. While several tools
exist for the simulation and inverse design of photonic devices [30], more
development is required on the application of inverse design to magnetized
PMMs to incorporate such physics. Ceviche also includes a finite difference
time domain solver that can be used to compute numerical gradients for
optimization. Operating in the time domain would allow for a number of
interesting designs to be pursued, such as analog recurrent neural networks
which take advantage of nonlinearities in the plasma permittivity. The plasma
column model could also be further improved via the use of kinetic particle-
in-cell simulations for each bulb given an operating voltage, current, and
fill pressure/composition. The particle density functions obtained from these
codes could be used to provide a more accurate density profile than the
polynomial profile used in this study, and more detailed collision modeling
can be incorporated. Finally, should experimental realization directly from
simulation results prove difficult (we consider this to be a possible
outcome), the device training could actually be carried out in-situ, with
plasma element currents adjusted while the objective (transmission at a
certain point, for example) is determined in real-time. In effect, the Ceviche
forward-mode differentiation scheme that we use here could be carried out in
the laboratory with real field measurements instead of simulated fields. This
method would automatically take all non-ideal factors into account.

In conclusion, we consider the application of inverse design methods to the
creation of highly optimized and realistic 2D PMM devices – including
waveguides and a demultiplexer – for the TM propagation mode. The optimization
algorithm is shown to take advantage of the rich physics inherent to PMMs
including tunable reflection and refraction in both the epsilon-near-zero
regime and the negative epsilon regime. In both the waveguiding and
demultiplexer cases, simulated devices are shown to perform well with various
non-ideal factors such as the presence of quartz envelopes,
collisionality/loss, non-uniform density profiles, and experimental
error/noise. The optimal parameters of each of the devices are experimentally
realizable with few necessary changes to existing experimental infrastructure.
However, further work is needed, both in simulations and experimental
development, to enable device development for more complex objectives like
boolean logic [1] as well as other devices with complex wave propagation
structures.

###### Acknowledgements.

This research is partially supported by the Air Force Office of Scientific
Research through a Multi-University Research Initiative (MURI) with Dr. Mitat
Birkan as Program Manager. J.A.R. acknowledges support by the U.S. Department
of Energy, Office of Science, Office of Advanced Scientific Computing
Research, Department of Energy Computational Science Graduate Fellowship under
Award Number [DE-SC0019323].

## References

  * Rodríguez  _et al._ [2021] J. A. Rodríguez, A. I. Abdalla, B. Wang, B. Lou, S. Fan, and M. A. Cappelli, Inverse design of plasma metamaterial devices for optical computing, Phys. Rev. Applied 16, 014023 (2021).
  * Minkov  _et al._ [2020] M. Minkov, I. A. Williamson, L. C. Andreani, D. Gerace, B. Lou, A. Y. Song, T. W. Hughes, and S. Fan, Inverse design of photonic crystals through automatic differentiation, ACS Photonics 7, 1729 (2020).
  * Hughes  _et al._ [2019] T. W. Hughes, I. A. D. Williamson, M. Minkov, and S. Fan, Forward-Mode Differentiation of Maxwell’s Equations, ACS Photonics 6, 3010 (2019).
  * Su  _et al._ [2020] L. Su, D. Vercruysse, J. Skarda, N. V. Sapra, J. A. Petykiewicz, and J. Vučković, Nanophotonic inverse design with SPINS: Software architecture and practical considerations, Applied Physics Reviews 7, 11407 (2020).
  * Burger  _et al._ [2004] M. Burger, S. J. Osher, and E. Yablonovitch, Inverse problem techniques for the design of photonic crystals, IEICE transactions on electronics 87, 258 (2004).
  * Borel  _et al._ [2004] P. I. Borel, A. Harpøth, L. H. Frandsen, M. Kristensen, P. Shi, J. S. Jensen, and O. Sigmund, Topology optimization and fabrication of photonic crystal structures, Opt. Express 12, 1996 (2004).
  * Lin  _et al._ [2016] Z. Lin, A. Pick, M. Lončar, and A. W. Rodriguez, Enhanced spontaneous emission at third-order dirac exceptional points in inverse-designed photonic crystals, Phys. Rev. Lett. 117, 107402 (2016).
  * Miller [2013] O. D. Miller, Photonic design: From fundamental solar cell physics to computational inverse design, arXiv preprint arXiv:1308.0212 (2013).
  * Christiansen  _et al._ [2020] R. E. Christiansen, Z. Lin, C. Roques-Carmes, Y. Salamin, S. E. Kooi, J. D. Joannopoulos, M. Soljačić, and S. G. Johnson, Fullwave maxwell inverse design of axisymmetric, tunable, and multi-scale multi-wavelength metalenses, Opt. Express 28, 33854 (2020).
  * Pestourie  _et al._ [2018] R. Pestourie, C. Pérez-Arancibia, Z. Lin, W. Shin, F. Capasso, and S. G. Johnson, Inverse design of large-area metasurfaces, Opt. Express 26, 33732 (2018).
  * Chung and Miller [2020] H. Chung and O. D. Miller, Tunable metasurface inverse design for 80% switching efficiencies and 144° angular deflection, ACS Photonics 7, 2236 (2020), https://doi.org/10.1021/acsphotonics.0c00787 .
  * Meem  _et al._ [2020] M. Meem, S. Banerji, C. Pies, T. Oberbiermann, A. Majumder, B. Sensale-Rodriguez, and R. Menon, Large-area, high-numerical-aperture multi-level diffractive lens via inverse design, Optica 7, 252 (2020).
  * Piggott  _et al._ [2015] 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 9, 374 (2015).
  * Liu  _et al._ [2013] D. Liu, L. H. Gabrielli, M. Lipson, and S. G. Johnson, Transformation inverse design, Optics express 21, 14223 (2013).
  * Hughes  _et al._ [2018] T. W. Hughes, M. Minkov, I. A. Williamson, and S. Fan, Adjoint method and inverse design for nonlinear nanophotonic devices, ACS Photonics 5, 4781 (2018).
  * Molesky  _et al._ [2018] S. Molesky, Z. Lin, A. Y. Piggott, W. Jin, J. Vucković, and A. W. Rodriguez, Inverse design in nanophotonics, Nature Photonics 12, 659 (2018).
  * Christiansen and Sigmund [2021] R. E. Christiansen and O. Sigmund, Inverse design in photonics by topology optimization: tutorial, JOSA B 38, 496 (2021).
  * Jiang  _et al._ [2020a] J. Jiang, M.-K. Chen, and J. A. Fan, Deep neural networks for the evaluation and design of photonic devices, ArXiv abs/2007.0 (2020a).
  * Andrade  _et al._ [2019] N. M. Andrade, S. Hooten, S. A. Fortuna, K. Han, E. Yablonovitch, and M. C. Wu, Inverse design optimization for efficient coupling of an electrically injected optical antenna-led to a single-mode waveguide, Opt. Express 27, 19802 (2019).
  * Iwai  _et al._ [2020] A. Iwai, F. Righetti, B. Wang, O. Sakai, and M. A. Cappelli, A tunable double negative device consisting of a plasma array and a negative-permeability metamaterial, Physics of Plasmas 27, 23511 (2020).
  * Wang  _et al._ [2019] B. Wang, J. A. Rodríguez, and M. A. Cappelli, 3D woodpile structure tunable plasma photonic crystal, Plasma Sources Science and Technology 28, 02LT01 (2019).
  * Wang and Cappelli [2016a] B. Wang and M. A. Cappelli, Waveguiding and bending modes in a plasma photonic crystal bandgap device, AIP Advances 6, 65015 (2016a).
  * Wang and Cappelli [2016b] B. Wang and M. A. Cappelli, A plasma photonic crystal bandgap device, Applied Physics Letters 108, 161101 (2016b).
  * Wang  _et al._ [2021] B. Wang, J. A. Rodríguez, O. Miller, and M. A. Cappelli, Reconfigurable plasma-dielectric hybrid photonic crystal as a platform for electromagnetic wave manipulation and computing, Physics of Plasmas 28, 043502 (2021), https://doi.org/10.1063/5.0043336 .
  * Williamson and Minkov [2019] I. Williamson and M. Minkov, Optical inverse design with automatic differentiation, https://github.com/fancompute/workshop-invdesign (2019).
  * Righetti  _et al._ [2018] F. Righetti, B. Wang, and M. Cappelli, Enhanced attenuation due to lattice resonances in a two-dimensional plasma photonic crystal, Physics of Plasmas 25, 124502 (2018).
  * Kingma and Ba [2014] D. P. Kingma and J. Ba, Adam: A method for stochastic optimization (2014), arXiv:1412.6980 [cs.LG] .
  * Thomas [2019] W. Thomas,  _Development and validation of computational tools for the characterization and analysis of non-linear plasma photonic crystals_ , Master’s thesis, University of Washington ResearchWorks Archive (2019).
  * Xu  _et al._ [2021] X. Xu, M. Tan, B. Corcoran, J. Wu, A. Boes, T. G. Nguyen, S. T. Chu, B. E. Little, D. G. Hicks, R. Morandotti, A. Mitchell, and D. J. Moss, 11 tops photonic convolutional accelerator for optical neural networks, Nature 589, 44 (2021).
  * Jiang  _et al._ [2020b] J. Jiang, R. Lupoiu, E. W. Wang, D. Sell, J. P. Hugonin, P. Lalanne, and J. A. Fan, MetaNet: a new paradigm for data sharing in photonics research, Opt. Express 28, 13670 (2020b).


theoremTheorem[section] lemma[theorem]Lemma corollary[theorem]Corollary
University of California, Davis, USA and https://web.cs.ucdavis.edu/~doty/
<EMAIL_ADDRESS>https://orcid.org/0000-0002-3922-172X University of
California, Davis, USA<EMAIL_ADDRESS>David Doty, Aaron Ong [100]Theory of
computation Models of computation simulator: https://github.com/UC-Davis-
molecular-computing/simd-dna

###### Acknowledgements.

The authors were supported by NSF grants 1900931 and 1844976.

# Simulating 3-symbol Turing machines with SIMD||DNA

David Doty  Aaron Ong

###### Abstract

SIMD||DNA [wang2019simd] is a model of DNA strand displacement allowing
parallel in-memory computation on DNA storage. We show how to simulate an
arbitrary 3-symbol space-bounded Turing machine with a SIMD||DNA program,
giving a more direct and efficient route to general-purpose information
manipulation on DNA storage than the Rule 110 simulation of Wang, Chalk, and
Soloveichik [wang2019simd]. We also develop software [simd_simulator] that can
simulate SIMD||DNA programs and produce SVG figures.

###### keywords:

DNA storage, strand displacement, parallel computation

###### category:

##  1 Introduction

DNA storage typically encodes information in the choice of DNA sequences
[church2012next, organick2018random, bornholt2016dna], so that reading and
writing require expensive sequencing (reading DNA) and synthesis (writing DNA)
steps. An alternative “nicked storage” scheme of Tabatabaei et al.
[tabatabaei2020dna] uses a single long strand called a _register_ , with a
fixed sequence. Information is stored in the choice of short complementary
strands to bind to the register. This gives the potential to process the
stored information using _DNA strand displacement_ (see Figure 1), which
reconfigures which DNA strands are bound, without changing their sequences.
Thus manipulation of the stored information (i.e., computation) can
potentially be done _in vitro_ with simpler lab steps than DNA sequencing or
synthesis.

Figure 1: DNA strand displacement (see [seelig2006enzyme] for more details).
An input DNA strand (turquoise) binds to the short toehold region of a
complementary strand (pink) and displaces the output strand (purple). The
toehold region is so-called because, although too short to bind stably, it
allows temporary binding of the input, giving it a “foot in the door” to begin
the displacement process.

The SIMD||DNA model of Wang, Chalk, and Soloveichik [wang2019simd] is an
abstract model of such a system. It allows parallel in-memory computation on
several copies of the register; each register may store different data. In the
experimental implementation, each register strand is attached to a magnetic
bead, enabling _elution_ : washing away strands not bound to a register, while
keeping the registers (and their bound strands) in the solution due to their
attachment to the bead. This motivates the “multi-stage” SIMD||DNA model of
DNA strand displacement, which at a high level works as follows. Each stage is
called an _instruction_ , consisting of a set of strands to add to the
solution. It is assumed that strand displacement reactions proceed until the
solution reaches equilibrium, at which point all strands and complexes not
attached to a register are washed away. The strands for the next instruction
are then added. A key aspect of the model is that the wash step can constrain
what strand displacement reactions are possible afterward, compared to “one-
pot” strand displacement schemes that mix all strands from the start. This
gives the SIMD||DNA model potentially more power than one-pot DNA strand
displacement. Wang, Chalk, and Soloveichik [wang2019simd] showed SIMD||DNA
programs for binary counting and simulating cellular automata Rule 110, and
Chen, Solanki, and Riedel [chen2021parallel] showed SIMD||DNA programs for
sorting, shifting and searching in parallel. See Section 2 for a formal
definition and [wang2019simd] for more details and motivation for the model.

A major theoretical result of [wang2019simd] is a SIMD||DNA program that
simulates a space-bounded version of cellular automata Rule 110. When space is
unbounded, Rule 110 is known to be efficiently Turing universal, i.e., able to
simulate any single-tape Turing machine [cook2004universality] with only a
polynomial-time slowdown [neary2006p], though by an awkward indirect
construction and encoding with very large constant factors. We show how to
simulate an arbitrary 3-symbol space-bounded single-tape Turing machine
_directly_ with a SIMD||DNA program. Since custom manipulation of bits is much
easier to program in a Turing machine than Rule 110, this gives a more direct,
efficient, and conceptually simple method of general-purpose information
processing on nicked DNA storage. Although we have not worked out the details,
it seems likely that the construction can be extended straightforwardly to
Turing machines with alphabet sizes larger than 3. However, it is
straightforward to simulate a larger-alphabet Turing machine $M$ with a
3-symbol Turing machine $S$, for example representing each of 16 non-blank
symbols of $M$ by 4 consecutive bits of $S$.

Our construction was designed and tested using software we developed
[simd_simulator] for simulating the SIMD||DNA model. It is able to take a
description of an arbitrary SIMD||DNA program: a list of instructions, where
each instruction is a set of DNA strands to add. It produces figures
indicating visually how the steps work, both with text printed on the command
line (for quickly testing ideas) and SVG figures, such as most of those in
this paper.

##  2 Model

In this section we define the model of SIMD||DNA [wang2019simd].

Figure 2:  Notational conventions and reactions in the SIMD||DNA model. The
register strand is on the bottom in each subfigure, with a yellow round
“magnetic bead” depicted on the left (not depicted in subsequent figures).
Lightly shaded gray or pink regions denote bonds (double-stranded regions),
but later figures omit this and simply draw a forward strand (one with
$5^{\prime}$ end on left and $3^{\prime}$ end on right) immediately above the
domains to which it is bound on the register strand. (a) Conventions for
domain names of strands. Domains are numbered $1,\ldots,d$ within each cell;
$d=6$ in Figure 2(a) and $d=18$ in subsequent figures. The register strand has
the starred versions of these domains. If a top strand is horizontal over
domain $i$, it has domain $i$. If it is diagonal over the whole domain, it has
an unlabelled domain distinct from all register domains (used as a toehold
overhang for detachment, see subfigure (c)). If two strands both partially
cover a domain then they both have that domain. (b) A forward instruction
strand can attach if at least two complementary consecutive domains are
unbound on the register. (c) Reverse instruction strands can bind to toehold
overhangs on forward bound strands to detach them from the register. The fact
that the (unlabelled) toeholds are complementary is indicated by a diagonal
bend in the reverse strand matching. (d) Forward instruction strands can do
toehold-mediated strand displacement, one-way if the displacing strand
contains all the domains of the displaced strand. (e) If the displacing strand
is missing the last domain of the displaced, displacement can also happen,
known as _toehold exchange_. This is often called “reversible” since it
conserves the number of bound domains, but in the SIMD||DNA model, instruction
strands are added in large excess over registers, making it effectively
irreversible due to the entropic bias toward binding the instruction strand.
Thus it is depicted with irreversible arrows in the figure. (f) Two forward
strands can cooperate to displace a single bound top strand, even if neither
has enough domains to displace on its own.

See Figure 2 for notational conventions in the SIMD||DNA model and an
explanation of the basic strand displacement reactions. The register strand is
on the bottom in each sub-figure, with a yellow round “magnetic bead” depicted
on the left (bead not depicted in subsequent figures). A DNA strand has an
orientation, with one end called the $5^{\prime}$ end and the other called the
$3^{\prime}$ end; by convention strands are drawn as arrows with the arrowhead
on the $3^{\prime}$ end. The register strand has its $3^{\prime}$ end on the
left and $5^{\prime}$ end on the right. The model allows multiple registers to
be present in solution at once, each possibly configured differently. However,
it is assumed that register strands are sufficiently dilute that they do not
interact with each other or with strands that have been displaced from other
registers. Thus all figures depict only a single register and instruction
strands that interact with it.

The register strand is divided into _cells_ , which are further divided into
_domains_. Each domain can be thought of as a fixed-length DNA sequence with
relatively weak binding (e.g., 5-7 bases). A strand is stably attached to the
register strand only if it is bound by at least two domains, but one domain is
sufficiently long to act as a “toehold” to help initiate strand displacement
(Figure 2(c-f)). Within a cell with $d$ domains, each domain is unique and
assumed to be named $1,2,\ldots,d$. The register strand has the starred
version of these domains, e.g.,
$1^{*},2^{*},3^{*},4^{*},5^{*},6^{*},1^{*},2^{*},3^{*},4^{*},5^{*},6^{*}$
reading from the register’s $3^{\prime}$ to $5^{\prime}$ end (left to right)
in Figure 2(a). All cells have the same ordered list of domains, so for
example in Figure 2(a), domain 5 in cell 1 is the same DNA sequence as domain
5 in cell 2.

An _instruction_ is a set of strands that are added to the solution at once.
Figure 2(b-f) shows the various reactions that these strands might conduct to
change the configuration of strands attached to the register. Multiple
reactions can occur in a cascade in a single instruction.111 See for example
instruction 39 in Figure 8. In the right cell, an orange instruction strand
displaces an orange strand bound to the register via toehold exchange. This
opens a toehold for a blue instruction strand to displace the bound blue
strand, resulting in the configuration shown at the beginning of instruction
41.  In particular, the model is nondeterministic, and in general multiple
reactions might be possible. It is the job of the system designer to ensure
that only one final configuration can result no matter the order of reactions.
Instruction strands can either be _forward_ ($3^{\prime}$ arrow on right) or
_reverse_ ($3^{\prime}$ arrow on left). Forward instruction strands can do
attachment and displacement reactions (Figure 2(b,d-f)) and reverse
instruction strands can detach forward strands previously bound to the
register (Figure 2(c)).

Crucially, instruction strands are added in large excess over the register
strands. Thus even the toehold exchange displacement, which is often
considered reversible due to being enthalpically balanced (same number of
domains bound before and after), is actually irreversible in the SIMD||DNA
model, due to the entropic bias toward binding the new instruction strand with
much larger concentration than the strand it displaces.

The notation of [wang2019simd] uses dashed lines for reverse strands used for
detachment, as a visual reminder that they do not bind to the register. We
leave reverse strands as solid lines and rely on the $3^{\prime}$ arrow to
denote that the strand is reversed. We reserve the dashed line notation for
later figures to depict _inert_ instruction strands: instruction strands that
are shown above the register where they would bind if possible, but where no
reaction allows them to do so in the current configuration. We also have a
slightly different notation for strands with domains mismatching the register:
in [wang2019simd], these are depicted by writing an explicit domain name. In
our convention, the drawing of that part of the strand as diagonal and lying
entirely above the register domain indicates that the top strand domain and
register domain are not complementary (Figure 2(a), cell 2, domain 3). To
denote that two adjacent top strands share the same domain, both of which can
bind to the register (so they dynamically compete with strand displacement),
we draw both strands partially horizontal over the domain, and partially
diagonal (Figure 2(a), cell 1, domain 6).

Although these rules allow for nondeterministically competing reactions, our
construction is deterministic in the sense that there is only one sequence of
reactions possible in any instruction step.

After instruction strands are added and the described reactions go to
completion, the _wash_ step removes all strands not bound to the register.
This includes excess instruction strands that never reacted, as well as
strands that were displaced or complexes formed in a detachment reaction.

##  3 Simulation of Turing machine in SIMD||DNA

In this section we describe how to simulate an arbitrary 3-symbol single-tape
Turing machine with SIMD||DNA instructions.

###  3.1 High-level overview of construction

Since the SIMD||DNA model as defined has no mechanism to grow the register
strand, it can only simulate a fixed-space-bound Turing machine (a.k.a.,
linear-bounded automaton), which starts with $s$ total tape cells and never
moves the tape head off of them. A 3-symbol, space-$s$ Turing machine has
three tape symbols: $0,1,\llcorner\negthinspace\lrcorner$. The binary input
$x\in\\{0,1\\}^{<s}$ is represented by string
$x\llcorner\negthinspace\lrcorner^{s-|x|}$ on the tape in the initial
configuration, i.e. $x$ padded with enough blank symbols to make $s$ total
tape cells. We use as a running example the 5-transition Turing machine in
Figure 3, which increments a binary number.

Figure 3:  Turing machine (start state $a$) that increments an integer
represented in binary, with the least significant bit on the right. This
example is simulated in all subsequent figures.

Each cell of the register represents a tape cell of the Turing machine. If the
Turing machine has $t$ total transitions, then each cell uses $d=2t+8$
domains.

For each Turing machine, there is a fixed sequence of instructions that, after
executing, will update the register to represent the next configuration of the
Turing machine.

The cell with the tape head is the only cell with uncovered register domains.
Which domains are uncovered (known as a _transition region_) represents both
the current state of the Turing machine and the symbol written on that tape
cell. For all other cells, a disjoint region (the _symbol region_) represents
the symbol on that cell through its pattern of nicks. On the cell with the
tape head, the symbol region has no nicks (and represents no symbol) since it
is covered by a longer 8-domain strand.

###  3.2 Representation of Turing machine tape cell as a register cell

In the SIMD||DNA representation of a Turing machine, each register cell
represents a single Turing machine tape cell. We represent each Turing machine
with tape alphabet $\Gamma=\\{0,1,\llcorner\negthinspace\lrcorner\\}$, state
set $Q$, and halt state $h$, as a set of _transitions_ , where each transition
$(q,b)\to(r,c,m)$ means that if the Turing machine is in state $q\in
Q\setminus\\{h\\}$ reading symbol $b\in\Gamma$, it changes to state $r$,
writes symbol $c$, and moves one cell by $m\in\\{L,R\\}$ (left or right).
Since the Turing machine is deterministic, for each state-symbol pair, there
is at most one transition with that pair on the left. (But some such pairs
could be undefined, e.g., there is no
$(b,\llcorner\negthinspace\lrcorner)\to\ldots$ transition in Figure 3.)

#### Representation of tape cell with tape head

Figure 4:  A SIMD||DNA cell where the tape head is presently located. The
$(a,0)$ region is fully exposed, indicating that the Turing machine is in
state $a$ and that the cell contains the symbol $0$. The other transition
regions are fully covered, and the symbol region (rightmost 8 domains of the
cell) is covered by a single long strand, not encoding any symbol (which is
encoded by the uncovered transition region).

We take every state-symbol pair $(q,\sigma)\in(Q\setminus\\{h\\})\times\Gamma$
(each possible left side of a transition) and represent each as two
consecutive domains in a SIMD||DNA register cell. See Figure 4. Recall the
binary incrementing Turing machine of Figure 3. It has five transitions:
$(a,0)\to(a,0,R),\quad(a,1)\to(a,1,R),\quad(a,\llcorner\negthinspace\lrcorner)\to(b,\llcorner\negthinspace\lrcorner,L),\quad(b,0)\to(1,h,L),\quad(b,1)\to(0,b,L)$.
We call the pair on the left the _transition input_. Each of the given
transition inputs is represented in the SIMD||DNA cell using two domains,
requiring ten domains total for our example. Since each register cell
represents a cell in $M$, we must denote the presence of the tape head on one
of the cells. If the tape head is present on a given cell and if the current
Turing machine configuration has a valid transition, then the two domains that
represent that transition will have no top strand attached to them, leaving
them exposed. For example, if the tape head is on a cell with the $0$ symbol,
and the Turing machine is currently in state $a$, then the region that
represents $(a,0)$ in that cell will be exposed to serve as a toehold for
strand displacement. The other transition regions are fully covered by
2-domain strands.

#### Representation of tape cell without tape head

If the tape head is not present on a cell, or if no valid transitions exist
for the current configuration,222 For example, if the machine has halted; see
the bottom register configuration of Figure 9 for a case where the state is
non-halting but no valid transition exists.  then every transition region is
covered by 2-domain strands. Eight additional domains at the rightmost part of
the cell, called the symbol region represent the current symbol written on
that cell.

Figure 5:  High-level overview of construction. A Turing machine register
currently in state $a$, with the tape head on the second cell. The second cell
contains the symbol 1. The leftmost cell contains symbol 0. The inset above
shows encodings for 1 and $\llcorner\negthinspace\lrcorner$. All the
transition regions are fully covered in cells lacking the tape head. After the
full list of instructions in the SIMD||DNA program are complete, the register
represents the Turing machine configuration with state $a$ and the tape head
moved to the rightmost cell with the $\llcorner\negthinspace\lrcorner$. The
same full list of instructions updates the register again, now representing
the Turing machine configuration in state $b$ with the tape head back on the
middle cell.  Figure 6:  An overview of the first 19 instructions of the
construction, which represent the $(a,0)\to(a,0,R)$ transition of the Turing
machine shown in Figure 3. Instructions marked with a red cross are fully
inert, meant for cases not exhibited by this register.

Whenever the tape head is present on a cell, the symbol region is covered by a
single 8-domain strand that does not encode any symbol, since the symbol
information is already encoded in the transition region with an open toehold.

###  3.3 Detailed description of SIMD||DNA instructions simulating a Turing
machine

We designed an algorithm that converts Turing machine specifications from
https://turingmachine.io into SIMD||DNA representations, along with the
equivalent instructions. Each transition $\tau_{i}$ has an associated sublist
of instructions $L_{i}$, and, not knowing which transition is applicable to
the current configuration, we simply add instruction strands in order from
$L_{1},L_{2},\ldots$. For $i\neq j$, to ensure that $L_{j}$ instructions have
no effect when the current applicable transition is $\tau_{i}$, we “plug” the
open domains of other transition regions with a strand and remove the plug
strand once it’s time to process that transition. Because the SIMD||DNA model
allows parallel computation among multiple registers in the same solution,
this prevents instructions meant for one configuration from affecting
registers currently not in that configuration. In the beginning, all
transition regions are plugged, where the order of processing for the
transitions is arbitrary.

Figure 7:  On the left, the first instruction in the whole SIMD||DNA program
covers the transition region representing the next applicable transition. The
two horizontal rows of strands have the following interpretation: Bottom are
strands bound to register, top are instruction strands. Dashed instruction
strands will not have an effect on the current cell (but to help verify
correctness, they are shown above where they would bind to the register). On
the right, the last instruction in the whole SIMD||DNA program, which removes
the post-plug strands in each register. In the above example, the post-plug
strand covers the $(a,0)$ transition region, indicating the cell’s next Turing
machine transition. After this, the entire register is updated to appear as a
configuration similar to those in Figure 5.

#### Pre-plug and post-plug strands to protect instructions for inapplicable
transitions from affecting configuration

The full list of instructions to simulate a Turing machine transition works as
follows. Recall that in the “clean” configurations shown in Figure 5, the only
exposed register domains are on the cell representing the tape head. The first
instruction in the entire list contains _pre-plug_ strands for each transition
region. At the end of each instruction sublist $L_{i}$, a _post-plug strand_
is also placed on the transition region that represents the Turing machine’s
next applicable transition. Examples of both can be seen in Figure 7. These
strands act like a “chemical protecting group” that prevents instruction
sublists $L_{i}$ from modifying the register unless they apply to the intended
transition. The difference between the pre-plug and post-plug strands is that
pre-plug strands protect the configuration when using instructions strands
_before_ the applicable transition, whereas post-plug strands protect the
configuration when using instructions strands _after_ the applicable
transition. In the left part of Figure 7, because the register has a pre-plug
strand in $(a,\llcorner\negthinspace\lrcorner)$, it means that the instruction
sublist $L_{(a,\llcorner\negthinspace\lrcorner)}$ has not been applied to it
yet. Instruction sublists for the other transitions
$(a,0)\to\ldots,\quad(a,1)\to\ldots,\quad(b,0)\to\ldots,\quad(b,1)\to\ldots$
will be inert, not affecting the register. The first instruction of
$L_{(a,\llcorner\negthinspace\lrcorner)}$ will remove this pre-plug strand so
that any register in the $(a,\llcorner\negthinspace\lrcorner)$ configuration
can be processed. The instruction sublists will result in a configuration like
that of the bottom of Figures 8 and 9. Figure 9 shows instructions that affect
the cell where the tape head _was_ (right cell), not where it will be _next_
(left cell), which is why the left cell is the same in both Figures 8 and 9.
This almost represents the next Turing machine configuration, but with the
appropriate transition region covered by a _post-plug strand_. The final
instruction in the entire list (Figure 7) removes this post-plug strand,
restoring the register configuration to be as shown in Figure 5.

The post-plug strand placed on the transition region at the end of simulating
a transition has a different purpose from the pre-plug strand placed in
instruction 1. Its purpose is to prevent the register from updating its state
multiple times in the same instruction iteration. For example, if a register
has a post-plug strand on the $(b,1)$ region (indicating that it has been
processed and that its next transition is $(b,1)$), and the instruction
sublist that processes $(b,1)$ comes after, the register will be unaffected by
$(b,1)$’s deprotecting instruction, keeping it inert throughout. The final
instruction in the entire iteration removes these post-plug strands from the
registers, as seen in Figure 7, preparing the registers for the next iteration
of the instruction set.

Note the duality between pre-plug and post-plug instructions. _All_ pre-plug
strands are included in the first instruction, though only one of them will
bind (the one matching the applicable transition), and instruction strands
removing _all_ post-plug strands are included as part of the last instruction,
though only one will find its complementary post-plug strand to remove. On the
other hand, each pre-plug instruction is removed more specifically, by adding
a single complementary strand to remove it just prior to the sublist of
instructions corresponding to the applicable transition. Similarly, each post-
plug strand is added by itself, at the end of the instruction sublist
corresponding to the applicable transition.

In the next section, we will describe the details of the instruction sublists
that represent the Turing machine transitions.

#### Sublist of instructions representing a single Turing machine transition

_Figure conventions._ For the figures explaining SIMD||DNA instructions that
simulate a single transition of the Turing machine (Figure 8 and beyond), we
use the following conventions in figures. Several register configurations are
shown, but they are not necessarily consecutive. Each is numbered with its
absolute index in the list of all 86 instructions implementing the Turing
machine of Figure 3. An example with all instructions shown can be found in
Appendix C. If two adjacent configurations have non-consecutive instruction
indices, this means that the instructions not shown are inert: their strands
do not affect the register in that configuration. The instruction strands that
have just been added are always shown above the register, with a solid line if
they will do a reaction as in Figure 2, and with a dashed line if that
instruction strand is inert for that configuration. The final configuration in
each figure does not show any instruction strands, but for all figures there
is a followup figure showing what happens next from that configuration
(possibly the followup is Figure 7, the final instruction in the entire
program, removing the post-plug strand from the next applicable transition
region).

Each Turing machine transition is individually processed by a sublist of
instructions. The pre-plug strand is first removed by an instruction
containing its complementary strand, so that its corresponding transition
region in the cell can be used as a toehold, such as instruction 38 in Figure
8. The next instructions then update the contents of the current cell to
encode the symbol that the tape head writes. For example, in Figure 9, the
transition
$(a,\llcorner\negthinspace\lrcorner)\to(b,\llcorner\negthinspace\lrcorner,L)$
is represented, and the strand encoding of $\llcorner\negthinspace\lrcorner$
is placed in the right cell after the tape head writes on it and moves left.
After that, the instructions check the contents of the tape head’s new
location and determine the Turing machine’s next configuration. In Figure 8,
the tape head moves to the left cell and finds a 1, and the Turing machine
goes to state $b$, so it leaves a post-plug strand on $(b,1)$’s transition
region to show that the register has been processed for that instruction
iteration, as seen in instruction 46.

#### Left versus right tape head moves

In the SIMD||DNA instructions implementing a single transition, there are two
sublists: next-cell instructions and previous-cell instructions. As their
names indicate, next-cell instructions update the contents of the tape head’s
destination, while previous-cell instructions update the contents of the tape
head’s former location. Other factors such as the symbol to be written on the
current cell and the next applicable transition region only introduce minor
variations in the instruction strands.

For transitions moving the tape head _left_ , the next-cell instructions
precede the previous-cell instructions (see Figures 8, 10, 11 and 9). For
transitions moving the tape head _right_ , this order is reversed (see Figures
12, 13, 14 and 15).

#### Left tape head moves

Figure 8 shows the next-cell instructions for transition
$(a,\llcorner\negthinspace\lrcorner)\to(b,\llcorner\negthinspace\lrcorner,L)$,
for the special case when the cell to the left of the tape head has the symbol
1. Figure 10 shows the next-cell instructions for the same transition when the
cell to the left of the tape head has the symbol 0, and Figure 11 shows the
next-cell instructions when the cell to the left has the symbol
$\llcorner\negthinspace\lrcorner$. Note that in any given configuration, the
same instructions will result in exactly one of the situations depicted in
Figures 8, 10 and 11. Once these next-cell instructions are applied, the
leftmost domain of the previous cell will serve as a toehold for the previous-
cell instructions that follow in Figure 9.

Figure 8:  First half of instructions (next-cell instructions) to implement
transition $a,\llcorner\negthinspace\lrcorner\to
b,\llcorner\negthinspace\lrcorner,L$, in the case that the cell to the left
(where the tape head will move) has a symbol 1. After instruction 46, the left
cell (where the tape head will move next) now encodes the next state $b$ and
the symbol 1 on the new tape cell. Rather than show a red cross, inert
instructions are instead omitted in this figure and the following ones. Inert
instructions (40,42,47,48) are used when the cell to the left has a symbol 0
or $\llcorner\negthinspace\lrcorner$ on it instead. Figures 10 and 11
respectively show these cases. Figure 9 shows instructions completing the
transition by writing $\llcorner\negthinspace\lrcorner$ over the right cell.
In instruction 38, the transition region is first unplugged to expose the
toehold. The strands in instruction 39 cascade until the left cell’s symbol
region, allowing the instructions to branch out depending on the tape content
of the left cell. The remaining instructions process the left cell so that it
encodes the next configuration. In instruction 46, a special post-plug strand
whose leftmost domain is orthogonal is attached to the region that represents
the next configuration; this strand will be removed once all transitions have
been processed. The angled dashed strands to the right of the register
indicate that no cells are present to the right of the rightmost cell; if
there were another cell, then one more domain of each of these strands would
be horizontal, bound to the leftmost domain of the cell to the right (just as
with their solid counterparts to the left). Figure 9:  Second half (previous-
cell instructions) of transition $a,\llcorner\negthinspace\lrcorner\to
b,\llcorner\negthinspace\lrcorner,L$ whose first 11 instructions are shown in
Figure 8, Figure 10, and Figure 11; these instructions write
$\llcorner\negthinspace\lrcorner$ over the old cell (right of Figure 8) where
the tape head was at the start of the transition. Slight variations of
instruction 15 write 0 or 1 instead of $\llcorner\negthinspace\lrcorner$.
Figure 10:  The next-cell instructions of the transition in Figure 8 that are
applicable when the cell to the left of the tape head has a 0. Figure 11:  The
next-cell instructions of the transition in Figure 8 that are applicable when
the cell to the left of the tape head has a $\llcorner\negthinspace\lrcorner$.
Because no transition exists for the state-symbol pair
$(b,\text{$\llcorner\negthinspace\lrcorner$})$, the left cell is left
unchanged.

#### Right tape head moves

For right transitions, the first three instructions are the previous-cell
instructions, as shown in Figure 12. The next-cell instructions follow, where
Figure 13 shows what happens when the cell to the right of the tape head has a
$\llcorner\negthinspace\lrcorner$, Figure 14 shows what happens when the cell
to the right has a 0, and Figure 15 shows what happens when the cell to the
right has a 1. As with the left-cell instructions, the next-cell situations
depicted in Figures 13, 14 and 15 are disjoint.

Figure 12:  First half (previous-cell instructions) of transition $a,0\to
a,0,R$. Instruction 3’s strands cascade to the right side of the current cell,
while instruction 4 removes the previously introduced strands. Instruction 5
then covers up all the transition regions and adds the strands of the symbol
to be written on the current cell (0 in this example), but leaving the
rightmost domain exposed to act as a toehold for the next instructions.

#### Final deprotecting instruction

Figure 7 shows the final deprotecting instruction, which removes the post-plug
strand put in place during the last instruction of the non-inert instruction
sublists. This puts the register back into a “clean” configuration
representing a Turing machine configuration, such as those shown in Figure 5,
opening up new toeholds for the next iteration of instructions.

###  3.4 Complexity of construction

A common metric of “complexity” of DNA systems is the number of unique domains
they require. Fewer is better because it is a nontrivial task to design
_orthogonal_ domains: domains that, if they are not perfectly complementary,
will have low binding affinity. Low domain complexity is particular important
in SIMD||DNA, where each domain is considered “toehold-length”: sufficiently
short (5-7 bases) that the off-rate of a strand bound by a single domain is
large enough to detach in a short amount of time. There are only $4^{7}=16384$
DNA sequences of length $7$. In practice even fewer are available: half are
complementary to the other half, leaving only $8192$ available to assign to
the unstarred versions of each domain. DNA sequence design heuristics such as
the “3-letter code” (using only A,T,C for forward strands, thus only A,T,G for
the register and reverse strands) reduce this number to $3^{7}=2187$.
Reasonable design constraints, e.g., avoiding almost equal domains such as
$5^{\prime}$-AAAAAAG-$3^{\prime}$ and $5^{\prime}$-AAAAAAA-$3^{\prime}$, which
both bind almost equally strongly to $3^{\prime}$-TTTTTTT-$5^{\prime}$,
further limit the set of available domains.

Our construction uses $d=2t+8$ total unique DNA domains (which repeat
throughout the register), where $t$ is the number of transitions of the
simulated Turing machine. Each transition is represented by 2 domains, assumed
to be bound strongly enough not to spontaneously dissociate. (The construction
can be altered to use more domains, in case this assumption is overly ideal.)
To simulate a Turing machine with space bound $s$, the register has $s$
“cells”, where each cell is simply a copy of one each of the $d$ domains
$1,\ldots d$. Thus, if each domain consists of $k$ nucleotides, the register
strand has $k\cdot s\cdot(2t+8)$ total nucleotides. There is an 11-state,
3-symbol universal Turing machine (directly simulating another Turing machine)
with 32 transitions [neary2006small], giving $2\cdot 32+8=72$ total domains
required in the worst case. However, specialized non-universal Turing machines
with a smaller number of transitions (for example the 5-transition binary
incrementor of Figure 5) could accomplish many computationally sophisticated
tasks.

Each Turing machine transition can be represented by approximately 16
SIMD||DNA instructions. The exact number varies depending on the specific
properties of the transition in question, such as the the direction of the
tape head, the symbol to be written, and whether any of the possible next
configurations are halting or not. This range has constant upper and lower
bounds, however, so the total number of instructions is in $O(t)$, where $t$
is the number of transitions in the Turing machine. Because the size of the
cell increases in $O(t)$ due to the transition regions, the number of DNA
strands present in some instructions also scales by a factor of $O(t)$, most
notably the instructions that cause a cascade of toehold exchanges, such as
instructions 3 to 7 in Figures 12 and 13.

##  4 Conclusion

Our construction, like the Rule 110 simulation of Wang, Chalk, and Soloveichik
[wang2019simd], is not Turing universal because it simulates a _space-bounded_
Turing machine. Truly universal computation should be possible without
advanced knowledge of the space requirement. An interesting question (raised
also in [wang2019simd]) is whether a suitable augmentation of the SIMD||DNA
model could allow Turing-universal computation. This would require unbounded
polymers such as those used in the two-stack machine DNA implementation of
Qian, Soloveichik, and Winfree [qian2010efficient]. That paper showed Turing
universal computation in the case where only a single copy of certain strands
are permitted to exist in solution, simulating only a single stack machine at
a time, in contrast to the SIMD||DNA model, where we can operate on many
registers, each representing their own Turing machine, in parallel.

Technically the strand displacement reactions of the SIMD||DNA model are
currently powerful enough to grow arbitrarily large polymers from a fixed set
of strands, as in [qian2010efficient], by alternating top and bottom strands,
making a long double-helix with nicks on the top _and_ bottom. However, it is
difficult to see how to use the ability of SIMD||DNA to exploit this to
simulate a Turing machine represented in this way. If there are multiple
bottom strands, i.e., there is a nick on the bottom, then any strand
displacement of top strands, upon reaching this nick, would separate the
polymer into two complexes to the left and right of this nick, and the right
polymer would be lost in the wash step. One could imagine, however, augmenting
the model to allow, for example, 3-arm junctions, which could be used to do
strand displacement that crosses over the boundary between two bottom strands
without separating them (since they would be joined to each other by a strong
domain representing the third arm “below” the main helix).

Although some of the toehold exchanges in the SIMD||DNA model are reversible
based on the principles of DNA strand displacement, we make the assumption
that the applied instructions are not undone by the displaced strands. This is
based on the assumption that the instruction strands are present in a
sufficiently high concentration that reversal is unlikely. Because multiple
registers can be present in the same solution, another possibility to consider
is a displaced DNA strand from Register A binding to an open toehold in
Register B, such that the attachment is irreversible even with the presence of
a high concentration of instruction strands. One open question is to design a
system that factors in these possibilities, reducing the likelihood of
unexpected strand displacement results.

Another open question is whether a more domain-efficient encoding exists for
the Turing machine construction. Given $n$ transitions, $2n$ domains are
required to represent them, which has $O(n)$ complexity. However, given $d$
domains, there are $2^{d-1}$ possible nick patterns among the attached
strands, which makes $O(\log n)$ domain complexity possible in theory.

## References

  * [1] James Bornholt, Randolph Lopez, Douglas M Carmean, Luis Ceze, Georg Seelig, and Karin Strauss.  A DNA-based archival storage system.  In ASPLOS 2016: Proceedings of the Twenty-First International Conference on Architectural Support for Programming Languages and Operating Systems, pages 637–649, 2016.
  * [2] Tonglin Chen, Arnav Solanki, and Marc Riedel.  Parallel Pairwise Operations on Data Stored in DNA: Sorting, Shifting, and Searching.  In DNA 27: 27th International Conference on DNA Computing and Molecular Programming, volume 205, pages 11:1–11:21, 2021.  doi:10.4230/LIPIcs.DNA.27.11.
  * [3] George M Church, Yuan Gao, and Sriram Kosuri.  Next-generation digital information storage in DNA.  Science, 337(6102):1628–1628, 2012.
  * [4] Matthew Cook.  Universality in elementary cellular automata.  Complex systems, 15(1):1–40, 2004.
  * [5] Turlough Neary and Damien Woods.  P-completeness of cellular automaton Rule 110.  In ICALP 2006: International Colloquium on Automata, Languages, and Programming, pages 132–143. Springer, 2006.
  * [6] Turlough Neary and Damien Woods.  Small fast universal Turing machines.  Theoretical Computer Science, 362(1-3):171–195, 2006.
  * [7] Lee Organick, Siena Dumas Ang, Yuan-Jyue Chen, Randolph Lopez, Sergey Yekhanin, Konstantin Makarychev, Miklos Z Racz, Govinda Kamath, Parikshit Gopalan, Bichlien Nguyen, Christopher N Takahashi, Sharon Newman, Hsing-Yeh Parker, Cyrus Rashtchian, Kendall Stewart, Gagan Gupta, Robert Carlson, John Mulligan, Douglas Carmean, Georg Seelig, Luis Ceze, and Karin Strauss.  Random access in large-scale DNA data storage.  Nature Biotechnology, 36(3):242, 2018.
  * [8] Lulu Qian, David Soloveichik, and Erik Winfree.  Efficient Turing-universal computation with DNA polymers.  In International Workshop on DNA-Based Computers, pages 123–140. Springer, 2010.
  * [9] Georg Seelig, David Soloveichik, David Yu Zhang, and Erik Winfree.  Enzyme-free nucleic acid logic circuits.  Science, 314(5805):1585–1588, 2006.
  * [10] SIMD||DNA simulator.
source code: https://github.com/UC-Davis-molecular-computing/simd-dna, 2021.

  * [11] S Kasra Tabatabaei, Boya Wang, Nagendra Bala Murali Athreya, Behnam Enghiad, Alvaro Gonzalo Hernandez, Christopher J Fields, Jean-Pierre Leburton, David Soloveichik, Huimin Zhao, and Olgica Milenkovic.  DNA punch cards for storing data on native DNA sequences via enzymatic nicking.  Nature Communications, 11(1):1–10, 2020.
  * [12] Boya Wang, Cameron Chalk, and David Soloveichik.  SIMD||DNA: Single instruction, multiple data computation with DNA strand displacement cascades.  In DNA 2019: International Conference on DNA Computing and Molecular Programming, pages 219–235, 2019.

##  Appendix A Appendix: Next-cell Instructions

Figure 13: The next-cell instructions for when the right cell contains a
blank. The Turing machine enters the
$(a,\text{$\llcorner\negthinspace\lrcorner$})$ configuration afterward. Figure
14: The next-cell instructions for when the right cell contains the symbol 0,
entering the $(a,0)$ configuration. Figure 15: The next-cell instructions for
when the right cell contains the symbol 1, entering the $(a,1)$ configuration.

##  Appendix B Appendix: $(a,1)$ Transition Overview

Figure 16:  An overview of instructions 20 to 37 of the construction, which
represent the $(a,1)\to(a,1,R)$ transition of the Turing machine shown in
Figure 3. Because the register in this example is not in the $(a,1)$
configuration, the instructions in this section do not affect the register, as
no toehold is present.

##  Appendix C Appendix: Full Instruction Overview

The following figures show the full instruction set on a register whose
current value is 01$\llcorner\negthinspace\lrcorner$, where the tape head is
on 1 and the current state is $a$. Because the $(a,1)$ transition is
represented in this configuration, the instruction sublists that represent the
other transitions (Figures 17, 19, 20 and 21) are fully inert.

Figure 17: The initial pre-plug instruction, as well as the $(a,0)$
instruction sublist, for a register whose next transition is $(a,1)$. Note all
instructions after the first are inert, because the pre-plug strand blocks all
subsequent instructions (as intended since $(a,0)$ is not the applicable
transition.) Figure 18: The $(a,1)$ instruction sublist for a register whose
next transition is $(a,1)$. Instructions 30 to 37 are inert, because they
apply to the cases when the cell to the right of $(a,1)$ contains a 0 or a 1.
Figure 19: The $(a,\llcorner\negthinspace\lrcorner)$ instruction sublist after
the register was processed by the $(a,1)$ instruction sublist. Even though the
register is now in the $(a,\llcorner\negthinspace\lrcorner)$ configuration,
the post-plug strand introduced in instruction 25 prevents the register from
being updated a second time using the $(a,\llcorner\negthinspace\lrcorner)$
instructions. Figure 20: The $(b,0)$ instruction sublist after the register
was processed by the $(a,1)$ instruction sublist. Figure 21: The $(b,1)$
instruction sublist after the register was processed by the $(a,1)$
instruction sublist, as well as the final deprotecting instruction that
removes the post-plug strand.