# Probing Hidden Leptonic Scalar Portals using the NA64 Experiment at CERN

A. Ponten  0009-0008-2463-2930 Uppsala University, Department of Physics and
Astronomy, High Energy Physics, Ångströmlaboratoriet, Lägerhyddsvägen 1, 752
36 Uppsala, Sweden  H. Sieber  0000-0003-1476-4258<EMAIL_ADDRESS>ETH Zürich, Institute for Particle Physics and Astrophysics, CH-8093 Zürich,
Switzerland  B. Banto Oberhauser  0009-0006-4795-1008 ETH Zürich, Institute
for Particle Physics and Astrophysics, CH-8093 Zürich, Switzerland  P.
Crivelli  0000-0001-5430-9394<EMAIL_ADDRESS>ETH Zürich, Institute for
Particle Physics and Astrophysics, CH-8093 Zürich, Switzerland  D.
Kirpichnikov  0000-0002-7177-077X Institute for Nuclear Research of the
Russian Academy of Sciences, 117312 Moscow, Russia  S. N. Gninenko
0000-0001-6495-7619 Institute for Nuclear Research of the Russian Academy of
Sciences, 117312 Moscow, Russia  Center for Theoretical and Experimental
Particle Physics, Facultad de Ciencias Exactas, Universidad Andres Bello,
Fernandez Concha 700, Santiago, Chile  M. Hösgen  Universität Bonn, Helmholtz-
Institut für Strahlen-und Kernphysik, 53115 Bonn, Germany  L. Molina Bueno
0000-0001-9720-9764<EMAIL_ADDRESS>Instituto de Fisica Corpuscular
(CSIC/UV), Carrer del Catedratic Jose Beltran Martinez, 2, 46980 Paterna,
Valencia, Spain  M. Mongillo  0009-0000-7331-4076 ETH Zürich, Institute for
Particle Physics and Astrophysics, CH-8093 Zürich, Switzerland  A. Zhevlakov
0000-0002-7775-5917 Bogoliubov Laboratory of Theoretical Physics, JINR, 141980
Dubna, Russia  Matrosov Institute for System Dynamics and Control Theory SB
RAS,
Lermontov str., 134, 664033, Irkutsk, Russia

###### Abstract

In this study, we demonstrate the potential of the NA64 experiment at CERN SPS
to search for New Physics processes involving $e\rightarrow\mu$ transitions
after the collision of 100 GeV electrons with target nuclei. A new Dark Sector
leptonic portal in which a scalar boson $\varphi$ could be produced in the
lepton-flavor-changing bremsstrahlung-like reaction, $eN\rightarrow\mu
N\varphi$, is used as benchmark process. In this work, we develop a realistic
Monte Carlo simulation of the NA64 experimental setup implementing the
differential and total production cross-section computed at exact tree-level
and applying the Weiszäcker-Williams phase space approximation. Using this
framework, we investigate the main background sources and calculate the
expected sensitivity of the experiment. The results indicate that with minor
setup optimization, NA64 can probe a large fraction of the available parameter
space compatible with the muon $g-2$ anomaly and the Dark Matter relic
predictions in the context of a new Dark Sector leptonic portal with $10^{11}$
EOT. This result paves the way to the exploration of lepton-flavour-changing
transitions in NA64.

##  I Introduction

The known shreds of evidence for the existence of Dark Matter (DM) are based
on its gravitational behaviour inferred from astronomical, astrophysical, and
cosmological observations. Those notably include galaxy rotational curves [1,
2, 3, 4], gravitational lensing [5, 6, 7], Cosmic Microwave Background (CMB)
anisotropies [8, 9], or large-scale structure (LSS) formation [10, 11, 12],
suggesting that the DM relic density abundance, $\Omega_{\text{DM}}\simeq
0.27$, exceeds that of atoms by a factor
$\Omega_{\text{DM}}/\Omega_{\text{baryon}}\simeq 5$.
Nonetheless, fundamental questions remain open, such as the microscopic nature
of DM or the origin of its observed relic density. Because of the absence of
sizeable interactions with Standard Model (SM) particles, a popular framework,
the so-called Dark Sector (DS), suggests that DM particles are neutral under
SM forces but charged under a new force thus interacting with SM particles
through mediators [13, 14, 15, 16, 17].
Different types of particles can intertwin these interactions opening a broad
parameter space and motivating new and complementary experimental avenues [18,
19]. Portal interactions involving lepton-flavor-changing (LFC) transitions
mediated by a scalar particle could explain the origin of Dark Matter and at
the same time accommodate the long-standing muon $g-2$ puzzle [20]. In this
case, a new scalar boson can act as a mediator between the two sectors
carrying dark and leptonic quantum numbers. As a result, non-diagonal flavour
interactions are allowed involving $e\leftrightarrow\mu$ transitions. The
effective Lagrangian governing the interaction is

$L_{\text{int}}=-h\bar{e}P_{R}\mu\varphi+h.c.,$ (1)

with $P_{R}=1/2(1+\gamma_{5})$, is the right-handed chiral projection
operator, and $h$ is the coupling strength. Within this model, it is assumed
that the new boson has a mass above the muon mass and would predominantly
decay invisibly avoiding existing experimental constraints from muon decay
experiments [21, 20, 22]. Given the thermally averaged annihilation cross-
section, $\langle\sigma v\rangle\simeq\mathcal{O}(1\ \text{pb})$, the DM relic
abundance associated with fermionic DS particles, $\psi_{1}$ and $\psi_{2}$,
is given through [20]

$\frac{h^{2}g^{2}_{\varphi\psi_{1}\psi_{2}}m_{1}^{2}}{4\pi\left\\{m_{\varphi}^{2}-\left(m_{1}+m_{2}\right)^{2}\right\\}^{2}}=\mathcal{O}\left(10\,\text{pb}\right),$
(2)

where $m_{i}$ is the mass of $\psi_{i}$ fermion and $g_{\phi\psi_{1}\psi_{2}}$
is the coupling strength of the scalar boson $\varphi$ to the DS particles.
Besides the possibility of probing the particle nature of DM, such a model
conveniently provides a solution to the muon anomalous magnetic moment $g-2$,
which problem is associated to the discrepancy between the experimental
measured value [23] and the prediction from the Standard Model (SM) [24, 25,
26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43]. In
this case, the contribution to the muon $g-2$ from a hypothetical scalar is:

$\Delta
a_{\mu}^{\varphi}=\frac{h^{2}}{16\pi^{2}}\frac{m_{\mu}^{2}}{m_{\varphi}^{2}}\int_{0}^{1}dx\frac{x^{2}\left(1-x\right)}{\left(1-x\right)\left(1-\frac{m_{\mu}^{2}}{m_{\varphi}^{2}}x\right)+\frac{m_{e}^{2}}{m_{\varphi}^{2}}x}.$
(3)

where $m_{e}$ and $m_{\mu}$ are the masses of electron and muon respectively.
The potential of fixed target experiments to search for such a boson was
described in [20]. In this work, we focus on implementing such a model in a
realistic GEANT4-based [44, 45] framework to derive the sensitivity of the
NA64 experiment [46, 47]. In addition, in this work, the lepton flavour
changing process $e(\mu)+N\rightarrow\mu(e)+N+\varphi$ is described through
its differential and total production cross-sections. Those later are derived
in section II both at exact tree-level (ETL) and in the Weiszäcker-Williams
(WW) phase space approximation. The accuracies of the two approaches are
compared through numerical studies in section II.3. The underlying LFC physics
is implemented in a fully realistic Monte Carlo (MC) simulation to study the
potential reach of the NA64 experiment, as presented in section III. We
discuss in detail the signal yield estimate through the GEANT4-based DMG4 [48,
49] package for the simulation of Dark Matter. Besides, we estimate the main
background sources for a statistics of $10^{11}$ electrons on target (EOT).
The final projected sensitivity for such a model is evaluated in section IV.
Finally, our results are summarised in section V.

##  II Cross-sections calculations

In the following, the production cross-section is derived both at exact tree-
level and using the Weiszäcker-Williams equivalent photon flux phase space
approximation. We follow closely the notations of [50, 51, 52, 53]. For the
sake of generality, the cross-sections are derived for generic incoming and
outgoing leptons, denoted respectively by $l_{i}$ and $l_{f}$, with $l=e,\
\mu$.

###  II.1 The exact tree-level

The $\varphi$ boson production associated with the LFC bremsstrahlung-like
reaction is defined by the kinematics

$l_{i}(p)+N(P_{i})\rightarrow l_{f}(p^{\prime})+N(P_{f})+\varphi(k),$ (4)

where we defined $p=(E_{i},\mathbf{p})$ and
$p^{\prime}=(E^{\prime},\mathbf{p}^{\prime})$ to be respectively the incoming
and outgoing four-momentum of the initial and final state lepton,
$k=(E_{\varphi},\mathbf{k})$ is the $\varphi$ boson momentum, and $M$ the mass
of the target (N). The nucleus has initial and final four-momenta
$P_{i}=(M,\mathbf{0})$ and $P_{f}=(P_{f}^{0},\mathbf{P}_{f})$, and is treated
as a spineless boson to which the associated photon vertex rule reads

$ie\mathcal{P}^{\mu}F(t)=ie(P_{i}^{\mu}+P_{f}^{\mu})F(t),$ (5)

where $F^{2}(t)\equiv G_{2}^{\text{el}}(t)$ is the squared elastic form factor
[51, 52]

$F^{2}(t)=Z^{2}\left(\frac{a^{2}t}{1+a^{2}t}\right)^{2}\left(\frac{1}{1+t/d}\right)^{2},$
(6)

with $a=111Z^{-1/3}/m_{e}$ and $d=0.164A^{-2/3}\ \text{GeV}^{-2}$ and
$P_{\mu}=P_{i\mu}+P_{f\mu}$. For the typical lead target properties of NA64
with $Z=82$ and $A=207$, the inelastic term contributing to Eq. (6) can be
neglected [52]. The relevant matrix element to the $2\rightarrow 3$ LFC
process is then

$\begin{split}i\mathcal{M}_{\varphi}^{2\rightarrow
3}&=ihe^{2}\frac{F(q^{2})}{q^{2}}\mathcal{P}_{\mu}\bar{u}_{f}^{(s^{\prime})}(p^{\prime})\bigg{\\{}P_{R}\frac{(\not{p^{\prime}}+\not{k})+m_{i}}{\tilde{s}}\gamma^{\mu}\\\
&+\,\gamma^{\mu}\frac{(\not{p}-\not{k})+m_{f}}{\tilde{u}}P_{R}\bigg{\\}}u_{i}^{(s)}(p),\end{split}$
(7)

