††thanks:  0000-0002-6122-4986††thanks:  0000-0002-9584-8877††thanks:
0000-0002-1013-3953††thanks:  0000-0002-2958-456X

# Semi-Visible Dark Photon Phenomenology at the GeV Scale

Asli M. Abdullahi<EMAIL_ADDRESS>Theoretical Physics Department, Fermi
National Accelerator Laboratory, Batavia, IL 60510, USA  Matheus Hostert
<EMAIL_ADDRESS>Perimeter Institute for Theoretical Physics, Waterloo, ON N2J
2W9, Canada  School of Physics and Astronomy, University of Minnesota,
Minneapolis, MN 55455, USA  William I. Fine Theoretical Physics Institute,
School of Physics and Astronomy, University of Minnesota, Minneapolis, MN
55455, USA  Daniele Massaro<EMAIL_ADDRESS>Dipartimento di Fisica
e Astronomia, Università di Bologna, via Irnerio 46, 40126 Bologna, Italy
INFN, Sezione di Bologna, viale Berti Pichat 6/2, 40127 Bologna, Italy  Centre
for Cosmology, Particle Physics and Phenomenology (CP3), Université Catholique
de Louvain, B-1348 Louvain-la-Neuve, Belgium  Silvia Pascoli
<EMAIL_ADDRESS>Dipartimento di Fisica e Astronomia, Università di
Bologna, via Irnerio 46, 40126 Bologna, Italy  INFN, Sezione di Bologna, viale
Berti Pichat 6/2, 40127 Bologna, Italy  CERN, Theoretical Physics Department,
Geneva, Switzerland

###### Abstract

In rich dark sector models, dark photons heavier than tens of MeV can behave
as semi-visible particles: their decays contain both visible and invisible
final states. We present models containing multiple dark fermions which allow
for such decays and inscribe them in the context of inelastic dark matter and
heavy neutral leptons scenarios. Our models represent a generalization of the
traditional inelastic dark matter model by means of a charge conjugation
symmetry. We revisit constraints on dark photons from $e^{+}e^{-}$ colliders
and fixed target experiments, including the effect of analysis vetoes on semi-
visible decays, $A^{\prime}\to\psi_{i}(\psi_{j}\to\psi_{k}\ell^{+}\ell^{-})$.
We find that in some cases the BaBar and NA64 experiments no longer exclude
large kinetic mixing, $\varepsilon\sim 10^{-2}$, and, specifically, the
related explanation of the discrepancy in the muon $(g-2)$. This reopens an
interesting window in parameter space for dark photons with exciting discovery
prospects. We point out that a modified missing-energy search at NA64 can
target short-lived $A^{\prime}$ decays and directly probe the newly-open
parameter space.

††preprint: CERN-TH-2023-019, FERMILAB-PUB-23-054-T, IRMP-CP3-23-08

###### Contents

  1. I Introduction
  2. II Semi-visible dark photons
    1. II.1 Two heavy neutral fermions (HNFs)
    2. II.2 Three HNFs
    3. II.3 Four HNFs
    4. II.4 Mixing with light neutrinos
  3. III Model-independent limits
  4. IV Reinterpreting constraints on invisible $A^{\prime}$
    1. IV.1 BaBar monophoton search
    2. IV.2 NA64 dark photon searches
    3. IV.3 Other limits from beam dumps
  5. V Results
  6. VI Prospects and constraints on HNL and DM interpretations
    1. VI.1 HNFs as dark matter candidates
    2. VI.2 HNFs as heavy neutral leptons
    3. VI.3 Prospects for detection
  7. VII Conclusions
  8. A Simulations and event distributions
  9. B New physics in the anomalous magnetic moment of the muon

##  I Introduction

The existence of hidden sectors containing light and feebly-interacting
particles offers a promising avenue to address the shortcomings of the
Standard Model (SM). With dynamics and structures of their own, these dark
sectors (DS) can contain stable dark matter (DM) particles, new gauge
symmetries, new fundamental scales, and additional sources of $C$ and $P$
violation. While the search for new heavy particles at the LHC continues, this
possibility offers a paradigm-shifting framework that is testable and provides
fertile ground for model building. In particular, if the DS contains light
dark matter particles, their production typically requires the existence of
new mediators that interact with both the DS and SM particles [1, 2, 3, 4].
These can be new scalars (Higgs portal), heavy neutral leptons (neutrino
portal), or new gauge bosons (vector portal). In the case of heavy neutral
leptons, the connection with neutrino masses and leptonic mixing provides
additional theoretical motivation. The presence of new gauge symmetries with
associated Higgs-like breaking mechanisms is also a natural possibility,
appearing in many breaking patterns of grand-unified theories.

Targeting experimental searches for these portal mediators can be an efficient
way to test the DS framework since they couple to both the SM and DS. While
tremendous experimental progress has been achieved for the three portal cases
above (see Refs. [5, 6, 7]), the focus has often been on scenarios with
minimal new particle content. While this is a sensible starting assumption,
relaxing the stringent conditions on minimality can help us to uncover rich DS
theories [8].

In this work, we focus our attention on the dark photon $A^{\prime}$, the
vector portal mediator. Unless explicitly forbidden by new symmetries, kinetic
mixing between the dark $U(1)_{D}$ gauge boson and the SM hypercharge [9],
$\frac{\varepsilon}{2c_{\rm W}}X_{\mu\nu}B^{\mu\nu}$, is expected to be
sizeable, providing a clear target for detection. The one-loop expectation for
kinetic mixing is

$\varepsilon\sim\frac{g^{\prime}g_{D}}{16\pi^{2}}\sum_{i}Q_{i}^{Y}Q_{i}^{X}\log\left(\frac{M_{i}^{2}}{\mu^{2}}\right)\sim\mathcal{O}\left(10^{-3}-10^{-2}\right)$
(1)

where $M_{i}$ and $Q_{i}^{(Y,\,X)}$ are the masses and charges of the heavy
new fermions that run in the loop, $\mu$ the renormalization scale, and
$g_{D}$ the gauge coupling, taken to be of the same order as the SM couplings.
So far, a kinetic mixing of this size has not been experimentally observed for
dark photons in minimal scenarios. Experimental limits have focused primarily
on models in which the $A^{\prime}$ decays to the SM as a fully visible
resonance, or decays invisibly to, e.g., DM particles [10, 11, 12, 13]. At
first sight, the naive expectation for kinetic mixing looks too strongly
constrained to remain a viable possibility, suggesting small dark couplings or
a higher-order origin for $\varepsilon$ [14]. It is, however, possible to
avoid laboratory constraints and remain compatible with the naive one-loop
estimation of $\varepsilon$. One possibility is that the dark photon decays
_semi-visibly_.

Semi-visible decays of the dark photon contain both visible and invisible
particles in the final state, precluding a full reconstruction of the dark
photon mass through its decay products and avoiding constraints from resonance
searches. As we will show, this is a natural prediction of multi-generational
DS models, such as inelastic dark matter models [15]. Nevertheless, the dark
photon may still appear as an invisible particle when produced in an
experiment due to the detector geometry and limited experimental resolution.
The strongest constraints on the GeV-scale invisible dark photon come from the
$e^{+}e^{-}$ colliders and fixed-target experiments [3, 16, 17, 18, 19, 20,
21], and so we revisit the leading constraints from BaBar [12] and NA64 [22,
13, 23, 24].

Another motivation for this work is the measurement of the anomalous magnetic
moment of the muon, $a_{\mu}=(g-2)_{\mu}/2$, at Fermilab (FNAL) [25, 26].
Through the same mechanism understood by Schwinger in the early days of
Quantum Electrodynamics (QED), kinetically-mixed dark photons contribute to
$a_{\mu}$ at the one loop-level with a positive sign, providing an elegant and
simple solution to the discrepancy between experiment and theoretical
predictions, $\Delta a_{\mu}$ [27, 28]. This solution is only possible for
light mediators, $m_{A^{\prime}}\lesssim 3$ GeV, and requires large values of
kinetic mixing, $\varepsilon\sim 10^{-3}-10^{-2}$, compatible with the naive
one-loop expectation. While this explanation is excluded in fully visible or
invisible dark photon models, it remains viable for semi-visible dark photons
in the mass region of $m_{A^{\prime}}\sim 0.6-1$ GeV, as first proposed in
Ref. [29], but later disputed in Refs. [30, 31]. We provide a detailed
analysis of this option, proposing new models that overcome the problems of
the minimal model of Refs. [29, 30, 31]. We also note that dark photons can be
semi-visible due if their decays contain dark showers [32, 33] (see Ref. [34]
for a recent study) and lepton jets [35, 36].

We focus on DS models with multiple fermions, interpreted as either thermal DM
models or seesaw neutrino mass models. In the DM interpretation, our
phenomenological considerations point towards models with dark-photon-mediated
DM coannihilations, automatically satisfying strong Cosmic Microwave
Background (CMB) constraints on light thermal DM. The seesaw interpretation,
albeit less predictive, has several striking predictions for neutrino
experiments [37]. As we will show, the original proposal for a semi-visible
$A^{\prime}$, based on a minimal inelastic DM (iDM) model [29], is strongly
constrained by collider, fixed target, and indirect searches and can only
explain $\Delta a_{\mu}$ in a very narrow region of parameter space. Our
collection of semi-visible $A^{\prime}$ models constitutes a viable
explanation of $\Delta a_{\mu}$ by adding new heavy neutral fermions (HNFs)
with several hundred MeV masses and sizeable mass splittings. Due to the
conservation of a charge-conjugation symmetry, $C$, many of our scenarios
ensure that $A^{\prime}$ couples only off-diagonally to HNF generations,
generalizing the popular iDM scenario.

Testing the allowed parameter space to exclude or discover these models is a
tangible task for current collider and fixed-target experiments. We discuss
some strategies to isolate the distinct semi-visible signatures in these
experiments. In particular, the newly-open parameter space can be explored
with displaced vertices and monophoton-like events at Belle-II, such as those
studied by the BaBar $e^{+}e^{-}$ collider. At fixed-target experiments, we
point out that invisible-$A^{\prime}$ searches can be adapted to be sensitive
to the missing energy in semi-visible dark photon decays. This strategy is
pursued in an accompanying paper to derive new experimental limits using NA64
data [38].

The paper is organized as follows. We introduce the canonical semi-visible
dark photon model in Section II and discuss specific realizations with varying
fermionic content. The set of relevant constraints in the parameter space is
discussed in Section III, and in Section IV, we detail our recasting procedure
to obtain the revised constraints from BaBar and NA64. Our results are then
presented in Section V. We discuss the implications of our results for models
of dark matter and heavy neutral leptons in Section VI, concluding with
Section VII.

##  II Semi-visible dark photons

We are interested in a kinetically-mixed, massive dark photon. In general
terms, the initial Lagrangian is given by

$\displaystyle\mathscr{L}$ $\displaystyle=\mathscr{L}_{\rm
SM}-\frac{\varepsilon}{2c_{\rm
W}}F_{\mu\nu}X^{\mu\nu}-\frac{1}{4}X_{\mu\nu}X^{\mu\nu}$ (2)
$\displaystyle+g_{D}X_{\mu}\mathcal{J}^{\mu}_{D}+\frac{m_{X}^{2}}{2}X_{\mu}X^{\mu},$

where $X_{\mu\nu}$ is the field strength tensor of the dark photon, and
$\mathcal{J}^{\mu}_{D}$ is the dark current containing new fermionic or scalar
degrees of freedom. The origin of the dark photon mass is thus far unspecified
111While introducing a dark Higgs is compelling from a model building point of
view, it also comes with additional assumptions and decreased predictivity.
Therefore, we proceed assuming a Stückelberg mass (see Ref. [39] for a recent
discussion).. After the diagonalization of the gauge kinetic terms, the dark
photon mass eigenstate $A^{\prime}$ with mass $m_{A^{\prime}}\simeq m_{X}$
couples to both the SM electromagnetic (EM) current and the SM weak neutral
current (NC),

$\displaystyle\mathscr{L}_{\rm int.}$ $\displaystyle\supset
A^{\prime}_{\mu}\left(g_{D}\mathcal{J}^{\mu}_{D}-e\varepsilon\mathcal{J}^{\mu}_{\rm
EM}-\varepsilon t_{\rm W}\frac{m_{A^{\prime}}^{2}}{m_{Z}^{2}}\frac{g}{2c_{\rm
W}}\mathcal{J}^{\mu}_{\rm NC}\right)$ (3)
$\displaystyle+Z_{\mu}\left(\frac{g}{2c_{\rm W}}\mathcal{J}^{\mu}_{\rm
NC}+g_{D}t_{\rm W}\varepsilon\mathcal{J}_{\rm
D}^{\mu}\right)+\mathcal{O}(\varepsilon^{2}),$

where $t_{\rm W}\equiv\tan\theta_{\rm W}$, with $\theta_{\rm W}$ the Standard
Model Weak angle. While the SM photon does not couple to the dark sector, the
SM Z boson mass eigenstate can.

Let us now discuss the particle content in the dark sector and how it can
render $A^{\prime}$ semi-visible. To appear as semi-visible, a dark photon has
to decay predominantly into dark particles that cascade-decay into SM states
plus missing energy, avoiding limits on both visible and invisible dark
photons, see Section IV. A particularly simple choice would be a dark complex
scalar $\Phi$, charged under the $U(1)_{D}$. The dark current is
$\mathcal{J}_{D}^{\mu}=\Phi^{*}i\overset{\leftrightarrow}{\partial^{\mu}}\Phi$,
and the hermiticity of the Lagrangian ensures that the dark photon
interactions must always take place between the real components of $\Phi$,
$\varphi_{1}$ and $\varphi_{2}$. A soft $U(1)_{D}$-breaking term,
$\mu\Phi^{2}$ can then split the masses of the real scalars, and render the
dark photon semi-visible due to the decay cascade
$A^{\prime}\to\varphi_{1}(\varphi_{2}\to\varphi_{1}e^{+}e^{-})$. For dark
fermions, the idea is similar, but a much richer structure can arise due to
the half-integer spin. We will consider models with $n$ fermions, where

$\mathcal{J}_{D}^{\mu}\equiv\sum_{i,j=1}^{n}V_{ij}\overline{\psi_{i}}\gamma^{\mu}\psi_{j}$
(4)

with $V_{ij}$ the model-dependent coupling vertices. If each of the dark
fermions has a dark charge $Q_{i}$, the vertices are constrained to
$\sum_{i,j}^{n}|V_{ij}|^{2}\leq\sum_{i}^{n}|Q_{i}|$. While we do not pursue
this possibility, we note that non-renormalizable interactions between
$\psi_{i}$ could be considered, including electric and magnetic moments for
the dark fermions [40, 41, 42, 43]. Like in the scalar case, the latter has
the advantage of being automatically off-diagonal in the $i,j$ fermion indices
(if the dark fermions are Majorana particles), and would also generate semi-
visible decays. For an analogous discussion of semi-visible dark sector
models, see Ref. [44]. In what follows, we focus purely on a fermionic dark
sector with the couplings in Eq. 4.

If the HNFs mix with neutrinos, they are usually called heavy neutral leptons
(HNL) and are typically labeled in the literature as $N_{4}$, $N_{5}$, and so
on. If the lightest is stable, they may be a DM candidate and would typically
be denoted as $\chi_{1}$, with the heavier particles in the spectrum
$\chi_{j=2,\dots,n}$. We retain here a more general notation,
$\psi_{i=1,\dots,n}$, which encompasses both options. The heavier states
$\psi_{j=2,\dots,n}$, are expected to decay in cascades down to the lightest
HNF, emitting two charged particles at each step 222If the dark photon is
heavier than all HNFs, then only three-body decays are allowed. If not
forbidden by the $C$-symmetry, decays into three HNFs are assumed to be
kinematically forbidden.. These models were initially linked to iDM [29, 30],
where the dark photon decays into a heavy and short-lived HNF, and a lighter
HNF, $\psi_{1}$, that is stable and constitutes the DM candidate. In this
work, we generalize this idea to richer dark sectors.

Figure 1:  A summary of the benchmark models we consider, including the HNF
spectra and their interaction vertices with the dark photon, $A^{\prime}$.
Blue lines indicate $C$-even states, and pink lines $C$-odd states. Pseudo-
Dirac states with negligible mass splitting are denoted by two opposite $C$
lines close together. The vertex matrices are defined in the basis
$(\psi_{1},\psi_{2},\dots)$, where the mixing angles are small
$\alpha,\beta\ll 1$, and the entries denoted by $\times$ are arbitrary. Mixed-
color states have no definite $C$ property.

In what follows, we systematically discuss the fermionic content of the models
we study. We start from the case of two Majorana HNFs, as for iDM. Then,
following Ref. [37], we further extend the model to include 3 HNFs, organized
into a pseudo-Dirac pair and one Majorana HNF, or 3 Majorana HNFs, depending
on the choice of parameters. The former option is akin to the iDM model, while
the latter is more general and can accommodate the neutrino mass model of Ref.
[37]. Finally, we present the case of four Majorana HNFs and its limit of two
Dirac particles, recently studied in [45].

In all cases, the dark photon decay into two $\psi_{1}$ needs to be suppressed
as this contributes to its invisible branching ratio, which is severely
constrained by missing-energy searches. We will show how a $C$-symmetry can
satisfy such a requirement.

###  II.1 Two heavy neutral fermions (HNFs)

