11institutetext: Department of Astronomy and Theoretical Physics, Lund
University, Sölvegatan 14 A, 223 62 Lund, Sweden 22institutetext: Albert
Einstein Center for Fundamental Physics, Institute for Theoretical Physics,
Universität Bern, Sidlerstrasse 5, 3012 Bern, Switzerland 33institutetext:
Scuola Internazionale Superiore di Studi Avanzati (SISSA), Via Bonomea 265,
34136, Trieste, Italy

# Short-distance constraints on the hadronic light-by-light

Johan Bijnens 11<EMAIL_ADDRESS>Nils Hermansson-Truedsson ,
_Speaker_ 1122<EMAIL_ADDRESS>Antonio Rodríguez-Sánchez
33<EMAIL_ADDRESS>
###### Abstract

The muon anomalous magnetic moment continues to attract interest due to the
potential tension between experimental measurement Muong-2:2006rrc ;
Muong-2:2021ojo  and the Standard Model prediction Aoyama:2020ynm . The
hadronic light-by-light contribution to the magnetic moment is one of the two
diagrammatic topologies currently saturating the theoretical uncertainty. With
the aim of improving precision on the hadronic light-by-light in a data-driven
approach founded on dispersion theory Colangelo:2015ama ; Colangelo:2017fiz ,
we derive various short-distance constraints of the underlying correlation
function of four electromagnetic currents. Here, we present our previous
progress in the purely short-distance regime and current efforts in the so-
called Melnikov-Vainshtein limit.

##  1 Introduction

The muon anomalous magnetic moment, or, $a_{\mu}=(g-2)_{\mu}/2$, is a
potential probe for new physics beyond the Standard Model (SM) Aoyama:2020ynm
. Comparing the most recent SM prediction from the White Paper Aoyama:2020ynm
with the experimental average from Brookhaven National Laboratory
Muong-2:2006rrc  and Fermilab National Laboratory Muong-2:2021ojo  yields a
discrepancy of $4.2\sigma$. The corresponding numbers are

$\displaystyle a_{\mu}^{\textrm{SM}}=116\,591\,810(43)\times 10^{-11}\,,$ (1)
$\displaystyle a_{\mu}^{\textrm{exp}}=116\,592\,061(41)\times 10^{-11}\,.$ (2)

The SM uncertainty is dominated by hadronic contributions, namely the hadronic
vacuum polarisation (HVP) and the hadronic light-by-light (HLbL), which can be
calculated either in lattice field theory or using a data-driven approach
Aoyama:2020ynm . Although the lattice and data-driven predictions for the HLbL
agree within uncertainties, there has since the publication of the White Paper
arisen a tension for the HVP Borsanyi:2020mff . The updated lattice prediction
for the HVP predicts a significantly reduced discrepancy for the muon magnetic
moment Borsanyi:2020mff . There is currently much effort in trying to resolve
the tension between the two approaches to calculate the HVP.

In the following, we will be concerned with improving the precision on the
HLbL contribution, whose White Paper average is Aoyama:2020ynm

$\displaystyle a_{\mu}^{\textrm{HLbL}}=92(18)\times 10^{-11}\,.$ (3)

Diagrammatically, this contribution is given by the topology in Fig. 1 where
the photons have momenta $q_{1}$, $q_{2}$, $q_{3}$ and $q_{4}$. The external
momentum $q_{4}$ is for the $(g-2)_{\mu}$ soft, i.e. in the kinematical limit
$q_{4}\rightarrow 0$. The remaining three photons have virtual momenta
integrated over, which means that in the evaluation of the HLbL there is an
intricate mixing between different energy scales. In the dispersive data-
driven approach Colangelo:2015ama ; Colangelo:2017fiz , one uses analyticity
and unitarity to systematically decompose the HLbL into a sum over hadronic
contributions. The appearing form-factors in this decomposition require input
from e.g. experiments or the lattice Aoyama:2020ynm . The unknown hadronic
contributions can, however, be controlled with the help of so-called short-
distance constraints (SDCs) on the underlying correlation function of four
electromagnetic currents, see Refs. Melnikov:2003xd ; Bijnens:2019ghy ;
Bijnens:2020xnl ; Bijnens:2021jqo . These kinds of SDCs can be obtained from
operator product expansion (OPE) techniques and have been used for model
calculations Bijnens:1995xf ; Colangelo:2019lpu ; Colangelo:2019uex ;
Cappiello:2019hwh ; Leutgeb:2019gbz ; Colangelo:2021nkr  and a complementary
approach building on interpolation techniques Ludtke:2020moa .

Figure 1: The HLbL topology. The grey blob connecting to the external muon
line contains all hadronic contributions.

Defining the Euclidean virtualities $Q_{i}^{2}=-q_{i}^{2}$ for $i=1,2,3$, we
here consider SDCs in the purely short-distance region
$Q_{i}^{2}\gg\Lambda_{\textrm{QCD}}^{2}$ (based on Refs. Bijnens:2019ghy ;
Bijnens:2020xnl ; Bijnens:2021jqo ), and preliminary work in the Melnikov-
Vainshtein limit $Q_{i}^{2},Q_{j}^{2}\gg
Q_{k}^{2},\Lambda^{2}_{\textrm{QCD}}$. These two expansions respectively
correspond to three and two of the electromagnetic currents in the underlying
correlation function being close.

##  2 Some generalities

The HLbL is defined in terms of a correlation function of four electromagnetic
currents. The corresponding so-called HLbL tensor is given by

$\displaystyle\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}(q_{1},q_{2},q_{3})=-i\int\frac{d^{4}q_{4}}{(2\pi)^{4}}\left(\prod_{i=1}^{4}\int
d^{4}x_{i}\,e^{-iq_{i}x_{i}}\right)\,\langle
0|T\left(\prod_{j=1}^{4}J^{\mu_{j}}(x_{j})\right)|0\rangle\,.$ (4)

Here, the currents are given by $J^{\mu}(x)=\bar{q}\,Q_{q}\gamma^{\mu}q$ where
$q=(u,d,s)$ and $Q_{q}=\textrm{diag}(e_{q})=\textrm{diag}(2/3,-1/3,-1/3)$ is
the associated charge matrix for light quarks. Here we use the convention
$q_{1}+q_{2}+q_{3}+q_{4}=0$. The above tensor satisfies the Ward identities

$q_{i,\,\mu_{i}}\,\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}(q_{1},q_{2},q_{3})=0\,.$
(5)

This implies the relation Aldins:1970id

$\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}(q_{1},q_{2},q_{3})=-q_{4,\,\nu_{4}}\frac{\partial\Pi^{\mu_{1}\mu_{2}\mu_{3}\nu_{4}}}{\partial
q_{4,\,\mu_{4}}}(q_{1},q_{2},q_{3})\,,$ (6)

which means that one can access $a_{\mu}^{\textrm{HLbL}}$ from the derivative
tensor on the right-hand side. The derivative can be Lorentz decomposed into a
set of 54 scalar functions $\Pi_{i}$ Colangelo:2015ama ; Colangelo:2017fiz .
The contribution to the magnetic moment can be written in terms of six linear
combinations $\hat{\Pi}_{1,4,7,17,39,54}$ of the $\Pi_{i}$, namely

$\displaystyle
a_{\mu}^{\mathrm{HLbL}}=\frac{2\alpha^{3}}{3\pi^{2}}\int_{0}^{\infty}dQ_{1}\int_{0}^{\infty}dQ_{2}\int_{-1}^{1}d\tau\,\sqrt{1-\tau^{2}}\,Q_{1}^{3}Q_{2}^{3}\sum_{i=1}^{12}T_{i}(Q_{1},Q_{2},\tau)\,\overline{\Pi}_{i}(Q_{1},Q_{2},\tau)\,,$
(7)

where $\overline{\Pi}_{i}$ are functions of the six $\hat{\Pi}_{j}$, and the
kernels $T_{i}$ are known. Note that the integration in (7) is for the full
HLbL contribution, whereas the kinematic regions for the two types of SDCs
considered here correspond to restricted integration domains. Common for the
two limits $Q_{i}^{2}\gg\Lambda_{\textrm{QCD}}^{2}$ and
$Q_{i}^{2},Q_{j}^{2}\gg Q_{k}^{2},\Lambda^{2}_{\textrm{QCD}}$ is the need to
define an onset $Q_{\textrm{min}}$ of the asymptotic limits:
$Q_{i}^{2}>Q_{\textrm{min}}^{2}$ and
$Q_{i}^{2},Q_{j}^{2}>Q_{\mathrm{min}}^{2}$, respectively. Since there
currently is no preferred choice of $Q_{\textrm{min}}$ we keep it variable,
and note that this ambiguity induces a source of uncertainty in
$a_{\mu}^{\mathrm{HLbL}}$.

##  3 An OPE for three currents

In the purely short-distance limit $Q_{i}^{2}\gg\Lambda_{\textrm{QCD}}^{2}$,
three of the currents in the four-point function defined in (4) are close.
Naively performing an OPE of these currents works only to leading order, since
the next-to-leading order contribution in the systematic expansion is ill-
defined in the static limit $q_{4}\rightarrow 0$, which is the limit relevant
for the magnetic moment Bijnens:2019ghy . The problem arises from propagators
of the momentum $q_{4}$ which clearly diverge. However, it was shown in Ref.
Bijnens:2019ghy  that the problem can be circumvented by treating the static
photon as an external electromagnetic field and instead of the four-point
function $\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}$ consider a three-point
correlator in the presence of that field. For this reason we study