where $\bar{u}_{f}$ and $u_{i}$ are the spinors of the initial and final state
leptons with spins $s,s^{\prime}$, $P_{R}$ the right-handed projection
operator and $q=P_{i}-P_{f}$ the momentum transfer to the nucleus. We also
define the following Mandelstam variables

$\begin{split}\@add@centering\centering\tilde{s}&=(p^{\prime}+k)^{2}-m_{i}^{2},\\\
\tilde{u}&=(p-k)^{2}-m_{f}^{2},\\\ t_{2}&=(p^{\prime}-p)^{2}\\\ t&=-q^{2},\\\
m_{\varphi}^{2}&=\tilde{s}+t_{2}+\tilde{u}+t.\end{split}$ (8)

where $m_{i}$ and $m_{f}$ are the leptons masses for either combinations of
$(m_{e},m_{\mu})$. The average squared amplitude of the process is calculated
using the FeynCalc package [54, 55, 56] from the Wolfram Mathematica program
[57], such that from Eq. (7)

$\overline{|\mathcal{M}_{\varphi}^{2\rightarrow
3}|^{2}}=e^{4}h^{2}\frac{F^{2}(t)}{t^{2}}|\mathcal{A}_{\varphi}^{2\rightarrow
3}|^{2},$ (9)

where

$\begin{split}|\mathcal{A}_{\varphi}^{2\rightarrow
3}|^{2}&=\frac{(\tilde{s}+\tilde{u})^{2}}{2\tilde{s}\tilde{u}}\mathcal{P}^{2}-\frac{2t}{\tilde{s}\tilde{u}}(\mathcal{P}\cdot
k)^{2}\\\
+\frac{(\tilde{s}+\tilde{u})^{2}}{2\tilde{s}^{2}\tilde{u}^{2}}&\Delta{m^{2}}\bigg{\\{}\mathcal{P}^{2}t-4\bigg{(}\frac{\tilde{u}(\mathcal{P}\cdot
p)+\tilde{s}(\mathcal{P}\cdot
p^{\prime})}{\tilde{s}+\tilde{u}}\bigg{)}^{2}\bigg{\\}}.\end{split}$ (10)

with $\Delta m^{2}\equiv m_{\varphi}^{2}-m_{i}^{2}-m_{f}^{2}$ and

$\begin{split}\@add@centering\centering\mathcal{P}^{2}&=4M^{2}+t,\\\
\mathcal{P}\cdot p&=2ME_{i}-(\tilde{s}+t)/2,\\\ \mathcal{P}\cdot
p^{\prime}&=2M(E_{i}-E_{\varphi})+(\tilde{u}+t)/2,\\\ \mathcal{P}\cdot
k&=\mathcal{P}\cdot p-\mathcal{P}\cdot p^{\prime}.\end{split}$ (11)

The double-differential cross-section can then be expressed following [51, 58]

$\begin{split}\left(\frac{d^{2}\sigma_{2\rightarrow
3}}{dxd\cos{\theta_{k}}}\right)_{\text{ETL}}=\frac{\alpha^{2}h^{2}}{4\pi}\frac{E_{i}|\textbf{k}|}{|\textbf{p}|V}&\int_{t_{\text{min}}}^{t_{\text{max}}}\,dt\,\frac{F^{2}(t)}{t^{2}}\\\
&\times\int_{0}^{2\pi}\,\frac{d\phi_{q}}{2\pi}\,\frac{|A_{\varphi}^{2\rightarrow
3}|^{2}}{8M^{2}}.\end{split}$ (12)

with $\alpha$ the fine-structure constant, $t_{\text{min}}$ and
$t_{\text{max}}$ the minimum and maximum momentum transfer respectively as in
[51] and $x=E_{\varphi}/E_{i}$ the fractional energy transferred to $\phi$
dark scalar. The integration is performed in the reference frame where
$\mathbf{V}=\mathbf{p}-\mathbf{k}$ is parallel to $\hat{z}$ axis and
$\mathbf{k}$ is in the $\widehat{xz}$ plane. As such, $\phi_{\mathbf{q}}$ is
the azimuthal angle of $\mathbf{q}$ defined in the polar coordinates system of
the aforementioned reference frame (see [51] for more details). The
$\theta_{k}$ is the radiation angle of the produced $\varphi$ scalar boson
such that

$V=|\textbf{p}-\textbf{k}|=\sqrt{\textbf{p}^{2}+\textbf{k}^{2}-2|\textbf{p}||\textbf{k}|\cos{\theta_{k}}}.$
(13)

###  II.2 The Weiszäcker-Williams approach

In this section, the production cross-sections are derived assuming that the
incoming lepton energy, $E_{i}$, is much larger than the mass of its lepton,
$m_{i}$, and of the scalar boson $m_{\varphi}$. As such we use the Weiszäcker-
Williams approach and we closely follow the procedure of [52, 53]. Therefore,
we can define the equivalent photon flux from nuclear as

$\chi^{\text{WW}}=\int_{t_{\text{min}}}^{t_{\text{max}}}\,dt\
\frac{t-t_{\text{min}}}{t^{2}}F(t)^{2},$ (14)

where $t_{min}\simeq U^{2}/4E_{i}^{2}(1-x)^{2}$ and
$t_{max}=m_{\varphi}^{2}+m_{i}^{2}$. The definition of the function $U$ is
given below in Eq. (18). It is worth noting that Eq. (14) can be expressed
analytically as presented in [52]

$\chi=Z^{2}\left[\tilde{\chi}(t_{\text{max}})-\tilde{\chi}(t_{\text{min}})\right].$
(15)

where we define $t_{d}=d$ and $t_{a}=1/a^{2}$ and express the antiderivative
as

$\begin{split}&\tilde{\chi}(t)=\frac{t_{d}^{2}}{(t_{a}-t_{d})^{3}}\bigg{[}\frac{(t_{a}-t_{d})(t_{a}+t_{min})}{t+t_{a}}\\\
&+\frac{(t_{a}-t_{d})(t_{d}+t_{min})}{t+t_{d}}+(t_{a}+t_{d}+2t_{min})\log{\frac{t+t_{d}}{t+t_{a}}}\bigg{]}.\\\
&\\\ \end{split}$ (16)

Thus, the double differential cross-section in terms of $(x,\ \theta_{k})$ can
be expressed as [52, 53]

$\left(\frac{d^{2}\sigma_{2\rightarrow
3}}{dxd\cos{\theta_{k}}}\right)_{\text{WW}}=\frac{\alpha\chi^{\text{WW}}}{\pi(1-x)}E_{i}|\textbf{k}|\frac{d\sigma_{2\rightarrow
2}}{d(p\cdot k)}\biggr{\rvert}_{t=t_{\text{min}}},$ (17)

with it being evaluated at $t=t_{\text{min}}$. In this approximation, we can
rewrite the Mandelstam variables in the form

$\begin{split}\@add@centering\centering\tilde{s}&\simeq-\frac{\tilde{u}}{(1-x)}\\\
t_{2}&\simeq\frac{x\tilde{u}}{(1-x)}+m_{\varphi}^{2},\\\ U&=-\tilde{u},\\\
\tilde{u}&\simeq-
xE_{i}^{2}\theta_{k}^{2}-m_{\varphi}^{2}\frac{1-x}{x}+m_{i}^{2}(1-x)-m_{f}^{2},\\\
m_{\varphi}^{2}&\simeq\tilde{s}+t_{2}+\tilde{u},\end{split}$ (18)

where we used the definition of $\tilde{u}$ from Eq. (8) and only kept terms
up to order $\mathcal{O}(\theta_{k}^{2})$, $\mathcal{O}(m_{i}^{2}/E^{\prime
2})$ and $\mathcal{O}(m_{i}^{2}/x^{2}E_{i}^{2})$.
The differential cross-section on the right-hand side of Eq. (17) for the
$2\rightarrow 2$ process can be expressed as [59]

$\frac{d\sigma_{2\rightarrow 2}}{d(p\cdot k)}=2\frac{d\sigma_{2\rightarrow
2}}{dt_{2}}=\frac{\alpha
h^{2}}{2\tilde{s}^{2}}|\mathcal{A}_{\varphi}^{2\rightarrow 2}|^{2}.$ (19)

In particular, the spin-averaged amplitude squared,
$|\mathcal{A}_{\varphi}^{2\rightarrow 2}|^{2}$, for $l_{i}+\gamma\rightarrow
l_{f}+a$ of the process reads

$\begin{split}|\mathcal{A}_{\varphi}^{2\rightarrow
2}|^{2}&=-\frac{\left(\tilde{s}+\tilde{u}\right)^{2}}{2\tilde{s}\tilde{u}}\\\
&+\Delta
m^{2}\left\\{\frac{(\tilde{s}+\tilde{u})(\tilde{s}m_{f}^{2}+\tilde{u}m_{i}^{2})}{\tilde{s}^{2}\tilde{u}^{2}}-\frac{t_{2}}{\tilde{s}\tilde{u}}\right\\}.\end{split}$
(20)

Evaluated at $t=t_{\text{min}}$, the above amplitude can be written as

$\begin{split}&|\mathcal{A}_{\varphi}^{2\rightarrow
2}|^{2}_{t_{\text{min}}}=\frac{x^{2}}{2(1-x)}\\\ &+\Delta
m^{2}\cdot\frac{\tilde{u}x+m_{\varphi}^{2}(1-x)+x\left[m_{f}^{2}+m_{i}^{2}(x-1)\right]}{\tilde{u}^{2}}.\\\
\end{split}$ (21)

Plugging Eqs. (21) and (19) into Eq. (17) gives us the following result

$\begin{split}\left(\frac{d^{2}\sigma_{2\rightarrow
3}}{dxd\cos{\theta_{k}}}\right)_{\text{WW}}&=\frac{\alpha^{2}h^{2}}{2\pi}\chi^{\text{WW}}E_{i}^{2}x\beta_{k}\frac{(1-x)}{\tilde{u}^{2}}|\mathcal{A}_{\varphi}^{2\rightarrow
2}|^{2}_{t_{\text{min}}},\end{split}$ (22)

