Evolving Programmable Computational Metamaterials
Atoosa Parsa

Dong Wang

Corey S. O’Hern

University of Vermont
Burlington, Vermont, USA
atoosa.parsa@uvm.edu

Yale University
New Haven, Connecticut, USA
dong.wang@yale.edu

Yale University
New Haven, Connecticut, USA
corey.ohern@yale.edu

Mark D. Shattuck

Rebecca Kramer-Bottiglio

Josh Bongard

City College of New YorkNew York, New York, USA
shattuck@ccny.cuny.edu

Yale University
New Haven, Connecticut, USArebecca.kramer@yale.edu

University of Vermont
Burlington, Vermont, USA
josh.bongard@uvm.edu

Figure 1: Overview. (a): Multiobjective optimization (MOO) evolves increasing amounts of computational potential into granularmetamatetrials. Each individual in the population is a granular metamaterial composed of two particle types. Different
configurations of particles confer different material behaviours. In the approach reported here, vibrations are supplied as input
and vibration (if any) is recorded as output. Materials have been found that act as an AND gate (a) or an XOR gate (b). We report
here how MOO can discover a single material (d) that, at one frequency acts as an AND gate, and at another frequency acts as
an XOR gate. Thus, the superposition of input waves supplied to the system will result in emergent behaviors other than what
the material was originally designed to perform. This suggest future materials amenable to reprogramming using increasingly
sophisticated programming languages expressed in the frequency domain.

ABSTRACT

the designed material can be “reprogrammed” via frequency modulation. As metamaterial design advances, more computationallydense materials may be evolved, amenable to reprogramming byincreasingly sophisticated programming languages written in thefrequency domain.

Digital signal processors are widely used in today’s computers to
perform advanced computational tasks. But, the selection of digitalelectronics as the physical substrate for computation a hundredyears ago was influenced more by technological limitations thansubstrate appropriateness. In recent decades, advances in chemical,
physical and material sciences have provided new options. Granularmetamaterials are one such promising target for realizing mechanical computing devices. However, their high-dimensional design
space and the unintuitive relationship between microstructure anddesired macroscale behavior makes the inverse design problem
formidable. In this paper, we use multiobjective evolutionary optimization to solve this inverse problem: we demonstrate the designof basic logic gates embedded in a granular metamaterial, and thatCCS CONCEPTS• Computing methodologies → Randomized search; • Theory of computation → Models of computation; • Hardware →Emerging technologies.

KEYWORDS
Granular Metamaterials, Mechanical Computing, Logic Gates, Inverse Design Problem, Multiobjective Optimization

Permission to make digital or hard copies of all or part of this work for personal orclassroom use is granted without fee provided that copies are not made or distributedfor profit or commercial advantage and that copies bear this notice and the full citationon the first page. Copyrights for components of this work owned by others than ACM
must be honored. Abstracting with credit is permitted. To copy otherwise, or republish,
to post on servers or to redistribute to lists, requires prior specific permission and/or afee. Request permissions from permissions@acm.org.
GECCO ’22, July 9–13, 2022, Boston, MA, USA© 2022 Association for Computing Machinery.
ACM ISBN 978-1-4503-9237-2/22/07. . . $15.00
https://doi.org/10.1145/3512290.3528861

ACM Reference Format:
Atoosa Parsa, Dong Wang, Corey S. O’Hern, Mark D. Shattuck, Rebecca
Kramer-Bottiglio, and Josh Bongard. 2022. Evolving Programmable Computational Metamaterials. In Genetic and Evolutionary Computation Conference(GECCO ’22), July 9–13, 2022, Boston, MA, USA. ACM, New York, NY, USA,
8 pages. https://doi.org/10.1145/3512290.3528861

122

GECCO ’22, July 9–13, 2022, Boston, MA, USA

1

Parsa et al.

INTRODUCTION

on the intuitions of a human designer. In this paper, we propose
using evolutionary algorithms to automatically optimize the material properties for dense computation. Moreover, the mentioned
works are mostly conducted using continuous rather than granularmetamaterials, or in one-dimensional particle chains, while ourwork is performed with two-dimensional sheets of granular metamaterial. Due to advantages of this type of material for dynamic
programmability (as we will show), we envision a wider potential
in expanding our work to more complex computational elements,
as explained in our conclusions.
Granular metamaterials possess many parameters that affecttheir behavioral response to stimuli. For example, the particle properties (size, shape, stiffness, mass, and so on) and their placements
can affect the eigenfrequencies of the system and consequently the
localized propagation or suppression of acoustic waves throughthe material. With so many design parameters, deciding on the
optimal micro-structure to achieve a desired macro-behavior (i.e. a
logic gate) is a challenge. Evolutionary algorithms have long beenshown capable of automated optimization in such design spaces[15], [18]. Given this, in previous work [19] we demonstrated that
we can evolve the configurations of granular metamaterials suchthat they act as logic gate. Here, we show that a MOO can densely
pack multiple operations into the same piece of material: at one
frequency the material acts as an AND gate; at another frequency
it acts as an XOR gate. Finally, we show that combining two different input frequencies can enable the material to simultaneously
compute both logical functions. This suggest future materials may
be amenable to reprogramming using increasingly sophisticatedprogramming languages expressed in the frequency domain.