The simplest model we consider is that of two Majorana HNFs. In the
interaction basis, we take $\chi_{L}$ and $\chi_{R}$ with charges $Q_{L}$ and
$Q_{R}$. The most general Lagrangian for their mass reads

$\displaystyle\mathscr{L}_{\chi}$
$\displaystyle=\overline{\chi_{L}}i\left(\not{\partial}-ig_{D}Q_{L}\not{A}^{\prime}\right)\chi_{L}+\overline{\chi_{R}^{c}}i\left(\not{\partial}+ig_{D}Q_{R}\not{A}^{\prime}\right)\chi_{R}^{c}-\frac{1}{2}\left[\left(\begin{matrix}\overline{\chi_{L}^{c}}&\overline{\chi_{R}}\end{matrix}\right)\left(\begin{matrix}\mu_{L}&m_{D}\\\
m_{D}&\mu_{R}\end{matrix}\right)\left(\begin{matrix}\chi_{L}\\\
\chi_{R}^{c}\end{matrix}\right)+{\rm h.c.}\right]$ (5)

where the Majorana masses $\mu_{L}$ and $\mu_{R}$ break the $U(1)_{D}$ softly
— they can be generated by the vacuum expectation value of a dark Higgs
$\Phi_{2}$ with $Q_{\Phi_{2}}=2Q_{L,R}$. The diagonalization of the symmetric
and complex mass matrix $M$ can be achieved with the usual Takagi
diagonalization, $\text{diag}(m_{1},m_{2})=U^{T}MU$, where
$U=R(\theta)\text{diag}(e^{i\varphi},1)$ with $R(\theta)$ the matrix rotation
by an angle,

$\tan 2\theta=\frac{m_{D}}{\Delta\mu}.$ (6)

We define $\Delta\mu=(\mu_{R}-\mu_{L})/2$ and $\mu=(\mu_{L}+\mu_{R})/2$. $CP$
conservation is ensured when the Majorana phase is $\varphi=0$ or $\pi/2$. In
terms of the Majorana mass eigenstates, the dark current is given by

$\mathcal{J}_{D}^{\mu}=\frac{Q_{A}-Q_{V}\cos{2\theta}}{2}\overline{\psi_{2}}\gamma^{\mu}\gamma^{5}\psi_{2}+\frac{Q_{A}+Q_{V}\cos{2\theta}}{2}\overline{\psi_{1}}\gamma^{\mu}\gamma^{5}\psi_{1}+iQ_{V}\sin{2\theta}\sin\varphi\,\overline{\psi_{2}}\gamma^{\mu}\psi_{1}+Q_{V}\sin{2\theta}\cos\varphi\,\overline{\psi_{2}}\gamma^{\mu}\gamma^{5}\psi_{1},$
(7)

where $Q_{V}\equiv(Q_{L}+Q_{R})/2$ and $Q_{A}\equiv(Q_{L}-Q_{R})/2$. Gauge
anomaly cancellation fixes $Q_{L}=Q_{R}$ ($Q_{A}=0$).

The $C$ symmetry — For $\Delta\mu\to 0$, the mixing angle $\theta$ is maximal,
and the dark photon couples only off-diagonally to the mass eigenstates. The
smallness of the on-diagonal couplings can be understood thanks to a $C$
symmetry. The $C$ operator $U_{c}$ acts on Weyl fermions as

$U_{c}\chi_{L}U_{c}^{-1}=\eta_{c}\psi_{R}^{c},\quad
U_{c}\chi_{R}U_{c}^{-1}=\eta_{c}\psi_{L}^{c},$ (8)

where $\psi_{L}^{c}=C\overline{\psi_{L}}^{T}$ and we choose the phase factor
$\eta_{c}=+1$, for simplicity. In the case of a Dirac fermion, the $C$
operation is achieved by charge conjugating the field,
$\chi=(\zeta,\,-i\sigma^{2}\xi^{*})^{T}\mapsto\chi^{c}=(\xi,\,-i\sigma^{2}\zeta^{*})^{T}$.
The $C$ symmetry is then respected when the Lagrangian is invariant under the
exchange $\xi\leftrightarrow\zeta$ (i.e.,
$\chi_{L}\leftrightarrow\chi_{R}^{c}$).

The left-handed fermions

$\displaystyle\chi_{+}=\frac{\chi_{L}+\chi_{R}^{c}}{\sqrt{2}},\quad\chi_{-}=e^{i\varphi}\left(\frac{\chi_{L}-\chi_{R}^{c}}{\sqrt{2}}\right),$
(9)

constitute the $C$ eigenbasis, where $\varphi$ is the same Majorana phase as
before. Given that $C(A^{\prime}_{\mu})=-1$, the intrinsic $C$-parity of the
fermions can be fixed as $C(\chi_{\pm})=\pm 1$ without loss of generality. The
Lagrangian in Eq. 5 in this basis reads

$\displaystyle\mathscr{L}_{\chi}$
$\displaystyle=\overline{\chi}_{+}i\not{\partial}\chi_{+}+\overline{\chi}_{-}i\not{\partial}\chi_{-}+g_{D}A_{\mu}^{\prime}\left[\frac{Q_{A}}{2}\left(\overline{\chi_{+}}\gamma^{\mu}\gamma^{5}\chi_{+}+\overline{\chi_{-}}\gamma^{\mu}\gamma^{5}\chi_{-}\right)+iQ_{V}\overline{\chi}_{+}\gamma^{\mu}\chi_{-}\right]$
(10)
$\displaystyle-\left[\frac{1}{2}\left(\begin{matrix}\overline{\chi_{-}^{c}}&\overline{\chi_{+}^{c}}\end{matrix}\right)\left(\begin{matrix}m_{D}-\mu&i\Delta\mu\\\
i\Delta\mu&m_{D}+\mu\end{matrix}\right)\left(\begin{matrix}\chi_{-}\\\
\chi_{+}\end{matrix}\right)+{\rm h.c.}\right],$

where we took $\varphi\to\pi/2$ to ensure that the mass terms are positive for
$m_{D}>\mu$. This signals that the two fermions have opposite $CP$ parities.
This basis is identified with the physical basis when $\Delta\mu\to 0$. As
expected, $\Delta\mu$ and $Q_{A}$ are the only parameters that break $C$ in
this model.

In the $C$-symmetric limit, $\chi_{\pm}$ behaves like the components of a
pseudo-Dirac particle with a mass gap $2\mu$. We can also conclude that if the
interactions with the dark photon are off-diagonal in the $C$-conserving
limit, then interactions with a $C$-even dark Higgs in the same limit would be
purely diagonal. In what follows, we assume that if any such scalar degree of
freedom is part of the spectrum, it is heavier than $A^{\prime}_{\mu}$ and has
negligible mixing with the SM Higgs.

Finally, the $C$ symmetry cannot be preserved in the Standard Model due to the
different hypercharges of left- and right-handed fermions. Therefore, we
consider the breaking of $C$ to be stronger in the SM than in the DS. If there
are two $C$ symmetries in the theory, one in the SM, $C_{\rm SM}$, and one in
the DS, $C_{\rm DS}$, then the conservation of $C_{\rm DS}$, but not $C_{\rm
SM}$, would forbid kinetic mixing, $F^{\mu\nu}X_{\mu\nu}\xrightarrow{C_{\rm
DS}}-F^{\mu\nu}X_{\mu\nu}$.

Inelastic dark matter (iDM) — In the $C$ symmetric limit and with an anomaly-
free charge assignment, $Q_{A}=0$, we recover the well-known iDM model [15,
46]. Taking $\varphi=\pi/2$ and $Q_{V}=1$, the dark current is simply

$\displaystyle\mathcal{J}_{\rm iDM}^{\mu}$
$\displaystyle=i\overline{\psi_{2}}\gamma^{\mu}\psi_{1}+\text{h.c.}$ (11)

This phenomenological model is defined by the following five parameters

$m_{1},\,\,\Delta_{21}\equiv\frac{m_{2}-m_{1}}{m_{1}},\,\,r\equiv\frac{m_{1}}{m_{A^{\prime}}},\,\,\alpha_{D}\equiv\frac{g_{D}^{2}}{4\pi},\text{
and }\varepsilon.$ (12)

Historically, the interest in these models stemmed from explanations of the
DAMA observation and the fact that the energy threshold for direct DM
detection, as induced by $\Delta_{21}$, varies between Sodium-Iodine and Xenon
experiments [15]. This explanation has since been ruled out by other direct
detection experiments [47, 48, 49]. Still, the interest in iDM has persisted,
especially in the context of accelerator experiments [50, 51, 52, 53], where
the co-annihilator can be searched for through its displaced decays.

Self-interactions of $\psi_{1}$ can only proceed through scalars in this
model. If the scalar mixes with the Higgs, $\lambda|\Phi|^{2}|H|^{2}$, direct
annihilation to SM fermions can take place. The mixing with the Higgs induces
couplings smaller than the Higgs’ Yukawa couplings with SM fermions;
annihilations are suppressed by the fermion masses, in addition to the Higgs
mixing parameter. The mixing should, therefore, be small enough for the scalar
contribution to the self-annihilation of DM to be sub-dominant to the
exponentially suppressed $A^{\prime}$-mediated coannihilations. In addition,
we assume the dark higgs is heavier than the HNFs, as otherwise secluded
annihilation would dominate 333This scenario was explored in Ref. [31]. The
authors also find that BaBar does not rule out the entire $\Delta a_{\mu}$
region of preference but did not study constraints from NA64 and Higgs
decays..

As we will show in Section V, the explanation of $\Delta a_{\mu}$ in this
model is in tension with invisible dark photon limits. Each dark photon
produced can only be accompanied by a single semi-visible decay,
$\psi_{2}\to\psi_{1}f^{+}f^{-}$, with $f$ a SM particle. Even small losses of
acceptance in the detector can lead to a missed $e^{+}e^{-}$ pair. With this
limitation in mind, we consider new models where multiple unstable fermions
accompany dark photon production.

There are two ways to achieve this:  1. in models where $\psi_{2}$ can be
produced in pairs, $A^{\prime}\to\psi_{2}\psi_{2}$, and  2. in models of three
or more HNFs, where the dark photon couples predominantly to the heaviest and
most short-lived states, e.g. $A^{\prime}\to\psi_{3}\psi_{2}$.

We explore these possibilities in models with three and four Weyl fermions,
always reducing the phenomenological model to, at most, three distinguishable
states in the spectrum. This allows for better compatibility between the
$\Delta a_{\mu}$ anomaly and a dark photon explanation. We will prove this in
Section V with a detailed analysis.

###  II.2 Three HNFs

We start by extending the two HNF model by a single fully sterile Weyl
fermion. We keep the two dark fermions, $\chi_{R}$ and $\chi_{L}$, with the
same charges $Q_{L}=Q_{R}=1$, and introduce a new singlet fermion, $\eta_{L}$.
This content is a simplified version of the three-portal model of Ref. [37].

The Lagrangian is given by

$\displaystyle\mathscr{L}_{\rm 3-HNF}$
$\displaystyle=\mathscr{L}_{\chi}+\overline{\eta_{L}}i\not{\partial}\eta_{L}\
$ (13)
$\displaystyle-\left[\frac{\mu_{L}^{\prime}}{2}\overline{\eta_{L}^{c}}\eta_{L}+\Lambda_{L}\overline{\eta_{L}^{c}}\chi_{L}+\Lambda_{R}\overline{\eta_{L}^{c}}\chi_{R}^{c}+\text{
h.c.}\right].$

The mixing terms break the $U(1)_{D}$ and can be generated by the vacuum
expectation value of a scalar particle $\Phi_{1}$, which carries charge
$Q_{\Phi_{1}}=1$. In that case, $\Lambda_{L,R}\equiv
Y_{L,R}\,v_{\Phi_{1}}/\sqrt{2}$ where $Y_{L,R}$ are the Yukawa couplings.
Another dark Higgs $\Phi_{2}$ with $Q_{\Phi_{2}}=2$ could generate the
Majorana masses of the two dark fermions after symmetry breaking.

Because $\eta_{L}$ is completely neutral, it can couple to the SM lepton
doublets via the Yukawa coupling $\overline{L}\tilde{H}\eta_{L}^{c}$. While
these terms play an essential role in the mass generation of light neutrinos,
they only give a small correction to the HNF masses. The neutrino Yukawa
coupling is constrained to be small and will have no impact on the collider
and fixed-target phenomenology we discuss. We will consider the impact of this
coupling on neutrino mass generation in Section II.4. These terms allow the
lightest HNF to decay into SM neutrinos. To ensure the stability of the dark
matter candidate, we forbid these terms by charging all DS fermions,
$\chi_{L},\chi_{R},\eta_{L}$, under a dark parity, e.g. $Z_{2}$ symmetry. This
dark parity can also be attributed to the conservation of lepton number if
$L(\eta_{L})=0$, which would forbid the neutrino Yukawa coupling [54].

We use the left-handed dark fermion basis $\chi_{+}$ and $\chi_{-}$ introduced
in Eq. 9, and set the Majorana phases to be such that $CP$ is conserved and
the mass terms are positive when $M_{X}>\mu$. In that case, the DS fermion
mass matrix is

$\displaystyle-\mathscr{L_{\rm 3-HNF}}$
$\displaystyle\supset\frac{1}{2}\left(\begin{matrix}\overline{\eta_{L}^{c}}&\overline{\chi_{-}^{c}}&\overline{\chi_{+}^{c}}\end{matrix}\right)\left(\begin{matrix}\mu_{L}^{\prime}&\Delta\Lambda&\Lambda\\\
\Delta\Lambda&M_{X}-\mu&\Delta\mu\\\
\Lambda&\Delta\mu&M_{X}+\mu\end{matrix}\right)\left(\begin{matrix}\eta_{L}\\\
\chi_{-}\\\ \chi_{+}\end{matrix}\right)+{\rm h.c.},$ (14)

where $\Lambda=(\Lambda_{L}+\Lambda_{R})/\sqrt{2}$ and
$\Delta\Lambda=(\Lambda_{R}-\Lambda_{L})/\sqrt{2}$. Imposing the $C$ symmetry
in the $\chi$ sector 444Just like in the SM, the $C$ symmetry is broken in
this model due to the odd number of Weyl fermions. This would indicate that
only the $U(1)_{D}$-charged sector respects $C$., we recover the limit where
$\Delta\mu=\Delta\Lambda=0$. We find an analogous situation to the two HNF
cases, with the difference that $\chi_{+}$ can now mix with a sterile state.
Indeed, since C$(\eta_{L})=+1$, $C$-conservation implies that only the
$C$-even fermion can mix with $\eta_{L}$. As we will see, the spectrum can
consist of one Dirac and one Majorana particle or three Majorana states.

The $C$-odd state $\chi_{-}\equiv\psi_{2}$ decouples, and $\eta$ and
$\chi_{+}$ mix. A single rotation in the $C$-even sector leads to the mass
basis,

$\displaystyle\psi_{1}$ $\displaystyle=c_{\alpha}\eta+s_{\alpha}\chi_{+},$
$\displaystyle m_{1}=\mu_{L}^{\prime}-M\frac{\sin^{2}\alpha}{\cos{2\alpha}}$
(15) $\displaystyle\psi_{2}$ $\displaystyle=\chi_{-},$ $\displaystyle
m_{2}=M_{X}-\mu_{L}^{\prime}$ (16) $\displaystyle\psi_{3}$
$\displaystyle=-s_{\alpha}\eta+c_{\alpha}\chi_{+},$ $\displaystyle
m_{3}=\mu_{L}^{\prime}+M\frac{\cos^{2}\alpha}{\cos{2\alpha}}$ (17)

with $\tan 2\alpha=2\Lambda/M$ and $M=M_{X}+\mu-\mu_{L}^{\prime}$. If
$\tan{2\alpha}\ll 1$, a seesaw mechanism is in place, and $\psi_{2}$ and
$\psi_{3}$ form a pseudo-Dirac pair. The other possibility, $\tan 2\alpha\gg
1$, does not preserve the pseudo-Dirac limit. The splittings in the model are
then given by $m_{3}-m_{1}\sim M+2\Lambda^{2}/M$ and
$m_{3}-m_{2}\sim\Lambda^{2}/M+2\mu$.

In the $C$ symmetric case, the dark current is also fully off-diagonal and is
given by

$\mathcal{J}_{\rm 3-HNF}^{\mu}\supset
s_{\alpha}\overline{\psi_{2}}\gamma^{\mu}\psi_{1}+c_{\alpha}\overline{\psi_{2}}\gamma^{\mu}\psi_{3}+\text{
h.c.}$ (18)

In the limit of small $\alpha$, the dark photon interacts more strongly with
the pseudo-Dirac pair. When the $\psi_{1}$ HNF is a dark matter particle, its
relic abundance is set exclusively through coannihilation with the heavier
pseudo-Dirac partner. Similarly to the minimal iDM model, this mechanism
evades constraints from the CMB and direct detection experiments. Below, we
highlight the two types of phenomenological models that can be derived from
Eq. 13 above.

Mixed inelastic dark matter (mixed-iDM) — The first phenomenological scenario
we can consider is the limit where two of the Weyl fermions make up a mostly-
dark, pseudo-Dirac particle, while the third Weyl fermion remains a Majorana
particle, mostly in the direction of the sterile state. In this case, the
lighter, mostly-sterile Majorana fermion would constitute dark matter, while
the mostly-dark pseudo-Dirac fermion plays the role of the co-annihilator. In
this case, the self-annihilation of dark matter via $A^{\prime}$ interactions
is forbidden by the $C$ symmetry, and not constrained by CMB limits. The model
is a trivial extension of the iDM model and invokes the same $C$ symmetry used
there.