$\displaystyle\Pi^{\mu_{1}\mu_{2}\mu_{3}}(q_{1},q_{2})=-\frac{1}{e}\int\frac{d^{4}q_{3}}{(2\pi)^{4}}\left(\prod_{i=1}^{3}\int
d^{4}x_{i}\,e^{-iq_{i}x_{i}}\right)\,\langle
0|T\left(\prod_{j=1}^{3}J^{\mu_{j}}(x_{j})\right)|\gamma(q_{4})\rangle\,,$ (8)

where the static background field is captured in the external state. This is
related to the full HLbL tensor through Bijnens:2019ghy

$\displaystyle\Pi^{\mu_{1}\mu_{2}\mu_{3}}(q_{1},q_{2})=\epsilon_{\mu_{4}}(q_{4})\,\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}(q_{1},q_{2},q_{3})=-\epsilon_{\mu_{4}}(q_{4})\,q_{4,\,\nu_{4}}\frac{\partial\Pi^{\mu_{1}\mu_{2}\mu_{3}\nu_{4}}}{\partial
q_{4,\,\mu_{4}}}\,(q_{1},q_{2},q_{3})\,,$ (9)

where $\epsilon_{\mu_{4}}(q_{4})$ is the polarisation vector of the external
photon and the Ward identity (6) was used. By performing a systematic
background field OPE between the three currents in (8) then yields the HLbL
derivative through (9). The background field induces non-perturbative
condensates in the expansion which are different than those appearing in
vacuum OPEs Shifman:1978bx . The condensates appear from diagrammatic
contributions to the correlator in (8) where some of the fields are left
uncontracted in the Wick expansion. Background field OPEs were first used for
nucleon magnetic moments in Refs. Balitsky:1983xk ; Ioffe:1983ju , but were in
fact later also employed for the electroweak contribution to $a_{\mu}$ in Ref.
Czarnecki:2002nt .

###  3.1 Leading order perturbative and non-perturbative corrections

The systematic OPE is done in powers of $1/Q_{i}$ and the strong coupling
$\alpha_{s}$. Through order $1/Q_{i}^{6}$ and $\alpha_{s}$ the leading order
term and non-perturbative corrections arise from the operators Bijnens:2019ghy
; Bijnens:2020xnl

$\displaystyle S_{1,\,\mu\nu}$ $\displaystyle=e\,e_{q}F_{\mu\nu}\,,\qquad
S_{2,\,\mu\nu}=\bar{q}\,\sigma_{\mu\nu}q\,,\hskip 17.0pt\quad
S_{3,\,\mu\nu}=i\,\,\bar{q}\,G_{\mu\nu}q\,,$ (10) $\displaystyle
S_{4,\,\mu\nu}$ $\displaystyle=i\,\bar{q}\,\bar{G}_{\mu\nu}\gamma_{5}q\,,\quad
S_{5,\,\mu\nu}=\bar{q}q\;e\,e_{q}F_{\mu\nu}\,,\qquad
S_{6,\,\mu\nu}=\frac{\alpha_{s}}{\pi}\,G_{a}^{\alpha\beta}G^{a}_{\alpha\beta}\;e\,e_{q}F_{\mu\nu}\,,$
(11) $\displaystyle S_{7,\,\mu\nu}$
$\displaystyle=\bar{q}\,\left(G_{\mu\lambda}D_{\nu}+D_{\nu}G_{\mu\lambda}\right)\gamma^{\lambda}q-(\mu\leftrightarrow\nu)\,,$
(12) $\displaystyle S_{\\{8\\},\,\mu\nu}$
$\displaystyle=\alpha_{s}\,\left(\bar{q}\,\Gamma q\;\bar{q}\,\Gamma
q\right)_{\mu\nu}\,.$ (13)

Here, the gluon field is given by $G_{\mu\nu}=ig_{S}\lambda^{a}G^{a}_{\mu\nu}$
as well as its dual
$\bar{G}^{\mu\nu}=\frac{i}{2}\epsilon^{\mu\nu\lambda\rho}G_{\lambda\rho}$. The
matrix $\Gamma$ is in flavour, spin and colour space . For more detail, we
refer the reader to Ref. Bijnens:2020xnl . The class of four-quark operators
$S_{\\{8\\},\,\mu\nu}$ yields flavour mixing, and in the chiral limit there
are only 12 independent operators. In Fig. 2 we show the various diagrammatic
contributions appearing in the OPE. The leading order term is given by the
perturbative quark-loop in Fig. 2(a), which in fact coincides with the leading
order term in a naive OPE of the four-point function in (4). This finding in
Ref. Bijnens:2019ghy  confirmed that indeed the perturbative quark loop
describes the leading behaviour in the short-distance limit Aoyama:2020ynm .

Figure 2: Various contributions to the OPE. The leading order contribution is
$(a)$ the quark loop. The crossed vertex indicates the interaction of the
external field on the perturbative quark line. At next-to leading order the
non-perturbative condensates $(b_{1,2})$ appear, namely
$\langle\overline{q}\,\sigma_{\mu\nu}q\rangle$, induced by $F_{\mu\nu}$, and
$\langle\overline{q}q\rangle$. At next-to-next-to leading order diagram $(c)$
is a four-quark operator with condensate
$\langle\overline{q}\,\Gamma_{1}q\,\overline{q}\,\Gamma_{2}q\rangle$, and
diagram $(d)$ contains the gluon condensate $\langle\alpha_{s}GG\rangle$.

The obtained systematic expansion allows us to quantify for the first time
from where the short-distance representation of the HLbL is valid. We study
first the impact of the non-perturbative corrections. We note in passing that
one also has to renormalise the condensates, which is described in detail in
Ref. Bijnens:2020xnl . It suffices here to write the renormalised result on
the form

$\Pi^{\mu_{1}\mu_{2}\mu_{3}}(q_{1},q_{2})=\frac{1}{e}\,\vec{C}^{\,T,\,\mu_{1}\mu_{2}\mu_{3}\mu_{4}\nu_{4}}_{\overline{MS}}(q_{1},q_{2})\,\langle\vec{Q}_{\overline{MS},\,\mu_{4}\nu_{4}}(\mu)\rangle\,,$
(14)

where the perturbative short-distance Wilson coefficients are contained in the
vector $\vec{C}^{\,T,\,\mu_{1}\mu_{2}\mu_{3}\mu_{4}\nu_{4}}_{\overline{MS}}$
and the condensates at renormalisation scale $\mu$ in
$\langle\vec{Q}_{\overline{MS},\mu_{4}\nu_{4}}(\mu)\rangle$. For convenience,
we write

$\langle\vec{Q}_{\overline{MS},\,\mu\nu}(\mu)\rangle=e\,\vec{X}\,\langle
e_{q}F_{\mu\nu}\rangle\,.$ (15)

Numerical estimates for the condensates are obtained from Refs.
Bijnens:2020xnl ; Shifman:1978bx ; Aoki:2019cca ; Belyaev:1982sa .

Analytical formulae for the $\hat{\Pi}_{i}$ needed to numerically calculate
$a_{\mu}^{\mathrm{HLbL}}$ in (7) in the short-distance region
$Q_{i}^{2}>Q_{\textrm{min}}^{2}$ can be found in Refs. Bijnens:2019ghy ;
Bijnens:2020xnl . Performing the numerical integration yields Figs. 4–5. The
leading order perturbative quark loop corresponds to $X_{1,0}$, which as can
be seen is completely dominant by one to two orders of magnitude. This shows
that the non-perturbative corrections through next-to-next-to leading order
are negligible as compared to the perturbative quark loop.

Figure 3: Contributions from $X_{1,0}$, $X_{1,m^{2}}$, $X_{2,m}$ and
$X_{2,m^{3}}$ to $a_{\mu}^{\mathrm{HLbL}}$.

Figure 4: Contributions from $X_{1,0}$, $X_{6}$, $X_{7}$, $X_{8,1}$ and
$X_{8,2}$ to $a_{\mu}^{\mathrm{HLbL}}$.

Figure 5: Contributions from $X_{1,0}$, $X_{3}$, $X_{4}$ and $X_{5}$ to
$a_{\mu}^{\mathrm{HLbL}}$.

###  3.2 Gluonic corrections

From the above numerical results it appears that the perturbative quark loop
is a good description of the short-distance behaviour of the HLbL. However,
perturbative order $\alpha_{s}$ corrections to the quark loop can also be
important, which is what we studied in Ref. Bijnens:2021jqo . Including two
vertices from the Dyson series with gluon interactions yields two-loop
topologies like those in Fig. 7. Details of the calculation are given in Ref.
Bijnens:2021jqo  together with the $\hat{\Pi}_{i}$, and we report here only
the numerical comparison to the quark loop contribution to
$a_{\mu}^{\mathrm{HLbL}}$ in Fig. 7. As can be seen, the shift in going from
the leading order quark loop to including the two-loop corrections is small.
Numerically, it corresponds to a shift around $-10\%$ compared to the quark
loop. This indeed shows that the quark loop describes the short-distance
dynamics of the HLbL to within 10% in the region
$Q_{i}^{2}\gg\Lambda_{\textrm{QCD}}^{2}$, and can be used by dispersive and
model studies as in Ref. Colangelo:2021nkr .

Figure 6: Two-loop gluonic corrections to the perturbative quark loop.