Metamaterials are an emerging class of engineered composite materials that exhibit properties different from their constituent materials and behave in ways not observed in nature [12]. In the same
way that metamaterials exhibit non-intuitive and exotic material
properties, computational metamaterials have potential to perform
computation in useful ways as well [23, 28]. Here, we explore various ways to evolve a particular class of metamaterial—granularmetamaterials—to increase its computational capabilities. Granularmetamaterials are metamaterials consisting of discrete particles thatexhibit increased plasticity compared to continuous metamaterials
because they can be dynamically reconfigured: particle properties
can change in response to external stimuli [26]. In this paper, we
propose using granular metamaterials as a physical substrate toperform mechanical computation.
Considering logic gates as the basic computational blocks uponwhich more complex units can be built, we demonstrate here thedesign of metamaterials that can act as basic acoustic logic gates:
low or high vibration frequencies or amplitudes can serve as zeros
or ones. There are several advantages of this type of computationcompared to digital computation where logic gates are built fromelectrical transistors. First, by moving to a mechanical substrate
we can avoid analogue to digital conversion, sidestepping all of
the limitations of abstract representations and discretizations necessary for a digital computing system [28]. Second, by supplying
vibrations composed of multiple frequencies, multiple computations may be performed simultaneously using the same patch ofmaterial, suggesting increasingly computationally dense materials
may be evolved in future. Third, acoustically driven computational
metamaterials could serve as useful building blocks for more complex machines such as robots: sensing, control and actuation could
all respond to and/or produce vibration, rendering these otherwise
separate robot components as just different regions within a singlematerial robot. Finally, by exploiting the natural dynamics of themetamaterial, energy efficient computation and higher robustness
and stability may be achieved. Moreover, because computational
components do not need to be physically separate modules (as weshow below), this raises the possibility of bottom-up design of computational architecture where the exact form of computation is notpredetermined [17].
In recent years there has been some research on embedding mechanical computation into material. In [22], a universal logic gate is
implemented as a nonlinear mass-spring-damper model. In [20], asoft bistable building block is designed and used in the implementation of soft mechanical diodes and logic gates. [14] and [3] present
examples of acoustic gate design in a 1D chain of elastic particles.
Computational metamaterials have been introduced to perform
specialized computing tasks such as integration, differentiation andconvolution [28]. Perhaps closest to our work is [26], in which theauthors show the potential of a hand-designed 2D granular system
to act as an acoustic switch modulated by the packing pressure toswitch between on and off states. Despite these recent advances, innone of the aforementioned works is the computational unit automatically optimized to perform computation, let alone how best to
densely pack computation in new ways into materials is explored.
Instead, computational building blocks are hand-designed based

2

METHODS

In order to build a computing system, we need to choose a physicalsubstrate within which to embed the architecture of the computational model [27]. The conventional choice is digital electronic
devices. As discussed in the introduction, this choice was heavily influenced by the manufacturing technology of the time: the promiseof the semiconductor industry and Gordon Moore’s repeatedly validated prediction about the doubling of transistors on integrated
circuits every two years [16]. This promise led to the abandonmentof the other computational substrates and huge investments in digital electronics. General purpose digital computers in effect won
the “Hardware Lottery” and enjoyed significant advances during
the last 50 years [10]. However, the physical constraints of miniaturization are now causing Moore’s Law to slow significantly [24].
This has triggered a renewed interest in exploring other possible
physical substrates for computation. In addition to photonics [25],DNA computing [1] and quantum computing [8], Metamaterials
are one such promising computational substrate.
In this paper, we explore the computational capabilities of granular metamaterials. Granular metamaterials have been studied toa great extent in the material sciences literature [2, 11, 13]. They
are made from many individual grains with possibly different material properties. They can exhibit complex interesting behaviors
but there are (simplified) computational models that can accurately
predict them [21]. The basic principle that we use in this work is
the ability of a granular assembly to have extinguishable responses

123Evolving Programmable Computational Metamaterials