Considering the Majorana state $\psi_{1}$ and a (pseudo-)Dirac state
$\Psi_{2}$, the dark current is

$\mathcal{J}_{\rm mixed-iDM}^{\mu}\supset
s_{\alpha}\overline{\Psi_{2}}\gamma^{\mu}\psi_{1}+c_{\alpha}\overline{\Psi_{2}}\gamma^{\mu}\Psi_{2}+\text{h.c.},$
(19)

and so this model is fully specified by Eq. 12 and $\alpha$. To make use of
the model above, one must guarantee the coherence of the Majorana states
$\psi_{2}$ and $\psi_{3}$ in the pseudo-Dirac state $\Psi_{2}$, so that we are
justified in treating them as a single Dirac particle in the phenomenological
work. Note that $\Delta_{32}$ will only play a minor role in the dark matter
hypothesis since the relevant splitting for coannihilations is the one between
$\psi_{1}$, the dark matter candidate, and $\Psi_{2}$, its interaction
partner. We can express $M_{X}$ and $\Lambda$ in terms of $\tan 2\alpha$,
which controls the decay rate $\Psi_{2}\rightarrow\psi_{1}$, and the splitting
$\Delta_{21}\equiv(m_{2}-m_{1})/m_{1}$, which has an important impact on the
coannihilation rate for dark matter. For $\mu<M_{X}$, we find

$\displaystyle m_{1}$
$\displaystyle\simeq\mu_{L}^{\prime}-\frac{1}{4}\mu\Delta_{21}\tan^{2}2\alpha~{},$
(20) $\displaystyle m_{2}$ $\displaystyle=\mu_{L}^{\prime}(1+\Delta_{21})~{},$
(21) $\displaystyle m_{3}$
$\displaystyle\simeq\mu_{L}^{\prime}(1+\Delta_{21})+\frac{1}{4}\mu\Delta_{21}\tan^{2}2\alpha~{}.$
(22)

We notice that
$\Delta_{32}=\frac{1}{4}\frac{\Delta_{21}}{1+\Delta_{21}}\tan^{2}2\alpha$ and
is small as far as the condition $\tan^{2}2\alpha\ll 1$ holds. For
$\mu_{L}^{\prime}=\mu=0$, the splitting of the pseudo-Dirac pair is
$\Delta_{32}\propto\alpha^{2}$, so that in the limit of $\alpha\to 0$, we
recover an exact Dirac state. In summary, provided the mixing angle is small,
the $\Delta_{32}$ splitting is negligible. The decay rate for three-body
decays like $\psi_{3}\to\psi_{2}+\dots$ are suppressed by
$\Delta_{32}^{5}\propto\alpha^{10}$, and so, can be safely neglected for the
mixing angles considered here.

Three Majorana fermions — Relaxing the condition on $\alpha\ll 1$ and the $C$
symmetry in the dark sector, it is possible to split the mass eigenstates away
from a heavy Pseudo-Dirac pair and therefore have three hierarchical Majorana
HNFs. This structure enhances the semi-visible decay rates of $\psi_{3}$ and
$\psi_{2}$ while suppressing the on-diagonal terms in the dark sector current.
The benchmark points exemplify this in Ref. [37]. This case is of interest for
providing both a viable inelastic DM model, compatible with CMB bounds, or,
alternatively, a heavy neutral lepton interpretation with interesting
phenomenological consequences, e.g. Ref. [37].

We can also obtain some useful approximate formulas. For the $CP$ conserving
case, a mild hierarchy can be obtained for large values of $\tan 2\alpha$. For
$0<\mu\simeq M_{X}$, we have $(M_{X}-\mu)/M_{X}\tan^{2}2\alpha\simeq
2\Delta_{21}/(1+\Delta_{21})$ and
$\Delta_{32}\simeq\Delta_{21}/(1+\Delta_{21})$, implying a mildly hierarchical
HFN spectrum. Moving away from the $C$ symmetric limit, a sizeable
$\Delta\Lambda$ can lead to a stronger hierarchy of masses, as, for instance,
in benchmark BP5.

Depending on the mass hierarchy, it would also be possible for the heavy
states to decay into multiple lighter HNFs, in particular into 3 $\psi_{1}$.
In the case under consideration, such decays are kinematically forbidden. We
avoid this possibility as these decay channels could easily dominate and
enhance the invisible branching ratio of the dark photon if $\psi_{1}$ is
stable or long-lived, as is typical in these models.

###  II.3 Four HNFs

If we further enlarge the fermionic sector, it is possible to recover a
2-Dirac fermion picture. Two families of HNF exist: one neutral under all
gauge symmetries, $\eta$, and one charged under the dark gauge symmetry,
$\chi$. Our Lagrangian in this case reads

$\mathscr{L}=\mathscr{L}_{\chi}+\overline{\eta}i\not{\partial}\eta-
M_{\eta}\overline{\eta}\eta-\left[\frac{\mu^{\prime}_{R}}{2}\overline{\eta_{R}}\eta_{R}^{c}+\Lambda_{L}^{\prime}\overline{\eta_{R}}\chi_{L}+\Lambda_{R}^{\prime}\overline{\eta_{R}}\chi_{R}^{c}+\frac{\mu_{L}^{\prime}}{2}\overline{\eta_{L}^{c}}\eta_{L}+\Lambda_{L}\overline{\eta_{L}^{c}}\chi_{L}+\Lambda_{R}\overline{\eta_{L}^{c}}\chi_{R}^{c}+\text{h.c.}\right],$
(23)

where again we have omitted potential Yukawa couplings between SM neutrinos
and the sterile fermions, $\eta_{L}$ and $\eta_{R}$ (cf. Section II.2). In the
$C$ symmetric limit, one can show that $\Lambda_{L}=\Lambda_{R}^{\prime}$ and
$\Lambda_{L}=\Lambda_{R}^{\prime}$. In this limit, for an appropriate choice
of Majorana phases, the mass matrix in the $C$ eigenbasis is,

$\displaystyle-\mathscr{L_{\rm 4-HNF}}$
$\displaystyle\supset\frac{1}{2}\left(\begin{matrix}\overline{\eta_{-}^{c}}&\overline{\eta_{+}^{c}}&\overline{\chi_{-}^{c}}\overline{\chi_{+}^{c}}\end{matrix}\right)\left(\begin{matrix}M_{\eta}-\mu^{\prime}&0&\Lambda_{-}&0\\\
0&M_{\eta}+\mu^{\prime}&0&\Lambda_{+}\\\ \Lambda_{-}&0&M_{X}-\mu&0\\\
0&\Lambda_{+}&0&M_{X}+\mu\end{matrix}\right)\left(\begin{matrix}\eta_{-}\\\
\eta_{+}\\\ \chi_{-}\\\ \chi_{+}\end{matrix}\right)+{\rm h.c.},$ (24)

where
$\Lambda_{\pm}\equiv(\Lambda_{L}^{\prime}+\Lambda_{R})/2\pm(\Lambda_{L}+\Lambda_{R}^{\prime})/2$.
The $C$-even and $C$-odd sectors decouple. We introduce two commuting
rotations defined by the mixing angles $\tan
2\beta_{\pm}=2\Lambda_{\pm}/\Delta_{\pm}$, where
$\Delta_{\pm}=\pm(M_{X}-M_{\eta})+\mu-\mu^{\prime}$. The spectrum is then
given by

$\displaystyle\psi_{1}$
$\displaystyle=c_{\beta_{-}}\eta_{-}+s_{\beta_{-}}\chi_{-},$ $\displaystyle
m_{1}=M_{\eta}-\mu^{\prime}+\Delta_{-}\frac{\sin^{2}\beta_{-}}{\cos{2\beta_{-}}}$
$\displaystyle\psi_{2}$
$\displaystyle=c_{\beta_{+}}\eta_{+}+s_{\beta_{+}}\chi_{+},$ $\displaystyle
m_{2}=M_{\eta}+\mu^{\prime}-\Delta_{+}\frac{\sin^{2}\beta_{+}}{\cos
2\beta_{+}}$ $\displaystyle\psi_{3}$
$\displaystyle=s_{\beta_{-}}\chi_{-}-c_{\beta_{-}}\eta_{-},$ $\displaystyle
m_{3}=M_{X}-\mu-\Delta_{-}\frac{\sin^{2}\beta_{-}}{\cos{2\beta_{-}}}$
$\displaystyle\psi_{4}$
$\displaystyle=s_{\beta_{+}}\eta_{+}-c_{\beta_{+}}\chi_{+},$ $\displaystyle
m_{4}=M_{X}+\mu+\Delta_{+}\frac{\sin^{2}\beta_{+}}{\cos 2\beta_{+}}$

When $\mu,\mu^{\prime},\Lambda\ll M_{X},M_{\eta}$, the spectrum is composed of
two pseudo-Dirac particles, split by the $U(1)_{D}$-breaking terms.

BP | model | $r$ | $\Delta_{21}$ | $\Delta_{32}$ | $\alpha_{D}$ | $V_{11}$ | $V_{21}$ | $V_{22}$ | $V_{31}$ | $V_{32}$ | $V_{33}$ | Comment
---|---|---|---|---|---|---|---|---|---|---|---|---
$/10^{-2}$
1a | iDM | $1/3$ | $0.5$ | $-$ | 0.5 | $0$ | $1$ | $0$ | $0$ | $0$ | $0$ | -
1b | iDM | $1/3$ | $0.4$ | $-$ | 0.1 | $0$ | $1$ | $0$ | $0$ | $0$ | $0$ | same as [29]
2a | mixed-iDM | $1/3$ | $0.3$ | $-$ | 0.5 | $0$ | $s_{\alpha}c_{\alpha}$ | $c_{\alpha}^{2}$ | $-$ | $-$ | $-$ | $\alpha=8^{\circ}$
2b | mixed-iDM | $1/3$ | $0.3$ | $-$ | 0.5 | $0$ | $s_{\alpha}c_{\alpha}$ | $c_{\alpha}^{2}$ | $-$ | $-$ | $-$ | $\alpha=4^{\circ}$
3a | i2DM | $1/3$ | $0.4$ | $-$ | 0.5 | $s_{\beta}^{2}$ | $s_{\beta}c_{\beta}$ | $c_{\beta}^{2}$ | $-$ | $-$ | $-$ | $\beta=8.6^{\circ}$
3b | i2DM | $1/3$ | $0.4$ | $-$ | 0.5 | $s_{\beta}^{2}$ | $s_{\beta}c_{\beta}$ | $c_{\beta}^{2}$ | $-$ | $-$ | $-$ | $\beta=4.6^{\circ}$
3c | i2DM | $1/3$ | $0.4$ | $-$ | 0.5 | $s_{\beta}^{2}$ | $s_{\beta}c_{\beta}$ | $c_{\beta}^{2}$ | $-$ | $-$ | $-$ | $\beta=2.3^{\circ}$
3d | i2DM | $1/3$ | $0.4$ | $-$ | 0.5 | $s_{\beta}^{2}$ | $s_{\beta}c_{\beta}$ | $c_{\beta}^{2}$ | $-$ | $-$ | $-$ | $\beta=1.1^{\circ}$
4a | 3 HNFs | $0.11$ | $2.44$ | $0.54$ | $0.3$ | $0$ | $3.9$ | $0$ | $0$ | $99$ | $0$ | same as [37]
4b | 3 HNFs | $0.16$ | $2.44$ | $0.54$ | $0.3$ | $0$ | $3.9$ | $0$ | $0$ | $99$ | $0$ | same as [37]
4c | 3 HNFs | $0.15$ | $0.85$ | $0.77$ | $0.3$ | $0$ | $0.10$ | $0$ | $0$ | $99$ | $0$ | same as [37]
5 | 3 HNFs | $0.16$ | $0.573$ | $0.586$ | $0.3$ | $0.40$ | $7.8$ | $8.3$ | $2.8$ | $98$ | $69$ | same as [37]
Table 1: The benchmark models for semi-visible dark photons used in this work.
In the second column, we specify the type of model considered. Here,
$r=m_{1}/m_{A^{\prime}}$ and $\Delta_{ij}=(m_{i}-m_{j})/m_{j}$. The dark
photon coupling vertices $V_{ij}$ are defined in Eq. 4.

One more limit of interest is considering the $U(1)_{D}$ to be exclusively
broken by one unit, such that $\mu=\mu^{\prime}=0$ and
$\Delta_{+}=-\Delta_{-}\equiv\Delta$. In that case, the Dirac pairs are split
by

$\Delta_{43}\sim\Delta_{12}\sim\Delta(\beta_{+}^{2}-\beta_{-}^{2}),$ (25)

which is small for small mixing angles and vanishes when
$\beta_{+}=\beta_{-}$. This last regime is the limit where the two pairs
compose exact Dirac fermions, achieved in two cases: i)
$\Lambda_{+}=-\Lambda_{-}$ ($\Lambda_{L}^{\prime}=\Lambda_{R}=0$), or ii)
$\Lambda_{+}=\Lambda_{-}$ ($\Lambda_{L}=\Lambda_{R}^{\prime}=0$).

In terms of the mass eigenstates, the dark current takes the simple form,

$\displaystyle\mathcal{J}_{X}^{\mu}$
$\displaystyle=c_{\beta_{+}}c_{\beta_{-}}\overline{\psi_{4}}\gamma^{\mu}\psi_{3}+s_{\beta_{+}}c_{\beta_{-}}\overline{\psi_{4}}\gamma^{\mu}\psi_{1}$
(26)
$\displaystyle+s_{\beta_{-}}c_{\beta_{+}}\overline{\psi_{3}}\gamma^{\mu}\psi_{2}+s_{\beta_{+}}s_{\beta_{-}}\overline{\psi_{2}}\gamma^{\mu}\psi_{1}+\text{h.c.},$

where only interactions between $C$-odd and $C$-even states are allowed, and
where the heaviest pseudo-Dirac pair couples most strongly to the dark photon.
Decays of the type $\psi_{4,2}\to\psi_{3,1}+\dots$ are suppressed with respect
to the dominant $\psi_{4,3}\to\psi_{1,2}+\dots$ by factors of
$(\beta_{+}^{2}-\beta_{-}^{2})^{5}\beta_{\pm}^{2}$. In fact, for the typical
mixing angles we consider, the particle $\psi_{2}$ is semi-stable, as
$m_{2}-m_{1}<2m_{e}$.

In summary, in the $C$ symmetric limit, we find a pair of pseudo-Dirac
particles, each split by a small gap, proportional to
$\Delta(\beta_{+}^{2}-\beta_{-}^{2})$, where $\beta_{+}$ and $\beta_{-}$ are
the mixing angles in the $C$-even and $C$-odd sectors, respectively. The
individual splittings are only relevant for large mixings, and vanish exactly
in the limit $\beta_{+}=\beta_{-}$.

Inelastic Dirac dark matter (i2DM) — In the exact Dirac limit, a simple
phenomenological model arises with only two particles in the spectrum. A
light, mostly-neutral Dirac fermion $\Psi_{1}$, constituting a dark matter
candidate, and a heavier, mostly-dark Dirac fermion $\Psi_{2}$. In terms of
these degrees of freedom, the dark current is given by

$\mathcal{J}_{\rm
i2DM}^{\mu}=s_{\beta}^{2}\overline{\Psi_{1}}\gamma^{\mu}\Psi_{1}+s_{\beta}c_{\beta}\left(\overline{\Psi_{2}}\gamma^{\mu}\Psi_{1}+\text{h.c.}\right)+c_{\beta}^{2}\overline{\Psi_{2}}\Psi_{2}.$
(27)

As expected, $\Psi_{2}$ is coupled more strongly to the dark photon, and,
therefore, will aid in the depletion of $\Psi_{1}$ particles in the freeze-out
mechanism through its coannihilations. The mixing-suppressed self-interactions
of the dark matter particle weaken the CMB limits on $\Psi_{1}\Psi_{1}\to
e^{+}e^{-}$. The relic density is typically set by the self-annihilation of
the heavy partner, $\Psi_{2}\Psi_{2}\to f^{+}f^{-}$, and coscattering,
$\Psi_{2}\Psi_{2}\leftrightarrow\Psi_{2,1}\Psi_{1}$ and
$\Psi_{2}f\leftrightarrow\Psi_{1}f$ [45]. This model differs from the mixed-
iDM scenario mostly in that the off-diagonal interactions between dark matter
and its co-annihilator are suppressed with respect to the self-interactions of
the co-annihilator, $\Psi_{2}$. The phenomenology is fully determined by the
parameters in Eq. 12, in addition to $\beta$. Similar ideas of a sterile dark
matter particle co-annihilating with heavier dark partners have been explored
before in the context of a toy model with two Majorana fermions [55].

With regards to the accelerator phenomenology, the branching ratios of the
dark photon to the lighter fermions will be hierarchical, approximately
following a proportion of $(1:\beta^{2}:\beta^{4})$ for decays into
$(\Psi_{2}\Psi_{2},\Psi_{2}\Psi_{1},\Psi_{1}\Psi_{1})$ final states. The
dominance of $A^{\prime}\to\Psi_{2}\Psi_{2}$ decays guarantees a large number
of events with two semi-visible particles, further relaxing constraints on
kinetic mixing coming from invisible dark photon searches. The presence of
more visible final states enhances the prospects for discovery.

