# Adaptive control for hindlimb locomotion in a simulated mouse through
temporal cerebellar learning

T. P. Jensen<EMAIL_ADDRESS>Automation and Control, Department of
Electrical Engineering, Technical University of DenmarkElektrovej, Building
326Kgs. Lyngby2800Denmark ,  S. Tata<EMAIL_ADDRESS>Biorobotics Laboratory, School of Engineering, École Polytechnique Fédérale de
LausanneStation 14Lausanne1015Switzerland ,  A. J. Ijspeert
<EMAIL_ADDRESS>Biorobotics Laboratory, School of Engineering, École
Polytechnique Fédérale de LausanneStation 14Lausanne1015Switzerland and  S.
Tolu<EMAIL_ADDRESS>Automation and Control, Department of Electrical
Engineering, Technical University of DenmarkElektrovej, Building 326Kgs.
Lyngby2800Denmark

###### Abstract.

Human beings and other vertebrates show remarkable performance and efficiency
in locomotion, but the functioning of their biological control systems for
locomotion is still only partially understood. The basic patterns and timing
for locomotion are provided by a central pattern generator (CPG) in the spinal
cord. The cerebellum is known to play an important role in adaptive
locomotion. Recent studies have given insights into the error signals
responsible for driving the cerebellar adaptation in locomotion. However, the
question of how the cerebellar output influences the gait remains unanswered.
We hypothesize that the cerebellar correction is applied to the pattern
formation part of the CPG. Here, a bio-inspired control system for adaptive
locomotion of the musculoskeletal system of the mouse is presented, where a
cerebellar-like module adapts the step time by using the double support
interlimb asymmetry as a temporal teaching signal. The control system is
tested on a simulated mouse in a split-belt treadmill setup similar to those
used in experiments with real mice. The results show adaptive locomotion
behavior in the interlimb parameters similar to that seen in humans and mice.
The control system adaptively decreases the double support asymmetry that
occurs due to environmental perturbations in the split-belt protocol.

Bio control, learning algorithms, adaptive locomotion, brain models, motor
control

††copyright: acmcopyright††conference: Neuro-inspired Computational Elements
Workshop; March 17–20, 2020; Heidelberg, Germany

##  1\. Introduction

The locomotion of vertebrates has been studied for a very long time in
research areas of biology and neurology, and it has also attracted interest
from roboticists because of the inherent ability of biological control systems
to adapt and being robust to disturbances within the environment. The
consensus in the neurological literature is that the spinal cord contains a
central pattern generator (CPG) responsible for generating the timing and
patterns of muscle activation signals for locomotion (Dietz, 2003; Kiehn,
2006; Takakusaki, 2013). Examples of robotic control systems with CPG s
include quadrupedal robots (Espinal et al., 2016), robot snakes (Wu and Ma,
2010), salamanders (Ijspeert et al., 2007) and bipedal robots (Matsubara et
al., 2006; Fujiki et al., 2013).

Within the brain, the cerebellum has a key role for adaptive locomotion
(Morton and Bastian, 2004). While the CPG provides the basic control patterns
to generate the gait, the cerebellum is required to adapt the patterns in case
of changes or perturbations in the environment (Morton and Bastian, 2006).
Evidence suggests (Malone et al., 2012; Mawase et al., 2013) that the
cerebellum works as an adaptive feed-forward controller for the locomotion
task since an after-effect is present when the perturbation is removed. The
fundamental research question is how adaptive locomotion is achieved in
vertebrates, and whether we can apply it to robots. In this work, we aim to
identify the error signal responsible for driving the adaptation, and to
locate where the correction is given by the cerebellum. It is debated whether
the cerebellum gives corrections in terms of a reference trajectory for the
limbs or rather with some higher objectives in mind. Considering that the
spinal CPG is a complex network of interneurons, with groups that are
responsible for rhythm generation, pattern generation, and motor output
respectively (Takakusaki, 2013), it is natural to speculate to which group the
cerebellar correction applies. In this paper, we hypothesize that for rodents
and humans in locomotion, the double support asymmetry is used as an error
signal in the cerebellum to give corrections to the pattern formation part of
the CPG.

Some researchers performed trials with humans on a split-belt that consist of
a treadmill with two belts that are individually controlled (Reisman et al.,
2005; Morton and Bastian, 2006; Malone et al., 2012), to analyze the adaptive
behavior. Morton describe the adaptation of gait parameters (Morton and
Bastian, 2006). The results showed two types of adaptation, described as
reactive and predictive. The reactive adaptation happened rapidly and returned
promptly to the baseline values in the post-adaptation period. The predictive
feed-forward adaptation occurred more slowly than the reactive one and it had
an after-effect in the post-adaptation period. Their findings for patients
with cerebellar damage showed no loss of the reactive adaptation. The
predictive adaptation was, however, much less (or absent) compared to healthy
patients. They identified the reactive adaptation as impacting the intralimb
parameters, stance time, and length. The predictive adaptation was seen in
interlimb parameters, step length, and double support time.

In a series of articles, it has been proposed that the locomotor cerebellar
adaptation in humans consists of a spatial and a temporal component.
Experiments during split-belt trials show no impact on a temporal adaptation,
when the subjects are distracted (Malone and Bastian, 2010). The spatial
adaptation, however, suffers from the distraction, which leads to the
conclusion that two different neural circuits are involved in the learning
process to reduce the step length asymmetry. These were shown to adapt at
different rates and in an experiment it was possible to circumvent the
adaptation of the temporal component (Malone et al., 2012). For the temporal
component, the authors proposed that the double support asymmetry is utilized
as a temporal error signal for adaptation in human locomotion. It has been
shown that cerebellar adaptations can be separated into spatial and temporal
contributions also in mice (Darmohray et al., 2019). Their results suggest
that the adaptation is primarily on the front limbs and that the hindlimbs
might simply adjust to match the patterns of the front limbs. It is argued
that the cerebellum adapts to reduce the step length asymmetry, which they
decompose into adaptations in spatial and temporal domains. The spatial
adaptation affects the asymmetry in the center of oscillation and the temporal
adaptation affects the double support asymmetry. Neither of those asymmetries
is reduced to zero in their experimental results, rather they are both reduced
to give almost no step length asymmetry. The cerebellar adaptations of humans
and mice carry many similarities, according to (Darmohray et al., 2019). This
hints that the structure of the involved cerebellar circuits might be very
similar for rodents and humans, even considering the differences between
bipedal and quadrupedal locomotion.

To study the biological control systems for locomotion, Fujiki et al.
presented a CPG model for controlling the hindlimbs in a simulated rat (Fujiki
et al., 2018). The model used a hip reflex for phase resetting in the neural
oscillators to achieve reactive adaptation in the split-belt protocol.
Cerebellar models are less common in locomotion control. In an earlier study
by Fujiki (Fujiki et al., 2015), a bipedal robot is controlled by means of a
CPG with a similar oscillator model, and a cerebellar-like model to adapt the
phase resetting with a foot contact sensor. The cerebellar model used a simple
equation for learning based on a phase prediction error. This is the only
previous example of a cerebellar-like model for adaptive locomotion in a
robot, where the cerebellar correction is not applied to the motor output. For
example, in (Massi et al., 2019) a cerebellar-like controller was used to give
joint torque corrections for locomotion of a quadrupedal robot, by utilizing a
reference limb trajectory to generate the error signals.

In this work, we sought to model the cerebellar temporal adaptation in mice
locomotion. To do so, a control architecture was designed by including a
modified version of the CPG by Fujiki et al. (Fujiki et al., 2018). A
cerebellar learning rule similar to the one described in (Fujiki et al., 2015)
was applied, however, the error signal was based on the theories of Malone et
al. (Malone et al., 2012). In contrast to (Fujiki et al., 2015), the
cerebellar output is applied as a correction signal to the pattern formation
part of the CPG. To investigate how adaptive locomotion is achieved for the
mouse, the musculoskeletal system was simulated using muscle models and
controlled by the developed bio-inspired control architecture.

The organization of the paper is as follows: an overview of the mouse hindlimb
simulation is given in section 2 together with a definition of the
experimental protocol and a presentation of the bio-inspired control system;
in section 3 the results are presented; finally, in section 4, the results are
discussed and compared with the ones from real mice and humans to highlight
the effect of the selected error signal on the adaptive locomotion behavior.

##  2\. Methods

###  2.1. Experimental setup

The experimental setup consists of a simulation of a mouse on a split-belt
treadmill in the Webots simulation environment (Michel, 2004). Only hindlimb
locomotion is considered. The hip, knee and ankle joints of each hindlimb are
controlled by a total of 8 muscles consisting of both monoarticular (spanning
across one joint) and biarticular (spanning across two joints) muscles. The
muscles that are included in the mouse model are shown in Figure 1.

The mouse model was based on a high-resolution 3-dimensional scan of a mouse
skeleton and the muscles are modeled using Geyer muscle equations (Geyer et
al., 2003; Geyer and Herr, 2010) with muscle parameters from (Charles et al.,
2016).

Figure 1. Illustration of the mouse hindlimb model, not showing the exact
attachments of the muscles. The muscles that have been modeled was selected
from the full set of muscles due to their relevance in forward locomotion. For
a more accurate picture see (Charles et al., 2016). PMA: psoas major, CF:
caudofemoralis, SM: semimembranosus, POP: popliteus, RF: rectus femoris, LG:
lateral gastrocnemius, SOL: soleus, TA: tibialis anterior.

The experimental split-belt protocol has been used extensively on real mice,
cats, and humans in the literature (Fujiki et al., 2018; D’Angelo et al.,
2014; Helm and Reisman, 2015). The setup consists of two separate treadmills,
one for each side of the body. By changing the velocity of one of the belts,
the gait is perturbed. An example of an experimental trial is illustrated in
Figure 2. Typically, it consists of a baseline period, where the belt
velocities are tied, a split period where the velocity of one of the belts is
increased, and finally a period of tied belt velocities. Here, the base belt
velocity is selected to
$6\text{\,}\mathrm{m}\mathrm{/}\mathrm{m}\mathrm{i}\mathrm{n}$ and speeds of
9, 10.2 and $12\text{\,}\mathrm{m}\mathrm{/}\mathrm{m}\mathrm{i}\mathrm{n}$
are used for the fast belt (1.5x, 1.7x and 2.0x speed ratio). This base
velocity was selected by matching belt speed to the CPG frequency and the step
length in a single belt experiment. To describe the adaptation, the split
period is divided into early and late adaptation periods, with the second
tied-belt period being divided similarly.

Figure 2. The split-belt protocol (2.0x speed ratio), showing the velocity of
each belt for one experiment.

Because the mouse is only walking on the hindlimbs, the spine is fixed in
space by a physics plugin in the simulator. This setup is similar to that in
(Fujiki et al., 2018), where the rats wore a harness, and the front limbs were
supported by a bar.

For the analysis, the parameters related to the gait are illustrated in Figure
3 and 4 for the spatial and temporal parameters, respectively. These are
parameters that can be measured in every stride. We define a stride as one
full locomotion cycle, consisting of a step on each limb with double-support
periods in between. For the analysis of the gait, the joint angles, toe
positions, and ground contact sensor data of the feet is logged. All of the
gait parameters are then derived from this data, using the mathematical
definitions found in Appendix A.