GECCO ’22, July 9–13, 2022, Boston, MA, USA

next section. The truth table in Fig. 1 shows our bit representation.
The output 𝑂𝑖 𝑗 (𝑖, 𝑗 ∈ 0, 1) is the amplitude of the displacement of
the particle at the output port.
In each of the four cases shown in the truth table in Fig. 1, thegain of the system is defined as the amplitude of the fast Fouriertransform ( 𝑓ˆ) at the driving frequency (𝜔) in the output, divided
by the sum of the amplitudes of the fast Fourier transform at thedriving frequency (𝜔) in the inputs. In other words:
Figure 2: The system setup. A logic gate with two inputs(green and blue particles) and one output (red particle) embedded in a 2D granular assembly with two types of particles(indicated with different shades of grey). 𝐴1 and 𝐴2 are theamplitudes of oscillations applied in 𝑥 direction to the input
ports. 𝜔 is the input frequency. The truth table is shown on
the right: 𝐷 represents the diameter of particles. 𝑂𝑖 𝑗 (𝑖, 𝑗 ∈ 0, 1)
is the magnitude of the vibrations in the output.

𝐺𝑖 𝑗 (𝜔) =𝑖, 𝑗 ∈ 0, 1

(1)

In our experiments, we can measure the gain (𝐺𝑖 𝑗 ) for each of
the four input cases. Then, in order for the material to act as alogic gate, the relative magnitude of the gain in each case must be
consistent with desired functionality of the gate. Note that because
it’s experimentally impractical to get an absolute zero or one, weare looking at the relative magnitudes of the two cases. For example,
for an AND gate when both input ports are driven with a sinusoidal
wave (𝑖𝑛𝑝𝑢𝑡 = ‘11‘), we expect to see a high amplitude of oscillation
at the output (𝑜𝑢𝑡𝑝𝑢𝑡 = ‘1‘) and therefore we expect a high gain
(𝐺 11 ). But in the other remaining three cases (𝑖𝑛𝑝𝑢𝑡 = ‘00‘, ‘01‘,and ‘10‘) we expect a low amplitude of vibration (𝑜𝑢𝑡𝑝𝑢𝑡 = ‘0‘)
and thus a low gain (𝐺 00, 𝐺 01 and 𝐺 10 ). Based on this formulation,
we can devise a metric to measure the similarity of the material’sfunctionality to a desired logic gate. We will show this metric in
the next section.

to mechanical vibrations with different frequencies. In the next
section, we expand this idea and formalize the problem setup.

2.1

𝑓ˆ(𝑂𝑖 𝑗 )𝑓ˆ(in𝑖 ) + 𝑓ˆ(in 𝑗 )Embedding Computation

In a typical electronic system, the inputs and outputs are electricalsignals. Here, we evolve an analogous mechanical system where
inputs and outputs are acoustic signals. Our system is a two dimensional assembly of two types of circular particles placed on ahexagonal lattice. The setup is shown in Fig. 2. As we mentioned
before, such granular system can be widely tuned by changing theparticles’ properties to achieve different responses [11]. The response is dependent on a large number of variables, such as particle
arrangements, mass, modulus, shape, interactions between particles, and boundary conditions of the system. These properties canaffect the frequency spectrum of the material, normal modes of thesystem, and consequently position of the band gap: a contiguous
region of the input frequency spectrum muffled by the material.
Our goal is to regulate the intrinsic material properties to obtain ametamaterial with a desired vibrational response. We are interested
in embedding computation in the granular metamaterial and thefirst step towards this goal is to see if we can design the basic logic
gates. In this system, the inputs and outputs are acoustic waves. So
we choose two particles on one side of the material as the input
ports (the particles with green and blue markers in Fig. 2) and one
particle on the other side as the output port (the particle with thered marker). The input signal is a sinusoidal wave with amplitude
𝐴𝑖 and frequency 𝜔 in 𝑥 direction. This signal is applied to a particle,
causing it to move from its initial position (𝑥𝑖0 ). We can look at theamplitude of the displacement signal in time (𝑥𝑖 (𝑡)) as the abstract
bit representation: small or no displacement means a ‘0‘ and highmagnitude of displacement means a ‘1‘. Going back to our analogy
for the input signal, applying a sinusoidal wave with amplitude zero
(𝐴𝑖 = 0) to an input port means a ‘0‘ bit at that port of the logic gate.
In our simulations, we fixed the non-zero amplitude to 1 × 10−2(which is 10% of the diameter of a particle.). The frequency of theapplied signal (𝜔) is also a fix value, chosen based on the frequency
spectrum of a regular configuration which will be discussed in the