###  II.4 Mixing with light neutrinos

So far we have considered a secluded sector that feebly interacts with the SM
only via the vector (and possibly scalar) portal. Generically, in the presence
of sterile fermions, Yukawa couplings with both the SM leptonic doublet and
the DS are allowed, and, after symmetry breaking, will lead to mixing between
neutrinos and HNFs. Conventionally, the HNFs are called HNLs in this scenario.

The HNLs are unstable, as they can always decay to e.g. neutrinos, and the
lightest particle in the spectrum cannot constitute DM. In order to recover a
stable candidate for DM, it is necessary to advocate an additional symmetry
which distinguishes the HNLs from the light neutrinos and forbids Yukawa
couplings with the leptonic doublets. The simplest such symmetry is a $Z_{2}$
and would guarantee the stability of the lightest HNL.

The HNL scenario is most easily realized in models with three or more HNFs, in
which the neutral fermions $\eta$ are free to mix with the SM neutrinos, in
the absence of any additional symmetries. For the model in Section II.2, one
can add the following Yukawa interaction

$\mathscr{L}=\mathscr{L}_{3{\rm-
HNF}}-\sum_{\alpha=e,\mu,\tau}\left(y_{\alpha}\overline{L_{\alpha}}\widetilde{H}\eta_{L}^{c}+\text{h.c.}\right),$
(28)

where $L_{\alpha}$ is the $SU(2)$ leptonic doublet of the SM, and
$\widetilde{H}=i\sigma_{2}H^{*}$ the conjugate of the Higgs doublet. Similar
terms involving $\eta_{R}$ could be added to the four HNF models of Section
II.3. After EW and $U(1)_{D}$ symmetry breaking, the HNFs mix amongst
themselves and we can justifiably call them HNLs, $N_{i}\equiv\psi_{i}$. With
the addition of a dark scalar $\Phi$ with a dark charge $Q_{\Phi}=1$, the
Lagrangian above is identical to the model discussed in Ref. [37].

The mixing of active SM neutrinos and HNLs is constrained to be small by
direct laboratory searches. For the values of kinetic mixing and $A^{\prime}$
mass considered in this paper, the lightest HNL, $N_{4}$, will decay via
$N_{4}\to\nu\ell^{+}\ell^{-}$ or $N_{4}\to\nu\pi^{+}\pi^{-}$. Due to the
suppression of the small mixing with the neutrinos, $N_{4}$ is usually long-
lived ($c\tau_{N_{4}}\gg 1$ m) and constitutes missing energy at $e^{+}e^{-}$
colliders and fixed-target experiments. The experimental consequences of the
mixing with neutrinos is discussed in Section VI.2. For a review of the
phenomenology of HNLs, see Ref. [56].

A key consequence of this setup is light neutrino mass generation. In fact, a
GeV-scale seesaw mechanism as the origin of the observed neutrino masses and
mixing angles has been extensively studied in the literature [57, 58, 59, 60,
61, 62, 63, 64, 65]. For a review of this topic, see Ref. [5] and references
therein.

It is also interesting to consider the role of lepton number in these
scenarios. Charging $\eta_{L}$ so that the Yukawa coupling is lepton number
conserving implies that the Majorana mass term $\mu^{\prime}_{L}$ breaks it by
two units. In the dark sector, the charge assignment of $\chi_{L}$ and
$\chi_{R}$ is arbitrary and, depending on the specific choice, lepton number
will be broken by $\Lambda_{L,R}$, $M_{X}$ and $\mu_{L,R}$, or by
$\Lambda_{L}$ and $\mu_{L,R}$, for $L(\chi_{L})=L(\chi_{R}^{c})=0$,
$L(\chi_{L})=L(\chi_{R}^{c})=1$ or $L(\chi_{L})=-L(\chi_{R}^{c})=1$,
respectively. Both $\Lambda_{L,R}$ and $\mu_{L,R}$ terms also break $U(1)_{D}$
by one and two units, respectively, and can arise once multiple scalars,
carrying $U(1)_{D}$ charges, develop a vacuum expectation value. In the most
minimal case of one scalar $\Phi$ with dark charge $Q_{\Phi}=1$, either the
resulting $\Lambda_{L}$, or $\Lambda_{R}$, term breaks the lepton number
explicitly by 2 units. We leave further theoretical considerations to future
work. We also note that the C-symmetry introduced earlier is not compatible
with $U(1)_{L}$ in this minimal realization if $\eta_{L}$ is charged under
$U(1)_{L}$.

Light neutrino masses need to depend on all the $U(1)_{L}$-breaking
parameters. As an interesting example, let us consider the case in which the
$\chi_{L,R}$ do not carry lepton number and only one scalar is included in the
theory, so that $\mu_{L}=\mu_{R}=0$. We assume that lepton number violating
terms are small, implying $\Lambda_{L,R}\ll M_{X}$ after $U(1)_{D}$ breaking.
Another choice for the charges consists in having all new fields neutral. In
this case, the lepton number violating term is the Yukawa coupling itself,
explaining naturally its smallness. For negligible $\mu^{\prime}_{L}$ and
$\Lambda_{L,R}\ll M_{X}$, we have $m_{1}\simeq\Lambda^{2}/M$, $m_{2}\simeq
m_{3}=M$ and light neutrino masses arise

$m_{\nu}\simeq\frac{y^{2}v_{H}^{2}}{\Lambda^{2}}M.$ (29)

Note that this is a one-generation estimate, but a full flavor analysis is
needed to determine the values of the three light masses and mixing
parameters.

It should be pointed out that additional contributions can also come from
loops, especially if $\mu^{\prime}_{L}$ is large, and they can be significant
owing to the fact that the scale of symmetry breaking in the dark sector is
smaller than the electroweak scale [66].

The case with four HNFs is even richer in possibilities owing to the enlarged
fermionic sector. In this case, it is possible to add to Yukawa interactions
with the SM

$\mathscr{L}=\mathscr{L}_{4{\rm-
HNF}}-\sum_{\alpha=e,\mu,\tau}\left(y_{\alpha}\overline{L_{\alpha}}\widetilde{H}\eta_{R}+y_{\alpha}^{\prime}\overline{L_{\alpha}}\widetilde{H}\eta_{L}^{c}+\text{h.c.}\right).$
(30)

One option is not to charge either of $\eta_{L,R}$ implying that both Yukawa
coupings are suppressed being lepton number violating. On the contrary, if
$L(\eta_{L})=L(\eta_{R})=1$, $M_{\eta}$ is allowed while $L$-conservation
implies $y_{\alpha}^{\prime}$ to be very small. Depending on the lepton charge
assignment of the $\chi_{L,R}$ fields, the different terms in the full
Lagrangian are $L$-violating, in addition to being $U(1)_{D}$ violating, and
will be either forbidden or can be taken to be naturally small, if the L
symmetry is just approximate. A full analysis of the different cases is beyond
the scope of the current discussion. We highlight one specific case which is
of particular interest: the case in which $L(\chi_{L})=L(\chi_{R})=1$ 555The
case $L(\chi_{L})=L(\chi_{R})=-1$ is equivalent.. This choice is compatible
with the $C$-symmetry discussed earlier in order to avoid diagonal dark photon
vertices and forbids the terms $\Lambda^{\prime}_{R}$ and $\Lambda_{L}$. We
notice that in this case, the lightest neutrino mass is zero as it is
protected by the accidental lepton number symmetry. Small neutrino masses can
then be controlled by lepton number breaking terms either introduced directly
in the Lagrangian as technically natural, as in standard extended seesaw
models, or induced by additional scalars which take a $U(1)_{L}$-breaking vev.

##  III Model-independent limits

Our region of interest for dark photons includes $10\text{
MeV}<m_{A^{\prime}}<10$ GeV and kinetic mixing of order
$10^{-4}<\varepsilon<0.1$. In this region, it has long been known that
colliders, fixed-target, and beam dump experiments provide the best limits on
dark photons [28, 67, 68, 19, 20]. For a compilation of constraints on dark
photons, see Refs. [69, 5, 7]. A discussion of the phenomenology at colliders
is given in Ref. [70].

Semi-visible dark photons decay into visible particles and missing energy,
modifying both bounds on visible and invisible $A^{\prime}$ models. The
dominant branching ratio (BR) of semi-visible $A^{\prime}$ is into the HNFs,
which subsequently produce both visible and invisible particles. This BR
cannot be reconstructed as a visible resonance due to the missing energy, and
it also does not satisfy the criteria of missing energy searches when the
visible products are picked up by the detector. We leave a detailed discussion
of the reinterpretation of searches for missing energy to Section IV.

Visible resonance searches — In principle, resonance searches at $e^{+}e^{-}$
colliders [71, 10] and the LHC [72] can still constrain the direct decays of
$A^{\prime}$ into SM particles,
$A^{\prime}\to\ell^{+}\ell^{-},\,\pi^{+}\pi^{-},\pi^{+}\pi^{-}\pi^{0}$. These
BRs, while still present, are typically much smaller than the BRs into HNFs,
as they are of the order of $\varepsilon^{2}\alpha/\alpha_{D}$. In addition,
when $\alpha_{D}$ is large, the dark photon will be a much wider resonance,
somewhat decreasing the effectiveness of the bump hunt method. We do not show
the rescaled limits from visible searches in our plots, as they are typically
much weaker than the model-independent constraints discussed below. We come
back to the importance of visible searches in Section VI.

Constraints on kinetic mixing that are independent of the BRs of $A^{\prime}$
can be obtained from processes that are sensitive to the exchange of virtual
dark photons. Barring fine-tuning from other new-physics contributions to
these observables, the derived limits on kinetic mixing can be regarded as
model-independent. We show these constraints in Fig. 2, comparing them with
the limits on fully invisible dark photons, shown in thin purple and dotted
black lines.

Figure 2: Model independent limits on kinetic mixing $\varepsilon$ alongside
limits on fully invisible dark photons. The limits shown as gray regions are
independent of the decay channels of the dark photon. Navy colors indicate
constraints from meson decays, $\pi^{0}\to\gamma A^{\prime}$ and
$K^{+}\to\pi^{+}A^{\prime}$. Both assume $A^{\prime}$ is invisible. In dashed
curves, we show the limits from BaBar [12], NA64 [23], and our recast of
Belle-II [73]. To obtain the latter, we neglect the interference between
initial and final state radiation of $A^{\prime}$ (see text). In green, we
show the preferred region to explain $\Delta a_{\mu}^{\rm disp}$, at
$3\sigma$, and in solid orange, we show the constraints obtained from the
lattice results for $\Delta a_{\mu}^{\rm lattice}$. In dashed orange, we also
show the region preferred by $\Delta a_{\mu}^{\rm lattice}$ at $1\sigma$.

Deep-inelastic scattering — A dark photon contributes to the deep-inelastic-
scattering (DIS) of charged leptons on nuclei via t-channel exchange,
impacting the extracted values of Parton distribution functions (PDF) [74, 75,
76, 77]. The authors of Ref. [74] set a limit for the first time, fixing the
PDF to the best-fit values of the HERA measurement, finding
$\varepsilon\lesssim 0.015$ for $m_{A^{\prime}}\lesssim 2$ GeV at $95\%$ C.L.
These limits are slightly relaxed when accounting for the effect of new
physics on the extraction of PDFs, as discussed in Refs. [78, 75]. In this
work, we use the limits of Ref. [75], where $\varepsilon\lesssim 0.034$ for
$m_{A^{\prime}}<1$ GeV at $95\%$ C.L.

Electroweak precision observables — Among the Electroweak precision
observables (EWPO) modified by kinetic mixing, the most important is
$M_{Z}^{2}\sim M_{Z^{0}}^{2}-\varepsilon^{2}M_{A^{\prime}}$ and the
corresponding shift in the mass of the $W$ boson. In Ref. [79], the global fit
to EWPO uses the W-mass measurements of LEP, finding $\varepsilon^{2}_{\rm
EWPO}<7.3\times 10^{-4}$ at $95\%$ C.L for $M_{A^{\prime}}\ll 10$ GeV. Since
then, new $W$ mass measurements have been performed by ATLAS [80] and LHCb
[81]. In addition, a recent analysis by the CDF collaboration reported a
significant deviation from the previous measurements [82]. In view of these
discrepancies with the SM, and the fact that light dark photons decrease
$M_{W}$ in the EW fit, we proceed by showing the limits from Ref. [79] with
the caveat that limits may turn to regions of preference, depending on future
developments with the $W$-mass measurement.

Meson decays — Direct production of the dark photon in meson decays provides
robust constraints on various dark photon models [28]. For invisible dark
photons, there are searches for $\pi^{0}\to\gamma A^{\prime}$ [83] and
$K\to\pi A^{\prime}$, with $A^{\prime}$ invisible [84]. We update the latter
using the latest NA62 measurement of $K\to\pi\nu\overline{\nu}$ [85],
including also the dedicated search for $\pi^{0}\to{\rm inv}$ in Ref. [86].
Since these limits assume the new vector to be invisible, they would be
modified in the semi-visible models of interest, especially in hermetic
detectors like NA62. In the rest of the paper, we follow the aggressive
strategy of showing these limits without modifications in all our plots.

Electron $(g-2)$ — Precision measurements of the electron anomalous magnetic
moment provide model-independent limits on $\varepsilon$ due to the exchange
of virtual dark photons. The dark photon contribution, in this case, comes
with a negative sign and acts to decrease $a_{e}\equiv(g-2)_{e}/2$. The two
most recent measurements of $a_{e}$ include the one in 2008 [87] and a recent
update with 2.2 times more precision in 2022 [88]. To use these results to
constrain new physics, it is necessary to compare them with high-precision SM
predictions [89]. The predictions, however, are not robust due to the
inconsistencies in the experimental determination of the fine structure
constant, $\alpha$. A group at Berkeley [90] measures the fine structure
constant using cesium-133 atoms to 120 parts per million. Another technique
employed by a group in Paris, referred here to as LKM, measures $\alpha$ to 81
parts per million [91]. These measurements are in disagreement at more than
$5~{}\sigma$, indicating that more experimental progress is needed for a
meaningful constraint to be derived. Hereafter, we follow the conservative
approach of quoting the most stringent limits, provided by Ref. [90], to set
In Fig. 2, we show both limits, as well as the region of preference that would
explain the measurement of Ref. [91]. In general, these constraints exclude
the $\Delta a_{\mu}$-explanation for dark photon masses below
$m_{A^{\prime}}\sim 30$ MeV (see below).

Muon $(g-2)$ — In the case of the muon anomalous magnetic moment, the
theoretical and experimental uncertainties on $a_{\mu}\equiv(g-2)_{\mu}/2$ are
much larger than the uncertainty in $\alpha$. Therefore, it is not subject to
the ambiguities discussed above and can still be more sensitive to new physics
due to the $a_{\mu}/a_{e}\sim m_{\mu}^{2}/m_{e}^{2}$ enhancement. The
theoretical predictions for $a_{\mu}$ in the SM (see Refs. [92, 93, 94] for a
review) have differed from the experimental measurement at E821 at the
Brookhaven National Laboratory (BNL) with a significance of 3.7$\sigma$ [95].
The E989 experiment at FNAL [96], running since 2018, has now confirmed the
central value measured at BNL, reporting $a_{\mu}^{\rm
FNAL}=116\,592\,040\,(54)\times 10^{-11}$, where the error in parenthesis is a
sum in quadrature of systematical and statistical errors. This result combined
with the BNL measurements, $a_{\mu}^{\rm BNL}=116\,592\,920\,(63)\times
10^{-11}$, provides an experimental average

$a_{\mu}^{\rm comb}=116\,592\,061(41)\times 10^{-11}.$ (31)

Eventually, with the five data-taking stages at FNAL completed, the FNAL
measurement can achieve 20 times more statistics than the BNL experiment [96]
and improve the precision on $a_{\mu}$ by a factor of 4. Eventually, this
value can also be tested by next-generation experiments, such as at the J-PARC
muon facility [97], which plans to achieve a similar precision to the BNL
measurement using a complementary technique with lower muon momenta,
$p_{\mu}\sim 300$ MeV.

The SM predictions are obtained by combining QED corrections up to 5 loops
[98, 89], electroweak ($a_{\mu}^{\rm EW}$) corrections up to 2 loops [99,
100], and hadronic contributions including the hadronic vacuum polarization
($a_{\mu}^{\rm HVP}$) [101, 102, 103, 104, 105, 106] and hadronic light-by-
light scattering ($a_{\mu}^{\rm HLbL}$) diagrams [107, 108, 109, 110, 111,
112, 113, 114, 115, 116, 117, 118, 119, 120, 121], with $a_{\mu}^{\rm HVP}$
dominating the uncertainty in the overall prediction and $a_{\mu}^{\rm HLbL}$
well below the value needed to explain the discrepancy.

Using a data-driven dispersive calculation for the hadronic contributions
$a_{\mu}^{\rm HVP}$ and $a_{\mu}^{\rm HLbL}$, the authors of Ref. [94]
converge on the prediction

$a_{\mu}^{\rm disp}=116\,591\,810\,(43)\times 10^{-11}.$ (32)

If the dispersive value holds, the tension between experiment and theory
reaches $4.2\sigma$, where $\Delta a_{\mu}^{\rm disp}=a_{\mu}^{\rm
disp}-a_{\mu}^{\rm comb}$ is