Figure 7: Contributions to $a_{\mu}^{\mathrm{HLbL}}$ from gluonic corrections
to the perturbative quark loop. The uncertainty band is associated to
$\alpha_{s}\left(\mu=Q_{\textrm{min}}\right)$.

##  4 An OPE for two currents

Having concluded the study for the purely short-distance region
$Q_{i}^{2}\gg\Lambda_{\textrm{QCD}}^{2}$ we now turn to the Melnikov-
Vainshtein limit $Q_{i}^{2},Q_{j}^{2}\gg
Q_{k}^{2},\Lambda^{2}_{\textrm{QCD}}$. This was studied first in Ref.
Melnikov:2003xd . We emphasise that results here are still preliminary, and we
exemplify by considering the case $Q_{1}^{2},Q_{2}^{2}\gg
Q_{3}^{2},\Lambda^{2}_{\textrm{QCD}}$. Recall that for three currents we
studied the three-point function in (8). In the Melnikov-Vainshtein limit we
can instead use the two-point function

$\Pi^{\mu_{1}\mu_{2}}=\frac{i}{e^{2}}\int\frac{d^{4}q_{4}}{(2\pi)^{4}}\int
d^{4}x_{1}\int d^{4}x_{2}e^{-i(q_{1}x_{1}+q_{2}x_{2})}\langle
0|T\left\\{J^{\mu_{1}}(x_{1})J^{\mu_{2}}(x_{2})\right\\}|\gamma^{*}(q_{3})\gamma(q_{4})\rangle\,.$
(16)

The relation to the four-point HLbL tensor is then (cf. (9))

$\Pi^{\mu_{1}\mu_{2}}=\epsilon_{\mu_{3}}(q_{3})\,\epsilon_{\mu_{4}}(q_{4})\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}=-\epsilon_{\mu_{3}}(q_{3})\,\epsilon_{\mu_{4}}(q_{4})\,q_{4,\,\nu_{4}}\frac{\partial\Pi^{\mu_{1}\mu_{2}\mu_{3}\nu_{4}}}{\partial
q_{4,\,\mu_{4}}}\,(q_{1},q_{2},q_{3})\,.$ (17)

In this OPE we will expand in powers of the large momentum
$\hat{q}=(q_{1}-q_{2})/2$. We further define a Dirac matrix in terms of the
massless quark propagator according to

$\displaystyle\Gamma^{\mu\nu}(k)=\gamma^{\mu}S(k)\gamma^{\nu}\,,\qquad
S(k)=\frac{\not{k}}{k^{2}}\,.$ (18)

Denoting the perturbative quark loop contribution to the HLbL tensor in (4) by
$\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}_{\mathrm{quark-loop}}$, one finds in the
chiral limit through next-to-leading order the result

$\displaystyle\Pi^{\mu_{1}\mu_{2}}$
$\displaystyle\approx{-\frac{e_{q}^{2}}{e^{2}}\,\left\langle\bar{q}(0)\left[\Gamma^{\mu_{1}\mu_{2}}(-\hat{q})-\Gamma^{\mu_{2}\mu_{1}}(-\hat{q})\right]q(0)\,\big{|}\gamma^{*}(q_{3})\gamma(q_{4})\right\rangle}$
$\displaystyle-\frac{ie_{q}^{2}}{e^{2}\hat{q}^{2}}\left(g^{\mu_{1}\delta}g^{\mu_{2}}_{\beta}+g^{\mu_{2}\delta}g^{\mu_{1}}_{\beta}-g^{\mu_{1}\mu_{2}}g^{\delta}_{\beta}\right)\left(g_{\alpha\delta}-2\frac{\hat{q}_{\delta}\hat{q}_{\alpha}}{\hat{q}^{2}}\right)\left\langle\bar{q}(0)\left[\overrightarrow{\partial}^{\alpha}-\overleftarrow{\partial}^{\alpha}\right]\gamma^{\beta}q(0)\big{|}\gamma^{*}(q_{3})\gamma(q_{4})\,\right\rangle$
$\displaystyle-\frac{1}{4}{F_{\nu_{3}\mu_{3}}F_{\nu_{4}\mu_{4}}}\left.\frac{\partial}{\partial
q_{3\,\nu_{3}}}\right|_{q_{3}\rightarrow 0}\left.\frac{\partial}{\partial
q_{4\,\nu_{4}}}\right|_{q_{4}\rightarrow
0}\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}_{\mathrm{quark-loop}}\,.$ (19)

We thus see that the two external photons in (16) induce new contributions as
compared to (8). There are several technical issues in obtaining the
$\hat{\Pi}_{i}$ in the Melnikov-Vainshtein limit which will be explained in
our upcoming paper. Note that we have made no ordering on the size of
$Q_{3}^{2}$ and $\Lambda_{\mathrm{QCD}}^{2}$ here, except for an omitted gluon
operator contribution to the flavour singlet piece when
$\Lambda_{\mathrm{QCD}}<Q_{3}$. However, the perturbative $Q_{3}^{2}$ regime
$Q_{1}^{2},Q_{2}^{2}\gg Q_{3}^{2}\gg\Lambda^{2}_{\textrm{QCD}}$ corresponds to
a special case of the purely short-distance limit
$Q_{i}^{2}\gg\Lambda_{\textrm{QCD}}^{2}$ discussed above. This means that for
$Q_{3}^{2}\gg\Lambda^{2}_{\textrm{QCD}}$ in the Melnikov-Vainshtein limit we
should reproduce our previous results for the perturbative quark loop as well
as the leading order result of Ref. Melnikov:2003xd . Evaluating the appearing
matrix elements in (4) perturbatively, this is indeed the case.

In the non-perturbative region $Q_{3}^{2}\not{\gg}\Lambda^{2}_{\textrm{QCD}}$,
the matrix elements in (4) must be form-factor decomposed. The leading order
term (first row in (4)) depends on two form-factors, i.e. the longitudinal and
transversal functions discussed in Ref. Melnikov:2003xd . At next-to-leading
order there are six form-factors, and the exact decomposition will be given in
our future paper.

We are currently working on gluonic corrections to the OPE in (4). As argued
in Ref. Ludtke:2020moa , the leading perturbative $\alpha_{s}$ correction is
expected to be given by the Melnikov-Vainshtein result Melnikov:2003xd  with
an overall factor $-\alpha_{s}/\pi$. In our OPE, through leading order in
$\alpha_{s}$ the first line in (4) instead yields

$\displaystyle\Pi^{\mu_{1}\mu_{2}}_{\mathrm{LO},\,\alpha_{s}}$
$\displaystyle\approx{-\frac{e_{q}^{2}}{e^{2}}\,\left(1-\frac{7}{3}\,\frac{\alpha_{s}}{\pi}\right)\,\left\langle\bar{q}(0)\left[\Gamma^{\mu_{1}\mu_{2}}(-\hat{q})-\Gamma^{\mu_{2}\mu_{1}}(-\hat{q})\right]q(0)\,\big{|}\gamma^{*}(q_{3})\gamma(q_{4})\right\rangle}\,.$
(20)

Although this might seem to be in tension, the Melnikov-Vainshtein result
refers to the axial current that preserves chiral symmetry within perturbation
theory, related to ours by a finite counterterm Trueman:1979en . We are
currently working to complete the calculation, taking these subtleties into
account.

SDCs in the Melnikov-Vainshtein limit are highly relevant for data-driven
approaches to evaluate the HLbL, see e.g. Refs. Melnikov:2003xd ;
Colangelo:2019lpu ; Colangelo:2019uex ; Ludtke:2020moa ; Colangelo:2021nkr .
Our anticipated results can be used to further improve the data-driven
prediction.

##  5 Conclusions

In these proceedings we have reported on recent progress on the derivation of
short-distance constraints for the data-driven evaluation of the HLbL. This
hadronic contribution to the muon magnetic moment is particularly complicated
since it mixes long- and short-distance dynamics. Constraints in the purely
short-distance as well as the Melnikov-Vainshtein limit are important to
control the systematic uncertainty in the data-driven approach to the HLbL,
both of which we have considered herein. Due to the soft external magnetic
field, the constraints are derived using background field operator product
expansion techniques. Our main result is a precise description for the HLbL in
the purely short-distance regime. In the Melnikov-Vainshtein limit, which is
still work in progress, we have reproduced important results in certain
limits, and are currently evaluating also novel higher-order non-perturbative
as well as gluonic corrections.

## Acknowledgements

N. H.-T. is funded in part by the Albert Einstein Center for Fundamental
Physics at the University of Bern, and in part by the Swedish Research
Council, project number 2021-06638. J. B. is supported by the Swedish Research
Council grants contract numbers 2016-05996 and 2019-03779. A. R.-S. is funded
in part by MIUR contract number 2017L5W2PT.

