# Cubic singularities in binary linear electromechanical oscillators

Xin Zhou  College of Intelligence Science and Technology, NUDT, Changsha
410073, China  These authors contributed equally: Xin Zhou, Hui Jing, Xingjing
Ren.  Hui Jing  Key Laboratory of Low-Dimensional Quantum Structures and
Quantum Control of Ministry of Education, Department of Physics and Synergetic
Innovation Center for Quantum Effects and Applications, Hunan Normal
University, Changsha 410081, China  These authors contributed equally: Xin
Zhou, Hui Jing, Xingjing Ren.  Xingjing Ren  College of Intelligence Science
and Technology, NUDT, Changsha 410073, China  These authors contributed
equally: Xin Zhou, Hui Jing, Xingjing Ren.  Jianqi Zhang  State Key Laboratory
of Magnetic Resonance and Atomic and Molecular Physics, Wuhan Institute of
Physics and Mathematics, Innovation Academy of Precision Measurement Science
and Technology, Chinese Academy of Sciences, Wuhan 430071, China  Ran Huang
Theoretical Quantum Physics Laboratory, Cluster for Pioneering Research,
RIKEN, Wako-shi, Saitama 351-0198, Japan  Zhipeng Li  Department of Electrical
and Computer Engineering, National University of Singapore, Singapore 117576,
Singapore.  Xiaopeng Sun  College of Intelligence Science and Technology,
NUDT, Changsha 410073, China  Xuezhong Wu  College of Intelligence Science and
Technology, NUDT, Changsha 410073, China  Cheng-Wei Qiu  Department of
Electrical and Computer Engineering, National University of Singapore,
Singapore 117576, Singapore.  Franco Nori  Theoretical Quantum Physics
Laboratory, Cluster for Pioneering Research, RIKEN, Wako-shi, Saitama
351-0198, Japan  Dingbang Xiao  College of Intelligence Science and
Technology, NUDT, Changsha 410073, China

###### Abstract

Singularities arise in diverse disciplines and play a key role in both
exploring fundamental laws of physics and making highly-sensitive sensors [1,
2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 7, 8]. Higher-order ($\geq 3$)
singularities, with further improved performance [15, 16, 17, 18, 19],
however, usually require exquisite tuning of multiple ($\geq 3$) coupled
degrees of freedom [15, 16, 17] or nonlinear control [18, 19], thus severely
limiting their applications in practice. Here we propose theoretically and
confirm using mechanics experiments that, cubic singularities can be realized
in a coupled binary system without any nonlinearity, only by observing the
phase tomography of the driven response. By steering the cubic phase-
tomographic singularities in an electrostatically-tunable micromechanical
system, enhanced cubic-root response to frequency perturbation and voltage-
controlled nonreciprocity are demonstrated. Our work opens up a new phase-
tomographic method for interacted-system research, and sheds new light on
building and engineering advanced singular devices with simple and well-
controllable elements, with a wide range of applications including precision
metrology, portable nonreciprocal devices, and on-chip mechanical computing.

††preprint: NC/123-Singularities

Singularities, sometimes referred to as catastrophes, arise in diverse
disciplines and play an essential role in describing how the properties of an
object, that are dependent on certain controlling parameters, change
qualitatively even if the controlling parameters vary minimally [1]. The
unusual landscapes near the singularities are very useful for enhancing the
sensitivities of detection [2, 3, 4, 5, 6, 7, 8] as well as generating
nonreciprocity [9, 10, 11, 12, 13, 14, 7, 8]. Recently, higher-order
singularities have increasingly attracted attention, which have the potential
to provide higher performance and engender richer physics [15, 16, 17, 18,
19]. However, what prevents the higher-order singularities from being well
explored or exploited thus far is their difficulty in practical realization
and control. Usually, higher-order singularities call for multiple ($\geq 3$)
degrees of freedom [15, 16, 17] and are very difficult to construct and
adjust. Lately, it was theoretically predicted [18] and experimentally
observed [19] that introducing nonlinearity to binary non-Hermitian systems
may also realize higher-order singularities. These studies point out the
possibility of exploring higher-dimensional ($\geq 3$) physics using binary
systems. However, nonlinearity is still a stringent condition that will bring
intrinsic power consumption and reliability limits.

Micro- and nanomechanical resonators, with broad applications [20, 21, 22, 23,
24, 25], excellent in-situ controllability [26, 27, 28, 29], and rich
interactive phenomena [26, 27, 30, 31, 32, 28, 33, 10], provide an ideal
platform for exploring and exploiting singularities. Here we demonstrate, both
theoretically and experimentally, that cubic singularity arcs and nexus can be
realized in a binary linear micromechanical system, by observing the phase-
locked-loop (PLL) enabled tomographic dynamics of the coherent-coupling phase
responses. The cubic singularity arcs are featured by a series of stable
boundaries on a partially folded geometry made by the closed-loop oscillation
frequency of a pair of coherently coupled modes. The intersection of two
singularity arcs makes a singularity nexus. Behind the singularities lies
interesting polarization dynamics. In an electrically tunable micromechanical
system, we experimentally observe the cubic singularities. We confirm that the
singularity nexus can provide improved detecting sensitivity with a cubic-root
response, surpassing the binary singularities. Moreover, we demonstrate a
single-parameter-controlled nonreciprocity by traversing the projected
sensitivity arcs. Using the electrostatic tuning method, the nonreciprocity is
steered electrically.

Figure 1:  Cubic singularity in phase-tomographic dynamics: Concept. a Setup
for the phase-tomographic singularity. Mode 1 is driven and coherently coupled
to mode 2. A PLL is used to obtain a phase-tomographic closed-loop
oscillation. In this study, coherent coupling is produced by rotation
$\vec{\Omega}$, giving a coupling strength of $g=2\kappa\Omega$ with
$\kappa=0.85$. b Open-loop phase-frequency response ($\theta_{1}$, colored
surface) of the driven mode 1 as a function of coupling strength. PLL adjusts
the drive frequency to track the phase $\theta_{1}=-\pi/2$. The phase
tomography shows a “pitchfork” bifurcation. c The $\omega_{\text{d}}^{*}$ as a
function of degeneracy condition $\Delta\omega$ and coupling strength $g$.
Singularity arcs connected by a singularity nexus are formed by the stability
boundaries. Colors on the surface represent the relative phase of the coupled
modes. d Polarization dynamics of the balanced “pitchfork” bifurcation on the
classical Bloch sphere. e Stabilities and singularities on the Bloch sphere.
Figure 2:  Experimental realization of the phase-tomographic singularity. a
Micromechanical resonator with degenerate modes that are coherently coupled by
Coriolis effect. b Mode natural frequencies $\omega_{1,2}$ and their
difference $\Delta\omega$ vs tuning voltage $V_{\text{t}}$. c Spontaneous
chiral symmetry breaking and d the corresponding order parameter $\mathcal{N}$
illustrating a second-order phase transition. Error bars are the standard
deviation. e Measured $\omega_{\text{d}}^{*}$ if the tuning voltage
$V_{\text{t}}$ is swept adiabatically at constant angular velocities. The blue
dashed (red solid) curves indicate the upward (downward) sweeps. Singularities
and hysteresis are shown above $\Omega_{0}$. The gray surface is simulation. f
Cubic singularities projected to $V_{\text{t}}$-$\Omega$ plane. White-faced
points (black curves) are experimental (theoretical) data. The red and blue
curves are the corresponding $\mathcal{N}$ data in e revealing first-order
phase transitions.

## Results

### Concept:

We consider a pair of coherently coupled mechanics modes with tunable natural
frequencies $\omega_{1,2}$ and identical dissipation rate $\gamma$, as shown
in Fig. 1a. In this study, the coherent coupling is produced by the rotation-
induced Coriolis effect [34] (see Supplementary note 2), giving an angular
velocity $\Omega$-controlled coupling strength $g=2\kappa\Omega$, where
$\kappa\approx 0.85$ is the Coriolis coupling coefficient. Mode 1 is linearly
driven by an external sinusoidal force with frequency $\omega_{\text{d}}$
while mode 2 is set free. The linear displacement responses of the two modes
are labelled by $q_{1,2}=\cos(\omega_{\text{d}}t+\theta_{1,2})$, where
$|q_{1,2}|$ ($\theta_{1,2}$) are the amplitude (phase) responses.

If the two modes are at degeneracy
($\Delta\omega\equiv\omega_{2}-\omega_{1}=0$), the open-loop amplitude-
frequency response $|q_{1}|$ of mode 1 as a function of the coupling strength
$g$ displays normal-mode splitting [26] (see Supplementary Fig. 4 or 5). The
accompanied phase response $\theta_{1}$ of mode 1 is shown by the colored
surface in Fig. 1b. Here, we consider the tomography of the driven-mode phase
response with a constant oscillation phase, which is $-\pi/2$ for ideal
oscillators. The phase tomography is realized by imposing a PLL, which
produces a stable closed-loop oscillation (Fig. 1a). As shown by the red
contour in Fig. 1b, the controlled closed-loop frequency (denoted by
$\omega_{\text{d}}^{*}$, that fulfills the phase-tomographic condition
$\theta_{1}=-\pi/2$) as a function of coupling strength shows a “pitchfork”
bifurcation. It is noteworthy that this bifurcation is only related to the
landscape of the linear phase response $\theta_{1}$, and is different from its
counterparts in nonlinear dynamics [35]. The degenerate bifurcation point is
exactly the threshold between weak and strong coupling, $g=\gamma$. The middle
branch of the bifurcation (dotted curves) is unstable because the
corresponding amplitude response is in an antiresonance valley (see
Supplementary Fig. 4 or 5), which is unlikely to be detected by the PLL. The
other two branches (solid curves) are stable because the corresponding
amplitude responses are close to the resonant peaks.

The “pitchfork” bifurcation of $\omega_{\text{d}}^{*}$ changes as the
degeneracy condition $\Delta\omega$ varies. As a function of the coupling
strength $g$ and degeneracy condition $\Delta\omega$, $\omega_{\text{d}}^{*}$
makes a 3D surface (Fig. 1c), which is simulated by the cubic equation (see
Supplementary note 3)

$\displaystyle(\omega_{\text{d}}^{*}-\omega_{1})(\omega_{\text{d}}^{*}-\omega_{2}+\frac{\text{i}}{2}\gamma)(\omega_{\text{d}}^{*}-\omega_{2}-\frac{\text{i}}{2}\gamma)-\frac{1}{4}g^{2}(\omega_{\text{d}}^{*}-\omega_{2})=0.$
(1)