$\Delta a_{\mu}^{\rm disp}=251(59)\times 10^{-11}.$ (33)

A strong ongoing effort aims at reducing the dominant source of uncertainties
in these hadronic contributions, but so far, it has not been demonstrated that
hadronic uncertainties alone can reconcile the discrepancy [122, 123, 124].
The agreement between EW precision measurements and the $e^{+}e^{-}\to$
hadrons data corroborates this hypothesis. A deviation in
$\sigma(e^{+}e^{-}\to\text{hadrons},s)$ to explain the observed value of
$\Delta a_{\mu}$ was shown to be ruled out for $e^{+}e^{-}$ center-of-mass
energies of $\sqrt{s}>0.7$ GeV [125], implying that any missed contributions
ought to be mostly concentrated in the $\pi^{+}\pi^{-}$ region [126]. It has
also been suggested that new physics could be hiding in this data [127, 128,
129, 130].

The discrepancy of Eq. 33 is still not conclusive, however. In particular, a
calculation of $a_{\mu}^{\rm HVP}$ on the lattice by the BMW collaboration
[131] finds a $2.1\sigma$ significant disagreement with the value obtained
using the data-driven dispersive method of Ref. [94]. Using the BMW result for
$a_{\mu}^{\rm HVP}$ reduces the disagreement between theory and the experiment
average, $a_{\mu}^{\rm comb}$, down to $1.5\sigma$,

$\Delta a_{\mu}^{\rm lattice}=107(69)\times 10^{-11}.$ (34)

This lattice result has been increasingly scrutinized in the search for
additional systematic uncertainties that are specific to discretization and
finite-volume effects of the lattice. Consistency checks of the BMW results
have been performed by other collaborations using “Euclidean window
observables”, namely, observables calculated in Euclidean time windows that
enhance or suppress specific systematic uncertainties [132]. One of these
observables isolates contributions to $a_{\mu}^{\rm HVP}$ into short,
$a_{\mu}^{\rm SD}$, and intermediate, $a_{\mu}^{W}$, time-distance pieces. The
data-driven dispersive method [133] and lattice calculations of $a_{\mu}^{\rm
SD}$ are in good agreement. However, a $3.8\sigma$ significant tension exists
between the intermediate time-distance observable, $a_{\mu}^{\rm W}$, and all
lattice results [131, 132, 134, 135, 136], suggesting the disagreement with
$e^{+}e^{-}\to$ Hadrons data is, in fact, largest in the energy region of
$\sqrt{s}\sim 1-3$ GeV [133]. While the nature of this discrepancy is not
identified, in this study, we proceed to entertain BSM explanations to both
$\Delta a_{\mu}^{\rm disp}$ and $\Delta a_{\mu}^{\rm lattice}$. The $3\sigma$
preference region for $\Delta a_{\mu}^{\rm disp}$ is shown as a green band in
Fig. 2, alongside the $1\sigma$ preference band and the $3\sigma$ exclusion
limit from $\Delta a_{\mu}^{\rm lattice}$. We summarize other new-physics
explanations to $\Delta a_{\mu}^{\rm disp}$ in Appendix B.

Belle II — As of now, no dedicated search for invisible dark photons has been
released by the Belle II $e^{+}e^{-}$ collider. Nevertheless, the
collaboration has performed a search for $L_{\mu}-L_{\tau}$ gauge bosons,
focusing on final state radiation (FSR) of dark photons,
$e^{+}e^{-}\to\mu^{+}\mu^{-}Z^{\prime}_{\mu-\tau}$. The search is based on the
missing energy carried by the gauge boson. A similar signature can take place
for dark photons, with the difference that $A^{\prime}$ can be emitted as
either initial state radiation (ISR) or FSR. In the analogous QED processes,
$e^{+}e^{-}\to\mu^{+}\mu^{-}\gamma$, the interference between ISR and FSR can
lead to significant charge asymmetries [137], so the differential event rate
will be different for a dark photon. This interference in QED vanishes,
however, if integrated over the total phase space. While a dedicated study of
the efficiencies is needed to recast the Belle-II limit on kinetic mixing, we
provide a simple estimate by neglecting the interference term and simply
adding the ISR and FSR pieces together. The estimated limit is shown in Fig. 2
with an asterisk to emphasize the approximation in the recast method. The
semi-visible decay of the dark photon can also lead to additional energy in
the final state, and relax this constraint. Therefore, we show it as a dashed
line in Fig. 2. Because of the similar geometries of BaBar and Belle II, we do
not include this limit in our recast analyses: a strong relaxation of BaBar
would lead to a similar effect in Belle II, which is, in addition, a more
hermetic detector.

##  IV Reinterpreting constraints on invisible $A^{\prime}$

Searches for invisibly decaying dark photons at $e^{+}e^{-}$ colliders and
fixed-target experiments provide the strongest constraints on models of semi-
visible dark photons. At collider experiments, the dark photon is produced
directly alongside initial state radiation (ISR) in the $e^{+}e^{-}\to\gamma
A^{\prime}$ process. Prompt decays of $A^{\prime}$ to a pair of HNFs,
$A^{\prime}\to\psi_{i}\psi_{j}$, in which the HNFs are:

  1. 1.

long-lived and decay outside the detector, or

  2. 2.

short-lived and decay inside to $\psi_{i}\to\psi_{j}\,\ell^{+}\ell^{-}$, with
final state leptons pairs whose energies fall below detector thresholds,

would appear identical to a monophoton signature accompanied by missing
$\not{E}_{T}$. The strongest bound of this kind is obtained by the BaBar
experiment, excluding explanations of $(g-2)_{\mu}$ with invisibly decaying
dark photons [12]. Dark photons can also be produced in bremsstrahlung at
fixed-target experiments such as NA64. [13, 23, 24]. In this type of
experiment, the dark photon signature constitutes a large amount of missing
energy in the primary electron beam that scatters on the target. Below we
discuss the reinterpretation of these two leading constraints on the parameter
space of the models discussed in Section II. In the following, the vector
$\bm{x}$ represents the set of free model parameters,
$\bm{x}=(\varepsilon,\,m_{A^{\prime}},\,\Delta_{21},\,\dots)$, that have been
varied in the analysis performed.

###  IV.1 BaBar monophoton search

Figure 3:  The signatures of semi-visible dark photons at $e^{+}e^{-}$
colliders. On the right, an inner view of the BaBar detector with the
displaced, semi-visible decay of the HNFs into charged lepton pairs. In this
example, the decay cascade is
$A^{\prime}\to(\psi_{2}\to\psi_{1}\mu^{+}\mu^{-})(\psi_{3}\to\psi_{2}e^{+}e^{-}\to\psi_{1}e^{+}e^{-}e^{+}e^{-})$,
where $\psi_{3}$ decayed promptly.

Based on the PEP-II asymmetric $e^{+}e^{-}$ collider at the SLAC National
Accelerator Laboratory, the BaBar experiment searched for single photons
(monophotons) accompanied by missing energy and momentum in the process
$e^{+}e^{-}\to\gamma A^{\prime}$. The search was conducted in the $53$ fb-1
dataset collected between 2007-2008 at the center-of-mass (COM) energies
$\Upsilon(2S),\Upsilon(3S)$ and $\Upsilon(4S)$. The components of the BaBar
detector relevant for our analysis are a charged-particle-tracking system
provided by a five-layer, double-sided silicon vertex tracker (SVT) and
$40$-layer drift chamber (DCH); an electromagnetic calorimeter (EMC) of $6580$
CsI(Tl) crystals. These systems are all contained within a $1.5$ T
superconducting solenoid magnet. Beyond the superconducting coil is located an
instrumented flux return (IFR) barrel that provides muon and neutral hadron
identification. An illustration of the detector and the HNF signatures is
shown in Fig. 3.

Event generation — Using our own MC event generator, we simulate the
production of dark photons at BaBar in the process $e^{+}e^{-}\to\gamma
A^{\prime}$ at a COM energy $\sqrt{s}\approx 10$ GeV. We boost and rotate to
the laboratory frame, considering that the $e^{+}e^{-}$ collision frame is
itself already boosted with respect to the laboratory by $\beta_{z}\approx
0.5$. In the laboratory frame, we take the z-direction to be aligned with the
direction of $e^{+}e^{-}$ collision. The experiment employs a primary
selection cut on the photon COM angle to suppress SM backgrounds,
$|\cos\theta_{\gamma}^{*}|<0.6$. To a good approximation, the $A^{\prime}$
decay to HNF pairs is prompt upon production. However, the HNF can travel
before decaying, modeled by random sampling an exponential distribution
according to its decay width. We simulate the decays $A^{\prime}\to
e^{+}e^{-},\mu^{+}\mu^{-}$, and $\pi^{+}\pi^{-}$ according to their
differential decay rates, taking into account the differences in the decay
kinematics of Majorana and Dirac particles.

As the original analysis searched for single photons and vetoed additional
activity in the detector above a certain energy, we introduced a set of veto
criteria to dispense those events that would not have passed the event
selection. We show kinematical distributions of $e^{+}$ and $e^{-}$ decay
products in Appendix A, also showing the impact of the analysis selection
criteria discussed next.

Figure 4:  The BaBar limit on the kinetic mixing parameter, $\varepsilon$. On
the left panel, we show the limit as a function of our analysis’s individual
$e^{\pm}$ energy threshold. In solid (dashed) lines, we use an analysis with
(without) a cut on the angle between the leptons and the pipeline. On the
right, we show the limit as a function of the mass splitting for BP1b.

Monophoton selection — An $e^{+}e^{-}\to\gamma A^{\prime}$ event passing the
initial monophoton selection is vetoed if, anywhere along the decay chain, a
charged particle is produced in an instrumented region of the detector, i.e.,
within the SVT, DCH, EMC, or IFR regions, and satisfies both of the following
conditions:

  1. 1.

For $e^{+}e^{-}$ pairs with angular separation
$\theta_{\text{sep.}}>10^{\circ}$, the energy of each electron exceeds BaBar’s
energy threshold to resolve charged particle tracks, $E_{\pm}>100$ MeV. For
overlapping $e^{+}e^{-}$ pairs with $\theta_{\text{sep.}}<10^{\circ}$, we
require $(E_{+}+E_{-})>100$ MeV.

  2. 2.

The polar angles, $\theta_{\text{pol.}}$, of the electrons individually, or as
a pair, are sufficiently wide that the electrons do not escape along the beam
pipeline, $17^{\circ}<\theta_{\text{pol.}}<142^{\circ}$.

The criteria above amount to a statement that all HNF decays that occur inside
the detector and produce charged lepton final states that leave visible tracks
are vetoed.

The threshold is assumed to be a step function with $100\%$ detection
efficiency above $E_{\rm threshold}$ and $0\%$ below it. Realistically, the
final state leptons can escape detection even if their energy is large, as
leptons can escape between the active materials in the detector. This effect
requires a more detailed description of the geometry and particle propagation
model, which is beyond the capabilities of our simulation. Nevertheless, we
expect this will not significantly change our conclusions, as the leptons are
always produced in pairs and follow bent trajectories due to the magnetic
field, especially at low energies.

We show the impact of varying energy thresholds used in our analysis on the
left panel of Fig. 4. We do this for a few benchmark points, demonstrating the
dependence on the threshold assumptions and confirming that the dominant
source of invisible events at BaBar comes from soft leptons. The effect of
omitting the pipeline is small, as seen in the comparison between the solid
lines (considering the pipeline effect) and dashed lines (neglecting it). In
addition, in varying the mass splitting, we vary the total energy emitted in
SM particles, observing a strong effect on the relaxation of constraints for
larger $\Delta_{21}$. The band in each curve represents the uncertainty
associated with Monte-Carlo statistics.

Recasting the bound — To derive their bound, BaBar assumed an invisibly
decaying dark photon

$\sigma^{\text{BaBar}}_{e^{+}e^{-}\to\gamma+\text{inv.}}\left(\bm{x}\right)=\sigma_{e^{+}e^{-}\to\gamma
A^{\prime}}\left(\bm{x}\right)\times\mathcal{B}\left(A^{\prime}\to\text{inv.}\right),$
(35)

with $\mathcal{B}\left(A^{\prime}\to\text{inv.}\right)=1$. To re-interpret the
bound for a semi-visible dark photon, we introduce a factor
$\text{P}^{\text{inv.}}$ that accounts for the probability that decays of
semi-visible dark photons produced alongside ISR appear invisible and
contribute to the monophoton dataset:

$\sigma^{\text{BaBar}}_{e^{+}e^{-}\to\gamma+\text{inv.}}\left(\bm{x}\right)=\sigma_{e^{+}e^{-}\to\gamma
X}\left(\bm{x}\right)\times\text{P}^{\text{inv.}}\left(\bm{x}\right),$ (36)

where $X$ is the semi-visible dark photon. From Eq. 35 and Eq. 36, we may
obtain the relation

$\varepsilon^{\text{BaBar}}=\varepsilon\,g_{D}\sqrt{\text{P}^{\text{inv.}}\left(\bm{x}\right)},$
(37)

where $\varepsilon^{\text{BaBar}}$ is the bound on the kinetic mixing
parameter obtained by BaBar. We may define the function
$\text{P}^{\text{inv.}}$ as follows

$\text{P}^{\text{inv.}}\left(\bm{x}\right)=1-\text{P}^{\text{veto}}\left(\bm{x}\right)=1-\frac{N^{\text{veto}}\left(\bm{x}\right)}{N\left(\bm{x}\right)},$
(38)

with $N$ being the total number of events that pass the initial monophoton
selection, and $N^{\text{veto}}$ the subset of events that are vetoed
according to the criteria set out above. To recast BaBar’s monophoton limit,
we solve Eq. 37 for $\varepsilon$ at each value of $\left(\bm{x}\right)$. The
function $\text{P}^{\text{inv.}}$ contains all model dependencies, including
the HNF masses and any mixing parameters. The results of our recast are given
in Section V for all benchmark points of interest.

Pseudo-monophotons — It was noticed that BaBar has a mild excess of mono-
photon events [37]. If one of the two HNFs decays with a lifetime of a few cm,
a significant fraction of them will decay within the EMC of BaBar, mimicking a
monophoton signature. These events have a relatively broad spectrum in missing
energy $M_{\textrm{miss}}^{2}\equiv s-2E_{\gamma}\sqrt{s}$ where a mild excess
has been observed in the region
$24~{}\mathrm{GeV}^{2}<M_{\textup{miss}}^{2}<50~{}\mathrm{GeV}^{2}$. This
explanation can fit the events well, explaining the $\sim 2.5\sigma$ excess
[37].

###  IV.2 NA64 dark photon searches

Figure 5: The NA64 setup and signatures considered in this work. Left panel:
production and decay processes of the heavy neutral fermions (HNF). Top right
panel: the first kind of signature with a displaced vertex already considered
in [24]. Bottom right panel: the prompt decay signature that NA64 can use to
constrain regions of parameter space where the HNFs are promptly decaying.

The fixed-target experiment NA64 searches for dark sector particles at the
CERN SPS, employing a $100~{}\mathrm{GeV}$ electron beam. The main search
strategy relies on the fact that invisible dark photons can carry away a large
amount of missing energy in hard electron-nucleus bremsstrahlung events [18,
20, 22, 13, 23],

$e^{-}\mathcal{N}\to
e^{-}\mathcal{N}A^{\prime},\,A^{\prime}\to~{}\mathrm{invisible}\,,$ (39)

where $\mathcal{N}$ is the target nucleus in the fixed target. Like monophoton
searches, the bremsstrahlung signal is proportional to $\varepsilon^{2}$.

The main parts of the detector relevant to our analysis are:

  * •

the electromagnetic calorimeter (ECAL), which also works as the active beam
dump, made of Pb+Sc layers, with an average photon conversion length of
$X_{0}=1.175~{}\mathrm{cm}$;

  * •

a large high-efficiency veto counter downstream of the ECAL;

  * •

a hadronic calorimeter (HCAL), made of three different modules.

The search for invisible decays of $A^{\prime}$ [23] was conducted on the
total sample of electrons on target (EOT) collected during the period
2016–2018, $n_{\textup{EOT}}=2.84\times 10^{11}$.

NA64 also searched for semi-visible dark photon decays assuming the iDM model
[24]. In this search, they considered the same data collected in the period
2016–2018, performing a recast-based analysis resembling that of their search
for axion-like particles [138], targeting visible final states coming from the
dark photon decay chain and putting a model-dependent constraint on the
kinetic mixing parameter.

Event generation — We simulated the production of dark photons and the
detection of $e^{+}e^{-}$ pairs in the NA64 detector with a fast MC generator.
From a comparison of the sensitivity curves with the limits obtained in Ref.
[38], we find that our projections are comparable to the limits obtained using
a proper GEANT4 detector simulation, with minor discrepancies that can be
attributed to the complete description of the detector geometry, shower
development, and energy collection within the different detectors. The
complete GEANT4 detector simulation, along with a discussion of the latest
NA64 constraints on several semi-visible dark photon models, can be found in
the accompanying [38].

We simulate bremsstrahlung events, producing an $A^{\prime}$ with an electron
beam at energy $E_{\textup{beam}}=100~{}\mathrm{GeV}$. We consider the
electron beam energy to be unaffected by any energy losses happening when
entering the ECAL, so that $E_{\textup{beam}}$ can be considered constant. The
beam is considered to impact the ECAL at coordinate $x=y=z=0$. The energy of
the $A^{\prime}$ is distributed according to the following formula, obtained
by applying the improved Weizsaker–Williams (IWW) approximation [139]:

$\frac{\mathrm{d}\sigma}{\mathrm{d}x}\propto\biggl{(}m_{A^{\prime}}^{2}\frac{1-x}{x}+m_{e}^{2}x\biggr{)}^{-1}\biggl{(}1-x+\frac{x^{2}}{3}\biggr{)}\,,$
(40)

where $x=E_{A^{\prime}}/E_{\textup{beam}}$, and the $A^{\prime}$ is in the $z$
direction. After radiating an $A^{\prime}$, the beam electron (or main
electron) will have energy $E_{e}=E_{\textup{beam}}-E_{A^{\prime}}$, and will
shower inside the ECAL completely. The $A^{\prime}$ decays promptly into a
pair of HNFs, which are boosted and rotated to the lab frame according to the
$A^{\prime}$ energy. Each HNF will then decay according to its decay modes.
The simulation automatically handles the decay of the secondary HNFs in the
same way. The lightest HNF from each model is considered stable with respect
to the size of the detector. We assume a simplified shower development for the
$e^{+}e^{-}$ pairs produced inside the ECAL. The energy loss can be computed
assuming an exponential law:

$\begin{split}\frac{\mathrm{d}E}{\mathrm{d}z}=-\frac{E}{X_{0}}\Rightarrow\Delta
E(z)&=E(z)-E_{0}\\\
{}&=E_{0}\biggl{[}1-\exp\biggl(-\frac{z-z_{0}}{X_{0}}\biggr{missing})\biggr{]}\,,\end{split}$
(41)

where $z_{0}$ is the production point of the pair. For pairs detected inside
the HCAL, we assumed each pair is able to shower completely inside it. Given
the high energies of the $A^{\prime}$, the $e^{+}e^{-}$ pairs are highly
boosted and collimated. We may then treat each pair as a single particle with
energy $E(z_{0})=E_{e^{+}}(z_{0})+E_{e^{-}}(z_{0})$. The kinematical
considerations on the $e^{+}e^{-}$ pairs are corroborated by the distributions
shown in Appendix A.

Each event is made of invisible final states (the stable HNFs) and visible
energy (the $e^{+}e^{-}$ pairs):

$\displaystyle e^{-}\mathcal{N}\to e^{-}\mathcal{N}A^{\prime}$ $\displaystyle
A^{\prime}\to(\psi_{i}\to(\psi_{k}\to\dots)\,e^{+}e^{-})(\psi_{j}\to(\psi_{l}\to\dots)\,e^{+}e^{-})$

We check the visible energy collected in each event against the veto criteria
applied to the different regions of the detector to see if the event is
recorded as signal for the experiment.

Semi-visible selection (S1) — This signature corresponds to the one employed
by NA64 to study semi-visible dark photon decays in the framework of the iDM
model [24]. In this work, we recast their limits to the more complex models
presented in Section II.

The following selection criteria were applied, according to the expected
signal yield coming from this model, relying on the decay of the heaviest HNF
$\psi_{i}=\psi_{2}$ or $\psi_{3}$ into a single $e^{+}e^{-}$ pair:

  1. 1.

$\psi_{i}$ is expected to decay inside the HCAL: in particular, the analysis
targets a fiducial HCAL volume composed by the last two modules of the HCAL
detector in which a consistent amount of energy coming from a single
$e^{+}e^{-}$ pair is expected to be detected.

  2. 2.

events with any other activity happening before the fiducial HCAL volume are
vetoed;

  3. 3.

events with more than one visible decay vertex in the fiducial HCAL volume are
vetoed.

Additionally, in case $\psi_{i}$ decays beyond the HCAL, no other significant
energy deposits are expected in the full detector, and the event resembles an
invisible dark photon event.

This proposed analysis is tailored to the iDM model, focusing in particular in
the parameter space where the $\psi_{2}$ lifetime is comparable to the size of
the detector. It rapidly loses sensitivity in the high mass region where
$\psi_{2}$ decays promptly. Other models involving large mass splittings
between HNFs are also characterized by a shorter lifetime, due to
$c\tau\propto\Delta^{-5}$, so this search will also be less sensitive to the
high mass region.

Invisible selection (S2) — The analysis focused on detecting invisible decays
of $A^{\prime}$ by constraining the amount of visible energy collected. The
set of energy cuts relies on the search strategy applied by the NA64
experiment to detect dark photon events containing missing energy, while
limiting the possible background [23]:

  1. 1.

pre-shower ECAL. The total energy collected in the first layers of the ECAL
should be compatible with the deposit expected for a primary electron.

  2. 2.

ECAL. The total energy collected, $E_{\textup{ECAL}}$, including both the main
electron and the $e^{+}e^{-}$ pairs, is compared against the missing energy
threshold chosen by the experiment:

$E_{\textup{ECAL}}>50~{}\mathrm{GeV}\,.$ (42)

  3. 3.

Veto counter. We do not expect any activity in the veto, in order to record an
invisible event. In the case of semi-visible dark photon models, a particle
can reach the veto in two cases:

     * •

if they are produced inside the ECAL, they can shower until its end, releasing
the remaining energy in the veto;

     * •

if they are produced between the end of the ECAL and the veto, they will
release their energy inside the veto.

We further assume that if a pair happens to reach the veto, it is able to
release all of its energy inside it. This can be justified by the fact that
the imposed veto threshold is sufficiently low that even the softest
$e^{+}e^{-}$ pairs we expect to produce would be able to trigger an event.
Moreover, the thickness of the veto counter is large enough to guarantee a
consistent energy deposit by the $e^{+}e^{-}$ pairs.

  4. 4.

HCAL: For an invisible event, no activity is expected inside the HCAL.
Particles created between the veto and the HCAL or in the empty space between
the three HCAL modules will be intercepted by the HCAL, eventually. The HCAL
detector of NA64 is sufficiently long so that we can assume that any pair
created inside it has enough space to shower completely, depositing its entire
energy inside the HCAL. This approximation may not apply only to the small
number of particles created at the very end of the HCAL.

Using our simulation, we present a few distributions of the kinematics of the
$e^{+}e^{-}$ pairs in Appendix A. It should be noted that in performing the
recast of the invisible selection, we have assumed that NA64 is not able to
distinguish the $e^{+}e^{-}$ showers due to short-lived HNF decays in the ECAL
from the main electron beam. We assume that the total energy deposition in the
ECAL for these cases is the aggregate of the $e^{+}e^{-}$ energy and beam
electron energy. Given the NA64 sensitivity to these prompt decays has not
been previously studied, we take this recast constraint as a projection of
NA64’s potential sensitivity to decays of this kind. The latest constraints
from NA64, Ref. [38], tackle this region and show good agreement with our
projections.

Recasting the bound — The derivation of the recast bound is done in a similar
fashion to the BaBar simulation. We start by considering the bound obtained by
the invisible dark photon search performed in [23], which we call
$\varepsilon^{\textup{NA64}}$. In addition, we extrapolated the bound above
$1$ GeV, through a 2-degree polynomial fit. Our toy Monte-Carlo yields the
probability of obtaining an invisible event assuming the semi-visible dark
photon decay $P^{\textup{inv}}$. The recast bound can be found by solving the
following equation, which matches Eq. 37 discussed for the BaBar recast, with
$\bm{x}$ being the set of model parameters:

$\varepsilon^{\textup{NA64}}=\varepsilon\,g_{D}\sqrt{\text{P}^{\text{inv.}}\left(\bm{x}\right)}.$
(43)

###  IV.3 Other limits from beam dumps

We also recast beam dump limits on displaced decays of dark sector particles.
In particular, the NuCal and E137 limits, originally recasted on the parameter
space of iDM models [140, 29, 52]. We do not simulate the experimental setup
but rescale the upper and lower bounds under a few approximations. For the
upper bound (small coupling regime), the rescaling is done by matching our
models’ production and detection rates with the iDM model. In this region, the
decay-in-flight rate is proportional to the number of particles produced times
the decay rate into the respective signal channels, in our case,
$\psi_{i}\to\psi_{j}\ell^{+}\ell^{-}$. In some cases, more production channels
are available than in the original study, such as when recasting iDM bounds
onto the parameter space of mixed-iDM and i2DM, where $\psi_{2}$ can be
produced in pairs in meson decays or bremsstrahlung. We multiply the signal
rate by the respective $A^{\prime}$ vertices and the corresponding long-lived
particle decay rate to account for that. For instance, to rescale a limit on
$\varepsilon$ from an iDM model to a mixed-iDM or i2DM model, we take

$\frac{\varepsilon^{4}_{\rm iDM}}{\varepsilon^{4}_{\rm
new}}=\frac{\alpha_{D}^{\rm new}}{\alpha_{D}^{\rm iDM}}\sum_{i\geq
j}|V_{ij}|^{2}\frac{\Gamma_{\psi_{i}\to\psi_{i-1}e^{+}e^{-}}^{\rm
new}+\Gamma_{\psi_{j}\to\psi_{j-1}e^{+}e^{-}}^{\rm
new}}{\Gamma_{\psi_{2}\to\psi_{1}e^{+}e^{-}}^{\rm iDM}},$ (44)

where $\hat{\Gamma}$ are the decay rates deprived of all coupling pre-factors.
The $\alpha_{D}|V_{ij}|^{2}$ pre-factor takes care of rescaling the production
of different dark states at the source or target, and the decay rate ratios
account for the different mass splitting and mediator masses. Note that we
neglect decay cascades and that no cases require taking
$\psi_{i}\to\psi_{i-2}e^{+}e^{-}$ decays into account. We also neglect the
effect of the mass splitting on the kinematics, which can lead to a different
acceptance for the new models. For the lower-bound (upper line), we match the
total lifetimes of the lightest unstable particles of the original bound to
the one in our models. That is, the new limit is found by setting
$\Gamma_{\psi_{2}}^{\rm new}=\Gamma_{\psi_{2}}^{\rm iDM}$. This is also an
approximate procedure but relies on the fact that in this lower-limit regime,
the number of particles produced is large, but the probability of the unstable
state decaying within the detector is going to zero exponentially. Our
procedure also does not increase the mass reach of the constraints and
therefore provides only an estimate of the excluded parameter space. For E137
(scatter), the recast is similar to Eq. 44, but we rescale the scattering rate
of dark particles instead, neglecting differences in mass splitting.

##  V Results

We have studied a series of constraints on the theoretical models presented in
Sec. II, recasting the limits from NA64 and BaBar. As the parameter space is
very large, we have fixed representative values for the gauge coupling $g_{D}$
and the masses of the HNFs, parameterized in terms of $m_{1}/m_{A^{\prime}}$
(known as $r$ in the literature) and $\Delta_{ij}\equiv(m_{i}-m_{j})/m_{j}$.
We have taken $\alpha_{D}\equiv g_{D}^{2}/(4\pi)$ to be in the range $0.1-0.5$
as we are interested in the fast decays of the HNFs while maintaining
perturbativity. Except for the 3 HNF models, we fix the mass ratio
$r=m_{1}/m_{A^{\prime}}$ to be $1/3$ for most of the benchmarks, a standard
value in the literature. The values $\Delta_{ij}$ have been set to minimize
the lifetime of the HNFs and maximize the amount of visible energy deposited
by their decay products in the BaBar and NA64 detectors. The constraints from
NA64 are labeled according to the dark photon signature.

  * •

_NA64 (S1)_ (solid line) — described in [24], a model-dependent search for iDM
was performed for the semi-visible dark photon signature of the model.

  * •

_NA64 proj. (S2)_ (dash-dotted line) — described in [23], this constraint is
on the invisible dark photon. Our recast expresses the potential future
sensitivity of the experiment towards a dedicated semi-visible dark photon
search, showing the capability to constrain the parameter space in a model-
independent way.

In addition, for benchmarks BP1, BP2, and BP3, the lightest HNF can be a dark
matter candidate. In this case, $\Delta_{ij}$ cannot be too large, to minimize
the Boltzmann suppression in coannihilations.

We show the model-dependent bounds from NuCal [141, 142, 52], E137 (scatter)
[140, 29, 52], and the model-independent limits discussed in Section II.

Inelastic dark matter (BP1a/b) — The results for the iDM benchmark are shown
in Fig. 6, expressed in the $\varepsilon$/$m_{A^{\prime}}$-plane.

Figure 6:  The kinetic mixing $\varepsilon$ as a function of the dark photon
mass $m_{A^{\prime}}$ for BP1a (left) and BP1b (right) in the inelastic dark
matter (iDM) model. We show the $\Delta a_{\mu}$-preferred $1$, $2$, and
$3\sigma$ regions in shades of green. The recast constraints from BaBar and
NA64 are shown in blue and orange, respectively. The NA64 curves show the
recast constraints using the dedicated _semi-visible_ search, corresponding to
those derived using displaced decays (S1 in Fig. 5), and the projected
sensitivity of a search for invisible and promptly-decaying particles (S2 in
Fig. 5). Assuming the lightest HNF to be dark matter, the relic density line
is shown in black. Other constraints from $(g-2)_{e}$, EWPO, DIS, and NA62 are
shown with thin gray lines and light gray regions and are referred to as model
independent constraints (see Section III). The constraint imposed by NuCal and
E137 (scatter) are shown with the same style. The masses of vector meson
resonances are shown as vertical grey dashed lines.

In BP1a, we see a significant relaxation of the NA64 and BaBar bounds, with a
sizeable $\Delta a_{\mu}$ preference region now open. Due to a large DS
coupling, $\alpha_{D}=0.5$, the decay rates of the HNFs are enhanced, allowing
for more semi-visible events in the detectors. Benchmark BP1b corresponds to
the choice of parameters used in Ref. [29]. Assuming the lightest HNL,
$\psi_{1}$, is a dark matter candidate, we find that the correct dark matter
abundance can be achieved in both benchmarks BP1a and BP1b.

In BP1b we find a much less significant relaxation of the bounds compared to
BP1a, leaving very little open parameter space for a $\Delta a_{\mu}$
explanation. In particular, we find the BaBar bound to be much more
constraining than in Ref. [29]. This is predominantly due to the difference in
selection criteria used in the two analyses. In Ref. [29], an energy cut of
$60$ MeV is applied to charged particles produced in semi-visible $A^{\prime}$
decays. In this work, we take the larger value of $100$ MeV, corresponding to
the energy threshold used in the analyses to veto additional tracks in the
BaBar drift chamber [143, 12]. The impact of the higher energy threshold is to
veto a greater proportion of the semi-visible decays, leading to a stronger
constraint from BaBar. In addition, we cut on the polar angle of the charged
lepton pairs, which allows us to exclude events in which the leptons are
produced in the direction of the beam pipeline. We find the fraction of these
“pipeline” events to be small and that the relative strength of the BaBar
bound is predominantly a consequence of the energy threshold. Contrary to the
more conservative analysis of Ref. [30], in which a threshold energy of $150$
MeV is taken, a very small region of preference for $\Delta a_{\mu}$ remains
open at the $2\sigma$ level for $m_{A^{\prime}}\sim 300-500$ MeV.

In both BP1a and BP1b, we find that a search for promptly decaying HNFs at
NA64, with signatures of the type S2, can cover the newly open $\Delta
a_{\mu}$ parameter space (see also the companion paper in Ref. [38]). This
region of parameter space is also accessible to other lower-energy
$e^{+}e^{-}$ colliders, including KLOE/KLOE-2 [144, 145] and BES-III [146].

In the BaBar signal selection, events with charged tracks above
$100~{}\mathrm{MeV}$ were vetoed, justifying our choice for the energy
threshold [12]. However, the probability of detection of such low-energy
tracks may not be exactly unity, setting an effective energy threshold.
Therefore, our hard requirement on the energy of leptons is an important
source of uncertainty. In Fig. 7, we show the effects of changing this value
on the BaBar limits for the iDM model in BP1a and BP1b.

Figure 7: Same as Fig. 6, with the addition of multiple BaBar recasts obtained
under different assumptions for the energy threshold for the $e^{+}/e^{-}$
tracks. The different energy thresholds are represented by shades of blue with
the adopted values shown in the legend. The threshold chosen for the analysis
performed in this work is marked with a $\star$ symbol.

Increasing the energy threshold makes it more difficult to veto events,
leading to a stronger bound. For a threshold larger than $250~{}\mathrm{MeV}$,
BaBar would cover the entire parameter space allowed by $\Delta a_{\mu}$ in
the case of BP1a, while a threshold of $150~{}\mathrm{MeV}$ is already
sufficient to rule it out completely in BP1b. We note, however, that
thresholds exceeding $200~{}\mathrm{MeV}$ are likely unrealistic as the
probability of missing such tracks within the inner BaBar detectors is very
small.

To understand whether it is possible to accommodate the different constraints
and the relic density, we can inspect the parameter space in
$\alpha_{D}$/$m_{A^{\prime}}$ and $\Delta_{21}$/$m_{A^{\prime}}$ plots, shown
respectively in Fig. 8 and in Fig. 9. In Fig. 8, we fix the value of
$\varepsilon$ to that required to explain $\Delta a_{\mu}$ and we vary
$\alpha_{D}$ and $m_{A^{\prime}}$. Increasing $\alpha_{D}$ affects the
lifetime of $\psi_{2}$, making them short-lived and allowing their decays to
happen inside the detector, increasing the amount of semi-visible events that
can be identified. Both scenarios are strongly constrained by BaBar, NA64,
E137, and NuCal, with only a small part of the parameter space allowed in the
left panel. That region is also able to accommodate the DM relic density. The
two scenarios differ only by the value of the mass splitting $\Delta_{21}$:
the right panel is characterized by a lower value. Decreasing $\Delta_{21}$
means that the two HNFs become more degenerate and the lifetime of $\psi_{2}$
increases, becoming larger than the size of the detector. This effect reduces
the amount of $\psi_{2}$ decays happening inside the detector and makes the
bound more constraining. The right panel is indeed entirely constrained.

