CUPID-Mo Collaboration

††thanks: Now at: Northwestern University, Evanston, IL 60208, USA ††thanks:
Now at: Northwestern University, Evanston, IL 60208, USA

# A new limit for neutrinoless double-beta decay of 100Mo from the CUPID-Mo
experiment

E. Armengaud  IRFU, CEA, Université Paris-Saclay, F-91191 Gif-sur-Yvette,
France  C. Augier  Univ Lyon, Université Lyon 1, CNRS/IN2P3, IP2I-Lyon,
F-69622, Villeurbanne, France  A. S. Barabash  National Research Centre
Kurchatov Institute, Institute of Theoretical and Experimental Physics, 117218
Moscow, Russia  F. Bellini  Dipartimento di Fisica, Sapienza Università di
Roma, P.le Aldo Moro 2, I-00185, Rome, Italy  INFN, Sezione di Roma, P.le Aldo
Moro 2, I-00185, Rome, Italy  G. Benato  INFN, Laboratori Nazionali del Gran
Sasso, I-67100 Assergi (AQ), Italy  A. Benoît  CNRS-Néel, 38042 Grenoble Cedex
9, France  M. Beretta  University of California, Berkeley, California 94720,
USA  L. Bergé  Université Paris-Saclay, CNRS/IN2P3, IJCLab, 91405 Orsay,
France  J. Billard  Univ Lyon, Université Lyon 1, CNRS/IN2P3, IP2I-Lyon,
F-69622, Villeurbanne, France  Yu. A. Borovlev  Nikolaev Institute of
Inorganic Chemistry, 630090 Novosibirsk, Russia  Ch. Bourgeois  Université
Paris-Saclay, CNRS/IN2P3, IJCLab, 91405 Orsay, France  V. B. Brudanin
Laboratory of Nuclear Problems, JINR, 141980 Dubna, Moscow region, Russia  P.
Camus  CNRS-Néel, 38042 Grenoble Cedex 9, France  L. Cardani  INFN, Sezione di
Roma, P.le Aldo Moro 2, I-00185, Rome, Italy  N. Casali  INFN, Sezione di
Roma, P.le Aldo Moro 2, I-00185, Rome, Italy  A. Cazes  Univ Lyon, Université
Lyon 1, CNRS/IN2P3, IP2I-Lyon, F-69622, Villeurbanne, France  M. Chapellier
Université Paris-Saclay, CNRS/IN2P3, IJCLab, 91405 Orsay, France  F. Charlieux
Univ Lyon, Université Lyon 1, CNRS/IN2P3, IP2I-Lyon, F-69622, Villeurbanne,
France  D. Chiesa  Dipartimento di Fisica, Università di Milano-Bicocca,
I-20126 Milano, Italy  INFN, Sezione di Milano-Bicocca, I-20126 Milano, Italy
M. de Combarieu  IRAMIS, CEA, Université Paris-Saclay, F-91191 Gif-sur-Yvette,
France  I. Dafinei  INFN, Sezione di Roma, P.le Aldo Moro 2, I-00185, Rome,
Italy  F. A. Danevich  Institute for Nuclear Research of NASU, 03028 Kyiv,
Ukraine  M. De Jesus  Univ Lyon, Université Lyon 1, CNRS/IN2P3, IP2I-Lyon,
F-69622, Villeurbanne, France  T. Dixon  University of California, Berkeley,
California 94720, USA  L. Dumoulin  Université Paris-Saclay, CNRS/IN2P3,
IJCLab, 91405 Orsay, France  K. Eitel  Karlsruhe Institute of Technology,
Institute for Astroparticle Physics, 76021 Karlsruhe, Germany  F. Ferri  IRFU,
CEA, Université Paris-Saclay, F-91191 Gif-sur-Yvette, France  B. K. Fujikawa
Lawrence Berkeley National Laboratory, Berkeley, California 94720, USA  J.
Gascon  Univ Lyon, Université Lyon 1, CNRS/IN2P3, IP2I-Lyon, F-69622,
Villeurbanne, France  L. Gironi  Dipartimento di Fisica, Università di Milano-
Bicocca, I-20126 Milano, Italy  INFN, Sezione di Milano-Bicocca, I-20126
Milano, Italy  A. Giuliani  e-mail<EMAIL_ADDRESS>Université
Paris-Saclay, CNRS/IN2P3, IJCLab, 91405 Orsay, France  V. D. Grigorieva
Nikolaev Institute of Inorganic Chemistry, 630090 Novosibirsk, Russia  M. Gros
IRFU, CEA, Université Paris-Saclay, F-91191 Gif-sur-Yvette, France  E. Guerard
Université Paris-Saclay, CNRS/IN2P3, IJCLab, 91405 Orsay, France  D. L. Helis
IRFU, CEA, Université Paris-Saclay, F-91191 Gif-sur-Yvette, France  H. Z.
Huang  Key Laboratory of Nuclear Physics and Ion-beam Application (MOE), Fudan
University, Shanghai 200433, PR China  R. Huang  University of California,
Berkeley, California 94720, USA  J. Johnston  Massachusetts Institute of
Technology, Cambridge, MA 02139, USA  A. Juillard  Univ Lyon, Université Lyon
1, CNRS/IN2P3, IP2I-Lyon, F-69622, Villeurbanne, France  H. Khalife
Université Paris-Saclay, CNRS/IN2P3, IJCLab, 91405 Orsay, France  M. Kleifges
Karlsruhe Institute of Technology, Institute for Data Processing and
Electronics, 76021 Karlsruhe, Germany  V. V. Kobychev  Institute for Nuclear
Research of NASU, 03028 Kyiv, Ukraine  Yu. G. Kolomensky  University of
California, Berkeley, California 94720, USA  Lawrence Berkeley National
Laboratory, Berkeley, California 94720, USA  S.I. Konovalov  National Research
Centre Kurchatov Institute, Institute of Theoretical and Experimental Physics,
117218 Moscow, Russia  A. Leder  Massachusetts Institute of Technology,
Cambridge, MA 02139, USA  P. Loaiza  Université Paris-Saclay, CNRS/IN2P3,
IJCLab, 91405 Orsay, France  L. Ma  Key Laboratory of Nuclear Physics and Ion-
beam Application (MOE), Fudan University, Shanghai 200433, PR China  E. P.
Makarov  Nikolaev Institute of Inorganic Chemistry, 630090 Novosibirsk, Russia
P. de Marcillac  Université Paris-Saclay, CNRS/IN2P3, IJCLab, 91405 Orsay,
France  R. Mariam  Université Paris-Saclay, CNRS/IN2P3, IJCLab, 91405 Orsay,
France  L. Marini  University of California, Berkeley, California 94720, USA
Lawrence Berkeley National Laboratory, Berkeley, California 94720, USA  INFN,
Laboratori Nazionali del Gran Sasso, I-67100 Assergi (AQ), Italy  S. Marnieros
Université Paris-Saclay, CNRS/IN2P3, IJCLab, 91405 Orsay, France  D. Misiak
Univ Lyon, Université Lyon 1, CNRS/IN2P3, IP2I-Lyon, F-69622, Villeurbanne,
France  X.-F. Navick  IRFU, CEA, Université Paris-Saclay, F-91191 Gif-sur-
Yvette, France  C. Nones  IRFU, CEA, Université Paris-Saclay, F-91191 Gif-sur-
Yvette, France  E. B. Norman  University of California, Berkeley, California
94720, USA  V. Novati  Université Paris-Saclay, CNRS/IN2P3, IJCLab, 91405
Orsay, France  E. Olivieri  Université Paris-Saclay, CNRS/IN2P3, IJCLab, 91405
Orsay, France  J. L. Ouellet  Massachusetts Institute of Technology,
Cambridge, MA 02139, USA  L. Pagnanini  INFN, Gran Sasso Science Institute,
I-67100 L’Aquila, Italy  INFN, Laboratori Nazionali del Gran Sasso, I-67100
Assergi (AQ), Italy  P. Pari  IRAMIS, CEA, Université Paris-Saclay, F-91191
Gif-sur-Yvette, France  L. Pattavina  INFN, Laboratori Nazionali del Gran
Sasso, I-67100 Assergi (AQ), Italy  Physik Department, Technische Universität
München, Garching D-85748, Germany  B. Paul  IRFU, CEA, Université Paris-
Saclay, F-91191 Gif-sur-Yvette, France  M. Pavan  Dipartimento di Fisica,
Università di Milano-Bicocca, I-20126 Milano, Italy  INFN, Sezione di Milano-
Bicocca, I-20126 Milano, Italy  H. Peng  Department of Modern Physics,
University of Science and Technology of China, Hefei 230027, PR China  G.
Pessina  INFN, Sezione di Milano-Bicocca, I-20126 Milano, Italy  S. Pirro
INFN, Laboratori Nazionali del Gran Sasso, I-67100 Assergi (AQ), Italy  D. V.
Poda  Université Paris-Saclay, CNRS/IN2P3, IJCLab, 91405 Orsay, France  O. G.
Polischuk  Institute for Nuclear Research of NASU, 03028 Kyiv, Ukraine  S.
Pozzi  INFN, Sezione di Milano-Bicocca, I-20126 Milano, Italy  E. Previtali
Dipartimento di Fisica, Università di Milano-Bicocca, I-20126 Milano, Italy
INFN, Sezione di Milano-Bicocca, I-20126 Milano, Italy  Th. Redon  Université
Paris-Saclay, CNRS/IN2P3, IJCLab, 91405 Orsay, France  A. Rojas  LSM,
Laboratoire Souterrain de Modane, 73500 Modane, France  S. Rozov  Laboratory
of Nuclear Problems, JINR, 141980 Dubna, Moscow region, Russia  C. Rusconi
Department of Physics and Astronomy, University of South Carolina, SC 29208,
Columbia, USA  V. Sanglard  Univ Lyon, Université Lyon 1, CNRS/IN2P3,
IP2I-Lyon, F-69622, Villeurbanne, France  J. A. Scarpaci  Université Paris-
Saclay, CNRS/IN2P3, IJCLab, 91405 Orsay, France  K. Schäffner  INFN,
Laboratori Nazionali del Gran Sasso, I-67100 Assergi (AQ), Italy  B. Schmidt
Lawrence Berkeley National Laboratory, Berkeley, California 94720, USA  Y.
Shen  Key Laboratory of Nuclear Physics and Ion-beam Application (MOE), Fudan
University, Shanghai 200433, PR China  V. N. Shlegel  Nikolaev Institute of
Inorganic Chemistry, 630090 Novosibirsk, Russia  B. Siebenborn  Karlsruhe
Institute of Technology, Institute for Astroparticle Physics, 76021 Karlsruhe,
Germany  V. Singh  University of California, Berkeley, California 94720, USA
C. Tomei  INFN, Sezione di Roma, P.le Aldo Moro 2, I-00185, Rome, Italy  V. I.
Tretyak  Institute for Nuclear Research of NASU, 03028 Kyiv, Ukraine  V. I.
Umatov  National Research Centre Kurchatov Institute, Institute of Theoretical
and Experimental Physics, 117218 Moscow, Russia  L. Vagneron  Univ Lyon,
Université Lyon 1, CNRS/IN2P3, IP2I-Lyon, F-69622, Villeurbanne, France  M.
Velázquez  Université Grenoble Alpes, CNRS, Grenoble INP, SIMAP, 38402 Saint
Martin d’Héres, France  B. Welliver  Lawrence Berkeley National Laboratory,
Berkeley, California 94720, USA  L. Winslow  Massachusetts Institute of
Technology, Cambridge, MA 02139, USA  M. Xue  Department of Modern Physics,
University of Science and Technology of China, Hefei 230027, PR China  E.
Yakushev  Laboratory of Nuclear Problems, JINR, 141980 Dubna, Moscow region,
Russia  M. Zarytskyy  Institute for Nuclear Research of NASU, 03028 Kyiv,
Ukraine  A. S. Zolotarova  Université Paris-Saclay, CNRS/IN2P3, IJCLab, 91405
Orsay, France

###### Abstract