## References

  * (1) G.W. Bennett et al. (Muon g-2), Phys. Rev. D 73, 072003 (2006), hep-ex/0602035
  * (2) B. Abi et al. (Muon g-2), Phys. Rev. Lett. 126, 141801 (2021), 2104.03281
  * (3) T. Aoyama et al., Phys. Rept. 887, 1 (2020), 2006.04822
  * (4) G. Colangelo, M. Hoferichter, M. Procura, P. Stoffer, JHEP 09, 074 (2015), 1506.01386
  * (5) G. Colangelo, M. Hoferichter, M. Procura, P. Stoffer, JHEP 04, 161 (2017), 1702.07347
  * (6) S. Borsanyi et al., Nature 593, 51 (2021), 2002.12347
  * (7) K. Melnikov, A. Vainshtein, Phys. Rev. D70, 113006 (2004), hep-ph/0312226
  * (8) J. Bijnens, N. Hermansson-Truedsson, A. Rodríguez-Sánchez, Phys. Lett. B798, 134994 (2019), 1908.03331
  * (9) J. Bijnens, N. Hermansson-Truedsson, L. Laub, A. Rodríguez-Sánchez, JHEP 10, 203 (2020), 2008.13487
  * (10) J. Bijnens, N. Hermansson-Truedsson, L. Laub, A. Rodríguez-Sánchez, JHEP 04, 240 (2021), 2101.09169
  * (11) J. Bijnens, E. Pallante, J. Prades, Nucl. Phys. B 474, 379 (1996), hep-ph/9511388
  * (12) G. Colangelo, F. Hagelstein, M. Hoferichter, L. Laub, P. Stoffer, Phys. Rev. D 101, 051501 (2020), 1910.11881
  * (13) G. Colangelo, F. Hagelstein, M. Hoferichter, L. Laub, P. Stoffer, JHEP 03, 101 (2020), 1910.13432
  * (14) L. Cappiello, O. CatÃ , G. D’Ambrosio, D. Greynat, A. Iyer, Phys. Rev. D 102, 016009 (2020), 1912.02779
  * (15) J. Leutgeb, A. Rebhan, Phys. Rev. D 101, 114015 (2020), 1912.01596
  * (16) G. Colangelo, F. Hagelstein, M. Hoferichter, L. Laub, P. Stoffer, Eur. Phys. J. C 81, 702 (2021), 2106.13222
  * (17) J. Lüdtke, M. Procura, Eur. Phys. J. C 80, 1108 (2020), 2006.00007
  * (18) J. Aldins, T. Kinoshita, S.J. Brodsky, A.J. Dufner, Phys. Rev. D1, 2378 (1970)
  * (19) M.A. Shifman, A.I. Vainshtein, V.I. Zakharov, Nucl. Phys. B147, 385 (1979)
  * (20) I.I. Balitsky, A.V. Yung, Phys. Lett. 129B, 328 (1983)
  * (21) B.L. Ioffe, A.V. Smilga, Nucl. Phys. B232, 109 (1984)
  * (22) A. Czarnecki, W.J. Marciano, A. Vainshtein, Phys. Rev. D67, 073006 (2003), [Erratum: Phys. Rev. D73, 119901 (2006)], hep-ph/0212229
  * (23) S. Aoki et al. (Flavour Lattice Averaging Group), Eur. Phys. J. C80, 113 (2020), 1902.08191
  * (24) V. Belyaev, B. Ioffe, Sov. Phys. JETP 56, 493 (1982)
  * (25) T.L. Trueman, Phys. Lett. B 88, 331 (1979)


# ST-Mamba: Spatial-Temporal Mamba for Traffic Flow Estimation Recovery using
Limited Data

Doncheng Yuan1, Jianzhe Xue1, Jinshan Su2, Wenchao Xu3, and Haibo Zhou1
1School of Electronic Science and Engineering, Nanjing University, Nanjing,
China, 210023.
2Key Laboratory of Vibration Signal Capture and Intelligent Processing, Yili
Normal University, Yining, China, 835000.
3Department of Computing, the Hong Kong Polytechnic University.
Email<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
and<EMAIL_ADDRESS>
###### Abstract

Traffic flow estimation (TFE) is crucial for urban intelligent traffic
systems. While traditional on-road detectors are hindered by limited coverage
and high costs, cloud computing and data mining of vehicular network data,
such as driving speeds and GPS coordinates, present a promising and cost-
effective alternative. Furthermore, minimizing data collection can
significantly reduce overhead. However, limited data can lead to inaccuracies
and instability in TFE. To address this, we introduce the spatial-temporal
Mamba (ST-Mamba), a deep learning model combining a convolutional neural
network (CNN) with a Mamba framework. ST-Mamba is designed to enhance TFE
accuracy and stability by effectively capturing the spatial-temporal patterns
within traffic flow. Our model aims to achieve results comparable to those
from extensive data sets while only utilizing minimal data. Simulations using
real-world datasets have validated our model’s ability to deliver precise and
stable TFE across an urban landscape based on limited data, establishing a
cost-efficient solution for TFE.

###### Index Terms:

ST-Mamba, traffic flow estimation, recovery method, limited data.

##  I Introduction

Traffic flow estimation (TFE) assumes paramount importance in the operation of
traffic systems, particularly in light of the escalating volume of vehicles
and the expansion of transportation networks [1] [2]. The conventional methods
of TFE, employing on-road detectors like cameras and inductive loops, exhibit
limited coverage and substantial expenses associated with infrastructure
installation and maintenance [3]. Utilizing data from vehicular network can
expand the scope of TFE, while the collection and transmission of massive data
also encounter technical and economic challenges. Considering the constraints
of existing methods, it is cost-effective to utilize limited data from
vehicular network to obtain real-time TFE by cloud computing [4][5]. However,
when the amount of traffic data decreases, the resulting TFE tends to be
inaccurate and unstable [6]. With the advancement of artificial intelligence,
employing deep learning methods to recover TFE from limited data emerges as a
feasible and promising approach.

Lately, deep learning methods have been developed to explore features of
traffic flow [7] [8]. The convolutional neural network (CNN) is proficient in
extracting spatial features of input through convolutional operations [9]. The
recurrent neural network (RNN) and its variants such as long short-term memory
(LSTM) are designed to capture sequential or temporal features of input by
utilizing recurrent connections [10]. Mamba, an advanced neural network based
on state space model (SSM), is effective at processing sequential data and
modeling temporal correlation via selection mechanism [11]. Moreover, some
methods combine the functions of CNN and RNN. For instance, the Conv-LSTM
module is proposed to forecast short-term traffic flow via learning spatial-
temporal dependencies [12]. However, existing works primarily focus on traffic
flow prediction rather than estimating real-time traffic flow, and are
provided with sufficient data instead of limited data. To enhance estimation
accuracy with limited data, there is a need to design a specialized spatial-
temporal model for traffic flow estimation recovery.

In this paper, we investigate the TFE utilizing limited data and propose a
spatial-temporal deep learning method to recover both the accuracy and
stability of TFE. The framework of TFE with limited vehicular network data is
illustrated in Fig. 1. We divide the city map evenly into grids, each of which
represent a specific region in the city. The data limitation in our scheme
refers to the practice of randomly sampling a subset of vehicles across the
city with equal probability of recruitment at each time interval [6].
Following the collection of limited data, the traffic information comprising
vehicle speeds and GPS coordinates are matched to corresponding grids. The
condition of traffic flow within each grid is represented by the average speed
of vehicles over a defined period. For the recovery of TFE, we propose a novel
deep learning model, named as spatial-temporal Mamba (ST-Mamba), to generate
stable and accurate TFE via effectively leveraging the spatial-temporal
features of traffic flow. Specifically, we employ the CNN to capture the local
spatial correlation of traffic flow, and utilize the Mamba to model the
temporal correlation of each grid in the map. The main contributions of this
paper are listed as the following:

  * •

We present a novel approach for cost-effective TFE, utilizing only a small
fraction of vehicular network data to alleviate the burden on network
communication and mitigate the expenses of traffic systems.

  * •

We design a spatial-temporal deep learning model named as ST-Mamba, consisting
of CNN and Mamba, to recover the TFE with limited vehicular network data.

  * •

We validate the effectiveness of the proposed method through comprehensive
simulations on the real-world dataset, varying in data limitation from 10% to
50%. The results demonstrate that our method can provide accurate and stable
city-wide TFE.

The rest of this paper is organized as follows. Section II demonstrates the
framework and problem analysis with limited vehicular network data. The
structure of ST-Mamba is detailed in Section III. Then, section IV showcases
the simulation performance of our proposed model. At last, we conclude this
paper in Section V.

##  II Framework and Problem Analysis

This section consists of two parts. First, the framework for TFE with limited
network data is delineated. Then, the problem of recovering estimation from
limited network data is formulated.

###  II-A Framework for TFE with Limited Data

Figure 1: An illustration of TFE with limited data.

The process of recovering TFE from limited data is illustrated in Fig. 1.
First, we divide the city map evenly into grids. Meanwhile, the limited data
of vehicle mobility information is gathered from vehicular network by traffic
data collector. Instead of removing data from some specific regions, we
acquire the limited data by randomly selecting a small fraction of vehicles as
data sources uniformly distributed across the city, with equal probability of
recruitment for each vehicle. Following the limited data collection, we design
a matching process between the data and the grid-structured city map. Then,
the original estimation is acquired via computing the average speeds of
vehicles in each grid to construct the traffic flow image. Lastly, the
recovery algorithm is applied to generate precise TFE based on the limited
data.

As to the specific scheme of our work, we first obtain the city map consisting
of grids with equal size and number of $H$*$W$. Corresponding to reality, each
grid in the map represents a particular region of the city. We construct the
grid-structured city map as,

$\mathrm{M}=\begin{bmatrix}\mathrm{m}^{1,1}&\cdots&\mathrm{m}^{1,W}\\\
\vdots&\mathrm{m}^{\text{h,w}}&\vdots\\\
\mathrm{m}^{H,1}&\cdots&\mathrm{m}^{H,W}\end{bmatrix},$ (1)

where $\mathrm{m}^{h,w}$ indicates the grid of map which is located at
$(h,w)$, $H$ and $W$ can be viewed as the height and width of the map.

