# Autoregressive Attention Neural Networks for Non-Line-of-Sight User Tracking
with
Dynamic Metasurface Antennas

Kyriakos Stylianopoulos1, Murat Bayraktar2, Nuria González-Prelcic2, and
George C. Alexandropoulos1
1Department of Informatics and Telecommunications, National and Kapodistrian
University of Athens
Panepistimiopolis Ilissia, 15784 Athens, Greece
2Department of Electrical and Computer Engineering, North Carolina State
University, 27695 Raleigh, USA
e-mails: {kstylianop<EMAIL_ADDRESS>{mbayrak<EMAIL_ADDRESS>

###### Abstract

User localization and tracking in the upcoming generation of wireless networks
have the potential to be revolutionized by technologies such as the Dynamic
Metasurface Antennas (DMAs). Commonly proposed algorithmic approaches rely on
assumptions about relatively dominant Line-of-Sight (LoS) paths, or require
pilot transmission sequences whose length is comparable to the number of DMA
elements, thus, leading to limited effectiveness and considerable measurement
overheads in blocked LoS and dynamic multipath environments. In this paper, we
present a two-stage machine-learning-based approach for user tracking,
specifically designed for non-LoS multipath settings. A newly proposed
attention-based Neural Network (NN) is first trained to map noisy channel
responses to potential user positions, regardless of user mobility patterns.
This architecture constitutes a modification of the prominent vision
transformer, specifically modified for extracting information from high-
dimensional frequency response signals. As a second stage, the NN’s
predictions for the past user positions are passed through a learnable
autoregressive model to exploit the time-correlated channel information and
obtain the final position predictions. The channel estimation procedure
leverages a DMA receive architecture with partially-connected radio frequency
chains, which results to reduced numbers of pilots. The numerical evaluation
over an outdoor ray-tracing scenario illustrates that despite LoS blockage,
this methodology is capable of achieving high position accuracy across various
multipath settings.

###### Index Terms:

Localization, tracking, dynamic metasurface antennas, deep learning,
autoregressive attention networks.

##  I Introduction

††footnotetext: This work has been supported by the EU H2020 RISE-6G project
under grant number 10101701.

Metasurface-based antennas, such as Reconfigurable Intelligent Surfaces (RISs)
[1] [2] and Dynamic Metasurface Antenna arrays (DMAs) [3, 4] are key enablers
of smart radio environments [5], and can offer benefits across various network
objectives, such as localization, sensing, and Radio Frequency (RF) mapping.
In many recent cases, RISs have been endowed with receive RF Chains (RFCs) and
play the role of large-aperture receivers, or sensors, of uplink signals,
either by deploying multiple single-RF metasurfaces [6] or partially-connected
multi-RFC architectures with reflecting [7] and hybrid reflecting/sensing [8]
meta-elements.

When solely reflecting surfaces are considered in bistatic sensing
architectures [9, 10, 11, 12], a two-stage procedure is most commonly adopted
in estimating channel or environment components (e.g., Angles of Arrivals
(AoAs)), and then, performing the position estimation. Typically, the
information extraction over the signals is carried out via subspace based
algorithms (e.g., MUSIC [13] and its extensions [14, 15]) or Compressed
Sensing (CS) approaches (e.g. [11, 9]), while the position estimation can be
obtained via minimizing the Cramér-Rao Bound (CRB) [16]. Recent methodologies
stemming from Machine Learning (ML) are either designed to augment
conventional approaches, such as the Deep-Augmented MUSIC approach [15], or
may focus on active sensing setups where the algorithms try to find favorable
metasurface configurations that lead to sufficient estimations. To that end,
very recently, a Long-Short-Term Memory (LSTM) neural network was proposed in
[17] for beamforming toward the static user position, and a similar strategy
based on reinforcement learning was used in [18] to localize passive objects.

While subspace/CS-based and ML approaches have demonstrated great potential,
they often rely on extensive signal sampling to perform estimations (i.e., in
the order of the number of RIS elements [19, 20, 6, 12, 21]) that may well
exceed the coherence time of the channel. This problem is only exacerbated by
the fact that, switching configurations in metasurfaces introduce non-
negligible delays (which may be even comparable to the channel coherence frame
[22]), and that most of the available techniques have non-linear computational
complexity [1]. Besides, the vast majority of localization approaches, even
the ones that follow similar ML-based methodologies [23, 24], assume Line-of-
Sight (LoS) propagation environments between the receiver/metasurface and the
user [23], or may entail large sensing and preprocessing overheads [24, 25].

Contributions: Taking the above considerations into account, in this paper, we
design a monostatic sensing system for user tracking in multipath non-LoS
environments, where a DMA with a limited number of RFCs observes an equal
number of pilot signals to perform a noisy estimation of the channel. A novel
modified Multi-Head Self-Attention (mMHSA) Neural Network (NN) architecture is
proposed that learns a mapping between the channel vectors and the (static)
user positions within the considered environment. To perform tracking, we
furthermore train an AutoRegressive (AR) model that receives as input the
time-agnostic predictions of the mMHSA network along sampled user
trajectories, and outputs a smoothed version of the predictions that better
fit the mobility patterns of the users.

Notation: Vectors (matrices) are given in bold lowercase (bold uppercase)
typeface and calligraphic letters denote sets. $(\cdot)^{T}$ and $(\cdot)^{H}$
return matrix transposition and hermitian transposition, respectively, while
$\|\cdot\|$ is the Euclidean norm and $|\cdot|$ provides a set’s cardinality.
$\mathcal{CN}(\cdot,\cdot)$ denotes the complex normal distribution and
$\jmath\triangleq\sqrt{-1}$. For a vector $\boldsymbol{v}$,
$\boldsymbol{v}_{i:j}$ denotes the sliced-vector constructed by taking the
elements of $\boldsymbol{v}$ from the $i$-th position up to, and including,
the $j$-th position.

##  II System Model and Objective

###  II-A System Model

We consider the uplink of a communication system, in which a single-antenna
user moves along pseudorandom 2D trajectories and periodically transmits
localization requests to a fixed-position planar DMA-based access point
comprising $N$ meta-elements, each with a tunable electromagnetic response
(effectively, a phase shift). It is assumed that the DMA meta-elements are
attached in groups in $N_{\rm RF}$ RFCs, so that all $N_{E}\triangleq N/N_{\rm
RF}$ elements of each DMA row are connected to a common RFC. By modeling the
signal propagation inside the DMA microstrip (assuming negligible attenuation)
and adopting a Lorentzian-constrained phase model for its unit elements, the
frequency response $w_{i,j}$ of the element located in the $i$-th row
($i=1,\dots,N_{\rm RF})$ and $j$-th column ($j=1,\dots,N_{E})$ of the
metasurface’s grid can be compactly expressed as follows [26]:

$\displaystyle w_{i,j}\triangleq
0.5\left(\jmath+\exp{(\jmath\phi_{i,j})}\right)\exp{(\jmath\rho_{j}\beta)},$
(1)

where $\phi_{i,j}\in[-\pi/2,\pi/2]$ denotes the controllable phase shift
induced by the DMA for the respected element, $\rho_{j}$ is the horizontal
position of the $j$-th element in its microstrip, and $\beta$ is the
wavenumber associated with the internal propagation, which is expressed as
$\beta\triangleq(2\pi/\lambda)\sqrt{\epsilon}$ with $\lambda$ being the
wavelength of the carrier signal and $\epsilon$ denotes the permittivity of
the microstrip (usually, $\epsilon=6$ [26]). For notation purposes, let us
denote the concatenated vector of phase shifts as:

$\boldsymbol{\phi}\triangleq[\phi_{1,1},\phi_{1,2},\dots,\phi_{2,1},\dots,\phi_{N_{\rm
RF},N_{E}}]^{T},$ (2)

and the corresponding vector of the DMA configuration as

$\boldsymbol{w}\triangleq[w_{1,1},w_{1,2},\dots,w_{2,1},\dots,w_{N_{\rm
RF},N_{E}}]^{T}.$ (3)

The phase shift configuration of the overall DMA satisfies the mean power
constraint $\mathbb{E}\left[\|\boldsymbol{w}\|^{2}\right]\leq P_{\rm max}$.
From now on in this paper, we will be referring to the elements of
$\boldsymbol{\phi}$ and $\boldsymbol{w}$ as $\phi_{n}$ and $w_{n}$
($n=1,\dots,N$), respectively, where $n=(i-1)N_{E}+j$.

At each localization request time slot $t$, the user is located at a 2D
position $\boldsymbol{p}^{t}$ and transmits a fixed pilot symbol $x$ with unit
power under a multipath channel. Orthogonal Frequency-Division Multiplexing
(OFDM) with $L$ subcarriers is assumed and the channel at each $t$-th time
slot and $l$-th subcarrier with $l=1,2,\ldots,L$ is modeled by the vector
$\boldsymbol{h}^{t,l}\in\mathbb{C}^{N\times 1}$. By using the notation
$\boldsymbol{\tilde{w}}_{i}^{t}\triangleq\boldsymbol{w}^{H}_{(i-1)N_{E}+1:iN_{E}}\in\mathbb{C}^{1\times
N_{E}}$ for the sub-vector of the configuration associated with the $i$-th DMA
RFC at the $t$-th time slot, and representing the associated channel sub-
vector as
$\boldsymbol{\tilde{h}}_{i}^{t,l}\triangleq\boldsymbol{h}_{(i-1)N_{E}+1:iN_{E}}^{t,l}\in\mathbb{C}^{N_{E}\times
1}$, the received signal at the $i$-th RFC of the DMA at the access point at
the $l$-th subcarrier and $t$-th time slot can be expressed as:

$\displaystyle
y_{i}^{t,l}\triangleq\boldsymbol{\tilde{w}}_{i}^{t}\boldsymbol{\tilde{h}}_{i}^{t,l}x+\boldsymbol{\tilde{w}}_{i}^{t}\boldsymbol{\tilde{n}}^{t,l},$
(4)

where
$\boldsymbol{\tilde{n}}^{t,l}\sim\mathcal{CN}(\mathbf{0}_{N_{E}},\sigma^{2}\mathbf{I}_{N_{E}})$
denotes the thermal noise vector at the receiver with $\mathbf{0}_{N_{E}}$ and
$\mathbf{I}_{N_{E}}$ being the $N_{E}$-element zero vector and $N_{E}\times
N_{E}$ identity matrix, respectively. Let us finally denote the concatenated
receive vector from all $N_{\rm RF}$ DMA’s RFCs as
$\boldsymbol{y}^{t,l}\triangleq[y^{t,l}_{1},y^{t,l}_{2},\ldots,y^{t,l}_{N_{\rm
RF}}]^{T}$.

To accurately model the non-LoS and multipath nature of the considered
localization setting, we adopt the urban outdoor “O1 - Blocked” scenario of
the DeepMIMO simulator [27] for the carrier frequency $f_{c}=28$ GHz and $500$
KHz bandwidth. This simulator provides ray-tracing-derived path details (see
[28] for a relevant survey with similar simulators) for user positions along a
grid with $20$ cm separation. A top view of the simulation setup is given in
Fig 1.

Figure 1: Top view of the non-LoS tracking scenario “O1” with one blocking and
two main reflecting wall surfaces. The receiving DMA is positioned at
$(235,485)$. The figure is colored by the received signal strength values, as
provided via ray-tracing from [27] using $3$ non-LOS channel paths. Values
shown in normalized units. The red curve illustrates a randomly generated user
trajectory of $100$ location points of $20$ cm separation.

###  II-B Tracking Objective

Let us denote as $\mathcal{T}$ the trajectory instance within a predefined
space containing the $|\mathcal{T}|$ sequential user positions. In this paper,
we make use of the Bézier curves [29] with $5$ control points for each
trajectory by further discretizing each random user position to the closest
one appearing on the DeepMIMO data set. The applied localization schemes
output an estimated user position $\boldsymbol{\hat{p}}^{t}$ at every time
slot $t$, therefore, the instantaneous localization error is given by $e_{\rm
pos}^{t}\triangleq\sqrt{\|\boldsymbol{p}^{t}-\boldsymbol{\hat{p}}^{t}\|^{2}}$
and our proposed tracking objective entails minimizing the mean of this error
within a trajectory. In mathematical terms, we consider the optimization
problem:

$\displaystyle\mathcal{OP}:$
$\displaystyle\underset{\boldsymbol{\\{\hat{p}}^{t}\\}_{t=1}^{|\mathcal{T}|}}{\min}~{}\frac{1}{|\mathcal{T}|}\sum_{t=1}^{|\mathcal{T}|}e_{\rm
pos}^{t}.$ (5)

##  III Proposed User Tracking Methodology

In this section, a two-stage approach for online user tracking with a DMA-
based receiver is presented. We first describe the channel estimation process,
followed by the overall user tracking procedure. Then, the details of the
proposed mMHSA NN that maps channel estimates with position estimations, and
the considered AR model that leverages the latter estimations to better fit
user mobility patterns, are provided.

###  III-A Channel Estimation

By setting $\phi_{i,j}=-\pi/2\triangleq\phi^{0}$ in (1), the response of each
($i,j$)-th DMA element gets zeroed. Using this setting for all indices $i$ and
all $j\neq j^{\prime}$, while simultaneously letting
$\phi_{i,j^{\prime}}\neq\phi^{0}$, the received signal in (4) at each $i$-th
DMA RFC at the $l$-th subcarrier and $t$-th time slot can be expressed as:

$y_{i}^{t,l}=\tilde{w}_{i,j^{\prime}}^{t}\tilde{h}_{i,j^{\prime}}^{t,l}x+\tilde{w}_{i,j^{\prime}}^{t}\tilde{n}_{i,j^{\prime}}^{t,l},$
(6)

which indicates that only the $j^{\prime}$-the element of the $i$-th
microstrip contributes to the received signal. By digitally removing the known
effects of $\tilde{w}_{i,j^{\prime}}^{t}$ and $x$ from $y_{i}^{t,l}$, an
estimate for $\tilde{h}_{i,j^{\prime}}^{t,l}$ can be obtained with an
associated estimation error given by
$\tilde{w}_{i,j^{\prime}}^{t}x\tilde{n}_{i,j^{\prime}}^{t,l}$. By sequentially
setting $j^{\prime}=1,\dots,N_{E}$, an estimate of the complete channel vector
is attained using $N_{E}$ repeated pilot signals. We note that for practical
purposes, $N_{E}$ is expected to be $o(N)$, thus, the overhead of this pilot-
assisted channel estimation procedure is limited.

###  III-B Tracking Procedure

Assume a user moving along its predefined trajectory at some arbitrary time
step $t$. The first stage of the proposed user tracking approach involves
training a regression NN, termed modified Multi-Head Self-Attention (mMHSA),
that receives a channel measurement matrix as described above and outputs an
initial position estimate
$\boldsymbol{\tilde{p}}^{t}=f_{\boldsymbol{\theta}}^{\rm
mMHSA}(\boldsymbol{H}^{t})$, which is parameterized by an arbitrary vector of
weights, $\boldsymbol{\theta}$. The mMHSA network disregards any temporal
correlation appearing in tracking problems. To that end, the second stage of
the procedure involves the refinement of those initial estimates by accounting
for the user’s learned mobility patterns. Concretely, assume that, at time
$t$, mMHSA has provided estimates for all past and current positions
$\\{\boldsymbol{\tilde{p}}^{k}\\}_{k=1}^{t}$. An AR model is then applied on
those initial predictions to obtain the final output position estimate
$\boldsymbol{\hat{p}}^{t}$, as
$\boldsymbol{\hat{p}}^{t}=f_{\gamma,\boldsymbol{z}}^{\rm
AR}(\\{\boldsymbol{\tilde{p}}^{k}\\}_{k=1}^{t})$, so that correlations between
past estimates can be learned.

Since this approach falls under supervised learning, both models require a
training procedure with known user positions. To that end, assume that $D$
trajectories of length $\mathcal{|T|}$ are collected as
$\mathcal{D}=\\{\\{(\boldsymbol{H}^{t,d},\boldsymbol{p}^{t,d})\\}_{t=1}^{|\mathcal{T}|}\\}_{d=1}^{D}$
and a shuffled version of $\mathcal{D}$ is compiled as
$\mathcal{D^{\prime}}=\\{(\boldsymbol{H}^{q},\boldsymbol{p}^{q})\\}_{q=1}^{|\mathcal{T}|D}$.
The mMHSA is, thus, trained via standard Stochastic Gradient Descent (SGD)
[30] over $\mathcal{D^{\prime}}$ on the loss function
$J_{\mathcal{D^{\prime}}}(\boldsymbol{\theta})=1/({|\mathcal{T}|D})\sum_{q=1}^{|\mathcal{T}|D}\sqrt{\|\boldsymbol{p}^{q}-f_{\boldsymbol{\theta}}^{\rm
mMHSA}(\boldsymbol{H}^{q})\|^{2}}$ which minimizes the tracking objective (5).
Once mMHSA is trained, it is used to predict the (known) positions for all
channel instances in $\mathcal{D}$, to obtain
$\\{\\{\boldsymbol{\tilde{p}}^{t,d}\\}_{t=1}^{|\mathcal{T}|}\\}_{d=1}^{D}$.
Therefore, the AR model is trained via SGD to minimize the error in the
estimation between the ground-truth positions and the initial estimates of
mMHSA as
$J_{\mathcal{D}}(\gamma,\boldsymbol{z})=1/({|\mathcal{T}|D})\sum_{d=1}^{D}\sum_{t=1}^{|\mathcal{T}|}\sqrt{\|\boldsymbol{p}^{t,d}-f_{\gamma,\boldsymbol{z}}^{\rm
AR}(\\{\boldsymbol{\tilde{p}}^{k,d}\\}_{k=1}^{t})\|^{2}}$.

###  III-C Modified Multi-Head Self-Attention Network (mMHSA)

Figure 2: The proposed mMHSA architecture for position estimation based on
noisy channel estimates, inspired by ViT [31]. The design of the patches,
estimation token, and positional embeddings have been modified from the
original architecture. Red rectangles signify learnable blocks.

The mMHSA step aims to map noisy channel estimates to position estimations.
Accounting for the $N$ DMA meta-elements and the $L$ OFDM subcarriers, and by
separating the real and imaginary channel components along a matrix dimension,
the estimated channel at each $t$-th time slot can be expressed as a 3D matrix
$\boldsymbol{H}^{t}\in\mathbb{R}^{2\times N\times L}$, which resembles the
matrix dimensions commonly appearing in computer vision applications. It is
noted that, due to the geometric nature of the considered channel model, the
$\boldsymbol{H}^{t}$ elements are expected to exhibit high spatial correlation
[32]. Under the latter observations, our mMHSA architecture uses components of
the state-of-the-art Vision Transformer (ViT), which was originally proposed
for image classification tasks [31]. The NN architecture is included in Fig. 2
following the ViT naming and components. In particular, each input
$\boldsymbol{H}^{t}$ is first rearranged to two dimensions of size $(N,2L)$,
which is then sliced into $\mathtt{p}\times\mathtt{p}$ smaller patches. The
patches are then concatenated and passed through a linear mapping layer to
attain the initial embeddings of dimension $d_{\rm hidden}$. Two learnable
tokens are appended in the embeddings instead of a single class token, since
our architecture is intended for prediction of 2D vectors. Positional
embeddings are then added, similar to [31], although the frequency of the
embeddings is set to $2N$ instead of $10000$ in our architecture. The
embeddings are then forwarded to $\mathtt{n}_{\rm blocks}$ sequential
transformer encoder blocks. Each encoder is comprised an mMHSA layer followed
by a Multi-Layer Perceptron (MLP), as proposed in [33]. Finally, the 2D
estimation tokens are extracted along each embedding dimension and an MLP head
is responsible for extracting the position estimation
$\boldsymbol{\tilde{p}}^{t}$ out of them.

###  III-D AR-Model-Based Position Estimation

This stage receives the estimates
$\boldsymbol{\\{\tilde{p}}^{t}\\}_{t=1}^{|\mathcal{T}|}$ of the user position
along a trajectory of $\mathcal{T}$ sequential user positions (i.e.,
$|\mathcal{T}|$ mMHSA outputs) and feeds an AR model with an exponentially
decaying weighting factor to refine the position estimate at each $t$-th time
slot as:

$\boldsymbol{\hat{p}}^{t}=f_{\gamma,\boldsymbol{z}}^{\rm
AR}(\\{\boldsymbol{\tilde{p}}^{k}\\}_{k=1}^{t})=\sum_{k=1}^{t}{\rm
pow}({\rm\sigma}(\gamma),t-k+1)\boldsymbol{z}\circ\boldsymbol{\tilde{p}}^{k},$
(7)

where $\gamma$ and $\boldsymbol{z}\in\mathbb{R}^{2}$ denote the model’s
parameters, which are optimized through SGD over a data set of collected user
trajectories and corresponding mMHSA predictions as described in Section
III-B. In this expression, $\sigma(\cdot)$ denotes the sigmoid function that
ensures that the decaying factor resides in $(0,1)$, $\circ$ indicates
element-wise multiplication, and ${\rm pow}(\cdot,\cdot)$ denotes
exponentiation.

##  IV Numerical Evaluation

For the numerical evaluation, we adopt DeepMIMO’s “O1 - Blocked” scenario
[27], where the base station with index $3$ plays the role of the receiving
DMA positioned at $(235,485,6)$. The geometry is given in Fig. 1, where
approximately $10^{6}$ user positions are sampled at a spacing of $20$ cm
inside a $110\times 40$ m outdoor urban area. The user travels along $75000$
randomly generated 2D trajectories of $|\mathcal{T}|=100$ positions at a mean
separation of $2.5$ m, remaining at a fixed height of $2$ m. We keep $N=256$,
$L=4$, and $N_{\rm RF}=16$ as fixed values and we set the estimation noise
floor power per symbol to $-60$ dBm. The mMHSA network is built with $3$
Transformer Encoder blocks with $2$ heads and a hidden embedding dimension of
$16$. A learning rate of $0.002$ and a training period of $60$ epochs were
used for both models. Only half the data set was used for training the mMHSA,
and the rest was used for testing to limit the number of random trajectory
positions that are seen by the network during training, regardless of the
internal regularization components.

Since typical localization and tracking works assume LoS [9, 17, 15, 16] or
more involved systems and procedures that don’t facilitate real time user
tracking [11, 18], we have implemented a fingerprinting localization scheme
[34] as a baseline for position estimation. Disregarding trajectories and
estimation errors, the Received Signal Strength Indicators (RSSIs) for all
$N_{\rm RF}$ RFCs and $L$ subcarriers corresponding to the data points of the
training set are stored in a data base. To make the position estimation, this
baseline observes the current RSSI from the user’s pilots and outputs the
position of the closest RSSI match in the data base (i.e. nearest neighbor
heuristic).

The position estimation errors along $25000$ randomly sampled test
trajectories are given in Fig. 3 over increasing numbers of signal paths, both
for the initial mMHSA estimations, the two-step mMHSA and AR process, and the
fingerprinting baseline. It can be inferred that the proposed architecture has
substantial overheads over the baseline, especially when fewer reflection
paths are present. The performance difference is decreasing as the effects of
multipath get more pronounced, with the mMHSA-AR methodology resulting in
smaller position errors nonetheless. The addition of the AR step is beneficial
to the overall scheme by reducing the estimation error roughly $1$-$2$ m in
all cases. When a single reflection path is considered, this approach is able
to attain a mean position error of $2.47$ m. Finally, it is worth stating that
while the presence of multipath has an adverse effect in the estimation
accuracy, the induced errors nearly plateau after a certain threshold of $5$
paths.

Figure 3: Localization error of the proposed methodology and baseline over
different number of signal paths.

##  V Conclusion

In this paper, the problem of online user tracking via a receiving DMA has
been considered, dealing specifically with mobile users in non-LoS multipath
environments. The DMA architecture was used to provide noisy estimates of the
channel with low pilot overheads, which are then fed into a variation of the
ViT architecture to map frequency responses to user position estimates. The
tracking problem was handled by an AR model over initial estimates along
sampled user trajectories to further reduce the estimation error based on
mobility patterns. The numerical evaluation has illustrated that this approach
is able to attain high localization performance in a large area non-LoS
scenario with multipath channels derived from the DeepMIMO ray-tracing
simulator. This approach was able to attain a mean position error below $2.5$
m when a single channel path is present, consistently outperforming the
fingerprinting baseline in all cases, while the effects of the increasing
number of paths are shown to be bounded.

## References

  * [1] M. Jian, G. C. Alexandropoulos, E. Basar, C. Huang, R. Liu, Y. Liu, and C. Yuen, “Reconfigurable intelligent surfaces for wireless communications: Overview of hardware designs, channel models, and estimation techniques,” _Intell. Converged Netw._ , vol. 3, no. 1, pp. 1–32, Mar. 2022.
  * [2] C. Huang, A. Zappone, G. C. Alexandropoulos, M. Debbah, and C. Yuen, “Reconfigurable intelligent surfaces for energy efficiency in wireless communication,” _IEEE Trans. Wireless Commun._ , vol. 18, no. 8, pp. 4157–4170, Aug. 2019.
  * [3] N. Shlezinger, G. C. Alexandropoulos, M. F. Imani, Y. C. Eldar, and D. R. Smith, “Dynamic metasurface antennas for 6G extreme massive MIMO communications,” _IEEE Wireless Commun._ , vol. 28, no. 2, pp. 106–113, Apr. 2021.
  * [4] I. Gavras, M. A. Islam, B. Smida, and G. C. Alexandropoulos, “Full duplex holographic MIMO for near-field integrated sensing and communications,” in _Proc. European Signal Proces. Conf._ , Helsinki, Finland, Sep. 2023.
  * [5] E. Calvanese Strinati, G. C. Alexandropoulos, H. Wymeersch, B. Denis, V. Sciancalepore, R. D’Errico, A. Clemente, D.-T. Phan-Huy, E. D. Carvalho, and P. Popovski, “Reconfigurable, intelligent, and sustainable wireless environments for 6G smart connectivity,” _IEEE Commun. Mag._ , vol. 59, no. 10, pp. 99–105, Oct. 2021.
  * [6] G. C. Alexandropoulos, I. Vinieratou, and H. Wymeersch, “Localization via multiple reconfigurable intelligent surfaces equipped with single receive RF chains,” _IEEE Wireless Commun. Lett._ , vol. 11, no. 5, pp. 1072–1076, May 2022.
  * [7] J. He, A. Fakhreddine, C. Vanwynsberghe, H. Wymeersch, and G. C. Alexandropoulos, “3D localization with a single partially-connected receiving RIS: Positioning error analysis and algorithmic design,” _IEEE Trans. Veh. Technol._ , vol. 72, no. 10, pp. 13 190–13 202, Oct. 2023.
  * [8] G. C. Alexandropoulos, N. Shlezinger, I. Alamzadeh, M. F. Imani, H. Zhang, and Y. C. Eldar, “Hybrid reconfigurable intelligent metasurfaces: Enabling simultaneous tunable reflections and sensing for 6G wireless communications,” _arXiv preprint arXiv:2104.04690_ , 2021.
  * [9] J. He, H. Wymeersch, L. Kong, O. Silvén, and M. Juntti, “Large intelligent surface for positioning in millimeter wave MIMO systems,” in _Proc. IEEE VTC-Spring_ , 2020, pp. 1–5.
  * [10] J. Yuan, M. Wen, Q. Li, E. Basar, G. C. Alexandropoulos, and G. Chen, “Receive quadrature reflecting modulation for RIS-empowered wireless communications,” _IEEE Trans. Veh. Technol._ , vol. 70, no. 5, pp. 5121–5125, May 2021.
  * [11] J. He, A. Fakhreddine, H. Wymeersch, and G. C. Alexandropoulos, “Compressed-sensing-based 3D localization with distributed passive reconfigurable intelligent surfaces,” in _Proc. IEEE Int. Conf. Acoustics, Speech, and Signal Process._ , Rhodes, Greece, 2023, pp. 1–5.
  * [12] H. Chen, H. Kim, M. Ammous, G. Seco-Granados, G. C. Alexandropoulos, S. Valaee, and H. Wymeersch, “RISs and sidelink communications in smart cities: The key to seamless localization and sensing,” _IEEE Commun. Mag._ , vol. 61, no. 8, pp. 140–146, 2023.
  * [13] R. Schmidt, “Multiple emitter location and signal parameter estimation,” _IEEE Trans. Antennas Propag._ , vol. 34, no. 3, pp. 276–280, Mar. 1986.
  * [14] B. Friedlander, “The root-MUSIC algorithm for direction finding with interpolated arrays,” _Signal Process._ , vol. 30, no. 1, pp. 15–29, Jan. 1993.
  * [15] J. P. Merkofer, G. Revach, N. Shlezinger, T. Routtenberg, and R. J. G. van Sloun, “DA-MUSIC: Data-driven DoA estimation via deep augmented MUSIC algorithm,” _arXiv preprint, arXiv:2109.10581_ , 2023.
  * [16] Y. Liang, W. Liu, Q. Shen, W. Cui, and S. Wu, “A review of closed-form Cramér-Rao bounds for DOA estimation in the presence of Gaussian noise under a unified framework,” _IEEE Access_ , vol. 8, pp. 175 101–175 124, 2020.
  * [17] Z. Zhang, T. Jiang, and W. Yu, “Active sensing for localization with reconfigurable intelligent surface,” in _Proc. IEEE Int. Conf. Commun._ , Rome, Italy, Jun. 2023.
  * [18] J. Hu, H. Zhang, K. Bian, D. M. M., Z. Han, and L. Song, “MetaSensing: Intelligent metasurface assisted RF 3D sensing by deep reinforcement learning,” _IEEE J. Sel. Areas Commun._ , vol. 39, no. 7, pp. 2182–2197, Jul. 2021.
  * [19] G. C. Alexandropoulos, S. Samarakoon, M. Bennis, and M. Debbah, “Phase configuration learning in wireless networks with multiple reconfigurable intelligent surfaces,” in _Proc. IEEE Global Commun. Conf._ , Taipei, Taiwan, Dec. 2020.
  * [20] C. Huang, Z. Yang, G. C. Alexandropoulos, K. Xiong, L. Wei, C. Yuen, and Z. Zhang, “Hybrid beamforming for RIS-empowered multi-hop terahertz communications: A drl-based method,” in _Proc. IEEE Global Commun. Conf._ , Taipei, Taiwan, Dec. 2020.
  * [21] G. C. Alexandropoulos, K. Stylianopoulos, C. Huang, C. Yuen, M. Bennis, and M. Debbah, “Pervasive machine learning for smart radio environments enabled by reconfigurable intelligent surfaces,” _Proc. IEEE_ , vol. 110, no. 9, pp. 1494–1525, Sep. 2022.
  * [22] ETSI white paper, “Reconfigurable intelligent surfaces (RIS): Use cases, deployment scenarios and requirements,” online document, Apr. 2023, GR RIS 001 - V1.1.1.
  * [23] Y. Chen, N. González-Prelcic, T. Shimizu, H. Lu, and C. Mahabal, “Sparse recovery with attention: A hybrid data/model driven solution for high accuracy position and channel tracking at mmWave,” _arXiv preprint, arXiv:2308.13908_ , 2023.
  * [24] Y. Chen, N. González-Prelcic, T. Shimizu, and H. Lu, “Learning to localize with attention: From sparse mmwave channel estimates from a single BS to high accuracy 3D location,” _arXiv preprint, arXiv:2307.00167_ , 2023.
  * [25] J. Gante, G. Falcão, and L. Sousa, “Deep learning architectures for accurate millimeter wave positioning in 5G,” _Neural Process. Lett._ , vol. 51, no. 1, pp. 487–514, Feb 2020.
  * [26] N. Shlezinger, O. Dicker, Y. C. Eldar, I. Yoo, M. F. Imani, and D. R. Smith, “Dynamic metasurface antennas for uplink massive MIMO systems,” _IEEE Trans. Commun._ , vol. 67, no. 10, pp. 6829–6843, Oct. 2019.
  * [27] A. Alkhateeb, “DeepMIMO: A generic deep learning dataset for millimeter wave and massive MIMO applications,” _arXiv preprint, arXiv:1902.06435_ , 2019.
  * [28] H. Choi, J. Chung, J. Oh, G. C. Alexandropoulos, , and J. Choi, “WiThRay: A versatile 3D Ray-tracing simulator for smart wireless environments,” _IEEE Access_ , vol. 11, pp. 56 822–56 845, Jun. 2023.
  * [29] M. K. Agoston, _Computer Graphics and Geometric Modelling: Implementation & Algorithms_. Springer London, 2005.
  * [30] I. Goodfellow, Y. Bengio, and A. Courville, _Deep Learning_. MIT Press, 2016.
  * [31] A. Kolesnikov _et al._ , “An image is worth 16x16 words: Transformers for image recognition at scale,” in _Proc. Int. Conf. Learning Representations (ICLR)_ , May 2021, pp. 1–22.
  * [32] G. C. Alexandropoulos, P. T. Mathiopoulos, and N. C. Sagias, “Switch-and-examine diversity over arbitrary correlated Nakagami-$m$ fading channels,” _IEEE Trans. Veh. Technol._ , vol. 59, no. 4, pp. 2080–2087, May 2010.
  * [33] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, L. Kaiser, and I. Polosukhin, “Attention is all you need,” in _Proc. Conf. Neural Inf. Process. Syst. (NIPS)_ , Long Beach, USA, Dec. 2017, pp. 1–11.
  * [34] S. Yiu, M. Dashti, H. Claussen, and F. Perez-Cruz, “Wireless RSSI fingerprinting localization,” _Signal Process._ , vol. 131, pp. 235–244, Feb. 2017.


# Dependence of the kinetic energy absorption capacity of bistable mechanical
metamaterials on impactor mass and velocity

R. Fancher  I. Frankel  K. Chin  M. Abi Ghanem  B. MacNider  L. S. Shannahan
J. F. Berry  M. Fermen-Coker  A. J. Boydston  N. Boechler

###### Abstract

Using an alternative mechanism to dissipation or scattering, bistable
structures and mechanical metamaterials have shown promise for mitigating the
detrimental effects of impact by reversibly locking energy into strained
material. Herein, we extend prior works on impact absorption via bistable
metamaterials to computationally explore the dependence of kinetic energy
transmission on the velocity and mass of the impactor, with strain rates
exceeding $10^{2}$ s-1. We observe a large dependence on both impactor
parameters, ranging from significantly better to worse performance than a
comparative linear material. We then correlate the variability in performance
to solitary wave formation in the system and give analytical estimates of
idealized energy absorption capacity under dynamic loading. In addition, we
find a significant dependence on damping accompanied by a qualitative
difference in solitary wave propagation within the system. The complex
dynamics revealed in this study offer potential future guidance for the
application of bistable metamaterials to applications including human and
engineered system shock and impact protection devices.

[1]organization=Department of Mechanical and Aerospace Engineering, University
of California, San Diego, city=La Jolla, state=CA, country=US
[4]organization=Department of Chemistry, University of Wisconsin at Madison,
city=Madison, state=WI, country=US [5]organization=Institut Lumière Matière,
niversité Claude Bernard Lyon 1, CNRS, city=Villeurbanne, state=F-69622,
country=France [3]organization=Army Research Laboratory, city=Aberdeen Proving
Ground, state=MD, country=US [2]organization=Materials Science and
Engineering, University of California, San Diego, city=La Jolla, state=CA,
country=US

##  1 Introduction

The local magnification of mechanical forces as a result of a dynamic
collision (impact) and their detrimental effects on natural and engineered
systems have been studied extensively [1, 2, 3, 4]. A relatively new approach
for the mitigation of damage induced by impact is the use of bistable
structures, which in contrast to the more ubiquitous mechanisms of dissipation
and scattering [5, 6, 7, 8], reduces the effect of impact by reversibly
“locking” some of the energy imparted by a shock or impact into the form of
strain energy [9, 10]. Both the performance and reversibility of bistable
structures for impact mitigation are attractive, as the energy locking
mechanism could be used in conjunction with other mechanisms [11], such as
dissipation, and the structures can ostensibly be reset in a controllable
fashion for reuse.

In addition to studies of the impact absorption characteristics of single
bistable structures [9, 10], more recently, the energy absorption properties
of bistable structures configured into multi-unit-cell mechanical
metamaterials has been explored [10, 12, 13, 14, 15, 16, 17, 18, 11]. In such
studies, the material response has been studied in either relatively low rate
regimes [12, 13, 14, 15, 16, 17, 18], where in some cases the response was
suggested to be rate independent [12], while in other cases the impact
response was not related to the material’s wave dynamics [11]. This can be
placed in contrast with other studies of nonlinear solitary, or “transition”,
wave propagation in similar materials, which show significant dependence on
the system excitation [19, 20, 21, 22, 23, 24]. A notable work that brings
these two features (impact absorption via bistability and nonlinear waves)
together is that of Ref. [25], which described the effect that the shape of
the bistable potential has on energy trapping and solitary wave emission,
however, they did not study the dependence of absorption performance on impact
conditions.

Figure 1: Concept and modeling overview. (a) Photograph of a 3D printed model
of the bistable mechanical metamaterial (with fewer unit cells than modeled
herein). (b) Photograph of a 3D printed unit cell. (c) One beam in the unit
cell design, showing the displacement of the beam near $\epsilon_{2}$
simulated using FEM, where the colorscale denotes the von Mises stress (arb.).
(d) Non-dimensionalized effective stress and strain of the mechanical
metamaterial (for $r=0.14$ and $\theta=60$ degrees), based upon experimental
compression measurements (solid black line), and the 3rd order polynomial fit
of the experimentally measured curve (solid blue line), which is used in the
DEM model. The vertical dashed red and blue lines denote $\epsilon_{0}$ and
$\epsilon_{2}$, respectively, and the horizontal dotted black line denotes
zero stress. (e) Illustration of the DEM model.

In this work, we computationally study the dependence of kinetic energy
transmission in a bistable mechanical metamaterial (shown in Fig. 1) on the
velocity and mass of the impactor, in wave dominated regimes (wavelengths less
than the absorbing material’s reference length) with maximum strain rates of
approximately $195$ s-1. Kinetic energy (KE) transmission is used as a
performance metric herein as it has been previously shown to be closely
related to damaging effects, for instance, in the case of behind armor blunt
trauma [26, 27]. The computational model of our mechanical metamaterial is a
discrete element model (DEM) composed of a one-dimensional (1D) chain of
masses connected by bistable nonlinear springs, with linear intersite damping,
and a contact spring at the boundary that allows release of the impactor under
tension during rebound. We find a large dependence of KE transmission on both
impactor mass and velocity, ranging from significantly better to worse
performance than a comparative linear material. We correlate said performance
to solitary wave formation and give analytical estimates of idealized energy
absorption capacity under dynamic loading. In addition, a significant effect
is found from the inclusion of the intersite damping, accompanied by a
qualitative difference in solitary wave propagation within the system. The
complex dynamics revealed in this study offer potential future guidance for
applications including improved packaging to prevent damage during shipment
[12], personal protection equipment [7, 12, 28], and crash mitigation for
vehicles [12, 16, 29].

##  2 Description of metamaterial to be modeled

The conceptual setup described herein is the impact of a relatively rigid
“impactor” block onto an arbitrarily-sized “absorbing” material block. The
dimensions of the absorbing block were chosen as $2d$ (height) by $d$ (both
width and depth), with $d=5$ cm. Our designed absorbing material consists of a
bistable mechanical metamaterial composed of a periodic array of structured
unit cells based upon the geometry from Ref. [12]. The unit cell of our
bistable mechanical metamaterial is composed of two elastic beams surrounded
by a lower-aspect ratio monolithic “frame”. Figure 1(a,b) shows a 3D printed
representation (with fewer unit cells and layers than modeled herein) of the
bistable mechanical metamaterial, where the beams are composed of a rubber
simulant material and the frame composed of an acrylic simulant. The thickness
to length aspect ratio of the beam elements $r=T/L$ (where $T$ is the beam
thickness and $L$ is its length) and the angle of the beam $\theta$ was chosen
such that $r=0.14$ and $\theta=60$ degrees. Large deformation finite element
method (FEM) simulations of beam deformation under compression are shown in
Fig. 1(c), assuming a Neo-Hookean material model (Lamé Parameters:
$\lambda=47.191$ GPa, $\mu=25$ GPa), fixed boundaries on the bottom end (red),
with a prescribed y-direction displacement on the top end (purple, fixed in
the x-direction), and free boundaries otherwise. We model the beam separately
from the frame, because we approximate the frame itself to undergo minimal
deformation as it is much thicker (structurally stiffer) than the beams.
Furthermore, we assume that the two beams within the unit cell deform
symmetrically with respect to each other. As such the mechanical response of a
single beam can be extrapolated to determine the elasticity of the
metamaterial. The choice of the beam aspect ratio was chosen to obtain
moderately high energy locking without self-contact [12], but without too
small of an energy barrier before unsnapping. The size of the frame was chosen
to avoid self-contact while maintaining rigidity without excessive density per
unit volume of the unit cell. Both the frame and geometry are natural
candidates for future optimization studies.

##  3 Measurement of quasi-static bistable material response

Figure 1(d) shows the experimentally measured mechanical response of a single
3D printed layer of bistable metamaterial (composed of the same material and
geometry as in Fig. 1(a,b), but less depth in the direction going into the
page), along with a 3rd order polynomial fit of the measurement. The response
of the layer was measured using a mechanical test frame in displacement
control, where the rigid (acrylic simulant) top and bottom parts of the layer
were attached to grips in a mechanical testing frame, allowing measurement of
tension, as well as compression, following the snap through event. Herein a
“layer” is defined as a row of unit cells where the normal is in the direction
of the impactor velocity vector and along the long axis of the absorbing
block. In Fig. 1(d), $\sigma$ and $\varepsilon$ are the effective bistable
mechanical metamaterial stress and strain, respectively, and $E=8.42e5$ N/m2
is the measured small strain elastic modulus of the rubber simulant beam
material. The transition to negative stiffness can be seen to occur in the fit
at $\varepsilon_{0}=-0.13$ and the second stability point at
$\varepsilon_{2}=-0.46$.

##  4 Discrete element model

Figure 1(e) shows a visualization of the DEM used to simulate the dynamics of
the metamaterial undergoing impact. Each layer of unit cells in the absorbing
material were described as lumped mass layers of mass $m$ connected by
massless springs and dampers, with the impactor of mass $M$ interacting with
the absorber via a contact spring. A DEM was chosen as a reasonable model due
to the following key assumptions: i) That the model is designed to describe
uniaxial loading of the impactor on the absorber block with minimal off axis
loading effects; ii) near zero effective Poisson’s ratio of the lattice; and
iii) The lattice is composed of stiffer and larger masses lumped within the
material, separated by softer, lower mass elements. Regarding the this third
assumption, more precisely, the vibrational frequencies of the separate mass
and spring components must be much higher than the modal frequency of the two
elements combined, where the spring deforms as the mass moves like a rigid
body, such that the higher frequencies can be reasonably ignored. The sample
mass was divided into equal “mass layers” in number equal to the amount of
unit cells along the height of the material ($N$).