The CUPID-Mo experiment at the Laboratoire Souterrain de Modane (France) is a
demonstrator for CUPID, the next-generation ton-scale cryogenic
$0\nu\beta\beta$ experiment. It consists of a 4.2 kg array of 20 enriched
Li2100MoO4 scintillating bolometers to search for the lepton number violating
process of $0\nu\beta\beta$ decay in 100Mo. With more than one year of
operation (2.16 kg$\times$yr of physics data), no event in the region of
interest and hence no evidence for $0\nu\beta\beta$ is observed. We report a
new limit on the half-life of $0\nu\beta\beta$ decay in 100Mo of
$T_{1/2}>1.5\times 10^{24}\,$yr at 90% C.I. The limit corresponds to an
effective Majorana neutrino mass $\langle m_{\beta\beta}\rangle$ $<$ (0.3–0.5)
eV, dependent on the nuclear matrix element in the light Majorana neutrino
exchange interpretation.

††preprint: APS/123-QED

The discovery that neutrinos are massive particles through the evidence of
neutrino flavor oscillations [1] opens the question of neutrino mass
generation. Instead of having Dirac nature as charged leptons and quarks, the
scale of neutrino masses could be well motivated by the Majorana theory [2,
3]. In this scenario neutrinos could coincide with their antimatter partner
[4, 5] which would have a tremendous impact on our vision of Nature, implying
the violation of the total lepton number $L$ as well as for the matter-
antimatter asymmetry in the Universe [6], see also review [7].

The distinction between Dirac and Majorana behavior is an extreme experimental
challenge. Neutrinoless double-beta ($0\nu\beta\beta$) decay is the
traditional and the most sensitive tool to probe the Majorana nature of
neutrinos. This process is a nuclear transition consisting in the
transformation of an even-even nucleus into a lighter isobar containing two
more protons and accompanied by the emission of two electrons and no other
particles, with a change of the lepton number $L$ by two units [8, 9, 10, 11].
An observation of this rare process would establish that neutrinos are
Majorana particles [12]. The current most stringent limits on $0\nu\beta\beta$
decay half-lifes are at the level of $10^{26}$ yr in 136Xe and 76Ge [13, 14].

$0\nu\beta\beta$ decay can be induced by a variety of mechanisms [11, 9, 15,
16]. Among them, the so-called mass mechanism — consisting in the exchange of
a virtual light Majorana neutrino — represents a minimal extension of the
Standard Model. In this mechanism, the $0\nu\beta\beta$ decay rate is
proportional to the square of the effective Majorana neutrino mass $\langle
m_{\beta\beta}\rangle$, a linear combination of the three neutrino mass
eigenvalues which fixes the absolute neutrino mass scale. Present limits on
$\langle m_{\beta\beta}\rangle$ are in the range of $(0.06-0.6)$ eV [11],
assuming that the axial charge $g_{A}$ is not quenched and equal to the free
nucleon value of $\simeq 1.27$ [17, 18, 19].

The distinctive signal of $0\nu\beta\beta$ decay is a peak in the total energy
spectrum of the two emitted electrons at the total available energy
$Q_{\beta\beta}$ of the $0\nu\beta\beta$ transition, enlarged only by the
finite detector energy resolution. Besides $0\nu\beta\beta$ decay, candidate
nuclei undergo standard two-neutrino double-beta ($2\nu\beta\beta$) decay
where two anti-neutrinos are emitted along with the two electrons [20].
$2\nu\beta\beta$ conserves the lepton number $L$ and as anti-neutrinos escape
detection, the two-electron energy spectrum is a continuum extending from zero
to $Q_{\beta\beta}$. This process has been observed in several nuclei with
half-lives in the 1018–1024 yr range [21, 22]. Among the 35 natural double-
beta emitters [23], only a few of them are experimentally relevant because of
their high $Q_{\beta\beta}$ ($>2$ MeV) with high decay probability, typically
favorable backgrounds at higher energy and their high natural isotopic
abundance and/or technically feasible isotopic enrichment.

Low-temperature calorimeters, often named bolometers, are the detectors of
choice for several experimental efforts, including the one reported here.
Featuring high energy resolution, high efficiency, and flexibility in
detector-material choice [24, 25, 26], bolometers are perfectly tailored to
$0\nu\beta\beta$ search. These detectors consist in a single crystal that
contains the $0\nu\beta\beta$ source coupled to a temperature sensor. The
signal is collected at very low temperatures $\lesssim 20$ mK for large (0.1–1
kg) bolometers and consists of a thermal pulse registered by the sensor.

A detector embedding a candidate with ${Q_{\beta\beta}>2615\,}$keV is an
optimal choice in terms of background control, as the bulk of the $\gamma$
natural radioactivity ends at 2615 keV, corresponding to the energy of the
${}^{208}\mathrm{Tl}$ line in the 232Th decay chain. However, the energy
region above $\sim 2.6\,$MeV is dominated by events due to surface radioactive
contamination, especially energy-degraded $\alpha$ particles [27, 28], as
shown by the results of CUORE, the largest $0\nu\beta\beta$ bolometric
experiment currently under way.

A method to reject $\alpha$ events in bolometers exploits light emission —
scintillation or Cherenkov — from certain detector materials [29, 26, 30, 25,
31, 32, 33, 34]. Since the $\gamma$ and $\beta$ light yield is appreciably
different from the $\alpha$ light yield at equal deposited energy, while the
bolometric thermal response is almost equivalent, the simultaneous detection
of light and heat can effectively reject the surface $\alpha$ background.

The detector technology of CUPID-Mo has been developed by the LUMINEU
Collaboration [35, 36] and is described together with the experimental setup
in [37]. The isotope of interest 100Mo features a $Q_{\beta\beta}$ of
$(3034.40\pm 0.17)$ keV [38] and a natural abundance of 9.7% making large
scale enrichment viable by gas centrifuge isotopic separation. In CUPID-Mo, it
is embedded into scintillating Li2100MoO4 (LMO) crystals by a double low-
thermal-gradient Czochralski crystallization process [39] from enriched Mo
previously used in the NEMO-3 experiment [40]. A total of twenty cylindrical
$\sim 210\,$g crystals are stacked into 5 towers which results in a 100Mo mass
of $(2.258\pm 0.005)\,$kg with an average 100Mo isotopic abundance of
$(96.6\pm 0.2)\%$. Bolometric Ge light detectors (LDs) complement each LMO
detector module, such that all but the top LMO elements of the 4-crystal tower
see a LD at both the top and bottom. The LMO crystals as well as the LDs are
instrumented with Neutron-Transmutation-Doped (NTD)-Ge sensors [41]. The
towers are installed with a mechanical decoupling inside the EDELWEISS
cryogenic infrastructure [42, 43] at the Laboratoire Souterrain de Modane in
France.

The data of the present analysis have been acquired over a 380 day period
between March 2019 and April 2020 at operation temperatures of 20.7 and 22 mK.
$82\,\%$ of the time was devoted to the $0\nu\beta\beta$ search, split into
240 days of physics data and 73 days of calibration data.

We only consider 213 out of the 240 days of physics data, split into 7
datasets for which we also have calibration data with high statistics and
hence a precise knowledge of the energy scale.

We further exclude periods of temperature instabilities, disturbances in the
underground laboratory and periods of excessive noise on the individual
detectors reducing the physics exposure by 6%. We reject one of the twenty LMO
bolometers that shows an abnormal performance [37] and obtain a physics
exposure of 2.16 kg$\times$yr (Li2100MoO4).

All data are acquired as a continuous stream with $500\,$Hz sampling frequency
and analyzed with a software package developed by the CUORE [44] and CUPID-0
[45] Collaborations, first used in CUPID-Mo in [37, 46]. We use an optimum
filter based analysis chain with 3 s long pulse traces for both the LMO and LD
channels. The data were re-triggered offline using the optimum filter [47,
28]. This resulted in significantly lower trigger thresholds with 90%
efficiency values attained at typical (median) energies of 9.4 keV / 0.5 keV
for the LMOs / LDs. The LMOs analysis and coincidence thresholds have been set
at 45 keV, well above this efficiency turn on. For each signal on an LMO
detector we evaluate the resolution weighted average light signal of the two
(one) adjacent LDs to discriminate $\alpha$ events, exhibiting $\sim$20% of
the light yield of $\gamma$&$\beta$ events of the same energy [37].
Furthermore, the energy difference of the upper and lower LDs, corrected for
geometrical light collection efficiencies is computed to extend the
identification of near surface decays. We developed a new pulse shape analysis
procedure for the signal of the LMO detector to remove remaining spikes and
pulser-induced cross-talk events from the data. As described in detail in [48]
a principal component analysis (PCA) is used to establish a characteristic
shape of a pulse as the leading PCA component(s). We then compute the
reconstruction error defined as the square-root of the variance from a pulse
with respect to its projection onto the leading PCA component(s). The
reconstruction error shows a linear energy dependence and a constant spread
over our $\gamma$&$\beta$ data which makes it ideally suitable for the
$0\nu\beta\beta$ analysis where it needs to be extrapolated to
$Q_{\beta\beta}$ of 100Mo at $3034\,$keV. Together with this very simple
energy dependence we gain a $\sim 9\%$ efficiency increase compared to the
analysis presented in [46]. Lastly we calibrate the response of the LMO
detectors with a 2nd order polynomial using the U/Th calibration data as in
[46] and cross-calibrate the LD against the LMO signals.

We adopt a blinding strategy removing all events in a $\pm 50\,$keV window
around $Q_{\beta\beta}$ to avoid any bias in the optimization of our analysis
procedures and consider the following event selections. For events

  * •

to be contained in a single crystal and in anti-coincidence with a triple
module trigger and energy deposit in the muon-veto system [49] based on a 100
ms time window;

  * •

to have a single trigger in each 3-s pulse window;

  * •

to have a flat pre-trace with a slope of less than 15 median absolute
deviations;

  * •

to have a pulse shape compatible with the principal components established by
the PCA. This cut is optimized using calibration data by maximizing a
hypothetical discovery sensitivity for a $0\nu\beta\beta$ process equal in
half-live to the previous best limit [40], see also [48];

  * •

to have the expected light yield for $\gamma$&$\beta$ events and no difference
in top and bottom LDs. Both of these cuts are set to obtain close to full
coverage at $\pm 3\,\sigma$ each, based on a Gaussian fit of the light yield
in calibration data. The energy dependence of the cut was modeled with a
linear function after we observed an excess broadening of the recorded light
yield with respect to the photon statistics model discussed in [37]. The
excess width is associated with an under-sampling of the faster LD pulses and
is presently under further investigation.

Figure 1: Physics spectrum (blue) for 2.16 kg$\times$yr of data and
calibration spectrum (red) scaled to match the 2615 keV counts from
${}^{208}\mathrm{Tl}$. A $\pm\,50\,$keV region around $Q_{\beta\beta}$ has
been blinded (gray). The most prominent $\gamma$-peaks in the physics and
calibration data are the following: (1) ${}^{208}\mathrm{Tl}$ 2615 keV, (2)
${}^{40}\mathrm{K}$ 1461 keV, (3) ${}^{60}\mathrm{Co}$ 1173 and 1332 keV, (4)
${}^{214}\mathrm{Bi}$ 1764, 1120 and 609 keV, (5) ${}^{212}\mathrm{Pb}$ 239
keV, (6) 99Mo 142 keV.

The resulting spectra summed over 19 LMO detectors and the entire data taking
period are shown in Fig. 1.

The physics data exhibit an extended $2\nu\beta\beta$ component, associated
with the short $2\nu\beta\beta$ half-life (high rate) of 100Mo [50, 51], which
dominates the spectrum from 0.5 to 3 MeV.  Above this component only a limited
number of $\gamma$ peaks from remaining natural radioactivity in the
experimental setup are visible. The most prominent $\gamma$ peaks in our data
are ${}^{208}\mathrm{Tl}$, ${}^{40}\mathrm{K}$, ${}^{60}\mathrm{Co}$ and an
activation $\gamma$ peak from 99Mo, visible for a short time after a neutron
irradiation of the detectors. For more details we refer to a prior
characterization of the backgrounds in this setup by the EDELWEISS
Collaboration in [42]. The four labeled $\gamma$ peaks in the scaled
calibration spectrum in Fig. 1 are used for the calibration of each data-set.

We optimize the $0\nu\beta\beta$ search window for a Poisson counting process
in the low background regime. We consider the current exposure of 2.16
kg$\times$yr as a first unblinding of the full data of the now completed
CUPID-Mo experiment. For the definition of the signal region of interest (ROI)
we hence optimize using the full expected exposure and consider detector and
data-set based resolutions and a preliminary background index evaluated from
the current data release.