With the construction of grid-structured city map, the vehicular network data
are matched to respective grids based on the GPS coordinates of vehicles. We
first obtain the vehicle mobility dataset $V_{t}$ with sufficient data, which
is collected at time slot $t$ under ideal circumstances. Subsequently, with
the grid matching outcomes of $V_{t}$ and $\mathrm{M}$, the ideal estimation
of traffic flow is calculated as,

$z_{t}^{h,w}=\mathrm{Mean}(v^{h,w}|v^{h,w}\in V_{t}),$ (2)
$\boldsymbol{Z}_{t}=\begin{bmatrix}z_{t}^{1,1}&...&z_{t}^{1,W}\\\
\vdots&z_{t}^{\text{h,w}}&\vdots\\\ z_{t}^{H,1}&...&z_{t}^{H,W}\end{bmatrix},$
(3)

where the mobility information of vehicles in the grid located at $(h,w)$ is
described as $v^{h,w}$, $z_{t}^{h,w}$ is the average speed of vehicles
obtained from sufficient data in the grid located at $(h,w)$ at time slot $t$,
$Z_{t}$ represents the ideal estimation of traffic flow at time slot $t$.

Then, the acquisition of TFE from limited data is conducted as follows. We
denote the limited vehicle mobility information gathered at time slot $t$ as
$V_{t}^{\mathit{{}^{\prime}}}$, which is obtained through random sampling from
the sufficient dataset $V_{t}$, i.e., $V_{t}^{{}^{\prime}}\subseteq V_{t}$.
This sampling approach assigns equal selection probability to each piece of
vehicle mobility information. With a similar calculating and matching process,
the original estimation of traffic flow at time slot $t$ is denoted as,

$x_{t}^{h,w}=\mathrm{Mean}(v^{h,w}|v^{h,w}\in V_{t}^{\mathit{{}^{\prime}}}),$
(4) $\boldsymbol{X}_{t}=\begin{bmatrix}x_{t}^{1,1}&...&x_{t}^{1,W}\\\
\vdots&x_{t}^{\text{h,w}}&\vdots\\\ x_{t}^{H,1}&...&x_{t}^{H,W}\end{bmatrix},$
(5)

where $x_{t}^{h,w}$ is the average speed of vehicles obtained from limited
data in the grid located at $(h,w)$ at time slot $t$, $X_{t}$ represents the
city-wide original estimation of traffic flow at time slot $t$. Note that both
$X_{t}$ and $Z_{t}$ can be perceived as images with a resolution of $H*W$ and
four channels, where the amount of channels is decided by the types of driving
directions defined in our work, i.e., east, south, west, and north.

Finally, we apply the specifically designed deep learning model to recover the
original estimation $X_{t}$. The ideal estimation $Z_{t}$ not only serves as
the ground truth label for the learning of the model, but also is the target
data that the model expects to output. Following the training and inference of
the model, we obtain the recovered TFE from limited data, which is reliable
and cost-effective.

###  II-B Problem of Traffic Flow Estimation

Based on the limited data, we can obtain the original estimation of traffic
flow approximating the ideal one. However, the original estimation exhibits a
considerable degree of inaccuracy and instability regardless of the
similarity. In other words, we found that there exits abundant noise in each
grid of original TFE $\boldsymbol{X}_{t}$, which fails to serve as the
definitive estimation for the applications of urban transportation systems.
Therefore, the task of real-time TFE with limited network data is
conceptualized as a grid-structured data recovery problem.

We set the input of the recovery model to be the historical and current
estimations of traffic flow derived from limited data, represented as
$[\boldsymbol{X}_{t-L+1},...,\boldsymbol{X}_{t-1},\boldsymbol{X}_{t}]$. The
objective is to reform an estimation of current traffic flow with high
accuracy and stability, mirroring the ideal estimation $\boldsymbol{Z}_{t}$
derived from sufficient data. Therefore, the crux of the TFE problem lies in
ascertaining the mapping function $\mathcal{F}(\cdot)$,

$\boldsymbol{Z}_{t}=\mathcal{F}([\boldsymbol{X}_{t-L+1},...,\boldsymbol{X}_{t-1},\boldsymbol{X}_{t}]),$
(6)

which aims to fully explore the spatial-temporal dependencies of traffic flow,
mitigate the noise resulting from the data limitation, and transform the
original TFE derived from limited data into a precise estimation of current
traffic flow.

##  III ST-Mamba for Traffic Flow Recovery

In this section, we present the ST-Mamba designed to recover accurate and
stable TFE from limited network data by exploring and aggregating the features
of traffic flow. The network structures employed for capturing the spatial and
temporal correlations by the CNN and Mamba are discussed. Subsequently, we
provide an overview of the architecture belonging to the ST-Mamba.

###  III-A Spatial Correlation Modeling

In our framework, the estimation of traffic flow can be viewed as an image
with a resolution of $H*W$ and four channels. The spatial correlation of
traffic flow is inherent in the grid structure. Given that the condition of
traffic flow in one region is impacted by those in adjacent regions while the
influence wanes with increasing spatial distance, the aggregation of spatial
information from neighboring regions is imperative for accurate TFE. CNN is a
classic and concise neural network suitable for processing grid-structured
data, the principle of which is convolving input data with learnable filters
to detect local patterns [9]. In recognition of the features of traffic flow,
we employ CNN for capturing local spatial correlation of traffic flow from the
entire traffic flow image.

We first conduct convolutional operations by CNN separately on four channels
of the image, capturing spatial correlation of each driving direction. The
working process of the CNN in our model can be denoted as,

$X_{n}^{\mathit{{}^{\prime}}}=\mathrm{CNN_{e}}(X_{n}),$ (7)

where subscript $n$ ranges from $t-L+1$ to $t$, $X_{n}\in\mathbb{R}^{4\times
H\times W}$ is the original TFE derived from limited data at time slot $n$,
$X_{n}^{\mathit{{}^{\prime}}}\in\mathbb{R}^{4\times H\times W}$ contains the
spatial correlation captured by CNN. Then, we utilize a fully connected (FC)
layer to lift data dimension and integrate spatial correlations of different
driving directions. The FC layer is operated as,

$X_{n}^{\mathit{{}^{\prime\prime}}}=\mathrm{FC_{e}}(X_{n}^{\mathit{{}^{\prime}}}),$
(8)

where $X_{n}^{\mathit{{}^{\prime\prime}}}\in\mathbb{R}^{K\times H\times W}$ is
the high-dimensional feature representation of $X_{n}^{\mathit{{}^{\prime}}}$
in latent space. Subsequently, one reshape layer is constructed to serialize
the format of traffic data. The sequential traffic data of each grid can be
obtained as,

$I^{\mathit{h,w}}=\mathrm{Concat}(x_{n}^{h,w}|n=t-L+1,...t-1,t),$ (9)

where $\mathrm{Concat}(\cdot)$ is the concatenating operation, $x_{n}$ is the
spatial feature representation in one grid of
$X_{n}^{\mathit{{}^{\prime\prime}}}$, $I\in\mathbb{R}^{K\times L}$ is the
sequence comprising spatial correlation of traffic flow in one grid of
$X_{n}^{\mathit{{}^{\prime\prime}}}$.

Figure 2: ST-Mamba framework.

###  III-B Temporal Correlation Modeling

The temporal correlation of traffic flow is hidden in the sequential
structure. The original TFE at time slot $t$, combined with those at preceding
$(L-1)$ time slots, are utilized to construct sequential traffic data for
recovering TFE at current time slot $t$. Mamba is an advanced neural network
especially effective at processing sequential data and modeling temporal
correlation. Integrating the selection mechanism which associates the
parameters of SSM with the input, Mamba is able to focus on relevant inputs
and ignore irrelevant ones, executing linear-time processing of sequences
while maintaining high performance [13]. It is appropriate to employ Mamba for
learning the complex temporal correlation of traffic flow.

We aim to explore temporal correlation in each grid of the original TFE, which
directs us to set sequence $I\in\mathbb{R}^{K\times L}$ as the input of Mamba.
As the pre-foundation of Mamba, we first define the continuous-time SSM to
describe the state representation of the sequence at each time step and
predict the next state based on the input, with $\mathrm{\mathbf{A}}$ as the
evolution parameter and $\mathbf{B}$, $\mathbf{C}$ as the projection
parameters. The continuous-time SSM is formulated as,

$H^{\mathit{{}^{\prime}}}(t)=\mathbf{A}H(t)+\mathbf{B}I(t),$ (10)
$O(t)=\mathbf{C}H(t),$ (11)

where $H(t)$ is the implicit hidden state at time slot $t$, $O(t)$ is the
output of the SSM and is obtained via the projection of the latest latent
state.

Then, we construct Mamba with four matrix format parameters $\mathbf{\Delta}$,
$\mathrm{\mathbf{A}}$, $\mathrm{\mathbf{B}}$, $\mathrm{\mathbf{C}}$, which is
the discrete version of the continuous-time SSM. The transformation method for
discretization in Mamba is zero-order hold, which can be written as,

$\mathbf{\overline{A}}=\mathrm{exp}(\mathbf{\Delta}\mathbf{A}),$ (12)
$\mathbf{\overline{\mathbf{B}}}=(\mathbf{\Delta}\mathbf{A})^{-1}(\mathrm{exp}(\mathbf{\Delta}\mathbf{A})-\mathbf{E})\cdot\mathbf{\Delta}\mathbf{B},$
(13)