Figure 8:  Parameter space of $\alpha_{D}$ versus $m_{A^{\prime}}$ for the iDM
model, fixing to the value needed to explain the $\Delta a_{\mu}$,
$\varepsilon=\varepsilon_{g-2}$. The mass splitting is fixed at
$\Delta_{21}=0.5$ (left) and $\Delta_{21}=0.4$ (right). The constraints are
shown with the same style used in Fig. 6.

Figure 9:  Parameter space of $\Delta_{21}$ versus $m_{A^{\prime}}$ for the
iDM model, fixing to the value needed to explain the $\Delta a_{\mu}$,
$\varepsilon=\varepsilon_{g-2}$. The dark coupling is fixed at $\alpha=0.5$
(left) and $\alpha=0.1$ (right). The constraints are shown with the same style
used in Fig. 6.

In Fig. 9 we fix the value of $\varepsilon$ to that required to explain
$\Delta a_{\mu}$ and we vary $\Delta_{21}$ and $m_{A^{\prime}}$. In this case,
the E137 constraint shown with a thin gray line has been extrapolated at large
$\Delta_{21}$ value. As discussed, previously, decreasing $\Delta_{21}$ makes
most of the $\psi_{2}$ decays to happen outside of the detector. No semi-
visible event would be detected in this case, and the bound would resemble the
original invisible $A^{\prime}$ bound, covering the region in which
$\varepsilon$ can explain $(g-2)_{\mu}$. In this case, the NA64 (S1)
constraint follows a similar trend as for the BaBar constraint, becoming
weaker for larger mass splitting. In the case of NA64 (S2), at large
$m_{A^{\prime}}$, the experiment loses sensitivity, and the bound becomes
naturally weaker, independently of the value of $\Delta_{21}$. This
corresponds to the loss of sensitivity in the original invisible bound posed
by NA64, caused by a lack of event rate.

Nevertheless, increasing the mass splitting between the two HNFs to
accommodate the constraints affects the dark matter relic abundance of
$\psi_{1}$. A larger mass splitting increases the Boltzmann suppression of
$\psi_{2}$ number density in the early universe, depleting it faster and
suppressing the coannihilation contribution to the cross section. This results
in an overabundant scenario, which can be controlled only assuming secluded
annihilations within the dark sector, and it is expressed by the parameter
space above the relic density line in Fig. 9. The smaller $\alpha_{D}$ value
in BP1b translates into a shift of the relic density line towards lower
$\Delta_{21}$, because of its effect in decreasing the annihilation cross
section. A smaller $\Delta_{21}$ ensures a smaller Boltzmann suppression of
$\psi_{2}$ number density and a larger coannihilation cross section.

It is interesting to notice that the projections shown by the NA64 (S2) line
could constrain the free parameter space. The search for promptly decaying
HNFs in the detector can address whether this minimal model can simultaneously
explain $\Delta a_{\mu}$ and dark matter.

Mixed Inelastic Dark Matter (BP2a/b) — We show the constraints for the mixed-
iDM model in Fig. 10, expressed as $\varepsilon$/$m_{A^{\prime}}$ constraints.

Figure 10:  Same as Fig. 6 but for BP2a (left) and BP2b (right), corresponding
to the mixed inelastic Dark Matter (mixed-iDM) model. The two panels represent
two different realizations of the model, obtained by varying the $\alpha$
mixing angle.

This model main feature, which is also expected in BP3, is an enhanced
$A^{\prime}$ decay BR into $\psi_{2}\psi_{2}$. The branching ratio to
$\psi_{2}\psi_{1}$ is suppressed by a factor of the mixing angle $\alpha$, and
the one into two $\psi_{1}\psi_{1}$ is forbidden by the $C$ symmetry. Because
most dark photon events come accompanied by two unstable particles, the
additional energy deposition is missed even less often, relaxing the bounds
further.

The relic density of $\psi_{1}$ for this model depends strongly on the
efficiency of coannihilations and co-scattering processes. In the realization
shown in Fig. 10, we find that a simultaneous explanation of $\Delta a_{\mu}$
and dark matter relic density, along with the constraints discussed, can be
achieved in the region $0.9~{}\mathrm{GeV}\lesssim m_{A^{\prime}}\lesssim
1.2~{}\mathrm{GeV}$ for BP2a. In the case of BP2b, the coannihilation
processes are inefficient due to the smaller $\alpha$, so that $\psi_{1}$ is
overabundant in the region of parameter space that can explain $\Delta
a_{\mu}$.

In addition, we report an analysis on the $\Delta_{21}$ and $\alpha$
parameters, showing the constraints in a $\Delta_{21}/\alpha$ in Fig. 11,
fixing the mass of the dark photon to
$m_{A^{\prime}}=1,\,1.25,\,2~{}\mathrm{GeV}$ and
$\varepsilon=0.01,\,0.02,\,\varepsilon_{(g-2)}$.

Figure 11: Parameter space of the mixed-iDM in the plane of $\Delta_{21}$ and
the mixing angle $\alpha$. The parameter $m_{A^{\prime}}$ has been fixed to
$1.0$, $1.25$ and $2.0~{}\mathrm{GeV}$ (column-wise), while the kinetic mixing
$\varepsilon$ has been chosen among $0.01$, $0.02$, and value
$\varepsilon_{g-2}$ (row-wise), corresponding to the central value of the
$\Delta a_{\mu}$ explanation. BaBar and NA64 recast bounds are shown with the
same style used in Fig. 6.

For some combinations of the parameters, the NA64 constraints are not present
because they are too weak, given the efficient relaxation that this model can
provide. The dependence on the angle $\alpha$ is expressed by the branching
ratio: a larger value favors a larger branching ratio to the
$\psi_{2}\psi_{1}$ channel, with respect to $\psi_{2}\psi_{2}$; it has the
effect of decreasing the amount of visible energy, and ultimately the
possibility to detect a semi-visible event. On the contrary, a smaller
$\alpha$ affects the decay rate of $\psi_{2}$, suppressing its decay, and
recovering the original invisible bound. The behavior of the constraints is
similar to BP1 for what concert $\Delta_{21}$ dependence: a larger
$\Delta_{21}$ means a shorter $\psi_{2}$ lifetime, and the energy of its decay
is released inside the detector. On the other hand, a lower $\Delta_{21}$
resembles the case of fully invisible dark photon decays, as the lifetime
becomes larger than the size of the detector.

Figure 12:  The relic density (solid) in the parameter space of mixed-
inelastic DM model in the plane of kinetic mixing $\varepsilon$ versus dark
photon mass $m_{A^{\prime}}$ for fixed values of $\alpha=0.5$,
$r=m_{1}/m_{A^{\prime}}=1/3$, and $\Delta_{21}$. CMB limits are not applicable
as the dark matter self-annihilation $\psi_{1}\psi_{1}\to f^{+}f^{-}$ is
forbidden in the $C$ symmetric limit.

The relic density lines in Fig. 11 identify the overabundant region,
corresponding to high $\Delta_{21}$ and low $\alpha$. In that case,
coannihilations and coscattering processes become too inefficient. For small
$\Delta_{21}$, the choice of mixing angle does not have a strong impact, as
coscattering remains efficient for longer, and the self-annihilation of
$\psi_{2}$ sets the relic abundance of $\psi_{1}$. For large $\alpha$ the
dependence on $\Delta_{21}$ on the relic density is relaxed: the enhanced
coscattering ratio obtained with a larger $\alpha$ allows to afford a larger
$\Delta_{21}$ value before ending up in an underabundant scenario. The kink
present in the dark matter relic density for
$m_{A^{\prime}}=1.0~{}\mathrm{GeV}$ is due to the presence of a resonance
region, that can be observed also in Fig. 12. In that figure we show the trend
of the relic density line for different choice of the parameters of the model,
along with the $3\sigma$ region accounting for the $\Delta a_{\mu}$
explanation.

Inelastic Dirac Dark Matter (BP3a–d) — We show the constraints for the mixed-
iDM model in Fig. 13, expressed as $\varepsilon$/$m_{A^{\prime}}$ constraints.

Figure 13:  Same as Fig. 6 but for BP3a (top left) and BP3b (top right), BP3c
(bottom left), and BP3d (bottom right), corresponding to the inelastic Dirac
dark matter (i2DM) model. The four panels represent the same choice of
parameters, varying solely the mixing angle, $\beta$.

Similarly to the previous model, the constraints are relaxed, and a new region
of the parameter space opens up. This model is characterized by an enhanced
dark photon decay rate to the channel $\psi_{2}\psi_{2}$ as it happens in the
case of the mixed-iDM model. The rate into the channels $\psi_{2}\psi_{1}$ and
$\psi_{1}\psi_{1}$ is suppressed by respectively a factor $\beta$ and
$\beta^{2}$. Differently from the mixed-iDM case, this model allows for dark
photon decays to $\psi_{1}\psi_{1}$ channel. A larger branching ratio to the
heaviest HNF increases the possibility to detect a semi-visible event in the
detector, given the larger abundance of those particles releasing $e^{+}e^{-}$
pairs after their decay.

The relic density of $\psi_{1}$ for this model depends on the efficiency of
coannihilations and coscattering processes. The main difference with respect
to the previous model is that it is not possible to evade the CMB bounds,
because of the possibility for the dark matter candidate $\psi_{1}$ to
annihilate through the vertex $\psi_{1}\psi_{1}$. Even though this vertex is
suppressed, it can have a sizable contribution to late time annihilations,
injecting additional energy into the CMB. In the realizations shown in Fig.
13, we find that, despite the relaxation of the main constraints, the CMB
bounds are unavoidable and can exclude large parts of the parameter space,
with the only exception being the choice of a small $\beta$ parameter, as
represented by benchmarks BP3c and BP3d. This choice corresponds to suppress
further the channels $\psi_{2}\psi_{1}$ and $\psi_{1}\psi_{1}$. However, it
has an impact on the relic abundance of $\psi_{1}$, because it suppresses the
contributions of coscattering and annihilations, resulting in an overabundant
scenario, which can be set under control only assuming secluded annihilations
within the DS.

Additionally, the relic density depends also on $\Delta_{21}$, because it
affects the Boltzmann suppression of the coannihilation $\psi_{2}$. It is
interesting to understand the interplay of $\Delta_{21}$ and $\beta$ in the
different constraints. In Fig. 14, the bounds are shown in a
$\Delta_{21}/\beta$, fixing the mass of the dark photon to
$m_{A^{\prime}}=1,\,1.25,\,2~{}\mathrm{GeV}$ and
$\varepsilon=0.01,\,0.02,\,\varepsilon_{(g-2)}$.

Figure 14:  Same as Fig. 6 but for BP3, corresponding to the inelastic Dirac
dark matter (i2DM) model. The parameter space is shown as a function of
$\Delta_{21}$ and the mixing angle $\beta$. The region dotted in black is
excluded by CMB limits, providing a full exclusion of these slices of
parameter space.

We can draw similar conclusions as the ones discussed for the Mixed-iDM model
in Fig. 11, with the difference that the $x$-axis now represents the parameter
$\beta$. In addition, the CMB constraints are shown: their exclusion region
corresponds to large $\beta$ values, because of the enhancement of
$\psi_{1}\psi_{1}$ annihilation rate. Regarding the relic density, we can draw
a similar conclusion as for benchmarks BP2a and b based on Fig. 11 and Fig.
15. The overabundant region corresponds to large $\Delta_{21}$ and low
$\beta$, due to both inefficient coannihilations and suppressed $\psi_{1}$
self-annihilations. However, the dependence on $\beta$ is stronger due to the
presence of $\psi_{1}$ self-annihilations, which can dominate over
coannihilations in depleting the DM density. For this reason, the relic
density is underabundant for large $\beta$ independently on the value of
$\Delta_{21}$: the coannihilator does not play a crucial role anymore in the
determination of the dark matter relic density. Nevertheless, the small
regions in which dark matter density is underabundant and compatible with
other constraints are excluded by the CMB bounds.

Figure 15:  The relic density (solid) and CMB limits (dashed) in the parameter
space of inelastic Dirac Dark Matter (i2DM) in the plane of kinetic mixing
$\varepsilon$ versus dark photon mass $m_{A^{\prime}}$ for fixed values of
$\alpha_{D}=0.5$, $r=m_{1}/m_{A^{\prime}}=1/3$, and $\Delta_{21}$. Each curve
corresponds to a different value of the mixing angle $\beta$. The CMB limits
exclude large values of $\varepsilon$. The $\Delta a_{\mu}$ $3\sigma$
preferred region is shown as a green band.

Three Heavy Neutral Fermions (BP4a–c, BP5) — We show the constraints for the
Three Heavy Neutral Fermions models in Fig. 16, expressed as
$\varepsilon$/$m_{A^{\prime}}$ constraints.

Figure 16: Same as Fig. 6 but for BP4a (top left), BP4b (top right), BP4c
(bottom left), BP5 (bottom right), corresponding to the models with 3 heavy
neutral fermions. The dark photon branching ratios are dominated by
$A^{\prime}\to N_{2}N_{3}$ decays.

Similarly to the previous cases, the constraints coming from both BaBar and
NA64 are relaxed, and a new region of the parameter space opens up. All
benchmarks are characterized by having three HNFs, and by a sizable $V_{32}$
coupling, which enhances the dark photon decay rate to $N_{3}N_{2}$ final
states. Furthermore, the produced HNFs can promptly decay, releasing
$e^{+}e^{-}$ pairs in the detectors. The presence of a new fermion, and the
enhanced annihilation rate to the heaviest HNFs make it possible to have a
larger number of $e^{+}e^{-}$ pairs, and, consequently, a larger probability
of detecting a semi-visible event. The main difference between the models is
that BP4a–c are characterized by only off-diagonal couplings, with the only
possible decay chains being:

$\displaystyle A^{\prime}$ $\displaystyle\to(N_{2}\to
N_{1}e^{+}e^{-})\,(N_{3}\to(N_{2}\to N_{1}e^{+}e^{-})e^{+}e^{-}),$
$\displaystyle A^{\prime}$ $\displaystyle\to(N_{2}\to N_{1}e^{+}e^{-})N_{1}.$

Differently, BP5 allows for any possible coupling among the HNFs.

The downward shift of the BaBar bound happening between BP4a and BP4b
benchmarks is due to the different values assumed by the parameter
$r=m_{1}/m_{A^{\prime}}$. This parameter affects both the HNF lifetime and the
$A^{\prime}$ branching ratios. The HNF lifetime depends on $r$ according to
$c\tau\sim m_{A^{\prime}}^{-1}r^{-5}$, so a larger value would imply a shorter
lifetime, translating into a more relaxed bound, because the potential larger
fraction of events releasing energy inside the instrumented regions of the
detector. However, as it can be observed, the bound becomes stronger from BP4a
to BP4b, even with a larger $r$. The new value modifies the branching ratio of
the dark photon decay and forbids the channel $A^{\prime}\to N_{3}N_{2}$,
because the kinematics requires:

$m_{A^{\prime}}>m_{2}+m_{3}\Rightarrow
r<\frac{1}{(\Delta_{32}+2)(\Delta_{21}+1)}\approx 0.116,$ (45)

which is satisfied for BP4a, but not for BP4b. Being the decay to the two
heavy HNFs forbidden means that the potential production of $e^{+}e^{-}$ pairs
is suppressed, because $A^{\prime}$ can decay only to $N_{2}N_{1}$, and only
$N_{2}$ can decay further.

The constraints show that the NA64 (S2) projected bound has the capability of
excluding new regions of the parameter space, demonstrating the capability of
the experiment to be sensitive to promptly decaying HNFs.

##  VI Prospects and constraints on HNL and DM interpretations

###  VI.1 HNFs as dark matter candidates

The U$(1)_{D}$ symmetry can be responsible for the stability of the lightest
HNF, and, therefore, provide a dark matter candidate. Light dark matter models
with self-annihilating Dirac fermions are excluded by CMB data due to the
s-wave, velocity-unsuppressed annihilation. Majorana fermions or scalar
particles have p-wave, velocity-suppressed, and self-annihilations; however,
in this case, the required values of self-interactions render the dark photon
fully invisible, and, therefore, excluded as an explanation of $\Delta
a_{\mu}$. Self-annihilation near the $A^{\prime}$ resonance,
$r=m_{1}/m_{A^{\prime}}\lesssim 1/2$, can significantly enhance cross
sections, but such mass spectrum would leave no room for semi-visible,
promptly-decaying fermions.

Coannihilations, $\psi_{1}\psi_{2,3}\to(A^{\prime})^{*}\to f^{+}f^{-}$, are
therefore the most natural possibility to achieve freeze-out. The
coannihilation cross section of opposite $C$ states, $\psi_{i}$ and
$\psi_{j}$, is given by