2.2

The Simulator

The simulator 1 is a simplified granular system based on the model
used in [26]. It’s a two dimensional system, made of frictionlesscircular disks with a fixed diameter (𝐷). The particles can havedifferent material properties such as different masses or different
stiffnesses (this is shown as dark/light colors in the figures.). They
are placed on a 5 by 6 hexagonal lattice and so there is a total
number of 30 particles. The system is periodic in 𝑥 direction andhas a fixed boundary in 𝑦 direction. Gravity is ignored, so the only
forces acting on the particles are the result of a purely repulsive
linear spring potential between the disks which can be formalized
as a Lennard-Jones potential.
Our system is simulated using discrete element method (DEM).
Starting from the initial positions of the particles (placed on ahexagonal lattice), the repulsive forces are calculated based on the
distance between the overlapping particles. Next, the accelerations,velocities and positions are updated using the Verlet integration.
But there is a pre-processing step to make sure that the system is
at an equilibrium at 𝑡 = 0 (meaning that the sum of total forces
between particles is (near) zero) and thus the particle packing isstatistically stable. This is done by calculating the total force andupdating the initial positions of the particles using the steepestdescent method. Here, we use Fast Inertial Relaxation Engine (FIRE)
to reduce the processing load.
Table 1 includes the main simulation parameters and their assigned values in our experiments. 𝑁𝑡 is total simulation time.
1 https://github.com/AtoosaParsa/gecco-2022 contains the source code necessary

for reproducing our results.

124

GECCO ’22, July 9–13, 2022, Boston, MA, USA

Parsa et al.

Table 1: Simulator Parameters
Parameter

Value

Particle Diameter
Particle Mass
Stiffness Ratio
Packing Fraction
𝑁𝑡

0.1
1
10
0.91
1𝑒4

2.4

Table 2: Optimization Parameters

As was mentioned briefly before, one of the important propertiesof a granular assembly is the existence of gaps in their vibrational
density of states [5]. To find the band gap (the biggest gap), themass-weighted dynamical matrix is calculated using the Hessian of
the total potential energy. The eigenvalues of this matrix are theeigenfrequencies of the system. We can plot the eigenfrequency
spectrum by sorting the frequencies in an increasing order, thenwe see the gaps in the spectrum, where We call the biggest gap
the band gap. When this granular system is excited at a specificfrequency, depending on where that frequency falls within the
spectrum, the acoustic signal will propagate or will be filtered. Thisis the underlying process in the granular metamaterial that enables
it to act as a logic gate.

2.3

Performance Measures3(𝐺 10 (𝜔 2 ) + 𝐺 01 (𝜔 2 ))/2
𝐺 11 (𝜔 2 )Value

Population Size
Generations
Runs
Mutation Probability
Crossover Probability
Bit-flip Probability

50
250
5
0.8
0.2
0.05

RESULTSIn this section, we first present our experimental setup using themethods and metrics developed in the previous section. Then, weshow the optimization results and take a closer look at the Pareto
optimal solutions by investigating their functionality as logic gates.

3.1

Experimental Setup

Our goal is to find a granular assembly made of two types of particles that can act as an AND gate at 𝜔 1 = 7 and as an XOR gate
at 𝜔 2 = 10. These frequencies have been chosen arbitrarily andby looking at the frequency spectrum of random configurations.
We formalized the problem as an optimization problem with two
objectives (equations 2 and 3) and proposed to use the NSGA-II
algorithm to solve the problem. The simulator and optimization
parameters were also presented in the previous section.

(2)

On the other hand, to have an XOR gate at frequency 𝜔 2 , we
expect to have a significant vibration at the output, if only one of
the input ports is being driven by a sinusoidal displacement (‘01‘and ‘10‘ cases). So the fitness in this case can be defined as follows:
𝑓“XOR-ness” =Parameter

The individuals are granular assemblies made of two types ofparticles (soft and stiff particles). We use a direct encoding, so the
genome is a binary string of length 30 (the lattice is 5 × 6). In thisbinary encoding, ‘1‘ means a stiff particle (dark grey particles in
the figures) and ‘0‘ means a soft particle (light grey particles in the
figures). In order to add variation to the population, both mutationand crossover operators are implemented. The mutation operator
is a bit-flip operator with probability of 0.05. The crossover is asingle-point crossover operator. We implement a (𝜆 + 𝜇) evolutionary algorithm where both 𝜆 and 𝜇 are equal to the population size.
Selection of the Pareto nondominated front is based on the Generalized Reduced Run-Time Complexity Non-Dominated Sorting
algorithm presented in [7]. We perform 5 independent runs, eachwith a different random initial population.
All of our code are written in Python. We use the DEAP library[6] for the evolutionary optimization. SCOOP [9] is used to parallelize the code to speed up the population evaluation. All of theexperiments ran on our computing cluster.