where $\text{($\cdot$)}^{-1}$ represents the inverse matrix operation,
$\mathbf{E}$ stands for the identity matrix, $\mathbf{\Delta}$ is the
timescale parameter targeted at transforming the continuous parameters
$\mathbf{A}$, $\mathbf{B}$ to discrete ones $\mathbf{\overline{A}}$,
$\mathbf{\overline{B}}$.

Subsequently, we realize the selection mechanism of Mamba by configuring the
weight matrix $\mathrm{\mathbf{B}}$, $\mathrm{\mathbf{C}}$, $\mathbf{\Delta}$
to be input-dependent, which filters out useless information, compresses input
sequence selectively into the hidden state, and handles contextual information
effectively. Th selection mechanism is formulated as,

$\mathbf{B}=\mathrm{Linear}_{N}(I),$ (14) $\mathbf{C}=\mathrm{Linear}_{N}(I),$
(15)
$\mathbf{\Delta}=\mathrm{softplus}(\mathrm{Broadcast}_{K}(\mathrm{Linear}_{1}(I))+\mathbf{D}),$
(16)

where $\mathbf{\mathrm{Linear}}_{d}$ is a parameterized projection to
dimension $d$, $N$ is the dimension of hidden state, $\mathrm{Broadcast}_{K}$
is targeted at broadcasting one-dimensional data into a $K$-dimensional space,
$\mathbf{D}$ is a constant weight matrix with dimension $K$,
$\mathrm{softplus}$ ensures numerical stability as an activation function.
Note that $\mathbf{A}$ comprises constant parameters for simplicity, while
$\mathbf{\overline{A}}$ is input-dependent via the process of discretization
by $\mathbf{\Delta}$. Meanwhile, although $\mathbf{B}$ and $\mathbf{C}$ are
computed in a similar way, they are independent sets of parameters and
represent different dynamics in Mamba.

Following the discretization of parameters and the implementation of selection
mechanism, we conduct the computation of global convolution to obtain the
output of Mamba as,

$\mathbf{\overline{S}}=(\mathbf{C\overline{B}},\mathbf{C\overline{A}\overline{B}},...,\mathbf{C}\mathbf{\overline{A}}^{L-1}\mathbf{\overline{B}}),$
(17) $O=I\ast\mathbf{\overline{S}},$ (18)

where $\mathbf{\overline{S}}$ is a structured convolutional kernel for
implementing the selective scan algorithm to speed up the model,
$O\in\mathbb{R}^{K\times L}$ is the final output of Mamba which contains
spatial and temporal correlations of traffic flow in each grid.

###  III-C Outcome Decoding

The final outcome is decoded from the output of Mamba. The decoder block
consists of one FC layer and one CNN layer. In the FC layer, we combine the
sequential traffic data $O$ to form an image format with dimension
$\mathbb{R}^{(K\times L)\times H\times W}$. Meanwhile, we reduce data
dimension to aggregate spatial-temporal correlations previously modeled. The
FC layer in decoder block is operated as,

$G=\mathrm{FC_{d}}(O),$ (19)

where $G\in\mathbb{R}^{4\times H\times W}$ can be viewed as the image of TFE
containing spatial and temporal correlations of traffic flow. Then, we conduct
convolutional operations in the CNN layer to decode the outcome as,

$\boldsymbol{Y}=\mathrm{CNN_{d}}(G),$ (20)

where $\boldsymbol{Y}\in\mathbb{R}^{4\times H\times W}$ is the final outcome
of the model and can be perceived as the recovered TFE approximating the ideal
estimation $\boldsymbol{Z}_{t}$.

###  III-D ST-Mamba Model

The main architecture of ST-Mamba model is illustrated in Fig. 2. We
specifically design three parts for ST-Mamba, which are the encoder block, the
Mamba block, and the decoder block, respectively.

In specific, the encoder block consists of one CNN layer and one FC layer in
turn. Given that vehicles driving in different directions have minimal
influence on each other, four CNNs with identical parameters are employed in
the CNN layer to independently perform convolutional operations in each
channel of input, extracting spatial correlation of traffic flow under four
distinct driving directions. Then, we utilize one FC layer to lift the
dimension of input data to explore much more detailed information and exchange
spatial features of different driving directions. For the Mamba block, we
stack one reshape layer before the Mamba layer, serializing the input so that
Mamba can process the limited data of traffic flow effectively. Then, the
generated sequences are fed into the Mamba layer in parallel to explore the
temporal correlation in each gird of the original traffic flow estimations.
For the decoder block which is similar to encoder one, the role of the FC
layer is to aggregate spatial-temporal correlations captured previously and
project the processed data into the image format. Then, convolutional
operations are done by the CNN layer to decode the final outcome while form a
formal symmetry of the model. Finally, ST-Mamba exports the recovered traffic
flow estimation at time slot $\mathrm{\mathit{t}}$ from limited data, with
high accuracy and stability which can meet the demands of traffic
transportation systems.

(a)  Figure 3: The ideal TFE on Friday 5:00 PM. (a)  Figure 4: The original
TFE on Friday 5:00 PM at 10% limitation. (a)  Figure 5: The recovered TFE on
Friday 5:00 PM at 10% limitation.

##  IV Simulation and Results

In this section, we assess our method utilizing the real traffic data from
Beijing. We begin by the experimental settings comprising the dataset,
baselines, and evaluation metrics. Then, we present and analyze the estimation
performance.

###  IV-A Experimental Setting

The simulation utilizes the real traffic data sourced from the Fourth Ring
Road in Beijing. Our dataset includes 6 days in 2012, and the time period for
TFE spans from 7:30 AM to 10:30 PM [6]. The resolution of TFE image is
100$\times$100\. The limited data is obtained by random sampling from
sufficient data, where the limitation degree indicates the proportion of
sampled data. Our approach is compared with the following four baselines:

  * •

Original: The original estimation derived from limited data, and the input of
deep learning models.

  * •

CNN: A convolutional neural network designed to capture spatial correlation by
convolutional operations [9].

  * •

PredRNN: A recurrent neural network proposed to capture spatial-temporal
dependencies for coherent future frame forecasting [14].

  * •

SimVP: A simple yet powerful architecture for video prediction, which makes
the most of CNN and its variants to explore spatial-temporal correlations
[15].

Meanwhile, we select three metrics to evaluate the performance of models. Root
mean squared error (RMSE) assesses the recovery accuracy by penalizing larger
errors more heavily. Improved percentage (IP) demonstrates the recovery
improvement brought by the model. Mean absolute error (MAE) calculates the
average absolute difference between the recovered and ideal estimation.

Figure 6: The CDFs of estimation error on Sunday.

###  IV-B Performance Analysis

TABLE I: The estimation recovery comparison. Limitation | Model | Metric
---|---|---
RMSE |  IP |  MAE
10% | Original |  13.269 |  * |  5.379
CNN |  9.955 |  24.979% |  5.359
PredRNN |  8.067 |  39.207% |  3.965
SimVP |  8.224 |  40.864% |  4.561
ST-Mamba |  $\mathbf{7.504}$ |  $\mathbf{43.449\%}$ |  $\mathbf{3.531}$
20% | Original |  11.072 |  * |  4.073
CNN |  8.786 |  20.638% |  5.124
PredRNN |  6.966 |  37.082% |  3.231
SimVP |  7.049 |  38.415% |  3.680
ST-Mamba |  $\mathbf{6.379}$ |  $\mathbf{42.386\%}$ |  $\mathbf{2.851}$
30% | Original |  9.481 |  * |  3.204
CNN |  7.957 |  16.065% |  4.067
PredRNN |  6.443 |  32.043% |  2.843
SimVP |  6.126 |  37.042% |  3.001
ST-Mamba |  $\mathbf{5.516}$ |  $\mathbf{41.813\%}$ |  $\mathbf{2.245}$
40% | Original |  8.187 |  * |  2.547
CNN |  7.134 |  12.854% |  3.529
PredRNN |  5.967 |  27.118% |  2.615
SimVP |  5.341 |  36.183% |  2.537
ST-Mamba |  $\mathbf{4.857}$ |  $\mathbf{40.670\%}$ |  $\mathbf{1.876}$
50% | Original |  7.056 |  * |  2.015
CNN |  6.631 |  6.008% |  3.077
PredRNN |  5.817 |  17.549% |  2.472
SimVP |  4.703 |  34.664% |  1.989
ST-Mamba |  $\mathbf{4.288}$ |  $\mathbf{39.220\%}$ |  $\mathbf{1.540}$

Table I presents a comparative analysis of recovery performance by the ST-
Mamba against other baselines. Models are evaluated across five degrees of
data limitation from 10% to 50%. ST-Mamba demonstrates superior recovery
capability compared to the other models, consistently attaining the lowest
RMSE and MAE, alongside the highest IP. These results indicate that the
traffic flow estimations recovered by ST-Mamba have the least error. In
addition, considering the significant errors associated with original TFE, the
errors of ST-Mamba are deemed acceptable for practical applications.

As to the recovery performance within space domain, Fig. 3, Fig. 4, and Fig. 5
show the TFE under different conditions. Fig. 3 is the ideal TFE, calculating
the average vehicle speeds from 100% data and clearly demonstrating the
framework of Fourth Ring Road. Fig. 4 is the original TFE at the data
limitation degree of 10%, where considerable number of data missing and errors
can be witnessed. Fig. 5 is the recovered TFE obtained by ST-Mamba and is
similar to the ideal one in terms of both grid color and image structure,
proving excellent recovery performance of ST-Mamba in space domain.