Figure 3. Illustration of the spatial gait parameters, stance length, and step
length. Blue dots denote the joints. Figure 4. Illustration showing the
temporal gait parameters, step time, and stance time. The horizontal direction
indicates time, increasing from left to right. A stride, consisting of a step
and stance period of each leg, is illustrated by the region between the two
dotted lines. The red and blue horizontal lines correspond to periods of
contact between the belts and the fast and slow limb, respectively, and the
empty areas between the lines of same color correspond to the swing time
periods. LO: time of lift-off, TD: touchdown, DS: double support.

###  2.2. Control architecture

The control architecture consists of a spinal CPG and a cerebellar-like module
for adaptation (see Fig. 5). The CPG is the same as described by Fujiki et al.
(Fujiki et al., 2018), but adapted and extended for this mouse model. It is
comprised of the rhythm generation and pattern formation layers, with the
latter generating the motor neuron activation patterns made of a series of
rectangular pulses. In this work, the periodic activation pattern was extended
to contain a fourth activation pulse. The addition of a fourth activation
pulse was inspired by kinematic data from real mice (Leblond et al., 2003).
Here the four gait stages are enumerated as flexion (F), contact (E1),
extension (E2), and push off (E3). Furthermore, the parameters of the CPG were
adjusted to this mouse model.

Figure 5. The bio-inspired control architecture used for adaptive locomotion
of the mouse hindlimbs. The spinal CPG provides timing and patterns of muscle
activation $u$ while the cerebellar model provides feed-forward corrections to
the step time for symmetric double support times. $\bm{\theta}$ is the vector
of all hindlimb joint angles, $\theta_{h}$ are the hip angles, and $c_{r}$ and
$c_{l}$ are the contact sensor values of the right and left hindlimb,
respectively.

In the CPG, the rhythm is generated by two oscillators, one for each hindlimb,
giving a phase $\phi\in[0,2\pi]$ for each. For the rhythm generator
parameters, the frequency was chosen to
$\omega=$19\text{\,}\mathrm{r}\mathrm{a}\mathrm{d}\mathrm{/}\mathrm{s}$$ and
the oscillator coupling strength to $K_{\phi}$ = 7.5. A phase-resetting
mechanism keeps the gait in sync with the belt. The phase of each oscillator
is reset when the corresponding hip angle reaches a threshold.

The motor neuron command is formed in the pattern formation layer. It consists
of four activation pulses, corresponding to each gait stage. From the phase
signal, pulses of unit magnitude are generated by the following equation:

(1)
$u_{j}^{\textnormal{CPG}}(\phi)=\begin{cases}1,&\phi_{j}^{\textnormal{start}}\leq\phi\leq\phi_{j}^{\textnormal{start}}+\Delta\phi_{j}\\\
0,&\textnormal{otherwise}\end{cases},j={1\dots 4},$

where $\phi_{j}^{\textnormal{start}}$ is the pulse phase onset,
$\Delta\phi_{j}$ is the duration, and $j$ denotes the activation pulse, E1,
E2, E3, or F.

The motor neuron command $u_{m}$ is then given as a weighted sum of the
pulses, with weighing coefficients $w_{m,j}$:

(2) $u_{m}=\sum_{j=1}^{4}w_{m,j}\cdot
u_{j}^{\textnormal{CPG}}(\phi),\;m=1\dots 8,$

where $m$ denotes each of the muscles. For more details on the CPG model, see
(Fujiki et al., 2018).

In the cerebellar microcircuit, the inferior olive (IO) gives the teaching or
error signal responsible for driving the adaptation. In the proposed
cerebellar-like module, the IO estimates the double support asymmetry from the
contact sensor, to give the temporal error signal $e_{t}$. The asymmetry is
given by the difference in the estimated double support times:

(3) $e_{t}=DS_{s}-DS_{f},$

where $DS_{s}$ and $DS_{f}$ are the slow and fast double support time,
respectively. Fast double support time corresponds to the duration of the
double support period which ends with liftoff of the limb on the fast belt,
and similarly for slow double support time (see Fig. 4). The double support
times are estimated from the contact sensor data using a state machine with
timers.

In the cerebellar microcircuit, the Purkinje cell (PC) is the central learning
component. In the proposed cerebellar-like model, the output from the PC is a
correction to the duration of the flexion activation pulse in the pattern
formation structure of the CPG:

(4) $\Delta\phi_{F}=\Delta\phi_{F}^{0}+y_{\textnormal{cerebellum}},$

where $y_{\textnormal{cerebellum}}$ is the correction and $\Delta\phi_{F}^{0}$
is the baseline value of the flexion pulse duration. This adjusts the duration
of the flexion pulse to adjust the step time, and was chosen to allow
correction of the double support time in accordance with Figure 4. The pattern
formation part of the CPG was selected as target for the correction because it
contains parameters related to the separate stages of the gait, thus the step
time is more easily influenced. As the error signal is only updated after each
step, a correction to the rhythm generation part would be constant in one
locomotion cycle, and thus affect all stages of the gait equally.

To ensure that the flexion (F) pulse and the early extension (E1) pulse have
the same phase spacing, the cerebellum also corrects the onset phase of the E1
pulse:

(5) $\phi_{E1}^{Start}=\phi_{E1}^{Start,0}+y_{\textnormal{cerebellum}},$

where $\phi_{E1}^{Start,0}$ is the baseline value for the onset phase of the
E1 pulse.

The cerebellar learning rule is a gradient-descent type with learning rate
$\alpha^{*}$ which is applied at the end of every period of double support.
The ($i+1$)’th output values for each limb is given by

(6) $y_{\textnormal{cerebellum,
s}}^{(i+1)}=y_{\textnormal{cerebellum}}^{(i)}-\alpha^{*}e_{t},$

and

(7) $y_{\textnormal{cerebellum,
f}}^{(i+1)}=y_{\textnormal{cerebellum}}^{(i)}+\alpha^{*}e_{t},$

for the slow and fast limb, respectively. The derivation of the signs for the
update rules can be found in Appendix B.

##  3\. Results

The gait on a single treadmill is shown in Figure 6 where the control
architecture has been applied without the cerebellar model.

Figure 6. Series of images showing the base gait on a single belt. The
locomotion direction is to the right, thus the belt moves to the left.

We now apply the full control architecture in the split-belt trial. Figure 7
shows that the cerebellar-like adaptation is able to reduce the double support
asymmetry in the split-belt protocol by making corrections to the pattern
formation part of the CPG. As the output signals are symmetric around zero,
only the output for the fast limb, $y_{\textnormal{cerebellum,f}}$, is shown.

(a) Error.

(b) Output.

Figure 7. The plots show the $e_{t}$ online error estimate and the
$y_{\textnormal{cerebellum,f}}$ output of the cerebellar-like model for the
fast limb in the split-belt protocol with a belt velocity ratio of 1.5. The
belt speed is perturbed at $t=10s$ as shown by the grayed out region.

The experimental protocol was carried out with and without cerebellar
learning. For the experiment with a velocity ratio of 1.5, the time-evolution
of the temporal and spatial parameters presented in Section II, throughout the
protocol, are compared for the mouse with and without cerebellar learning in
Figures 8 to 12.

It is seen in Figure 9(b) that the fast step time is increased as the
cerebellar output for that limb is increased (which is not surprising given
the derivation in Appendix B). However, the remainder of the gait parameters
also changes, as a result thereof. When the belt speed is perturbed, all of
the gait parameters change immediately in the early adaptation period. Without
cerebellar learning, the gait parameters still change due to the hip reflex,
which attempts to keep up to the faster belt speed, but they are not adapted
during the early- and late adaptation periods. This reflex effect happens
already at the first stride after the perturbation and does not change after
that. Without learning, there is a constant double support asymmetry (Fig.
12(a)) during the perturbation, which returns to zero immediately when the
perturbation is removed in post-adaptation.

The cerebellar adaptation does not only affect the step time, but also the
remaining gait parameters to a varying degree. While no effect is seen on the
stance length (Fig. 10) and stance time (Fig. 11) of the fast leg at all, the
stance length and time for the slow leg increases during the adaptation. When
the step time of the slow limb is decreased, the stance time automatically
increases because the total stride time is approximately constant. The total
stride time increases marginally when the effective oscillator frequency
increases slightly, due to the hip reflex. The reason for no increase in the
fast stance length could be that the hip reflex is activating, initiating the
lift-off of the limb.

The slow stance time is increased by a considerable amount (Fig. 11(b)), which
is working against the adaptation. To decrease the double support asymmetry,
the cerebellum is shortening the slow double support time and lengthening the
fast double support time. To shorten the slow double support time, the fast
step time is increased (see Fig. 4). But Figure 11(b) shows that the slow
stance time is increasing, which is working against decreasing the slow double
support time.

The step length (Fig. 8(b)) seems directly correlated to step time, as one
would expect; when the step time is increased due to the adaptation, the step
length increases as well.

With cerebellar learning, the double support asymmetry is reduced from the
early to late adaptation period, but when the perturbation is removed, a
double support asymmetry with the opposite sign arises in post-adaptation
(Fig. 7(a)). In late post-adaptation, this error has been reduced by the
learning rule. This after-effect happens because the cerebellar output has to
adapt back to the baseline value, as seen in Figure 7(b).

(a) Without cerebellum.

(b) With cerebellum.

Figure 8. Step length in the split-belt trial with a 1.5x speed ratio.

(a) Without cerebellum.

(b) With cerebellum.

Figure 9. Step time in the split-belt trial with a 1.5x speed ratio.

(a) Without cerebellum.

(b) With cerebellum.

Figure 10. Stance length in the split-belt trial with a 1.5x speed ratio.

(a) Without cerebellum.

(b) With cerebellum.

Figure 11. Stance time in the split-belt trial with a 1.5x speed ratio.

(a) Without cerebellum.

(b) With cerebellum.

Figure 12. Double support time in the split-belt trial with a 1.5x speed
ratio.

The double support adaptation is shown in Figure 12(b) and 13 for each of the
relative belt velocities (1.5x, 1.7x, and 2.0x). It can be seen that the
initial asymmetry is larger for higher velocity ratios, and for 2.0x there is
initially no period of fast double support. Still, in that case, the
cerebellar adaptation is able to reduce the asymmetry. When adapting, both
slow and fast double support times are changing, and the final late-adaptation
value is lower for higher velocity ratios. When the fast limb has to take
longer steps due to the faster moving belt, more time is spent in the swing
stage, and there is less time for double support. The adaptation time is
longer for higher velocity ratios, having adapted after $\approx$30,
$\approx$40, and $\approx$50 strides. The magnitude of the after-effect and
the time to de-adapt is also larger with a higher velocity ratio.

(a) 1.7x

(b) 2.0x

Figure 13. Double support time in the split-belt trial with two higher belt
velocity ratios.

##  4\. Discussion

This paper focuses on the input and output of the cerebellum rather than on
the evaluation of its learning mechanism. The cerebellar learning rule was
derived based on the internal mechanism of the CPG model, which enabled the
cerebellum to efficiently reduce the double support time asymmetry. This shows
that the assumption used to derive the learning rule (see Appendix B) is
coherent. The cerebellar-like module does not capture the complexity of the
real cerebellum, however, it is used as a phenomenological model to determine
the effect of applying the selected input (double support asymmetry) and
output (pattern formation part of the CPG) for the cerebellar-like model. We
will now compare these results with the ones from the literature, to determine
whether the chosen input and output gives rise to similar adaptive behavior.