The goal of the optimization is to find a configuration of particles
that will act as an AND gate at one frequency (𝜔 1 ) and an XOR
at another frequency (𝜔 2 ). Thus, we need to define metrics foreach of these two cases to measure the amount of “AND-ness”
and “XOR-ness” in a candidate solution. We use the same metrics
that we defined in our previous work [19]. As we discussed at
the start of this section, we can measure the gain of the system(𝐺𝑖 𝑗 (𝜔) in equation 1) for each of the four possible input cases(‘00‘, ‘01‘, ‘10‘, ‘11‘). 𝐺 00 is trivial: if the input is ‘00‘, meaning thatno vibration is applied to either of the input particles, the output
will remain ‘0‘ as well. But for the remaining three cases, we havespecific expectations in a logic gate. In an AND gate, we only wantto have a significant gain when both of the inputs are activatedat 𝜔 1 frequency (high 𝐺 11 (𝜔 1 ) is desired). In order to achieve thisgoal, one option is to define a single fitness function as follows:
𝐺 11 (𝜔 1 )𝑓“AND-ness” =(𝐺 10 (𝜔 1 ) + 𝐺 01 (𝜔 1 ))/2

Evolutionary Search

We employ a standard multiobjective optimization algorithm, Nondominated Sorting Genetic Algorithm II (NSGA-II) [4]. The hyperparameters used in the simulations are shown in Table 2 below:

(3)

3.2

Our goal is to find a configuration of particles that maximizes
both of these objectives, so we can formulate the problem as amultiobjective optimization problem. We present the details of the
optimization in the next section.

Evolved Solutions

Fig. 4 shows the evolutionary progress. Panel (a) of this figure plots
the evolution of the first fitness (“AND-ness” (Eqn. 2)) on top andthe second fitness (“XOR-ness” (Eqn. 3)) on the bottom. Solid blueline in these two plots shows the average fitness of the population,

125

Evolving Programmable Computational Metamaterials

GECCO ’22, July 9–13, 2022, Boston, MA, USA

averaged over 5 independent evolutionary runs, along with its
standard deviation in light blue. We see that the average fitness
fluctuates during the evolution and even decreases at a few points.
Since the problem is a multiobjective optimization, sometimes new
individuals are added to the population that might have lower
values in one fitness but higher values in the other. This can cause
a decrease in the the average fitness of the whole population.
In order to show the effectiveness of the evolutionary search, we
also performed a random search for each of the two objectives: onefor the AND-ness metric and the other for the XOR-ness metric.
Since each one of the 30 particles can either be soft or stiff, there
are a total number of 230 = 1073741824 possible configurations.
In each experiment, we drew 5000 random configurations from a
uniform distribution and evaluated them based on the objective.
Fig. 4c shows the histograms of these two random distributions. For
comparison, we’ve also plotted the best and worst found solutions ineach figure. The histograms show how rugged the fitness landscapeis: for the AND-ness metric, the histogram spreads from 0.026
to 7.88 with a mean value around 0.760 and for the XOR-ness, it
spreads from 0.21 to 60.24 with a mean at 1.622.
Fig. 4b shows the Pareto front at three different stages during theevolution for one of the 5 independent evolutionary runs. We cansee that at the start of the optimization (generation 0), individuals
in the population have objective values close to the mean value ofa randomly drawn solution. As evolution proceeds, the population
progresses towards higher fitness values and also spreads across theaxis to the extreme solutions where one objective function (ANDness/XOR-ness) dominates the other one (see red markers showing
the Pareto front at generation 250). The magenta square on thisplot shows one of the solutions in the middle of the Pareto front
which has moderately high fitness value in both objectives. the plot
on the bottom of this panel shows the chosen particle configuration
along with its fitness values.
All the solutions of the last generation of the evolutionary search
are Pareto-optimal. These solutions are shown in Fig. 3. It’s interesting to investigate these configurations visually, to see if one canfind a regular pattern or a specific order in the particle placementsthat correlates with high fitness values. We’ll talk about this in
more detail in the Discussion section. The magenta star in thisfigure shows a candidate solution with moderately high fitness
values in both objectives. This is the configuration that we choseto investigate further in Fig. 5.
In order to examine the performance of the evolved solution, we
tested the response of the configuration when different vibrations
are applied to the input ports. Fig. 5a shows the functionality of thesolution as an AND gate at 𝜔 1 = 7. Fig. 5b shows its function as an
XOR gate at 𝜔 1 = 10. For each gate, we supply the four different
input signals: ‘00‘, ‘01‘, ‘10‘ and ‘11‘. For the AND gate, we seea high amplitude of oscillation when both of the input ports areactivated at 𝜔 1 = 7 (see (𝑣) in panel (a)). The response is plottedboth in frequency and in time space. In the XOR case in panel (b),only when one of the inputs is activated we see a high magnitude ofdisplacement at the output. The plots confirm the multifunctionalityof the designed granular assembly at different frequencies. In otherwords, we can “program” the material by modulating the frequencyof the input signal. This means that without needing to reconfigure