The inflectional part of the folded $\omega_{\text{d}}^{*}$ surface is
unstable. If the stability boundaries (magenta curves) are crossed by changing
the control parameters $g$ and $\Delta\omega$ in an adiabatic manner,
catastrophic jumps of the oscillation state takes place. The stability
boundaries constitute singularity arcs. The degenerate bifurcation point
(white star) is actually a singularity nexus that connects two singularity
arcs, giving a severely twisted $\omega_{\text{d}}^{*}$ geometry. The
singularity arcs and nexus, given the discriminant of the cubic Eq. 1,
construct cubic cusp catastrophes, the projection of which to the
$\Delta\omega$-$g$ parameter plane forms two cusp-connected parabola loci [1,
35] (see Supplementary note 4).

The phase-tomographic closed-loop oscillation of the system is described by a
state vector
$|\psi\rangle=\cos\frac{\phi}{2}|1\rangle+\text{e}^{i\vartheta}\sin\frac{\phi}{2}|2\rangle$,
where $\\{|1\rangle,|2\rangle\\}$ is the orthonormal basis of modes 1 and 2,
$\phi=2\arctan(\lvert q_{2}\rvert/\lvert q_{1}\rvert)$ is the polar angle, and
$\vartheta\equiv\theta_{2}-\theta_{1}$ is the relative phase between the
coupled modes, which is revealed by the colors of the surface in Fig. 1c. The
state vectors can be projected to a classical Bloch sphere with polar and
azimuthal angles given by $\phi$ and $\vartheta$, respectively (see
Supplementary note 5). The red trajectories on the Bloch sphere in Fig. 1d
illustrate the dynamics that underlies the degenerate “pitchfork” bifurcation
in Fig. 1b. The arrows indicate the $g$-increasing direction. There is a
corresponding polarization pattern in the $q_{1}$-$q_{2}$ plane for every
state vector on the Bloch sphere. At the initial state $|1\rangle$ without
coupling, $q_{2}=0$ makes a horizontal linear polarization. When
$\Omega\to\infty$, the states corresponding to the upper or lower stable
branches of $\omega_{\text{d}}^{*}$ approach the counter-clockwise circular
polarization state $|\text{CCW}\rangle=(|1\rangle-i|2\rangle)/\sqrt{2}$, or
the clockwise circular polarization state
$|\text{CW}\rangle=(|1\rangle+i|2\rangle)/\sqrt{2}$, respectively. While the
state corresponding to the middle unstable branch approaches $|2\rangle$ with
vertical linear polarization, because the antiresonance valley makes $q_{1}\to
0$. At the singularity nexus, $|q_{1}|=|q_{2}|$ and $\vartheta=0$ make a
$45^{\circ}$ linear polarization. The singularity arcs are marked by the
magenta curves on the Bloch sphere in Fig. 1e. The unstable, bistable, and
monostable regions on the Bloch sphere are filled with light red, light green,
and grey, respectively. The front (back) hemisphere corresponds to the
oscillations with positive (negative) angular velocities. The binary
polarization dynamics is a projection of a cubic singular dynamics.

Figure 3:  High sensitivity near the singularity nexus. a Phase-tomographic
frequency $\omega_{\text{d}}^{*}$ as a function of angular velocity $\Omega$
and degeneracy condition $\Delta\omega$. The contours of $\Omega=\Omega_{0}$
(blue curve) and $\Delta\omega=0$ (green curves) portray the sharp variation
of $\omega_{\text{d}}^{*}$ at the singularity nexus. b Phase-tomographic
frequency at the nexus angular velocity $\omega_{\text{d}}^{*}(\Omega_{0})$
and its shift from $\omega_{1}$,
$\delta\omega_{\text{X}}=\omega_{\text{d}}^{*}(\Omega_{0})-\omega_{1}$ as
functions of $\Delta\omega$. $\omega_{0}$ denotes $\omega_{1}$ at
$\Delta\omega=0$. In the range of $-0.25\gamma\leq\Delta\omega\leq
0.25\gamma$, $\delta\omega_{\text{X}}$ decreases monotonically to
$\Delta\omega$. c Frequency output $\delta\omega_{\text{X}}$ near the
singularity nexus versus the natural-frequency perturbation
$\epsilon=\Delta\omega$ from simulation (red solid curve) and experiment
(points). Eigenfrequency splits near an EP (blue dashed curve) and a DP (black
dot-dashed curve) are simulated as well. Error bars are the standard
deviation. d Logarithmic plot of the absolute data in c. The singularity nexus
has a cubic-root output, which provides higher sensitivity than the EP and DP.
Figure 4:  Voltage controlled nonreciprocity. A closed 1D parameter trajectory
that crosses both projected singularity loci in down-up-down a and up-down-up
b $V_{\text{t}}$ traversing process. The numbers indicate the traversing
order. c Adiabatic evolutions start at high (upper panel) or low (lower panel)
bistable branches following the down-up-down $V_{\text{t}}$ traversing
process. d Same as c for the up-down-up traversing process. Two processes that
start at an identical location and travel the same path in opposite directions
will reach distinct destinations. The nonreciprocity is ensured by crossing
both sides of the projected singularity loci.

### Realization:

To experimentally demonstrate the cubic singularities, we realize the scheme
in Fig. 1a using a capacitive microelectromechanical system [28] (see
Supplementary note 1). As shown in Fig. 2a, two micromechanical modes with
near-degenerate natural frequencies $\omega_{1,2}/2\pi\approx 3.85$ kHz and
equal dissipation rates $\gamma=2\pi\times 55.8$ mHz are coherently coupled by
the Coriolis effect. The degeneracy condition $\Delta\omega$ can be adjusted
electrostatically by a tuning voltage $V_{\text{t}}$, as shown in Fig. 2b (see
Methods). Mode 1 is driven externally and the linear displacement responses of
the two modes are detected using the Homodyne method (see Methods). A stable
$\theta_{1}=-\pi/2$ phase-tomographic closed-loop oscillation is maintained by
applying a PLL to mode 1.

First, we adjust $V_{\text{t}}$ to make $\Delta\omega\approx 0$ and
adiabatically sweep $\Omega$ from zero to $40^{\circ}/\text{s}$ to observe the
singularity nexus. The calculated and measured phase-tomographic closed-loop
frequency $\omega_{\text{d}}^{*}$ is shown in Fig. 2c. If the angular velocity
is below the strong-coupling threshold $\Omega<\gamma/(2\kappa)$,
$\omega_{\text{d}}^{*}$ is latched to $\omega_{1}$. At the threshold
(singularity nexus), $\Omega_{0}=\gamma/(2\kappa)=11.83^{\circ}/\text{s}$, the
$\omega_{\text{d}}^{*}$ transfers to one of the two stable bifurcation
branches randomly, revealing a spontaneous breaking of chiral symmetry. To
better describe the chirality, we re-expand the state vector in the
$\\{|\text{CCW}\rangle,|\text{CW}\rangle\\}$ basis,
$|\psi\rangle=c_{\text{ccw}}|\text{CCW}\rangle+c_{\text{cw}}|\text{CW}\rangle$,
and define the relative population of the $|\text{CCW}\rangle$ and
$|\text{CW}\rangle$ states as the order parameter
$\mathcal{N}=(|c_{\text{ccw}}|^{2}-|c_{\text{cw}}|^{2})/(|c_{\text{ccw}}|^{2}+|c_{\text{cw}}|^{2})=-\sin\phi\sin\vartheta$
(see Supplementary note 6). The order-parameter evolution in the symmetry-
breaking process of Fig. 2c is plotted in Fig. 2d, showing a second-order
phase transition. Below the singularity nexus, the oscillation is in the
linear-polarization phase. At the nexus, the oscillation transfers to one of
the degenerate $|\text{CCW}\rangle$ and $|\text{CW}\rangle$ dominant phases
randomly. The oscillation frequency of the $|\text{CCW}\rangle$
($|\text{CW}\rangle$) dominant phase increases (decreases) if $\Omega$ is
further increased from $\Omega_{0}$, because of the rotational Doppler effect
[36, 37]. The experimental data in Fig. 2c and d indicate that we
experimentally observed a linear-polarization-to-$|\text{CCW}\rangle$ phase
transition.

Next, we change $\Delta\omega$ by tailoring $V_{\text{t}}$ adiabatically while
keeping $\Omega$ unchanged at some specific values to observe the singularity
arcs. The measured $\omega_{\text{d}}^{*}$ is shown in Fig. 2e. The blue
dashed (red solid) curves illustrate the experimental results for upward
(downward) $V_{\text{t}}$ sweeps. If $\Omega>\Omega_{0}$, the sweeping curves
encounter discontinuous jumps at some specific values of $V_{\text{t}}$,
referred to as catastrophes or singularities [1]. The two upward and downward
curves at identical $\Omega$ form a hysteresis loop. The area of the
hysteresis loop decreases when reducing $\Omega$, which vanishes if
$\Omega\leq\Omega_{0}$. The experimentally detected singularities mapped to
the $V_{\text{t}}$-$\Omega$ parameter plane develop the predicted cusp-
connected parabola loci in Fig. 1c, as shown in Fig. 2f (see Supplementary
note 4). The order parameters $\mathcal{N}$ corresponding to the upward
(downward) sweeps in Fig. 2e are shown by the blue dashed (red solid) curves
in Fig. 2f, which indicate first-order transitions from the
$|\text{CCW}\rangle$ ($|\text{CW}\rangle$) dominant phase to the
$|\text{CW}\rangle$ ($|\text{CCW}\rangle$) dominant phase.

### Cubic-root sensitivity:

We now demonstrate the enhanced cubic-root sensitivity of the singularity
nexus. As shown in Fig. 3a, at the nexus ($\Omega=\Omega_{0}$ and
$\Delta\omega=0$), the closed-loop oscillation frequency
$\omega_{\text{d}}^{*}(\Omega_{0})$ is latched to the driven-mode natural
frequency $\omega_{1}$. Otherwise, if the degeneracy is broken,
$\Delta\omega\neq 0$, $\omega_{\text{d}}^{*}(\Omega_{0})$ will suddenly but
continuously deviate from $\omega_{1}$. The deviation
$\delta\omega_{\text{X}}=\omega_{\text{d}}^{*}(\Omega_{0})-\omega_{1}$ changes
sharply if $\Delta\omega$ shift from the nexus (Fig. 3b). Here, we consider
the perturbation $\epsilon$ that can affect the degeneracy condition,
$\epsilon\sim\Delta\omega$, and regard $\delta\omega_{\text{X}}$ in the
vicinity of the nexus as the sensing output of $\epsilon$, as shown by the red
curve in Fig. 3c (see Supplementary note 7). When plotted on a logarithmic
scale, it gives a cubic-root response near the nexus:
$\delta\omega_{\text{X}}\sim\epsilon^{1/3}$ (Fig. 3d), confirming the cubic
nature of the singularity nexus. To experimentally verify the cubic-root
behavior, we maintain an $\Omega_{0}$ rotation and introduce a fine-tuning
voltage $V_{\text{t}}$ to sweep across the nexus. By transforming
$V_{\text{t}}$ to $\epsilon$, the experimental input-output data are shown by
the red circles in Fig. 3c and d, which coincide well with the cubic-root
simulation.

We compare the sensitivities produced by the singularity nexus and a binary
exceptional-point (EP) singularity that is generated by a passive parity-time-
symmetric system whose damping difference is chosen to be equal to the
dissipation of our system (see Supplementary note 7). The blue dashed curves
in Fig. 3c and d demonstrate an $\epsilon^{1/2}$ dependency of the
eigenfrequency split $\delta\omega_{\text{EP}}$ near the EP. The sensitivity
produced by the singularity nexus is greater than that of the binary EP [2, 3,
4, 5] and is on par with that of the third-order EP [15]. When compared to the
standard output $\delta\omega_{\text{DP}}\sim\epsilon$ of the diabolic-point
(DP) system shown by the black dot-dashed curves in Fig. 3c and d, both
$\delta\omega_{\text{EP}}$ and $\delta\omega_{\text{X}}$ are much improved.

### Voltage-controlled nonreciprocity:

Lastly, we show that the phase-tomographic cubic singularity can produce a
voltage controlled nonreciprocity. We consider a closed 1D trajectory along
the tuning-voltage $V_{\text{t}}$ direction in the parameter plane, which
starts and ends at the bistable region and crosses both projected singularity
loci. If the trajectory is traversed in the down-up-down (up-down-up)
direction, as shown in Fig. 4a (b), the system ends at the low (high) branch
no matter which branch it is started, as shown in Fig. 4c (d). Even if it is
started at an identical location, opposite traversal directions lead to
different ending branches, as shown in Fig. 4c or d. The ending location only
depends on the traversal direction, not the starting place.

In the experiments shown in Fig. 4, the tuning voltage $V_{\text{t}}$ is the
only steering knob of the nonreciprocal state transfer, while the angular
velocity is set to be a constant value, $\Omega=60^{\circ}/\text{s}$. In fact,
the $\Omega$ value of the $V_{\text{t}}$-controlled nonreciprocal process can
be changed almost arbitrarily, as long as the start/end points are located at
the bistable region and both projected singularity loci are crossed. This
single-parameter-controlled nonreciprocity is more desirable than that of the
binary EP singularity, which is guaranteed by two-parameter encircling [38, 9,
10, 11, 12, 13, 14]. Benefiting from the electrostatic tunability of our
device, the phase-tomographic cubic singularity provides a voltage-controlled
nonreciprocity.

## Discussion

Although we experimentally demonstrate the phase-tomographic closed-loop cubic
singularity based on the Coriolis coupling, in principle, it can also be
realized using ordinary linear coherent coupling (see Supplementary
discussion). This study may open up a new tomographic dimension for phase-
related interactive dynamics study, which is applicable for a wide range of
disciplines such as optics, optomechanics, or hybrid quantum systems. Our
discovery makes it possible to construct advanced singularities using highly
controllable elements. It also enhances the understanding of the closed-loop
oscillation dynamics, and extends coherent control into the singularity
region. Potential applications of the closed-loop singularity include precise
sensing, deep-sub-linewidth mode matching, rapid mode switching, and
generating portable nonreciprocity. Moreover, the hysteresis in the closed-
loop oscillation is promising for mechanical computing [39]. The bit
abstraction of the closed-loop oscillations is independent of vibration
amplitudes, which may provide potential advantages in power consumption and
lifetime. Future studies can also investigate phase-tomographic singularities
originating from different kinds of coupling, the interplay with other kinds
of singularities, and the phase-tomographic dynamics in many-body systems with
more degrees of freedom [40, 17, 41].

## Methods

### Electrostatic frequency tuning:

The tuning voltage $V_{\text{t}}$ introduces electrostatic negative stiffness
to both modes 1 and 2,
$\omega^{2}_{1,2}(V_{\text{t}})={\omega^{\prime}}_{1,2}^{2}(0)-T_{1,2}(V_{0}-V_{\text{t}})^{2}$,
where $\omega^{\prime}_{1,2}(0)$ denotes the natural frequencies of the bare
mechanical modes. The electrostatic tuning factors $T_{1,2}$ are proportional
to the capacitive area, the inverse of the modal mass, and the inverse of the
cubic of the capacitive gap. The stiffness perturbation induced by $V_{0}$
exists even if the tuning voltage $V_{\text{t}}$ is absent, so it can be
included in the intrinsic natural frequencies. By defining
$\omega_{1,2}^{2}(0)={\omega^{\prime}}^{2}_{1,2}(0)-T_{1,2}V_{0}^{2}$, and
assuming that the electrostatic stiffness perturbation is small relative to
the intrinsic stiffness, we have

$\displaystyle\omega_{1,2}(V_{\text{t}})$
$\displaystyle=\sqrt{\omega_{1,2}^{2}(0)+T_{1,2}(2V_{0}V_{\text{t}}-V_{\text{t}}^{2})}$
$\displaystyle\approx\omega_{1,2}(0)+K_{1,2}(2V_{0}V_{\text{t}}-V_{\text{t}}^{2}),$
(2)

where the tuning coefficients are defined by
$K_{1,2}=T_{1,2}/[2\omega_{1,2}(0)]$.

The experimentally measured natural frequencies $\omega_{1,2}$ as functions of
$V_{\text{t}}$ are shown by the red and blue points in the inset of Fig. 2b,
which are fitted (curves) to the model (2) with parameters
$\omega_{1}(0)=2\pi\times 3,852.92$ Hz, $\omega_{2}(0)=2\pi\times 3,856.43$
Hz, and $V_{0}=2.5~{}\text{V}$. The tuning coefficients are fitted to be
$K_{1}=1.29\times 10^{-2}~{}\text{rad}~{}\text{s}^{-1}~{}\text{V}^{-2}$ and
$K_{2}=6.40\times 10^{-2}~{}\text{rad}~{}\text{s}^{-1}~{}\text{V}^{-2}$. The
relationship between the difference of natural frequencies (degeneracy
condition) $\Delta\omega=\omega_{2}-\omega_{1}$ and the tuning voltage
$V_{\text{t}}$ is further given by

$\displaystyle\Delta\omega(V_{\text{t}})$
$\displaystyle=\Delta\omega(0)+(K_{2}-K_{1})(2V_{0}V_{\text{t}}-V_{\text{t}}^{2}),$
(3)

where $\Delta\omega(0)=\omega_{2}(0)-\omega_{1}(0)$. The experimentally
measured $\Delta\omega$ data (green points in Fig. 2b) coincide well with the
tuning model (3) (green curve in Fig. 2b).

### Homodyne measurement:

The antinodal displacements of the two micromechanical modes $q_{1,2}=\lvert
q_{1,2}\rvert\cos(\omega_{\text{d}}t+\theta_{1,2})$ are picked up by
capacitive transducers, transformed into voltage signals by two charge
amplifiers integrated into a printed circuit board, and then recorded by a
two-channel lock-in amplifier (Zurich Instruments HF2LI). The amplitudes
$\lvert q_{1,2}\rvert$ and phases $\theta_{1,2}$ relative to the driving
signal are obtained by dual-phase demodulation techniques. In this process,
$q_{j}(\omega_{\text{d}},t)$ is split and separately mixed with the driving
reference signal $\cos{\omega_{\text{d}}t}$ and a $\pi/2$-phase-shifted copy
of it,

$\displaystyle\lvert
q_{j}\rvert\cos(\omega_{\text{d}}t+\theta_{j})\times\cos{\omega_{\text{d}}t}$
$\displaystyle=$ $\displaystyle\frac{\lvert
q_{j}\rvert}{2}\left[\cos(2\omega_{\text{d}}t+\theta_{j})+\cos(\theta_{j})\right],$
$\displaystyle\lvert
q_{j}\rvert\cos(\omega_{\text{d}}t+\theta_{j})\times\cos{(\omega_{\text{d}}t+\frac{\pi}{2})}$
$\displaystyle=$ $\displaystyle\frac{\lvert
q_{j}\rvert}{2}\left[-\sin(2\omega_{\text{d}}t+\theta_{j})+\sin(\theta_{j})\right].$

The high-harmonic components of the mixed signals are removed using low-pass
filters and the remaining in-phase component $X_{j}=\frac{\lvert
q_{j}\rvert}{2}\cos(\theta_{j})$ and quadrature component $Y_{j}=\frac{\lvert
q_{j}\rvert}{2}\sin(\theta_{j})$ are obtained. By transforming into the polar
coordinates, we can obtain the amplitude and phase,

$\displaystyle\lvert q_{j}\rvert$ $\displaystyle=\sqrt{X^{2}+Y^{2}},$
$\displaystyle\theta_{j}$ $\displaystyle=\arctan\frac{Y}{X}.$

## Data availability

Data relevant to the figures and conclusions of this manuscript are available
at https://doi.org/10.6084/m9.figshare.19609350.