Fujiki et al. (Fujiki et al., 2018) ran simulations of hindlimb locomotion of
a rat and compared their results with data from real rats, to determine the
importance of phase resetting by a hip reflex in a split-belt trial. The CPG
model implemented here is a further development of the one from (Fujiki et
al., 2018). By comparing our results with theirs, without the cerebellar
model, it is clear that the effects of the reflex are qualitatively the same,
in that the gait parameters change in same direction due to the split-belt
perturbation. This hip reflex, or one with a similar effect, is very important
for a stable gait when the belt speed is perturbed.

The comparison of mice walking on hindlimbs and human bipedal locomotion in
split-belt trials gives many interesting insights. After all, the only
difference is that mice do not have to keep balance, where humans do (although
in most experiments they are assisted by a handrail). The results in (Morton
and Bastian, 2006) correspond very well with the results we obtained by
applying the simplistic control architecture on the musculoskeletal model. For
example, we also found that the faster reactive adaptations were provided by
the reflexes. The correction that is seen in the stance length (Fig. 10) is
also longer on the fast limb, and slightly shorter on the slow limb in their
case, however, the slow stance time does not show any adaptation in the
measurements made by Morton et al. (Morton and Bastian, 2006). The reactive
change in the step length in early adaptation is in the same direction, and
the adaptation decreases the asymmetry. For the healthy subjects, the step
length asymmetry was almost completely removed through the adaptation. The
magnitudes of the gait parameters are not expected to be comparable between
mice and humans. In (Reisman et al., 2005), healthy human subjects are
subjected to the split-belt, where the results show that the adaptations take
more time, for a larger speed ratio. The same is the case here, as seen in
Figure 12(b) and 13.

In (Darmohray et al., 2019) it was shown that cerebellar adaptations in mice
can be separated into spatial and temporal contributions, with the temporal
contribution adapting on a shorter time scale. Also, it was suggested that the
adaptation is primarily learned for the forelimbs. Their results showed that
the double support asymmetry is reduced in the split-belt period, albeit only
for the front limbs. In this paper, only the temporal component was modeled
and the simulated mouse is only walking on the hindlimbs. Thus, the results
here cannot be compared directly with the ones from the aforementioned paper,
however, the adaptive behavior carries a strong resemblance. When comparing
the results presented here for the hindlimbs, with their forelimb results, a
lot of similarities can be seen. The reactive change in stance length is in
their results almost completely symmetrical around the baseline value, which
is not the case here, however, both cases show no adaption of stance length in
the split period. At the same time, the step length shows very similar
adaptation, but the reactive change is not the same. The differences are not
surprising when comparing fore- and hindlimbs. In general, the comparison with
results from the literature shows a good agreement in that interlimb
parameters (step length and duration) are predominantly affected by the
cerebellar adaptation while the intralimb parameters (stance length and
duration) are almost unaffected.

From the results, it is clear that applying the double support asymmetry as an
error signal for adapting step timing in the pattern formation part of the CPG
does result in adaptive behavior that is similar to that observed for mice and
humans. It is unclear whether the real cerebellum does use double support
asymmetry as an error signal, or whether a reduction of the double support
asymmetry is only an effect of optimization of a higher objective, for example
to keep balance (Houk and Miller, 2001), or to save energy (Finley et al.,
2013).

##  5\. Conclusion

In this paper, we presented a bio-inspired adaptive control architecture for
hindlimb locomotion in the mouse to test the hypothesis that the cerebellum
gives corrections to the pattern formation part of the CPG, using the double
support asymmetry as a temporal teaching signal. In our simulated split-belt
experiment the musculoskeletal system of a mouse is able to adapt to
perturbations coming from the unequal belt velocities, by adjusting the step
timing. The comparison with results from real mice and humans in the
literature shows that our control architecture captures many of the same
adaptive features. This contributes to validate the double support asymmetry
as the teaching signal used for gait adaptation in humans and mice. At the
same time, it gives confidence to the cerebellar output being applied to
different parts of the CPG network, and not only as a correction to the motor
neuron output. However, it is clear that in order to capture all of the
dynamics, a model of spatial adaptation is needed. For this reason, we
hypothesize that cerebellar internal models of limb dynamics can be used to
accurately achieve desired foot positions to reduce a spatial error signal,
such as the asymmetry in the center of oscillation.

In a future work, it will be of interest to further develop the cerebellar-
like model for adaptive quadrupedal locomotion in the mouse, and to
investigate if both spatial and temporal errors will allow the musculoskeletal
system to keep the balance.

## References

  * (1)
  * Charles et al. (2016) James P Charles, Ornella Cappellari, Andrew J Spence, Dominic J Wells, and John R Hutchinson. 2016\.  Muscle moment arms and sensitivity analysis of a mouse hindlimb musculoskeletal model.  Journal of anatomy 229, 4 (2016), 514–535.
  * D’Angelo et al. (2014) Giuseppe D’Angelo, Yann Thibaudier, Alessandro Telonio, Marie-France Hurteau, Victoria Kuczynski, Charline Dambreville, and Alain Frigon. 2014.  Modulation of phase durations, phase variations, and temporal coordination of the four limbs during quadrupedal split-belt locomotion in intact adult cats.  Journal of Neurophysiology 112, 8 (2014), 1825–1837.  https://doi.org/10.1152/jn.00160.2014
  * Darmohray et al. (2019) Dana M Darmohray, Jovin R Jacobs, Hugo G Marques, and Megan R Carey. 2019.  Spatial and Temporal Locomotor Learning in Mouse Cerebellum.  Neuron 102 (2019), 217–231.e4.  https://doi.org/10.1016/j.neuron.2019.01.038
  * Dietz (2003) V. Dietz. 2003.  Spinal cord pattern generators for locomotion.  Clinical Neurophysiology 114, 8 (2003), 1379–1389.  https://doi.org/10.1016/S1388-2457(03)00120-2
  * Espinal et al. (2016) A. Espinal, H. Rostro-Gonzalez, M. Carpio, E. I. Guerra-Hernandez, M. Ornelas-Rodriguez, H. J. Puga-Soberanes, M. A. Sotelo-Figueroa, and P. Melin. 2016\.  Quadrupedal Robot Locomotion: A Biologically Inspired Approach and Its Hardware Implementation.  Computational Intelligence and Neuroscience 2016 (2016), 1–13.  https://doi.org/10.1155/2016/5615618
  * Finley et al. (2013) James M. Finley, Amy J. Bastian, and Jinger S. Gottschall. 2013\.  Learning to be economical: the energy cost of walking tracks motor adaptation.  The Journal of Physiology 591, 4 (feb 2013), 1081–1095.  https://doi.org/10.1113/jphysiol.2012.245506
  * Fujiki et al. (2018) Soichiro Fujiki, Shinya Aoi, Tetsuro Funato, Yota Sato, Kazuo Tsuchiya, and Dai Yanagihara. 2018\.  Adaptive hindlimb split-belt treadmill walking in rats by controlling basic muscle activation patterns via phase resetting.  Scientific Reports 8, 1 (2018), 17341.  https://doi.org/10.1038/s41598-018-35714-8
  * Fujiki et al. (2015) Soichiro Fujiki, Shinya Aoi, Tetsuro Funato, Nozomi Tomita, Kei Senda, and Kazuo Tsuchiya. 2015\.  Adaptation mechanism of interlimb coordination in human split-belt treadmill walking through learning of foot contact timing: a robotics study.  Journal of The Royal Society Interface 12, 110 (sep 2015), 20150542.  https://doi.org/10.1098/rsif.2015.0542
  * Fujiki et al. (2013) Soichiro Fujiki, Shinya Aoi, Tsuyoshi Yamashita, Tetsuro Funato, Nozomi Tomita, Kei Senda, and Kazuo Tsuchiya. 2013.  Adaptive splitbelt treadmill walking of a biped robot using nonlinear oscillators with phase resetting.  Autonomous Robots 35, 1 (2013), 15–26.  https://doi.org/10.1007/s10514-013-9331-6
  * Geyer and Herr (2010) Hartmut Geyer and Hugh Herr. 2010.  A Muscle-Reflex Model That Encodes Principles of Legged Mechanics Produces Human Walking Dynamics and Muscle Activities.  IEEE Transactions on Neural Systems and Rehabilitation Engineering 18, 3 (jun 2010), 263–273.  https://doi.org/10.1109/TNSRE.2010.2047592
  * Geyer et al. (2003) Hartmut Geyer, Andre Seyfarth, and Reinhard Blickhan. 2003\.  Positive force feedback in bouncing gaits?  Proceedings of the Royal Society of London. Series B: Biological Sciences 270, 1529 (oct 2003), 2173–2183.  https://doi.org/10.1098/rspb.2003.2454
  * Helm and Reisman (2015) Erin E. Helm and Darcy S. Reisman. 2015.  The Split-Belt Walking Paradigm.  Physical Medicine and Rehabilitation Clinics of North America 26, 4 (nov 2015), 703–713.  https://doi.org/10.1016/j.pmr.2015.06.010
  * Houk and Miller (2001) James C Houk and Lee E Miller. 2001.  Cerebellum: Movement Regulation and Cognitive Functions.  In Encyclopedia of Life Sciences. John Wiley & Sons, Ltd, Chichester.  https://doi.org/10.1038/npg.els.0000036
  * Ijspeert et al. (2007) Auke Jan Ijspeert, Alessandro Crespi, Dimitri Ryczko, and J.-M. Cabelguen. 2007.  From Swimming to Walking with a Salamander Robot Driven by a Spinal Cord Model.  Science 315, 5817 (2007), 1416–1420.  https://doi.org/10.1126/science.1138353
  * Kiehn (2006) Ole Kiehn. 2006\.  Locomotor Circuits in the Mammalian Spinal Cord.  Annual Review of Neuroscience 29, 1 (2006), 279–306.  https://doi.org/10.1146/annurev.neuro.29.051605.112910
  * Leblond et al. (2003) Hugues Leblond, Marion L’Espérance, Didier Orsal, and Serge Rossignol. 2003.  Treadmill Locomotion in the Intact and Spinal Mouse.  The Journal of Neuroscience 23, 36 (dec 2003), 11411–11419.  https://doi.org/10.1523/JNEUROSCI.23-36-11411.2003
  * Malone and Bastian (2010) Laura A. Malone and Amy J. Bastian. 2010.  Thinking About Walking: Effects of Conscious Correction Versus Distraction on Locomotor Adaptation.  Journal of Neurophysiology 103, 4 (2010), 1954–1962.  https://doi.org/10.1152/jn.00832.2009
  * Malone et al. (2012) Laura A. Malone, Amy J. Bastian, and Gelsy Torres-Oviedo. 2012\.  How does the motor system correct for errors in time and space during locomotor adaptation?  Journal of Neurophysiology 108, 2 (2012), 672–683.  https://doi.org/10.1152/jn.00391.2011
  * Massi et al. (2019) Elisa Massi, Lorenzo Vannucci, Ugo Albanese, Marie Claire Capolei, Alexander Vandesompele, Gabriel Urbain, Angelo Maria Sabatini, Joni Dambre, Cecilia Laschi, Silvia Tolu, and Egidio Falotico. 2019.  Combining Evolutionary and Adaptive Control Strategies for Quadruped Robotic Locomotion.  Frontiers in Neurorobotics 13, August (2019), 1–19.  https://doi.org/10.3389/fnbot.2019.00071
  * Matsubara et al. (2006) Takamitsu Matsubara, Jun Morimoto, Jun Nakanishi, Masa-aki Sato, and Kenji Doya. 2006.  Learning CPG-based biped locomotion with a policy gradient method.  Robotics and Autonomous Systems 54, 11 (nov 2006), 911–920.  https://doi.org/10.1016/j.robot.2006.05.012
  * Mawase et al. (2013) F. Mawase, T. Haizler, S. Bar-Haim, and A. Karniel. 2013\.  Kinetic adaptation during locomotion on a split-belt treadmill.  Journal of Neurophysiology 109, 8 (2013), 2216–2227.  https://doi.org/10.1152/jn.00938.2012
  * Michel (2004) Olivier Michel. 2004\.  Webots: professional mobile robot simulation.  International Journal of Advanced Robotic Systems 1, 1 (2004), 39–42.  https://doi.org/10.5772/5618 arXiv:cs/0412052
  * Morton and Bastian (2004) Susanne M. Morton and Amy J. Bastian. 2004.  Cerebellar Control of Balance and Locomotion.  The Neuroscientist 10, 3 (jun 2004), 247–259.  https://doi.org/10.1177/1073858404263517
  * Morton and Bastian (2006) S. M. Morton and Amy J. Bastian. 2006.  Cerebellar Contributions to Locomotor Adaptations during Splitbelt Treadmill Walking.  Journal of Neuroscience 26, 36 (2006), 9107–9116.  https://doi.org/10.1523/JNEUROSCI.2622-06.2006
  * Reisman et al. (2005) Darcy S. Reisman, Hannah J. Block, and Amy J. Bastian. 2005\.  Interlimb Coordination During Locomotion: What Can be Adapted and Stored?  Journal of Neurophysiology 94, 4 (oct 2005), 2403–2415.  https://doi.org/10.1152/jn.00089.2005
  * Takakusaki (2013) Kaoru Takakusaki. 2013\.  Neurophysiology of gait: From the spinal cord to the frontal lobe.  Movement Disorders 28, 11 (2013), 1483–1491.  https://doi.org/10.1002/mds.25669
  * Wu and Ma (2010) Xiaodong Wu and Shugen Ma. 2010.  CPG-based control of serpentine locomotion of a snake-like robot.  Mechatronics 20, 2 (mar 2010), 326–334.  https://doi.org/10.1016/j.mechatronics.2010.01.006