A “contact spring” with nonlinear stiffness, based on the Hertzian contact
model [30], between the top unit cell and the impactor mass was modeled to
describe the impactor hitting the top of the sample and allowing the impactor
to freely bounce rather than stick to the top of the lattice after impact. The
equations of motion for the impactor particle (i.e. particle $n=N+1$) is thus:

$M\ddot{y}_{N+1}={C_{1}}([y_{N}-y_{N+1}]_{+})^{3/2},$ (1)

where the $[]_{+}$ denotes spring’s inability to support tension (if the value
in brackets is negative it equals zero), $y_{n}$ is the displacement of the
$n$th layer, and $C_{1}$ is a fitting parameter (set to $-8.51e8$). As such,
the contact stiffness (shown in Fig. 1(e)) is
$c(u)=-\frac{3}{2}C_{1}([-u]_{+})^{1/2}$, where $u$ is the spring stretch
(positive in tension for both $u$ and $c(u)$). In addition to allowing the
impactor to rebound, the use of a contact spring allowed for a better
estimation of real impact conditions, including roughly describing the
impactor coming into contact with the sample at a slight relative angle, or
having asperities on the two surfaces.

We then define the force-displacement relation of an individual layer:

$F_{L}(u)=\beta_{3}(u)^{3}+\beta_{2}(u)^{2}+\beta_{1}(u),$ (2)