## References

  * Arnold [1984] V. I. Arnold,  _Catastrophe Theory_ (Springer-Verlag, 1984).
  * Liu  _et al._ [2016] Z.-P. Liu, J. Zhang, S. K. Ozdemir, B. Peng, H. Jing, X.-Y. Lü, C.-W. Li, L. Yang, F. Nori, and Y.-x. Liu, Metrology with $\mathcal{PT}$-symmetric cavities: Enhanced sensitivity near the $\mathcal{PT}$-phase transition, Physical Review Letters 117, 110802 (2016).
  * Chen  _et al._ [2017] W. Chen, S. K. Ozdemir, G. Zhao, J. Wiersig, and L. Yang, Exceptional points enhance sensing in an optical microcavity, Nature 548, 192 (2017).
  * Lai  _et al._ [2019] Y.-H. Lai, Y.-K. Lu, M.-G. Suh, Z. Yuan, and K. Vahala, Observation of the exceptional-point-enhanced Sagnac effect, Nature 576, 65 (2019).
  * Hokmabadi  _et al._ [2019] M. P. Hokmabadi, A. Schumer, D. N. Christodoulides, and M. Khajavikhan, Non-Hermitian ring laser gyroscopes with enhanced Sagnac sensitivity, Nature 576, 70 (2019).
  * Kononchuk  _et al._ [2022] R. Kononchuk, J. Cai, F. Ellis, R. Thevamaran, and T. Kottos, Exceptional-point-based accelerometers with enhanced signal-to-noise ratio, Nature 607, 697 (2022).
  * Miri and Alù [2019] M.-A. Miri and A. Alù, Exceptional points in optics and photonics, Science 363, eaar7709 (2019).
  * Ozdemir  _et al._ [2019] S. K. Ozdemir, S. Rotter, F. Nori, and L. Yang, Parity-time symmetry and exceptional points in photonics, Nature Materials 18, 783 (2019).
  * Doppler  _et al._ [2016] J. Doppler, A. A. Mailybaev, J. Böhm, U. Kuhl, A. Girschik, F. Libisch, T. J. Milburn, P. Rabl, N. Moiseyev, and S. Rotter, Dynamically encircling an exceptional point for asymmetric mode switching, Nature 537, 76 (2016).
  * Xu  _et al._ [2016] H. Xu, D. Mason, L. Jiang, and J. G. E. Harris, Topological energy transfer in an optomechanical system with exceptional points, Nature 537, 80 (2016).
  * Hassan  _et al._ [2017] A. U. Hassan, B. Zhen, M. Soljačić, M. Khajavikhan, and D. N. Christodoulides, Dynamically encircling exceptional points: Exact evolution and polarization state conversion, Physical Review Letters 118, 093002 (2017).
  * Yoon  _et al._ [2018] J. W. Yoon, Y. Choi, C. Hahn, G. Kim, S. H. Song, K.-Y. Yang, J. Y. Lee, Y. Kim, C. S. Lee, J. K. Shin, H.-S. Lee, and P. Berini, Time-asymmetric loop around an exceptional point over the full optical communications band, Nature 562, 86 (2018).
  * Zhang  _et al._ [2018] X.-L. Zhang, S. Wang, B. Hou, and C. T. Chan, Dynamically encircling exceptional points: In situ control of encircling loops and the role of the starting point, Physical Review X 8, 021066 (2018).
  * Nasari  _et al._ [2022] H. Nasari, G. Lopez-Galmiche, H. E. Lopez-Aviles, A. Schumer, A. U. Hassan, Q. Zhong, S. Rotter, P. LiKamWa, D. N. Christodoulides, and M. Khajavikhan, Observation of chiral state transfer without encircling an exceptional point, Nature 605, 256 (2022).
  * Hodaei  _et al._ [2017] H. Hodaei, A. U. Hassan, S. Wittek, H. Garcia-Gracia, R. El-Ganainy, D. N. Christodoulides, and M. Khajavikhan, Enhanced sensitivity at higher-order exceptional points, Nature 548, 187 (2017).
  * Tang  _et al._ [2020] W. Tang, X. Jiang, K. Ding, Y.-X. Xiao, Z.-Q. Zhang, C. T. Chan, and G. Ma, Exceptional nexus with a hybrid topological invariant, Science 370, 1077 (2020).
  * del Pino  _et al._ [2022] J. del Pino, J. J. Slim, and E. Verhagen, Non-Hermitian chiral phononics through optomechanically induced squeezing, Nature 606, 82 (2022).
  * Wang  _et al._ [2019] H. Wang, S. Assawaworrarit, and S. Fan, Dynamics for encircling an exceptional point in a nonlinear non-Hermitian system, Optics Letters 44, 638 (2019).
  * Bai  _et al._ [2022] K. Bai, L. Fang, T.-R. Liu, J.-Z. Li, D. Wan, and M. Xiao, Nonlinearity enabled higher-order exceptional singularities with ultra-enhanced signal-to-noise ratio, National Science Review 11, nwac259 (2022).
  * Ng  _et al._ [2015] E. Ng, Y. Yang, V. A. Hong, C. H. Ahn, D. B. Heinz, I. Flader, Y. Chen, C. L. M. Everhart, B. Kim, R. Melamud, R. N. Candler, M. A. Hopcroft, J. C. Salvia, S. Yoneoka, A. B. Graham, M. Agarwal, M. W. Messana, K. L. Chen, H. K. Lee, S. Wang, G. Bahl, V. Qu, C. F. Chiang, T. W. Kenny, A. Partridge, M. Lutz, G. Yama, and G. J. O’Brien, The long path from MEMS resonators to timing products, in  _IEEE International Conference on Micro Electro Mechanical Systems_ (2015) pp. 1–2.
  * Nguyen  _et al._ [1998] C.-C. Nguyen, L. Katehi, and G. Rebeiz, Micromachined devices for wireless communications, Proceedings of the IEEE 86, 1756 (1998).
  * Hanay  _et al._ [2012] M. S. Hanay, S. Kelber, A. K. Naik, D. Chi, S. Hentz, E. C. Bullard, E. Colinet, L. Duraffourg, and M. L. Roukes, Single-protein nanomechanical mass spectrometry in real time, Nature Nanotechnology 7, 602 (2012).
  * Cleland and Roukes [1998] A. N. Cleland and M. L. Roukes, A nanometre-scale mechanical electrometer, Nature 392, 160 (1998).
  * Middlemiss  _et al._ [2016] R. P. Middlemiss, A. Samarelli, D. J. Paul, J. Hough, S. Rowan, and G. D. Hammond, Measurement of the Earth tides with a MEMS gravimeter, Nature 531, 614 (2016).
  * Ayazi and Najafi [2001] F. Ayazi and K. Najafi, A HARPSS polysilicon vibrating ring gyroscope, Journal of Microelectromechanical Systems 10, 169 (2001).
  * Faust  _et al._ [2012] T. Faust, J. Rieger, M. J. Seitner, P. Krenn, J. P. Kotthaus, and E. M. Weig, Nonadiabatic dynamics of two strongly coupled nanomechanical resonator modes, Physical Review Letters 109, 037205 (2012).
  * Mahboob  _et al._ [2012] I. Mahboob, K. Nishiguchi, H. Okamoto, and H. Yamaguchi, Phonon-cavity electromechanics, Nature Physics 8, 387 (2012).
  * Zhou  _et al._ [2019] X. Zhou, C. Zhao, D. Xiao, J. Sun, G. Sobreviela, D. D. Gerrard, Y. Chen, I. Flader, T. W. Kenny, X. Wu, and A. A. Seshia, Dynamic modulation of modal coupling in microelectromechanical gyroscopic ring resonators, Nature Communications 10, 4980 (2019).
  * Eriksson  _et al._ [2023] A. M. Eriksson, O. Shoshani, D. López, S. W. Shaw, and D. A. Czaplewski, Controllable branching of robust response patterns in nonlinear mechanical resonators, Nature Communications 14, 161 (2023).
  * Okamoto  _et al._ [2013] H. Okamoto, A. Gourgout, C.-Y. Chang, K. Onomitsu, I. Mahboob, E. Y. Chang, and H. Yamaguchi, Coherent phonon manipulation in coupled mechanical resonators, Nature Physics 9, 598 (2013).
  * Faust  _et al._ [2013] T. Faust, J. Rieger, M. J. Seitner, J. P. Kotthaus, and E. M. Weig, Coherent control of a classical nanomechanical two-level system, Nature Physics 9, 485 (2013).
  * Sun  _et al._ [2016] F. Sun, X. Dong, J. Zou, M. I. Dykman, and H. B. Chan, Correlated anomalous phase diffusion of coupled phononic modes in a sideband-driven resonator, Nature Communications 7, 12694 (2016).
  * Miao  _et al._ [2022] T. Miao, X. Zhou, X. Wu, Q. Li, Z. Hou, X. Hu, Z. Wang, and D. Xiao, Nonlinearity-mediated digitization and amplification in electromechanical phonon-cavity systems, Nature Communications 13, 2352 (2022).
  * Li  _et al._ [2018] K. Li, H. Fu, and Y. Li, Coriolis-force-induced coupling between two modes of a mechanical resonator for detection of angular velocity, Physical Review A 98, 023862 (2018).
  * Strogatz [2015] S. H. Strogatz,  _Nonlinear dynamics and chaos: with applications to physics, biology, chemistry, and engineering_ , 2nd ed. (CRC Press, 2015).
  * Fleury  _et al._ [2014] R. Fleury, D. L. Sounas, C. F. Sieck, M. R. Haberman, and A. Alù, Sound isolation and giant linear nonreciprocity in a compact acoustic circulator, Science 343, 516 (2014).
  * Pan  _et al._ [2021] D. Pan, H. Xu, and F. J. G. de Abajo, Rotational Doppler cooling and heating, Science Advances 7, eabd6705 (2021).
  * Gao  _et al._ [2015] T. Gao, E. Estrecho, K. Y. Bliokh, T. C. H. Liew, M. D. Fraser, S. Brodbeck, M. Kamp, C. Schneider, S. Höfling, Y. Yamamoto, F. Nori, Y. S. Kivshar, A. G. Truscott, R. G. Dall, and E. A. Ostrovskaya, Observation of non-Hermitian degeneracies in a chaotic exciton-polariton billiard, Nature 526, 554 (2015).
  * Yasuda  _et al._ [2021] H. Yasuda, P. R. Buskohl, A. Gillman, T. D. Murphey, S. Stepney, R. A. Vaia, and J. R. Raney, Mechanical computing, Nature 598, 39 (2021).
  * Matheny  _et al._ [2019] M. H. Matheny, J. Emenheiser, W. Fon, A. Chapman, A. Salova, M. Rohden, J. Li, M. Hudoba de Badyn, M. Pósfai, L. Duenas-Osorio, M. Mesbahi, J. P. Crutchfield, M. C. Cross, R. M. D’Souza, and M. L. Roukes, Exotic states in a simple network of nanoelectromechanical oscillators, Science 363, eaav7932 (2019).
  * Brown  _et al._ [2022] C. D. Brown, S.-W. Chang, M. N. Schwarz, T.-H. Leung, V. Kozii, A. Avdoshkin, J. E. Moore, and D. Stamper-Kurn, Direct geometric probe of singularities in band structure, Science 377, 1319 (2022).