##  Appendix A Gait parameters

All gait parameters are defined in terms of the slow or fast limb, where the
slow limb refers to the limb on the slow belt, and vice versa for the fast
limb.

We define by $x_{s}(t)$ and $x_{f}(t)$ the horizontal toe position of the slow
and fast limb with respect to time. The stance length is the intralimb
distance between the toe positions of touchdown and liftoff. Given as, for the
fast limb:

(8)
$s_{\textnormal{stance},f}=x_{f}(t_{\textnormal{touchdown},f})-x_{f}(t_{\textnormal{liftoff},f})$

and for the slow limb:

(9)
$s_{\textnormal{stance},s}=x_{s}(t_{\textnormal{touchdown},s})-x_{s}(t_{\textnormal{liftoff},s})$

The step length is the interlimb distance between the toes, at touchdown of
one of the limbs:

(10)
$s_{\textnormal{step},f}=x_{f}(t_{\textnormal{touchdown},f})-x_{s}(t_{\textnormal{touchdown},f})$
(11)
$s_{\textnormal{step},s}=x_{s}(t_{\textnormal{touchdown},s})-x_{f}(t_{\textnormal{touchdown},s}).$

Similarly, the step time is an interlimb parameter:

(12)
$t_{\textnormal{step},s}=t_{\textnormal{touchdown},s}-t_{\textnormal{touchdown},f}$
(13)
$t_{\textnormal{step},f}=t_{\textnormal{touchdown},f}-t_{\textnormal{touchdown},s},$

and the stance time is an intralimb parameter:

(14)
$t_{\textnormal{stance},f}=t_{\textnormal{liftoff},f}-t_{\textnormal{touchdown},f}$
(15)
$t_{\textnormal{stance},s}=t_{\textnormal{liftoff},s}-t_{\textnormal{touchdown},s}.$

Finally, we define the double support time, as the time when both limbs have
ground contact. We define the double support time for both limbs, as shown by
the grayed out area in Figure 4:

(16) $DS_{f}=t_{\textnormal{stance},f}-t_{\textnormal{step},s}\;,\hskip
8.5359ptDS_{s}=t_{\textnormal{stance},s}-t_{\textnormal{step},f}.$

##  Appendix B Derivation of the learning rule

If the loss, or evaluation function, is defined as:

(17) $V=\frac{1}{2}e_{t}^{2}.$

Then an error-based update rule can be derived for the cerebellum:

(18)
$y_{\textnormal{cerebellum}}^{(i+1)}=y_{\textnormal{cerebellum}}^{(i)}-\alpha\frac{\partial
V}{\partial y_{\textnormal{cerebellum}}}$

If we assume the step time to be a linear function of the phase duration, with
a positive correlation constant $k_{corr}$:

(19)
$t_{\textnormal{step}}=t_{\textnormal{step}}^{0}+k_{corr}(\Delta\phi_{F}-\Delta\phi_{F}^{0}),$

where $t_{\textnormal{step}}^{0}$ is the base step time given when
$\Delta\phi_{F}=\Delta\phi_{F}^{0}$, then the gradient will be:

(20) $\displaystyle\frac{\partial V}{\partial y_{\textnormal{cerebellum}}}$
$\displaystyle=\frac{\partial V}{\partial e_{t}}\frac{\partial e_{t}}{\partial
t_{\textnormal{step}}}\frac{\partial
t_{\textnormal{step}}}{\partial\Delta\phi_{F}}\frac{\partial\Delta\phi_{F}}{\partial
y_{\textnormal{cerebellum}}}$ $\displaystyle=e_{t}\cdot\frac{\partial
e_{t}}{\partial t_{\textnormal{step}}}\cdot k_{corr}\cdot 1$

The remaining derivative is different for the slow and the fast limb. For the
slow limb:

(21) $\frac{\partial e_{t}}{\partial t_{\textnormal{step},s}}=\frac{\partial
e_{t}}{\partial DS_{f}}\frac{\partial DS_{f}}{\partial
t_{\textnormal{step},s}}=(-1)\cdot(-1)=1,$

and for the fast limb:

(22) $\frac{\partial e_{t}}{\partial t_{\textnormal{step},f}}=\frac{\partial
e_{t}}{\partial DS_{s}}\frac{\partial DS_{s}}{\partial
t_{\textnormal{step},f}}=1\cdot(-1)=-1.$

The final cerebellar adaptation rules are then

(23) $y_{\textnormal{cerebellum,
s}}^{(i+1)}=y_{\textnormal{cerebellum}}^{(i)}-\alpha^{*}e_{t},$

and

(24) $y_{\textnormal{cerebellum,
f}}^{(i+1)}=y_{\textnormal{cerebellum}}^{(i)}+\alpha^{*}e_{t},$

where $\alpha^{*}=\alpha k_{corr}$ is the learning rate. This is the only
parameter that has to be chosen for the algorithm.


aainstitutetext: Authors affiliated with an institute covered by a cooperation
agreement with CERNbbinstitutetext: INFN, Sezione di Genova, 16147 Genova,
Italiaccinstitutetext: CERN, European Organization for Nuclear Research,
CH-1211 Geneva, Switzerlandddinstitutetext: ETH Zürich, Institute for Particle
Physics and Astrophysics, CH-8093 Zürich, Switzerlandeeinstitutetext:
Università degli Studi di Genova, 16126 Genova, Italiaffinstitutetext: UCL
Departement of Physics and Astronomy, University College London, Gower St.
London WC1E 6BT, United Kingdomgginstitutetext: Authors affiliated with an
international laboratory covered by a cooperation agreement with
CERNhhinstitutetext: Physics Department, University of Patras, 265 04 Patras,
Greeceiiinstitutetext: Universität Bonn, Helmholtz-Institut für Strahlen-und
Kernphysik, 53115 Bonn, Germanyllinstitutetext: Center for Theoretical and
Experimental Particle Physics, Facultad de Ciencias Exactas, Universidad
Andres Bello, Fernandez Concha 700, Santiago, Chilemminstitutetext: Millennium
Institute for Subatomic Physics at High-Energy Frontier (SAPHIR), Fernandez
Concha 700, Santiago, Chilenninstitutetext: Universidad Técnica Federico Santa
María and CCTVal, 2390123 Valparaíso, Chileooinstitutetext: Instituto de
Fisica Corpuscular (CSIC/UV), Carrer del Catedratic Jose Beltran Martinez, 2,
46980 Paterna, Valencia, Spainppinstitutetext: Department of Physics and
Astronomy, York University, Toronto, ON, Canadaqqinstitutetext: Departamento
de Fisica, Facultad de Ciencias, Universidad de La Serena, Avenida Cisternas
1200, La Serena, Chile

# First constraints on the $L_{\mu}-L_{\tau}$ explanation of the muon g-2
anomaly from NA64-$e$ at CERN

Yu. M. Andreev  0000-0002-7397-9665 b A. Antonov  0000-0003-1238-5158 c D.
Banerjee  0000-0003-0531-1679 d B. Banto Oberhauser  0009-0006-4795-1008 c J.
Bernhard  0000-0001-9256-971X b,e P. Bisio  /0009-0006-8677-7495 b A.
Celentano  0000-0002-7104-2983 c N. Charitonidis  0000-0001-9506-1022 f D.
Cooke d P. Crivelli  0000-0001-5430-9394 d E. Depero  0000-0003-2239-1746 a A.
V. Dermenev  0000-0001-5619-376X a S. V. Donskov  0000-0002-3988-7687 a R. R.
Dusaev  0000-0002-6147-8038 g T. Enik  0000-0002-2761-9730 g V. N. Frolov h A.
Gardikiotis  0000-0002-4435-2695 a S. N. Gninenko  0000-0001-6495-7619 i M.
Hösgen a V. A. Kachanov  0000-0002-3062-010X g Y. Kambar  0009-0000-9185-2353
a A. E. Karneyeu  0000-0001-9983-1004 g G. Kekelidze  0000-0002-5393-9199 i B.
Ketzer  0000-0002-3493-3891 a D. V. Kirpichnikov  0000-0002-7177-077X a M. M.
Kirsanov  0000-0002-8879-6538 a V. N. Kolosov g S. V. Gertsenberger
0009-0006-1640-9443 c S. Girod g E. A. Kasianova a,g V. A. Kramarenko
0000-0002-8625-5586 a L. V. Kravchuk  0000-0001-8631-4200 a,g N. V. Krasnikov
0000-0002-8717-6492 l,m S. V. Kuleshov  0000-0002-3065-326X a,m,n V. E.
Lyubovitskij  0000-0001-7467-572X g V. Lysan  0009-0004-1795-1651 b A. Marini
0000-0002-6778-2161 b L. Marsicano  0000-0002-8931-7498 g V. A. Matveev
0000-0002-2745-5908 m R. Mena Fredes m,n R. Mena Yanssen o L. Molina Bueno
0000-0001-9720-9764 d M. Mongillo  0009-0000-7331-4076 g D. V. Peshekhonov
0009-0008-9018-5884 a V. A. Polyakov  0000-0001-5989-0990 p B. Radics
0000-0002-8978-1725 g K. Salamatin  0000-0001-6287-8685 a V. D. Samoylenko d
H. Sieber  0000-0003-1476-4258 a D. Shchukin  0009-0007-5508-3615 m,q O. Soto
a V. O. Tikhomirov  0000-0002-9634-0581 a I. Tlisova  0000-0003-1552-2015 a A.
N. Toropin  0000-0002-2106-4041 o M. Tuzi  0009-0000-6276-1401 l P. Ulloa
0000-0002-0789-7581 a,g P. V. Volkov  0000-0002-7668-3691 a V. Yu. Volkov
0009-0005-3500-5121 a I. V. Voronchikhin  0000-0003-3037-636X l,m J. Zamora-
Saá  0000-0002-5030-7516 g A. S. Zhevlakov  0000-0002-7775-5917
<EMAIL_ADDRESS>
###### Abstract