where $F_{L}$ is positive in tension, $\beta_{1}=2.257e4$ N/m,
$\beta_{2}=1.187e7$ N/m, and $\beta_{3}=1.524e9$ N/m are coefficients from the
fit of the mechanical response measured for the 3D printed bistable layer
(shown in Fig. 1(d)). The stiffness is thus defined as $K(u)=\partial
F_{L}(u)/\partial u$ (shown in Fig. 1(e)). Equation 2 is used in all other
elements of the DEM (corresponding to the metamaterial layers), where the top
layer of the absorber ($n=N$) is:

$\begin{split}{m}\ddot{y}_{N}&={C}_{1}(y_{N}-y_{N+1})^{3/2}-F_{L}(y_{N}-y_{N-1})\\\
&-{\eta}(\dot{y}_{N}-\dot{y}_{N-1}),\end{split}$ (3)

where $\dot{y}_{n}$ and $\eta$ are the velocity of the $n$th particle and
$\eta$ the damping coefficient, respectively.

The equations of motion for layers from $n=2$ to $n=N-1$ are given by:

$\begin{split}{m}\ddot{y}_{n}=&F_{L}(y_{n+1}-y_{n})-F_{L}(y_{n}-y_{n-1})\\\
+&{\eta}(\dot{y}_{n+1}-\dot{y}_{n})-{\eta}(\dot{y}_{n}-\dot{y}_{n-1}).\end{split}$
(4)