Figure 3: Configurations in the non-dominated Pareto front
at the last generation of the evolutionary process. The magenta star shows the chosen solution, the green and red stars
show the extreme solutions with highest AND-ness or highest XOR-ness. The numbers beside each configuration show
the fitness values as (AND-ness, XOR-ness). These are the candidate solution of our optimization problem: can you findany specific patterns common in the particle placements
shown in this figure?
the particles or change any material properties, we can observe
different functionalities from the material.

4

DISCUSSION AND CONCLUSION

In our previous paper [19], we introduced the possibility of performing computation in a granular substrate. We showed that we candesign a granular assembly made of two types of particles which
can act as an AND gate or an XOR gate depending on the frequencyof the input vibrations. But one of the reasons that mechanical computing devices were abandoned at the outset of the computer age
and overshadowed by their digital electronic counterparts was theminiaturization capability of electronic transistors and the promiseof the compact computational power predicted by Moore’s Law.
To address this concern in the application of granular metamaterials as a new generation of computational devices, we decided toinvestigate how much computational power can be packed into
one granular configuration. Because of the unique properties ofa granular metamaterial under vibrations with different frequencies, we decided to work on modulating the input frequency as away to program the material to perform different computationaltasks. In the Results section, we demonstrated the success of ourapproach in finding particle configurations that can exhibit two
different functionalities without any changes to the material after
fabrication.
The input frequency is the key to regulate the functionality of thedesigned metamaterial and potentially increase its computational
density. To explore this potential, we investigated the incorporationof a more complex computational block into metamaterials: a half
adder. The half adder is a logical computational block that adds
two input bits and generates a carry (C) and a sum (S) signal, withthe carry signal representing the overflow to the next digit. The

126

GECCO ’22, July 9–13, 2022, Boston, MA, USA

Parsa et al.

Figure 4: Evolutionary progress. The progress of each of the objective functions during the evolution (averaged over 5 independent runs) is plotted in (a) for AND-ness (top) and XOR-ness (bottom). (b) shows the Pareto front of one of the evolutionary
runs at three different stages during the optimization. The magenta square in this plot indicates the chosen solution from thefinal population and its corresponding particle configuration. (c) shows histograms of randomly drawn configurations. They
demonstrate that a random configuration will on average have an AND-ness of 0.760 and XOR-ness of 1.622.
simplest half adder design is made of one AND gate and one XOR
gate and is shown in Fig. 6b. Since the half adder can be implemented
using just AND and XOR gates, it is an interesting test case for ourdesigned metamaterial. As shown in Fig. 6a, we can change the input
vibration to the sum of two sinusoidal waves with two different
frequencies (𝜔 1𝑎𝑛𝑑𝜔 2 ). Since the output of the computational blockis defined as the magnitude of the oscillation at the excitation
frequency (refer to the Methods section) in this system, we will
have two (temporal) outputs: one at 𝜔 1 and the other at 𝜔 2 . The
frequency responses in Fig. 6d demonstrates this idea: there are
two spikes in each plot, corresponding to the system response at
two different frequencies: one at 𝜔 1 = 7 which produces the carrysignal (C) and the other at 𝜔 2 = 10 which produces the sum (S).
One of the important performance metrics in this system is thevariation between the magnitude of vibrations in ‘0‘ and ‘1‘ output
cases: the higher variation results in a more distinguishable on andoff states. We notice that in the third plot of Fig. 6d, when the input
is ‘10‘, the outputs are expected to be ‘0‘ at 𝐶 and ‘1‘ at 𝑆. But the two
spikes at 7 and 10 frequencies are almost of the same magnitude andthe difference is not significant. One way to fix this issue is to choose
different excitation frequencies, meaning that instead of choosing𝜔 1 = 7 and 𝜔 1 = 10 we can have frequencies that are further awayfrom each other in the spectrum. This might reduce the possibility
that the two phases of the input signal interrupt each other and thus
result in a better multifunctional metamaterial. The other option is
to increase the resolution of the granular metamaterial. By adding
more particles in the system, the frequency spectrum of the material