The recovery performance within time domain is illustrated in Fig. 6, where we
display the cumulative distribution function (CDF) of the estimation error per
minute from 7:30 to 22:30 on Sunday at limitation degree of 30%. It can be
observed the CDF of ST-Mamba converges the fastest, which means the errors of
ST-Mamba are concentrated on small values and the error range is relatively
narrow compared to baselines. In other words, ST-Mamba is able to provide
reliable recovery for TFE throughout a day.

In general, ST-Mamba is able to capture spatial-temporal correlations of
traffic flow effectively, providing accurate and stable TFE for the entire
city on various days.

##  V Conclusion

In this paper, we studied the framework for TFE utilizing limited data from
the vehicular network. We analyzed the TFE recovery problem caused by the
limitation of traffic data availability. To address this issue, we designed a
spatial-temporal deep learning model named as ST-Mamba to improve the
estimation accuracy and stability. Comprehensive simulation results based on
the real-world dataset have validated the effectiveness of ST-Mamba in
processing sequential grid-structured data. Our approach provides a cost-
effective solution for TFE in scenarios with limited data. For future work, we
intend to explore the traffic flow prediction based on limited data.

## Acknowledgment

This work is supported in part by the National Natural Science Foundation of
China under Grant 623B2052, 62271244, the Natural Science Fund for
Distinguished Young Scholars of Jiangsu Province under Grant BK20220067.

## References

  * [1] Y. Zhang, Q. Cheng, Y. Liu, and Z. Liu, “Full-Scale Spatio-Temporal Traffic Flow Estimation for City-Wide Networks: A Transfer Learning Based Approach,” _Transportmetrica B: Transport Dynamics_ , vol. 11, no. 1, pp. 869–895, 2023.
  * [2] P. Wang, J. Lai, Z. Huang, Q. Tan, and T. Lin, “Estimating Traffic Flow in Large Road Networks Based on Multi-Source Traffic Data,” _IEEE Transactions on Intelligent Transportation Systems_ , vol. 22, no. 9, pp. 5672–5683, 2020.
  * [3] N. Caceres, L. M. Romero, F. G. Benitez, and J. M. del Castillo, “Traffic Flow Estimation Models Using Cellular Phone Data,” _IEEE Transactions on Intelligent Transportation Systems_ , vol. 13, no. 3, pp. 1430–1441, 2012.
  * [4] A. Abadi, T. Rajabioun, and P. A. Ioannou, “Traffic Flow Prediction for Road Transportation Networks With Limited Traffic Data,” _IEEE Transactions on Intelligent Transportation Systems_ , vol. 16, no. 2, pp. 653–662, 2015.
  * [5] C. Zha, J. Xue, Q. Shen, and H. Zhou, “Low-Cost Traffic Perception for Road Detector Data Estimation: A Deep Learning Approach,” in _the Proceedings of International Conference on Wireless Communications and Signal Processing (WCSP)_. IEEE, 2022, pp. 1–5.
  * [6] J. Xue, Y. Xu, W. Wu, T. Zhang, Q. Shen, H. Zhou, and W. Zhuang, “Sparse Mobile Crowdsensing for Cost-Effective Traffic State Estimation With Spatio-Temporal Transformer Graph Neural Network,” _IEEE Internet of Things Journal_ , vol. 11, no. 9, pp. 16 227–16 242, 2024\.
  * [7] A. Miglani and N. Kumar, “Deep Learning Models for Traffic Flow Prediction in Autonomous Vehicles: A Review, Solutions, and Challenges,” _Vehicular Communications_ , vol. 20, p. 100184, 2019.
  * [8] J. Xue, K. Yu, T. Zhang, H. Zhou, L. Zhao, and X. Shen, “Cooperative Deep Reinforcement Learning Enabled Power Allocation for Packet Duplication URLLC in Multi-Connectivity Vehicular Networks,” _IEEE Transactions on Mobile Computing_ , vol. 23, no. 8, pp. 8143–8157, 2024\.
  * [9] L. Alzubaidi, J. Zhang, A. J. Humaidi, A. Al-Dujaili, Y. Duan, O. Al-Shamma, J. Santamaría, M. A. Fadhel, M. Al-Amidie, and L. Farhan, “Review of Deep Learning: Concepts, CNN Architectures, Challenges, Applications, Future Directions,” _Journal of Big Data_ , vol. 8, pp. 1–74, 2021.
  * [10] A. Sherstinsky, “Fundamentals of Recurrent Neural Network (rnn) and Long Short-Term Memory (LSTM) Network,” _Physica D: Nonlinear Phenomena_ , vol. 404, p. 132306, 2020.
  * [11] L. Li, H. Wang, W. Zhang, and A. Coster, “STG-Mamba: Spatial-Temporal Graph Learning via Selective State Space Model,” _arXiv preprint arXiv:2403.12418_ , 2024.
  * [12] Y. Liu, H. Zheng, X. Feng, and Z. Chen, “Short-Term Traffic Flow Prediction with Conv-LSTM,” in _the Proceedings of International Conference on Wireless Communications and Signal Processing (WCSP)_ , 2017, pp. 1–6.
  * [13] A. Gu and T. Dao, “Mamba: Linear-Time Sequence Modeling with Selective State Spaces,” _arXiv preprint arXiv:2312.00752_ , 2023.
  * [14] Y. Wang, H. Wu, J. Zhang, Z. Gao, J. Wang, S. Y. Philip, and M. Long, “PredRNN: A Recurrent Neural Network for Spatiotemporal Predictive Learning,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , vol. 45, no. 2, pp. 2208–2225, 2022.
  * [15] Z. Gao, C. Tan, L. Wu, and S. Z. Li, “SimVP: Simpler yet Better Video Prediction,” in _the Proceedings of the IEEE/CVF conference on computer vision and pattern recognition (CVPR)_ , 2022, pp. 3170–3180.


# Multi-Contact Force-Sensing Guitar for Training and Therapy

Zhiyi Ren1, Chun-Cheng Hsu1, Can Kocabalkanli1, Khanh Nguyen1, Iulian I.
Iordachita1,
Serap Bastepe-Gray23, Nathan Scott1 3Department of Neurology, School of
Medicine, Johns Hopkins University, Baltimore, MD, USA
Email: zren5, chsu40, ckocaba1, knguye44<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>1Department of Mechanical Engineering,
Johns Hopkins University, Baltimore, MD, USA  2The Peabody Institute, Johns
Hopkins University, Baltimore, MD, USA

###### Abstract

Hand injuries from repetitive high-strain and physical overload can hamper or
even end a musician’s career. To help musicians develop safer playing habits,
we developed a multiple-contact force-sensing array that can substitute as a
guitar fretboard. The system consists of 72 individual force sensing modules,
each containing a flexure and a photointerrupter that measures the
corresponding deflection when forces are applied. The system is capable of
measuring forces between 0-25 N applied anywhere within the first 12 frets at
a rate of 20 Hz with an average accuracy of $\scriptstyle\pm$0.4 N and a
resolution of 0.1 N. Accompanied with a GUI, the resulting prototype was
received positively as a useful tool for learning and injury prevention by
novice and expert musicians.

###### Keywords:

Force sensing, photointerrupter, sensing array, flexure deflection, injury
prevention

##  I Introduction

Many musicians never fully recover to play at their highest level after
suffering from injuries. For guitarists, injuries are common even in experts
due to the repetitive overloading of finger joints developed through years of
harmful playing habits [1, 2]. By measuring the dynamic forces applied by the
players on the frets and strings of the guitar, we aim (1) to warn musicians
in real time during practice when their grip is too strong, (2) and to enforce
correct force patterns.

In a preliminary experiment involving one of the senior authors, two 6-axis
ATI Nano 17 force/torque sensors were mounted on a guitar. Despite valuable
insight, the system modified the original shape of the guitar significantly
and only measures the gross force on the fretboard, thus unable to distinguish
individual finger forces. Similar work in force sensing on musical instruments
suffers the same issues [3, 4]. In our work, we address both shortcomings by
developing a modular sensing array that can be packed into a customized, 4mm
thick fretboard (Fig. 1). It can measure forces at 72 intersections of frets
and strings, and can be mounted onto a classical guitar easily with little
change in appearance.

Typical sensing elements used in force sensing arrays include piezoresistive
sensors [5], capacitive sensors [6], and resistive composites like Velostat
[7]. Both piezoresistive sensors and Velostat are known for their inherent
drift and hysteresis. Capacitive sensors require complicated circuitry and
customization processes [8, 9]. Instead, we adopted optical proximity sensors
as the sensing elements, inspired by [10]. The sensors are low-cost and
lightweight, and require minimal space for operation. The nature of the
technology also eliminates the possibility of drift or hysteresis. In our
work, an array of photointerrupters (Sharp GP2S60) measures the distance to
flexures that deflect proportionally to the force applied. While light
isolation is necessary to mitigate interference from adjacent modules,
specific design considerations detailed in Section II were developed to
resolve this issue. Though 3-axis force measurement would be ideal and
beneficial for our purposes (to measure forces from vibratos, for instance),
we focus only on the vertical component of the force in this preliminary work.
The preliminary experiment using ATI Nano 17 sensors reveals that a sensing
range of 0-25N is sufficient. After calibrating and validating to the full
sensing range, we achieved an average error of <0.4N RMSE (root mean squared
error) and a worst error of <5% FSO (full scale output) at a resolution of
0.1N on most sensing modules.