###### Acknowledgements.

X.Z. thank Prof. Ashwin Seshia from University of Cambridge and Prof. Zenghui
Wang from University of Electronic Science and Technology of China for helpful
discussions. This work is partly supported by the National Natural Science
Foundation of China (NSFC) (U21A20505 and 51905539), the Young Elite Scientist
Sponsorship Program by CAST (YESS20200127), and the Natural Science Foundation
of Hunan Province for Excellent Young Scientists (2021JJ20049). This work is
primarily supported the National Key R&D Program of China (NKPs)
(2022YFB3204901).

## Author contributions

X.Z. conceived the idea and designed the research. X.Z. and X.R. performed the
experiments with assistance of X.S. X.Z. designed the device. X.Z. and D.X.
fabricated the device. X.Z., X.R., X.S., D.X., and X.W. developed the test
equipments. X.Z., H.J., and J.Z. conducted the theory with advice from F.N.
and C.-W.Q. X.Z. and H.J. wrote the manuscript with inputs from all authors.
The project was jointly supervised by X.Z. and H.J.

## Competing interests

The authors declare no competing interests.


11institutetext: Drexel University, Philadelphia PA 19104, USA
11email<EMAIL_ADDRESS>
# Thermal-Aware Compilation of Spiking Neural Networks to Neuromorphic
Hardware

Twisha Titirsha  00000-0002-2142-2283 Anup Das  0000-0002-5673-2636

###### Abstract

Hardware implementation of neuromorphic computing can significantly improve
performance and energy efficiency of machine learning tasks implemented with
spiking neural networks (SNNs), making these hardware platforms particularly
suitable for embedded systems and other energy-constrained environments.  We
observe that the long bitlines and wordlines in a crossbar of the hardware
create significant current variations when propagating spikes through its
synaptic elements, which are typically designed with non-volatile memory
(NVM). Such current variations create a thermal gradient within each crossbar
of the hardware, depending on the machine learning workload and the mapping of
neurons and synapses of the workload to these crossbars. This thermal gradient
becomes significant at scaled technology nodes and it increases the leakage
power in the hardware leading to an increase in the energy consumption. We
propose a novel technique to map neurons and synapses of SNN-based machine
learning workloads to neuromorphic hardware. We make two novel contributions.
First, we formulate a detailed thermal model for a crossbar in a neuromorphic
hardware incorporating workload dependency, where the temperature of each NVM-
based synaptic cell is computed considering the thermal contributions from its
neighboring cells. Second, we incorporate this thermal model in the mapping of
neurons and synapses of SNN-based workloads using a hill-climbing heuristic.
The objective is to reduce the thermal gradient in crossbars. We evaluate our
neuron and synapse mapping technique using 10 machine learning workloads for a
state-of-the-art neuromorphic hardware. We demonstrate an average 11.4K
reduction in the average temperature of each crossbar in the hardware, leading
to a 52% reduction in the leakage power consumption (11% lower total energy
consumption) compared to a performance-oriented SNN mapping technique.

###### Keywords:

Neuromorphic computing Spiking Neural Network Non-Volatile Memory (NVM) Phase-
Change Memory (PCM) Temperature Leakage power consumption Crossbar.

##  1 Introduction

Spiking Neural Networks (SNNs) are machine learning models designed with
spike-based computations and bio-inspired learning algorithms [37]. Neurons
communicate information using spikes via synapses. SNNs are used to implement
both supervised and unsupervised machine learning approaches. Our focus is on
supervised approaches, where a machine learning model is first trained using
training data, and then deployed for inference with in-field data.

Neuromorphic hardware such as TrueNorth [24], Loihi [23], and DYNAP-SE [39]
can significantly improve the energy efficiency of SNNs, thanks to their
event-driven computations, efficient implementations of biological neurons
using CMOS and FinFET technologies, and the use of Non-Volatile Memory (NVM)
such as Phase-Change Memory (PCM) [10, 50, 46, 48, 49], Oxide-base Resistive
RAM (OxRRAM) [38], and Spin-Transfer Torque Magnetic or Spin-Orbit-Torque RAM
(STT- and SoT-MRAM) [41] for high density synaptic storage. Therefore,
neuromorphic hardware can be used to implement machine learning tasks on
power-constrained environments such as embedded systems, and sensor and edge
devices of the Internet-of-Things (IoT) [27].

A neuromorphic hardware is implemented as a tile-based architecture [11] with
a shared interconnect in the form of Networks-on-Chip (NoC) or Segmented Bus
[9] (see Figure 1a). A tile in a neuromorphic hardware is designed as a
crossbar, which is an organization of top electrodes (wordlines) and bottom
electrodes (bitlines), with NVM-based synaptic elements at their intersections
(Figure 1b). A synaptic element is connected to a bitline and a wordline using
an access transistor (Figure 1c). Within a crossbar, the pre-synaptic neurons
are mapped on the wordlines, while the post-synpatic neurons are mapped along
the bitlines. An $n\times n$ crossbar has $n$ pre-synaptic neurons, $n$ post-
synaptic neurons, and $n^{2}$ NVM cells. A pre-synaptic neuron’s spike voltage
from a wordline is multiplied with the conductance of the NVM to generate a
current. Currents from multiple wordlines are integrated on a bitline,
implementing forward propagation of neuron excitation. This is illustrated in
Figure 1b.

Figure 1: (a) Tile-based neuromorphic hardware. (b) A crossbar of a
neuromorphic tile. (c) An NVM-based synaptic cell consisting of an access
transistor and an NVM.

We investigate the internal architecture of a crossbar and observe that the
bitlines and wordlines of a crossbar consist of parasitic elements, which
consist of capacitance and resistance of the metal interconnect as shown in
Figure 2. These parasitic elements create variation in current propagating
along different paths in the crossbar. The figure illustrates the shortest and
the longest current paths in a crossbar, where the length of a path is
measured in terms of the number of parasitic components that are present on
the path. Current differences create variation in access speed of the
different synaptic elements in the crossbar [26, 55, 51]. A conservative
design practice is to use a common spike voltage to obtain the required access
speed of the synaptic element on the longest current path.

Figure 2: Parasitc components on the bitlines and wordlines in a crossbar.

We argue that this conservative approach creates current differences in a
crossbar, leading to a wide thermal gradient. Figure 3 illustrates the current
and thermal variations in a 128x128 PCM crossbar at 65nm technology node.
Accessing the synaptic cells on shorter current paths (bottom left corner of
Figure 3b) generate higher temperatures than those on longer current paths
(top right corner). Due to the exponential dependency of leakage current on
temperature [36], the leakage current through cells with higher temperature is
much higher than the current through cells with lower temperature. So,
frequently accessing the cells on shorter current paths when executing a
workload can lead to higher leakage power consumption in the crossbar.

(a) Current variation for PCM access operations in a 128x128 crossbar. (b)
Temperature gradient in a 128x128 crossbar.

Figure 3: Current variation and temperature gradient in a 128x128 crossbar at
$65nm$ process node with $T_{amb}=298K$. The PCM crystallization point is
360K.

Existing techniques to map neurons and synapses of SNNs to neuromorphic
hardware have mostly focused on improving performance and circuit aging [22,
17, 5, 6, 1, 44, 47, 45, 7, 8]. These techniques do not consider the thermal
gradient in a crossbar and therefore, they can increase the leakage power
significantly. We build the case for one such mapping technique – SpiNeMap
[6]. The leakage energy using this technique constitute between 20% to 30% of
the total energy consumption for the typical machine learning workloads (see
Section 4), where the total energy of a neuromorphic hardware includes the
energy to generate spikes, the energy to communicate spikes, and the leakage
energy. Therefore, reducing the leakage power (which we demonstrate in this
work) will lead to a significant reduction of the total energy consumption.

Our goal is to minimize the leakage power consumption. We achieve this goal by
lowering the average temperature of each crossbar using the proposed mapping
technique. To this end, we make the following two key contributions.

  * •

Contribution 1: We propose a new comprehensive thermal model of a crossbar
designed with phase-change memory (PCM). Our model incorporates 1) workload
dependency, i.e., the temperature obtained in processing spike trains from a
given SNN-based machine learning workload, and 2) spatial thermal
dependencies, i.e., the temperature contributions from the neighboring cells
based on their synaptic excitation in the workload.

  * •

Contribution 2: We propose a novel neuron and synapse mapping approach
incorporating the thermal model using a hill climbing heuristic. The objective
of the heuristic is to allocate the neurons and synapses of an SNN to the
crossbars of the hardware such that the maximum average temperature of all
crossbars is minimized, which lowers its leakage power consumption.

We evaluate the proposed technique with 10 machine learning applications from
three most commonly-used neural network topology – convolution neural network
(CNN), multilayer perceptron (MLP), and recurrent neural network (RNN).
Evaluation for DYNAP-SE [39], a state-of-the-art neuromorphic hardware
demonstrates the reduction of temperature, leading to a significant reduction
in the leakage current.

##  2 Workload-dependant Thermal Model of Crossbars

In this section, we develop a workload-dependent thermal model of crossbars in
a neuromorphic hardware, considering PCM-based synaptic elements. We start by
reviewing the internals of a PCM device. The proposed thermal model can be
generalized to other NVMs such as OxRRAM and SOT-/STT-MRAM exploiting their
specific structures.

Figure 4(a) illustrates how a chalcogenide semiconductor alloy is used to
build a PCM cell. The amorphous phase (logic ‘0’) in this alloy has higher
resistance than the crystalline phase (logic ‘1’). Ge2Sb2Te5 (GST) is the most
commonly used alloy for PCM [54] due to its high amorphous-to-crystalline
resistance ratio, fast switching between phases, and high endurance. However,
other chalcogenide alloys are also explored due to their better data retention
properties [40]. Phase changes in a PCM cell are induced by injecting current
into the resistor-chalcogenide junction and heating the chalcogenide alloy.