The detector resolutions are extracted from a data-set based, simultaneous
unbinned extended maximum likelihood fit of the ${}^{208}\mathrm{Tl}$ 2615 keV
peak in calibration data. The fit is performed with individual parameters for
the detector resolutions, peak amplitudes and position and with common
parameters for the peak-background ratio [46]. Limited by the statistics in
physics data we estimate a global scaling factor (common to all data-sets and
detectors) that is used to project the individual calibration resolutions to
$Q_{\beta\beta}$. We report the more conservative value of an extrapolation
considering the preferred scaling described in [46] as well as the polynomial
fit used in [28, 27]. This results in a characteristic resolution of $(7.6\pm
0.7)\,\mathrm{keV}$ FWHM at $Q_{\beta\beta}$ for the overall data taking. We
note that the estimation of the scaling factor does have some sensitivity to
non-Gaussian tails in calibration data, leading to an excess systematic error
of 2%.

The accuracy of the energy scale is further investigated based on the mean
position of the 352, 583, 609, 1461 and 2615 keV peaks in our physics spectrum
(Fig. 1). A 2nd order polynomial fit in reconstructed-to-expected peak
position results in an energy bias at $Q_{\beta\beta}$ compatible with zero:
${E_{bias}^{Q_{\beta\beta}}=(-0.2\pm 0.4)\,\mathrm{keV}}$. We note that we
excluded the prominent peaks from ${}^{60}\mathrm{Co}$, due to a known spatial
non-uniformity, and 99Mo, due to its time dependence.

The background index has been evaluated from the still blinded data with a
phenomenological fit model that contains an exponential to approximate both
the high energy part of the $2\nu\beta\beta$ spectrum as well as tails from
U/Th contaminants in the setup, and a constant as a conservative estimate for
the coincident detection of two $2\nu\beta\beta$ events in the same crystal,
remaining un-vetoed muon events and close contamination from the high energy
beta decays in the natural U/Th chains. The result of an unbinned extended
maximum likelihood fit is strongly dependent on the initial and final fit
range.  For a fit with the start range varied from 2.65–2.9 MeV and the final
fit range from the upper end of the blinded region to 4 MeV we obtain a
background index of $2\times 10^{-3}\,$counts/(keV$\times$kg$\times$yr) to
$6\times 10^{-3}\,$counts/(keV$\times$kg$\times$yr) in a 10 keV window around
$Q_{\beta\beta}$.  This background estimate is consistent with the results
from a preliminary Geant4-based reconstruction of the $\gamma$&$\beta$
backgrounds, similar to the ones described in [52, 53, 32].

Considering the large remaining uncertainty we round the background index for
the ROI optimization to ${b=5\times
10^{-3}\,}$counts/(keV$\times$kg$\times$yr). We model the background as
locally flat, consider detector and data-set based resolutions, and simulate
the $0\nu\beta\beta$ peak containment in our Geant4 Monte Carlo model. We
maximize the mean limit setting sensitivity of a Poisson process

$\overline{S_{90}}=\sum_{i=0}^{\infty}P(i,b)\cdot S_{90}(i)$

with the sum running over the product of the Poisson probability $P(i,b)$ of
obtaining $i$ background events for a background index $b$ times the expected
90% confidence exclusion limit $S_{90}(i)$. We transfer this maximization from
the optimization of the energy range for a peak search in 19 (detectors) times
7 (data-sets) to the optimization of a single parameter by splitting the
simulated smeared $0\nu\beta\beta$ peaks into 0.1 keV bins and ranking each
bin associated with a triplet (detector, data-set, energy) in Signal-
Background (B/S) likelihood space. The optimal cutoff parameter (B/S)cutoff
results in a ROI that is on average (exposure weighted) 17.9 keV wide. It has
a mean signal containment of $75.8\%$ with a spread of $\pm 1.0\%$. The ROI
width corresponds to an average 2.7$\,\sigma$ Gaussian coverage with the loss
of $0\nu\beta\beta$ decay events in the full energy peak dominated by
Bremsstrahlung loss from events close to the surface of the crystals. The
optimization exhibits only a mild dependence on the background index or the
knowledge of the resolution, with the overall containment changing by $\pm
0.7\%$ for a $50\%$ change in $b$ ($2.2\,$keV wider, $1.5\,$keV narrower ROI).
We truncate the computation of the mean limit setting sensitivity after the
first three terms as the probability of 3 or more background events is
negligible for the considered ROI.

We used two independent analyses to extract a limit or potential signal on the
rate of $0\nu\beta\beta$ events. We present results from a likelihood analysis
of the 100 keV wide blinded analysis region. This analysis has been
implemented in the framework of the Bayesian Analysis Toolkit (BAT) [54], but
we note that we obtained compatible results from a Poisson counting analysis.
In the BAT analysis we consider both the signal region as well as the
sidebands in our likelihood function

$L=\prod_{i=1}^{3}\frac{e^{\lambda_{i}}\lambda_{i}^{n_{i}}}{n_{i}!},$

which is the product over three Poisson terms with observed events $n_{i}$ and
expected events $\lambda_{i}$. The mean number of expected events
$\lambda_{i}$ is computed considering the phenomenological background model
described above and a Gaussian signal contribution in which we leave the
strength of the signal and flat background component free by using
uninformative flat priors. After defining all analysis steps we unblind and
obtain the spectrum in Fig. 2. We observe no event in the signal region and a
single event (cyan) in the right-hand side region. The corresponding
marginalized posterior distribution for the number of signal events has a most
probable value of zero with an upper limit of 2.4 events at 90% C.I.,
resulting in a half-life limit for $0\nu\beta\beta$ decay in 100Mo of
$T_{1/2}^{0\nu}>1.4\times 10^{24}\,\mathrm{yr}\quad\mathrm{(90\%~{}C.I.)}.$
The posterior for the flat background is non-zero with a 1$\,\sigma$ interval
of $3\begin{subarray}{c}+7\\\ -3\end{subarray}\times
10^{-3}\,$counts/(keV$\times$kg$\times$yr) and the posterior distributions for
the parameters of the exponential are compatible with priors from a fit of the
$2\nu\beta\beta$ spectrum in the 2650–2980 keV interval. We repeat the same
fit for the approximation of a Gaussian signal with a locally flat background
over the 100 keV analysis region. The limit on $0\nu\beta\beta$ decay of 100Mo
is unaffected.

Figure 2: Physics spectrum for 2.16 kg$\times$yr of data after unblinding. No event is observed in the detector and dataset based ROI. A single event, highlighted in cyan has been observed in the analysis region. In a further refinement of the analysis it was identified as a $\beta$-candidate out of the 212Bi$\rightarrow^{208}$Tl$\rightarrow^{208}$Pb part of the natural decay chain (see text). For visualization the exposure weighted mean ROI for $0\nu\beta\beta$ decay (17.9 keV wide) has been indicated with solid black lines.  Table 1: Nuisance parameters included in the analysis and their implementation with Flat or Gaussian prior in the Bayesian fit. Parameters 2 and 5 are multiplicative scaling factors instead of absolute uncertainties, see text for details.  Systematic | Index | Value | Prior
---|---|---|---
Isotopic enrichment | 1 | $0.966\pm 0.002$ | Gauss
$0\nu\beta\beta$ containment MC | 2 | $1.000\pm 0.015$ | Gauss
$0\nu\beta\beta$ cont. detector response | 3 | 0.95–1.00 | Flat
Analysis efficiency a | 4 | $0.906\pm 0.004$ | Gauss
Light yield selection 111Data-set dependent; exposure weighted mean value presented. | 5 | 0.998–1.008 | Flat

Nuisance parameters are summarized in Table 1. Uncertainties on the detector
response, in particular the energy scale and resolutions, are included in the
simulation of $0\nu\beta\beta$ events. They are hence covered in the resulting
containment in the optimized central ROI on a detector and data-set basis and
not considered independently. The only remaining uncertainty for the detector
response (Index 3) is based on a potential non-gaussianity of the
$0\nu\beta\beta$ peak. In this analysis we estimate this contribution based on
the shape of the 2615 keV calibration peak. We observe evidence for non-
gaussian tails, which are dominated by unrejected pile-up events caused by the
high trigger-rate in calibration data. Given the much lower trigger-rate in
physics data we believe that the containment reduction of up to 5%, obtained
for pulse shape cuts that are representative of the ones used in physics data,
is a very conservative estimate.

The second nuisance parameter on the containment (Index 2) accounts for the
Geant4 modeling uncertainty of Bremsstrahlung events. Reported accuracies for
the Geant4 Bremsstrahlung production of a few MeV electrons in thick targets
[55, 56] of $\sim$10% result in a systematic uncertainty in the containment of
$1.5\%$ for our crystal geometry, which we include as a multiplicative factor
with a Gaussian prior.

The inclusion of the analysis efficiency

$\epsilon=(90.6\pm 0.4\,\mathrm{(stat.)}\,\begin{subarray}{c}+0.8\\\
-0.2\end{subarray}\,\mathrm{(syst.)})\%$

is split into two parts. For the evaluation of the mean value and its
statistical uncertainty (Index 4), we make use of the two independent signals
in the LDs and LMOs to evaluate cut efficiencies on a clean sample of signal
events in the 1.3 MeV to 2 MeV $2\nu\beta\beta$ spectrum or from the
${}^{210}\mathrm{Po}$ peak [46]. Energy independent cuts are evaluated
directly from the ratio between passed and total events with binomial
uncertainty. The pulse shape analysis efficiency is extracted from a linear
fit extrapolated to $Q_{\beta\beta}$ in order to account for the energy
dependence in the reconstruction error.

The systematic uncertainty associated with the excess broadening of the light
yield cut has been evaluated considering either the applied linear cut and
linear energy dependence or a modified statistical behavior with arbitrary
$a_{0}$ where the width of the light yield band is
$\sigma(E)=\sqrt{\sigma_{0}^{2}+(a_{0}\sqrt{E})^{2}}$. The quantification with
a set of pseudo-experiments results in a likely underestimate of the signal
acceptance which is reflected in our limit setting as a multiplicative factor
with uniform prior in 0.998–1.008 (Index 5). Lastly, we include the
subdominant uncertainty in the number of 100Mo atoms of 0.2% (Index 1). All of
these nuisance parameters are considered in the limit of our blind analysis
${T_{1/2}^{0\nu}>1.4\times 10^{24}\,\mathrm{yr}~{}\mathrm{(90\%~{}C.I.)}}$.

We further refined our analysis after unblinding, implementing a cut designed
to reject high energy $\beta$ events from the
212Bi$\xrightarrow{\alpha}$${}^{208}\mathrm{Tl}$ $\xrightarrow{\beta}$ 208Pb
branch in the thorium chain ($T_{1/2}^{{}^{208}\textrm{Tl}}=183\,$s, 5 MeV
$Q$-value). Similar to previous analyses with scintillating bolometers [29,
45, 32] we tag 212Bi $\alpha$ candidates with energies in the 6.0–6.3 MeV
range, and veto any decay in the same crystal in a 10 half-life period (1832
s). This cut has a negligible impact on the life-time (0.02%) accidentally
rejecting 2:10000 events, while being fully efficient in rejecting any tagged
true 212Bi$\xrightarrow{\alpha}$${}^{208}\mathrm{Tl}$ $\xrightarrow{\beta}$
208Pb event. It rejects the event highlighted in cyan in Fig. 2. The energy of
the preceding $\alpha$ candidate is consistent with the $Q$-value of 212Bi
within 10 keV and the time difference between the events is 113 s.

We report a final $0\nu\beta\beta$ limit that is 1.3% stronger and rounds to

${T_{1/2}^{0\nu}>1.5\times 10^{24}\,\mathrm{yr}~{}\mathrm{(90\%~{}C.I.)}}.$

The posterior for the flat background of the bayesian fit in this case is
peaked at zero with a 90% C.I. of $1.1\times
10^{-2}\,$counts/(keV$\times$kg$\times$yr).