will expand (because the number of degrees of freedom increases).
Therefore we will have more normal modes in the system and havemore freedom in designing the material.
This paper is meant to confirm the computational power of a
granular metamaterial and demonstrate one possible method toprogram such systems. We are not trying to imply that they havehigher performance with regards to traditional digital devices. Much
work remains before a fair comparison of the computational powerof such a system with today’s general purpose digital computers canbe performed as the digital electronics and semiconductor industry
have benefited from billion dollar investments and huge research
efforts. However, it is our hope that this work might help revive
the idea of analogue mechanical computing by highlighting their
computational potential.
In the end, we would also like to point out the potential of granular metamaterials. In this work, we showcased an overly simplified
model made of frictionless circular particles with just one different
material property (stiffness). There are other parameters such as
particle shapes, masses, modulus, friction and other inter-particle
interactions that can be incorporated in building a more sophisticated model. There is also the possibility of utilizing particles thatcan change their physical properties dynamically in response tosome stimuli such as temperature. Such particles can be used tomake a “reconfigurable” granular assembly with the ability to shift
from one configuration to another without the need to reassemble
the whole material. There are many possible ways to reconfigure127Evolving Programmable Computational Metamaterials

GECCO ’22, July 9–13, 2022, Boston, MA, USA

Figure 5: Evaluating the functionality of the candidate solution. (a) examines its performance as an AND gate at 𝜔 = 7. Panel (𝑖)shows the particle configuration on the left and the frequency spectrum on the right. The horizontal lines mark the excitation
frequencies at 𝜔 1 and 𝜔 2 . Panels (𝑖𝑖) through (𝑣) show the response of the system in frequency and time space for the four
input cases. (b) examines the performance of the candidate configuration as an XOR gate at 𝜔 = 10.
metamaterials to exhibit more complex behaviour yet to be explored, some of which may allow these materials to assume muchcomputational power.

5[3] Osama R Bilal, André Foehr, and Chiara Daraio. 2017. Bistable metamaterial forswitching and cascading elastic vibrations. Proceedings of the National Academyof Sciences 114, 18 (2017), 4603–4606.
[4] Kalyanmoy Deb, Amrit Pratap, Sameer Agarwal, and TAMT Meyarivan. 2002. A
fast and elitist multiobjective genetic algorithm: NSGA-II. IEEE transactions on
evolutionary computation 6, 2 (2002), 182–197.
[5] Pierre A Deymier. 2013. Acoustic metamaterials and phononic crystals. Vol. 173.
Springer Science & Business Media.
[6] Félix-Antoine Fortin, François-Michel De Rainville, Marc-André Gardner, Marc
Parizeau, and Christian Gagné. 2012. DEAP: Evolutionary Algorithms Made Easy.
Journal of Machine Learning Research 13 (jul 2012), 2171–2175.
[7] Félix-Antoine Fortin, Simon Grenier, and Marc Parizeau. 2013. Generalizing
the improved run-time complexity algorithm for non-dominated sorting. In
Proceedings of the 15th annual conference on Genetic and evolutionary computation.
615–622.
[8] Tony Hey. 1999. Quantum computing: an introduction. Computing & Control
Engineering Journal 10, 3 (1999), 105–112.
[9] Yannick Hold-Geoffroy, Olivier Gagnon, and Marc Parizeau. 2014. Once you
SCOOP, no need to fork. In Proceedings of the 2014 Annual Conference on Extreme
Science and Engineering Discovery Environment. ACM, 60.

ACKNOWLEDGMENTS

We would like to acknowledge financial support from National
Science Foundation under the DMREF program (award number:
2118810). We also acknowledge computation provided by the Vermont Advanced Computing Core.

REFERENCES[1] Martyn Amos, Gheorghe Păun, Grzegorz Rozenberg, and Arto Salomaa. 2002.
Topics in the theory of DNA computing. Theoretical computer science 287, 1(2002), 3–38.
[2] Robert P Behringer and Bulbul Chakraborty. 2018. The physics of jamming for
granular materials: a review. Reports on Progress in Physics 82, 1 (2018), 012601.

128

GECCO ’22, July 9–13, 2022, Boston, MA, USA

Parsa et al.