Figure 4 (b) shows the different current profiles needed to program and read
in a PCM device. To RESET a PCM cell, a high power pulse of short duration is
applied and quickly terminated. This first raises the temperature of the
chalcogenide alloy to 650$\degree$C, above its melting point. The melted alloy
subsequently cools extremely quickly, locking into an amorphous phase. To SET
a PCM cell, the chalcogenide alloy is heated above its crystallization
temperature, but below its melting point for a sufficient amount of time.
Finally, to read the content (i.e., know the phase) of a PCM cell, a small
electrical pulse is applied that is sufficiently low so as not to induce phase
change in the PCM cell. We focus on PCM read for the inference of supervised
machine learning approaches

Figure 4: (a) A PCM cell and (b) Current needed to operate a PCM cell.

Many prior works have developed thermal models for PCM devices [52, 12, 32].
However, these models are developed for individual PCM cell considering the
effect of crystallization and amorphization (synaptic weight updates in the
context of machine learning). In other words, these models have the following
two key limitations for their use in the context of neuromorphic computing.
First, they do not consider spatial dependencies, i.e., the thermal
contributions from neighboring PCM cells considering their utilization in a
machine learning workload. Second, the thermal impact due to PCM reads (as
required for machine learning inference) is not modeled. Figure 5a shows the
thermal interactions in a crossbar. When a cell is accessed repeatedly within
a short time window, there remains very little scope for heat inside the cell
to be dissipated. As a result temperature keeps rising on every access,
building on the undissipated components, and dissipating the heat to its
neighboring cells, raising their temperature.

(a) Thermal interactions in a PCM-based crossbar. (b) Iterative thermal
computations.

Figure 5: Building thermal model of a PCM-based crossbar.

Figure 5b shows the proposed iterative approach of computing the temperature
of a crossbar. The model computes the temperature of a single PCM cell
incorporating 1) thermal contributions from its neighbors and 2) its
activation within a workload.

The temperature of a single PCM cell is computed using Joule heating, $W_{j}$
and heat dissipation, $W_{d}$ which is given by the following equation [35],

$T_{PCM}=\int\frac{W_{j}-W_{d}}{C\times V}dt$ (1)

where C and V are heat capacity of GST and volume of the active region of the
cell respectively. The heat generation in the PCM cell is given by,

$W_{j}=I_{PCM}^{2}\times R_{PCM}$ (2)

where $I_{PCM}$ is the current through the PCM cell and $R_{PCM}$ is the
effective resistance of the cell. We use $R_{PCM}$ = 10K$\Omega$ in the low
resistance (SET) state and 200K$\Omega$ in the high resistance (RESET) state.
A part of this generated heat is dissipated to the surrounding and this heat
dissipation is given by the Equation [34],

$W_{d}=-k\sum\Delta T$ (3)

where $\Delta T$ represents the temperature dispersion around the active
region and expressed as,

$\Delta T=\frac{\partial T_{PCM}}{\partial x}+\frac{\partial T_{PCM}}{\partial
y}+\frac{\partial T_{PCM}}{\partial z}$ (4)

For simplicity we assume that the heat is mainly dispersed along the thickness
of the cell and the temperature outside the dispersion region is close to the
temperature surrounding the cell. Therefore, Equation 3 can be written as
[31], [53],

$W_{d}=\frac{kV}{l^{2}}(T_{PCM}-T_{surrounding})$ (5)

where $l$ is the thickness of the GST material and $k$ is the thermal
conductivity. Substitution of Equations 2 and 5 in Equation 1 yields,

$\frac{dT_{PCM}}{dt}=\frac{W_{j}-W_{d}}{C\times V}$ (6)

Solving this ODE gives,

$T_{PCM}=\frac{I_{PCM}^{2}R_{PCM}l^{2}}{kV}-C_{1}exp\left(-\frac{kt}{l^{2}C}\right)+T_{surrounding}$
(7)

Initially the PCM cell’s temperature is assumed to be the same as its
surrounding temperature. This boundary condition is used to determine the
constant $C_{1}$. Finally the cell temperature is modeled as,

$T_{PCM}=\frac{I_{PCM}^{2}R_{PCM}l^{2}}{kV}-\left[1-exp\left(-\frac{kt}{l^{2}C}\right)\right]+T_{surrounding}$
(8)

The surrounding temperature $T_{surrounding}$ is computed as

$T_{surrounding}=T_{amb}+\sum_{j}k\cdot T_{PCM_{j}}/D_{j}$ (9)

where $D_{j}$ is the thermal distance of the PCM cell from its neighboring
cell $j$, $T_{PCM_{j}}$ is the temperature of the neighboring cell, and
$T_{amb}$ is the ambient temperature of the neuromorphic hardware.

Equations 8 and 9 combine the following effects — 1) temporal thermal effect
of accessing a PCM cell in a machine learning workload, 2) the spatial thermal
contributions from the neighboring cell based on their activation.

Finally, we use the PCM temperature $T_{PCM}$ to compute the leakage current
through the access transistor of the PCM cell using Equation 10, where the
fitting parameters $A$ and $\eta$, and the nominal parameters $I_{nominal}$
and $T_{nominal}$ are obtained using [36, 42, 21, 19, 18, 3, 15].

$I_{leakage}\approx A\cdot
I_{nominal}\left(T_{PCM}-T_{nominal}\right){}^{\eta}$ (10)

##  3 Proposed Neuron and Synapse Mapping Technique

Figure 6 shows an overview of the proposed neuron and synapse mapping
approach. A machine learning application is first simulated using PyCARL [2],
a framework for simulating SNN-based applications. PyCARL internally uses the
CARLsim [13] simulator to extract the precise spike times on every synaptic
element in the SNN for representative training data. These spike times,
together with the neuron and synapse information constitute the SNN workload
for the machine learning application.

Figure 6: Overview of the proposed technique.

Next, the SNN workload is clustered using a greedy clustering approach,
roughly based on the Kernighan-Lin Graph Partitioning algorithm of SpiNeMap
[29].  Each cluster is a collection of pre- and post-synaptic neurons,
synapses connecting these neurons, and the spike times on these synapses. From
the mapping perspective, each cluster maps to a crossbar in the hardware,
while the inter-cluster communication channels are mapped on the shared
interconnect of the hardware. Therefore, the clustering technique ensures that
the neurons and synapses of a cluster can fit onto the resources of the
crossbar. PyCARL clusters an SNN to minimize the inter-cluster communication.
This reduces the spike congestion on the shared interconnect which improves
application latency.

The final step in our approach is the cluster mapping to the hardware. To
describe this step, let ${G}({C,S})$ be the machine learning workload with set
${C}$ of clusters and a set ${S}$ of connections between the clusters. The
workload is to be executed on the hardware $H(\mathcal{T},L)$ with a set
$\mathcal{T}$ of tiles (each tile has one crossbar) and a set $L$ of links
between the tiles. The mapping of the application $G$ to the hardware $H$,
$\mathcal{M}=\\{m_{x,y}\\}$ is defined as

$\footnotesize m_{x,y}=\begin{cases}1&\text{if cluster }{c}_{x}\in{C}\text{ is
mapped to tile }{t}_{y}\in{T}\\\ 0&\text{otherwise}\end{cases}$ (11)

Input: $G,H$

Output: $\mathcal{M}$

1 for  _$i$ in $MaxIter$_ do

2 $\mathcal{M}_{\text{init}}$ = allocate clusters to crossbars randomly;

3 $T_{\text{init}}$ = CalculateAvgTemperature($\mathcal{M}_{\text{init}}$);

4 do

5 for  _$c_{x}\in C$_ do

6 for  _$t_{y}\in\mathcal{T}$_ do

$\mathcal{M}_{y}=\mathcal{M}_{\text{init}}\big{|}m_{x,z}=\begin{cases}1&\text{if
}z=y\\\ 0&\text{otherwise}\end{cases}$ /* Move $c_{x}$ to tile $t_{y}$ and
generate the new mapping $\mathcal{M}_{y}$. */

7 $T_{y}$ = CalculateAvgTemperature($\mathcal{M}_{y}$);

8

9 end for

$x_{\text{idx}}$ = argmin $\\{T_{y}\big{|}y\in 1,2,\cdots,|\mathcal{T}|\\}$/*
Find the index of the mapping with the minimum temperature. */

10 if  _$T_{y} <T_{\text{min}}$_ then

$T_{\text{min}}=T_{y}$ and $\mathcal{M}_{\text{min}}=\mathcal{M}_{y}$/* Update
the mapping if the average temperature reduces. */

11

12 end if

13

14 end for

15

16 while  _$T_{\text{min}} <T_{\text{init}}$_;

17 end for

Return $\mathcal{M}_{\text{init}}$

Algorithm 1 Generate neuron and synapse mapping $\mathcal{M}$ to minimize the
average temperature of crossbars.

Algorithm 1 provides the pseudo-code of the hill-climbing based average
temperature minimization algorithm. The algorithm takes the clustered
application $G$ and the neuromorphic hardware $H$ as input. The algorithm
returns the mapping of $G$ to $H$, which minimizes the average temperature of
the crossbars. The algorithm is iterated for $MaxIter$ iterations (outer loop
lines 1-16). For each iteration of the outer loop, the algorithm generates a
random allocation of the clusters to the tiles (line 2) and calculate the
average temperature (line 3). The routine CalculateAvgTemperature calculates
the temperature of each crossbar for a mapping $\mathcal{M}$ using the
iterative approach of Figure 5b, specifically utilizing Equations 8 & 9, and
return the maximum average temperature of all crossbars in the neuromorphic
hardware.

At each iteration of the Algorithm 1, a cluster is moved to one of the tiles
(line 7), computing the average temperature of this new mapping (line 8). The
one mapping that leads to reduction of the average temperature is retained as
the new mapping (lines 10-13) and the process is repeated for the next cluster
(5-14). Once every cluster is analyzed, the iteration is repeated (lines 4-15)
to check if the clusters can be remapped again to reduce the average
temperature. The user-defined parameter $MaxIter$ governs the convergence of
the algorithm.

Algorithm Complexity: The complexity of Algorithm 1 is calculated as follows.
Let the inner loop (lines 4-15) be executed $\zeta$ times on average. At each
of these iterations, the algorithm performs $|C|\times|\mathcal{T}|$
operations. Therefore, the complexity of Algorithm 1 is
$O(MaxIter\times\zeta\times|C|\times|\mathcal{T}|)$.

##  4 Evaluation