We interpret the obtained half-life limit in the framework of light Majorana
neutrino exchange using ${g_{A}=1.27}$, phase space factors from [57, 58] and
nuclear matrix element calculations from [59, 60, 61, 62, 63, 64, 65, 66]. The
resulting limit on the effective Majorana neutrino mass of $\langle
m_{\beta\beta}\rangle<$(0.3–0.5) eV is the fourth most stringent limit world
wide, obtained with a modest 100Mo exposure of 1.19 kg$\times$yr. It is the
leading constraint for 100Mo, exceeding the previous best limit from NEMO-3
[40] by $30\%$ with almost 30 times lower 100Mo exposure. The technology of
CUPID-Mo has proven that it can be operated reliably, reaches high efficiency
for $0\nu\beta\beta$ search of $68.6\%$ (containment $\times$ analysis
efficiency) and a resolution of ${0.11\%~{}(1\sigma)}$ at $Q_{\beta\beta}$.
Together with a preliminary estimate of the background in the ROI at the few
$10^{-3}\,$counts/(keV$\times$kg$\times$yr) level, this results in an
excellent sensitivity to $0\nu\beta\beta$ decay for the present and future
experiments. Further analyses will be focused on precisely reconstructing
remaining backgrounds, comparing to the best reported background index for a
cryogenic $0\nu\beta\beta$ search [32] and to optimally design and use the
technology of the CUPID-Mo experiment in CUPID [67].

The help of the technical staff of the Laboratoire Souterrain de Modane and of
the other participant laboratories is gratefully acknowledged. This work has
been partially performed in the framework of the LUMINEU program, funded by
ANR (France). This work was also supported by CEA and IN2P3, by P2IO LabEx
with the BSM-Nu project, by Istituto Nazionale di Fisica Nucleare, the Russian
Science Foundation (Russia), the National Research Foundation (Ukraine), the
US Department of Energy (DOE) Office of Science, the DOE Office of Science,
Office of Nuclear Physics, the National Science Foundation (USA), by the
France-Berkeley fund, the MISTI-France fund, and by the Chateaubriand
Fellowship of the Office for Science & Technology of the Embassy of France in
the United States. Individuals have received support from the National Academy
of Sciences of Ukraine and P2IO LabEx managed by the ANR (France). This work
makes use of the DIANA data analysis software which has been developed by the
Cuoricino, CUORE, LUCIFER, and CUPID-0 Collaborations.

## References

  * Esteban  _et al._ [2019] I. Esteban, M. C. Gonzalez-Garcia, A. Hernandez-Cabezudo, M. Maltoni, and T. Schwetz, J. High Energy Phys. 01, 106 (2019).
  * Mohapatra and Senjanović [1980] R. N. Mohapatra and G. Senjanović, Phys. Rev. Lett. 44, 912 (1980).
  * Schechter and Valle [1980] J. Schechter and J. W. F. Valle, Phys. Rev. D 22, 2227 (1980).
  * Majorana [1937] E. Majorana, Nuovo Cim. 14, 171 (1937).
  * Racah [1937] G. Racah, Nuovo Cim. 14, 322 (1937).
  * Fukugita and Yanagida [1986] M. Fukugita and T. Yanagida, Phys. Lett. B 174, 45 (1986).
  * Davidson  _et al._ [2008] S. Davidson, E. Nardi, and Y. Nir, Phys. Rept. 466, 105 (2008).
  * Furry [1939] W. H. Furry, Phys. Rev. 56, 1184 (1939).
  * Bilenky and Giunti [2015] S. M. Bilenky and C. Giunti, Int. J. Mod. Phys. A 30, 1530001 (2015).
  * Dell’Oro  _et al._ [2016] S. Dell’Oro, S. Marcocci, M. Viel, and F. Vissani, Adv. High Energy Phys. 2016, 2162659 (2016).
  * Dolinski  _et al._ [2019] M. J. Dolinski, A. W. Poon, and W. Rodejohann, Annu. Rev. Nucl. Part. Sci. 69, 219 (2019).
  * Schechter and Valle [1982] J. Schechter and J. W. F. Valle, Phys. Rev. D 25, 2951 (1982).
  * Gando  _et al._ [2016] A. Gando _et al._ (KamLAND-Zen Collaboration), Phys. Rev. Lett. 117, 082503 (2016).
  * Agostini  _et al._ [2019] M. Agostini _et al._ , Science 365, 1445 (2019).
  * Deppisch  _et al._ [2012] F. F. Deppisch, M. Hirsch, and H. P$\mathrm{\ddot{a}}$s, J. Phys. G 39, 124007 (2012).
  * Rodejohann [2012] W. Rodejohann, J. Phys. G 39, 124008 (2012).
  * Suhonen and Kostensalo [2019] J. Suhonen and J. Kostensalo, Front. in Phys. 7, 29 (2019).
  * Gysbers  _et al._ [2019] P. Gysbers _et al._ , Nat. Phys. 15, 428 (2019).
  * Šimkovic  _et al._ [2018a] F. Šimkovic, R. Dvornický, and D. Štefánik, EPJ Web Conf. 194, 02002 (2018a).
  * Goeppert-Mayer [1935] M. Goeppert-Mayer, Phys. Rev. 48, 512 (1935).
  * Saakyan [2013] R. Saakyan, Annu. Rev. Nucl. Part. Sci. 63, 503 (2013).
  * Barabash [2020] A. S. Barabash, Universe 6, 159 (2020).
  * Tretyak and Zdesenko [2002] V. I. Tretyak and Y. G. Zdesenko, At. Data Nucl. Data Tables 80, 83 (2002).
  * Fiorini and Niinikoski [1984] E. Fiorini and T. Niinikoski, Nucl. Instrum. Meth. A 224, 83 (1984).
  * Giuliani [2012] A. Giuliani, J. Low Temp. Phys. 167, 991 (2012).
  * Poda and Giuliani [2017] D. Poda and A. Giuliani, Int. J. Mod. Phys. A 32, 1743012 (2017).
  * Alduino  _et al._ [2018] C. Alduino _et al._ (CUORE), Phys. Rev. Lett. 120, 132501 (2018).
  * Adams  _et al._ [2020] D. Q. Adams _et al._ (CUORE), Phys. Rev. Lett. 124, 122501 (2020).
  * Pirro  _et al._ [2006] S. Pirro _et al._ , Phys. Atom. Nucl. 69, 2109 (2006).
  * Tabarelli de Fatis [2010] T. Tabarelli de Fatis, Eur. Phys. J. C 65, 359 (2010).
  * Giuliani  _et al._ [2018] A. Giuliani _et al._ , Eur. Phys. J. C 78, 272 (2018).
  * Azzolini  _et al._ [2019a] O. Azzolini _et al._ (CUPID-0), Eur. Phys. J. C 79, 583 (2019a).
  * Azzolini  _et al._ [2019b] O. Azzolini _et al._ (CUPID-0), Phys. Rev. Lett. 123, 032501 (2019b).
  * Alenkov  _et al._ [2019] V. Alenkov _et al._ (AMoRE), Eur. Phys. J. C 79, 791 (2019).
  * Armengaud  _et al._ [2017] E. Armengaud _et al._ (LUMINEU), Eur. Phys. J. C 77, 785 (2017).
  * Poda  _et al._ [2017] D. V. Poda _et al._ (LUMINEU), AIP Conf. Proc. 1894, 020017 (2017).
  * Armengaud  _et al._ [2020a] E. Armengaud _et al._ (CUPID-Mo), Eur. Phys. J. C 80, 44 (2020a).
  * Rahaman  _et al._ [2008] S. Rahaman _et al._ , Phys. Lett. B 662, 111 (2008).
  * Grigorieva  _et al._ [2017] V. D. Grigorieva _et al._ , J. Mat. Sci. Eng. B 7, 63 (2017).
  * Arnold  _et al._ [2015] R. Arnold _et al._ (NEMO-3), Phys. Rev. D 92, 072011 (2015).
  * Haller [1994] E. E. Haller, Infrared Phys. Techn. 35, 127 (1994).
  * [42] E. Armengaud _et al._ (EDELWEISS), J. Instrum. 12, P08010 (2019).
  * Hehn  _et al._ [2016] L. Hehn _et al._ (EDELWEISS), Eur. Phys. J. C 76, 548 (2016).
  * Alduino  _et al._ [2016] C. Alduino _et al._ (CUORE), Phys. Rev. C 93, 045503 (2016).
  * Azzolini  _et al._ [2018] O. Azzolini _et al._ (CUPID-0), Eur. Phys. J. C 78, 734 (2018).
  * Schmidt  _et al._ [2020] B. Schmidt _et al._ (CUPID-Mo), J. Phys.: Conf. Ser. 1468, 012129 (2020).
  * Di Domizio  _et al._ [2011] S. Di Domizio, F. Orio, and M. Vignati, J. Instrum. 6, P02007 (2011).
  * [48] R. Huang _et al._ (CUPID-Mo), arXiv:2010.04033 .
  * Schmidt  _et al._ [2013] B. Schmidt _et al._ (EDELWEISS), Astropart. Phys. 44, 28 (2013).
  * Arnold  _et al._ [2019] R. Arnold _et al._ (NEMO-3), Eur. Phys. J. C 79, 440 (2019).
  * Armengaud  _et al._ [2020b] E. Armengaud _et al._ (CUPID-Mo), Eur. Phys. J. C 80, 674 (2020b).
  * Alduino  _et al._ [2017a] C. Alduino _et al._ (CUORE-0), Eur. Phys. J. C 77, 13 (2017a).
  * Alduino  _et al._ [2017b] C. Alduino _et al._ (CUORE), Eur. Phys. J. C 77, 543 (2017b).
  * Caldwell  _et al._ [2009] A. Caldwell, D. Kollár, and K. Kröninger, Comput. Phys. Commun. 180, 2197 (2009).
  * Faddegon  _et al._ [2008] B. A. Faddegon _et al._ , Medical Physics 35, 4308 (2008).
  * Pandola  _et al._ [2015] L. Pandola, C. Andenna, and B. Caccia, Nucl. Instrum. Meth. A 350, 41 (2015).
  * Kotila and Iachello [2012] J. Kotila and F. Iachello, Phys. Rev. C 85, 034316 (2012).
  * Mirea  _et al._ [2015] M. Mirea, T. Pahomi, and S. Stoica, Rom. Rep. Phys. 67, 872 (2015).
  * Rath  _et al._ [2013] P. K. Rath _et al._ , Phys. Rev. C 88, 064322 (2013).
  * Šimkovic  _et al._ [2013] F. Šimkovic, V. Rodin, A. Faessler, and P. Vogel, Phys. Rev. C 87, 045501 (2013).
  * Vaquero  _et al._ [2013] N. L. Vaquero, T. R. Rodríguez, and J. L. Egido, Phys. Rev. Lett. 111, 142501 (2013).
  * Barea  _et al._ [2015] J. Barea, J. Kotila, and F. Iachello, Phys. Rev. C 91, 034304 (2015).
  * Hyvärinen and Suhonen [2015] J. Hyvärinen and J. Suhonen, Phys. Rev. C 91, 024613 (2015).
  * Song  _et al._ [2017] L. S. Song, J. M. Yao, P. Ring, and J. Meng, Phys. Rev. C 95, 024305 (2017).
  * Šimkovic  _et al._ [2018b] F. Šimkovic, A. Smetana, and P. Vogel, Phys. Rev. C 98, 064325 (2018b).
  * Rath  _et al._ [2019] P. K. Rath, R. Chandra, K. Chaturvedi, and P. K. Raina, Front. in Phys. 7, 64 (2019).
  * [67] W. R. Armstrong _et al._ (CUPID Interest Group), CUPID pre-CDR, arXiv:1907.09376 .


These authors contributed equally to this work]

These authors contributed equally to this work]

Corresponding author<EMAIL_ADDRESS>
# Mechanical intelligence via fully reconfigurable elastic neuromorphic
metasurfaces