Figure 1: (a) GUI showing forces at each string-fret intersection during user
testing. Tiles turn red when user presses above a specified force threshold.
(b) Fretboard terminology. Also note the contortion of the hand.

##  II Implementation

###  II-A Sensing module characteristics

The system is comprised of 72 sensing modules, each containing a pair of
flexure and photointerrupter. Each photointerrupter has an emitter and a
detector: infrared light is emitted by the emitter, reflected off from a fixed
surface, and received by the detector. The output is proportional to the
amount of light received, which corresponds to the distance to the surface
when the environment is dark. As the force is applied onto the flexure, the
distance from the flexure to photointerrrupter changes.

###  II-B Flexure design

The purpose of the flexures is twofold: to provide a measurable deflection and
a fixed boundary condition for which the string can vibrate on to create the
standing wave we hear as musical notes. Through our finite element analysis in
ANSYS [11], we designed the thickness of the flexure such that it elastically
deflects 0.2mm under 25N of force. To create the fixed boundary condition, the
flexures are designed with the exact same crown profile as a regular fret,
with the additional ability to independently deflect under each of the 6
strings (Fig. 2). A compromise has been made to replace the bottom 7 frets
(less often used) with non-functional ones in order to house electronic
components.

Figure 2: (a) A single flexure unit with the profile of a regular fret. (b)
FEA analysis under static load. (c) Deflection process when force is applied.

###  II-C Electronics design

All electronics components are mounted on a single piece of PCB spanning the
whole fretboard. A chain of shift registers (Texas Instruments SN74HC164)
activate each of the 72 photointerrupters in sequence. The output signal of
each photointerrupter, passing through a unity gain amplifier, is then fed
into a differential amplifier (Texas Instruments LM324). There are 6 output
lines corresponding to the 6 strings (Fig 3). A digital-to-analog converter
(Texas Instruments TLV5638) provides the reference voltage for each
differential amplifier. The whole system is controlled by a microprocessor
(Texas Instruments MSP430F5342), which has an internal, 12-bit analog-to-
digital converter that reads the sensor signal.

The system is powered through a USB connection through a Micro-USB port
installed at the bottom of the guitar. By connecting the fretboard to a
computer, the user can visualize real-time force measurements at all locations
at a rate of 20Hz and save them in text files with a custom GUI made in
Processing language [12].

Figure 3: Simplified diagram of the circuits. The red dotted lines in the
photointerrupter array indicate light isolation between frets. Figure 4: Top
view of the fretboard and the PCB underneath. The top 12 frets measure forces
with photointerrupters. The lower 7 dummy frets house the rest of the
electronics.

###  II-D Assembly

The aluminum flexures are glued onto the PCB surface with Araldite glue
(Huntsman, The Woodlands, TX) such that each photointerrupter sits directly
under each flexure. The PCB fits into an aluminum case that houses the
flexures and covers the remaining area of the PCB. The case is screwed onto
the guitar neck with inserts, and a wooden veneer is glued onto the outermost
surface.

Figure 5: Exploded views of the sensing fretboard. The whole system is
contained in a custom, 4mm thick fretboard, which can be mounted on a
classical guitar with screws.

###  II-E Light isolation

The mechanical and electronic designs ensure that there is no light
interference between photointerrupters, which is critical for the performance
of a optical sensing array. The aluminum case fills the space between the
frets thus no light can transmit among them, as shown with the red dotted
lines in Fig. 3. This allows the 12 photointerrupters on the same string to
share a single output bus. Since each of the 6 photointerrupters within the
same fret are on separate output lines and are activated in sequence, there is
no light interference within the fret (among the 6 adjacent sensors) either.
In addition, the veneer on top of the case also protects the photointerrupters
from any external light source.

##  III Experiments and Results

###  III-A Testing setup

We developed a testing setup that can move freely along both the string and
fret directions to calibrate and validate each of the 72 sensing modules
efficiently. By precisely displacing a micron stage we apply varying forces
onto the the crown on the sensing unit (Fig. 2c). A Sparkfun TAL220B load cell
($\scriptstyle\pm$0.025 N accuracy), calibrated with known weights, is mounted
at the tip of the micron stage to provide the ground truth reading against the
output of the sensing module.

Figure 6: Testing rig used to apply force across all 72 modules. Figure 7:
Calibration results of five randomly chosen modules. The coefficient of
determination for each is higher than 0.99. Figure 8:  Histogram of average
error in RMSE (left) and histogram of worst error for all modules (right). Red
dotted lines show the desired 0.4 N average error in RMSE and 5% worst error
in FSO.

###  III-B Force-to-sensor calibration

We loaded and unloaded each of the 72 sensing modules between 0 and 25 N for 2
trials. Results for all modules exhibit a high linearity between sensor output
and true force measured by the load cell (Fig. 7). Since stiffness varies with
size of the flexures as fret width also varies, the calibration slopes of the
five modules shows some differences.

###  III-C Validation

To test the accuracy of the sensing modules, we apply random level of forces
to the sensors using the micron stage and calculate the errors based on
individual calibration curves. Fig. 8 shows the histograms for the average and
worst error for all 72 sensing modules. We find that 81% of the modules show
an average error of <0.4 N RMSE, and 90% of the modules show a worst error of
<5% FSO. The results satisfy the high accuracy desired for our purposes.

We believe that variances among sensing modules may be due to local
temperature differences. To combat this issue, a temperature compensation
system is installed within the current prototype: beyond the 6
photointerrupters at each fret, an extra one is placed in the middle and
measures against a fixed surface. Assuming that temperature affects all
photointerrupters in a fret the same way, the differential between each of the
6 photointerrupters with the fixed one should remain constant independent of
the temperature. This compensation feature was tested independently but was
not used during the user tests.

##  IV Conclusion and Future Work

This paper presents a novel compact, multi-contact force-sensing fretboard for
a classical guitar. The 4mm thick design addresses the previous prototype’s
inability to distinguish individual forces on different fingers and the
changes in the guitar’s shape and weight. It was evaluated using a custom-
designed test rig to demonstrate the desired accuracy and linearity. The
prototype weighs and looks virtually identical to a regular fretboard once
mounted. The functionality and feel of our system have also been evaluated by
musicians from the Peabody Institute of the Johns Hopkins University, using a
Likert scale questionnaire with a section for general comments. Users found
the force sensing capabilities and real-time visual feedback constructive to
their playing, despite noticing slight differences in sound quality.

In the future more data will be collected and analyzed through user testing.
Other directions include testing and evaluating the temperature compensation
system, improving the update rate of the system and reducing electronic noise.

## Acknowledgment

We thank Melissa Hullman and Dr. Jing Xu for their valuable insights, luthier
Garrett Lee for his help on building the final guitar, and Drs. Jae Kun Shim
and Hyun Joon Kwon (Department of Kinesiology, University of Maryland) for
earlier prototype development, preliminary data collection and analysis. The
project was funded by private philanthropy through the Peabody Conservatory.

## References

  * [1] I. Winspur and C. Wynn Parry, “The musician’s hand,” _Journal of Hand Surgery_ , vol. 22, no. 4, pp. 433–440, 1997.
  * [2] J. L. Rigg, R. Marrinan, and M. A. Thomas, “Playing-related injury in guitarists playing popular music,” _Medical Problems of Performing Artists_ , vol. 18, no. 4, pp. 150–152, 2003.
  * [3] M. Hori, M. Hosono, H. Takahashi, K. Matsumoto, and I. Shimoyama, “3-axis fingertip force during playing the string instrument,” in _2013 Transducers & Eurosensors XXVII: The 17th International Conference on Solid-State Sensors, Actuators and Microsystems (Transducers & Eurosensors XXVII)_. IEEE, 2013, pp. 2745–2748.
  * [4] H. Kinoshita and S. Obata, “Left hand finger force in violin playing: Tempo, loudness, and finger differences,” _The Journal of the Acoustical Society of America_ , vol. 126, no. 1, pp. 388–395, 2009.
  * [5] T. V. Papakostas, J. Lima, and M. Lowe, “A large area force sensor for smart skin applications,” in _Sensors, 2002 IEEE_ , vol. 2. IEEE, 2002, pp. 1620–1624.
  * [6] T. Hoshi and H. Shinoda, “A sensitive skin based on touch-area-evaluating tactile elements,” in _2006 14th Symposium on Haptic Interfaces for Virtual Environment and Teleoperator Systems_. IEEE, 2005, pp. 89–94.
  * [7] D. Giovanelli and E. Farella, “Force sensing resistor and evaluation of technology for wearable body pressure sensing,” _Journal of Sensors_ , vol. 2016, 2016.
  * [8] M. R. Cutkosky, R. D. Howe, and W. R. Provancher, “Force and tactile sensors,” _Springer Handbook of Robotics_ , pp. 455–476, 2008.
  * [9] X. A. Wu, J. V. Ulmen, and M. R. Cutkosky, “Capacitive force/torque sensor,” Mar. 15 2018, uS Patent App. 15/703,122.
  * [10] J. D. Goldberg and R. S. Fearing, “Force sensing shell using a planar sensor for miniature legged robots,” in _2015 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2015, pp. 1494–1500.
  * [11] Engineering simulation & 3d design software ansys. [Online]. Available: https://www.ansys.com/
  * [12] Processing. [Online]. Available: http://processing.org