The equation of motion for the $n=1$ mass, next to the fixed boundary is then:

$\begin{split}{m}\ddot{y}_{1}={}&F_{L}(y_{2}-y_{1})-F_{L}(y_{1})\\\
+&{\eta}(\dot{y}_{2}-\dot{y}_{1})-{\eta}\dot{y}_{1}.\end{split}$ (5)

The equations of motion were numerically integrated using the ODE45 integrator
in MATLAB [31], given an initial velocity $V$ applied to the impactor mass, to
solve for particle displacements and velocities as a function of time.

##  5 Analytical estimate of nominal impact conditions

The kinetic energy transmitted through the half-way-point of the absorbing
block (with respect to the impactor velocity vector), or unit cell $n=N/2$,
was chosen as the performance metric (so as to avoid boundary effects), which
was then compared against that of a linearly coupled (with the exception of
the contact spring) absorbing material.

We then endeavored to estimate the impactor parameters to minimize KE
transmission through $n=N/2$ of our chosen bistable mechanical metamaterial.
The nominal impact velocity was estimated as:

$V_{0}=2c_{0}\varepsilon_{0},$ (6)

where $c_{0}$ is the long wavelength linear soundspeed of the lattice [2]. For
all simulations, a unit cell length of $a=1$ mm was used. Equation 6 is meant
to estimate the minimum velocity threshold such that the snap through to the
second stable state is induced. In an ideal scenario, the entire material
between the $n=N/2$ and the impactor would then snap to the second stable
state, and have locked the entire KE from the impactor into stored potential
energy (PE). The PE absorbed from the first half of the absorber is then:

$\text{PE}=\frac{N}{2}\int_{o}^{a\epsilon_{2}}F_{L}(u)du.$ (7)

The nominal impactor mass $M_{0}$ was solved as the only remaining unknown
when impactor KE (i.e. $MV^{2}/2$) was set equal to PE.

Figure 2: Undamped simulations. (a) KE ratio as a function of impactor
conditions. Non-dimensionalized KE, $\kappa$, for bistable (b) and linear (d)
materials at nominal impact conditions ($M_{0}$, $V_{0}$). The colorbar is
saturated at $\text{KE}_{I}/10$. Layer strain, $\epsilon$, for bistable (c)
and linear (e) materials at nominal impact conditions ($M_{0}$, $V_{0}$).

##  6 Simulated dependence on impact conditions for a conservative material

To assess the performance of the bistable mechanical metamaterial, the linear
“control” material was simulated using a linear stiffness corresponding to the
slope between zero force and the force at $\varepsilon_{0}$ for the bistable
material (equal to $1.77\beta_{1}$), identical mass, contact spring and
damping parameters to the bistable sample. Simulations with varied impactor
mass and velocity using an undamped linear material showed negligible
difference in KE transmission at the $n=N/2$ (approximately $1\%$ maximum),
between linear materials of different stiffnesses, when analyzed from impact
time until the time the first wavefront hits the bottom of the sample. We thus
define a “KE ratio”, as the maximum KE of the $N/2$ mass of the control
lattice divided by the maximum KE of the $N/2$ mass of the bistable lattice,
where a value greater than one indicates the bistable material is
outperforming the linear material in terms of minimizing KE transmission.