M. Moghaddaszadeh  [  Department of Mechanical and Aerospace Engineering,
University at Buffalo (SUNY),
Buffalo, NY 14260-4400, USA  Department of Civil, Structural and Environmental
Engineering, University at Buffalo (SUNY), Buffalo, NY 14260-4300, USA  M.
Mousa  [  Department of Mechanical and Aerospace Engineering, University at
Buffalo (SUNY),
Buffalo, NY 14260-4400, USA  A. Aref  Department of Civil, Structural and
Environmental Engineering, University at Buffalo (SUNY), Buffalo, NY
14260-4300, USA  M. Nouh  [  Department of Mechanical and Aerospace
Engineering, University at Buffalo (SUNY),
Buffalo, NY 14260-4400, USA  Department of Civil, Structural and Environmental
Engineering, University at Buffalo (SUNY), Buffalo, NY 14260-4300, USA

###### Abstract

The ability of mechanical systems to perform basic computations has gained
traction over recent years, providing an unconventional alternative to digital
computing in off grid, low power, and severe environments which render the
majority of electronic components inoperable. However, much of the work in
mechanical computing has focused on logic operations via quasi-static
prescribed displacements in origami, bistable, and soft deformable matter. In
here, we present a first attempt to describe the fundamental framework of an
elastic neuromorphic metasurface that performs distinct classification tasks,
providing a new set of challenges given the complex nature of elastic waves
with respect to scattering and manipulation. Multiple layers of reconfigurable
waveguides are phase-trained via constant weights and trainable activation
functions in a manner that enables the resultant wave scattering at the
readout location to focus on the correct class within the detection plane. We
further demonstrate the neuromorphic system’s reconfigurability in performing
two distinct tasks, eliminating the need for costly remanufacturing.

##  I INTRODUCTION

Mechanical computing Yasuda  _et al._ (2021), a research field older than
electronic computing, has gained a lot of interest over the past few years.
Despite the limited capabilities of mechanical computers, the ability to
integrate accurate computational and morphological capabilities with minimal
energy requirements in a self-contained mechanical structure remains
invaluable Song  _et al._ (2019). The recent surge in mechanical computing
research has been motivated by new and transformative technologies requiring
rapid (yet simple) data processing, and directly encoded autonomy and
intelligence Riley  _et al._ (2022). The latter presents an opportunity to
exploit the dynamic behavior of structures, and capitalize on certain
synergies associated with physical phenomena such as material response,
deformation, and scattering to accomplish effective computation with minimal
resources DAR (2019). In keeping with these needs, several mechanical concepts
have recently emerged with the goal of performing combinational logic El Helou
_et al._ (2022) and basic mathematics Zuo  _et al._ (2018, 2022), using for
example conductive polymers El Helou  _et al._ (2021) or bistable spring-mass
chains Ion  _et al._ (2017); Bilal, Foehr, and Daraio (2017) that propagate
mechanical signals when triggered. In tandem, there have been reports of
origami-inspired metamaterials in which shear and expansion responses along
different directions change the configuration state Treml  _et al._ (2018);
Meng  _et al._ (2021); Liu  _et al._ (2023a), and snap-through mechanisms
where elastic instabilities are exploited to enable the creation of logic
gates Mei  _et al._ (2021).

In pursuit of systems capable of learning and adaptation, structures employing
physical neural architectures and elements of artificial intelligence in the
mechanical domain have arguably gained more traction Lee, Mulder, and Hopkins
(2022). Ranging from binary-stiffness beams which learn desired shape-morphing
behaviors Hopkins, Lee, and Sainaghi (2023), to Hopf oscillators Shougat  _et
al._ (2021); Shougat, Li, and Perkins (2022) , phononic Zhang, Deshmukh, and
Wang (2023) and multi-stable metamaterials Liu  _et al._ (2023b), there has
been a spurt of activity demonstrating mechanical learning via programmable
structures at different length scales.  Additionally, recent investigations
have demonstrated the utility of origami structures in physical reservoir
computing, leveraging the multistability and reconfigurability of the Muira-
ori folding technique Liu  _et al._ (2023a); Wang and Li (2023), and
exploiting the nonlinear folding dynamics of dynamic truss-frame models Bhovad
and Li (2021), to perform intelligent tasks such as input recognition,
classification, and even robotic-related functions such as emulation and
pattern generation. At their very core, neuromorphic computers (initially
proposed by Mead Mead and Ismail (1989); Mead (1990)) seek to mimic the
intricate workings of the brain’s nervous system and its vastly distributed
nature Sarkar  _et al._ (2022); Gkoupidenis (2022); van De Burgt  _et al._
(2018); Wright  _et al._ (2022); Prucnal and Shastri (2017); Shastri  _et al._
(2021); Merolla  _et al._ (2014). There are two main approaches to develop
neuromorphic systems. The first involves transferring and translating existing
digital neural architectures to physical substrates, while the second pertains
to creating new algorithms which more accurately emulate the functionality of
biological neurons and synapses Marković  _et al._ (2020).  At the hardware
level, the implementation of neuromorphic computing has traditionally been
achieved using various techniques ranging from diffractive optics Lin  _et
al._ (2018); Zhou  _et al._ (2021); Léonard  _et al._ (2022) to integrated
photonics Shen  _et al._ (2017); Ballarini  _et al._ (2020), memristors Li
_et al._ (2018); Xia and Yang (2019), and MEMS-based neural networks Barazani
_et al._ (2020); Nikfarjam  _et al._ (2023). Most recently, neuromorphic
metasurfaces have emerged as a novel approach to achieve neural functions in
optical Fu  _et al._ (2023); Léonard  _et al._ (2021); Wu  _et al._ (2020) and
acoustic Weng  _et al._ (2020) forms. In these systems, the neuromorphic
inference is rooted in the trained wave behavior of a set of metasurfaces
comprised of tunable unit cells (physical neurons). These metasurfaces, spaced
along the direction of wave propagation, form the basis of a neural network,
and the subsequent wave scattering of a physical agent (e.g., light or sound)
in each of the physical layers constitutes the interactions within these
layers. A readout mechanism is then used to interpret the results culminating
in a desired output.

While recent advances have been made in the development of photonic
neuromorphic metasurfaces, much less has been accomplished in the realm of
mechanics. In here, we present a first attempt to describe the fundamentals of
an elastic neuromorphic metasurface, which instills a level of cognitive
intelligence in a multi-layered classifier, enabling it to mechanically
execute neural network classification tasks while overcoming a new set of
challenges associated with the complex nature of elastic waves with respect to
scattering and manipulation. Despite the understandable trade-off in
computational pace and precision associated with elastoacoustic waves
(compared to their optical counterparts), the goal here is not to compete with
the performance metrics achieved via photonic systems, but rather provide a
minimal threshold of intelligence and computational function in core
mechanical parts and metallic structures which already constitute a core
component of several equipment and devices. This form of mechanical computing
is even more prudent in applications where the computational input is readily
available in mechanical form (e.g., vibrations, dynamic loads, noise, or
impinging acoustic waves) thus eliminating the need for encoding the input
features, or extreme environments which render electronic parts dysfunctional
(e.g., high magnetic fields or elevated temperatures).

In addition to the aforementioned challenges, we seek to resolve another
limitation of wave-based neuromorphic systems, namely, their inability to be
reconfigured to conduct a new task, differing from what it was initially
designed and trained for, unless the system is fully dismantled and/or
reconstructed. In this work, we demonstrate the inner-workings of multiple
layers of reconfigurable waveguides, which are phase-trained via constant
weights and trainable activation functions in a manner that enables the
resultant wave scattering at the readout location to focus on the correct
class label within the detection plane. A comprehensive framework detailing
the theory behind the unit cell, waveguide array, and overall metasurface
designs is provided, and is then used to reconstruct the classic MNIST problem
using custom tailored hyperparameters of the proposed system. Finally,
shifting from the MNIST to the Iris dataset, we demonstrate the neuromorphic
system’s ability to achieve full reconfigurability and perform distinct
classification tasks without the need for remanufacturing.

Figure 1: Waveguide design and transmission profile. (A) Reconfigurable unit
cell realized as a solid continuum. Blue and red dimensions denote fixed and
adjustable parameters. (B) Equivalent lumped parameter representation of the
unit cell. $m_{e}$ and $k_{e}$ are the effective mass and stiffness,
respectively, and $a$ is the lattice constant. (C) Effective mass and
stiffness ratios, $R_{m}$ and $R_{k}$, as functions of $t_{a}$ and $t_{m}$ at
$24$ kHz. (D) Transmission amplitude $|\bar{T}|$ and phase $\angle\bar{T}$ for
a waveguide comprised of $r=6$ unit cells. The dashed line in $|\bar{T}|$
indicates a full transmission scenario ($|\bar{T}|=1$), while the mask in
$\angle\bar{T}$ represents $|\bar{T}|>0.8$ regions. Schematic diagrams are not
drawn to scale and all dimensions are in mm.

##  II NEUROMORPHIC METASURFACE

###  II.1 Concept

We present an elastic neuromorphic metasurface composed of 3 sections: input
gates, metasurface neurons, and detection units. These are equivalent to the
input, hidden, and output layers of a digital neural network, respectively.
The system can be excited at the input gates via any excitational signal which
can be encoded into a vibrational waveform (e.g., noise, image, or dynamic
load), culminating in elastic waves which propagate through the medium
encompassing the metasurface layers. In current applications, the phase
profile of a metasurface can be designed using an appropriate delay law to
focus or steer a wavefront in a given direction. In this work, however, we
will show the ability of a neuromorphic metasurface to generate a desirable
output from a given input via intelligent, interlayer wave scattering.
Provided with a set of tunable degrees of freedom (in this case, phase
gradients between neighboring metasurface cells), the neuromorphic system
trains itself to attain an intricate configuration of meta-neurons which allow
the fully-assembled elastic structure to execute the required task.

While the model shown here takes inspiration from digital neural networks, its
architecture is notably different in important ways, in order to accommodate
the nuanced mechanics of elastic wave propagation. A conventional digital
neural network is comprised of trainable weights and biases between the
network layers. However, the activation function applied at each layer is
typically non-trainable. In a neuromorphic metasurface, we will show that such
weights represent wave scattering characteristics which are constrained
functions of the chosen geometry, inertial, and stiffness properties. Once
specified, these weights remain unchanged and are, therefore, non-trainable.
Instead, phase delays within the metasurface layers provide an alternative
tunable (and thus trainable) platform, shown to be analogous to a trainable
activation function in a neural network framework.

Figure 2: Waveguide design. (A) Displacement field for the selected waveguides
highlighted using a dark background in (B), obtained via finite element
modeling. The top panel depicts the details of the numerical simulations
showing locations of the continuity boundary conditions and the perfectly
matched layer (PML). Shown in the middle is a 3D physical realization of the
reconfigurable waveguide design, exploiting the change in bending stiffness
for rotating cross sections to achieve the desired tunability. (B) Waveguides
designed using $R_{k,m}$ and $\bar{T}$ maps in order to exhibit high
transmission at $20^{\circ}$ phase differences over a full $360^{\circ}$
range.

###  II.2 Unit cell

Driven by the need for a mechanically-tunable system, we utilize a
reconfigurable locally resonant unit cell which exhibits subwavelength
scattering. Since computational pace is directly proportional to the wave
speed in the medium, the focus here is on axial vibrations of the resonator
given the higher speed of longitudinal wave propagation compared to flexural
waves. The solid unit cell shown in Fig. 1A is designed to have effective
properties, $m_{e}$ and $k_{e}$, equivalent to these of the lumped parameter
representation (depicted in Fig. 1B), which can be retrieved from finite
element simulations of the solid cell via transfer matrix equivalence (see SI
Appendix Note 1). The unit cell shown in Fig. 1A is composed of an Aluminum
matrix ($E=70$ GPa, $\rho=2700$ kg/m3 and $\nu=0.33$), a Brass resonator
($E=97$ GPa, $\rho=8490$ kg/m3 and $\nu=0.31$) and two slender Aluminum beams
(resonator arms) which attach the resonator to the matrix. These arms are
placed along the $y$-axis, perpendicular to the longitudinal deformation of
the unit cell along the $x$-direction, enabling a coupling between the bending
modes of the resonator arms and the axial vibrations of the unit cell.

At the core of the neuromorphic metasurface theory is the ability to achieve
full phase tunability over a $2\pi$ range while maintaining high transmission
within the waveguide. The unit cell, therefore, needs to be reconfigurable in
order to admit different values of $m_{e}$ and $k_{e}$ as needed. The
effective widths of the resonator arm $t_{a}$ and the vertical section of the
unit cell matrix $t_{m}$, both indicated on Fig. 1A, play a central role in
the bending stiffness of both parts. Thus, changes in these two parameters
significantly alter the dynamics of the unit cell described by $m_{e}$ and
$k_{e}$. Control over these effective widths can be exercised via a 3D
realization of the unit cell with rotating arms having rectangular cross-
sections, as shown in Fig. 2A and successfully implemented in literature
Attarzadeh, Callanan, and Nouh (2020). We define the effective mass and
stiffness ratios as $R_{m}=m_{e}/m_{\text{slab}}$ and
$R_{k}=k_{e}/2k_{\text{slab}}$, respectively, where $m_{\text{slab}}=\rho Aa$
and $k_{\text{slab}}=EA/[a(1-\nu^{2})]$ denote the mass and axial stiffness,
respectively, of an aluminum slab of the same size. Figure 1C depicts the
variation of $R_{m}$ and $R_{k}$ as functions of $t_{a}$ and $t_{m}$ at a
frequency of $24$ kHz. The behavior shows that while both thicknesses affect
both effective ratios (i.e., mass and stiffness), $R_{m}$ and $R_{k}$ are
observed to be largely dependent on $t_{a}$ and $t_{m}$, respectively.

Figure 3: Reconfigurable elastic neuromorphic metasurface. (A) Schematic
illustration of the assembled neuromorphic system. Six vertical cross
sections, referred to as spines and labeled $1$ through $6$, define the
different mechanical components. Spines $1$ and $6$ mark the input and output
(readout) layers, with $K$ sources and $N$ detection units. Two metasurface
layers span the regions between spines $2\to 3$ and $4\to 5$ housing $J$ and
$L$ meta-neurons, respectively. Each meta-neuron consists of $6$ identical
unit cells and retains a unique design based on tuned values of $t_{a}$ and
$t_{m}$ (see Fig. 1). Regions before, between, and after the metasurface
layers are homogenous aluminum plates. (B) Logic diagram depicting the neural
architecture of the elastic neuromorphic system, emphasizing the non-trainable
weights $\mathbf{W}^{(1\rightarrow 2)}$, $\mathbf{W}^{(3\rightarrow 4)}$, and
$\mathbf{W}^{(5\rightarrow 6)}$, as well as the trainable activation functions
$\mathbf{m}^{(2\rightarrow 3)}$ and $\mathbf{m}^{(4\rightarrow 5)}$. All
layers are dense.

###  II.3 Scattering and waveguide design

The scattering properties of a waveguide comprising $r$ unit cells can be
evaluated starting with the dispersion relation, which can be derived as
follows Lee  _et al._ (2018):

$m_{e}\omega^{2}=2k_{e}\sin^{2}(\frac{\kappa a}{2})$ (1)

where $\kappa=\omega/c$ is the Bloch wavenumber and $\omega$ is the angular
frequency of the harmonic wave. The transmission coefficient $\bar{T}$ can be
computed from:

$\bar{T}=\frac{2}{2\cos{(r\kappa
a)}+\mathrm{i}(zz_{s}^{-1}+z_{s}z^{-1})\sin{(r\kappa a)}}$ (2)

where $\mathrm{i}$ is the imaginary unit, $z=\rho Ac$ is the mechanical
impedance, and $c=\sqrt{E/[\rho(1-\nu^{2})]}$ is the axial wave speed (See SI
Appendix Note 2 for detailed derivation). The system’s characteristic
impedance $z_{s}$ can be defined as the ratio of $\sqrt{\frac{m_{e}k_{e}}{2}}$
to $\cos{(\frac{\kappa a}{2})}$. For a waveguide consisting of $r=6$ unit
cells, shown in the bottom panel of Fig. 1B, we start by considering a wide
range of $m_{e}$ and $k_{e}$ values, and obtain $\kappa$ from Eq. (1). This is
then used to determine the corresponding $z_{s}$, allowing the transmission
amplitude $|\bar{T}|$ and phase $\angle\bar{T}$ to be evaluated by using Eq.
(2), as shown in Fig. 1D. A high transmission ($|\bar{T}|>0.8$) mask is
applied to the phase map to facilitate the selection of $R_{k}$ and $R_{m}$
combinations that maintain a minimum acceptable transmission at a given phase.

As indicated earlier, the neuromorphic metasurface needs to maintain
sufficient tunability for good training and subsequent performance. Making use
of $R_{k,m}$ (Fig. 1C) and the transmission amplitude and phase maps (Fig.
1D), eighteen arrays are generated with intervals of approximately
$20^{\circ}$ phase difference as summarized in the table provided in Fig. 2B.
These can be the building blocks of the neuromorphic metasurface. Finite
element modeling is used to validate the obtained transmission and phase
delays for six select cases from the table, and the resultant displacement
field is demonstrated in Fig. 2A. Finally, in the culminating neuromorphic
metasurface, each unit cell array in the metasurface is treated as a singular
point, neglecting lateral propagating waves. As such, the Aluminum slabs are
chosen to be $1$ mm thicker than the array itself to mitigate vertical
interactions that are unaccounted for and ensure isolation from adjacent
arrays.

##  III OPERATIONAL THEORY

###  III.1 Neural architecture

The neuromorphic metasurface’s neural architecture consists of four
consecutive layers: an input layer, two metasurface layers, and an output
layer, as illustrated in Fig. 3. These are equivalent to the input, hidden,
and output layers of a digital neural network, respectively. For a given
classification problem, training and test samples with $K$ features are
reshaped into a single dimensional vector with the same number of elements.
This vector is then fed to the elastic medium in the form of $K$ displacement
point sources at the input layer, with the input features manifesting
themselves in the amplitude of each excitation point. At the far side of the
neuromorphic metasurface (output layer), a set of $N$ detection units are
defined corresponding to the $N$ labels of the classification problem, as a
computational readout. Each waveguide is comprised of an array of $r=6$ unit
cells, henceforth referred to as a “meta-neuron”. For a successful mechanical
neural network, the goal of the training is to distinctly engineer the two
metasurface layers such that for a given input (i.e., test sample) generating
a distinct excitation pattern stemming from the $K$ point sources constituting
this input, the scattered elastic wavefronts would predominantly focus on the
appropriate detection unit spatially corresponding to the correct $N$ label.
While the model shown here takes inspiration from digital neural networks as
can be seen in Fig. 3B, its architecture is notably different in important
ways, in order to accommodate the nuanced mechanics of elastic wave
propagation. A conventional digital neural network is comprised of trainable
weights and biases between the network layers. However, the activation
function applied at each layer is typically non-trainable. In a neuromorphic
metasurface, we show that such weights represent wave scattering
characteristics which are constrained functions of the chosen geometry,
inertial, and stiffness properties. Once specified, these weights remain
unchanged and are, therefore, non-trainable. Instead, phase delays within the
metasurface layers provide an alternative tunable (and thus trainable)
platform, shown to be analogous to a trainable activation function in a neural
network framework. The analytical framework, detailed in SI Appendix Note 3,
tracks the trajectory of incident waves as they propagate through the
different components, evaluating wave amplitudes at the six spines labeled $1$
through $6$ in Fig. 3A.

A conventional neural network is mathematically described by:

$\mathbf{y}^{(i+1)}=f(\mathbf{b}^{(i+1)}+\mathbf{W}^{(i\rightarrow
i+1)}\mathbf{y}^{(i)})$ (3)

where $\mathbf{y}^{(i)}$ is a vector of values assigned to neurons in the
$i^{\text{th}}$ layer, $\mathbf{W}^{(i\rightarrow i+1)}$ is the weights
matrix, and $\mathbf{b}^{(i+1)}$ is the bias. The architecture of the elastic
neuromorphic system is similar to the above framework in that
$\mathbf{y}^{(i)}$ is analogous to the displacement vector for the
$i^{\text{th}}$ layer $\mathbf{u}^{(i)}$ and the weights matrix
$\mathbf{W}^{(i\rightarrow i+1)}$ is equivalent to the wave propagation matrix
from one spine to the next. Furthermore, the activation function $f$ is
analogous to the metasurface effect inflicted by $\mathbf{m}^{(i\rightarrow
i+1)}$ vector on the incident wave. However, the bias $\mathbf{b}^{(i+1)}$ is
not applicable in this study. More importantly, in the elastic neuromorphic
system, the wave propagation matrices do not change once dimensions and
material properties are finalized, rendering them hyperparameters, contrary to
the weights of a digital neural network which are the primary trainable
parameters throughout the learning process. Another stark contrast is the
activation function applied at each layer of a digital neural network, which
is typically identical for all neurons. In the elastic neuromorphic system,
however, phase shifts applied by each meta-neuron of the metasurface layers
(and consequently the $\mathbf{m}^{(i\rightarrow i+1)}$ vectors) are trainable
parameters in the learning algorithm. The neuromorphic metasurface’s wave
manipulation approach can be perceived as multiple beam-forming segments
combined together in one line and steering separate parts of the incident
wave.

###  III.2 Encoding of input data features

To illustrate the performance of the elastic neuromorphic metasurface, and
without loss of generality, we utilize the classical MNIST dataset which is
widely used for training and testing in the field of machine learning. While
the system presented here is ideally suited for inputs which take a mechanical
form, MNIST has served as a commonly-adopted benchmark in the context of
photonic Wu  _et al._ (2020) and even acoustic Weng  _et al._ (2020)
neuromorphic metasurfaces in a manner that goes beyond optical digit
recognition, owing to its simplicity and high diversity. As will be detailed
later, the MNIST data used here is eventually fed into the neuromorphic
metasurface merely as a set of vibrational excitations of varying amplitude
and, as such, are taken at their face values (raw numbers) which can be
interpreted differently for different applications.

Each sample in the MNIST dataset consists of a $28\times 28$ matrix of gray-
scale pixel intensities. The matrix is flattened into a 1D vector and the
pixel intensities are used as displacement amplitudes for $K=784$ sinusoidal
excitation point sources at the input layer which share the same operational
frequency $\omega$, as depicted in the top panel of Fig. 4. At the far end of
the neuromorphic metasurface (output layer), $N=10$ detection units are
defined corresponding to the possible input classifications: $0$ to $9$. A
properly trained and designed metasurface is capable of focusing most of the
energy in the detection layer at the correct output unit, indicating the
recognized digit as the classifier’s readout. The bottom panel of the figure
depicts a different interpretation of the same problem where the elastic
deformations at the $K$ locations correspond to pressure amplitudes reaching
each input node as a result of an active sound-emanating source which is a
euclidean distance $r_{i}$, $1\leq i\leq K$, away. This equivalency between
the two tasks demonstrates the system’s comparable performance for both the
pre-processed MNIST and acoustic tasks, and the system’s ability to leverage
mechanical environments to carry out a meaningful computational task on an
input which already exists in the system’s native physical domain.

Figure 4: Encoding of input data. The top row shows the pre-processing of a
single digit sample ($2$) from the MNIST dataset prior to being fed into the
neuromorphic metasurface. The bottom row depicts a different interpretation of
the same problem where the elastic deformations at the input excitation nodes
correspond to pressure amplitudes reaching each input node as a result of an
active sound-emanating source which is a given distance away, following the
acoustic inverse square law. The output layers in both examples consist of 10
detection units that are perceived as either digits or speakers, corresponding
to the possible classes of the classification problem.

###  III.3 Training

The metasurface layers are set to be $d=2$ m apart with $J=L=1,000$ meta-
neurons in each metasurface layer. The distance between neighboring point
sources at the input layer is set to $10$ mm, equal to the width of a single
meta-neuron (see Section V for an optimization of geometrical and
hyperparameters). The model is developed in TensorFlow using custom Keras
layers and is trained in Python. The flattened input samples are cast into a
vector of complex numbers $\mathbf{u}^{(1)}_{784\times 1}$ since the weights
and the phase shift have to be complex in order to retain the displacement
components from one layer to another. Two custom dense layers are then defined
in the model for the two metasurface layers, with non-trainable weights
($\mathbf{W}^{(1\rightarrow 2)}_{1000\times 784}$ and
$\mathbf{W}^{(3\rightarrow 4)}_{1000\times 1000}$ from Eqs. (S23) and (S27),
respectively), no biases, and custom activation functions
($\mathbf{m}^{(2\rightarrow 3)}_{1000\times 1}$ and $\mathbf{m}^{(4\rightarrow
5)}_{1000\times 1}$ from Eqs. (S24) and (S28), respectively). As detailed in
SI Appendix Note 4, the codependency between the transmission amplitude and
the phase delay in each meta-neuron is considered when applying the activation
functions in the training process. Lastly, an additional custom dense layer is
defined as the output layer with non-trainable weights
($\mathbf{W}^{(5\rightarrow 6)}_{1000\times 1000}$) and with no biases. In
this final layer, the well-known Softmax activation function is implemented.

While the inference step of the neuromorphic metasurface is purely mechanical,
the training takes place in a computer Marković  _et al._ (2020) equipped with
an Intel Xeon® Gold 6230 CPU @ $2.10$ GHz and a RAM of 128 GB, running on
Microsoft Windows 10. The learning is carried out using $60,000$ training
samples and a $10\%$ split for validation. The Sparse Categorical Cross-
entropy loss function, which is highly effective in multi-class classification
problems, is adopted, and the Adam optimizer is employed to improve the
training performance via a defined learning rate schedule. After training, the
model is evaluated using $10,000$ (blind) test samples to report its accuracy.
Upon ensuring that the training has been executed with high efficacy, the
trained metasurface vectors $\mathbf{m}^{(i\rightarrow i+1)}$ are used to
determine the corresponding design parameters $t_{a}$ and $t_{k}$ for each and
every meta-neuron.

Figure 5: Performance of the elastic neuromorphic metasurface. (A) Output of
the neuromorphic metasurface for two digit recognition test samples from the
MNIST dataset corresponding to the digits “$2$” (left) and “$5$” (right). The
top row shows the intensities (blue bars) and probabilities (red dashed lines)
at different detection units within the output (readout) layer. Bottom and
middle rows show the scattered wavefields for trained and randomly initialized
neuromorphic metasurfaces, respectively. (B) Classification accuracy (right)
and model loss (left) for the training dataset, and confusion matrix (bottom)
for the testing dataset. The training dataset is divided into 30 batches
processed at each epoch. (C) Examples of erroneous recognition: “$2$” and
“$9$” misinterpreted as “$6$” and “$4$”, respectively, due to subtle shape
similarities, corresponding to the E1 and E2 markings on the confusion matrix
in (B).

##  IV PERFORMANCE

Figure 5A showcases the output of the neuromorphic metasurface for two digit
recognition samples corresponding to the digits “$2$” and “$5$”. The top row
illustrates the distribution of intensity and probability over the output
layer at different detection units. The latter is computed as the ratio of the
summation of point intensities at the detection units to the total intensity
of the output layer. As evident from these plots, the neuromorphic metasurface
is able to successfully classify both samples. The shown wavefields at the
bottom of Fig. 5A further corroborate this by depicting the wave scattering
through the different layers of the neuromorphic metasurface from the input
location at spine $1$ to the output (readout) location at spine $6$. The
displacements in each of the three sections shown ($1\to 2$, $3\to 4$, and
$5\to 6$) are individually normalized with respect to the same section for
better visualization. The middle row of Fig. 5A displays the wavefield
resulting from a randomized, untrained design, further confirming the
effectiveness of the trained wave scattering.

The top panel of Fig. 5B summarizes the loss and accuracy as functions of the
epoch number. An accuracy of $87.5$% in $100$ epochs was achieved taking
approximately $2.6$ hours of training on the reported PC in Section III.3.
Wavefields showing successful detection of all numeric digits are provided in
SI Appendix Note 5. The overall performance is evaluated using a confusion
matrix for the trained model based on its performance on the entire testing
dataset, as demonstrated in the bottom panel of Fig. 5B. While this provides
important insights into the model’s capability of correctly classifying the
tested samples, it also helps identify potential areas of improvement, such as
recognizing specific digit shapes that may be challenging to accurately
predict. As a case in point, the labels E1 and E2 marked on the confusion
matrix correspond to erroneous recognition of the digits “$2$” and “$9$” as
$6$ and $4$, respectively, due to subtle shape similarities illustrated in
Fig. 5C.

##  V HYPERPARAMETERS

The size and dimensions of the elastic neuromorphic metasurface are scalable,
and can be downsized depending on the operational frequency and waveguide
design. Moreover, the computational accuracy of the elastic neuromorphic
metasurface is affected by several hyperparameters that, if adequately tuned,
can considerably improve several performance metrics. These are aspects which
are set at the design stage and control the learning process, as opposed to
being derived via training. The contribution of each of these to the overall
classification accuracy is evaluated by individually altering parameters of
interest while holding the rest unchanged from the values reported earlier, as
shown in Table 1. It is important to acknowledge that, for efficiency
purposes, the model used in these studies is trained for a predetermined
number of epochs, reaching a near-stabilized accuracy level, while the
training duration may be shorter compared to the main model.

Table 1: Effect of different hyperparameters stated in the left column of the table on the overall classification accuracy. Cases typed in boldface refer to the main model. No. of layers | 2 | 3 | 4 | 5
---|---|---|---|---
87.5% | 90.7 | 91.9% | 93.0%
Distance between layers | 0.3 m | 1 m | 2 m | 3 m
87.0% | 87.3% | 87.5% | 86.9%
No. of neurons | 100 | 250 | 784 | 1,000
45.1% | 70.7% | 86.7% | 87.5%
Width of detection units | 0.3 m | 0.5 m | 0.7 m | 1 m
89.2% | 88.7% | 88.3% | 87.5%
Training metric | Intensity | Displacement
87.5% | 82.5%
Digit processing | Gray-scale | Binary
87.5% | 85.8%

As can be inferred from the first row of Table 1, increasing the number of
metasurface layers improves the accuracy in a proportional manner (See SI
Appendix Note 6 for details). A similar trend is observed for both the
distance between layers $d$ and the number of neurons in each layer ($J$ or
$L$). Nonetheless, a trade-off between accuracy and practicality is
unavoidable. Accuracy gains become insignificant beyond a certain range while
the structural size multiplies, detrimentally influencing wave intensity at
the detection units, which, in turn, adversely affects the accuracy as noted
in the last column. Contrarily, reducing the dimensions of the detection units
marginally enhances the accuracy which, while appearing advantageous,
concurrently permits a substantial portion of the energy to be redirected
away, resulting in low intensities at the detection units. As evident from the
fifth row of Table 1, employing the intensity of the wavefield at the output
layer as the evaluation metric yields improved accuracy compared to the
displacement. Finally, although the utilization of binary simplifications for
the sample digit pixels at the input layer could offer more flexibility in
fabrication and testing, allowing on/off actuation, reducing the number of
sources through combination, and eliminating the need for amplitude adjustment
at the excitation location, it is observed that a gray-scale representation
leads to a slightly improved accuracy, as can be seen in the last row of Table
1.

Beyond the discussed hyperparameters, another important factor is the system’s
scalability, which directly influences practical implementation. This includes
the total number of meta-neurons and the overall footprint (length and width)
affected by the chosen number of layers, the number of meta-neurons per layer,
and the inter-layer spacing. To explore this aspect, we re-engineered the
system for significant downsizing using two design variations. Both designs
utilize five layers, with either 250 or 100 meta-neurons per layer, and
achieve reductions of $92.5\%$ and $98\%$, respectively of the overall area
(see SI Appendix Note 7). While significantly smaller in size and meta-neuron
count compared to the original model, the downsized designs achieved accuracy
levels of $79\%$ and $75\%$, respectively. While this confirms the inherent
trade-off between accuracy and practicality, it is clear that the system
maintains an acceptable performance even when substantially reduced in scale.

It is worth noting that in addition to the optimization of the aforementioned
parameters for enhanced accuracy and an overall better performance, there
exists the possibility of training the homogenous plate domains before and
after the metasurface layers. This could be accomplished through the
utilization of graded stiffness or thickness profiles or optimized topological
configurations. While these domains would not be re-trainable like the
reconfigurable metasurface layers, they could conceivably be fine-tuned to sit
around a handful of input-output datasets, such that the active domains can
efficiently map each of the individual tasks.

##  VI RECONFIGURABILITY

Figure 6: Reconfigurability. (A) Change in the phase profile of the 1,000
meta-neurons in the first (left) and second (right) metasurface layers
following training as the neuromorphic system switches from the MNIST (M) to
the Iris (I) classification task. $\phi_{j}^{(2\rightarrow 3)}$ and
$\phi_{l}^{(4\rightarrow 5)}$ represent phase delays imposed by the first and
second metasurface layers, respectively (see SI Appendix Note 3). The insets
reveal phase reconfiguration of the metasurface unit cell via resonator arm
rotation, effectively switching the trainable meta-neurons from (M) to (I).
(B) Confusion matrix generated from the testing data depicting the system’s
accuracy in the Iris problem. (C) Intensities (blue bars) and probabilities
(red dashed lines) at the three detection units within the output layer for
input samples corresponding to, from left to right, Setosa, Versicolor, and
Virginica.

It is critically important to note that metasurface-based neuromorphic systems
developed to date, even those deployed in non-mechanical domains have not yet
overcome the need to reconstruct their core wave-scattering components in
order to retrain and carry out a different task. Whether it is TiO2 pillars Wu
_et al._ (2020) or SiO2-filled silicon slots Fu  _et al._ (2023) which form
the waveguide array, these geometries are fixed once manufactured. Therefore,
owing to their inability to alter trainable nodes post design, the
computational capacity of such systems remains limited to executing the single
computational task that the constitutive cells were trained and physically
constructed for. Motivated by the need to overcome this drawback, we
demonstrate here the ability of the elastic neuromorphic metasurface to
exercise an unprecedented degree of reconfigurability, allowing it to cater to
distinct classification tasks, forgoing reconstruction and eliminating the
need for costly remanufacturing. To illustrate, we re-utilize the baseline
model generated for MNIST digit recognition and retrain it without changing
any of its hyperparameters to execute the Iris flower classification task
Fisher (1936). The new task involves the binning of Iris flower samples into
three distinct classes (Setosa, Versicolor, and Virginica) based on four
features: sepal length, sepal width, petal length, and petal width. The values
of these features represent the amplitudes of four excitation point sources
placed at the input layer (i.e., $K=4$), each two separated by $1.5$ m. The
number of meta-neurons in the two metasurface layers is kept at $J=L=1,000$,
while $N=3$ detection units corresponding to the three classes are defined at
the output (readout) layer. The model is trained on $70\%$ of the available
dataset ($105$ samples) and tested on the remaining $30\%$, yielding an
accuracy of $97\%$ and $96\%$ for the training and testing datasets,
respectively.

Figure 6A shows the precise phase profiles for all $1,000$ meta-neurons in the
first (Spine $2\to 3$) and second (Spine $4\to 5$) metasurface layers for both
the originally trained (MNIST, or M) and retrained (Iris, or I) datasets. The
small insets reveal the phase reconfiguration of the metasurface unit cell via
resonator arm rotation, effectively switching the trainable meta-neurons from
(M) to (I). Figure 6B shows the confusion matrix for the 3 classes of the Iris
dataset generated from the testing data. Furthermore, we show the intensity
(blue bars) and probability (red dashed lines) distribution over the three
detection units within the output (readout) layer in Fig. 6C, confirming the
neuromorphic system’s robust ability to relearn a diverse set of tasks within
the same mechanical platform.

##  VII CONCLUSIONS

In this work, the foundations of metasurface-based neuromorphic computing have
been introduced to the field of structural mechanics, utilizing guided in-
plane vibrational waves in an elastic substrate. Edge excitations depicting
the features of an input sample from two distinct multivariate datasets were
fed to the neuromorphic system in the form of an array of spatially
equidistant mechanical monopoles. These excitations were made to propagate
through multiple layers of trained waveguides, eventually focusing the bulk
portion of scattered wave intensity on the correct label identifier on a
prescribed detection plane. Exploiting analogies between the different
components of the proposed neuromorphic assembly and neural architectures, a
customized neural network was defined accounting for non-trainable, constant
wave propagation matrices (weights) and trainable phase delays (activation
functions).  While the proposed neuromorphic system exhibits a high level of
mechanical intelligence, it’s equally important to consider its advantages
against existing physical neural architectures. Diffractive neural networks,
for example, share the principle of wave propagation through layered
structures using passive optical components mimicking artificial neurons to
perform intelligent tasks. However, their bulky designs and limited
reconfigurability for different tasks pose challenges. MEMS-based neural
networks, on the other hand, offer appealing features like adaptability to
complex problems and low power consumption by mimicking biological neural
networks. However, their reliance on microscopic moving parts introduces
significant fabrication complexity and hinders miniaturization efforts. Owing
to the inherent reconfigurable nature of our physical meta-neurons within the
metasurface arrays, the proposed system represents an embodiment of
metasurface-based neuromorphic computers which are no longer constrained to a
single unaltered task, and one where neurons can mechanically recompose their
local phase profile to adapt to a new training scheme and execute a new task,
on demand while maintaining a high performance accuracy.

## Supplementary Material

See the supplementary material for a comprehensive understanding of the
theoretical background of this work, and additional results referred to in
this manuscript.

###### Acknowledgements.

This work was supported by the US Army Research Office (ARO) and the
Mechanical Behavior of Materials program under Grant W911NF-23-1-0078.

## Data Availability

The data that supports the findings of this study are available within the
article and its supplementary material.

## References

## References

  * Yasuda  _et al._ (2021) H. Yasuda, P. R. Buskohl, A. Gillman, T. D. Murphey, S. Stepney, R. A. Vaia, and J. R. Raney, “Mechanical computing,” Nature 598, 39–48 (2021).
  * Song  _et al._ (2019) Y. Song, R. M. Panas, S. Chizari, L. A. Shaw, J. A. Jackson, J. B. Hopkins, and A. J. Pascall, “Additively manufacturable micro-mechanical logic gates,” Nature communications 10, 882 (2019).
  * Riley  _et al._ (2022) K. S. Riley, S. Koner, J. C. Osorio, Y. Yu, H. Morgan, J. P. Udani, S. A. Sarles, and A. F. Arrieta, “Neuromorphic metamaterials for mechanosensing and perceptual associative learning,” Advanced Intelligent Systems 4, 2200158 (2022).
  * DAR (2019) “Darpa nature as computer (nac) program,” https://www.darpa.mil/program/nature-as-computer (2019).
  * El Helou  _et al._ (2022) C. El Helou, B. Grossmann, C. E. Tabor, P. R. Buskohl, and R. L. Harne, “Mechanical integrated circuit materials,” Nature 608, 699–703 (2022).
  * Zuo  _et al._ (2018) S.-Y. Zuo, Y. Tian, Q. Wei, Y. Cheng, and X.-J. Liu, “Acoustic analog computing based on a reflective metasurface with decoupled modulation of phase and amplitude,” Journal of Applied Physics 123, 091704 (2018).
  * Zuo  _et al._ (2022) S. Zuo, C. Cai, X. Li, Y. Tian, and E. Liang, “Performing broadband and tunable mathematical operations based on acoustic reconfigurable metasurfaces,” Journal of Physics D: Applied Physics 55, 354001 (2022).
  * El Helou  _et al._ (2021) C. El Helou, P. R. Buskohl, C. E. Tabor, and R. L. Harne, “Digital logic gates in soft, conductive mechanical metamaterials,” Nature Communications 12 (2021).
  * Ion  _et al._ (2017) A. Ion, L. Wall, R. Kovacs, and P. Baudisch, “Digital mechanical metamaterials,” in  _Proceedings of the 2017 CHI Conference on Human Factors in Computing Systems_ , CHI ’17 (Association for Computing Machinery, New York, NY, USA, 2017) p. 977–988.
  * Bilal, Foehr, and Daraio (2017) O. R. Bilal, A. Foehr, and C. Daraio, “Bistable metamaterial for switching and cascading elastic vibrations,” Proceedings of the National Academy of Sciences 114, 4603–4606 (2017).
  * Treml  _et al._ (2018) B. Treml, A. Gillman, P. Buskohl, and R. Vaia, “Origami mechanologic,” Proceedings of the National Academy of Sciences 115, 6916–6921 (2018).
  * Meng  _et al._ (2021) Z. Meng, W. Chen, T. Mei, Y. Lai, Y. Li, and C. Chen, “Bistability-based foldable origami mechanical logic gates,” Extreme Mechanics Letters 43, 101180 (2021).
  * Liu  _et al._ (2023a) Z. Liu, H. Fang, J. Xu, and K.-W. Wang, “Discriminative transition sequences of origami metamaterials for mechanologic,” Advanced Intelligent Systems 5, 2200146 (2023a).
  * Mei  _et al._ (2021) T. Mei, Z. Meng, K. Zhao, and C. Q. Chen, “A mechanical metamaterial with reprogrammable logical functions,” Nature Communications 12 (2021).
  * Lee, Mulder, and Hopkins (2022) R. H. Lee, E. A. Mulder, and J. B. Hopkins, “Mechanical neural networks: Architected materials that learn behaviors,” Science Robotics 7 (2022).
  * Hopkins, Lee, and Sainaghi (2023) J. B. Hopkins, R. H. Lee, and P. B. Sainaghi, “Using binary-stiffness beams within mechanical neural-network metamaterials to learn,” Smart Materials and Structures (2023).
  * Shougat  _et al._ (2021) M. R. E. U. Shougat, X. Li, T. Mollik, and E. Perkins, “A hopf physical reservoir computer,” Scientific Reports 11, 19465 (2021).
  * Shougat, Li, and Perkins (2022) M. R. E. U. Shougat, X. Li, and E. Perkins, “Dynamic effects on reservoir computing with a hopf oscillator,” Physical Review E 105, 044212 (2022).
  * Zhang, Deshmukh, and Wang (2023) Y. Zhang, A. Deshmukh, and K.-W. Wang, “Embodying multifunctional mechano-intelligence in and through phononic metastructures harnessing physical reservoir computing,” Advanced Science 10, 2305074 (2023).
  * Liu  _et al._ (2023b) Z. Liu, H. Fang, J. Xu, and K.-W. Wang, “Cellular automata inspired multistable origami metamaterials for mechanical learning,” Advanced Science 10, 2305146 (2023b).
  * Wang and Li (2023) J. Wang and S. Li, “Building intelligence in the mechanical domain—harvesting the reservoir computing power in origami to achieve information perception tasks,” Advanced Intelligent Systems 5, 2300086 (2023).
  * Bhovad and Li (2021) P. Bhovad and S. Li, “Physical reservoir computing with origami and its application to robotic crawling,” Scientific Reports 11, 13002 (2021).
  * Mead and Ismail (1989) C. Mead and M. Ismail,  _Analog VLSI implementation of neural systems_ , Vol. 80 (Springer Science & Business Media, 1989).
  * Mead (1990) C. Mead, “Neuromorphic electronic systems,” Proceedings of the IEEE 78, 1629–1636 (1990).
  * Sarkar  _et al._ (2022) T. Sarkar, K. Lieberth, A. Pavlou, T. Frank, V. Mailaender, I. McCulloch, P. W. Blom, F. Torricelli, and P. Gkoupidenis, “An organic artificial spiking neuron for in situ neuromorphic sensing and biointerfacing,” Nature Electronics 5, 774–783 (2022).
  * Gkoupidenis (2022) P. Gkoupidenis, “Artificial neurons emulate biological counterparts to enable synergetic operation,” Nature Electronics 5, 721–722 (2022).
  * van De Burgt  _et al._ (2018) Y. van De Burgt, A. Melianas, S. T. Keene, G. Malliaras, and A. Salleo, “Organic electronics for neuromorphic computing,” Nature Electronics 1, 386–397 (2018).
  * Wright  _et al._ (2022) L. G. Wright, T. Onodera, M. M. Stein, T. Wang, D. T. Schachter, Z. Hu, and P. L. McMahon, “Deep physical neural networks trained with backpropagation,” Nature 601, 549–555 (2022).
  * Prucnal and Shastri (2017) P. R. Prucnal and B. J. Shastri,  _Neuromorphic photonics_ (CRC press, 2017).
  * Shastri  _et al._ (2021) B. J. Shastri, A. N. Tait, T. Ferreira de Lima, W. H. Pernice, H. Bhaskaran, C. D. Wright, and P. R. Prucnal, “Photonics for artificial intelligence and neuromorphic computing,” Nature Photonics 15, 102–114 (2021).
  * Merolla  _et al._ (2014) P. A. Merolla, J. V. Arthur, R. Alvarez-Icaza, A. S. Cassidy, J. Sawada, F. Akopyan, B. L. Jackson, N. Imam, C. Guo, Y. Nakamura, _et al._ , “A million spiking-neuron integrated circuit with a scalable communication network and interface,” Science 345, 668–673 (2014).
  * Marković  _et al._ (2020) D. Marković, A. Mizrahi, D. Querlioz, and J. Grollier, “Physics for neuromorphic computing,” Nature Reviews Physics 2, 499–510 (2020).
  * Lin  _et al._ (2018) X. Lin, Y. Rivenson, N. T. Yardimci, M. Veli, Y. Luo, M. Jarrahi, and A. Ozcan, “All-optical machine learning using diffractive deep neural networks,” Science 361, 1004–1008 (2018).
  * Zhou  _et al._ (2021) T. Zhou, X. Lin, J. Wu, Y. Chen, H. Xie, Y. Li, J. Fan, H. Wu, L. Fang, and Q. Dai, “Large-scale neuromorphic optoelectronic computing with a reconfigurable diffractive processing unit,” Nature Photonics 15, 367–373 (2021).
  * Léonard  _et al._ (2022) F. Léonard, E. J. Fuller, C. M. Teeter, and C. M. Vineyard, “High accuracy single-layer free-space diffractive neuromorphic classifiers for spatially incoherent light,” Optics Express 30, 12510–12520 (2022).
  * Shen  _et al._ (2017) Y. Shen, N. C. Harris, S. Skirlo, M. Prabhu, T. Baehr-Jones, M. Hochberg, X. Sun, S. Zhao, H. Larochelle, D. Englund, _et al._ , “Deep learning with coherent nanophotonic circuits,” Nature photonics 11, 441–446 (2017).
  * Ballarini  _et al._ (2020) D. Ballarini, A. Gianfrate, R. Panico, A. Opala, S. Ghosh, L. Dominici, V. Ardizzone, M. De Giorgi, G. Lerario, G. Gigli, _et al._ , “Polaritonic neuromorphic computing outperforms linear classifiers,” Nano Letters 20, 3506–3512 (2020).
  * Li  _et al._ (2018) Y. Li, Z. Wang, R. Midya, Q. Xia, and J. J. Yang, “Review of memristor devices in neuromorphic computing: materials sciences and device challenges,” Journal of Physics D: Applied Physics 51, 503002 (2018).
  * Xia and Yang (2019) Q. Xia and J. J. Yang, “Memristive crossbar arrays for brain-inspired computing,” Nature materials 18, 309–323 (2019).
  * Barazani  _et al._ (2020) B. Barazani, G. Dion, J.-F. Morissette, L. Beaudoin, and J. Sylvestre, “Microfabricated neuroaccelerometer: integrating sensing and reservoir computing in mems,” Journal of Microelectromechanical Systems 29, 338–347 (2020).
  * Nikfarjam  _et al._ (2023) H. Nikfarjam, M. Megdadi, M. Okour, S. Pourkamali, and F. Alsaleem, “Energy efficient integrated mems neural network for simultaneous sensing and computing,” Communications Engineering 2, 19 (2023).
  * Fu  _et al._ (2023) T. Fu, Y. Zang, Y. Huang, Z. Du, H. Huang, C. Hu, M. Chen, S. Yang, and H. Chen, “Photonic machine learning with on-chip diffractive optics,” Nature Communications 14, 70 (2023).
  * Léonard  _et al._ (2021) F. Léonard, A. S. Backer, E. J. Fuller, C. Teeter, and C. M. Vineyard, “Co-design of free-space metasurface optical neuromorphic classifiers for high performance,” ACS Photonics 8, 2103–2111 (2021).
  * Wu  _et al._ (2020) Z. Wu, M. Zhou, E. Khoram, B. Liu, and Z. Yu, “Neuromorphic metasurface,” Photonics Research 8, 46–50 (2020).
  * Weng  _et al._ (2020) J. Weng, Y. Ding, C. Hu, X.-F. Zhu, B. Liang, J. Yang, and J. Cheng, “Meta-neural-network for real-time and passive deep-learning-based object recognition,” Nature communications 11, 6309 (2020).
  * Attarzadeh, Callanan, and Nouh (2020) M. Attarzadeh, J. Callanan, and M. Nouh, “Experimental observation of nonreciprocal waves in a resonant metamaterial beam,” Physical Review Applied 13, 021001 (2020).
  * Lee  _et al._ (2018) H. Lee, J. K. Lee, H. M. Seung, and Y. Y. Kim, “Mass-stiffness substructuring of an elastic metasurface for full transmission beam steering,” Journal of the Mechanics and Physics of Solids 112, 577–593 (2018).
  * Fisher (1936) R. A. Fisher, “The use of multiple measurements in taxonomic problems,” Annals of eugenics 7, 179–188 (1936).