###  4.1 Evaluated Applications

We evaluated 10 machine learning applications that are representative of three
most commonly used neural network classes — convolutional neural network
(CNN), multi-layer perceptron (MLP), and recurrent neural network (RNN). Table
1 summarizes the topology, the number of neurons and synapses of these
applications, and their baseline accuracy.

Table 1: Applications used to evaluate the proposed technique. Class | Applications | Synapses | Neurons | Topology | Accuracy
---|---|---|---|---|---
CNN | LeNet [33] | 282,936 | 20,602 | CNN | 85.1%
AlexNet [30] | 38,730,222 | 230,443 | CNN | 90.7%
VGG16 [43] | 99,080,704 | 554,059 | CNN | 69.8 %
HeartClass [4, 16] | 1,049,249 | 153,730 | CNN | 63.7%
MLP | DigitRecogMLP | 79,400 | 884 | FeedForward (784, 100, 10) | 91.6%
EdgeDet [13] | 114,057 | 6,120 | FeedForward (4096, 1024, 1024, 1024) | 100%
ImgSmooth [13] | 9,025 | 4,096 | FeedForward (4096, 1024) | 100%
RNN | HeartEstm [14] | 66,406 | 166 | Recurrent Reservoir | 100%
VisualPursuit [28] | 163,880 | 205 | Recurrent Reservoir | 47.3%
R-DigitRecog [25] | 11,442 | 567 | Recurrent Reservoir | 83.6%

###  4.2 Hardware Models

We model the DYNAP-SE neuromorphic hardware [39] with the following
configurations.

  * •

A tiled array of 4 tiles, each with a 128x128 crossbar. There are 65,536
memristors per crossbar.

  * •

Spikes are digitized and communicated between cores through a mesh routing
network using the Address Event Representation (AER) protocol.

  * •

Each synaptic element is a PCM-based memristor.

Table 2 reports the hardware parameters of DYNAP-SE.

Table 2: Major simulation parameters extracted from [39]. Neuron technology |  32nm FD-SOI
---|---
Synapse technology |  PCM
Supply voltage |  1.0V
Energy per spike |  50pJ at 30Hz spike frequency
Energy per routing |  147pJ
Switch bandwidth |  1.8G. Events/s

###  4.3 Evaluated Techniques

We evaluate the following two approaches.

  * •

SpiNeMap [6]: This is a performance-oriented approach to map SNN-based
applications to neuromorphic hardware. This approach first generates clusters
of neurons and synapses, where each cluster can fit on to the resources of a
tile in the hardware. Then, it uses an optimization algorithm to place these
clusters to the hardware, maximizing performance of the machine learning
application on the hardware. Temperature gradients are not incorporated in the
mapping process.

  * •

Proposed: In this technique the neurons and synapses of an SNN are mapped to
the hardware considering the thermal gradient. It uses the clustering
technique of SpiNeMap to generate clusters of neurons and synapses, where each
cluster can fit on to the resources of a tile. The clusters are mapped to the
crossbar using a hill-climbing approach to minimize the average temperature.
This reduces the leakage power consumption.

###  4.4 Evaluated Metrics

We evaluate the following metrics.

  * •

Average Temperature: This is the average temperature of each crossbar in the
hardware. We report the highest average temperatures of all crossbars.

  * •

Leakage Power: This is the total leakage power consumed in the hardware.

  * •

Performance: This is the latency, i.e., the time it takes to execute each
model on hardware.

  * •

Compilation Time: This is the time it takes to generate the minimum
temperature mapping of an application for the hardware.

##  5 Results and Discussion

###  5.1 Average Temperature

Figure 7 compares the maximum average temperature of the crossbars for each
evaluated application on DYNAP-SE using SpiNeMap and the proposed technique.
We make the following two key observations.

Figure 7: Maximum average temperature of the crossbars on DYNAP-SE.

First, the maximum average temperature increases with model size. VGG16, which
has more neurons and synapses than AlexNet (see Table 1), results in higher
average temperature than AlexNet for both SpiNeMap and the proposed technique.
MLP-MNIST, on the other hand, have lower temperature than both these models
due to its lower model complexity. Although R-DigitRecog has comparatively
fewer neurons and synapses, the average temperature is much higher. This is
because R-DigitRecog has higher activation, i.e., spikes in its workload,
which increases the temperature. These results clearly demonstrate the
workload-dependent nature of the temperature obtained on the hardware. Second,
the temperature obtained using the proposed mapping technique is lower than
SpiNeMap by an average 11.4K (between 6.4K and 17K) for these 10 applications.
This reduction is because of the proposed hill climbing algorithm (Algorithm
1), which incorporates the thermal gradient in optimizing the mapping of
neurons and synapses to the crossbars of the hardware.

###  5.2 Leakage Power

Figure 8 compares the leakage power on DYNAP-SE for each evaluated application
using SpiNeMap and the proposed technique. The leakage power constitute
between 20%–30% (average 22.8%) of the total energy consumption in the
hardware.  Results are normalized with respect to the leakage power obtained
on the hardware using SpiNeMap. We observe that the leakage power obtained
using the proposed technique is lower than SpiNeMap by an average 52%. This
significant improvement in the leakage power is due to the reduction of the
average temperature of the crossbars, which we analyzed in Section 5.1.  This
reduction in leakage power results in a reduction of the total energy
consumption by 11%.

Figure 8: Normalized leakage power on DYNAP-SE.

###  5.3 Performance

Figure 9 compares the latency of SpiNeMap and the proposed technique on DYNAP-
SE for the evaluated applications. We observe that the latency of the proposed
technique is only 5% higher (average) than SpiNeMap. Although the optimization
objective of SpiNeMap (which is performance) is different from the
optimization objective of the proposed technique (which is temperature), the
proposed technique uses the clustering technique of SpiNeMap to first generate
clusters, minimizing the spike communication on the shared interconnect of the
hardware. This results in lower spike latency. Therefore, in the next step
when the proposed technique optimizes for temperature during placement of the
clusters to crossbars of the hardware, the latency is not significantly higher
than SpiNeMap.

Figure 9: Normalized latency on DYNAP-SE.

###  5.4 Thermal Model Validation

We validate our thermal model against 1) the thermal model of [56], which
models the temperature of a single PCM cell and 2) the detailed model of [20],
which performs a detailed layout-based thermal simulations. The individual PCM
cell model is fast. However, it does not incorporate the thermal contributions
from neighboring PCM cells in a crossbar. Therefore, this model is not
accurate. On the other hand, the model in [20] is accurate because it
incorporates the spatial thermal contributions. However, it takes 30 minutes
of wall clock time to perform each thermal simulation for a 128x128 crossbar.
Therefore, incorporating this model in Algorithm 1 to evaluate the temperature
of a mapping makes the exploration time infeasible. Instead, we validated our
spatial formulation (Equation 9) by incorporating this equation into the
framework of [20].

Figure 10 plots the peak temperature obtained using the model of [56] and the
proposed model (Equation 8 & 9) for each evaluated application on DYNAP-SE. We
observe that existing models such as [56] lead to underestimation of the peak
temperature by an average 1.6K for these applications. This is because they do
not incorporate the spatial dependency. Underestimation of temperature leads
to an underestimation of the leakage power consumption of the hardware.

Figure 10: Comparison of Peak temperature.

Figure 11 plots the spatial contribution obtained using the model of [20] and
the proposed model (Equation 9) for 10 synthetic applications. We observe that
the accuracy of the proposed spatial model is close to that of the detailed
model [20]. The spatial contribution obtained using Equation 9 is on average
8.2% lower than [20] (0.3K in absolute terms).

Figure 11: Comparison of spatial contribution.

These results validate the thermal model proposed in this work.

###  5.5 Compilation Time and Solution Tradeoff

Table 3 reports the compilation time and the average temperature obtained for
three different settings of the variable $MaxIter$. We observe that as
$MaxIter$ is increased, the average temperature reduces for all applications.
This is because with increase in the number of iterations, Algorithm 1 is able
to find a better solution. However, the compilation time also increases.
Finally, we observe that increasing $MaxIter$ from 100 to 1000 results in a
significant increase in compilation time with a minimal improvement of the
average temperature. We conclude that setting $MaxIter=100$ gives the best
trade-off in terms of compilation time and the solution quality.  User can use
this $MaxIter$ parameter to set a limit on the compilation time of their
algorithm by analyzing the complexity of their model against the ones we
evaluate (see Table 1).

Table 3: Compilation time and solution tradeoff. Application | $\mathbf{MaxIter=10}$ | $\mathbf{MaxIter=100}$ | $\mathbf{MaxIter=1000}$
---|---|---|---
Compilation | Avg. | Compilation | Avg. | Compilation | Avg.
Time | Temperature | Time | Temperature | Time | Temperature
(sec) | (K) | (sec) | (K) | (sec) | (K)
LeNet | 26 | 326.3 | 259 | 323.2 | 2641 | 322.2
AlexNet | 114 | 330.1 | 1144 | 327.6 | 11480 | 326.0
VGG16 | 241 | 344.6 | 2413 | 335.8 | 24180 | 335.3
HeartClass | 96 | 315.1 | 965 | 311.3 | 9699 | 309.9
MLP-MNIST | 14 | 319.7 | 149 | 313.2 | 1520 | 311.6
EdgeDet | 12 | 323.5 | 132 | 321.5 | 1337 | 320.8
ImgSmooth | 26 | 327.11 | 268 | 324.7 | 2740 | 322.8
HeartEstm | 12 | 328.2 | 125 | 321.8 | 1255 | 320.4
VisualPursuit | 27 | 329.1 | 284 | 319.2 | 2883 | 318.7
R-DigitRecog | 15 | 336.3 | 159 | 328.5 | 1615 | 327.9

##  6 Conclusions

We propose a technique to map the neurons and synapses of SNN-based machine
learning applications to neuromorphic hardware. Prior work in this space have
focused extensively on performance, with no consideration of the thermal
aspects and the associated leakage power problem in the hardware. Our
technique is based on two key contributions. First, we propose a new thermal
model of a crossbar incorporating contributions from the adjacent cells.
Second, we incorporate this thermal model in a hill-climbing approach to
minimize the average temperature across the crossbars of the hardware. We
evaluate our approach using 10 machine learning applications and show the
significant reduction of the average temperature of the hardware. By lowering
the average temperature, we also show a reduction of leakage power
consumption.