Figure 6: Programming the granular metamaterial to act as a half adder. (a) shows the optimal configuration and the input
signals. The material was designed to act as an AND gate at 𝜔 = 7 and as an XOR gate at 𝜔 = 10. The inputs are the sum of two
sinusoidal waves with these two frequencies. (b) reports the logic diagram of a half adder along with its truth table. (c) shows
the time-series response of the programmable granular metamaterial in the four possible cases: sum of two sinusoidal waves
with frequencies 𝜔 1 = 7 and 𝜔 2 = 10 is applied to none of the input ports (‘00‘), input 1 (‘01‘), input 2 (‘10‘) or both (‘11‘). (d)demonstrates the response of the system in the frequency spectrum: the spikes in each signal indicate a ‘1‘ at that frequency.
The material acts as an AND gate at 𝜔 1 , producing the carry signal (C) in the half adder. At 𝜔 2 the material acts as an XOR gate
and so the output port produces the sum signal (S) in the half adder.
[10] Sara Hooker. 2021. The hardware lottery. Commun. ACM 64, 12 (2021), 58–65.
[11] Benjamin Jenett, Christopher Cameron, Filippos Tourlomousis, Alfonso Parra
Rubio, Megan Ochalek, and Neil Gershenfeld. 2020. Discretely assembled mechanical metamaterials. Science advances 6, 47 (2020), eabc9943.
[12] Muamer Kadic, Graeme W Milton, Martin van Hecke, and Martin Wegener. 2019.
3D metamaterials. Nature Reviews Physics 1, 3 (2019), 198–210.
[13] Eunho Kim and Jinkyu Yang. 2019. Wave propagation in granular metamaterials.
Functional Composites and Structures 1, 1 (2019), 012002.
[14] Feng Li, Paul Anzel, Jinkyu Yang, Panayotis G Kevrekidis, and Chiara Daraio.
2014. Granular acoustic switches and logic elements. Nature communications 5,
1 (2014), 1–6.
[15] Marc Z Miskin and Heinrich M Jaeger. 2013. Adapting granular materials throughartificial evolution. Nature materials 12, 4 (2013), 326–331.
[16] Gordon E Moore et al. 1965. Cramming more components onto integrated
circuits.
[17] Kohei Nakajima. 2020. Physical reservoir computing—an introductory perspective. Japanese Journal of Applied Physics 59, 6 (2020), 060501.
[18] Corey S O’Hern and Mark D Shattuck. 2013. Highly evolved grains. Naturematerials 12, 4 (2013), 287–288.
[19] Atoosa Parsa, Dong Wang, Corey S. O’Hern, Mark D. Shattuck, Rebecca KramerBottiglio, and Josh Bongard. 2022. Evolution of Acoustic Logic Gates in Granular
Metamaterials. In International Conference on the Applications of Evolutionary
Computation (Part of EvoStar). Springer. (in press).
[20] Jordan R Raney, Neel Nadkarni, Chiara Daraio, Dennis M Kochmann, Jennifer A
Lewis, and Katia Bertoldi. 2016. Stable propagation of mechanical signals insoft media using stored elastic energy. Proceedings of the National Academy of
Sciences 113, 35 (2016), 9722–9727.
[21] CF Schreck and CS O’Hern. 2010. Computational methods to study jammed
systems. Experimental and computational techniques in soft condensed matterphysics (2010), 25–61.
[22] Marc Serra-Garcia. 2019. Turing-complete mechanical processor via automatednonlinear system design. Physical Review E 100, 4 (2019), 042202.
[23] Alexandre Silva, Francesco Monticone, Giuseppe Castaldi, Vincenzo Galdi, Andrea Alù, and Nader Engheta. 2014. Performing mathematical operations withmetamaterials. Science 343, 6167 (2014), 160–163.
[24] Thomas N Theis and H-S Philip Wong. 2017. The end of moore’s law: A new
beginning for information technology. Computing in Science & Engineering 19, 2(2017), 41–50.
[25] Jiamin Wu, Xing Lin, Yuchen Guo, Junwei Liu, Lu Fang, Shuming Jiao, and Qionghai Dai. 2021. Analog Optical Computing for Artificial Intelligence. Engineering(2021).
[26] Qikai Wu, Chunyang Cui, Thibault Bertrand, Mark D Shattuck, and Corey S
O’Hern. 2019. Active acoustic switches using two-dimensional granular crystals.
Physical Review E 99, 6 (2019), 062901.
[27] Hiromi Yasuda, Philip R Buskohl, Andrew Gillman, Todd D Murphey, Susan
Stepney, Richard A Vaia, and Jordan R Raney. 2021. Mechanical computing.
Nature 598, 7879 (2021), 39–48.
[28] Farzad Zangeneh-Nejad, Dimitrios L Sounas, Andrea Alù, and Romain Fleury.
2021. Analogue computing with metamaterials. Nature Reviews Materials 6, 3
(2021), 207–225.

129