We initially simulated the response of $N=100$ layer conservative (i.e.
$\eta=0$) materials for duration $T_{s}=1.2\tau_{l}$ where $\tau_{l}$ is the
time for the linear wave to travel across the material once (based on
$c_{0}$). Figure 2(a) shows the KE ratio for varied $M$ and $V$, wherein the
bistable material outperforms the linear material by up to 21x, while
performing worse than the linear material particularly for impactor conditions
where $M>M_{0}$ and $V>V_{0}$. A reduction in performance above the nominal
impact conditions of $M_{0}$ and $V_{0}$ as they represent the point where the
kinetic energy of the impactor could be equally distributed across the first
half of the material and stably locked into strain energy. This is consistent
with the approximate iso-energy diagonal threshold that can be observed in
Fig. 2(a) that separates KE ratios above and below unity. This however opens
the question however, as to why the bistable material would perform less well
than the linear material without any energy locking capability. To address
this question we proceed to study the spatiotemporal response for both
systems.

To study the spatiotemporal response of the absorber materials at specific
impact conditions, we define a normalized KE
$\kappa=N(\text{KE}_{n}/\text{KE}_{I})$, where $\text{KE}_{n}$ is the KE of
the nth particle and $\text{KE}_{I}$ is the KE of the impactor, and strain
$\epsilon=u/a$. The simulation time $t$ is expressed in terms of
$\tau=t\sqrt{{\beta_{1}}/m}$. Figure 2(b,c) shows $\kappa$ and $\epsilon$ at
impactor conditions $M_{0}$, $V_{0}$, and Fig. 2(d,e) shows the linear sample
at the same impact conditions. We note that the impactor is not shown in the
plots of $\epsilon$, although it is included as the $n=N+1$ layer in the
$\kappa$ diagrams.

Figure 2(d,e) shows a minimally dispersive pulse that propagates through the
linear sample. In contrast, the bistable material (Fig. 2(b,c)) shows
nonlinear effects, demonstrating the richness of the bistable system. In Fig.
2(b,c), we see solitary, transition wave emission from the initial impact, as
shown by the straight, yellow lines that progress $\sim 25$ particles into the
sample from the point of impact. The identification of these as transition
waves is noted by strains exceeding $\varepsilon_{2}$. Looking first at each
of these emissions, we see that the first wave is the fastest, then the next
two are progressively slower, which is consistent with the amplitude dependent
wavespeed of many solitary waves [22]. Beyond $\tau\sim 4$, between the
$80-100$th particles, we see a complex oscillatory behavior, likely a
combination of low amplitude vibrations and dynamic snapping and unsnapping of
the bistable layers. Crucially, past $\sim 25$ particles from the impactor,
very little KE is transmitted further into the material.