where $\beta_{k}=\sqrt{1-(m_{\phi}^{2}/(x^{2}E_{\varphi}})$. Purely in terms
of the fractional energy $x$, one writes

$\begin{split}\left(\frac{d\sigma}{dx}\right)_{\text{WW}}&=\frac{h^{2}\alpha^{2}}{2\pi}E_{i}^{2}x\beta_{k}(1-x)\\\
&\times\int_{0}^{\theta_{\text{max}}}\,d\cos\theta_{k}\
\frac{|\mathcal{A}_{\varphi}^{2\rightarrow
2}|^{2}_{t_{\text{min}}}}{\tilde{u}^{2}}\chi^{\text{WW}}.\\\ \end{split}$ (23)

The value maximum angle $\theta_{\text{max}}$ is estimated from numerical
studies and is found to be $\theta_{\text{max}}\simeq 0.1$, given that larger
bound values do not contribute much to the integral. This result is in good
agreement with previous works such as [52, 53].

#### Final state lepton kinematics

In the case where the kinematics of the final state lepton are of interest, it
is instructive to perform similar computations based on its fractional energy
and emission angle, $(y,\ \psi)$, with $y=E_{f}/E_{i}$. The corresponding
double-differential cross-section then reads [52]

$\left(\frac{d^{2}\sigma}{d\cos{\psi}dy}\right)_{\text{WW}}=\frac{\alpha\chi}{\pi}\frac{E_{i}|\textbf{p}^{\prime}|}{1-y}\cdot\frac{d\sigma_{2\rightarrow
2}}{d(p\cdot p^{\prime})}\biggr{\rvert}_{t=t_{\text{min}}},$ (24)

for which we can express the Mandelstam variables in a similar fashion as in
Eqs. (18)

$\begin{split}t_{2}&\approx-
yE_{i}^{2}\psi^{2}-m_{f}^{2}\left(\frac{1-y}{y}\right)+m_{i}^{2}(1-y),\\\
\tilde{t}&\approx m_{\varphi}^{2}-t_{2},\\\
\tilde{s}&\approx\frac{\tilde{t}}{1-y},\\\
\tilde{u}&\approx-\frac{y\tilde{t}}{1-y}.\end{split}$ (25)

Plugging Eq. (25) into Eq. (20) gives

$\begin{split}|\mathcal{A}_{\varphi}^{2\rightarrow
2}|^{2}_{t_{\text{min}}}&=(1-y)^{2}\\\
\times\bigg{\\{}\frac{1}{2y}&-\frac{\Delta m^{2}}{yt}+\frac{\Delta
m^{2}[m_{f}^{2}+y(m_{i}^{2}y+\Delta
m^{2})]}{y^{2}t^{2}}\bigg{\\}}.\end{split}$ (26)

Combining Eqs. (19) and (26), the double-differential cross-section of Eq.
(24) reads

$\begin{split}\left(\frac{d^{2}\sigma}{d\cos{\psi}dy}\right)_{\text{WW}}&=\frac{\alpha^{2}h^{2}}{2\pi}E_{i}^{2}\beta_{f}\chi^{\text{WW}}(1-y)^{3}\\\
\times\bigg{(}\frac{1}{2\tilde{t}^{2}}-\frac{\Delta
m^{2}}{\tilde{t}^{3}}&+\frac{\Delta m^{2}[m_{f}^{2}+y(m_{i}^{2}y+\Delta
m^{2})]}{y\tilde{t}^{4}}\bigg{)},\end{split}$ (27)

where $\beta_{f}=(1-m_{f}^{2}/E_{f}^{2})^{1/2}$.

###  II.3 Numerical integration

The accuracy of the phase space approximation in the WW approach is estimated
through comparison with the ETL computations and shown in Fig. 1.

Figure 1: Top: Single-differential cross-sections as a function of the
fractional scalar boson $\varphi$ energy, evaluated at ETL (black line) and in
both the WW and IWW approaches (respectively green and blue dash lines). Eq.
(12) is evaluated with the Mathematica program [57]. Bottom: Single-
differential cross-section as a function of the outgoing lepton emission
angle, $\psi$, in both the WW and IWW approaches. The mass range of the scalar
boson spans from the $m_{\varphi}\geq m_{\mu}\simeq 110$ MeV up to 1 GeV, with
$h=10^{-4}$.

The relevant expressions, Eqs. (12) and (23), are integrated numerically using
methods from Mathematica. The primary lepton is chosen to be an electron with
energy $E_{i}=100$ GeV. The mass range of the $\varphi$ boson spans from a few
hundred MeV up to 1 GeV, the lower bound being justified by the kinematics
requiring $m_{\varphi}>m_{\mu}$. For completeness, the improved Weiszäcker-
Williams (IWW) approach is also shown, as it provides faster run-time
performance to the cost of loss of accuracy to the ETL [51, 52] (see Appendix
A). The relative error, expressed as
$(\mathcal{O}_{\text{exact}}-\mathcal{O}_{\text{approx.}})/\mathcal{O}_{\text{exact}}$,
is computed for both the IWW and WW approximations. In the WW approach, for
the whole $m_{\varphi}$ range, an error of $\lesssim 2\%$ is found with
respect to the ETL. In the IWW case, the relative error is significant, in
particular in the boundary regions, $x\rightarrow 0$ and $x\rightarrow 1$, due
to the flux integral simplifications.
Because of the importance of the final-state lepton kinematics in the LFC
process $e\rightarrow\mu$, the scattered angle of the muon is also shown in
Fig. 1 as a function of the $m_{\varphi}$ mass. It is found that in both
approximations the differential cross-sections peak at $\psi\sim
m_{\varphi}/E_{i}$. Because the ETL computations were not performed for the
$(y,\ \psi)$ kinematical variables, the IWW and WW approaches are not compared
to them.
For completeness, the total cross-section as a function of the incoming
electron energy is shown in Fig. 2 for two selected scalar boson masses,
$m_{\varphi}=110$ MeV and $m_{\varphi}=1$ GeV. While the WW reproduces
accurately the ETL with a precision of $\mathcal{O}(\leq 2\%)$, the IWW
approximation overestimates at typical NA64 electron beam energy ($\sim 100$
GeV) the tree-level computations by a factor $\sim 2$ and $\sim 3$
respectively for the lower and higher mass choices.

Figure 2: Total cross-sections evaluated at ETL (black line) and in both the
WW and IWW approaches (respectively green and blue dash lines) for different
initial electron energy. Top: $m_{\varphi}=110$ MeV. Bottom: $m_{\varphi}=1$
GeV.

##  III Search strategy at NA64

###  III.1 Experimental setup

NA64 is a fixed target experiment at the CERN Super Proton Synchrotron (SPS)
accelerator searching for Dark Sectors in the collisions of high energy lepton
beams with an active target [46]. The experiment runs in two modes, NA64e [60,
47, 61] and NA64$\mu$ [62, 63], with different detector setups, using the 100
GeV H4 beamline with electrons/positrons and the 160 GeV M2 muon beam
respectively. Since $\varphi$ mediates a non-diagonal interaction between a
muon and an electron, it could be produced in both NA64 configurations via the
bremsstrahlung-like process outlined in section II. This paper focuses only on
the sensitivity to $\varphi$ production using the H4 beamline.
The experimental technique for the search of $\varphi$ in the collision
$e^{-}N\rightarrow\mu N\varphi$ relies on the missing energy technique
developed by NA64 since $\varphi$ decays invisibly. Single $e^{-}$ originating
from the slow extracted beam are tagged through a set of scintillator and veto
counters. A magnet spectrometer consisting of two dipole magnets and tracking
detectors (Micromegas chambers (MMs), straw tubes (ST) and gas electron
multiplier (GEM)), measures the incoming momentum with a precision of $\delta
p/p\simeq 1\%$ [64]. The hadron contamination at the H4 beamline with 100 GeV
electrons is at the level of $\pi/e^{-}<10^{-2}$, with a $\sim 3\%$ kaon
contamination of the pion rate [65]. The synchrotron radiation detector (SRD)
identifies the electrons over the remaining hadrons. The SRD consists of an
array of a lead-scintillator sandwich calorimeter of a fine segmentation
placed after the bending magnet. It measures the synchrotron radiation of the
incoming high energy electrons through the magnetic field [66]. A 95$\%$
electron identification is achieved, thus suppressing the remaining hadron
contamination in the H4 beam down to a level of $\sim 2\times 10^{-5}$ [47].
The electrons collide with the active target, a 40 radiation length ($X_{0}$)
high-efficiency shashlik electromagnetic calorimeter (ECAL) made by sandwiched
lead-scintillator plates where $\varphi$ dark scalar would be produced. A
large high-efficiency VETO counter and three 7.5 nuclear interaction length
($\lambda_{\text{int}}$) iron hadronic calorimeters (HCALs) complete the setup
hermeticity, detecting charged and neutral secondaries produced from the
electron interaction in the target and measuring any energy leakage.
A key ingredient of this search is the addition of a second spectrometer
consisting of an MBPL dipole magnet together with four MMs trackers placed
after the last HCAL of the NA64 setup. A signal-like event will be missing
energy carried away by the $\varphi$ produced in the target (see Fig. 1) and a
muon in the final state. Such muon would need to traverse the three HCALs,
leaving $\sim 2.5$ GeV per module corresponding to the energy of a minimum
ionizing particle (MIP). The magnet spectrometer will allow us to identify it
and reconstruct its momenta. The missing energy is measured as the difference
between the incoming beam electron’s energy and the total energy carried away
by the final muon:

$E_{\text{miss}}=E_{e}^{i}-E_{\text{ECAL}}^{f}-E_{\text{HCAL}}^{f}-E_{\mu}^{f}.$
(28)

$E_{e}^{i}$ refers to the initial electron energy, $E_{\text{HCAL}}^{f}$ and
$E_{\text{ECAL}}^{f}$ are the energy in the HCAL and ECAL calorimeters
respectively, and $E_{\mu}^{f}$ is the energy of the final state muon
reconstructed in the magnet spectrometer.
A schematic representation of the NA64 setup focusing on the part after the
target is shown in Fig. 3.

Figure 3: Schematic view of the setup at NA64e optimized for
$e\rightarrow\mu+\varphi$ searches. The MBPL magnet combined with four MM
trackers constitute the magnet spectrometer. The whole set-up is not shown in
the diagram focusing on the relevant part for the $e$ to $\mu$ process. A
detailed description of the full setup can be found in [47].

###  III.2 Signal yield estimates

For a realistic estimate of the NA64 experiment sensitivity to the scalar
boson from the $e\rightarrow\mu$ LFC bremsstrahlung-like process, the full
physics of the underlying model is implemented in the MC using the setup
illustrated in Fig. 3. In particular, for the signal simulation, we use the
DMG4 package [48, 49] as it provides an application programming interface to
GEANT4 [44, 45]. The double- and single-differential cross-sections derived in
section II are implemented in the MC simulation framework. The related
distributions are used to perform the sampling of the final states’ kinematics
using a Von Neumann-based accept-reject method [67]. The mean free path of the
interaction is estimated at run-time for each particle step within the target
based on pre-computed tabulated ratios of the IWW approximation to the ETL
value of the total cross-section (see Fig. 4). This is done for a wide range
of $m_{\varphi}$ masses and initial electron energies.
Having implemented the underlying production mechanism in the simulations, the
total number of expected signal events, $\mathcal{N}_{\varphi}$, is given by

$\mathcal{N}_{\varphi}=N_{\text{EOT}}\times\frac{\rho\mathcal{N}_{A}}{A}\sum_{k}\sigma_{e\rightarrow\varphi\mu}(E_{k})\Delta
L_{k},$ (29)

where $N_{\text{EOT}}$ is the number of electrons on target (EOT), $\rho$ and
$A$ the target density and atomic weight, $\mathcal{N}_{A}$ the Avogadro
number, $\sigma_{e\rightarrow\varphi\mu}$ the production cross-section at the
$k-$step within the target evaluated at an electron energy $E_{k}$ and $\Delta
L_{k}$ the step length. Note that the total cross-section depends also on the
values of the mass of the $\varphi$ boson and coupling $h$.

Figure 4: Ratios of the IWW approximated total cross-section to the ETL value
for different initial electron energies and mass values.

Based on the simulation results, the signature of $\varphi$ in the
experimental setup was characterized. Signal-like events were identified
according to the following criteria: (i) One and only one particle compatible
with the trajectory expected for a $\mu^{-}$ reconstructed in the magnet
spectrometer. (ii) The energy deposit in each HCAL module and in the VETO is
consistent with that of a single minimum ionizing particle. The minimum
required energy for a muon to pass all three HCALs at NA64 is approximately
7.5 GeV. (iii) A missing energy larger than 20 GeV (see Eq. (28)). This value
was chosen to optimize the ratio of signal over background, as it will be
explained in more detail in the next sections.
Since $\varphi$ is likely to carry most of the initial electron’s energy (see
Fig. 1), many final state muons have too little energy to overcome the
calorimeters, and are either fully absorbed in the ECAL or stopped by the HCAL
modules. As the mass $m_{\varphi}$ increases, the differential cross section
is increasingly peaked towards $x=1$ (Fig. 1), meaning even fewer signal muons
are energetic enough to penetrate all HCALs. To be selected as candidate
signal events, the muons produced in the interaction will exit the last HCAL,
and will pass through the magnet gap hitting all MM trackers in the magnet
spectrometer. The signal efficiency as a function of the mass ($m_{\varphi}$),
is shown in Fig. 5 after applying the selection cuts (i-iii) described
previously.

Figure 5: Signal selection efficiency with respect to the mass of $\varphi$.
As the mass increases, the final state muon is likely to carry less of the
incoming electron energy (see Fig. 1), increasing the likelihood of the muon
to be absorbed in the ECAL or HCALs and not be reconstructed in the magnet
spectrometer.

###  III.3 Potential background sources

Standard Model events leading to missing energy and a final state muon in the
magnet spectrometer can mimic the signal signature. The two dominant
background sources in this type of processes are the dimuon production in the
target and muonic decays from hadronic beam contamination,
$h\rightarrow\mu+X$. Within this work, they have been studied in detail using
the MC framework described in the previous section.
Events with dimuon production $\gamma\rightarrow\mu^{+}\,\mu^{-}$ originating
from the shower of the primary electron in the ECAL can imitate the signal
characteristics if one and only one (anti-)muon of the pair is reconstructed
in the magnet spectrometer. The dimuon production at this energy is at the
level of $1.05\times 10^{-4}$ per electron on target (EOT). The HCAL signature
of signal-like events and such background events is shown in Fig. 6. The
signal spectrum is centered at the single MIP energy deposit of $\sim 2.5$
GeV, and is overall reduced as the muons are absorbed deeper into the HCAL
setup. The dimuon spectrum displays a larger peak deposit in the first HCAL
due to the double MIP deposit, but migrates towards a single MIP spectrum as
one of the muons is absorbed.

Figure 6: Deposited energy without any selection cuts in each HCAL modules 0-2
$(\mathrm{a})-(\mathrm{c})$ from simulated signal (blue circle) and dimuon
background events (red plus). With increasing HCAL module the signal peak is
reduced due to absorption. The background spectrum is likewise reduced and
shifted towards a single MIP deposit ($\sim 2.5$ GeV) as one of the dimuons is
absorbed.

Dimuon events where at least one muon is reconstructed in the magnet
spectrometer can generally be sorted into three categories: (I) One (anti)muon
is absorbed in the ECAL or HCALs while the other is reconstructed in the
magnet spectrometer. Such an event is not problematic since the energy is
fully contained and thus would not be tagged as a missing energy candidate
(II) Both muons pass through the HCALs and are reconstructed in the magnet
spectrometer. The requirements (i-iii) on single MIP energy deposit, missing
energy, and single track in the magnet remove such events. (III) Both muons
penetrate the HCALs but only one is reconstructed in the magnet. These events
have missing energy and only one track in the magnet but a double MIP energy
deposit in VETO and HCALs, and thus do not pass the selection criteria (ii).
After applying the selection criteria (i-iii), the missing energy variable
separates the dimuon background and signal sample as seen in Fig. 7. For a
simulated background data set equivalent to $1.5\times 10^{11}$ EOT, the event
selection shown in the previous section removes all dimuon events.
The other main background source arises from the muonic in-flight decays of
the remaining hadrons in the beam, $h\rightarrow\mu+X$. These events yield a
single muon and have also missing energy carried away by the accompanying
neutrino. A pion and kaon beam were simulated in the realistic MC simulation
framework described in section III to estimate the background level in these
searches considering the hadron admixture quoted in section III
($\pi/e^{-}<10^{-2}$). Muonic in-flight decays occurring before the SRD are
suppressed by the NA64 tagging system, so the dangerous events are those
occurring between the SRD and the ECAL. In this case, if the muon carries most
of the energy, the missing energy cut removes those events. However, events in
which the neutrino carries a sizeable missing energy, more than 30 GeV
according to simulations, will pass the selection being potential background.
After applying the event selection we get a survival rate of $9.6\times
10^{-5}$ for pions and $5.6\times 10^{-5}$ for kaons. The difference in the
rates associated with each hadron species is attributed to the underlying
decay process. The final states’ kinematics from pions and kaons result in
different angular acceptance within the magnet spectrometer. From one side,
the beam contamination rate of kaons is $3\%$ of that of pions. In addition,
although the probability for kaons to decay between the SRD and ECAL is higher
than for pions (0.32% compared to 0.04%), the opening angle of the kaon decay
products is larger than for pions. For that reason, the momentum
reconstruction in the Micromegas located after the bending will further
suppress those events. Finally, this result combined with the suppression of
the hadron component by the SRD gives a final background event rate of
$1.9\times 10^{-11}$ per EOT.

Figure 7: Missing energy
$E_{miss}=E_{e}^{i}-E_{HCAL}-E_{ECAL}-E_{\mu}^{\prime}$ for the simulated
dimuon background sample (left, green) and signal sample (right, blue) after
applying all selection criteria besides missing energy cut. The missing energy
for the dimuon background sample is centered around 0, meaning all those
events are fully contained in the detector.

##  IV Sensitivity

The sensitivity of the NA64 experiment to the new leptonic scalar boson
$\varphi$ is estimated based on the 90% confidence level (CL) upper limits
that would be obtained under the hypothesis that only SM events are observed.
Within this work, it is computed in the modified frequentist approach, using
the RooFit/RooStats-based [68, 69, 70] profile-likelihood ratio statistical
test in the asymptotic approximation [71]. The 90% CL projected sensitivity in
the $(m_{\varphi},\ h)$ plane based on the full signal simulation described in
section III.2 and considering the background level derived in section III.3 is
shown in Fig. 8. With a total statistics of $N_{\text{EOT}}=10^{11}$, NA64 has
the reach to probe a broad region of the parameter space compatible with a
fermionic DM relic abundance (see Eq. 2) for masses above the muon mass up to
few GeV. Additionally, it can also provide stringent constraints on the
allowed $(m_{\varphi},\ h)$ values in the explanation of the muon $g-2$ up to
masses $m_{\varphi}\lesssim 2m_{\mu}\sim\mathcal{O}(200\ \text{MeV})$ and
coupling $h\lesssim 7\times 10^{-4}$.

Figure 8: 90% CL projected sensitivity in mass-coupling $(m_{\varphi},\ h)$
parameter space of $\varphi$ given the background rate from section III.3. The
preferred band to solve the muon $g-2$ anomaly at $3\sigma$ from Eq. (3) is
shown in purple. The relic DM annihilation constraint is shown in black and is
computed through Eq. (2) with the parameter choices
$g_{\varphi\psi_{1}\psi_{2}}=0.1$ and $m_{1}=m_{2}=3m_{\varphi}$.

##  V Conclusion

In this paper, we presented the capabilities of fixed target experiments such
as NA64, to search for processes involving electron-to-muon transitions. In
particular, in this study, we have focused on using the 100 GeV electron beam
at H4 to search for a new lepton-flavor-changing scalar boson in the
Bremsstrahlung-like reaction $eN\rightarrow\mu N\varphi$. We have calculated
the exact-tree level differential and total cross-sections involved in such a
process and its WW phase-space approximation. The corresponding formulas have
been implemented in a realistic Monte Carlo framework using the GEANT4-based
DMG4 package as well as the NA64 experimental setup. The two main potential
background sources, dimuons production in the ECAL and the remaining hadron
decay in-flight have also been identified and studied. We have shown that the
most dangerous background arises from the remaining muons from the beam
hadrons decay-in-flight. This could be reduced further by improvements of the
SRD which are currently under study. In this paper, we have demonstrated that
adding a spectrometer after the last HCAL module to reconstruct the final
state muon momentum will allow to suppress the background at the required
level. The experiment projected sensitivity based on simulations indicates
that with $10^{11}$ EOT NA64 can probe the parameter space motivated by a
simultaneous solution of the muon g-2 anomaly and the DM problem for $0.2\leq
m_{\varphi}\leq 1$ GeV. The same process can be also probed in the NA64 muon
configuration searching for $\mu\rightarrow e$ and will be the subject of
future works. In addition, the positrons resulting from the electromagnetic
shower initiated by the incoming electron can also induce the reaction
$e^{+}N\rightarrow\mu^{+}N\varphi$. In this case, a signal-like event will be
identified as detecting an incoming $e^{-}$ and an outgoing $\mu^{+}$. This
study will also be covered in future works. Finally, this work opens the
possibility to further exploit other processes involving electron-to-muon
conversion using the NA64 missing energy and momentum techniques as ALP [72],
scalar[73] or dark vector [74, 75].

##  VI Acknowledgments

We acknowledge the members of the NA64 collaboration for fruitful discussions.
In particular, we thank N. V. Krasnikov for his helpful suggestions and
correspondence. The work of P. Crivelli, B. Banto Oberhauser, M. Mongillo, and
H. Sieber is supported by the SNSF Grants No. 186158, and No. 197346 No.
216602, No. 219485(Switzerland) and by ETH Grant No. 22-2 ETH-031. The work of
L. Molina-Bueno is supported by SNSF Grant No. 186158 (Switzerland),
RyC-030551-I, PID2021-123955NA-100 and CNS2022-135850 funded by MCIN/AEI/
10.13039/501100011033/FEDER, UE (Spain).

##  Appendix A Improved Weiszäcker-Williams approach

The computations performed in section II.2 can be further simplified under the
assumption that the lower bound of the effective photon flux is independent of
$x$ and $\theta_{k}$. As such we can pull
$\chi^{\text{WW}}\rightarrow\chi^{\text{IWW}}$ outside of the integral over
$\cos{\theta_{k}}$ in Eq. (23), permitting us to derive a fully analytical
differential cross section. The new bounds of $\chi^{\text{IWW}}$ are
$t_{\text{min}}=m_{\varphi}^{4}/(4E_{i}^{2})$ and
$t_{\text{max}}=m_{\varphi}^{2}+m_{i}^{2}$. The integral over $\theta_{k}$ in
Eq. (23) is analytically calculated by changing variables from
$\cos{\theta_{k}}$ to $\tilde{u}$, where we use $\tilde{u}$ evaluated at
$t_{\text{min}}$ just like in the WW approximation. This allows us to expand
$\sin{\theta_{k}}=\theta_{k}+\mathcal{O}(\theta_{k}^{3})$ leading to
$d\cos{\theta_{k}}\simeq-\theta_{k}\cdot d\theta_{k}$ and thus
$d\tilde{u}=2xE_{i}^{2}d\cos{\theta_{k}}$. After the substitution, the single-
differential cross-section takes the following form

$\left(\frac{d\sigma}{dx}\right)_{\text{IWW}}=\frac{\alpha^{2}h^{2}}{4\pi}\beta_{k}(1-x)\chi^{\text{IWW}}\int_{\tilde{u}_{\text{min}}}^{\tilde{u}_{\text{max}}}\,d\tilde{u}\
\frac{|\mathcal{A}_{\varphi}^{2\rightarrow
2}|^{2}_{t_{\text{min}}}}{\tilde{u}^{2}},$ (30)

where $\tilde{u}_{\text{max}}=\tilde{u}(0)$ and
$\tilde{u}_{\text{min}}=\tilde{u}(\theta_{\text{max}})$. Using the expression
for $|\mathcal{A}_{\varphi}^{2\rightarrow 2}|^{2}$ in Eq. (21) the integral in
Eq. (30) can finally be solved analytically as

$\begin{split}\int_{\tilde{u}_{\text{min}}}^{\tilde{u}_{\text{max}}}&\,d\tilde{u}\
\frac{|\mathcal{A}_{\varphi}^{2\rightarrow
2}|^{2}_{t_{\text{min}}}}{\tilde{u}^{2}}=\Bigg{\\{}\frac{x^{2}}{-2(1-x)\tilde{u}}+\frac{\Delta
m^{2}x}{-2\tilde{u}^{2}}\\\ &+\Delta
m^{2}\frac{m_{\varphi}^{2}(1-x)+x[m_{f}^{2}+m_{i}^{2}(x-1)]}{-3\tilde{u}^{3}}\Bigg{\\}}\biggr{\rvert}_{\tilde{u}_{\text{min}}}^{\tilde{u}_{\text{max}}}.\end{split}$
(31)

Similarly, the IWW approximation can be applied to the single-differential
cross-section with respect to the fractional lepton energy, $y$. Integrating
Eq. (27) over $\cos{\psi}$ gives the differential cross-section

$\begin{split}\left(\frac{d\sigma}{dy}\right)_{\text{IWW}}&=\frac{\alpha^{2}h^{2}}{2\pi}E_{i}^{2}\beta_{f}\chi^{\text{IWW}}(1-y)^{3}\int_{0}^{\psi_{\text{max}}}\,d\cos{\psi}\
\\\ \times\bigg{\\{}\frac{1}{2\tilde{t}^{2}}&-\frac{\Delta
m^{2}}{\tilde{t}^{3}}+\frac{\Delta m^{2}[m_{f}^{2}+y(m_{i}^{2}y+\Delta
m^{2})]}{y\tilde{t}^{4}}\bigg{\\}}.\end{split}$ (32)

As in section II.2, the maximal $\psi$ angle is estimated through numerical
integration and bound dependence on the contribution to the integral. It is
found $\psi_{max}=0.3$. We thus perform the integral in Eq. 32 analytically
through the change of variables from $\cos{\psi}\rightarrow\tilde{t}$. Since
we use the IWW approximation we can assume $d\cos{\psi}\simeq-\psi d\psi$. The
Jacobian of the transformation is $d\tilde{t}\simeq 2yE_{i}^{2}\psi
d\psi-2yE_{i}^{2}d\cos{\psi}$. Eq. (32) becomes

$\begin{split}\left(\frac{d\sigma}{dy}\right)_{\text{IWW}}&=\frac{\alpha^{2}h^{2}}{4\pi}\beta_{f}\chi^{\text{IWW}}\frac{(1-y)^{3}}{y}\\\
\times\Bigg{\\{}\frac{-1}{2\tilde{t}}&+\frac{\Delta
m^{2}}{2\tilde{t}^{2}}-\frac{\Delta m^{2}[m_{f}^{2}+y(m_{i}^{2}y+\Delta
m^{2})]}{3y\tilde{t}^{3}}\Bigg{\\}}\biggr{\rvert}_{\tilde{t}_{\text{min}}}^{\tilde{t}_{\text{max}}},\end{split}$
(33)

with $t_{\text{min}}$ and $t_{\text{max}}$ defined previously. For
completeness, it is worth noting that the IWW single-differential cross-
section for the outgoing lepton scattered angle is given by integrating Eq.
(27) over y.

## References

  * Rubin and Ford [1970] V. C. Rubin and W. K. Ford, Jr., Astrophys. J. 159, 379 (1970).
  * Freeman [1970] K. C. Freeman, Astrophys. J. 160, 811 (1970).
  * van Albada  _et al._ [1985] T. S. van Albada, J. N. Bahcall, K. Begeman, and R. Sancisi, Astrophys. J. 295, 305 (1985).
  * Begeman  _et al._ [1991] K. G. Begeman, A. H. Broeils, and R. H. Sanders, Mon. Not. Roy. Astron. Soc. 249, 523 (1991).
  * Zwicky [1937] F. Zwicky, Astrophys. J. 86, 217 (1937).
  * Tyson  _et al._ [1998] J. A. Tyson, G. P. Kochanski, and I. P. Dell’Antonio, Astrophys. J. Lett. 498, L107 (1998), arXiv:astro-ph/9801193 .
  * Clowe  _et al._ [2006] D. Clowe, M. Bradac, A. H. Gonzalez, M. Markevitch, S. W. Randall, C. Jones, and D. Zaritsky, Astrophys. J. Lett. 648, L109 (2006), arXiv:astro-ph/0608407 .
  * Ade  _et al._ [2014] P. A. R. Ade _et al._ (Planck), Astron. Astrophys. 571, A16 (2014), arXiv:1303.5076 [astro-ph.CO] .
  * Aghanim  _et al._ [2020] N. Aghanim _et al._ (Planck), Astron. Astrophys. 641, A6 (2020), [Erratum: Astron.Astrophys. 652, C4 (2021)], arXiv:1807.06209 [astro-ph.CO] .
  * Refregier [2003] A. Refregier, Ann. Rev. Astron. Astrophys. 41, 645 (2003), arXiv:astro-ph/0307212 .
  * Cole  _et al._ [2005] S. Cole _et al._ (2dFGRS), Mon. Not. Roy. Astron. Soc. 362, 505 (2005), arXiv:astro-ph/0501174 .
  * Abbott  _et al._ [2022] T. M. C. Abbott _et al._ (DES), Phys. Rev. D 105, 023520 (2022), arXiv:2105.13549 [astro-ph.CO] .
  * Kobzarev  _et al._ [1966] I. Y. Kobzarev, L. B. Okun, and I. Y. Pomeranchuk, Sov. J. Nucl. Phys. 3, 837 (1966).
  * Blinnikov and Khlopov [1982] S. I. Blinnikov and M. Y. Khlopov, Sov. J. Nucl. Phys. 36, 472 (1982).
  * Foot  _et al._ [1991] R. Foot, H. Lew, and R. R. Volkas, Phys. Lett. B 272, 67 (1991).
  * Hodges [1993] H. M. Hodges, Phys. Rev. D 47, 456 (1993).
  * Berezhiani  _et al._ [1996] Z. G. Berezhiani, A. D. Dolgov, and R. N. Mohapatra, Phys. Lett. B 375, 26 (1996), arXiv:hep-ph/9511221 .
  * Jaeckel  _et al._ [2020] J. Jaeckel, M. Lamont, and C. Vallée, Nature Phys. 16, 393 (2020).
  * Lanfranchi  _et al._ [2021] G. Lanfranchi, M. Pospelov, and P. Schuster, Ann. Rev. Nucl. Part. Sci. 71, 279 (2021), arXiv:2011.02157 [hep-ph] .
  * Gninenko and Krasnikov [2022] S. N. Gninenko and N. V. Krasnikov (NA64), Phys. Rev. D 106, 015003 (2022), arXiv:2202.04410 [hep-ph] .
  * Bayes  _et al._ [2015] R. Bayes _et al._ (TWIST), Phys. Rev. D 91, 052020 (2015), arXiv:1409.0638 [hep-ex] .
  * Workman and Others [2022] R. L. Workman and Others (Particle Data Group), PTEP 2022, 083C01 (2022).
  * Aguillard  _et al._ [2023] D. P. Aguillard _et al._ (Muon g-2), (2023), arXiv:2308.06230 [hep-ex] .
  * Davier  _et al._ [2017] M. Davier, A. Hoecker, B. Malaescu, and Z. Zhang, Eur. Phys. J. C 77, 827 (2017), arXiv:1706.09436 [hep-ph] .
  * Keshavarzi  _et al._ [2018] A. Keshavarzi, D. Nomura, and T. Teubner, Phys. Rev. D 97, 114025 (2018), arXiv:1802.02995 [hep-ph] .
  * Colangelo  _et al._ [2019] G. Colangelo, M. Hoferichter, and P. Stoffer, JHEP 02, 006 (2019), arXiv:1810.00007 [hep-ph] .
  * Hoferichter  _et al._ [2019] M. Hoferichter, B.-L. Hoid, and B. Kubis, JHEP 08, 137 (2019), arXiv:1907.01556 [hep-ph] .
  * Davier  _et al._ [2020] M. Davier, A. Hoecker, B. Malaescu, and Z. Zhang, Eur. Phys. J. C 80, 241 (2020), [Erratum: Eur.Phys.J.C 80, 410 (2020)], arXiv:1908.00921 [hep-ph] .
  * Keshavarzi  _et al._ [2020] A. Keshavarzi, D. Nomura, and T. Teubner, Phys. Rev. D 101, 014029 (2020), arXiv:1911.00367 [hep-ph] .
  * Kurz  _et al._ [2014] A. Kurz, T. Liu, P. Marquard, and M. Steinhauser, Phys. Lett. B 734, 144 (2014), arXiv:1403.6400 [hep-ph] .
  * Melnikov and Vainshtein [2004] K. Melnikov and A. Vainshtein, Phys. Rev. D 70, 113006 (2004), arXiv:hep-ph/0312226 .
  * Masjuan and Sanchez-Puertas [2017] P. Masjuan and P. Sanchez-Puertas, Phys. Rev. D 95, 054026 (2017), arXiv:1701.05829 [hep-ph] .
  * Colangelo  _et al._ [2017] G. Colangelo, M. Hoferichter, M. Procura, and P. Stoffer, JHEP 04, 161 (2017), arXiv:1702.07347 [hep-ph] .
  * Hoferichter  _et al._ [2018] M. Hoferichter, B.-L. Hoid, B. Kubis, S. Leupold, and S. P. Schneider, JHEP 10, 141 (2018), arXiv:1808.04823 [hep-ph] .
  * Gérardin  _et al._ [2019] A. Gérardin, H. B. Meyer, and A. Nyffeler, Phys. Rev. D 100, 034520 (2019), arXiv:1903.09471 [hep-lat] .
  * Bijnens  _et al._ [2019] J. Bijnens, N. Hermansson-Truedsson, and A. Rodríguez-Sánchez, Phys. Lett. B 798, 134994 (2019), arXiv:1908.03331 [hep-ph] .
  * Colangelo  _et al._ [2020] G. Colangelo, F. Hagelstein, M. Hoferichter, L. Laub, and P. Stoffer, JHEP 03, 101 (2020), arXiv:1910.13432 [hep-ph] .
  * Colangelo  _et al._ [2014] G. Colangelo, M. Hoferichter, A. Nyffeler, M. Passera, and P. Stoffer, Phys. Lett. B 735, 90 (2014), arXiv:1403.7512 [hep-ph] .
  * Blum  _et al._ [2020] T. Blum, N. Christ, M. Hayakawa, T. Izubuchi, L. Jin, C. Jung, and C. Lehner, Phys. Rev. Lett. 124, 132002 (2020), arXiv:1911.08123 [hep-lat] .
  * Aoyama  _et al._ [2012] T. Aoyama, M. Hayakawa, T. Kinoshita, and M. Nio, Phys. Rev. Lett. 109, 111808 (2012), arXiv:1205.5370 [hep-ph] .
  * Aoyama  _et al._ [2019] T. Aoyama, T. Kinoshita, and M. Nio, Atoms 7, 28 (2019).
  * Gnendiger  _et al._ [2013] C. Gnendiger, D. Stöckinger, and H. Stöckinger-Kim, Phys. Rev. D 88, 053005 (2013), arXiv:1306.5546 [hep-ph] .
  * Ignatov  _et al._ [2023] F. V. Ignatov _et al._ (CMD-3), (2023), arXiv:2302.08834 [hep-ex] .
  * Agostinelli  _et al._ [2003] S. Agostinelli _et al._ (GEANT4), Nucl. Instrum. Meth. A 506, 250 (2003).
  * Allison  _et al._ [2016] J. Allison _et al._ , Nucl. Instrum. Meth. A 835, 186 (2016).
  * Gninenko [2014] S. N. Gninenko, Phys. Rev. D 89, 075008 (2014), arXiv:1308.6521 [hep-ph] .
  * Andreev  _et al._ [2023a] Y. M. Andreev _et al._ (NA64), Phys. Rev. Lett. 131, 161801 (2023a), arXiv:2307.02404 [hep-ex] .
  * Bondi  _et al._ [2021] M. Bondi, A. Celentano, R. R. Dusaev, D. V. Kirpichnikov, M. M. Kirsanov, N. V. Krasnikov, L. Marsicano, and D. Shchukin, Comput. Phys. Commun. 269, 108129 (2021), arXiv:2101.12192 [hep-ph] .
  * Oberhauser  _et al._ [2024] B. B. Oberhauser _et al._ , (2024), arXiv:2401.12573 [hep-ph] .
  * Liu  _et al._ [2017] Y.-S. Liu, D. McKeen, and G. A. Miller, Phys. Rev. D 95, 036010 (2017), arXiv:1609.06781 [hep-ph] .
  * Liu and Miller [2017] Y.-S. Liu and G. A. Miller, Phys. Rev. D 96, 016004 (2017), arXiv:1705.01633 [hep-ph] .
  * Kirpichnikov  _et al._ [2021] D. V. Kirpichnikov, H. Sieber, L. M. Bueno, P. Crivelli, and M. M. Kirsanov, Phys. Rev. D 104, 076012 (2021), arXiv:2107.13297 [hep-ph] .
  * Sieber  _et al._ [2023] H. Sieber, D. V. Kirpichnikov, I. V. Voronchikhin, P. Crivelli, S. N. Gninenko, M. M. Kirsanov, N. V. Krasnikov, L. Molina-Bueno, and S. K. Sekatskii, Phys. Rev. D 108, 056018 (2023), arXiv:2305.09015 [hep-ph] .
  * Mertig  _et al._ [1991] R. Mertig, M. Bohm, and A. Denner, Comput. Phys. Commun. 64, 345 (1991).
  * Shtabovenko  _et al._ [2016] V. Shtabovenko, R. Mertig, and F. Orellana, Comput. Phys. Commun. 207, 432 (2016), arXiv:1601.01167 [hep-ph] .
  * Shtabovenko  _et al._ [2020] V. Shtabovenko, R. Mertig, and F. Orellana, Comput. Phys. Commun. 256, 107478 (2020), arXiv:2001.04407 [hep-ph] .
  * [57] W. R. Inc., “Mathematica, Version 13.2,” Champaign, IL, 2022.
  * Davoudiasl  _et al._ [2023] H. Davoudiasl, R. Marcarelli, and E. T. Neil, JHEP 02, 071 (2023), arXiv:2112.04513 [hep-ph] .
  * Bjorken  _et al._ [2009] J. D. Bjorken, R. Essig, P. Schuster, and N. Toro, Phys. Rev. D 80, 075018 (2009), arXiv:0906.0580 [hep-ph] .
  * Banerjee  _et al._ [2017] D. Banerjee _et al._ (NA64), Phys. Rev. Lett. 118, 011802 (2017), arXiv:1610.02988 [hep-ex] .
  * Andreev  _et al._ [2023b] Y. M. Andreev _et al._ (NA64), (2023b), arXiv:2308.15612 [hep-ex] .
  * Gninenko  _et al._ [2015] S. N. Gninenko, N. V. Krasnikov, and V. A. Matveev, Phys. Rev. D 91, 095015 (2015), arXiv:1412.1400 [hep-ph] .
  * Sieber  _et al._ [2022] H. Sieber, D. Banerjee, P. Crivelli, E. Depero, S. N. Gninenko, D. V. Kirpichnikov, M. M. Kirsanov, V. Poliakov, and L. Molina Bueno, Phys. Rev. D 105, 052006 (2022), arXiv:2110.15111 [hep-ex] .
  * Banerjee  _et al._ [2015] D. Banerjee, P. Crivelli, and A. Rubbia, Adv. High Energy Phys. 2015, 105730 (2015).
  * Andreev  _et al._ [2023c] Y. M. Andreev _et al._ , Nucl. Instrum. Meth. A 1057, 168776 (2023c), arXiv:2305.19411 [hep-ex] .
  * Depero  _et al._ [2017] E. Depero _et al._ , Nucl. Instrum. Meth. A 866, 196 (2017), arXiv:1703.05993 [physics.ins-det] .
  * Kalos and Whitlock [2008] M. H. Kalos and P. A. Whitlock,  _Monte Carlo methods_ (WILEY-VCH, Weinheim, 2008) pp. XII, 203 S.
  * Verkerke and Kirkby [2003] W. Verkerke and D. P. Kirkby, eConf C0303241, MOLT007 (2003), arXiv:physics/0306116 .
  * Wolffs  _et al._ [2022] Z. Wolffs, P. Bos, C. Burgard, E. Michalainas, L. Moneta, J. Rembser, and W. Verkerke, PoS ICHEP2022, 249 (2022).
  * Moneta  _et al._ [2010] L. Moneta, K. Belasco, K. S. Cranmer, S. Kreiss, A. Lazzaro, D. Piparo, G. Schott, W. Verkerke, and M. Wolf, PoS ACAT2010, 057 (2010), arXiv:1009.1003 [physics.data-an] .
  * Cowan  _et al._ [2011] G. Cowan, K. Cranmer, E. Gross, and O. Vitells, Eur. Phys. J. C 71, 1554 (2011), [Erratum: Eur.Phys.J.C 73, 2501 (2013)], arXiv:1007.1727 [physics.data-an] .
  * Jho  _et al._ [2022] Y. Jho, S. Knapen, and D. Redigolo, JHEP 10, 029 (2022), arXiv:2203.11222 [hep-ph] .
  * Araki  _et al._ [2021] T. Araki, K. Asai, and T. Shimomura, JHEP 11, 082 (2021), arXiv:2107.07487 [hep-ph] .
  * Zhevlakov  _et al._ [2024] A. S. Zhevlakov, D. V. Kirpichnikov, and V. E. Lyubovitskij, Phys. Rev. D 109, 015015 (2024), arXiv:2307.10771 [hep-ph] .
  * Heeck [2016] J. Heeck, Phys. Lett. B 758, 101 (2016), arXiv:1602.03810 [hep-ph] .


# Data-driven Thermal Anomaly Detection for Batteries using Unsupervised Shape
Clustering

Xiaojun Li*, Jianwei Li, Ali Abdollahi, Trevor Jones and Asif Habeebullah
This work has been submitted to the IEEE for possible publication. The Authors
are with Gotion Inc, 48660 Kato Rd, Fremont, CA, 94538 USA.Contact:{t.li,
j.li, a.abdollahi, t.jones, a<EMAIL_ADDRESS>
###### Abstract

For electric vehicles (EV) and energy storage (ES) batteries, thermal runaway
is a critical issue as it can lead to uncontrollable fires or even explosions.
Thermal anomaly detection can identify problematic battery packs that may
eventually undergo thermal runaway. However, there are common challenges like
data unavailability, environment variations, and battery aging. We propose a
data-driven method to detect battery thermal anomaly based on comparing shape-
similarity between thermal measurements. Based on their shapes, the
measurements are continuously being grouped into different clusters. Anomaly
is detected by monitoring deviations within the clusters. Unlike model-based
or other data-driven methods, the proposed method is robust to data loss and
requires minimal reference data for different pack configurations. As the
initial experimental results show, the method not only can be more accurate
than the onboard BMS, but also can detect unforeseen anomalies at the early
stage.

##  I Introduction

###  I-A Background

Conventional anomaly detection methods for batteries usually depend on
thresholds or lookup tables, often determined by lab-testing of sample
batteries, and may not apply to the individual battery that operates under
different conditions. On the other hand, advanced anomaly detection methods,
such as machine learning-based algorithms, require significant computational
resources. Traditional battery management systems (BMS) deployed on electric
vehicles or energy storage systems are based on embedded micro-controllers,
which lack computing power or memory to execute these algorithms on-board.

For a cloud-based BMS, the gathered data is transmitted to a data-center for
further analysis, during which advanced algorithms can be utilized. As a part
of the fault diagnosis process, anomaly/fault detection is the most critical
step. Based on timely detection, the BMS or vehicle controller (VCU) of the EV
can take proper actions, which prevent a relative small issue from developing
into a server problem. There are three mainstream methods for battery
fault/anomaly detection: knowledge-based, model-based, and data-driven[1]. The
threshold-based method, which is the industry’s standard practice, can be
categorized as knowledge-based. In general, every battery manufacturer has its
own heuristic ”recipe” combined with testing data and engineering known-hows,
which can not be generalized easily. Model-based approaches often include a
physical model and an estimator[2]. In [3], a lumped battery thermal model
with time-variant internal resistance, altogether with an adaptive observer,
are used to estimate the core battery temperature. [4] uses a cell difference
model with an extended Kalman filter to estimate the micro-shot-circuit
current, which can be used for short circuit diagnosis. Some of the data-
driven approaches depend on the recurrent neural network[5, 6]. For example,
[5] uses a long short-term memory neural network to create residual signals
for battery surface temperatures. A thermal fault is set when the residual is
over a certain threshold. For model-based approaches, training data is needed
for finding the optimal model parameters for different types of batteries,
making it difficult and time-consuming to be implemented. Besides, data
anomalies like long-time unavailability and shifting can cause the model to
malfunction. Relatively speaking, training recurrent neural networks requires
fewer efforts. However, the neural network also depends on the continuous
signal influx. For both methods, battery degradation over time poses a
challenge. In conclusion, the following common issues still pose significant
challenges for battery anomaly detection.

Figure 1: Typical issues of thermal measurements for cloud-BMS data. This
figure is based on fleet EV battery data.

  * •

Loss of data/invalid data: For cost-saving, some wireless transmission modules
in BMS or vehicles usually use less reliable networks such as 3G or even 2G.
Due to network issues, loss of data or invalid data is common. As depicted in
Fig. 1, it may happen just intermittently or for an extended period. Data is
also not available when the vehicle is shut down. Model-based approaches, such
as Kalman filters and recurrent-neural networks, have internal states that
rely on continuous data input. They may tolerate intermittent data loss, but
long-term data loss will cause the model to lose track and give inaccurate
results.

  * •

Asynchronous data: Because of communication or sensor issues, one or several
signals may be sampled with some time-delay. Residual-based methods, which
evaluate all signals at the exact moment, is likely to generate false positive
with asynchronous data.

  * •

Battery aging: As the battery ages and deteriorates, the thermal and cell
voltage measurement will deviate from their nominal range. Both residual and
model-based method needs to adjust the threshold and parameters for them to
work correctly. However, estimating the battery’s state-of-health (SOH) itself
is a challenging task[7].

  * •

Environment variations: The environment variations, such as different geo-
locations, sensor locations, and seasonal temperature changes, cause similar
difficulties as battery aging. But they are more hard to be factored in the
model.

  * •

Lack of training data: Electric vehicle and energy storage systems have
different configurations. It is challenging to train the model so that it can
adequately identify anomalies for all configurations. A tailored model is more
accurate but is challenging to implement. Besides, labeled data may not be
available for some configurations.

###  I-B The proposed method

To address some of the challenges mentioned above, this paper proposes a data-
driven, model-free approach that monitors the shape-similarities across the
measurements to detect battery thermal anomalies. The proposed method does not
require the data to be continuous, making it robust to data loss and invalid
data. The same reason also alleviates the effects of battery aging or
environment variations since cells are likely to deteriorate or be influenced
by these variations as a whole. Asynchronous data issue is handled by the
shape-based distance measurement, which is invariant to signal shifting.
Furthermore, this method can be applied to different configurations easily
since it needs very little reference data.

Clustering is an unsupervised learning algorithm that does not require
training with labeled data. Therefore, this method can detect unforeseen
anomalies. Unlike the existing cross-similarity approaches[8], the proposed
method does not need feature extraction. Also, comparing shapes(which is
scaling-invariant) instead of values allows it to capture abnormality at the
early stage. Therefore, this new method has the potential to give early
warnings.

##  II The Anomaly Detection Method

Figure 2: The Proposed Methodology

###  II-A The K-shape Clustering Algorithm

Proposed by[9], the K-shape clustering algorithm is intended to be used for
time-series analysis. It adopts a cross-correlation sequence
($CC_{\omega}\left(\mathbf{x},\mathbf{y}\right)$) to determine signal
similarity. The shape-based distance function is given as:

$SBD(\vec{x},\vec{y})\ =\
1-\underset{\omega}{max}{\left(\frac{CC_{\omega}\left(\mathbf{x},\mathbf{y}\right)}{\sqrt{R_{0}\left(\mathbf{x},\mathbf{x}\right)R_{0}\left(\mathbf{y},\mathbf{y}\right)}}\right)}$
(1)

where $\vec{x},\vec{y}$ are the normalized time-series measurements, and
$R_{0}$ is the Rayleigh Quotient. The overall K-shape algorithm has two steps
per iteration, repeated until convergence or max iteration is reached. In the
first step (assignment), measurements are assigned to each cluster based on
their similarity to the centroid. In the second step (update), the clusters’
centroids are updated. The overall time complexity
$\mathcal{O}\left(max\left\\{(ncm)log(m),nm^{2},cm^{3}\right\\}\right)$ scales
linearly with the number of measurements($n$) and the number of clusters ($c$)
but increase noticeably with the number of time-steps ($m$). K-shape has been
applied to time-series analysis and forecast[10, 11, 12], including battery
cell voltage monitoring[13].

###  II-B Assumptions and Limitations

Similar to the advantages are given above, the proposed method relies on some
assumptions and has its limitations.

  1. 1.

The battery pack initially operates normally: We use the original operational
data to acquire the reference cluster membership. It is worth mentioning that
the proposed method only needs a small size of measurement data to extract the
infomation. For example, in Section III, only 4 minutes of normal operation
data is used.

  2. 2.

The anomaly does not occur in all measurements in the same manner: In theory,
anomalies can cause no deviation in the measurement clustering if they affect
all measurements simultaneously and in the same way without disrupting the
existing cluster memberships. Under such circumstances, the proposed method
would not be effective. However, it is unlikely to happen in real-life. A
typical anomaly first appears only in one or a few measurements and can be
detected by the proposed method.

###  II-C Methodology

Input: membership and distance for the $kth$ segment, which has $c$ clusters
($\mathbf{mem}[k],\mathbf{dist}[k]$)

Output: anomaly indicator and confidence level for the $kth$ segments
($p[k],d[k]$)

// initialization

$a[k]=False$

$p[k]=0$

// initialization

for  _$mem[k][i]$ in $\mathbf{mem}[k]$_ do

if  _$mem[k][i]\neq mem[k-1][i]$_ then

$a[k]=True$

return _$(a[k],p[k])$_

end for

for  _$dist[k][i]$ in $\mathbf{dist}[k]$_ do

_ _ if  _$dist[k-1][i] >\epsilon$ // avoid zero-division _

then

$p[k]+=(1-\alpha)\cdot{sat}\left(\frac{{dist}[k]}{dist[k-1]}\right)+\alpha\cdot{sat}\left(\frac{{dist}[k]}{{dist}[0]}\right)$

else

$p[k]+=sat\left(\frac{{dist}[k]}{{dist}[0]}\right)$

end for

return _$a[k],p[k]$_

Algorithm 1 ${a[k],p[k]=AnmlyChk(\mathbf{mem}[k],\mathbf{dist}[k])}$

As the flowchart depicted in Fig. 2 , the proposed anomaly detection method
contains the following stages:

  1. 1.

Segmentation/buffering: For offline implementation, existing data is segmented
into smaller pieces, each of which will be clustered accordingly. For online
implementation, data should be continuously collected and buffered as a
segment, which is then processed. Because the K-shape’s time complexity grows
noticeably as the number of samples increase, segment size should be limited
based on the sampling rate and data type. For example, temperature
measurements tend to changes slower than voltage measurements and a thermal
fault typically takes a longer time to develop than a voltage anomaly does.
Therefore, temperature measurement may need larger segments than voltage
measurement.

  2. 2.

Cleaning: In this stage, invalid data points, like those that are out of the
sensor’s measurement range, are removed. However, there is no need to fill in
the missing data due to this method’s advantage.

  3. 3.

Preprocess: This stage includes filtering and data normalization. Segments
where all signals do not show noticeable dynamics (changes) are excluded from
the clustering process.

  4. 4.

K-shape: In this stage, the K-shape algorithm given in equation 1 is applied
to the segment.

  5. 5.

Anomaly Confirmation: Details of the confirmation algorithm are given in
Algorithm 1. We look at two criteria. Firstly, if one or multiple measurements
changed its membership, it indicates an anomaly($a[k]$). If there is no change
in cluster membership, we check if there is a significant increase in fitting
errors/distance($p[k]$). The following total distance function[14] is used:

$dist=\frac{1}{n}\sum_{i\in
n}\left(\underset{j}{min}\left({SBD^{2}(x_{i},c_{j})}\right)\right)$ (2)

As depicted in Fig 2, during an iteration for the $kth$ segment, each cluster
($i$) is compared to the reference cluster ($dist[0]$) and the previous
cluster ($dist[k-1]$). The first part is for capturing the accumulated changes
associated with anomalies that developed gradually. For example, thermal
anomaly due to battery’s increased internal resistances. The second part is
for capturing incremental changes associated with anomaly that developed
abruptly, such as a thermal anomaly caused by short-circuiting.

##  III Experimental results

###  III-A Database and data sources

The R&D team at Gotion has built a cloud-based battery data platform. The data
platform receives and cleans battery data from different sources such as fleet
vehicles, onboard tests, and lab tests. Afterward, it is uploaded to a time-
series database. There are also software modules for data-visualization,
battery simulation, and analysis. Most of the algorithms are implemented in
Python. The proposed method is applied to battery data collected from fleet
vehicle batteries between 2019 and 2021. The data sampling rate is about 0.1
Hz. There are different types of battery chemistry ($LiNiMnCo$ and
$LiFePO_{4}$), pack configurations, and vehicle types.

(a)

(b)

(c)

Figure 3: Case I: (a) Thermal measurement data with the over-temperature fault. (b) The zoom-in view of fault occurrence, notice a temperature anomaly is detected by the proposed method 90min before it surpasses the threshold. (c) The plot of measurement shapes for the segment in which the anomaly is detected. TABLE I: Description of the Vehicle Battery Data field names | data type | value range | source
---|---|---|---
cell voltage | double | (0,5) V | BMS
temperature | double | (-50,100) C | BMS
pack current | double | (-500,500) A | BMS
battery faults | bool | [0,1] | BMS
battery status | enum | - | BMS

###  III-B Initial results and discussions

In this section, the testing results of two different battery packs are
presented and compared to the fault detection based on the on-board BMS. The
proposed algorithm is implemented in Python and deployed on a cloud computing
platform (AWS). For the following testing results, the segment size is set to
25, equivalent to 4.1 min. Notice that both cases have intermittent and
longtime data losses. In the following figures, the data gaps are connected
for better visualization.

####  III-B1 Case I

In the first case, the proposed method is applied to a battery pack
($LiFePO_{4}$) that undergoes an over-temperature anomaly. As depicted in Fig
3(a), on Oct 30th, the battery’s temperature near sensor #13 rises
significantly to over 70 C. While both the proposed method and BMS are able to
flag the anomaly, the new method is more than 90 min earlier. The detailed
timing difference between the two methods is illustrated in Fig 3(b). As the
figure shows, the BMS reports an over-temperature fault at 3:45 PM when the
maximum temperature is over 55 C. On the other hand, the proposed method send
an anomaly warning around 2:15 pm, just when sensor #13 start to depart from
the rest of the measurements. Fig 3(c) shows a shape plot of the segment where
the anomaly is detected. Clearly, one of the signal’s (#13) rising shape
stands out from the rest. As signal #13 continues to grow, its shapes for the
following segments become less steep. Therefore, the confident level ($d[k]$)
is the largest at the beginning of the anomaly. This explains why the new
method can send early warnings. In conclusion, case I validates that the
proposed method can send an early warning for battery over-temperature faults.

####  III-B2 Case II

The second case is from an EV with a two-pack configuration, in which two
battery packs ($LiNiMnCo$) are installed in different locations inside the
vehicle. As a result, the thermal measurements behave very differently. As Fig
4 shows, the temperature difference between the two packs grows noticeably in
8:00-9:00 AM and in 1:00-2:00 PM, when a large current is discharged from the
battery. In both cases, the onboard BMS reports thermal fault despite that
there is no real anomaly. The reason is that the on-board BMS uses hard
thresholds from a look-up table. Also depicted in Fig 4, the proposed method
successfully recognizes two clustering groups, including $\\{1,3,5,7\\}$ and
$\\{2,4,6,8\\}$ and does not report any anomaly since there is no discrepancy
found with the cluster groups. This testing case shows that compared to BMS,
the proposed method is more robust to variations caused by the pack design and
sensor location.

Figure 4: Case II: (a) thermal measurements and (b) current for a two-pack
battery. Notice the temperature difference grows under large current loads.

##  IV Conclusion and Future work

In this paper, we identify common issues in the field of cloud-based battery
anomaly/fault detection. Then, a method based on unsupervised shape-clustering
is proposed for detecting battery thermal anomalies. The proposed method does
not depend on a model or large training data. It also has several unique
advantages, such as the resilience to data loss and the capability of early
detection. Two test cases based on real vehicle data are studied. In one case,
the new method is found to be more accurate than the BMS when applied to a
multi-pack vehicle. In the other case study where a battery over-temperature
fault occurs, the proposed method is capable of flagging the anomaly at the
very early stage, more than 90 min ahead of the BMS.

Future works include the following. Firstly, this method is being currently
tested on large data sets. The rate of false positives and false negatives
needs to be investigated and compared to the BMS. Secondly, detection accuracy
can be further improved by making use of another measurement signals like BMS
status and cell voltages.

## References

  * [1] X. Hu, K. Zhang, K. Liu, X. Lin, S. Dey, and S. Onori, “Advanced Fault Diagnosis for Lithium-Ion Battery Systems: A Review of Fault Mechanisms, Fault Features, and Diagnosis Procedures,” _IEEE Industrial Electronics Magazine_ , vol. 14, no. 3, pp. 65–91, 2020.
  * [2] G. L. Plett, _Battery Management Systems, Volume I: Battery Modeling_ , ser. Artech House power engineering series. Artech House, 2015.
  * [3] X. Lin, H. E. Perez, J. B. Siegel, A. G. Stefanopoulou, Y. Li, R. D. Anderson, Y. Ding, and M. P. Castanier, “Online parameterization of lumped thermal dynamics in cylindrical lithium ion batteries for core temperature estimation and health monitoring,” _IEEE Transactions on Control Systems Technology_ , vol. 21, no. 5, pp. 1745–1755, 2013.
  * [4] W. Gao, Y. Zheng, M. Ouyang, J. Li, X. Lai, and X. Hu, “Micro-short-circuit diagnosis for series-connected lithium-ion battery packs using mean-difference model,” _IEEE Transactions on Industrial Electronics_ , vol. 66, no. 3, 2019.
  * [5] O. Ojo, H. Lang, Y. Kim, X. Hu, B. Mu, and X. Lin, “A Neural Network-Based Method for Thermal Fault Detection in Lithium-Ion Batteries,” _IEEE Transactions on Industrial Electronics_ , vol. 68, no. 5, pp. 4068–4078, 2020\.
  * [6] D. Li, Z. Zhang, P. Liu, Z. Wang, and L. Zhang, “Battery Fault Diagnosis for Electric Vehicles Based on Voltage Abnormality by Combining the Long Short-Term Memory Neural Network and the Equivalent Circuit Model,” _IEEE Transactions on Power Electronics_ , vol. 36, no. 2, pp. 1303–1315, 2 2021.
  * [7] A. Nuhic, T. Terzimehic, T. Soczka-Guth, M. Buchholz, and K. Dietmayer, “Health diagnosis and remaining useful life prognostics of lithium-ion batteries using data-driven methods,” _Journal of Power Sources_ , vol. 239, pp. 680–688, 2013. [Online]. Available: http://dx.doi.org/10.1016/j.jpowsour.2012.11.146
  * [8] M. Schmid, H.-G. Kneidinger, and C. Endisch, “Data-Driven Fault Diagnosis in Battery Systems through Cross-Cell Monitoring,” _IEEE Sensors Journal_ , pp. 1–1, 8 2020.
  * [9] J. Paparrizos and L. Gravano, “K-shape: Efficient and accurate clustering of time series,” in _Proceedings of the ACM SIGMOD International Conference on Management of Data_ , vol. 2015-May. Association for Computing Machinery, 5 2015, pp. 1855–1870.
  * [10] D. Bega, M. Gramaglia, M. Fiore, A. Banchs, and X. Costa-Perez, “DeepCog: Cognitive Network Management in Sliced 5G Networks with Deep Learning,” in _Proceedings - IEEE INFOCOM_ , vol. 2019-April, 2019.
  * [11] P. Gianniou, X. Liu, A. Heller, P. S. Nielsen, and C. Rode, “Clustering-based analysis for residential district heating data,” _Energy Conversion and Management_ , vol. 165, pp. 840–850, 6 2018.
  * [12] E. Calikus, S. Nowaczyk, A. Sant’Anna, H. Gadd, and S. Werner, “A data-driven approach for discovering heat load patterns in district heating,” _Applied Energy_ , vol. 252, p. 113409, 10 2019.
  * [13] S. N. Haider, Q. Zhao, and X. Li, “Data driven battery anomaly detection based on shape based clustering for the data centers class,” _Journal of Energy Storage_ , vol. 29, 6 2020.
  * [14] R. Tavenard, J. Faouzi, G. Vandewiele, F. Divo, G. Androz, C. Holtz, M. Payne, R. Yurchak, M. Rußwurm, K. Kolar, and E. Woods, “Tslearn, a machine learning toolkit for time series data,” _Journal of Machine Learning Research_ , vol. 21, 2020.