[1]Emanuele Riva

[1]Department of Mechanical Engineering, Politecnico di Milano, Via La Masa,1,
Milano, 20156, Italy

2]Department of Civil, Chemical, Environmental and Materials Engineering,
University of Bologna, Bologna, 40136, Italy

# Controlling surface acoustic waves (SAWs) via temporally graded metasurfaces

Jonatha Santini<EMAIL_ADDRESS>Xingbo Pu<EMAIL_ADDRESS>Antonio Palermo<EMAIL_ADDRESS>Francesco Braghin
<EMAIL_ADDRESS><EMAIL_ADDRESS>*  [

###### Abstract

In this manuscript, the temporal rainbow effect for surface acoustic waves
(SAW) is herein illustrated through a temporal analog of space metagradings.
We show that a time-modulated array of mechanical resonators induces a
wavenumber-preserving frequency transformation which, in turn, dictates
Rayleigh-to-Shear wave conversion. The process is unfolded through the
adiabatic theorem, which allows us to delineate the transition between a
solely frequency-converted wave packet and a temporally-driven mode
conversion. In other words, our paper explores the role of time modulation in
the context of elastic metasurfaces, and we envision our implementation to be
suitable for designing a new family of SAW devices with frequency conversion,
mode conversion, and unusual transport capabilities.

###### keywords:

Acoustic metasurface, time-varying metamaterials, metagrading, surface-to-bulk
mode conversion, rainbow effect, Surface Acoustic Waves (SAW).

##  1 Introduction

Waves, whether they take the form of sound, light, or mechanical vibrations,
serve as fundamental carriers of information and energy, thus enabling a
plethora of applications across several disciplines. Not surprisingly, the
opportunity to manipulate wave motion has had far-reaching implications in
diverse technologies that we experience every day for communication, sensing,
and isolation purposes [1]. Within this framework, metamaterials have
blossomed a groundbreaking advancement in material engineering, whereby
intricate structural designs have promoted unconventional properties and
ensuing wave propagation functions, propelling the quest for new underpinning
behaviors. See, for example, the recent advances in the context of topological
phases of matter, waveguiding, and nonreciprocal materials [1, 2].

Among the available physical platforms, we herein operate with surface
acoustic waves (SAWs), which are known to be particularly attractive for
signal-processing applications [3]. Indeed, SAW devices are frequently
integrated with solid-state platforms as filters, sensors, and delay lines [4,
5, 6] and, in addition, have been conceived to operate with quantum systems
such as superconducting circuits, optical components, and optomechanical
devices [7, 8, 9, 10].

An invaluable tool for manipulating surface waves consists of mechanical
metasurfaces, which are made of a host medium equipped with periodic
tessellations of units placed atop an otherwise homogeneous solid, including
features such as holes, rigid cylinders, spheres, or resonant arrays [11, 12].
Non-periodic configurations are also employed to provide a diverse degree of
controllability, whereby the spatial distribution of a relevant parameter is
key in wave manipulation. The concept, known as space metagrading, has been
applied in several frameworks, including studies on 1D structures with rainbow
trapping capabilities [13, 14], which is a phenomenon enabling vibration
confinement as a result of the gradual variation of wave velocity in space.
Such variation dictates a frequency-preserving wavenumber transformation and
is responsible for the formation of “slow waves” [15, 16]. Spatially varying
metasurfaces have been explored also in the realm of SAWs where the rainbow
trapping and underlying wavenumber transformation are accompanied by mode
conversion between a wave sitting on the surface and a bulk wave propagating
toward the interior of the material [17, 18, 19].

Another promising way to manipulate wave motion leverages the temporal degree
of freedom, where the modulation does not occur in space, but over time
instead. Time-modulation schemes are known to produce frequency conversion,
unidirectional waves [20, 21, 22, 23], temporal compression [24] and pumping
[25, 26], which are behaviors more easily observable in 1D media. 2D media
with time-varying properties are more intriguing and challenging but, at the
same time, allow waveguiding in a controllable manner. A notable example,
recently explored both in phononics [27] and photonics [28], consists of
tailored time-modulation schemes able to modify the poynting vector which, in
turn, promotes wave steering. Other forms of time-modulation functions include
periodic drives for mode conversion of light [29].

Motivated by the intriguing waveguiding capabilities of time-varying materials
and by the recent discoveries in the context of space and time metagradings,
we herein merge these two concepts to produce a solid for frequency and mode
conversion. We show that a time-dependent array of mechanical resonators
placed atop an elastic half-space can be used to provide different degrees of
controllability of wave motion, which depends on the modulation scheme. The
dynamics of such a system is herein studied through a temporal analog of space
metagradings which, in contrast to its spatial counterpart, dictates a
wavenumber-preserving frequency conversion during wave motion. The underlying
frequency transformation is herein unfolded by the adiabatic theorem [27, 26,
25, 16, 30], delineating the limiting modulation speed for the transformation
to occur without energy scattering, which would otherwise imply leakage toward
the other wave modes populating the solid.

In this paper, we show how a time-modulated metasurface can be tailored to
produce: (i) pure frequency conversion; (ii) frequency conversion with
scattering toward other surface modes; (iii) frequency conversion accompanied
by surface-to-bulk mode conversion; (iv) energy transport between two
metasurfaces via adiabatic passage. All these aspects will be studied in light
of the adiabatic theorem, which provides a suitable mathematical tool for the
interpretation of the underlying dynamics.

Our work provides guidelines for designing time-varying metasurfaces with
novel functionalities, such as frequency conversion, mode conversion, and
unusual transport properties, and may find new exciting applications in the
framework of SAW technologies.

The paper is organized as follows. In section 2, we delineate the transition
between adiabatic (slow) and nonadiabtic (fast) transformations for frequency
conversion. In section 3, we study the interplay between frequency and
surface-to-bulk mode conversion, with emphasis on the role of the modulation
velocity. Concluding remarks are presented in section 4.

##  2 Frequency conversion induced by slow and fast time-modulation

Figure 1: Schematic of the system. A time-varying metasurface is realized by
way of time-dependent mechanical resonators rigidly connected atop an elastic
solid which, in a high-frequency regime, emulates a half-space. Time
modulation is employed to frequency-convert and mode-convert an impinging
surface wave into a bulk wave that propagates towards the bulk of the material
with an inclination angle $\gamma$ dictated by the modulation parameters.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 2: 2(a) Analytical dispersion relation obtained upon varying the
resonator stiffness from $k_{r1}$ (continuous black lines) to $k_{r2}$ (dashed
black lines). Stars are representative of pressure and shear waves. The
vertical dashed line denotes the incident wavenumber $\kappa_{x}^{*}$ employed
(and preserved) during numerical simulations, while the red dash-dotted line
represents the frequency of the impinging wave. 2(b) Zoomed view of Fig. 2(a).
2(c) Degree of adiabaticity as a function of the modulation speeds, the
modulation stiffness, and for the impinging wavenumber $\kappa_{x}^{*}$.
Velocity and stiffness values generating a curve below the grey horizontal
plane can be considered adiabatic. The white dashed line represents a set of
parameters used in this paper for simulation purposes. 2(d)-2(e) Numerical
dispersion relations computed in a COMSOL Multiphysics environment 2(d) before
and 2(e) after the modulation takes place. The colors are proportional to the
inverse participation ratio (IPR), which denotes localized (red) and bulk
modes (blue). The inset displays the unit cell with implied periodic
conditions employed in the simulations. 2(f) Evolution of incident state
$\mathbf{\psi}_{1}\left(\kappa_{x}^{*}\right)$ and all neighboring states
$\mathbf{\psi}_{2}\left(\kappa_{x}^{*}\right),\;\mathbf{\psi}_{3}\left(\kappa_{x}^{*}\right),\dots$
as time elapses. The grey-shaded region delimits the region populated by bulk
modes. Figure 3: Simulation results for a surface wave impinging an array of
time-modulated resonators, where the modulation parameters are compliant with
the adiabatic theorem. (a - b) Comparison between the expected and numerical
dispersions obtained (a) before ($t<T_{1}$) and (b) after the modulation ends
($t>T_{1}$). (c) Frequency spectrogram obtained at constant wavenumber
$\kappa_{x}^{*}$. The black lines show the expected central frequencies before
(continuous line), during (dot-dashed line), and after (dashed line) the
modulation. We observe that all the energy remains in correspondence of the
excited state, which is the blueprint of adiabaticity. (d) Comparison between
the frequency content before and after the modulation. (e) Wavenumber
spectrogram for $\kappa_{x}$. (f) Snapshots of the wavefield before and after
time modulation. The wave packet exhibits delocalization. (g) Numerical
(black) and analytical (red) surface wave profiles before (continuous line)
and after (dashed line) the modulation takes place.

We start the discussion with the schematic illustrated in Fig. 1, which
consists of a passive solid equipped with an active metasurface specifically
engineered to interact with surface acoustic waves (SAW) propagating in the
$x-y$ plane. The interaction is achieved through a periodic arrangement of
time-modulated mechanical resonators, where the time dependence is herein
induced in the form of stiffness-modulation. Note that a mass-modulation could
be also used. To elucidate the role of time modulation, we illustrate how an
impinging wave with frequency $\omega_{1}$ is transformed over time. The
output, dictated by a temporal version of Snell’s law, is two-fold and depends
on the modulation parameters. First and foremost, this transformation
preserves the wavenumber content along the propagation direction while, due to
time modulation, the frequency changes to $\omega_{2}$ [27]. The modulation is
hence reminiscent of a knob, which is screwed to vary the frequency content of
the output. When the knob is further rotated to produce a greater output
frequency $\omega_{3}$, the process is accompanied by mode conversion from a
wave sitting on the surface to a shear wave propagating toward the interior
with an inclination $\gamma$, which is tailorable to a certain extent.

The dynamics of such a system is herein studied under the assumption of a
quasi-static evolution of the wave modes and is accomplished starting from a
dispersion analysis. Details are reported in the supplementary material (SM)
[31]. The computation is carried out following the approach reported in Ref.
[21] and the results are illustrated in Fig. 2(a). Two stiffness levels
$k_{r1}$ and $k_{r2}$ are herein considered, which correspond to the initial
and final stiffness values represented with solid and dashed lines. Stars are
representative of pressure and shear waves with characteristic speeds $c_{L}$
and $c_{S}$, respectively.

We now consider a configuration where the effective stiffness is initially
$k_{r1}=k_{0}$ and smoothly evolves to $k_{r2}=1.57\;k_{0}$ following a linear
modulation $k(t)=k_{r1}+v_{m}(t-T_{1})$, where $v_{m}$ is the velocity of the
modulation.

The energy, initially injected in the neighborhood of the dispersion branch at
$f^{*}=100$ kHz, must follow the wavenumber-preserving transformation
illustrated with a vertical dashed line in Fig. 2(a) which, to ease
visualization, is also reported in the zoomed view (Fig. 2(b)). In both
diagrams, $\kappa_{x}^{*}$ is the impinging wavenumber corresponding to
$f^{*}$ and $\kappa_{x}^{*}l/\pi$ is its dimensionless form. $l$ is the
lattice size. Provided the modulation is slow enough, the spectral content of
the wave is expected to adiabatically follow the Rayleigh wave dispersion
curve and, hence, induce a frequency transformation from the solid line to the
dashed line. It is initially assumed that the stiffness modulation is
sufficient only to induce a small shift in frequency and, hence, without
triggering any interaction between Rayleigh and shear modes, which would be
present in case the Rayleigh dispersion approaches the shear dispersion [17].

We remark that the analysis herein discussed is relative to a wave propagating
along $x$ with evanescent wavenumber $\kappa_{y}$, which determines the degree
of penetration of the wave rather than an inclined propagation toward the
bulk. We also remark that the shear and pressure waves displayed in Figs. 2(a)
and 2(b) are characterized by $\kappa_{y}=0$, since waves with $\kappa_{y}\neq
0$ would populate another plane of the dispersion cone
$\omega\left(\kappa_{x},\kappa_{y}\right)$. This concept will be further
discussed in the following section and is key for mode conversion.

To shed light on the adiabaticity of this transformation, the dispersion
analysis of a finite strip of half-plane with implied $x$-periodic conditions
is developed in a COMSOL Multiphysics environment. The dimension along $y$ is
set sufficiently large to emulate an infinite medium. To ease readability,
details on the dispersion analysis, methods, and system parameters are
reported in the supplementary material (SM) [31]. Similar to previous studies
on elastic and acoustic lattices [27, 24], the framework herein discussed
allows us to discriminate between adiabatic and nonadiabatic modulations
which, depending on the modulation velocity, may drive the energy transfer
between the incident Rayleigh wave and all other bulk and surface wave modes
populating the solid for the impinging wavenumber.

We now focus on the numerical dispersion relation shown in Fig. 2(d) and Fig.
2(e) obtained with $k_{r1}=k_{0}$ and $k_{r2}=1.57\;k_{0}$, respectively. The
diagrams are color-coded: colors are proportional to the inverse participation
ratio (IPR), which quantitatively describes mode localization and is defined
as:

$IPR=\frac{\int_{S}\;\left<\mathbf{\psi}_{n}|\mathbf{\psi}_{n}\right>^{2}dS}{\left(\int_{S}\left<\mathbf{\psi}_{n}|\mathbf{\psi}_{n}\right>dS\right)^{2}}$
(1)

where $S$ is the unit cell domain. $\left|\mathbf{\psi}_{n}\right>$ is the
eigenvector relative to the eigenvalue $\omega_{n}$ evaluated through the
reduced eigenvalue problem
$j\omega\left|\mathbf{\psi}_{n}\right>=H\left(\kappa_{x},t\right)\left|\mathbf{\psi}_{n}\right>$,
where the Hamiltonian matrix $H\left(\kappa_{x}^{*},t\right)$ quasi-statically
varies over time. Red dots are representative of localized modes, e.g. a
surface wave, while blue dots denote a bulk wave. Note that the hybridized
Rayleigh wave is highly localized near resonance, while the upper branch
becomes less localized as the Rayleigh dispersion approaches the bulk one.
This is key for the following discussion.

We now discuss how the dispersion of the strip for the impinging wavenumber
$\kappa_{x}^{*}$ varies during time modulation. Results are reported in Fig.
2(f). The evolution of all wave modes is shown, including those populating the
bulk and herein delimited within the grey area. We observe that these modes do
not change their frequency and localization properties. In contrast, Rayleigh
waves experience a frequency transformation accompanied by a delocalization
mechanism whereby the degree of penetration changes over time. In this
section, we stop the knob before the surface mode enters the grey region. This
limitation will be later relaxed to show how a larger stiffness modulation
allows driving a conversion mechanism between the surface wave and the bulk
modes.

We remark that Fig. 2(f) is representative of a quasi-static (adiabatic)
transformation of an incident state, i.e. in case the modulation from $k_{r1}$
to $k_{r2}$ is sufficiently slow for the energy not to be scattered towards
other wave modes. We now assess the adiabaticity of this transformation under
the assumption of an arbitrary modulation speed and by following the procedure
in Ref. [27, 26, 16], which requires the elastodynamic equation to be written
in terms of first-order differential form akin to Schrödinger’s equation for
an impinging wavenumber $\kappa_{x}^{*}$:

$\left|\mathbf{\Psi}_{n}\right>_{\;,t}=H\left(\kappa_{x}^{*},t\right)\left|\mathbf{\Psi}_{n}\right>$
(2)

The energy, initially injected as a Rayleigh wave, is represented as
$\left|\mathbf{\Psi}_{1}\right>$ at $t=T_{1}$ in Fig. 2(f). The
transformation, which takes place over time, can be considered adiabatic if
there is negligible energy leakage toward the neighboring states
$\left|\mathbf{\Psi}_{2}\right>$, $\omega_{2}$ or, alternatively,
$\left|\mathbf{\Psi}_{3}\right>$, $\omega_{3}$. We found
$\left|\mathbf{\Psi}_{2}\right>$, $\omega_{2}$ more critical from the
adiabaticity viewpoint, which is justified by the smaller frequency separation
$\omega_{2}-\omega_{1}$. For adiabaticity, the following inequality must be
satisfied:

$p\;v_{m}<<1\hskip 28.45274ptwhere:\hskip
14.22636ptp=\left|\frac{\Big{<}\mathbf{\psi}_{2}^{L}\left|\displaystyle\frac{\partial
H}{\partial
k_{r}}\right|\mathbf{\psi}_{1}^{R}\Big{>}}{\left(\omega_{2}-\omega_{1}\right)^{2}}\right|$
(3)

The term on the left-hand side is computed for the impinging wavenumber
$\kappa_{x}^{*}$ and mapped in Figure 2(c) as a function of time and upon
varying the modulation speed $v_{m}$. For sufficiently small modulation
speeds, the transformation can be considered adiabatic, i.e. when the surface
is below the unitary threshold highlighted by the grey plane. In the diagram,
the white-dashed curve is relative to the adiabatic modulation employed for
simulation purposes and hereafter discussed in detail.

The theoretical discussion is now verified through numerical simulations. To
delineate the transition between adiabatic and non-adiabatic modulations, two
scenarios with the same starting and ending stiffness values, but different
modulation speeds, are reported. Fig. 3 illustrates results compliant with the
adiabatic limit, whereby the velocity of modulation is low. In contrast, Fig.
4, illustrates the results with a step-like modulation, which can be
considered nonadiabatic.

Figure 4: Simulation results for a surface wave impinging on an array of
resonators modulated with a non-adiabatic (very fast) transformation in time.
(a)-(b) Numerical dispersion before and after the step modulation centered at
$\tau=T$. (c) Frequency spectrogram at constant wavenumber $\kappa_{x}^{*}$.
Black lines show the expected central frequency before (continuous line) and
after (dashed line) the modulation. The energy, initially injected in
correspondence of the excited state, leaks toward all the states that populate
the dispersion relation after modulation. (d) Comparison between the frequency
content before and after the modulation. In the inset, a zoomed view of the
shaded area is used to highlight leakage towards the lower branch of the
dispersion. (e) Wavenumber spectrogram for $\kappa_{x}$. (f) Wavefield before
and after time modulation. The wave packet exhibits delocalization. An
additional mode, caused by the non-adiabaticity of the process, is present on
the surface. (g): Numerical (black) and analytical (red) wave profiles before
(continuous line) and after (dashed line) the modulation.

Two representative snapshots before and after time modulation are reported in
Fig. 3(f). As expected, since the modulation is adiabatic, the wave packet
experiences a transformation without energy leak toward other modes populating
the surface and the bulk material. As a result, a single wave packet is
obtained as output, which is however characterized by a greater degree of
penetration in the bulk of the material, as anticipated through the
numerically computed IPR in Figs. 2(f), 2(d) and 2(e). The analytical (red)
and numerical (black) envelopes before (continuous line) and after (dashed
line) the modulation are also reported in Fig. 3(g) for completeness.

The time history is further processed to shed more light on the results, as
shown in Fig. 3(a)-(b). The figure compares the analytical and numerical
(reconstructed from simulation) dispersions of the medium in reciprocal space
before ($t<T_{1}$) and after ($t>T_{2}$) the modulation takes place. The
numerical dispersion is computed by taking the 3D Fourier transform of wave
motion $\bm{u}\left(x,y,t\right)$ measured within the solid, which yields the
complex-valued coefficients
$\bm{\hat{u}}\left(\kappa_{x},\kappa_{y},\omega\right)$. The wavenumber
dependence $\kappa_{y}$ is eliminated by taking the RMS value along
$\kappa_{y}$, leading to the dispersion
$\bm{\hat{u}}\left(\kappa_{x},\omega\right)$ displayed as a contour plot. In
the first diagram, the energy content follows the solid line, which is
representative of the dispersion with $k_{r1}$. In the second diagram, the
energy content matches the dashed line, which is relative to the final
modulation value $k_{r2}$. Note that the amplitude is zero elsewhere and,
hence, we can conclude that a scattering-free frequency conversion is
achieved.

The frequency shift is more visible in the spectrogram reported in Fig. 3(c),
which is obtained by windowing the time history with a moving Gaussian
function $e^{-\left(t-t_{0}\right)^{2}/2c^{2}}$. Here, the time history is
Fourier-transformed to
$\hat{u}\left(\kappa_{x},\kappa_{y},\omega,t_{0}\right)$, and reduced to
$\hat{u}\left(\omega,t_{0}\right)$ by taking the RMS value along both
wavenumber dimensions. Also, the spectrum $\hat{u}\left(\omega,t<T1\right)$
and $\hat{u}\left(\omega,t>T2\right)$, evaluated before and after time
modulation, is reported alongside the spectrogram (see Fig. 3(d)), which
further shows the aforementioned frequency shift. We observe that, due to the
adiabatic nature of the modulation, the incident wave evolves sufficiently
slow for the transformation to avoid scattering toward the other available
wave modes, and, as such, the energy remains successfully located in the upper
state. Finally, the spectrogram reported in Fig. 3-V demonstrates the
conservation of the wavenumber content, which is preserved during time
modulation.

The outcome of a non-adiabatic transformation, obtained through a very fast
modulation of the stiffness, is displayed in Fig. 4. Interestingly, the
envelopes before and after time modulation (Fig. 4(g)) are very similar to the
adiabatic case, but there is an additional wave packet that nucleates due to
the temporal modulation, which is visible in Fig. 4(f) for $t>T$. To shed
light on the nature of this additional wave packet and the underlying
dynamics, the spectral content of the simulation is reported in Figs.4(a)-(c).
In contrast to the adiabatic case, we observe a bit of energy leak toward the
resonant branch of the dispersion. This is particularly visible in the
frequency spectrogram in Fig.4(c), where the frequency conversion is
accompanied by an energy jump, also present in Fig.4(d). This energy jump is
justified by the additional wave packet formed across the time interface.

##  3 Mode conversion induced by slow and fast time modulation

(a)

(b)

Figure 5: 5(a) Dispersion relation upon varying the resonator stiffness from
$k_{r1}$ (continuous black lines) to $k_{r3}$ (dashed black lines). White
stars represent the pressure and shear modes. The vertical dashed line
represents the incident wavenumber $\kappa_{x}^{*}$, which is preserved during
the modulation. 5(b) Zoomed view of Fig. 5(a) limited to the grey highlighted
area. The arrow illustrates that the wavenumber-preserving transformation is
tailored to frequency-convert the impinging energy above the shear cone.
Figure 6: Simulation results for a surface wave impinging on an array of
resonators modulated with a non-adiabatic (very fast) transformation in time.
(a)-(b) Numerical dispersion relation before and after the step modulation at
time $t=T$. (c) Frequency spectrogram at constant wavenumber $\kappa_{x}^{*}$.
Black lines show the expected central frequencies of the surface waves before
(continuous line) and after (dashed line) the modulation. White stars
represent the output frequency of the shear mode for the impinging wavenumber
$\kappa_{x}^{*}$ achieved after the modulation. (d) Comparison between the
frequency content before and after the modulation. The inset highlights
leakage towards the lower branch of the dispersion. (e)-(f) Wavenumber
spectrogram for (e) $\kappa_{x}$ and (f) $\kappa_{y}$. (g) Wavefield before
and after the modulation. The wave packet exhibits surface-to-bulk mode
conversion and additional modes nucleate on the surface due to the non-
adiabatic nature of the process. Figure 7: (a) Orientation of the planes
hosting the dispersion in case of mode conversion. The grey cone represents
shear modes and the blue plane hosts the Rayleigh dispersion.
$\kappa_{\gamma}-\omega$ is the plane after the modulation takes place. Waves
populating the $\kappa_{\gamma}-\omega$ plane point towards the bulk of the
material with an inclination $\gamma$. (b): Numerically reconstructed
wavenumber content before and after the modulation. There is a visible change
of $\kappa_{y}$ that occurs over time and is responsible for a non-null
$\gamma$. (c) Snapshots of the wavefield for $t>T$ evaluated for different
stiffness modulation values $k_{r}$, (I) $k_{r}=1.95\;k_{0}$, (II)
$k_{r}=2.40\;k_{0}$, (III) $k_{r}=2.91\;k_{0}$, and (IV) $k_{r}=3.25\;k_{0}$.
The angle at the outlet varies with $k_{r}$.

We now further rotate the knob to achieve mode conversion, whereby the
stiffness is modulated to shift the upper dispersion branch above the shear
cone for the impinging wavenumber $\kappa_{x}^{*}$, see Fig. 5. In contrast to
prior works on space metagradings, we remark that a gradual variation of the
properties would slowly transform the surface mode into a shear wave
populating the $\kappa_{x}-\omega$ plane, hence, without any directionality
$\gamma$ toward the interior, which would be dictated by non-evanescent and
non-null values of $\kappa_{y}$. This behavior is justified by the adiabatic
nature of the transformation, which implies energy conversion into a bulk mode
before a further increase of the frequency content takes place. This
assumption is numerically verified through a slow modulation from
$k_{r1}=k_{0}$ to $k_{r3}=2.4\;k_{0}$, thereby achieving full conversion
between a surface wave to a shear wave with $\kappa_{y}=0$ and, hence, with a
null propagation angle $\gamma=0$. Numerical results on this matter are
reported in the supplementary material (SM) [31].

A nonadiabatic step modulation is hence necessary to suddenly change the
frequency content $\omega_{1}\rightarrow\omega_{3}$ of the surface wave above
the shear cone which, due to momentum conservation along $x$, would determine
a complex-valued wavenumber
$\kappa_{y}=\sqrt{{\kappa_{x}^{*}}^{2}-\omega_{3}^{2}/c_{S}^{2}}$ at the
outlet. In contrast to adiabatic modulations with mode conversion, the energy
sitting on the surface is altered by the time duration of the drive, which is
very quick. Hence, we expect the frequency conversion illustrated by the arrow
in Fig. 5(b) to be accompanied by undesired mode scattering towards the
available states that populate the solid. The stiffness levels considered
herein are again $k_{r1}=k_{0}$ and $k_{r3}=2.4\;k_{0}$, corresponding to the
solid and dashed line in the figure, respectively, while stars are
representative of shear and pressure waves.

The concept is corroborated via numerical simulations. In contrast to the
previous results, three wave modes are visible in Fig. 6(g) for $t>T$, i.e.
after the modulation takes place, which consist of two counter-propagating
surface waves and a shear wave. As expected, the impinging wave is converted
into a shear wave in response to a frequency transformation
$\omega_{1}\rightarrow\omega_{3}$ and the non-adiabatic nature of the
modulation is revealed by the nucleation of additional wave packets on the
surface with minor energy content.

The frequency shift and the formation of such waves are better illustrated in
the frequency domain, see Fig. 6(a)-(c). First, we note that before and after
the modulation, the numerical and theoretical dispersion are in good
agreement. However, while Fig. 6(a) displays results in the
$\kappa_{x}-\omega$ plane, Fig. 6(b) illustrates that the energy content lies
in the $\kappa_{\gamma}-\omega$ plane, where $\gamma$ is coincident with the
expected propagation angle
$\gamma=\arctan{\left(\kappa_{y}/\kappa_{x}\right)}$. This means that the
frequency up-conversion is accompanied by a wavenumber transformation for
$\kappa_{y}$. To ease visualization, the two planes are reported schematically
in Fig. 7(a). In addition, due to the nonadiabatic nature of the modulation,
both the frequency spectrogram and spectrum before and after modulation in
Fig. 6(c)-(d) not only report frequency up-conversion, but also energy leakage
towards the resonant dispersion branch. This justifies the presence of
additional surface modes trapped within the array of resonators.

Now, the time history is processed to show the wavenumber content in
reciprocal space and relative spectrograms. The resulting wavenumber
spectrograms in Fig. 6(e)-(f) further demonstrate momentum conservation for
$\kappa_{x}^{*}$, while there is an evolution of wavenumber $\kappa_{y}$
toward a nonzero value. This de facto determines a nonzero angle $\gamma$
visible in the $\kappa_{x}-\kappa_{y}$ plane, see Fig. 7(b), which better
highlights the inclination angle achieved in the numerical simulation after
the modulation takes place. We remark that further increasing the final
stiffness value would modify the wavenumber
$\kappa_{y}=\sqrt{{\kappa_{x}^{*}}^{2}-\omega_{3}^{2}/c_{S}^{2}}$ and, hence,
the angle $\gamma$. This is shown through the temporal snapshots illustrated
in Fig. 7(c), which are obtained with different stiffness modulation values.
To ease readability, the spectral analysis is left in the supplementary
material (SM) [31].

To further illustrate the functionalities of the metasurface, we now employ a
smooth modulation that drives the energy transfer from the upper to the lower
surface of a solid with finite dimensions along $y$. In here, a longer solid
is considered, which is necessary to induce a slow transition between surface
and bulk waves. Results are reported in Fig. 8(a) in terms of snapshots
before, during, and after the modulation. To achieve this effect, a slow
temporal drive is applied to the top and bottom layer $k_{r}=k_{0}\rightarrow
2.4\;k_{0}$, in a way that two coincident surface states belonging to the
opposite edges form and, later, delocalize into a shear wave. In other words,
due to an adiabatic mechanism, the surface wave slowly transforms into a shear
wave. Then, once the shear mode is formed, the temporal drive
$k_{r}=2.4\;k_{0}\rightarrow k_{0}$ is applied only to the bottom metasurface,
thereby achieving localization in the form of a surface wave. The behavior is
well summarized by the RMS of the entire time history reported in Fig. 8(b),
which further confirms the edge-to-edge energy transfer. The corresponding
spectrogram with superimposed evolution of the states is displayed in Fig.
8(c). We observe that during the entire duration of the process, the energy,
initially injected in the surface, follows the expected evolution of the
surface states which, later, temporarily merge with the bulk. When the
temporal drive is applied to the bottom surface, such a mode back-transforms
into a surface wave traveling along the opposite surface, thereby achieving an
edge-to-edge energy transfer.

(a)

(b)

(c)

Figure 8: Results for a temporally driven edge-to-edge mode conversion, which
is achieved through a double array of modulated resonators. 8(a) Three
relevant snapshots before, during, and after time modulation. 8(b) RMS of the
time history. 8(c) Comparison between numerically reconstructed (colors) and
expected (dots) evolution of the modes. There is a mode that migrates toward
the bulk band, which is representative of a surface wave that transforms first
into a bulk mode and later back-transforms into a surface wave.

##  4 Conclusions

This work illustrates the dynamics of time-varying metasurfaces. The
implementation, which relies on time-modulated resonators, induces a
wavenumber-preserving frequency transformation which is suitable to manipulate
wave motion. The role of the modulation parameters in the frequency and mode
conversion capabilities of time-varying metasurfaces are critically discussed
and several functionalities inherently linked with the velocity and strength
of the modulation are revealed. We have shown: pure frequency conversion,
frequency conversion with energy scattering, concurrent frequency and mode
conversion, and surface-to-surface energy transfer across the solid. The
transition between these functionalities is delineated by the adiabatic
theorem, which dictates a condition for the transformation to occur with or
without scattering toward the wave modes populating the solid.

We envision this new family of metasurfaces to suit the design of SAW devices
with unusual transport capabilities, such as frequency conversion, mode
conversion, and edge-to-edge transport.

See pages - of Supplementary.pdf

## References

  *   * Oudich et al. [2023] Oudich, M., Gerard, N.J., Deng, Y., Jing, Y.: Tailoring structure-borne sound through bandgap engineering in phononic crystals and metamaterials: A comprehensive review. Advanced Functional Materials 33(2), 2206309 (2023)
  * Nassar et al. [2020] Nassar, H., Yousefzadeh, B., Fleury, R., Ruzzene, M., Alù, A., Daraio, C., Norris, A.N., Huang, G., Haberman, M.R.: Nonreciprocity in acoustic and elastic materials. Nature Reviews Materials 5(9), 667–685 (2020)
  * Campbell [1991] Campbell, C.: Surface acoustic wave devices and their signal processing applications. Acoustical Society of America (1991)
  * Shao et al. [2020] Shao, L., Mao, W., Maity, S., Sinclair, N., Hu, Y., Yang, L., Lončar, M.: Non-reciprocal transmission of microwave acoustic waves in nonlinear parity–time symmetric resonators. Nature Electronics 3(5), 267–272 (2020)
  * Ferreira et al. [2009] Ferreira, G.N., Da-Silva, A.-C., Tomé, B.: Acoustic wave biosensors: physical models and biological applications of quartz crystal microbalance. Trends in biotechnology 27(12), 689–697 (2009)
  * Weigel et al. [2002] Weigel, R., Morgan, D.P., Owens, J.M., Ballato, A., Lakin, K.M., Hashimoto, K.-Y., Ruppel, C.C.: Microwave acoustic materials, devices, and applications. IEEE Transactions on microwave theory and techniques 50(3), 738–749 (2002)
  * Manenti et al. [2017] Manenti, R., Kockum, A.F., Patterson, A., Behrle, T., Rahamim, J., Tancredi, G., Nori, F., Leek, P.J.: Circuit quantum acoustodynamics with surface acoustic waves. Nature communications 8(1), 975 (2017)
  * Forsch et al. [2020] Forsch, M., Stockill, R., Wallucks, A., Marinković, I., Gärtner, C., Norte, R.A., Otten, F., Fiore, A., Srinivasan, K., Gröblacher, S.: Microwave-to-optics conversion using a mechanical oscillator in its quantum ground state. Nature Physics 16(1), 69–74 (2020)
  * Navarro-Urrios et al. [2022] Navarro-Urrios, D., Colombano, M.F., Arregui, G., Madiot, G., Pitanti, A., Griol, A., Makkonen, T., Ahopelto, J., Sotomayor-Torres, C.M., Martínez, A.: Room-temperature silicon platform for ghz-frequency nanoelectro-opto-mechanical systems. ACS photonics 9(2), 413–419 (2022)
  * Satzinger et al. [2018] Satzinger, K.J., Zhong, Y., Chang, H.-S., Peairs, G.A., Bienfait, A., Chou, M.-H., Cleland, A., Conner, C.R., Dumur, É., Grebel, J., et al.: Quantum control of surface acoustic-wave phonons. Nature 563(7733), 661–665 (2018)
  * Palermo et al. [2016] Palermo, A., Krödel, S., Marzani, A., Daraio, C.: Engineered metabarrier as shield from seismic surface waves. Scientific reports 6(1), 39356 (2016)
  * Zeighami et al. [2021] Zeighami, F., Palermo, A., Marzani, A.: Rayleigh waves in locally resonant metamaterials. International Journal of Mechanical Sciences 195, 106250 (2021)
  * Santini et al. [2022] Santini, J., Sugino, C., Riva, E., Erturk, A.: Harnessing rainbow trapping via hybrid electromechanical metastructures for enhanced energy harvesting and vibration attenuation. Journal of Applied Physics 132(6) (2022)
  * Alshaqaq and Erturk [2020] Alshaqaq, M., Erturk, A.: Graded multifunctional piezoelectric metastructures for wideband vibration attenuation and energy harvesting. Smart Materials and Structures 30(1), 015029 (2020)
  * De Ponti et al. [2021] De Ponti, J.M., Iorio, L., Riva, E., Ardito, R., Braghin, F., Corigliano, A.: Selective mode conversion and rainbow trapping via graded elastic waveguides. Physical Review Applied 16(3), 034028 (2021)
  * Riva et al. [2023] Riva, E., De Ponti, J., Marconi, J., Braghin, F., Ardito, R., Corigliano, A.: Adiabatic limit for scattering-free waveguiding in space-graded arrays of micro-resonators. European Journal of Mechanics-A/Solids, 105162 (2023)
  * Colombi et al. [2016] Colombi, A., Colquitt, D., Roux, P., Guenneau, S., Craster, R.V.: A seismic metamaterial: The resonant metawedge. Scientific reports 6(1), 27717 (2016)
  * Colombi et al. [2017] Colombi, A., Ageeva, V., Smith, R.J., Clare, A., Patel, R., Clark, M., Colquitt, D., Roux, P., Guenneau, S., Craster, R.V.: Enhanced sensing and conversion of ultrasonic rayleigh waves by elastic metasurfaces. Scientific Reports 7(1), 6750 (2017)
  * Alan et al. [2019] Alan, S., Allam, A., Erturk, A.: Programmable mode conversion and bandgap formation for surface acoustic waves using piezoelectric metamaterials. Applied Physics Letters 115(9) (2019)
  * Marconi et al. [2020] Marconi, J., Riva, E., Di Ronco, M., Cazzulani, G., Braghin, F., Ruzzene, M.: Experimental observation of nonreciprocal band gaps in a space-time-modulated beam using a shunted piezoelectric array. Physical Review Applied 13(3), 031001 (2020)
  * Palermo et al. [2020] Palermo, A., Celli, P., Yousefzadeh, B., Daraio, C., Marzani, A.: Surface wave non-reciprocity via time-modulated metamaterials. Journal of the Mechanics and Physics of Solids 145, 104181 (2020)
  * Attarzadeh et al. [2020] Attarzadeh, M., Callanan, J., Nouh, M.: Experimental observation of nonreciprocal waves in a resonant metamaterial beam. Physical Review Applied 13(2), 021001 (2020)
  * Wu et al. [2021] Wu, Q., Chen, H., Nassar, H., Huang, G.: Non-reciprocal rayleigh wave propagation in space–time modulated surface. Journal of the Mechanics and Physics of Solids 146, 104196 (2021)
  * Riva et al. [2023] Riva, E., Rosa, M.I., Guo, Y., Ruzzene, M., et al.: Adiabatic sound transport in acoustic waveguides with time-varying helmholtz resonators. Frontiers in Acoustics 1, 1–8 (2023)
  * Riva et al. [2021] Riva, E., Castaldini, G., Braghin, F.: Adiabatic edge-to-edge transformations in time-modulated elastic lattices and non-hermitian shortcuts. New Journal of Physics 23(9), 093008 (2021)
  * Xia et al. [2021] Xia, Y., Riva, E., Rosa, M.I., Cazzulani, G., Erturk, A., Braghin, F., Ruzzene, M.: Experimental observation of temporal pumping in electromechanical waveguides. Physical Review Letters 126(9), 095501 (2021)
  * Santini and Riva [2023] Santini, J., Riva, E.: Elastic temporal waveguiding. New Journal of Physics 25(1), 013031 (2023)
  * Pacheco-Peña and Engheta [2020] Pacheco-Peña, V., Engheta, N.: Temporal aiming. Light: Science & Applications 9(1), 129 (2020)
  * Galiffi et al. [2020] Galiffi, E., Wang, Y.-T., Lim, Z., Pendry, J.B., Alù, A., Huidobro, P.A.: Wood anomalies and surface-wave excitation with a time grating. Physical Review Letters 125(12), 127403 (2020)
  * Bösch et al. [2023] Bösch, C., Fichtner, A., Garcia, M.S.: Differences between quantum and classical adiabatic evolution. arXiv preprint arXiv:2309.08510 (2023)
  * [31] See supplemental material at xxxx for more details on the simulation procedures and additional results.


# Reconfigurable Stochastic Neurons Based on Strain Engineered Low Barrier
Nanomagnets

Rahnuma Rahman, Samiran Ganguly and Supriyo Bandyopadhyay  Department of
Electrical and Computer Engineering, Virginia Commonwealth University,
Richmond, VA 23284, USA<EMAIL_ADDRESS>
###### Abstract

Stochastic neurons are efficient hardware accelerators for solving a large
variety of combinatorial optimization problems. “Binary” stochastic neurons
(BSN) are those whose states fluctuate randomly between two levels +1 and -1,
with the probability of being in either level determined by an external bias.
“Analog” stochastic neurons (ASNs), in contrast, can assume any state between
the two levels randomly (hence “analog”) and can perform analog signal
processing. They may be leveraged for such tasks as temporal sequence
learning, processing and prediction. Both BSNs and ASNs can be used to build
efficient and scalable neural networks. Both can be implemented with low
(potential energy) barrier nanomagnets (LBMs) whose random magnetization
orientations encode the binary or analog state variables. The difference
between them is that the potential energy barrier in a BSN LBM, albeit low, is
much higher than that in an ASN LBM. As a result, a BSN LBM has a clear double
well potential profile, which makes its magnetization orientation assume one
of two orientations at any time, resulting in the binary behavior. ASN
nanomagnets, on the other hand, hardly have any energy barrier at all and
hence lack the double well feature. That makes their magnetizations fluctuate
in an analog fashion. Hence, one can reconfigure an ASN to a BSN, and vice-
versa, by simply raising and lowering the energy barrier. If the LBM is
magnetostrictive, then this can be done with local (electrically generated)
strain. Such a reconfiguration capability heralds a powerful field
programmable architecture for a p-computer, and the energy cost for this type
of reconfiguration is miniscule. There are also other applications of strain
engineered barrier control, e.g., adaptive annealing in energy minimization
computing (Boltzmann or Ising machines), emulating memory hierarchy in a
dynamically reconfigurable fashion, and control over belief uncertainty in
analog stochastic neurons. Here, we present a study of this modality.

Keywords: binary stochastic neurons, analog stochastic neurons,
reconfigurability, low barrier nanomagnets, magnetostriction, strain

##  1 Introduction: Binary and analog stochastic neurons implemented with LBMs

Binary stochastic neurons (BSNs) are a well known route to implementing
“spins” in Ising machines [1, 2] and have been used to solve computationally
hard problems such as graph theoretic problems like Max Cut [3], factorization
[4], etc. very efficiently. A popular approach to realizing them is with a low
(energy) barrier nanomagnet (LBM), whose magnetization fluctuates randomly
between two preferred orientations representing the binary bits +1 and -1. The
probability of being in either bit can be altered by a “bias”, such as a spin-
polarized current injected into the LBM [1].

The LBM is usually a nanomagnet with in-plane anisotropy that is shaped like
an elliptical disk with small (but non-zero) eccentricity. The in-plane
potential energy profile (energy versus magnetization orientation) of such a
LBM is shown schematically in Fig. 1(a). Normally, there is a clear double-
well feature which can be discerned despite the low potential barrier. The two
ground states (or wells) correspond to the magnetization pointing along either
direction along the major axis (or easy axis) of the elliptical nanomagnet. At
room temperature, thermal energy can allow the magnetization to transcend the
energy barrier separating the wells, which will allow the magnetization to
fluctuate randomly between the two potential wells. If we take a snapshot in
time, we will usually find the magnetization in one of the two wells, i.e., it
will tend to point along one of the two directions along the major axis, which
encode the bits +1 and -1. This leads to the digital or “binary” behavior.

Figure 1: (a) Potential energy as a function of the in-plane magnetization
orientation in a magnetostrictive nanomagnet shaped like an elliptical disk
with small eccentricity. (b) Potential energy as a function of the in-plane
magnetization orientation when the nanomagnet is subjected to uniaxial stress
along the major axis such that the sign of the product of the stress and the
magnetostriction is negative. The inset shows the nanomagnet and the
magnetization orientation, with $\theta$ being the angle subtended by the
magnetization with the nanomagnet’s major axis.

One can depress the energy barrier separating the two wells with (electrically
generated) mechanical strain if the nanomagnet is made of a magnetostrictive
material like Co or FeGa or Terfenol-D. For this to happen, the sign of the
product of the magnetostriction and the strain has to be negative. When the
energy barrier is depressed sufficiently, it begins to lose the double-well
feature, as shown in Fig. 1(b). At that point, the magnetization has very
little tendency to settle into either of the two degenerate ground states
preferentially, and all orientations are almost equally likely. Consequently,
the magnetization will fluctuate among all orientations with the same
likelihood, leading to the “analog” behavior. We can, therefore, reconfigure a
stochastic neuron from binary to analog, and vice-versa, with strain.

##  2 Strain effects on energy barrier

The steady-state in-plane potential energy in a strained nanomagnet with in-
plane anisotropy depends on the magnetization orientation as [5]

$E=\left(\mu_{0}/2\right)M_{s}^{2}\Omega\left[N_{1}cos^{2}\theta+N_{2}sin^{2}\theta\right]-(3/2)\lambda_{s}Y\epsilon\Omega
cos^{2}\theta,$ (1)

where $\mu_{0}$ is the permeability of vacuum, $M_{s}$ is the saturation
magnetization of the nanomagnet’s material, $\Omega$ is the nanomagnet volume,
$\lambda_{s}$ is the saturation magnetostriction of the nanomagnet material,
$Y$ is the Young’s modulus of the nanomagnet, $\epsilon$ is the strain,
$\theta$ is the angle shown in the inset of Fig. 1 to denote the magnetization
orientation (it is the angle subtended by the in-plane component of the
magnetization with the major axis of the elliptical nanomagnet), and

$\displaystyle N_{1}$ $\displaystyle=$
$\displaystyle{{\pi\over{4}}}\left({{t}\over{a}}\right)\left[1-\frac{1}{4}\left(\frac{a-b}{a}\right)-\frac{3}{16}\left(\frac{a-b}{a}\right)^{2}\right]$
$\displaystyle N_{2}$ $\displaystyle=$
$\displaystyle{{\pi\over{4}}}\left({{t}\over{a}}\right)\left[1+\frac{5}{4}\left(\frac{a-b}{a}\right)+\frac{21}{16}\left(\frac{a-b}{a}\right)^{2}\right],$

where $a$ is the major axis, $b$ is the minor axis and $t$ is the thickness of
the nanomagnet [6].

Equation (1) clearly shows that if the product $\lambda_{s}\epsilon$ has a
negative sign, then application of strain will depress the energy barrier,
which is what we depicted schematically in Fig. 1. A material like Co has
negative magnetostriction and hence a tensile uniaxial strain along the major
axis of the nanomagnet will depress the energy barrier. A material like FeGa
or Terfenol-D has positive magnetostriction and hence a compressive uniaxial
strain along the major axis will depress the energy barrier. Fig. 2 shows the
potential energy $E-E_{min}$ ($E_{min}$ is the minimum value of $E$) as a
function of the magnetization orientation in a Co nanomagnet of major axis 100
nm, minor axis 99 nm and thickness 5 nm. The energy barrier is the maximum of
this quantity, and can be seen to decrease with increasing stress.

Figure 2: The potential energy profile in a Co nanomagnet shaped like an
elliptical disk as function of the angle $\theta$ subtended by the
magnetization with the major axis. The results are shown for different stress
values. The quantity $E$ is calculated from Equation (1) and $E_{min}$ is the
minimum value of $E$. The nanomagnet has major axis = 100 nm, minor axis = 99
nm and thickness = 5 nm.  Figure 3: Methodology to reconfigure a BSN into an
ASN and vice versa. Applying a gate voltage of the right polarity will
generate the right type of biaxial strain in the poled piezoelectric region
underneath the nanomagnet. This strain will be transferred to the nanomagnet
and it will lower the energy barrier in the latter, making its magnetization
fluctuate in an analog manner rather than a binary manner. We can build a
magnetic tunnel junction (MTJ) on top of the nanomagnet which will act as its
soft layer. This MTJ will be a fluctuating resistor that can be transformed
from a BSN to an ASN by turning on the gate voltage to generate strain.

One way to strain a nanomagnet electrically is to use the configuration shown
is Fig. 3. The nanomagnet is deposited on a poled piezoelectric film and gate
pads are delineated around it such that the line joining the pads passes
through the major axis. The two pads are shorted together and a voltage is
applied between the shorted pads and the grounded conducting substrate. The
substrate is made ‘conducting’ so that the applied gate voltage drops mostly
across the piezoelectric layer and not the substrate. If the resulting
electric field is parallel to the direction of the poling, then tensile strain
will appear along the major axis and compressive along the minor axis of the
elliptical nanomagnet. Reversing the polarity of the gate voltage will reverse
the signs of the strains. The dimensions of the nanomagnet and the electrodes,
the separation between the nanomagnet edge and the nearest electrode, and the
piezoelectric film thickness have to satisfy certain conditions for the
biaxial strain generation as described [7], but these conditions are
relatively easy to fulfill.

We can build a magnetic tunnel junction (MTJ) on top of the nanomagnet, with
the latter acting as the soft layer whose magnetization fluctuates. This will
transform the MTJ into a fluctuating resistance that acts as either a BSN (no
gate voltage applied to cause strain and lower the energy barrier) or an ASN
(gate voltage of the right polarity applied to cause strain that lowers the
energy barrier). This is the basis of a reconfigurable stochastic neuron
(RSN).

An important question now is whether the strain generated by the gate voltage
can be non-volatile. This will allow the reconfiguration to be non-volatile as
well. There are many reports of non-volatile remanent strain in piezoelectrics
at room temperature [8, 9, 10, 11, 12, 13, 14] although the strain’s longevity
has not been studied. If the strain remains non-volatile, we can reconfigure a
BSN to an ASN and the reconfiguration will survive subsequent removal of the
gate voltage. To revert the ASN back to a BSN, we can simply apply strain of
the opposite sign, which will raise the energy barrier back in the nanomagnet
and convert the ASN to a BSN.

##  3 Landau-Lifshitz-Gilbert simulations to study random magnetization
dynamics in an LBM under different strains

We carried out Landau-Lifshitz-Gilbert (LLG) simulations of the magnetization
dynamics in an LBM at room temperature under different strains to see how the
magnetization fluctuation behaves. The LBM we studied is an elliptical Co
nanomagnet of major axis 100 nm, minor axis 99 nm and thickness 5 nm. A
nanomagnet of these dimensions are likely to be monodomain and hence the
macrospin approximation holds. The saturation magnetization $M_{s}$ = 106 A/m,
the magnetostriction coefficient $\lambda_{s}$ = -35 ppm and the Gilbert
damping coefficient $\alpha$ = 0.01 correspond to a Co nanomagnet. The coupled
LLG equations governing the temporal evolutions of the scalar components of
the magnetization were solved with finite difference method [15, 16] with a
time step of 0.1 ps. We assumed positive (tensile) stress since Co has
negative magnetostriction. The initial condition was that the magnetization
was aligned along the minor axis of the nanomagnet.

The coupled LLG equations describing the temporal evolution of the three
components of the magnetization are:

$\displaystyle{{dm_{x}(t)}\over{dt}}$ $\displaystyle=$
$\displaystyle-\gamma\left[H_{z}(t)m_{y}(t)-H_{y}(t)m_{z}(t)\right]$
$\displaystyle-\alpha\gamma\left[H_{y}(t)m_{x}(t)m_{y}(t)-H_{x}(t)m_{y}^{2}(t)-H_{x}(t)m_{z}^{2}(t)+H_{z}(t)m_{x}(t)m_{z}(t)\right]$
$\displaystyle{{dm_{y}(t)}\over{dt}}$ $\displaystyle=$
$\displaystyle-\gamma\left[H_{x}(t)m_{z}(t)-H_{z}(t)m_{x}(t)\right]$
$\displaystyle-\alpha\gamma\left[H_{z}(t)m_{y}(t)m_{z}(t)-H_{y}(t)m_{z}^{2}(t)-H_{y}(t)m_{x}^{2}(t)+H_{x}(t)m_{x}(t)m_{y}(t)\right]$
$\displaystyle{{dm_{z}(t)}\over{dt}}$ $\displaystyle=$
$\displaystyle-\gamma\left[H_{y}(t)m_{x}(t)-H_{x}(t)m_{y}(t)\right]$
$\displaystyle-\alpha\gamma\left[H_{x}(t)m_{z}(t)m_{x}(t)-H_{z}(t)m_{x}^{2}(t)-H_{z}(t)m_{y}^{2}(t)+H_{y}(t)m_{y}(t)m_{z}(t)\right]$

where $\alpha$ is the Gilbert damping factor of the nanomagnet material,
$\gamma$ is the gyromagnetic factor (a constant), $m_{i}(t)$ is the $i$-th
component of the magnetization at time $t$, and $H_{i}(t)$ is the $i$-th
component of the effective magnetic field experienced by the nanomagnet at
time $t$.

The effective magnetic field components are given by

$\displaystyle H_{x}(t)$ $\displaystyle=$ $\displaystyle-
M_{s}N_{2}m_{x}(t)+h_{x}^{noise}(t)$ $\displaystyle H_{y}(t)$ $\displaystyle=$
$\displaystyle-M_{s}N_{1}m_{x}(t)+h_{y}^{noise}(t)$ $\displaystyle H_{z}(t)$
$\displaystyle=$ $\displaystyle-M_{s}N_{3}m_{x}(t)+h_{z}^{noise}(t)$ (3)

where $N_{3}$ = 1 - $N_{1}$ \- $N_{2}$ and $h_{i}^{noise}(t)=\sqrt{{{2\alpha
kT}\over{\gamma\left(1+\alpha^{2}\right)\mu_{0}M_{s}\Omega\Delta
t}}}G_{(0,1)}^{i}(t)$ with $G_{(0,1)}^{i}(t)$ ($i=x,y,z$) being three
uncorrelated Gaussians of zero mean and unit standard deviation, $\Omega$ is
the nanomagnet volume and $\Delta t$ is the attempt period which is the time
step of the simulation.

Figure 4: Temporal fluctuations in the magnetization component directed along
the major axis of the nanomagnet $\left(m_{y}\right)$ at different values of
stress. Note that the behavior gradually transitions from BSN to ASN with
increasing stress as the energy barrier within the nanomagnet is progressively
depressed.

Fig. 4 shows the time variations of the normalized magnetization component
along the major axis of the nanomagnet, i.e. $m_{y}$ (which is also
cos$\theta$) under different tensile stress. The magnetization is normalized
to the saturation magnetization. Clearly under no stress, the behavior is that
of a BSN where the magnetization fluctuates rail to rail and is mostly in the
state +1 or -1, and not in any intermediate state. As we increase the stress
(and depress the energy barrier), the behavior gradually transitions to that
of an ASN wherein the magnetization visits all states between -1 and +1 with
almost equal likelihood.

###  3.1 Energy cost of reconfiguration

We can make an order estimate of the energy cost of reconfiguration. This is
the quantity $(1/2)CV_{g}^{2}$, where $C$ is the capacitance of the two gate
pads in Fig. 3 and $V_{g}$ is the gate voltage needed to generate the required
stress. We see from Fig. 4 that 6.5 MPa of stress is enough to reconfigure a
BSN into an ASN. The gate voltage needed to generate a given stress $\sigma$
is given by $V_{g}\approx\sigma d/(Yd_{33})$, where $d_{33}$ is the diagonal
element of the piezoelectric tensor, $Y$ is the Young’s modulus of the
nanomagnet material and $d$ is the piezoelectric layer thickness. We will
assume that the piezoelectric is PMN-PT whose reported $d_{33}$ value is 2500
pC/N [17] and that $d$ = 300 nm. The Young’s modulus of Co is 209 GPa. This
will make the gate voltage needed to generate 6.5 MPa of stress = 3.6 mV. The
capacitance $C$ of the two gate pads is $C$ = $2\times\epsilon A/d$, where
$\epsilon$ is the dielectric constant of PMN-PT =
4000$\times$8.854$\times$10-12 F/m [18] and $A$ is the area of the gate pads =
100 nm $\times$ 100 nm. This makes $C$ = 2.4 fF. Hence the energy cost of
reconfiguration $(1/2)CV_{g}^{2}$ is only $\sim$ 8.5$\times$10-20 Joules,
which is miniscule. The gate voltage of 3.6 mV is above the noise voltage at
room temperature, which is $\sqrt{kT/C}$ = 1.3 mV ($kT$ = thermal energy). If
more noise resilience is desired, one can increase the gate voltage beyond 3.6
mV to obtain the desired noise margin, while still dissipating negligible
energy to reconfigure.

##  4 Application Space for Dynamically Reconfigurable Stochastic Neurons

Dynamic reconfigurability of the barrier height in a low barrier nanomagnet
through precise voltage (strain) control opens up some interesting
possibilities in neuromorphic hardware fabrics. We list a few potential
applications.

###  4.1 Precision and Adaptive Annealing Control in Energy-based Computation

One of the most important application of BSN is in solving binary optimization
problems, specifically finding a bit string that minimizes a complex Boolean
expression [19]. Often such problems are in computationally intractable NP-
Complete class and a Monte Carlo approach yields a computationally tractable
probabilistic solution method [20]. Stochastic neurons can be used to build
natural Monte Carlo hardware engines for solving such problems. The problem
expression can be cast as an adjacency matrix of a network spanning over the
bitstring basis set, and this adjacency matrix is then analogous to the
Hamiltonian of an equivalent Ising network, where the bitstring represents the
Ising spins [21]. In a purely software implementation, an annealing process is
simulated where the “temperature” of the system is slowly reduced and the
system settles in its ground state which encodes the solution bitstring.

In a hardware substrate, we directly build the Ising network and perform a
physical annealing. Typically, such a hardware approach will use low barrier
magnets (LBMs) and a current based scaling factor will control the “inverse”
temperature [22]. This current factor is global in nature and therefore
scaling this current up cools the whole network [23]. In contrast, the
approach presented here can seamlessly control the temperature via the control
over the barrier height through a gate voltage rather than through current
control. This is significantly easier to implement in a VLSI circuit.

Moreover, since this voltage control is directly over each individual neuron,
completely arbitrary annealing schedules can be implemented with only
polynomial increase in gate control circuit complexity and it brings precision
control over the annealing approach, i.e. a screwdriver instead of a hammer.

###  4.2 Control over Device-to-Device Variability and Data Retention Time

The gate control of barrier height via strain can be also used to control
variability in the barrier height in a large network. Deviation from the
designed shape of a LBM due to lack of precision in the fabrication method can
lead to significant variations in the natural barrier heights in LBM networks.
It has been shown that in a network with large barrier height variability,
significant error in computational results can be expected unless the compute
process is exponentially long on the magnitude of the barrier height
variability spread in the network, which in most practical applications is
unacceptable from the perspectives of both energy cost and throughput [24].
Therefore, gate control of barrier height (fig. 5 a) can enable both better
energy efficiency and throughput.

On a different note, magnetic random access memory (MRAM) has long been touted
as a solution to the memory hierarchy problem. However it’s use has thus far
been confined to replacement of solid-state non-volatile memory (NVM). The
retention time in a magnetic memory cell scales as $exp(\Delta/kT)$, where
$\Delta$ is the energy barrier height in the storage nanomagnet. The ability
to dynamically control the barrier heights in certain sections of a large
memory array may provide the ability to control memory retention times of
sections of the memory fabric (fig. 5 b). Since the vast majority of computing
architecture and programming models are optimized to a memory hierarchy of
speed/data retention rates, the ability to emulate the memory hierarchy in a
dynamically reconfigurable fashion can open up intriguing possibilities of
software defined hardware architectures, especially in the new age of
heterogeneous integration and chiplet based compute fabric design.

Figure 5: (a) Gate control over the barrier height (blue dashed line to red
solid line) of each individual memory cell to equalize memory retention and
hence reduce device-to-device variability in a memory array. (b) Gate control
over sections of memory fabric to emulate memory hierarchy in terms of
retention time-scales through barrier height modulation (color coded as light
to dark red boxes and corresponding barrier heights), e.g. $\sim\mu s$, $\sim
s$, $\sim$ years in a single integrated fabric

###  4.3 Control over Belief Uncertainty in Analog Stochastic Neurons

In analog-stochastic-neurons the noise is banded around the sigmoidal transfer
function and it’s magnitude is maximal at the mid point of the transfer curve
while minimal at the edges (fig. 6 a) [25], a feature which reflects the
belief uncertainty of the neuron in its state over the transfer curve:

$V_{out}(t)=\tanh(\beta*V_{in}(t))+\alpha(V_{in}(t))*V^{0}_{noise}(t)$ (4)

Where $\beta$, $\alpha$ are system transfer gain and noise profile function
respectively, $*$ is the convolution operator, $V^{0}_{noise}$ is the
normalized noise voltage (i.e. between $+V_{DD}/2$ and $-V_{DD}/2$), with
$V_{DD}$ being the power supply voltage. Through detailed simulations we have
observed that the empirical expression for $\alpha$ is given by the Gaussian
profile:

$\alpha(V_{in}(t))=\kappa\exp(-\frac{\nu V_{in}}{\sigma_{V_{n}}^{2}})$ (5)

where $\kappa$, $\nu$ are non-linear fitting functions dependent on the
barrier height, whereas $\sigma_{V_{n}}$ is the standard deviation of the
noise voltage profile, again a non-linear function of barrier height and the
transfer gain of the neuron cell.

Voltage control over the noise through barrier height modulation (fig. 6 b)
provides a clean way to control dynamically the magnitude of the belief and
its spread (fig. 6 c) . This can be used in applications such as reservoir
computing in continuous online-training mode where the network can reduce the
noise during the training phase, while incorporating it during the inference
stage [26]. This provides both robustness to the reservoir computing
operation, i.e. the training can be made accurate, and the inference can be
performed with significantly lower spread of the dynamic range of the
observation model weight spread.

Figure 6: (a) Illustrative circuit design of an ASN cell with a stochastic MTJ
(piezoelectric substrate omitted for clarity). (b) Gate control over the
barrier height modulates the stochasticity. (c) Three example input-output
characteristic curves of ASN transfer function (noisy signal is instantaneous
output, green smooth signal is expected time averaged output) over three
different barrier heights shows the degree of uncertainty on the belief of the
ASN over its transfer function.

##  5 Conclusion

We showed that strain mediated control of the energy barrier height in a low
barrier nanomagnet allows reconfiguring a BSN to an ASN and vice versa, thus
allowing a multitude of tasks to be performed in the same substrate. We also
discussed other applications of this modality. These few applications are only
a small subset of the potential functionalities that can be achieved in a gate
controlled stochastic neuron compute fabric. The dynamic and precision control
over the individual neurons through voltage control is well suited for
conventional VLSI design methodologies and fabrication practices. Further
exploration and development of this technology will provide a useful widget in
the toolkit of the rapidly expanding discipline of hardware neuromorphics.

## References

## References

  * [1] Hassan O, Faria R, Camsari K Y, Sun J Z and Datta S 2019 Low barrier magnet design for efficient hardware binary stochastic neurons. IEEE Magn. Lett., 10, 4502805.
  * [2] Hassan O, Datta S and Camsari K Y 2021 Quantitative evaluation of hardware binary stochastic neurons. Phys. Rev. Appl., 15, 064046.
  * [3] Khilwani D, Moghe V, Saraswat V, Kumbhare P, Baghini M J, Jandhyala S, Subramoney S and Ganguly U 2019 PrxCa1-xMnO3 based neuron for Boltzmann machine to solve “maximum cut” problem. APL Mater., 7 091112\.
  * [4] Borders W A, Pervaiz A Z, Fukami S, Camsari K Y, Ohno H and Datta S 2019 Nature, 573, 390.
  * [5] Salehi Fashami M, Roy K, Atulasimha J and Bandyopadhyay S 2011 Magnetization dynamics, Bennett clocking and associated energy dissipation in multiferroic logic. Nanotechnology, 22, 155201.
  * [6] Chikazumi S 1964 Physics of Magnetism (New York: Wiley).
  * [7] Cui J Z, Liang C Y, Paisley E A, Sepulveda A, Ihlefeld J F, Carman G P and Lynch C S 2015 Generation of localized strain in a thin film piezoelectric to control individual magnetoelectric heterostructures Appl. Phys. Lett. 107, 092903.
  * [8] Yang L, Zhao Y, Zhang S, Li P, Gao Y, Yang H, Hu H, Miao P, Liu Y, Chen A, Nan C W and Gao C 2014 Bipolar loop-like non-volatile strain in the (001)-oriented Pb(Mg1/3Nb2/3)O3-PbTiO3 single crystals Sci. Rep., 4, 4591.
  * [9] Yu L, Yu S W and Feng X Q 2007 Effects of electric fatigue on the butterfly curves of ferroelectric ceramics Mater. Sci. Eng. , 459, 273.
  * [10] Lupascu D C and Verdier C 2004 Fatigue anisotropy in lead zirconate titanate J. Eur. Ceram. Soc., 24, 1663.
  * [11] Wu T, et al. 2011 Domain-engineered switchable strain states in ferroelectric (011) $\left[Pb(Mg_{1/3}Nb_{2/3})O_{3}\right]_{1-x}-\left[PbTiO_{3}\right]_{x}$ (PMN-PT, x $\approx$ 0.32) single crystals J. Appl Phys., 109 124101\.
  * [12] Wu T, et al. 2011 Electrical control of reversible and permanent magnetization reorientation for magnetoelectric memory devices Appl. Phys. Lett., 98, 262504.
  * [13] Liu M., et al., 2013 Voltage impulse induced non-volatile ferroelastic switching of ferromagnetic resonance for reconfigurable magnetoelectric microwave devices Adv. Mater., 25, 4886.
  * [14] Bandyopadhyay S 2022 Comment on “Using dipole interaction to achieve nonvolatile voltage control of magnetism in multiferroic heterostructures”[Adv. Mater., 2021, 2105902] arXiv:2203.08236.
  * [15] Rahman R and Bandyopadhyay S 2023 Increasing flips per second and speed of p-computers by using dilute magnetic semiconductors IEEE Magnetics Letters, 14, 4500604.
  * [16] Rahman R and Bandyopadhyay S 2023 The strong sensitivity of the characteristics of binary stochastic neurons employing low barrier nanomagnets to small geometrical variations IEEE Transactions on Nanotechnology, 22, 112.
  * [17] Park S E and Shrout T R 1997 Ultrahigh strain and piezoelectric behavior ion relaxor based ferroelectric single crystals J. Appl. Phys., 82, 1804-1811.
  * [18] Pramanik R, Sahukar M K, Mohan Y, Praveenkumar B, Sangawar S R and Arockiarajan A 2019 Effect of grain size on piezoelectric, ferroelectric and dielectric properties of PMN-PT ceramics Ceramics International, 45, 5731.
  * [19] Glover Fred and Kochenberger Gary 2018 A tutorial on formulating and using QUBO models. arXiv preprint arXiv:1811.11538.
  * [20] Hochba Dorit S 1997 Approximation Algorithms for NP-Hard Problems. SIGACT News 28, 2 (June 1997), 40–52.
  * [21] Bian Z, Chudak F, Macready W G and Rose G 2010 The Ising model: teaching an old problem new tricks. D-wave systems, 2, 1-32.
  * [22] Camsari K Y, Faria R, Sutton B M and Datta S 2017 Stochastic p-bits for invertible logic. Physical Review.X, 7(3) 031014\.
  * [23] Sutton B, Camsari K Y, Behin-Aein B and Datta S 2017 Intrinsic optimization using stochastic nanomagnets. Sci. Rep., 7(1), 44370.
  * [24] Morshed M G, Ganguly S and Ghosh A W 2023. A Deep Dive into the Computational Fidelity of High Variability Low Energy Barrier Magnet Technology for Accelerating Optimization and Bayesian Problems. IEEE Magnetics Letters, 14, 6100405.
  * [25] Ganguly S, Camsari K Y and Ghosh A W 2021 Analog signal processing using stochastic magnets. IEEE Access, 9, 92640-92650.
  * [26] Morshed M G, Ganguly S and Ghosh A W 2023 Choose your tools carefully: a comparative evaluation of deterministic vs. stochastic and binary vs. analog neuron models for implementing emerging computing paradigms Frontiers in Nanotechnology, 5, 1146852.


* [30] V. V. Petrov Sums of independent random variables, Ergebnisse der Mathematik und ihrer Grenzgebiete 82 (1975) Springer, New York-Heidelberg, x+346 pp.
  * [31] V. V. Petrov, J. Robinson Large deviations for sums of independent non identically distributed random variables, Comm. Statist. Theory Methods 37 (2008) 2984–2990.
  * [32] Yu. V. Prokhorov, On a local limit theorem/or lattice distributions, Dokl. Akad. Nauk SSSR, 98 (1954) 535–538.
  * [33] Yu. A. Rozanov, On a local limit theorem for lattice distributions, Theory of Probability and Its Applications, 1957.
  * [34] L. Saulis, V. A. Statulevicius, Limit theorems for large deviations, Math. & Appl. (Soviet Series) 73 (1991) Kluwer, Dordrecht, viii+232 pp.
  * [35] S. Sethuraman and S.R.S Varadhan, A martingale proof of Dobrushin?s theorem for non-homogeneous Markov chains, Electron. J. Probab. 10 (2005) 1221–1235.
  * [36] Vo An Zung; Muhin, A. B.; To An Zung Certain local limit theorems for independent integer-valued random variables, (Russian) Izv. Akad. Nauk UzSSR Ser. Fiz.-Mat. Nauk 1980, no. 5, 9–15.