Figure 3: Normalized KE, $\kappa$, for the simulated undamped bistable
material shown in Fig. 2(a) for the (a) worst performance (low mass and high
velocity, $M/M_{0}=10^{-1}$ and $V/V_{0}=10^{0.75}$), (b) good performance
(close to nominal mass and velocity, $M/M_{0}=10^{-0.625}$ and
$V/V_{0}=10^{0.25}$), and (c) poor performance (high mass and high velocity,
$M/M_{0}=10^{0.75}$ and $V/V_{0}=10^{0.75}$).

In order to investigate the changing performance shown in Fig. 2(a), we study
three different impactor conditions. Figure 3(a) shows a location to the top
left of Fig. 2(a), which also corresponds to the worst performance of the
sweep. At these impactor conditions ($M/M_{0}=10^{-1}$ and
$V/V_{0}=10^{0.75}$), a high amplitude, short wavelength, transition wave
travels throughout the sample, at speeds faster than the waves in the linear
material, resulting in a KE ratio of $0.0848$x. Defining $M_{R}=M/m$,
$M_{R}=1.5$ and the formation of a single solitary wave in this instance
agrees with previous studies, for instance in granular chain impacts [4],
favoring single solitary wave formation when impactor and layer mass are
closely matched. Within the context of solitary waves, nonlinear self-
localization likely contributes to the poor performance seen in such cases.
Figure 3(b) shows impactor conditions: $M/M_{0}=10^{-0.625}$ and
$V/V_{0}=10^{0.25}$ representative of good performance (KE ratio $=4.44$x),
but not as good as the nominal conditions shown in Fig. 2(b). A side by side
comparison reveals qualitative differences in the behavior; namely longer
lasting and more initial solitary waves generated in Fig. 2(b), which
transition to an oscillatory phase encompassing more particles than Fig. 3(b).
Additionally, close inspection at the $n=N/2$ location in Fig. 3(b) show wave-
fronts that travel through particle $n=N/2$ of higher KE density than any of
the crossing wave-fronts in Fig. 2(b). We observe that these wavefronts occur
later in time than the time it would take for the first solitary wave to reach
$n=N/2$ and appear to stem from dynamic unsnapping of the unit cells. Figure
3(c) shows an example for high $M$ and $V$ ($M/M_{0}=10^{0.75}$ and
$V/V_{0}=10^{0.75}$), where $M_{R}=87$, and parallels could be drawn with
known “shock” cases in granular chains [32]). A “train” of solitary waves is
emitted after impact, and appear to continue to do so along the right side of
Fig. 3(c). The KE density magnitude of subsequent emitted waves appears to
decrease, and the highest KE density solitary wave, which is also the wave
that triggers the maximum $\kappa$, is the first wave emitted from the time of
initial impact. We note the fastest transition wavespeeds are shown under
these high-energy impact conditions. The KE ratio was $0.246$x, which we note,
does not occur at the first passage of the first transition wave but rather
after the reflection of the first wave off the fixed bottom (after
$\tau\approx 5$). Reflections and interference are representative of the
complications of longer simulation duration in undamped simulations.

##  7 Effect of damping on system response

The effect of the addition of a predetermined amount of damping to the
simulations is assessed by implementing an intersite damping value of
$\eta=0.164$ Ns/m, or a normalized value of
${\eta}{V_{0}}/{\beta_{1}}{\varepsilon_{2}}a=5.8e-3$ which represents the
ratio of viscous to elastic effects. Otherwise, the simulation setup is
identical to the prior section. Figure 4 shows the results of the damped
simulations in the same format as shown previously in Fig. 2. As before, Fig.
4(b,c) shows the bistable material response at impactor conditions $M_{0}$,
$V_{0}$, and Fig. 4(d,e) shows the linear material response at the same
conditions. The maximum KE ratio in the sweep in Fig. 4(a) was $34.58$x, which
occurred at impactor conditions: $M/M_{0}=10^{-.25}$ and $V/V_{0}=10^{.25}$,
noting a shift in the optimal impact conditions with the addition of damping.
The maximum KE ratio from the sweep in Fig. 4(a) is higher than any values
observed from the undamped sweep (Fig. 2(a)). In addition, Fig. 4(a) also
shows more regions where the bistable material outperforms the linear
material, as indicated by fewer regions of white color in the figure (noting
the threshold of white to blue in the plot in Fig. 2(a) and Fig. 4(a) is a KE
ratio of one).

Figure 4: Damped simulation. (a) KE ratio as a function of impactor
conditions. Non-dimensional kinetic energy, $\kappa$, for bistable (b) and
linear (d) materials at nominal impact conditions ($M_{0}$, $V_{0}$). The
colorbar axis is saturated at $\text{KE}_{I}/10$. Layer strain, $\varepsilon$
for bistable (c) and linear materials (e) at nominal impact conditions
($M_{0}$, $V_{0}$).

Figure 4(b) shows four solitary waves that propagate in succession after
impact, then transition to the oscillatory region by approximately $n=80$. In
contrast to the undamped case, the wavespeed of each transition wave decreases
after the point of impact, as could be expected by the damping-induced
reduction of amplitude coupled with the amplitude dependent wavespeed of the
solitary waves. Figure 4(c) shows regions with layers that have snapped to
their secondary stable state and remained there, ostensibly aided by damping.
Figure 4(d,e) show similar behavior of the damped linear sample as was seen in
the undamped linear sample in Fig. 2(d,e), with the exception of a lower
maximum $\kappa$ value at $n=N/2$ when damping is included.

Further comparison of $\kappa$ for three chosen impactor conditions picked
from within Fig. 4(a) are shown in Fig. 5. Figure 5(a) shows impactor
conditions corresponding to the top left corner and the worst performance
within Fig. 4(a), representing impactor mass $M/M_{0}=10^{-1}$ and velocity
$V/V_{0}=10^{0.75}$, and a KE ratio of $0.2073$x. It is also the same impactor
conditions that resulted in the worst performance from the undamped sweep.
Figure 5(b) shows the best performance, a KE ratio of $34.58$x, at impactor
conditions $M/M_{0}=10^{-0.25}$ and $V/V_{0}=10^{0.25}$. Qualitatively, some
similarities are observed when Fig. 5(b) is compared with the undamped high
performance cases (Fig. 2(b) and Fig. 3(b)). We see an initial solitary
wavefront that transitions (in time) to an oscillatory region prior to
crossing particle $n=N/2$ (noting again the slowing wavespeed). Crucially,
there appears to be an interaction between the solitary wave fronts and
damping such that when multiple solitary wave emissions occur, it is possible
that the waves slow as they reduce in amplitude (due to damping) such that
they approach, but do not pass, the $n=N/2$ location. This may further be a
synergistic effect in the sense that the damping-induced speed reduction gives
the solitary waves even more time to decay in amplitude, and eventually
vanish, before reaching the halfway point. Figure 5(c) shows impactor
conditions $M/M_{0}=10^{0.75}$ and $V/V_{0}=10^{0.75}$, which is a region in
the upper right corner of the sweep in Fig. 4(a), showing poor performance. A
comparison of Fig. 5(c) to the undamped case (Fig. 3(c)) at the same impactor
conditions shows similar behavior with the exception of the slowing transition
waves in the damped case.

Figure 5: Normalized KE, $\kappa$, for the simulated damped bistable material
shown in Fig. 4(a) for the (a) worst performance (low mass and high velocity,
$M/M_{0}=10^{-1}$ and $V/V_{0}=10^{0.75}$), (b) best performance (close to
nominal mass and velocity, $M/M_{0}=10^{-0.25}$ and $V/V_{0}=10^{0.25}$), and
(c) poor performance (high mass and high velocity, $M/M_{0}=10^{0.75}$ and
$V/V_{0}=10^{0.75}$).

##  8 Conclusions

The kinetic energy transmission performance in response to impact of a
simulated bistable mechanical metamaterial was found to be highly dependent on
the impactor conditions (mass and velocity). In the undamped simulations, the
performance of a $N=100$ sample ranged from worse ($0.08$x) to far superior
($20.9$x) in comparison to the linear control lattice. Similarly, in the
damped sweep, the performance again ranged from worse ($0.21$x) to far
superior ($34.6$x). The presence of damping was seen to have a beneficial,
potentially synergistic, effect on the performance of the system, as indicated
by both higher maximum KE ratio values and higher minimum KE ratio values for
the same set of impactor conditions. One likely reason for this is the higher
amount of “permanent” snapping of layers prior to the $n=N/2$ point. This
permanent snapping is aided by damping both by the viscous resistance to
unsnapping and reduction of traveling waves that may cause unsnapping at later
times after subsequent reflections and constructive interference.