##  7 Acknowledgment

This work is supported by the National Science Foundation Faculty Early Career
Development Award CCF-1942697 (CAREER: Facilitating Dependable Neuromorphic
Computing: Vision, Architecture, and Impact on Programmability).

## References

  * [1] Balaji, A., Song, S., Das, A., Dutt, N., Krichmar, J., Kandasamy, N., Catthoor, F.: A framework to explore workload-specific performance and lifetime trade-offs in neuromorphic computing. CAL (2019)
  * [2] Balaji, A., Adiraju, P., Kashyap, H.J., Das, A., Krichmar, J.L., Dutt, N.D., Catthoor, F.: PyCARL: A PyNN interface for hardware-software co-simulation of spiking neural network. In: IJCNN (2020)
  * [3] Balaji, A., Corradi, F., Das, A., Pande, S., Schaafsma, S., Catthoor, F.: Power-accuracy trade-offs for heartbeat classification on neural networks hardware. JOLPE (2018)
  * [4] Balaji, A., Corradi, F., Das, A., Pande, S., Schaafsma, S., et al.: Power-accuracy trade-offs for heartbeat classification on neural networks hardware. JOLPE (2018)
  * [5] Balaji, A., Das, A.: A framework for the analysis of throughput-constraints of snns on neuromorphic hardware. In: ISVLSI (2019)
  * [6] Balaji, A., Das, A., Wu, Y., Huynh, K., Dell’anna, F.G., Indiveri, G., Krichmar, J.L., Dutt, N.D., Schaafsma, S., Catthoor, F.: Mapping spiking neural networks to neuromorphic hardware. TVLSI (2020)
  * [7] Balaji, A., Marty, T., Das, A., Catthoor, F.: Run-time mapping of spiking neural networks to neuromorphic hardware. JSPS (2020)
  * [8] Balaji, A., Song, S., Das, A., Krichmar, J., Dutt, N., Shackleford, J., Kandasamy, N., Catthoor, F.: Enabling resource-aware mapping of spiking neural networks via spatial decomposition. ESL (2020)
  * [9] Balaji, A., Wu, Y., Das, A., Catthoor, F., et al.: Exploration of segmented bus as scalable global interconnect for neuromorphic computing. In: GLSVLSI (2019)
  * [10] Burr, G.W., Shelby, R.M., et al.: Neuromorphic computing using non-volatile memory. Advances in Physics: X (2017)
  * [11] Catthoor, F., Mitra, S., Das, A., Schaafsma, S.: Very large-scale neuromorphic systems for biological signal processing. In: CMOS Circuits for Biological Sensing and Processing (2018)
  * [12] Chen, I.R., Pop, E.: Compact thermal model for vertical nanowire phase-change memory cells. TED (2009)
  * [13] Chou, T.S., Kashyap, H.J., Xing, J., Listopad, S., Rounds, E.L., Beyeler, M., Dutt, N., et al.: CARLsim 4: An open source library for large scale, biologically detailed spiking neural network simulation using heterogeneous clusters. In: IJCNN (2018)
  * [14] Das, A., Pradhapan, P., Groenendaal, W., Adiraju, P., Rajan, R., Catthoor, F., Schaafsma, S., Krichmar, J., Dutt, N., et al.: Unsupervised heart-rate estimation in wearables with Liquid states and a probabilistic readout. Neural Networks (2018)
  * [15] Das, A., Catthoor, F., Bourdoux, A., Gyselinckx, B.: Energy-efficient mapping of lte-a phy signal processing tasks on microservers. TGCN (2018)
  * [16] Das, A., Catthoor, F., et al.: Heartbeat classification in wearables using multi-layer perceptron and time-frequency joint distribution of ECG. In: CHASE (2018)
  * [17] Das, A., Kumar, A.: Dataflow-based mapping of spiking neural networks on neuromorphic hardware. In: GLSVLSI (2018)
  * [18] Das, A., Kumar, A., Veeravalli, B.: Communication and migration energy aware design space exploration for multicore systems with intermittent faults. In: DATE (2013)
  * [19] Das, A., Kumar, A., Veeravalli, B.: Communication and migration energy aware task mapping for reliable multiprocessor systems. FGCS (2014)
  * [20] Das, A., Kumar, A., Veeravalli, B.: Reliability and energy-aware mapping and scheduling of multimedia applications on multiprocessor systems. TPDS (2015)
  * [21] Das, A., Walker, M.J., Hansson, A., Al-Hashimi, B.M., Merrett, G.V.: Hardware-software interaction for run-time power optimization: A case study of embedded linux on multicore smartphones. In: ISLPED (2015)
  * [22] Das, A., Wu, Y., Huynh, K., Dell’Anna, F., Catthoor, F., Schaafsma, S.: Mapping of local and global synapses on spiking neuromorphic hardware. In: DATE (2018)
  * [23] Davies, M., Srinivasa, N., et al.: Loihi: A neuromorphic manycore processor with on-chip learning. IEEE Micro (2018)
  * [24] Debole, M.V., Taba, B., et al.: TrueNorth: Accelerating from zero to 64 million neurons in 10 years. Computer (2019)
  * [25] Diehl, P.U., Cook, M.: Unsupervised learning of digit recognition using spike-timing-dependent plasticity. Frontiers in Computational Neuroscience (2015)
  * [26] Fouda, M.E., Eltawil, A.M., Kurdahi, F.: Modeling and analysis of passive switching crossbar arrays. TCAS I: Regular Papers (2017)
  * [27] Gubbi, J., Buyya, R., Marusic, S., Palaniswami, M.: Internet of Things (IoT): A vision, architectural elements, and future directions. FGCS (2013)
  * [28] Kashyap, H.J., et al.: A recurrent neural network based model of predictive smooth pursuit eye movement in primates. In: IJCNN (2018)
  * [29] Kernighan, B.W., Lin, S.: An efficient heuristic procedure for partitioning graphs. Bell system technical journal (1970)
  * [30] Krizhevsky, A., Sutskever, I., Hinton, G.E.: Imagenet classification with deep convolutional neural networks. In: NeurIPS (2012)
  * [31] Kwong, K.C., Li, L., He, J., Chan, M.: Verilog-A model for phase change memory simulation. ICSICT (2008)
  * [32] Le Gallo, M., Athmanathan, A., Krebs, D., Sebastian, A.: Evidence for thermally assisted threshold switching behavior in nanoscale phase-change memory cells. Journal of Applied Physics (2016)
  * [33] LeCun, Y., et al.: Lenet-5, convolutional neural networks (2015)
  * [34] Liao, Y.B., Chen, Y.K., Chiang, M.H.: An analytical compact PCM model accounting for partial crystallization. EDSSC (2007)
  * [35] Liao, Y.B., Lin, J.T., Chiang, M.H.: Temperature-based phase change memory model for pulsing scheme assessment. ICICDT (2008)
  * [36] Liu, Y., Dick, R.P., Shang, L., Yang, H.: Accurate temperature-dependent integrated circuit leakage power estimation is easy. In: DATE (2007)
  * [37] Maass, W.: Networks of spiking neurons: The third generation of neural network models. Neural Networks (1997)
  * [38] Mallik, A., Garbin, D., Fantini, A., Rodopoulos, et al.: Design-technology co-optimization for OxRRAM-based synaptic processing unit. In: VLSIT (2017)
  * [39] Moradi, S., Qiao, N., Stefanini, F., Indiveri, G.: A scalable multicore architecture with heterogeneous memory structures for dynamic neuromorphic asynchronous processors (DYNAPs). TBCAS (2017)
  * [40] Morikawa, T., Kurotsuchi, K., Kinoshita, M., et al.: Doped in-ge-te phase change memory featuring stable operation and good data retention. In: IEDM (2007)
  * [41] Ramasubramanian, S.G., Venkatesan, R., Sharad, M., et al.: Spindle: Spintronic deep learning engine for large-scale neuromorphic computing. In: ISLPED (2014)
  * [42] Shafik, R.A., Das, A., Yang, S., Merrett, G., Al-Hashimi, B.M.: Adaptive energy minimization of openmp parallel applications on many-core systems. In: PARMA-DITAM (2015)
  * [43] Simonyan, K., Zisserman, A.: Very deep convolutional networks for large-scale image recognition. arXiv (2014)
  * [44] Song, S., Balaji, A., Das, A., Kandasamy, N., Shackleford, J.: Compiling spiking neural networks to neuromorphic hardware. In: LCTES (2020)
  * [45] Song, S., Das, A.: A case for lifetime reliability-aware neuromorphic computing. In: MWSCAS (2020)
  * [46] Song, S., Das, A., Kandasamy, N.: Exploiting inter- and intra-memory asymmetries for data mapping in hybrid tiered-memories. In: ISMM (2020)
  * [47] Song, S., Das, A., Kandasamy, N.: Improving dependability of neuromorphic computing with non-volatile memory. In: EDCC (2020)
  * [48] Song, S., Das, A., Mutlu, O., Kandasamy, N.: Enabling and exploiting partition-level parallelism (PALP) in phase change memories. TECS (2019)
  * [49] Song, S., Das, A., Mutlu, O., Kandasamy, N.: Improving phase change memory performance with data content aware access. In: ISMM (2020)
  * [50] Song, S., Das, A., Mutlu, O., Kandasamy, N.: Aging aware request scheduling for non-volatile main memory. In: ASP-DAC (2021)
  * [51] Titirsha, T., Das, A.: Reliability-performance trade-offs in neuromorphic computing. In: CUT (2020)
  * [52] Warren, R., Reifenberg, J., Goodson, K.: Compact thermal model for phase change memory nanodevices. In: ICTTPES (2008)
  * [53] Wei, Y., Lin, X., Jia, Y., Cui, X., He, J., Zhang, X.: A SPICE model for a phase-change memory cell based on the analytical conductivity model. JOS (2012)
  * [54] Wong, H.S.P., Raoux, S., Kim, S., Liang, J., Reifenberg, J.P., Rajendran, B., Asheghi, M., et al.: Phase change memory. Proc. of the IEEE (2010)
  * [55] Woo, J., Yu, S.: Resistive memory-based analog synapse: The pursuit for linear and symmetric weight update. IEEE Nanotechnology Magazine (2018)
  * [56] Xi, L., Zhitang, S., Daolin, C., Xiaogang, C., Houpeng, C.: An spice model for phase-change memory simulations. JOS (2011)