The inclusion of an additional $U(1)$ gauge $L_{\mu}-L_{\tau}$ symmetry would
release the tension between the measured and the predicted value of the
anomalous muon magnetic moment: this paradigm assumes the existence of a new,
light $Z^{\prime}$ vector boson, with dominant coupling to $\mu$ and $\tau$
leptons and interacting with electrons via a loop mechanism. The
$L_{\mu}-L_{\tau}$ model can also explain the Dark Matter relic abundance, by
assuming that the ${Z^{\prime}}$ boson acts as a “portal” to a new Dark Sector
of particles in Nature, not charged under known interactions. In this work we
present the results of the ${Z^{\prime}}$ search performed by the NA64-$e$
experiment at CERN SPS, that collected $\sim 9\times 10^{11}$ 100 GeV
electrons impinging on an active thick target. Despite the suppressed
${Z^{\prime}}$ production yield with an electron beam, NA64-$e$ provides the
first accelerator-based results excluding the $g-2$ preferred band of the
${Z^{\prime}}$ parameter space in the 1 keV $<m_{Z^{\prime}}\lesssim 2$ MeV
range, in complementarity with the limits recently obtained by the NA64-$\mu$
experiment with a muon beam.

††arxiv: 2404.06982

##  1 Introduction

The Standard Model (SM) of fundamental interactions is one the greatest
successes of particle physics, explaining many phenomena at different energy
scales. Despite these results, the SM needs to be extended to account for
several experimentally observed anomalies or effects, currently not described
by the model. A remarkable example is provided by the Dark Matter (DM)
particle content puzzle. Nowadays, the existence of DM is confirmed by
multiple astrophysical and cosmological observations, but the SM does not
include any viable DM particle candidate Bertone  _et al._ (2005); Arcadi  _et
al._ (2018); Roszkowski  _et al._ (2018). This calls for an extension of the
SM, with new fields and forces not yet experimentally observed Zurek (2024).
Another significant example is provided by the measurement of the anomalous
muon magnetic moment $a_{\mu}\equiv(g_{\mu}-2)/2$. For this observable, the
most updated experimental average, mostly constrained by the latest Fermilab
Muon $g-2$ experiment, reads $a_{\mu}\mathrm{(Exp)}=116\,592\,059(22)\times
10^{-11}$ Aguillard  _et al._ (2023), to be compared with the latest
theoretical prediction from the Muon $g-2$ Theory Initiative,
$a_{\mu}\mathrm{(Teo)}=116\,591\,810(43)\times 10^{-11}$ Aoyama  _et al._
(2020). The leading order hadronic contribution for $a_{\mu}$ comes from the
cross section for the $e^{+}e^{-}\rightarrow$ hadrons process, measured by
many experiments. While new results from both the experimental side (e.g., the
latest CMD-3 result Ignatov  _et al._ (2023), still unpublished) and the
phenomenological one (e.g., a recent lattice calculation by the BMW
collaboration Borsanyi  _et al._ (2021)) tend to reduce the experiment-to-
theory discrepancy, still this difference motivates the investigation of new
physics scenarios with a preferred connection to SM second generation leptons.

One of the most attractive and simplest model explaining the $g-2$
discrepancy, the so-called $L_{\mu}-L_{\tau}$ scenario, predicts the existence
of a new $U(1)$ gauge boson, $Z^{\prime}$, that couples to second and third
generation leptons via a coupling $g_{Z^{\prime}}$ He  _et al._ (1991a, b),
and can also be connected to DM phenomenology. Although cosmological and
astrophysical considerations put stringent constraints on the model parameter
space, see e.g. Ref. Escudero  _et al._ (2019); Kamada and Yu (2015), so far
no laboratory-based experiment explored the preferred parameters region
releasing the $g-2$ tension in the sub-MeV ${Z^{\prime}}$ mass range.
Recently, NA64-$\mu$ at CERN set its first accelerator-based constraints using
a muon beam Andreev  _et al._ (2024), excluding the $L_{\mu}-L_{\tau}$
explanation of the $g-2$ anomaly for ${Z^{\prime}}$ masses heavier than $\sim
20$ MeV. In complementarity with this result, we present in this work the
first results from the electron-beam missing-energy experiment NA64-e, wich
exclude the $g-2$ band in the 1 keV $<m_{Z^{\prime}}\lesssim 2$ MeV mass
range. We note that, in principle, the NA64-$e$ experiment could be sensitive
also to a ${Z^{\prime}}$ lighter than 1 keV, but the production of such a
light particle would be affected by non trivial wave function coherence-loss
effects Demidov  _et al._ (2019), whose evaluation goes beyond the aim of this
work. These results are based on a $\sim$3-times larger accumulated statistics
if compared to the previous limits set by NA64-$e$ Andreev  _et al._ (2022),
resulting in a significant improvement of the explored parameter space.
Compared to the previous work, the analysis has been refined by performing a
full signal simulation via the DMG4 package Bondi  _et al._ (2021); Oberhauser
_et al._ (2024), including an effective description of the atomic electrons
motion in the target. In addition, the limits have been obtained with a
rigorous statistical approach, improving the uncertainty in the determination
of the systematic effects.

##  2 The $L_{\mu}-L_{\tau}$ model

This model introduces the following new lagrangian terms Bauer  _et al._
(2018):

$\displaystyle\mathcal{L}\subset$
$\displaystyle-\frac{1}{4}{Z^{\prime}}_{\mu\nu}{{Z^{\prime}}}^{\mu\nu}+\frac{1}{2}m^{2}_{Z^{\prime}}{Z^{\prime}}_{\mu}{Z^{\prime}}^{\mu}$
(1) $\displaystyle-
g_{Z^{\prime}}{Z^{\prime}}_{\mu}\left(\overline{\mu}\gamma^{\mu}\mu+\overline{\nu}_{\mu}\gamma^{\mu}P_{L}\nu_{\mu}-\overline{\tau}\gamma^{\mu}\tau-\overline{\nu}_{\tau}\gamma^{\mu}P_{L}\nu_{\tau}\right)\;\;,$

where
$Z^{\prime}_{\mu\nu}\equiv\partial_{\mu}{Z^{\prime}}_{\nu}-\partial_{\nu}{Z^{\prime}}_{\mu}$
is the ${Z^{\prime}}$ field strength, $m_{Z^{\prime}}$ is the ${Z^{\prime}}$
mass, and $P_{L}=(1-\gamma_{5})/2$. In this model, loop-order effects generate
an additional positive contribution to $a_{\mu}(\mathrm{Theo})$, bringing it
closer to $a_{\mu}(\mathrm{Exp})$. In particular, for values
$g_{Z^{\prime}}~{}\simeq~{}4.5\times 10^{-4}$ and $m_{Z^{\prime}}\ll m_{\mu}$,
the $Z^{\prime}$ contribution would actually solve the muon $g-2$ discrepancy.
Recently, results from the search of a ${Z^{\prime}}$ boson decaying to muons
have been reported by the BaBar Lees  _et al._ (2016), Belle Czank  _et al._
(2022), Belle-II Adachi  _et al._ (2024), and CMS Sirunyan  _et al._ (2019)
experiments, with null observations so far. The same result was reported by
Belle-II and BESIII from the search of an invisibly-decaying ${Z^{\prime}}$
Adachi  _et al._ (2023); Ablikim  _et al._ (2024).

The model can also be connected to the DM phenomenology in the context of Dark
Sector (DS) scenarios, by postulating the existence of a new ensemble of
particles in Nature, not charged under SM interactions. The lightest particle
of the DS (here denoted as $\chi$), if stable, can be the DM candidate.
Assuming, for illustration, a complex scalar nature for the $\chi$, the
Lagrangian is extended by the ${Z^{\prime}}-\chi$ interaction term ${\cal
L}_{D}=g_{D}{Z^{\prime}}_{\mu}J^{\mu}_{D}$, where
$J^{\mu}_{D}=\left(\chi^{*}i\partial^{\mu}\chi-\chi
i\partial^{\mu}\chi^{*}\right)$ is the DS current and
$g_{D}\equiv\sqrt{4\pi\alpha_{D}}$ is the coupling constant. In this picture,
considering the mass hierarchy $m_{Z^{\prime}}>m_{\chi}$, and far from the
resonance condition $m_{Z^{\prime}}=2m_{\chi}$, the presently observed DM
relic abundance is set by the velocity-averaged cross section for the process
$l\overline{l}\rightarrow Z^{\prime}\rightarrow\chi\chi$, where $l$ is a SM
lepton, if DM had a thermal origin in the early Universe. This results to a
preferred combination of the model parameters, related as Kahn  _et al._
(2018):

${g^{2}_{D}}g^{2}_{Z^{\prime}}\left(\frac{m_{\chi}}{m_{Z^{\prime}}}\right)^{4}\simeq
3\cdot 10^{-15}\left(\frac{m_{\chi}}{1\,\mathrm{MeV}}\right)^{2}\;\;.$ (2)

As discussed in Ref. Andreev  _et al._ (2022), even if the $L_{\mu}-L_{\tau}$
model does not explicitly include a tree-level coupling between electrons and
the $Z^{\prime}$, this is straightforwardly introduced via the one-loop level
kinetic mixing between the $Z^{\prime}$ and the photon, resulting to the
effective coupling $e\Pi(q^{2})$ dependent on the momentum squared $q^{2}$
carried by the ${Z^{\prime}}$ Araki  _et al._ (2017); Gninenko and Krasnikov
(2018); Zhang  _et al._ (2021)111The full analytic expression for $\Pi(q^{2})$
can be found in the Appendix of Ref. Andreev  _et al._ (2022).:

$\Pi(q^{2})=\frac{e\,g_{Z^{\prime}}}{2\pi^{2}}\int_{0}^{1}dx\,x(1-x)\ln\frac{m^{2}_{\tau}-x(1-x)q^{2}}{m^{2}_{\mu}-x(1-x)q^{2}}\;\;.$
(3)

In particular, for $q^{2}\ll m^{2}_{\mu}$, the loop function $\Pi(q^{2})$
approaches to
$\Pi(q^{2})\simeq\Pi(0)=\frac{e\,g_{Z^{\prime}}}{6\pi^{2}}\ln\frac{m_{\tau}}{m_{\mu}}\simeq
0.014\cdot g_{Z^{\prime}}$, and the phenomenology of the $L_{\mu}-L_{\tau}$
model is similar to that of a dark photon model Fabbrichesi  _et al._ (2020),
with the substitution $\varepsilon=0.0144\cdot g_{Z^{\prime}}$, with
$\varepsilon$ being the kinetic mixing parameter. At large $q^{2}$,
$|\Pi(q^{2})|$ grows, reaching a maximum value $\simeq 0.021g_{Z^{\prime}}$ at
$q^{2}=m^{2}_{\mu}$, and then smoothly goes to zero for $q^{2}\gg
m^{2}_{\mu}$.

Figure 1: Schematic representation of the NA64$-e$ detector, searching for an
invisibly decaying ${Z^{\prime}}$ produced by the interaction of 100 GeV
electrons with the material of the active ECAL target. See text for further
details.

##  3 The ${Z^{\prime}}$ search at NA64$-e$

Due to the loop-induced effective coupling, the ${Z^{\prime}}$ paradigm can be
explored at electron-beam experiments. In this work, we present the latest
results of the search for the ${Z^{\prime}}$ particle at the NA64$-e$
experiment at CERN SPS. In a collision of a high-energy electron beam with a
thick target, $Z^{\prime}$ particles can be produced by two main processes
involving the primary electron and the secondaries produced in the
electromagnetic cascade, namely the radiative emission
(${Z^{\prime}}$-strahlung) in the electromagnetic field of an atom by an
electron/positron, $e^{\pm}N\rightarrow e^{\pm}N{Z^{\prime}}$, and the
resonant annihilation of a secondary positron with an atomic electron,
$e^{+}e^{-}\rightarrow{Z^{\prime}}$. Depending on the mass and on the
couplings, the produced ${Z^{\prime}}$ will then decay to one of the allowed
dominant channels, i.e. $\mu^{+}\mu^{-}$, $\tau^{+}\tau^{-}$, $\nu\bar{\nu}$,
or $\chi\chi$, where the neutrino channel refers both to $\nu_{\mu}$ and
$\nu_{\tau}$. NA64$-e$ exploits the invisible decay channel
($\nu\overline{\nu}$ or $\chi\chi$) to investigate the ${Z^{\prime}}$
existence by relying on the missing-energy technique with a $E_{0}=$100 GeV
electron beam impinging on an active thick target. The signature for the
$Z^{\prime}$ production, followed by its invisible decay, is the observation
of events with a well identified impinging $e^{-}$ track, in coincidence with
a large missing energy $E_{miss}$, measured as the difference between the
nominal beam energy $E_{0}$ and the one deposited in the active target
$E_{ECAL}$, and no activity in the other downstream apparatus Gninenko  _et
al._ (2016a). A schematic representation of the NA64$-e$ detector assembly is
shown in Fig. 1. Incoming particles are detected by a set of three plastic
scintillator counters (S1, S2, S3) and two veto counters (V1, V2). A magnetic
spectrometer, consisting of tracking detectors (GEMs, MicroMegas, and Straw
Tubes) placed upstream and downstream two dipole magnets (MBPL 1, MBPL 2) with
a combined magnetic strength of approximately $7$ T$\cdot$m, is utilized for
measuring particle momentum, with a resolution $\delta p/p$ of about $1\%$
Banerjee  _et al._ (2018a). Particle identification is achieved by detecting
synchrotron radiation (SR) photons emitted by the electrons deflected by the
magnetic field through a compact Pb/Sc calorimeter (SRD) Depero  _et al._
(2017). The NA64 active target is a Pb/Sc calorimeter (ECAL) with a thickness
of 40 $X_{0}$, organized in a $5\times 6$ matrix of $3.82\times 3.82$ cm2
cells. Each cell has independent PMT readout and is longitudinally segmented
into a 4 $X_{0}$ pre-shower section (PRS) and a main section. Following the
ECAL, a hermetic Fe/Sc hadron calorimeter (HCAL) with three modules, totaling
approximately 21 $\lambda_{I}$ (nuclear interaction lengths), is installed to
detect secondary hadrons and muons produced in the ECAL and upstream beamline
elements. A fourth module is installed at zero degrees, to capture neutral
particles produced by the beam interacting with upstream materials. A high-
efficiency plastic scintillator counter (VETO) is placed between the ECAL and
HCAL to further reduce background signals. The trigger for the experiment
required the coincidence between signals from the scintillator counters, in
anti-coincidence with the veto detectors, paired with an in-time energy
deposition in the pre-shower $E_{PRS}\gtrsim 300$ MeV and in the whole ECAL
$E_{ECAL}\lesssim 90$ GeV Banerjee  _et al._ (2018b).

##  4 Data analysis

The results here presented are based on a total statistics of $(9.1\pm
0.5)\times 10^{11}$ electrons on target (EOT) accumulated in different runs by
NA64$-e$ during the 2016-2022 period, with beam intensity up to $\simeq
6\times 10^{6}$ electrons per SPS spill of 4.8 s. The data analysis was based
on the same dataset already scrutinized and unblinded for the recently
published invisibly-decaying dark photon search Andreev  _et al._ (2023a).
Therefore, supported by the fact that the phenomenology of an invisible-
decaying ${Z^{\prime}}$ is very similar to that of the ${A^{\prime}}$ search,
we decided to adopt the same reconstruction algorithms and event selection
criteria, resulting in zero observed events in the signal window Banerjee  _et
al._ (2018b). The criteria included the request for a well identified
impinging electron track with momentum $100\pm 10$ GeV, paired with an in-time
energy deposition in the SRD detector in the range $\simeq 1-100$ MeV. Events
with electro- and photo-nuclear interactions in the ECAL were suppressed by
requiring the compatibility of the measured energy deposition in the different
cells, in terms of longitudinal and transverse shape, with that expected from
an electron-induced electromagnetic shower Banerjee  _et al._ (2018b). In
addition, no activity must be observed in the VETO and in the HCAL detector;
for the latter, a 1 GeV energy threshold was set, just above the noise level.
Finally, we adopted the same value of the ECAL missing energy threshold
defining the signal box ($E_{ECAL}<E^{thr}_{ECAL}$) that was optimized
independently for each particular run period, accounting for small differences
in the detector response and in the observed background levels. Obtained
values were in the range $E^{thr}_{ECAL}\sim 47-50$ GeV. The most critical
background source originates from the interactions of the beam particles with
beamline materials upstream the ECAL. These interactions can produce secondary
hadrons emitted at large angles, outside the detector acceptance, while the
scattered low-energy electron may hit the ECAL with reduced energy deposition,
mimicking a signal Andreas  _et al._ (2013); Gninenko  _et al._ (2016b). To
suppress this background, a multiplicity cut on the measured number of hits in
the STRAW detectors was used. The residual yield was evaluated directly from
the data, considering events lying in the sideband region $E_{HCAL}\simeq 0$,
$E_{ECAL}<E_{0}$, and extrapolating the corresponding ECAL energy distribution
in the signal region through an exponential model. For the 2016-2018
(2021-2022) runs, the expected background is $0.4\pm 0.2\,(0.3\pm 0.2)$
events, where the error has been evaluated by varying the parameters of the
model within their fit uncertainty and repeating the procedure. Secondary
contributions arose from the in-flight decay of $\pi^{-}$ and $\mu^{-}$ beam
contaminants Andreev  _et al._ (2023b), estimated with Monte Carlo, as well as
from the punch-through of leading neutral hadrons $(n,K)$ produced in the
ECAL, valuated from a dedicated measurement Banerjee  _et al._ (2020).

##  5 Signal simulation

The expected yield $S$ of ${Z^{\prime}}$ events within the signal window for
given values of the model parameters was computed using the DMG4 software
Bondi  _et al._ (2021); Oberhauser  _et al._ (2024), integrated within the
full Geant4-based simulation of the NA64$-e$ detector Agostinelli  _et al._
(2003); Allison  _et al._ (2016). We considered two model variations: a
“vanilla” case in which $g_{D}=0$ and a “dark sector” one featuring non-zero
$Z^{\prime}-\chi$ DS coupling, focusing on the $1\,\mathrm{MeV}\lesssim
m_{Z^{\prime}}\lesssim 500\,\mathrm{MeV}$ mass range for $g_{Z^{\prime}}<0.1$.
For these parameters values, and in particular in the first scenario, the
${Z^{\prime}}$ width is very narrow Andreev  _et al._ (2022); this was
appropriately treated in the simulation, as discussed in the following. Both
the ${Z^{\prime}}$-strahlung and the resonant annihilation production
processes were implemented in DMG4, addressing the model peculiarities arising
from the loop-induced $e^{-}-{Z^{\prime}}$ coupling, as well the
${Z^{\prime}}$ width. We first considered the effects associated with the
$q^{2}$ dependency of the $\Pi$ function, starting from the
${Z^{\prime}}$-strahlung process. In the “dark sector” case and for the
${Z^{\prime}}\rightarrow\chi\chi$ decay channel, the cross-section expression
contains an average ${Z^{\prime}}-e^{-}$ kinetic mixing value $I_{|\Pi|^{2}}$
(see also Ref. Andreev  _et al._ (2022), Eq. A.5):

$I_{|\Pi|^{2}}=\int_{q^{2}_{min}}^{q^{2}_{max}}dq^{2}\frac{\left|\Pi(q^{2})\right|^{2}}{\pi}\frac{\sqrt{q^{2}}\,\Gamma_{Z^{\prime}}}{(q^{2}-m^{2}_{Z^{\prime}})^{2}+m^{2}_{Z^{\prime}}\Gamma^{2}_{Z^{\prime}}}\;,$
(4)

where $q^{2}_{min}=4m^{2}_{\chi}$ and $q^{2}_{max}\gg m^{2}_{Z^{\prime}}$ in
the NA64$-e$ kinematic regime. In Fig. 2, we compare $I_{|\Pi|^{2}}$ with the
on-shell value $|\Pi(m^{2}_{Z^{\prime}})|^{2}$, for $\alpha_{D}=0.1$ (black
curve) and $\alpha_{D}=0.02$ (red curve), by showing the relative difference
among the two as a function of the ${Z^{\prime}}$ mass222In the DS case, the
width $\Gamma_{Z^{\prime}}$ is dominated by the $\alpha_{D}$ contribution; we
therefore neglected its dependency on $g_{Z^{\prime}}$.. As expected, the
largest variation is seen for $m_{Z^{\prime}}\simeq 2m_{\mu}$, since here the
variation of the $\Pi$ function with respect to $q^{2}$ is larger. We observe
that, as a consequence of the narrow ${Z^{\prime}}$ width, the relative
difference is always smaller than $\simeq 1\%$. This justifies the use of an
on-shell approximation for the ${Z^{\prime}}-$strahlung cross section, i.e.
$I_{|\Pi|^{2}}\simeq|\Pi(m^{2}_{Z^{\prime}})|^{2}$. The same conclusion holds
for the “vanilla” scenario, given the even narrower ${Z^{\prime}}$ width.
Finally, for the $e^{+}e^{-}$ resonant production the full $q^{2}$ dependency
of the $\Pi$ function is naturally accounted for in the tree-level expression
for $\sigma(s)$, where $s$ is the $e^{+}e^{-}$ invariant mass squared, given
the identity $q^{2}=s$ – see also Ref. Andreev  _et al._ (2022), Eqs. (A.3)
and (A.4).