Two significant additional implications of the findings regarding impact
conditions and damping are the following. First, the bistable mechanism has
the potential to yield significant performance benefits in terms of KE
abatement if the material and impact conditions are well paired, but if not
well paired, the bistable material can underperform a more traditional
material, sometimes by a substantial margin. These high performance conditions
were found, via simulation, to match well with predictions made using simple
analytical estimates (namely, Eqs. 6 and 7). Second, while only a single
viscosity was used in the damped examples, the performance improvement as a
result of the addition of damping leads the authors to suspect that subsequent
“tuning” of the damping may significantly affect the performance of these
systems.

We further suggest several potential avenues for future study in the context
of both the fundamental understanding of energy transmission in bistable
systems as well as their application in impact mitigation systems. The
simulations considered herein idealized the response of a physically realized
bistable mechanical metamaterial by not considering self contact (or full
compaction) of the physical unit cells. The future inclusion of self contact
will be critical for accurately describing the response of such systems, and
can be expected to modify ranges of poor and superior impact performance. In
addition, we suggest that computational optimization has a key future role to
play, particularly given the strong nonlinearities involved herein. This may
take the form of shape optimization of the unit cell aimed at maximizing
energy absorption per unit mass density considering the metamaterial’s dynamic
response. One may also consider optimizing the unit cell to broaden the
envelope of impactor conditions (e.g. Fig. 2 and Fig. 4) wherein the bistable
material gives superior impact absorption response. Finally, we consider there
may exist additional rich dynamics and unique capabilities for impact
absorption with bistable media in higher dimensions, such as in the case of
localized point impacts, instead of the 1D (plate-impact-like) scenarios
considered herein.

##  9 Acknowledgments

This project was supported by the US Army Research Laboratory and Army
Research Office under grant no. W911NF-17-1-0595. I.F. acknowledges support
from the Department of Defense (DoD) through the National Defense Science &
Engineering Graduate (NDSEG) Fellowship Program. B.M. acknowledges support
from the U.S. Department of Energy (DOE) National Nuclear Security
Administration (NNSA) Laboratory Graduate Residency Fellowship (LRGF) under
Cooperative Agreement DE-NA0003960. The authors declare they have no competing
financial interests.

## References

  * [1] W. Schiehlen, “The long history of impact mechanics, rolling contact and multibody system dynamics,” PAMM, vol. 17, no. 1, pp. 165–166, 2017.
  * [2] M. A. Meyers, Dynamic behavior of materials.  John wiley & sons, 1994.
  * [3] C. M. Harris and A. G. Piersol, Harris’ shock and vibration handbook, vol. 5.  McGraw-Hill New York, 2002.
  * [4] V. Nesterenko, Dynamics of heterogeneous materials.  Springer Science & Business Media, 2013.
  * [5] R. Lakes and R. S. Lakes, Viscoelastic materials.  Cambridge university press, 2009.
  * [6] R. Bouferra, H. Pron, J.-F. Henry, C. Bissieux, and J. Beaudoin, “Study of the intrinsic dissipation associated to the plastic work induced by a ball impact,” International journal of thermal sciences, vol. 44, no. 2, pp. 115–119, 2005.
  * [7] Y. Li, H. Fan, and X.-L. Gao, “Ballistic helmets: Recent advances in materials, protection mechanisms, performance, and head injury mitigation,” Composites Part B: Engineering, p. 109890, 2022.
  * [8] W. Mason and H. McSkimin, “Energy losses of sound waves in metals due to scattering and diffusion,” Journal of Applied Physics, vol. 19, no. 10, pp. 940–946, 1948.
  * [9] K.-W. Wang and R. L. Harne, Harnessing bistable structural dynamics: for vibration control, energy harvesting and sensing.  John Wiley & Sons, 2017.
  * [10] Y. Cao, M. Derakhshani, Y. Fang, G. Huang, and C. Cao, “Bistable structures for advanced functional systems,” Advanced Functional Materials, vol. 31, no. 2106231, pp. 1–23, 2021.
  * [11] S. Jeon, B. Shen, N. Traugutt, Z. Zhu, L. Fang, C. Yakacki, T. Nguyen, and S. Kang, “Synergistic energy absorption mechanisms of architected liquid crystal elastomers,” Advanced Materials, vol. 34, no. 2200272, pp. 1–8, 2022.
  * [12] S. Shan, S. H. Kang, J. R. Raney, P. Wang, L. Fang, F. Candido, J. A. Lewis, and K. Bertoldi, “Multistable architected materials for trapping elastic strain energy,” Advanced Materials, vol. 27, no. 29, pp. 4296–4301, 2015\.
  * [13] T. Frenzel, C. Findeisen, M. Kadic, P. Gumbsch, and M. Wegener, “Tailored buckling microlattices as reusable light-weight shock absorbers,” Advanced Materials, vol. 28, no. 28, pp. 5865–5870, 2016.
  * [14] B. Haghpanah, A. Shirazi, L. Salari-Sharif, A. Guell Izard, and L. Valdevit, “Elastic architected materials with extreme damping capacity,” Extreme Mechanics Letters, vol. 17, pp. 56–61, 2017.
  * [15] A. Guell Izard, R. Fabian Alfonso, G. McKnight, and L. Valdevit, “Optimal design of a cellular material encompassing negative stiffness elements for unique combinations of stiffness and elastic hysteresis,” Materials and Design, vol. 135, pp. 37–50, 2017.
  * [16] L. Wu, X. Xi, B. Li, and J. Zhou, “Multi-stable mechanical structural materials,” Advanced Engineering Materials, vol. 20, no. 2, p. 1700599, 2018.
  * [17] C. S. Ha, R. S. Lakes, and M. E. Plesha, “Design, fabrication, and analysis of lattice exhibiting energy absorption via snap-through behavior,” Materials & Design, vol. 141, pp. 426–437, 2018.
  * [18] T. Raj Giri and R. Mailen, “Controlled snapping sequence and energy absorption in multistable mechanical metamaterial cylinders,” International Journal of Mechanical Sciences, vol. 204, no. 106541, pp. 1–12, 2021.
  * [19] B. Deng, P. Wang, V. Tournat, and K. Bertoldi, “Nonlinear transition waves in free-standing bistable chains,” Journal of the Mechanics and Physics of Solids, vol. 136, p. 103661, 2020.
  * [20] N. Nadkarni, A. F. Arrieta, C. Chong, D. M. Kochmann, and C. Daraio, “Unidirectional transition waves in bistable lattices,” Physical review letters, vol. 116, no. 24, p. 244501, 2016.
  * [21] J. R. Raney, N. Nadkarni, C. Daraio, D. M. Kochmann, J. A. Lewis, and K. Bertoldi, “Stable propagation of mechanical signals in soft media using stored elastic energy,” Proceedings of the National Academy of Sciences, vol. 113, no. 35, pp. 9722–9727, 2016.
  * [22] T. Dauxois and M. Peyrard, Physics of solitons.  Cambridge University Press, 2006.
  * [23] M. Frazier and D. Kochmann, “Band gap transmission in periodic bistable mechanical systems,” Journal of Sound and Vibration, vol. 388, pp. 315–326, 2017.
  * [24] S. Katz and S. Givli, “Solitary waves in a nonintegrable chain with double-well potentials,” Physical Review E, vol. 100, no. 3, p. 032209, 2019.
  * [25] S. Katz and S. Givli, “Solitary waves in a bistable lattice,” Extreme Mechanics Letters, vol. 22, pp. 106–111, 2018.
  * [26] N. R. Council et al., Opportunities in protection materials science and technology for future army applications.  National Academies Press, 2011.
  * [27] A. Sondén, D. Rocksén, L. Riddez, J. Davidsson, J. K. Persson, D. Gryth, J. Bursell, and U. P. Arborelius, “Trauma attenuating backing improves protection against behind armor blunt trauma,” Journal of Trauma and Acute Care Surgery, vol. 67, no. 6, pp. 1191–1199, 2009.
  * [28] X. Q. Wang, “The failure analysis and the innovative design of hammer head on hammer crusher based on triz theory,” in Applied Mechanics and Materials, vol. 741, pp. 85–90, Trans Tech Publ, 2015.
  * [29] J. A. Sherwood and C. C. Frost, “Constitutive modeling and simulation of energy absorbing polyurethane foam under impact loading,” Polymer Engineering & Science, vol. 32, no. 16, pp. 1138–1146, 1992.
  * [30] H. Hertz, “Uber die berührung fester elastischer körper. (on the contact of elastic solids),” J. reine angew. Math., vol. 92.
  * [31] MATLAB, version 9.7.0.1216025 (R2019b Update 1).  Natick, Massachusetts: The MathWorks Inc., 2022.
  * [32] A. Molinari and C. Daraio, “Stationary shocks in periodic highly nonlinear granular chains,” Physical Review E, vol. 80, no. 5, p. 056602, 2009.