$\sigma v(\psi_{i}\psi_{j}\to
f^{+}f^{-})=8\pi\alpha_{D}\alpha\varepsilon^{2}|V_{ij}|^{2}\frac{(2m_{ij}^{2}+m_{f}^{2})}{(4m_{ij}^{2}-m_{A^{\prime}}^{2})^{2}+m_{A^{\prime}}^{2}\Gamma_{A^{\prime}}^{2}}\sqrt{1-\frac{m_{f}^{2}}{m_{ij}^{2}}}+\mathcal{O}(v^{2}),$
(46)

where $m_{ij}=(m_{i}+m_{j})/2$. Just like the self-annihilation of Dirac
fermions, the cross section is velocity unsuppressed. Nevertheless, the
annihilation is exponentially suppressed at late times due to the mass
splitting between the two states and the subsequent decays of the co-
annihilator.

To calculate the DM relic density of $\psi_{1}$ in our benchmarks, we assume
that all dark sector fermions are in chemical equilibrium at the time of
freeze-out, and employ the formalism of Ref. [147]. We find good agreement
with the literature on iDM [30, 21] and i2DM [45]. We find a $50\%$
disagreement with the relic curves of Ref. [50] for $m_{1}\gtrsim 100$ MeV.

Direct detection — Direct detection of a dark matter particle of mass
$m_{\chi}\sim\mathcal{O}(100)$ MeV, with large kinetic mixing, would provide
strong evidence for the DM nature of the HNFs. For the parameter space we
consider in iDM and mixed-iDM models, low-energy direct detection can only
probe the loop-induced elastic scattering of DM. The tree-level upscattering
rates are exponentially suppressed as only the largest DM velocities can
overcome the kinematical threshold of the large mass splittings. Direct
detection prospects are instead dominated by the loop-induced, elastic DM-
quark coupling.

In the case of kinetic mixing, the scalar-current dominates,
$c_{5}^{q}(\overline{\chi}\chi)(\overline{q}q)$ [148]. In terms of coupling to
nucleons,

$\sigma_{\chi N}^{1{\rm-loop}}=\frac{(c_{1}^{N}\mu_{N})^{2}}{\pi}$ (47)

where $\mu_{N}=m_{\chi}m_{N}/(m_{\chi}+m_{N})$ is the reduced mass of the DM
and the nucleon $N$, and $c_{1}^{N}$ is the loop-induced, DM-nucleon coupling.
The matching to nucleon currents gives [149]

$c_{1}^{N}=4\frac{\alpha_{D}\alpha_{\rm
QED}\varepsilon^{2}}{m_{A^{\prime}}^{4}}m_{\chi}F_{3}(r^{2})\sum_{q}F_{S}^{q/N}Q_{q}^{2}\sum_{i}|V_{1i}|^{2},$
(48)

where $Q_{q}$ is the quark electric charge, $r=m_{\chi}/m_{A^{\prime}}$,
$F_{3}(x)$ is the loop function from Ref. [148], and $F_{S}^{q/N}(Q^{2})$ the
nucleon scalar form factors. Approximating the form factors to their $Q^{2}=0$
value and $F_{S}^{q/p}(0)\simeq F_{S}^{q/n}(0)\simeq(15,35,40)$ MeV for
$q=(u,d,s)$, we estimate the elastic DM-nucleon cross sections. At a typical
point of parameter space,

$\sigma_{\chi N}^{1-{\rm loop}}\simeq 1.4\times 10^{-3}\text{
pb}\times\left(\frac{\alpha_{D}}{0.3}\right)^{2}\left(\frac{\varepsilon}{10^{-2}}\right)^{4}$
(49)

where we assumed $m_{A^{\prime}}=3m_{\chi}=1$ GeV and
$\sum_{i}|V_{1i}|^{2}=1$. This value is not currently probed by any low-energy
direct detection experiments. In this mass region, CREST-III (2019) [150]
provides the leading limits on elastic DM-nucleus scattering using nuclear
recoils. Those limits are over two orders of magnitude above our estimate. The
next-generation SuperCDMS detectors at SNOLAB may be able to probe part of the
parameter space of interest using nuclear recoils [151, 152]. A more promising
avenue in sub-GeV DM direct detection, however, is the use of nuclear-
inelastic processes. In this case, DM can impart all of its kinetic energy
into excitating a target nucleus, which subsequently de-excited emitting an
electron through the Migdal effect, or a photon. The electron recoil, in this
case, can significantly improve the prospects for sub-GeV DM direct detection
[153, 154, 155, 156]. This method has been used by the LUX [157], SENSEI
[158], XENON1T [159], and DarkSide-50 [160] collaborations to set limits in
our mass region of interest. The best constraints come from DarkSide, where
$\sigma_{\chi N}\lesssim 0.1$ pb at $m_{\chi}=300$ MeV at 90% C.L.

Unfortunately, the loop-induced scattering on electrons is very suppressed in
iDM and mixed-iDM models due to the small electron mass.

In the i2DM model, direct detection is sensitive to the mixing-angle-
suppressed elastic scattering of $\chi\equiv\psi_{1}$ on electrons. For a
heavy dark photon, the total cross section is

$\displaystyle\overline{\sigma}_{e}$
$\displaystyle=\frac{16\pi\sin^{4}{\beta}\alpha_{D}\alpha_{\rm
QED}\varepsilon^{2}\mu_{e}^{2}}{m_{A^{\prime}}^{4}}$ (50) $\displaystyle\simeq
4\times
10^{-7}\text{\,pb}\times\left(\frac{\sin{\beta}}{0.14}\right)^{\\!\\!4}\\!\\!\left(\frac{\alpha_{D}}{0.3}\right)\\!\left(\frac{\varepsilon^{2}}{10^{-2}}\right)\\!\left(\frac{1\text{
GeV}}{m_{A^{\prime}}}\right)^{\\!\\!4}\\!\\!.$

The leading limits in this parameter space are from XENON1T [161, 162],
PandaX-II [163], and SENSEI [158]. At $m_{\chi}\sim 100$ MeV, XENON1T
constrains $\overline{\sigma}_{e}\lesssim 10^{-4}$ pb, still orders of
magnitude above our estimates for $\beta=8^{\circ}$. The prospects are more
interesting for scattering on protons, where the bounds discussed just above
apply as well. In this case, DarkSide-50 already probes the largest values of
kinetic mixing and $\beta$ for $m_{A^{\prime}}\gtrsim 1.5$ GeV. However, these
are already excluded by BaBar and CMB constraints.

Another possibility for direct detection is to search for a boosted DM
population [164, 165, 166, 167]. Cosmic rays can interact with the DM
background, upscatter $\chi\to\psi_{2,3,\dots}$, which subsequently decay to
fast DM particles. This cosmic-ray-boosted DM population can then be searched
for in direct detection and neutrino experiments. Refs. [168, 169] derive
limits on similar models using XENON1T data, from where we can conclude that
current limits are still too weak to constrain our parameter space, in all
models of interest. Large neutrino detectors can further enhance the
sensitivity thanks to their large mass and excellent detector performance
[165]. A more detailed study is needed to assess the flux of boosted DM
particles in our models and their testability via this strategy.

Cosmic Microwave Background — Precision measurements of the Cosmic Microwave
Background (CMB) also provide significant limits on the models we consider
when the HNFs are dark matter. If the dark matter fermions significantly
annihilate or decay to charged particles at the time of recombination, they
can inject additional energy into the SM plasma, re-ionize Hydrogen, and delay
the formation of the CMB [170, 171, 172, 173, 174]. The latest constraints
from Planck [175] rule-out light and thermal dark matter candidates with
s-wave annihilations for $m_{\chi}\simeq$. This constraint is much weaker and,
therefore, not significantly constraining for models with co-annihilating dark
matter candidates, like iDM and mixed-iDM, and in models with p-wave
annihilations, like the case of Majorana dark matter fermions.

Out of the models and mass splittings we consider, only the i2DM model is
subject to such constraints. This is because the light Dirac dark matter
fermion can undergo velocity-unsuppressed self-annihilations, even if
suppressed by the fourth power of a small mixing angle $\beta$. The self-
annihilation cross section to charged leptons in the i2DM model is given by
Eq. 46, with $m_{ij}\to m_{1}$ and $V_{ij}\to\beta^{2}$. The curves where the
correct relic density of DM can be achieved are compared with the limits from
CMB in Figure 10. For the typical lifetimes and mass splittings considered in
this work, the late-time annihilations involving $\psi_{2},\psi_{3},\dots$ can
be safely neglected.

###  VI.2 HNFs as heavy neutral leptons

Having discussed the theoretical aspects of a HNL interpretation in Section
II.4, we now turn to the phenomenological consequences. Searches of HNLs
require mixing with active neutrinos, which emerges from the Yukawa coupling
between the sterile neutrinos and the leptonic doublet after symmetry
breaking. While a successful $\Delta a_{\mu}$ explanation does not lead to any
constraint on this mixing, the latter will be constrained from below by BBN,
such that $\tau_{N_{4}}\lesssim 0.1$ s, and from above by laboratory searches.
As highlighted in Refs. [176, 37], the phenomenology of the models under
consideration can be very different from the minimal case.

For most of the parameter space of interest in this paper, the heavier HNLs
will decay very fast into lighter HNLs and dark photons, into 3 lighter HNLs,
and into HNLs and dilepton pairs, depending on which channels are
kinematically allowed. We focus on a hierarchy of HNF and dark photon masses
such that the latter decay dominates, allowing to evade BaBar bounds as
discussed in the previous section. The lightest HNL decays primarily into a
dilepton pair and missing energy.

Thanks to the presence of a light dark photon both HNL scattering and decays
are enhanced compared to the standard case in which HNL interact with the SM
via mixing with the neutrinos.

HNLs are tested experimentally mainly via peak searches and via visible
decays. In the former case, the emission of a HNL is a pion or kaon decay
leads to a small peak in the charged lepton spectrum at a lower energy. These
bounds are very robust as they rely uniquely on the kinematics of the meson
decay and pose some of the strongest constraints in the sub-GeV HNL mass
region [177]. Similarly, for HNLs coupled exclusively to the tau flavor, peak
searches in $\tau$ and $D$ decays, such the recent BaBar analysis [178],
provide strong limits. In the models we propose, even peak searches can be
affected as, for very fast decays, these events would be vetoed by the
requirement of no additional charged particles. A weakening of the bounds can
be expected in certain ranges of parameter space. A more detailed discussed is
provided in Ref. [37].

GeV-scale HNLs can be produced via mixing in meson decays and in neutrino
scatterings, typically in beam dump and neutrino experiments. In the first
setup, high energy protons impinge on a target producing copious amounts of
pions, kaons and, for sufficiently high energies, heavier mesons, which
subsequently decay producing HNLs. These travel some distance before decaying
into missing energy and visible particles that can be revealed in dedicated or
multi-purpose neutrino detectors. Due to the kinetic mixing and light dark
photon mass,$\Gamma^{A^{\prime}}_{N_{4}\to\nu e^{+}e^{-}}\gg\Gamma^{\rm
W,\,Z}_{N_{4}\to{SM}}$, where $\Gamma^{\textup{Md}}$ is the decay rate
mediated by the particle Md. There are two cases: in the long-lived regime,
$c\tau_{4}^{\rm LAB}>L$ with $L$ the baseline of the experiment, the event
rate of $N_{4}$ decays in DIF searches can be enhanced as it scales as
$\Gamma^{A^{\prime}}_{N_{4}\to\nu e^{+}e^{-}}/L$ and the bounds get stronger.
Alternatively, if $N_{4}$ is too short-lived, $c\tau_{4}^{\rm LAB}\ll d$, with
$d$ being the distance between the source and the detector, the limits do not
apply at all as the HNLs do not even reach the detector.

The strongest limits on $U_{e4}$ and $U_{\mu 4}$ are set by T2K [179] and
MicroBooNE [180] (see also PS191 limits [181] and the discussion in Refs.
[182, 183]), while $U_{\tau 4}$ is only constrained by higher-energy
experiments like CHARM [184] (see also [185, 186, 187]), NOMAD [188], and LEP
[189, 190, 191]. As discussed, these bounds need to be revisited in the light
of the considerations above and depend critically on the choice of parameters.
A more in depth discussion is available in Ref. [37].

In the second type of setups, HNLs are produced by a neutrino beam via
upscatterings in the detector itself and subsequently decay leading to visible
signatures. In this case, the upscattering cross-section can be very
significantly higher than the standard HNL one. This, combined with much
shorter decay lengths, can lead to striking signatures in neutrino experiments
with short baselines, such as at the SBN programme at Fermilab, and at near
detectors of long baseline accelerator neutrino experiments, including T2K,
NoVA, DUNE.

A particularly interesting signature is the electron-positron pairs from HNLs
decays produced by neutrino upscattering in the MiniBooNE experiment. It has
been shown that for suitable values of the parameters, this signature can
explain the anomalous excess events at MiniBooNE [37], as well as the
$(g-2)_{\mu}$ anomaly. The particle $\psi_{2}$ could be efficiently produced
and decay into a dilepton pair which can mimic an electron neutrino
scattering, as either the two leptons are very collimated or one of them is
not reconstructed being too soft [192]. This explanation critically relies on
the large values of kinetic mixing that are allowed in our models.
Specifically, in order to explain the MiniBooNE anomalous excess, it is
necessary for the HNLs to decay within the detector, that has a typical size
of few meters. Decay lengths $\tau_{HNL}\ll 1$ m cannot be obtained in the
standard HNL scenario and require light dark photons and large kinetic mixing.

###  VI.3 Prospects for detection

In this subsection, we discuss future prospects for the detection of semi-
visible dark photons and HNFs.

Low-energy $e^{+}e^{-}$ colliders — Direct tests of the parameter space
present in our work can be achieved with a dedicated semi-visible search at
BaBar, KLOE, BES-III, and Belle-II $e^{+}e^{-}$ colliders. These searches can
be divided into two categories. On-shell production of dark photon through
initial state radiation, $e^{+}e^{-}\to A^{\prime}\gamma$, or the production
of HNFs through off-shell dark photons,
$e^{+}e^{-}\to(A^{\prime})^{*}\to\psi_{i}\psi_{j}$. In particular, low-energy
machines like KLOE and BES-III ran with center-of-mass energies close to the
dark photon mass, significantly enhancing their prospects for direct
production of HNFs.

_Initial State Radiation —_ One advantage of keeping the ISR topology is the
kinematic imbalance in the photon-dark-photon center-of-mass system. Since
multiple invisible particles are emitted in the decay cascade, it is not
possible to reconstruct the dark photon mass with visible energy. However, by
isolating the photon, a resonance on $M_{X}^{2}=s-2E_{\gamma}\sqrt{s}$ would
still be visible. A detailed sensitivity study of the Belle-II reach to iDM
through this channel was carried out in Refs. [30, 31]. We expect the
sensitivity to be even better in models with two or more HNF decays and leave
a detailed study for mixed-iDM and HNL models to future literature.

_S-channel Production —_ Unlike the ISR channel, the s-channel production
cross section is proportional to the dark coupling $\alpha_{D}$, and so can be
large for models where $\alpha_{D}\gg\alpha$ [37, 193]. In terms of the cosine
of the angle of $\psi_{i}$ with respect to the beam in the COM, $c_{\theta}$,
the differential cross section is given by,

$\frac{\differential\sigma}{\differential
c_{\theta}}=|V_{ij}|^{2}\frac{\pi\alpha_{D}\alpha_{\rm
QED}\varepsilon^{2}}{(E_{\rm
CM}^{2}-m_{A^{\prime}}^{2})^{2}+\Gamma_{A^{\prime}}^{2}m_{A^{\prime}}^{2}}\frac{E_{\rm
CM}^{2}}{2}\left[1+\left(c_{\theta}-(1-c_{\theta})\frac{m_{i}^{2}\Delta_{ij}(2+\Delta_{ij})}{E_{\rm
CM}^{2}}\right)^{2}\right],$ (51)

where $\Delta_{ij}=(m_{i}-m_{j})/m_{j}$, $E_{\rm CM}=\sqrt{s}$ is the center-
of-mass collision energy, and $\Gamma_{A^{\prime}}$ is the total decay width
of the dark photon. Unlike ISR events, where the displaced HNF vertices would
depend on the direction of travel of $A^{\prime}$, s-channel production would
provide a source of back-to-back displaced vertices. The dileptons would not
point allow pointing back to the collision point due to the missing energy. In
addition, for secondary decays, like $\psi_{3}\to\psi_{2}\to\psi_{1}$, a
third, lower-energy dilepton pair could be visible, keeping the two primary
decay vertices and the collision point on the same line. In Ref. [193], the
authors have explored these events in iDM models, finding that Belle-II can
cover open regions of parameter space. The sensitivity reaches values of
$\varepsilon<10^{-3}$ for dark photon precisely in the region of interest for
$\Delta a_{\mu}$, $m_{A^{\prime}}\gtrsim 500$ MeV. The case of mixed-iDM and
i2DM have not been studied, but the additional semi-visible vertices can
provide additional discrimination from backgrounds, and extend its reach into
parameter space. A detailed study of these events is left to future
literature.

In the presence of a signal at Belle-II, both of the channels above would shed
light on the mass splittings and masses of the HNFs. Firstly, the ISR channel
could reconstruct the dark photon mass via
$m_{A^{\prime}}\overset{!}{=}M_{X}^{2}=\equiv s-2E_{\gamma}\sqrt{s}$. Then, in
both ISR and s-channel events, the dilepton invariant mass of