Figure 2: The relative difference between the $\Pi$ function on-shell value
$|\Pi|^{2}$ and the $q^{2}$-averaged one $I_{|\Pi|^{2}}$ as a function of the
${Z^{\prime}}$ mass, for different $\alpha_{D}$ values.

A second critical effect considered in the simulation is the modification of
the effective line shape for the $e^{+}e^{-}$ annihilation channel due to the
atomic electrons motion. This effect is well known in atomic physics, where it
induces, for example, a broadening of the photon energy spectrum in the
$e^{+}e^{-}\rightarrow\gamma\gamma$ reaction Ghosh  _et al._ (2000). In
particular, for in-flight annihilation atomic motions manifest as the
appearance of events with $E_{\gamma}>E_{max}$, where $E_{max}$ was the
maximum allowed photon energy if the atomic electron was at rest Hunt  _et
al._ (2001). To our knowledge, this effect has never been implemented in a
full Monte Carlo simulation of the resonant production of DS states, although
recently an analytical calculation was performed in Ref. Arias-Aragón  _et
al._ (2024) (see Appendix A). Indeed, for a given positron energy $E_{+}$ in
the reference frame of the target atom, the $e^{+}e^{-}$ invariant mass is
given by the expression $s=2m_{e}^{2}+2E_{+}(E_{-}-zP_{-})$, being $E_{-}$
($P_{-}$) the electron energy (momentum), and $z$ the cosine of the angle
between the two leptons. This implies that if the energy of a positron does
not align with the resonance mass when interacting with a stationary electron,
it can still do so when annihilating with an orbiting electron, and vice
versa.

Figure 3: Total atomic cross section for the process
$e^{+}e^{-}\rightarrow{Z^{\prime}}\rightarrow\nu\bar{\nu}$ on lead as a
function of the impinging positron energy, considering (red) or not (black)
the motion of atomic electrons.

To account for this, we computed an average cross section per atom, starting
from the tree-level expression for the annihilation cross section, given by
the formula

$\sigma=\frac{A(E_{+})}{(2m_{e}^{2}+2E_{+}(E_{-}-zP_{-})-m^{2}_{Z^{\prime}})^{2}+m^{2}_{Z^{\prime}}\Gamma_{Z^{\prime}}^{2}}\;\;,$
(5)

where $A(E_{+})$ is a smooth-varying function of the impinging positron
energy. We assumed a flat distribution for $z$ and integrated over it to
obtain an average value $\overline{\sigma}(E_{+},E_{-})$, depending solely on
the positron and electron energy. Then, to integrate over $E_{-}$, we
parameterized the kinetic energy distribution through a model inspired by the
virial theorem, $\overline{T}=E_{B}$, where $-E_{B}$ is the binding energy,
adopting an exponential ansatz $p(T_{-})=\frac{1}{E_{B}}\exp(-T_{-}/E_{B})$
Hunt  _et al._ (2001), and then we summed over the contributions from each
atomic shell Oberhauser  _et al._ (2024). For illustration, Fig. 3 shows the
cross section per atom as a function of the positron energy for the reaction
$e^{+}e^{-}\rightarrow{Z^{\prime}}\rightarrow\nu\bar{\nu}$ on lead (“vanilla”
case) when $m_{Z^{\prime}}=225$ MeV and $g_{A^{\prime}}=0.03$, comparing the
result obtained with (red) and without (black) the atomic effects. We observe
that, due to atomic effects, the cross section value at the peak is reduced by
many orders of magnitude, and simultaneously the shape is significantly
enlarged, keeping a constant value for the integral with respect to the
positron energy across the allowed kinematic range.

Figure 4: The new exclusion limits from the NA64$-e$ experiment in the
$(g_{Z^{\prime}},m_{Z^{\prime}})$ parameter space for the $L_{\mu}-L_{\tau}$
model, considering the “vanilla” scenario without any DS coupling (left) and
the full “dark sector one” (right), for $\alpha_{D}=0.1$ (continuous curve)
and $\alpha_{D}=0.02$ (dashed curve). The most competitive bounds reported by
other experiments are also shown (see text for further details). The green
band correspond to preferred combination of the model parameters that would
solve the muon $g-2$ tension. The blue dashed line defines the region excluded
by cosmological and astrophysical considerations. Finally, the black curves in
the right plot are the preferred combination of the parameters to explain the
observed dark matter relic density (so-called “thermal target”).

After the ${Z^{\prime}}$ yield evaluation via simulation, we estimated the
signal efficiency $\kappa$ directly from the experimental data, applying the
same selection cuts, except for the $E^{thr}_{ECAL}$ threshold, to pure
impinging electron events measured during calibration runs. During these runs,
no condition on the ECAL energy was included in the trigger, allowing to
collect events where the 100-GeV electrons fully deposit their energy in the
ECAL. The overall result reads $\kappa\simeq 50\%$, with a $\simeq 10\%$
variation between different runs mostly due to an increased beam intensity,
and thus to a larger pile-up probability, for 2016 and 2022 periods. To
account for the systematic uncertainty on the ECAL energy threshold, for each
$Z^{\prime}$ mass value the procedure was repeated multiple times by sampling
$E^{thr}_{ECAL}$ from a Gaussian PDF, and then re-evaluating $S$. The RMS
value of all results was taken as an estimate of $\Delta S$. Finally, an
additional 15$\%$ uncertainty on the absolute yield was also introduced to
account for the observed difference on the number of “di-muon” events between
data and Monte Carlo Andreev  _et al._ (2023a).

##  6 Results

Starting from the zero events measured in the signal window, the predicted
background level $B$, and the simulated ${Z^{\prime}}$ yield, we computed an
upper limit on the $g_{Z^{\prime}}$ coupling as a function of
$m_{Z^{\prime}}$, for the cases $\alpha_{D}=0$ (“vanilla”) and
$\alpha_{D}=0.02$ / $\alpha_{D}=0.1$ (“dark sector”). We used a frequentist
method, considering the $90\%$ Confidence Level (CL) of a one-sided profile-
likelihood test statistics Gross (2007). The likelihood model was built
assuming a Poisson PDF for the number of events in the signal region, with
mean value $\mu=S+B$. To handle the non-trivial dependency of $S$ on
$g_{Z^{\prime}}$, entering both in the cross-section multiplicative pre-factor
and in the ${Z^{\prime}}$ width Andreev  _et al._ (2022), we proceeded by
iteration, repeating at each time the computation of ${S}$ via Monte Carlo and
the extraction of the limit by adopting the $g_{Z^{\prime}}$ value returned
from the previous iteration. The first iteration was performed using the
$g_{Z^{\prime}}$ values reported in our previous work Andreev  _et al._
(2022). Convergence was observed already after three iterations.

The obtained results are shown in Fig. 4, reporting in red the $90\%$ CL limit
for the “vanilla” scenario (left plot) and for the “dark sector” one (right
plot), with $\alpha_{D}=0.1$ (continuous line) and $\alpha_{D}=0.02$ (dashed
line) and $m_{Z^{\prime}}/m_{\chi}=3$. In the same plots, we also show recent
exclusion limits from BaBar Lees  _et al._ (2016), Belle-II Adachi  _et al._
(2023) and BESIII Ablikim  _et al._ (2024), as well as the latest result
obtained by the NA64 collaboration through a dedicated missing-momentum
experiment with a muon beam, NA64-$\mu$Andreev  _et al._ (2024). Results from
previous neutrinos experiments (CCFR Altmannshofer  _et al._ (2014) and
Borexino Kamada and Yu (2015)) are also reported – as already underlined in
Ref. Gninenko and Gorbunov (2021), these should be taken with care, being them
a theoretical re-interpretation of the original experimental data. The blue
dashed line defines the region excluded by cosmological and astrophysical
considerations Escudero  _et al._ (2019); this has been reported only in the
“vanilla” case, since no connection with DM has been considered in Escudero
_et al._ (2019). In both figures, the green band correspond to the preferred
combination of the ${Z^{\prime}}$ parameters that would solve the muon $g-2$
tension. Finally, in the “dark sector” scenario, the black lines report the
predicted value of $g_{Z^{\prime}}$ as a function of the ${Z^{\prime}}$ mass
from Eq. 2. NA64-$e$ is the first accelerator-based experiment to
unambiguously probe the $g-2$ band in the 1 keV $<m_{Z^{\prime}}\lesssim 2$
MeV mass range, confirming the re-analysis of the Borexino results. We notice
that the limits set by NA64-$e$ tend to a constant value for a ${Z^{\prime}}$
lighter than $\sim 0.1$ MeV. This is due to the interplay between the cross
section of the ${Z^{\prime}}$ radiative production and the signal selection
cuts in NA64-$e$. Due to the signal window definition, the experiment is
sensitive only to ${Z^{\prime}}$s produced with a significant fraction of the
initial $e^{+}/e^{-}$ energy $x\geq 0.5$. In this kinematic range, if
$m_{Z^{\prime}}\ll m_{e}$, the differential cross section for the
${Z^{\prime}}$ radiative production is approximately independent from
$m_{Z^{\prime}}$ Gninenko  _et al._ (2018), resulting in a flat limit curve.

##  7 Conclusions

In conclusion, we performed the analysis of the data collected by NA64-$e$ in
the 2016-2022 period, searching for an invisibly decaying ${Z^{\prime}}$
boson. In the data analysis, we adopted the same selection criteria identified
in the search for an invisibly decaying Dark Photon, resulting in zero
observed events in the signal box. The signal yield was evaluated by means of
MC simulation: the dominant ${Z^{\prime}}$ production mechanisms where
implemented in the DMG4 software package, integrated in the Geant4-based MC
framework of NA64-$e$, including the momentum-transfer dependency of the
${Z^{\prime}}-e^{-}$ coupling and the effect of the motion of atomic
electrons. Despite the disfavoured ${Z^{\prime}}$ production processes with an
electron beam, the obtained limits prove to be complementary to the NA64-$\mu$
constraints, excluding for the the first time with an accelerator-based
experiment the $g-2$ band in the low (1 keV - 2 MeV) ${Z^{\prime}}$ mass
region.

##  Appendix A  Considerations on the effects of the motion of atomic
electrons

As described in the Sec. 5, an effect to be considered in the simulation of
the $e^{+}e^{-}\rightarrow{Z^{\prime}}$ annihilation process is the
modification of the effective annihilation line shape due to the motion of
atomic electrons in the target. An analytical calculation of this effect,
assuming Compton profiles (CP) and Roothan-Hartree-Fock (RHF) wave functions
for the modelization of the momentum distribution of the atomic shells, has
been performed in Arias-Aragón  _et al._ (2024). The exclusion limits
presented in this work have been obtained using the DMG4 simulation package,
considering an effective exponential momentum distribution, based on the
virial theorem, which allows a simple integration in a Monte Carlo simulation
tool.

Figure 5: The three curves represent the cross sections of the resonant
${Z^{\prime}}$ production on a thin diamond target as a function of the
positron beam nominal energy $E_{B}$, obtained analytically considering the
RHF (orange) and CP (green) model from Arias-Aragón  _et al._ (2024), and from
a simulation via DMG4 (gray). See the text for more details.

In order to check the agreement between the different approaches, we report in
Fig. 5 a comparison between the resonant cross section for a thin diamond
target as a function of the beam energy as reported in Fig.1 (left panel) of
Arias-Aragón  _et al._ (2024) and the same quantity evaluated with a
DMG4-based simulation. The considered mass of the new vector boson is 17 MeV,
while the beam energy spread is 0.5%. We observe that DMG4 (gray curve)
reproduces with very good accuracy the analytical calculation, with a maximum
peak magnitude between the CP (green) and RHF (orange) curves.

###### Acknowledgements.

We gratefully acknowledge the support of the CERN management and staff, and
the technical staffs of the participating institutions for their vital
contributions. This result is part of a project that has received funding from
the European Research Council (ERC) under the European Union’s Horizon 2020
research and innovation programme, Grant agreement No. 947715 (POKER). This
work was supported by the HISKP, University of Bonn (Germany), ETH Zurich and
SNSF Grant No. 186181, No. 186158, No. 197346, No. 216602 (Switzerland), and
FONDECYT (Chile) under Grant No. 1240066, and ANID - Millenium Science
Initiative Program - ICN2019 044 (Chile), and RyC-030551-I and
PID2021-123955NA-100 funded by MCIN/AEI/ 10.13039/501100011033/FEDER, UE
(Spain), and COST Action COSMIC WISPers CA21106, supported by COST (European
Cooperation in Science and Technology. This work is partially supported by
ICSC – Centro Nazionale di Ricerca in High Performance Computing, Big Data and
Quantum Computing, funded by European Union – NextGenerationEU.

## References

  * Bertone  _et al._ (2005) G. Bertone, D. Hooper, and J. Silk, Phys. Rept. 405, 279 (2005), arXiv:hep-ph/0404175 .
  * Arcadi  _et al._ (2018) G. Arcadi, M. Dutra, P. Ghosh, M. Lindner, Y. Mambrini, M. Pierre, S. Profumo, and F. S. Queiroz, Eur. Phys. J. C 78, 203 (2018), arXiv:1703.07364 [hep-ph] .
  * Roszkowski  _et al._ (2018) L. Roszkowski, E. M. Sessolo, and S. Trojanowski, Rept. Prog. Phys. 81, 066201 (2018), arXiv:1707.06277 [hep-ph] .
  * Zurek (2024) K. M. Zurek, Dark Matter Candidates of a Very Low Mass (2024), arXiv:2401.03025 [hep-ph] .
  * Aguillard  _et al._ (2023) D. P. Aguillard _et al._ (Muon g-2), Phys. Rev. Lett. 131, 161802 (2023), arXiv:2308.06230 [hep-ex] .
  * Aoyama  _et al._ (2020) T. Aoyama _et al._ , Phys. Rept. 887, 1 (2020), arXiv:2006.04822 [hep-ph] .
  * Ignatov  _et al._ (2023) F. V. Ignatov _et al._ (CMD-3), (2023), arXiv:2302.08834 [hep-ex] .
  * Borsanyi  _et al._ (2021) S. Borsanyi _et al._ , Nature 593, 51 (2021), arXiv:2002.12347 [hep-lat] .
  * He  _et al._ (1991a) X. G. He, G. C. Joshi, H. Lew, and R. R. Volkas, Phys. Rev. D 43, 22 (1991a).
  * He  _et al._ (1991b) X.-G. He, G. C. Joshi, H. Lew, and R. R. Volkas, Phys. Rev. D 44, 2118 (1991b).
  * Escudero  _et al._ (2019) M. Escudero, D. Hooper, G. Krnjaic, and M. Pierre, JHEP 03, 071, arXiv:1901.02010 [hep-ph] .
  * Kamada and Yu (2015) A. Kamada and H.-B. Yu, Phys. Rev. D 92, 113004 (2015), arXiv:1504.00711 [hep-ph] .
  * Andreev  _et al._ (2024) Y. M. Andreev _et al._ (NA64), Phys. Rev. Lett. 132, 211803 (2024), arXiv:2401.01708 [hep-ex] .
  * Demidov  _et al._ (2019) S. Demidov, S. Gninenko, and D. Gorbunov, Journal of High Energy Physics 2019, 10.1007/jhep07(2019)162 (2019).
  * Andreev  _et al._ (2022) Y. M. Andreev _et al._ (NA64), Phys. Rev. D 106, 032015 (2022), arXiv:2206.03101 [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] .
  * Bauer  _et al._ (2018) M. Bauer, P. Foldenauer, and J. Jaeckel, JHEP 07, 094, arXiv:1803.05466 [hep-ph] .
  * Lees  _et al._ (2016) J. P. Lees _et al._ (BaBar), Phys. Rev. D 94, 011102 (2016), arXiv:1606.03501 [hep-ex] .
  * Czank  _et al._ (2022) T. Czank _et al._ (Belle), Phys. Rev. D 106, 012003 (2022), arXiv:2109.08596 [hep-ex] .
  * Adachi  _et al._ (2024) I. Adachi _et al._ (Belle-II), Search for a $\mu^{+}\mu^{-}$ resonance in four-muon final states at Belle II (2024), arXiv:2403.02841 [hep-ex] .
  * Sirunyan  _et al._ (2019) A. M. Sirunyan _et al._ (CMS), Phys. Lett. B 792, 345 (2019), arXiv:1808.03684 [hep-ex] .
  * Adachi  _et al._ (2023) I. Adachi _et al._ (Belle-II), Phys. Rev. Lett. 130, 231801 (2023), arXiv:2212.03066 [hep-ex] .
  * Ablikim  _et al._ (2024) M. Ablikim _et al._ (BESIII), Phys. Rev. D 109, L031102 (2024), arXiv:2311.01076 [hep-ex] .
  * Kahn  _et al._ (2018) Y. Kahn, G. Krnjaic, N. Tran, and A. Whitbeck, JHEP 09, 153, arXiv:1804.03144 [hep-ph] .
  * Araki  _et al._ (2017) T. Araki, S. Hoshino, T. Ota, J. Sato, and T. Shimomura, Phys. Rev. D 95, 055006 (2017), arXiv:1702.01497 [hep-ph] .
  * Gninenko and Krasnikov (2018) S. N. Gninenko and N. V. Krasnikov, Phys. Lett. B 783, 24 (2018), arXiv:1801.10448 [hep-ph] .
  * Zhang  _et al._ (2021) Y. Zhang, Z. Yu, Q. Yang, M. Song, G. Li, and R. Ding, Phys. Rev. D 103, 015008 (2021), arXiv:2012.10893 [hep-ph] .
  * Fabbrichesi  _et al._ (2020) M. Fabbrichesi, E. Gabrielli, and G. Lanfranchi (2020) SpringerBrief in Physics, arXiv:2005.01515 [hep-ph] .
  * Gninenko  _et al._ (2016a) S. N. Gninenko, N. V. Krasnikov, M. M. Kirsanov, and D. V. Kirpichnikov, Phys. Rev. D 94, 095025 (2016a), arXiv:1604.08432 [hep-ph] .
  * Banerjee  _et al._ (2018a) D. Banerjee _et al._ , Nucl. Instrum. Meth. A 881, 72 (2018a), arXiv:1708.04087 [physics.ins-det] .
  * Depero  _et al._ (2017) E. Depero _et al._ , Nucl. Instrum. Meth. A 866, 196 (2017), arXiv:1703.05993 [physics.ins-det] .
  * Banerjee  _et al._ (2018b) D. Banerjee, V. E. Burtsev, A. G. Chumakov, D. Cooke, P. Crivelli, E. Depero, _et al._ (NA64), Phys. Rev. D 97, 072002 (2018b), arXiv:1710.00971 [hep-ex] .
  * Andreev  _et al._ (2023a) Y. M. Andreev _et al._ (NA64), Phys. Rev. Lett. 131, 161801 (2023a), arXiv:2307.02404 [hep-ex] .
  * Andreas  _et al._ (2013) S. Andreas, S. V. Donskov, P. Crivelli, A. Gardikiotis, S. N. Gninenko, N. A. Golubev, F. F. Guber, A. P. Ivashkin, M. M. Kirsanov, N. V. Krasnikov, V. A. Matveev, Y. V. Mikhailov, Y. V. Musienko, V. A. Polyakov, A. Ringwald, A. Rubbia, V. D. Samoylenko, Y. K. Semertzidis, and K. Zioutas, Proposal for an experiment to search for light dark matter at the sps (2013), arXiv:1312.3309 [hep-ex] .
  * Gninenko  _et al._ (2016b) S. N. Gninenko, N. V. Krasnikov, M. M. Kirsanov, and D. V. Kirpichnikov, Phys. Rev. D 94, 095025 (2016b).
  * Andreev  _et al._ (2023b) Y. M. Andreev _et al._ , Nucl. Instrum. Meth. A 1057, 168776 (2023b), arXiv:2305.19411 [hep-ex] .
  * Banerjee  _et al._ (2020) D. Banerjee _et al._ (NA64), Phys. Rev. Lett. 125, 081801 (2020), arXiv:2005.02710 [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).
  * Ghosh  _et al._ (2000) V. J. Ghosh, M. Alatalo, P. Asoka-Kumar, B. Nielsen, K. G. Lynn, A. C. Kruseman, and P. E. Mijnarends, Phys. Rev. B 61, 10092 (2000).
  * Hunt  _et al._ (2001) A. W. Hunt, D. B. Cassidy, P. A. Sterne, T. E. Cowan, R. H. Howell, K. G. Lynn, and J. A. Golovchenko, Phys. Rev. Lett. 86, 5612 (2001).
  * Arias-Aragón  _et al._ (2024) F. Arias-Aragón, L. Darmé, G. G. di Cortona, and E. Nardi, Production of dark sector particles via resonant positron annihilation on atomic electrons (2024), arXiv:2403.15387 [hep-ph] .
  * Gross (2007) E. Gross, in  _PHYSTAT-LHC Workshop on Statistical Issues for LHC Physics_ (2007).
  * Altmannshofer  _et al._ (2014) W. Altmannshofer, S. Gori, M. Pospelov, and I. Yavin, Phys. Rev. Lett. 113, 091801 (2014), arXiv:1406.2332 [hep-ph] .
  * Gninenko and Gorbunov (2021) S. Gninenko and D. Gorbunov, Phys. Lett. B 823, 136739 (2021), arXiv:2007.16098 [hep-ph] .
  * Gninenko  _et al._ (2018) S. N. Gninenko, D. V. Kirpichnikov, M. M. Kirsanov, and N. V. Krasnikov, Phys. Lett. B 782, 406 (2018), arXiv:1712.05706 [hep-ph] .