Semi-Visible Dark Photon Phenomenology at the GeV Scale

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.

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 ′ → ψi(ψj → ψ k ℓ + ℓ − ).We find that in some cases the BaBar and NA64 experiments no longer exclude large kinetic mixing, ε ∼ 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 ′ decays and directly probe the newly-open parameter space.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 parti-cles [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 ′ , 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], ε 2cW X µν B µν , is expected to be sizeable, providing a clear target for detection.The one-loop expectation for kinetic mixing is where M i and Q (Y, X) i are the masses and charges of the heavy new fermions that run in the loop, µ 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 ′ 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 ε [14].It is, however, possible to avoid laboratory constraints and remain compatible with the naive one-loop estimation of ε.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 [13,[22][23][24].
Another motivation for this work is the measurement of the anomalous magnetic moment of the muon, a µ = (g − 2) µ /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 µ at the one loop-level with a positive sign, providing an elegant and simple solution to the discrepancy between experiment and theoretical predictions, ∆a µ [27,28].This solution is only possible for light mediators, m A ′ ≲ 3 GeV, and requires large values of kinetic mixing, ε ∼ 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 ′ ∼ 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 ′ , based on a minimal inelastic DM (iDM) model [29], is strongly constrained by collider, fixed target, and indirect searches and can only explain ∆a µ in a very narrow region of parameter space.Our collection of semi-visible A ′ models constitutes a viable explanation of ∆a µ 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 ′ 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 monophotonlike events at Belle-II, such as those studied by the BaBar e + e − collider.At fixed-target experiments, we point out that invisible-A ′ 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 where X µν is the field strength tensor of the dark photon, and J µ D is the dark current containing new fermionic or scalar degrees of freedom.The origin of the dark photon mass is thus far unspecified 1 .After the diagonalization of the gauge kinetic terms, the dark photon mass eigenstate A ′ with mass m A ′ ≃ m X couples to both the SM electromagnetic (EM) current and the SM weak neutral current (NC), where t W ≡ tan θ W , with θ 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 ′ 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 Φ, charged under the U (1) D .The dark current is J µ D = Φ * i ↔ ∂ µ Φ, and the hermiticity of the Lagrangian ensures that the dark photon interactions must always take place between the real components of Φ, φ 1 and φ 2 .A soft U (1) D -breaking term, µΦ2 can then split the masses of the real scalars, and render the dark photon semi-visible due to the decay cascade A ′ → φ 1 (φ 2 → φ 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 V ij ψ i γ µ ψ 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 While we do not pursue this possibility, we note that non-renormalizable interactions between ψ 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 χ 1 , with the heavier particles in the spectrum χ j=2,...,n .We retain here a more general notation, ψ i=1,...,n , which encompasses both options.The heavier states ψ j=2,...,n , are expected to decay in cascades down to the lightest HNF, emitting two charged particles at each step 2 .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, ψ 1 , that is stable and constitutes the DM candidate.In this work, we generalize this idea to richer dark sectors.
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 ψ 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.

A. Two heavy neutral fermions (HNFs)
The simplest model we consider is that of two Majorana HNFs.In the interaction basis, we take χ L and χ R with charges Q L and Q R .The most general Lagrangian for their mass reads where the Majorana masses µ L and µ R break the U (1) D softly -they can be generated by the vacuum expectation value of a dark Higgs Φ 2 with Q Φ2 = 2Q L,R .The diagonalization of the symmetric and complex mass matrix M can be achieved with the usual Takagi diagonalization, diag(m 1 , m 2 ) = U T M U , where U = R(θ)diag(e iφ , 1) with R(θ) the matrix rotation by an angle, We define ∆µ = (µ R − µ L )/2 and µ = (µ L + µ R )/2.CP conservation is ensured when the Majorana phase is φ = 0 or π/2.In terms of the Majorana mass eigenstates, the dark current is given by The C symmetry -For ∆µ → 0, the mixing angle θ 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 where ψ c L = Cψ L T and we choose the phase factor η c = +1, for simplicity.In the case of a Dirac fermion, the C operation is achieved by charge conjugating the field, χ The C symmetry is then respected when the Lagrangian is invariant under the exchange ξ ↔ ζ (i.e., χ L ↔ χ c R ).The left-handed fermions constitute the C eigenbasis, where φ is the same Majorana phase as before.Given that C(A ′ µ ) = −1, the intrinsic C-parity of the fermions can be fixed as C(χ ± ) = ±1 without loss of generality.The Lagrangian in Eq. ( 5) in this basis reads where we took φ → π/2 to ensure that the mass terms are positive for m D > µ.This signals that the two fermions have opposite CP parities.This basis is identified with the physical basis when ∆µ → 0. As expected, ∆µ and Q A are the only parameters that break C in this model.In the C-symmetric limit, χ ± behaves like the components of a pseudo-Dirac particle with a mass gap 2µ.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 ′ µ 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 SM , and one in the DS, C DS , then the conservation of C DS , but not 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 φ = π/2 and Q V = 1, the dark current is simply This phenomenological model is defined by the following five parameters , and ε.
(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 ∆ 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 ψ 1 can only proceed through scalars in this model.If the scalar mixes with the Higgs, λ|Φ| 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 subdominant to the exponentially suppressed A ′ -mediated coannihilations.In addition, we assume the dark higgs is heavier than the HNFs, as otherwise secluded annihilation would dominate 3 .
As we will show in Section V, the explanation of ∆a µ 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, ψ 2 → ψ 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 ψ 2 can be produced in pairs, A ′ → ψ 2 ψ 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.
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 ∆a µ anomaly and a dark photon explanation.We will prove this in Section V with a detailed analysis.

B. Three HNFs
We start by extending the two HNF model by a single fully sterile Weyl fermion.We keep the two dark fermions, χ R and χ L , with the same charges Q L = Q R = 1, and introduce a new singlet fermion, η L .This content is a simplified version of the three-portal model of Ref. [37].
The Lagrangian is given by The mixing terms break the U (1) D and can be generated by the vacuum expectation value of a scalar particle Φ 1 , which carries charge where Y L,R are the Yukawa couplings.Another dark Higgs Φ 2 with Q Φ2 = 2 could generate the Majorana masses of the two dark fermions after symmetry breaking.
Because η L is completely neutral, it can couple to the SM lepton doublets via the Yukawa coupling L Hη c L .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 D. 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, χ L , χ R , η 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(η L ) = 0, which would forbid the neutrino Yukawa coupling [54].
We use the left-handed dark fermion basis χ + and χ − 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 > µ.In that case, the DS fermion mass matrix is where Imposing the C symmetry in the χ sector4 , we recover the limit where ∆µ = ∆Λ = 0. We find an analogous situation to the two HNF cases, with the difference that χ + can now mix with a sterile state.Indeed, since C(η L ) = +1, C-conservation implies that only the C-even fermion can mix with η L .As we will see, the spectrum can consist of one Dirac and one Majorana particle or three Majorana states.
The C-odd state χ − ≡ ψ 2 decouples, and η and χ + mix.A single rotation in the C-even sector leads to the mass basis, with tan 2α = 2Λ/M and M = M X + µ − µ ′ L .If tan 2α ≪ 1, a seesaw mechanism is in place, and ψ 2 and ψ 3 form a pseudo-Dirac pair.The other possibility, tan 2α ≫ 1, does not preserve the pseudo-Dirac limit.The splittings in the model are then given by m 3 − m 1 ∼ M + 2Λ 2 /M and m 3 − m 2 ∼ Λ 2 /M + 2µ.
In the C symmetric case, the dark current is also fully off-diagonal and is given by In the limit of small α, the dark photon interacts more strongly with the pseudo-Dirac pair.When the ψ 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 mostlydark, 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 ′ 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 ψ 1 and a (pseudo-)Dirac state Ψ 2 , the dark current is and so this model is fully specified by Eq. ( 12) and α.
To make use of the model above, one must guarantee the coherence of the Majorana states ψ 2 and ψ 3 in the pseudo-Dirac state Ψ 2 , so that we are justified in treating them as a single Dirac particle in the phenomenological work.Note that ∆ 32 will only play a minor role in the dark matter hypothesis since the relevant splitting for coannihilations is the one between ψ 1 , the dark matter candidate, and Ψ 2 , its interaction partner.We can express M X and Λ in terms of tan 2α, which controls the decay rate Ψ 2 → ψ 1 , and the splitting ∆ 21 ≡ (m 2 − m 1 )/m 1 , which has an important impact on the coannihilation rate for dark matter.For µ < M X , we find We notice that ∆ 32 = 1 4 ∆21 1+∆21 tan 2 2α and is small as far as the condition tan 2 2α ≪ 1 holds.For µ ′ L = µ = 0, the splitting of the pseudo-Dirac pair is ∆ 32 ∝ α 2 , so that in the limit of α → 0, we recover an exact Dirac state.In summary, provided the mixing angle is small, the ∆ 32 splitting is negligible.The decay rate for three-body decays like ψ 3 → ψ 2 + . . .are suppressed by ∆ 5  32 ∝ α 10 , and so, can be safely neglected for the mixing angles considered here.
Three Majorana fermions -Relaxing the condition on α ≪ 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 ψ 3 and ψ 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α.For 0 < µ ≃ M X , we have (M X − µ)/M X tan 2 2α ≃ 2∆ 21 /(1 + ∆ 21 ) and ∆ 32 ≃ ∆ 21 /(1 + ∆ 21 ), implying a mildly hierarchical HFN spectrum.Moving away from the C symmetric limit, a sizeable ∆Λ 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 ψ 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 ψ 1 is stable or long-lived, as is typical in these models.

C. 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, η, and one charged under the dark gauge symmetry, χ.Our Lagrangian in this case reads where again we have omitted potential Yukawa couplings between SM neutrinos and the sterile fermions, η L and η R (cf.Section II B).In the C symmetric limit, one can show that Λ L = Λ ′ R and Λ L = Λ ′ R .In this limit, for an appropriate choice of Majorana phases, the mass matrix in the C eigenbasis is, where The C-even and C-odd sectors decouple.
One more limit of interest is considering the U (1) D to be exclusively broken by one unit, such that µ = µ ′ = 0 and ∆ + = −∆ − ≡ ∆.In that case, the Dirac pairs are split by which is small for small mixing angles and vanishes when β + = β − .This last regime is the limit where the two pairs compose exact Dirac fermions, achieved in two cases: In terms of the mass eigenstates, the dark current takes the simple form, 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 ψ 4,2 → ψ 3,1 + . . .are suppressed with respect to the dominant ψ 4,3 → ψ 1,2 + . . .by factors of ( In fact, for the typical mixing angles we consider, the particle ψ 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 ∆(β 2 + − β 2 − ), where β + and β − 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 β + = β − .
Inelastic Dirac dark matter (i2DM) -In the exact Dirac limit, a simple phenomenological model arises with only two particles in the spectrum.A light, mostlyneutral Dirac fermion Ψ 1 , constituting a dark matter candidate, and a heavier, mostly-dark Dirac fermion Ψ 2 .In terms of these degrees of freedom, the dark current is given by As expected, Ψ 2 is coupled more strongly to the dark photon, and, therefore, will aid in the depletion of Ψ 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 Ψ 1 Ψ 1 → e + e − .The relic density is typically set by the self-annihilation of the heavy partner, Ψ 2 Ψ 2 → f + f − , and coscattering, . 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 selfinteractions of the co-annihilator, Ψ 2 .The phenomenology is fully determined by the parameters in Eq. ( 12), in addition to β.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 : 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.

D. 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 η are free to mix with the SM neutrinos, in the absence of any additional symmetries.For the model in Section II B, one can add the following Yukawa interaction where L α is the SU (2) leptonic doublet of the SM, and H = iσ 2 H * the conjugate of the Higgs doublet.Similar terms involving η R could be added to the four HNF models of Section II C.After EW and U (1) D symmetry breaking, the HNFs mix amongst themselves and we can justifiably call them HNLs, N i ≡ ψ i .With the addition of a dark scalar Φ with a dark charge Q Φ = 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 ′ mass considered in this paper, the lightest HNL, N 4 , will decay via N 4 → νℓ + ℓ − or N 4 → νπ + π − .Due to the suppression of the small mixing with the neutrinos, N 4 is usually long-lived (cτ N4 ≫ 1 m) and constitutes missing energy at e + e − colliders and fixedtarget experiments.The experimental consequences of the mixing with neutrinos is discussed in Section VI B. 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 η L so that the Yukawa coupling is lepton number conserving implies that the Majorana mass term µ ′ L breaks it by two units.In the dark sector, the charge assignment of χ L and χ R is arbitrary and, depending on the specific choice, lepton number will be broken by Λ L,R , M X and µ L,R , or by Λ L and µ L,R , for , respectively.Both Λ L,R and µ 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 Φ with dark charge Q Φ = 1, either the resulting Λ L , or Λ 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 Light neutrino masses need to depend on all the U (1) Lbreaking parameters.As an interesting example, let us consider the case in which the χ L,R do not carry lepton number and only one scalar is included in the theory, so that µ L = µ R = 0. We assume that lepton number violating terms are small, implying Λ L,R ≪ 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 µ ′ L and Λ L,R ≪ M X , we have 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 µ ′ 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 ) One option is not to charge either of η L,R implying that both Yukawa coupings are suppressed being lepton number violating.On the contrary, if α to be very small.Depending on the lepton charge assignment of the χ 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(χ L ) = L(χ R ) = 1 5 .This choice is compatible with the C-symmetry discussed earlier in order to avoid diagonal dark photon vertices and forbids the terms Λ ′ R and Λ 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 MeV < m A ′ < 10 GeV and kinetic mixing of order 10 −4 < ε < 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 [19,20,28,67,68].
For a compilation of constraints on dark photons, see Refs.[5,7,69].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 ′ models.The dominant branching ratio (BR) of semi-visible A ′ 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 [10,71] and the LHC [72] can still constrain the direct decays of A ′ into SM particles, A ′ → ℓ + ℓ − , π + π − , π + π − π 0 .These BRs, while still present, are typically much smaller than the BRs into HNFs, as they are of the order of ε 2 α/α D .In addition, when α 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 ′ 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.
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 ε ≲ 0.015 for m A ′ ≲ 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.[75,78].In this work, we use the limits of Ref. [75], where ε ≲ 0.034 for m A ′ < 1 GeV at 95% C.L.
Electroweak precision observables -Among the Electroweak precision observables (EWPO) modified by kinetic mixing, the most important is 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 ε 2 EWPO < 7.3 × 10 −4 at 95% C.L for M A ′ ≪ 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 π 0 → γA ′ [83] and K → πA ′ , with A ′ invisible [84].We update the latter using the latest NA62 measurement of K → πνν [85], including also the dedicated search for π 0 → 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 modelindependent limits on ε 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 ≡ (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, α.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 α to 81 parts per million [91].These measurements are in disagreement at more than 5 σ, 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 ∆a µ -explanation for dark photon masses below m A ′ ∼ 30 MeV (see below).
Muon (g − 2) -In the case of the muon anomalous magnetic moment, the theoretical and experimental uncertainties on a µ ≡ (g − 2) µ /2 are much larger than the uncertainty in α.Therefore, it is not subject to the ambiguities discussed above and can still be more sensitive to new physics due to the a µ /a e ∼ m 2 µ /m 2 e enhance- 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 ′ (see text).In green, we show the preferred region to explain ∆a disp µ , at 3σ, and in solid orange, we show the constraints obtained from the lattice results for ∆a lattice µ .In dashed orange, we also show the region preferred by ∆a lattice µ at 1σ. ment.The theoretical predictions for a µ 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σ [95].The E989 experiment at FNAL [96], running since 2018, has now confirmed the central value measured at BNL, reporting a FNAL 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 µ 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 µ ∼ 300 MeV.
The SM predictions are obtained by combining QED corrections up to 5 loops [89,98], electroweak (a EW µ ) corrections up to 2 loops [99,100], and hadronic contributions including the hadronic vacuum polarization (a HVP ( If the dispersive value holds, the tension between experiment and theory reaches 4.2σ, where ∆a 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 − → hadrons data corroborates this hypothesis.
A deviation in σ(e + e − → hadrons, s) to explain the observed value of ∆a µ was shown to be ruled out for e + e − center-of-mass energies of √ s > 0.7 GeV [125], implying that any missed contributions ought to be mostly concentrated in the π + π − 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 HVP µ on the lattice by the BMW collaboration [131] finds a 2.1σ significant disagreement with the value obtained using the data-driven dispersive method of Ref. [94].Using the BMW result for a HVP µ reduces the disagreement between theory and the experiment average, a comb µ , down to 1.5σ, 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 HVP µ into short, a SD µ , and intermediate, a W µ , time-distance pieces.The data-driven dispersive method [133] and lattice calculations of a SD µ are in good agreement.However, a 3.8σ significant tension exists between the intermediate time-distance observable, a W µ , and all lattice results [131,132,[134][135][136], suggesting the disagreement with e + e − → Hadrons data is, in fact, largest in the energy region of √ s ∼ 1 − 3 GeV [133].
While the nature of this discrepancy is not identified, in this study, we proceed to entertain BSM explanations to both ∆a disp µ and ∆a lattice µ .The 3σ preference region for ∆a disp µ is shown as a green band in Fig. 2, alongside the 1σ preference band and the 3σ exclusion limit from ∆a lattice µ .We summarize other new-physics explanations to ∆a 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 µ − L τ gauge bosons, focusing on final state radiation (FSR) of dark photons, e + e − → µ + µ − Z ′ µ−τ .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 ′ can be emitted as either initial state radiation (ISR) or FSR.In the analogous QED processes, e + e − → µ + µ − γ, 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 ′
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 − → γA ′ process.Prompt decays of A ′ to a pair of HNFs, A ′ → ψ i ψ j , in which the HNFs are: 1. long-lived and decay outside the detector, or 2. short-lived and decay inside to ψ i → ψ j ℓ + ℓ − , with final state leptons pairs whose energies fall below detector thresholds, would appear identical to a monophoton signature accompanied by missing / E T .The strongest bound of this kind is obtained by the BaBar experiment, excluding explanations of (g − 2) µ 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 x represents the set of free model parameters, x = (ε, m A ′ , ∆ 21 , . . .), that have been varied in the analysis performed.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 ), where ψ3 decayed promptly.
e + e − → γA ′ .The search was conducted in the 53 fb −1 dataset collected between 2007-2008 at the center-of-mass (COM) energies Υ(2S), Υ(3S) and Υ(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 − → γA ′ at a COM energy √ s ≈ 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 β z ≈ 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 θ * γ | < 0.6.To a good approximation, the A ′ 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 ′ → e + e − , µ + µ − , and π + π − 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.
Monophoton selection -An e + e − → γA ′ 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.For e + e − pairs with angular separation θ sep.> 10 • , the energy of each electron exceeds BaBar's energy threshold to resolve charged particle tracks, E ± > 100 MeV.For overlapping e + e − pairs with θ sep.< 10 • , we require (E + + E − ) > 100 MeV.
2. The polar angles, θ pol., of the electrons individually, or as a pair, are sufficiently wide that the electrons do not escape along the beam pipeline, 17 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.
FIG. 4. The BaBar limit on the kinetic mixing parameter, ε.On the left panel, we show the limit as a function of our analysis's individual e ± 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.
The threshold is assumed to be a step function with 100% detection efficiency above E 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 ∆ 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 σ BaBar e + e − →γ+inv.(x) = σ e + e − →γA ′ (x)×B (A ′ → inv.) , (35) with B (A ′ → inv.) = 1.To re-interpret the bound for a semi-visible dark photon, we introduce a factor P inv. that accounts for the probability that decays of semi-visible dark photons produced alongside ISR appear invisible and contribute to the monophoton dataset: σ BaBar e + e − →γ+inv.(x) = σ e + e − →γX (x) where X is the semi-visible dark photon.From Eq. ( 35) and Eq. ( 36), we may obtain the relation where ε BaBar is the bound on the kinetic mixing parameter obtained by BaBar.We may define the function P inv. as follows with N being the total number of events that pass the initial monophoton selection, and N 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 ε at each value of (x).The function P 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 2 miss ≡ s − 2E γ √ s where a mild excess has been observed in the region 24 GeV 2 < M 2 miss < 50 GeV 2 .This explanation can fit the events well, explaining the ∼ 2.5σ excess [37].

B. NA64 dark photon searches
The fixed-target experiment NA64 searches for dark sector particles at the CERN SPS, employing a 100 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 [13,18,20,22,23], where N is the target nucleus in the fixed target.Like monophoton searches, the bremsstrahlung signal is proportional to ε 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 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 ′ [23] was conducted on the total sample of electrons on target (EOT) collected during the period 2016-2018, n EOT = 2.84 × 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 ′ with an electron beam at energy E beam = 100 GeV.We consider the electron beam energy to be unaffected by any energy losses happening when entering the ECAL, so that E beam can be considered constant.The beam is considered to impact the ECAL at coordinate x = y = z = 0.The energy of the A ′ is distributed according to the following formula, obtained by applying the improved Weizsaker-Williams (IWW) approximation [139]: where x = E A ′ /E beam , and the A ′ is in the z direction.
After radiating an A ′ , the beam electron (or main electron) will have energy E e = E beam −E A ′ , and will shower inside the ECAL completely.The A ′ decays promptly into a pair of HNFs, which are boosted and rotated to the lab frame according to the A ′ 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: 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 ′ , 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): 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 semivisible 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, l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " / 8 p Q v h u F S w w o I w + K J D z + P l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " / 8 p Q v h u F S w w o I w + K J D z + P l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " / 8 p Q v h u F S w w o I w + K J D z + P l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " / 8 p Q v h u F S w w o I w + K J D z + P

NA64
GeV E e = 100 GeV E e = 100 FIG. 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.
relying on the decay of the heaviest HNF ψ i = ψ 2 or ψ 3 into a single e + e − pair: 1. ψ 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. events with any other activity happening before the fiducial HCAL volume are vetoed; 3. events with more than one visible decay vertex in the fiducial HCAL volume are vetoed.
Additionally, in case ψ 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 ψ 2 lifetime is comparable to the size of the detector.It rapidly loses sensitivity in the high mass region where ψ 2 decays promptly.Other models involving large mass splittings between HNFs are also characterized by a shorter lifetime, due to cτ ∝ ∆ −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 ′ 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. 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. ECAL.The total energy collected, E ECAL , including both the main electron and the e + e − pairs, is compared against the missing energy threshold chosen by the experiment: 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. 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 ε 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 inv .The recast bound can be found by solving the following equation, which matches Eq. ( 37) discussed for the BaBar recast, with x being the set of model parameters:

C. 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 [29,52,140].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, ψ i → ψ j ℓ + ℓ − .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 ψ 2 can be produced in pairs in meson decays or bremsstrahlung.We multiply the signal rate by the respective A ′ vertices and the corresponding long-lived particle decay rate to account for that.For instance, to rescale a limit on ε from an iDM model to a mixed-iDM or i2DM model, we take where Γ are the decay rates deprived of all coupling prefactors.The α 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 ψ i → ψ 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 Γ new ψ2 = Γ iDM ψ2 .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 ′ (known as r in the literature) and ∆ ij ≡ (m i − m j )/m j .We have taken α D ≡ g 2 D /(4π) 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 ′ to be 1/3 for most of the benchmarks, a standard value in the literature.The values ∆ 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 modeldependent search for iDM was performed for the semi-visible dark photon signature of the model.
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, ∆ ij cannot be too large, to minimize the Boltzmann suppression in coannihilations.
Inelastic dark matter (BP1a/b) -The results for the iDM benchmark are shown in Fig. 6, expressed in the ε/m A ′ -plane.In BP1a, we see a significant relaxation of the NA64 and BaBar bounds, with a sizeable ∆a µ preference region now open.Due to a large DS coupling, α 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, ψ 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 ∆a µ 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 ′ 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 [12,143].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 ∆a µ remains open at the 2σ level for m A ′ ∼ 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 ∆a µ 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 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.Increasing the energy threshold makes it more difficult to veto events, leading to a stronger bound.For a threshold larger than 250 MeV, BaBar would cover the entire parameter space allowed by ∆a µ in the case of BP1a, while a threshold of 150 MeV is already sufficient to rule it out completely in BP1b.We note, however, that thresholds exceeding 200 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 α D /m A ′ and ∆ 21 /m A ′ plots, shown respectively in Fig. 8 and in Fig. 9.In Fig. 8, we fix the value of ε to that required to explain ∆a µ and we vary α D and m A ′ .Increasing α D affects the lifetime of ψ 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 ∆ 21 : the right panel is characterized by a lower value.Decreasing ∆ 21 means that the two HNFs become more degenerate and the lifetime of ψ 2 increases, becoming larger than the size of the detector.This effect reduces the amount of ψ 2 decays happening inside the detector and makes FIG. 6.The kinetic mixing ε as a function of the dark photon mass m A ′ for BP1a (left) and BP1b (right) in the inelastic dark matter (iDM) model.We show the ∆aµ-preferred 1, 2, and 3σ 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.the bound more constraining.The right panel is indeed entirely constrained.In Fig. 9 we fix the value of ε to that required to explain ∆a µ and we vary ∆ 21 and m A ′ .In this case, the E137 constraint shown with a thin gray line has been extrapolated at large ∆ 21 value.As discussed, previously, decreasing ∆ 21 makes most of the ψ 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 ′ bound, covering the region in which ε can explain (g − 2) µ .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 ′ , the experiment loses sensitivity, and the bound becomes naturally weaker, independently of the value of ∆ 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 ψ 1 .A larger mass splitting increases the Boltzmann suppression of ψ 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 α D value in BP1b translates into a shift of the relic density line towards lower ∆ 21 , because of its effect in decreasing the annihilation cross section.A smaller ∆ 21 ensures a smaller Boltzmann suppression of ψ 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 ∆a µ and dark matter.
Mixed Inelastic Dark Matter (BP2a/b) -We show the constraints for the mixed-iDM model in Fig. 10, expressed as ε/m A ′ constraints.This model main feature, which is also expected in BP3, is an enhanced A ′ decay BR into ψ 2 ψ 2 .The branching ratio to ψ 2 ψ 1 is suppressed by a factor of the mixing angle α, and the one into two 1 ψ 1 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 ψ 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 ∆a µ and dark matter relic density, along with the constraints discussed, can be achieved in the region 0.9 GeV ≲ m A ′ ≲ 1.2 GeV for BP2a.In the case of BP2b, the coannihilation processes are inefficient due to the smaller α, so that ψ 1 is overabundant in the region of parameter space that can explain ∆a µ .
In addition, we report an analysis on the ∆ 21 and α parameters, showing the constraints in a ∆ 21 /α in Fig. 11, fixing the mass of the dark photon to m A ′ = 1, 1.25, 2 GeV and ε = 0.01, 0.02, ε (g−2) .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 α is expressed by the branching ratio: a larger value favors a larger branching ratio to the ψ 2 ψ 1 channel, with respect to ψ 2 ψ 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 α affects the decay rate of ψ 2 , suppressing its decay, and recovering the original invisible bound.The behavior the constraints is similar to BP1 for what concert ∆ 21 dependence: a larger ∆ 21 means a shorter ψ 2 lifetime, and the energy of its decay is released inside the detector.On the other hand, a lower ∆ 21 resembles the case of fully invisible dark photon decays, as the lifetime becomes larger than the size of the detector.
The relic density lines in Fig. 11 identify the overabundant region, corresponding to high ∆ 21 and low α.In that case, coannihilations and coscattering processes become too inefficient.For small ∆ 21 , the choice of mixing angle does not have a strong impact, as coscattering remains efficient for longer, and the self-annihilation of ψ 2 sets the relic abundance of ψ 1 .For large α the dependence on ∆ 21 on the relic density is relaxed: the enhanced coscattering ratio obtained with a larger α allows to afford a larger ∆ 21 value before ending up in an underabundant scenario.The kink present in the dark matter relic density for m A ′ = 1.0 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σ region accounting for the ∆a µ explanation.
Inelastic Dirac Dark Matter (BP3a-d) -We show the constraints for the mixed-iDM model in Fig. 13, expressed as ε/m A ′ constraints.
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 ψ 2 ψ 2 as it happens in the case of the mixed-iDM model.The rate into the channels ψ 2 ψ 1 and ψ 1 ψ 1 is suppressed by respectively a factor β and β 2 .Differently from the mixed-iDM case, this model allows for dark photon decays to ψ 1 ψ 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 ψ 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 ψ 1 to annihilate through the vertex ψ 1 ψ 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 β parameter, as represented by benchmarks BP3c and BP3d.This choice corresponds to suppress further the channels ψ 2 ψ 1 and ψ 1 ψ 1 .However, it has an impact on the relic abundance of ψ 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 ∆ 21 , because it affects the Boltzmann suppression of the coannihilation ψ 2 .It is interesting to understand the interplay of ∆ 21 and β in the different constraints.In Fig. 14, the bounds are shown in a ∆ 21 /β, fixing the mass of the dark photon to m A ′ = 1, 1.25, 2 GeV and ε = 0.01, 0.02, ε (g−2) .We can draw similar conclusions 11. Parameter space of the mixed-iDM in the plane of ∆21 and the mixing angle α.The parameter m A ′ has been fixed to 1.0, 1.25 and 2.0 GeV (column-wise), while the kinetic mixing ε has been chosen among 0.01, 0.02, and value εg−2 (row-wise), corresponding to the central value of the ∆aµ explanation.BaBar and NA64 recast bounds are shown with the same style used in Fig. 6.
as the ones discussed for the Mixed-iDM model in Fig. 11, with the difference that the x-axis now represents the parameter β.In addition, the CMB constraints are shown: their exclusion region corresponds to large β values, because of the enhancement of ψ 1 ψ 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 ∆ 21 and low β, due to both inefficient coannihilations and suppressed ψ 1 self-annihilations.However, the dependence on β is stronger due to the presence of ψ 1 selfannihilations, which can dominate over coannihilations in depleting the DM density.For this reason, the relic density is underabundant for large β independently on the value of ∆ 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.

Three Heavy Neutral Fermions (BP4a-c, BP5)
Relic density 12. The relic density (solid) in the parameter space of mixed-inelastic DM model in the plane of kinetic mixing ε versus dark photon mass m A ′ for fixed values of α = 0.5, r = m1/m A ′ = 1/3, and ∆21.CMB limits are not applicable as the dark matter self-annihilation ψ1ψ1 → f + f − is forbidden in the C symmetric limit.
-We show the constraints for the Three Heavy Neutral Fermions models in Fig. 16, expressed as ε/m A ′ constraints.
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: 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 dif-ferent values assumed by the parameter r = m 1 /m A ′ .This parameter affects both the HNF lifetime and the A ′ branching ratios.The HNF lifetime depends on r according to cτ ∼ m −1 A ′ 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 ′ → N 3 N 2 , because the kinematics requires: 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 ′ 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
A. 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 selfannihilating 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 ∆a µ .Self-annihilation near the A ′ resonance, r = m 1 /m A ′ ≲ 1/2, can significantly enhance cross sections, but such mass spectrum would leave no room for semi-visible, promptly-decaying fermions.
Coannihilations, ψ 1 ψ 2,3 → (A ′ ) * → f + f − , are therefore the most natural possibility to achieve freeze-out.The coannihilation cross section of opposite C states, ψ i and ψ j , is given by 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 ψ 1 in our bench- marks, 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 [21,30] and i2DM [45].We find a 50% disagreement with the relic curves of Ref. [50] for m 1 ≳ 100 MeV.
Direct detection -Direct detection of a dark matter particle of mass m χ ∼ 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 q 5 (χχ)(qq) [148].In terms of coupling to nucleons, where µ N = m χ m N /(m χ + m N ) is the reduced mass of the DM and the nucleon N , and c N 1 is the loop-induced, DM-nucleon coupling.The matching to nucleon currents gives [149] c where Q q is the quark electric charge, r = m χ /m A ′ , F 3 (x) is the loop function from Ref. [148], and F q/N S (Q 2 ) the nucleon scalar form factors. Approximating the form factors to their Q 2 = 0 value and F q/p S (0) ≃ F q/n S (0) ≃ (15,35,40) MeV for q = (u, d, s), we estimate the elastic DM-nucleon cross sections.At a typical point of parameter space, where we assumed m A ′ = 3m χ = 1 GeV and 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], SEN-SEI [158], XENON1T [159], and DarkSide-50 [160] collaborations to set limits in our mass region of interest.The best constraints come from DarkSide, where σ χN ≲ 0.1 pb at m χ = 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 χ ≡ ψ 1 on electrons.For a heavy dark photon, the total cross section ≃ 4 × 10 −7 pb × sin β 0.14 The leading limits in this parameter space are from XENON1T [161,162], PandaX-II [163], and SENSEI [158].At m χ ∼ 100 MeV, XENON1T constrains σ e ≲ 10 −4 pb, still orders of magnitude above our estimates for β = 8 • .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 β for m A ′ ≳ 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 χ → ψ 2,3,... , 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 χ ≃.This constraint is much weaker and, not significantly constraining for models with co-annihilating dark matter candidates, like iDM and mixed-iDM, and in mod-els 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 β.The self-annihilation cross section to charged leptons in the i2DM model is given by Eq. ( 46), with m ij → m 1 and V ij → β 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 ψ 2 , ψ 3 , . . .can be safely neglected.

B. HNFs as heavy neutral leptons
Having discussed the theoretical aspects of a HNL interpretation in Section II D, 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 ∆a µ explanation does not lead to any constraint on this mixing, the latter will be constrained from below by BBN, such that τ N4 ≲ 0.1 s, and from above by laboratory searches.As highlighted in Refs.[37,176], 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 τ 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,Γ A ′ N4→νe + e − ≫ Γ W, Z N4→SM , where Γ Md is the decay rate mediated by the particle Md.There are two cases: in the long-lived regime, cτ LAB 4 > 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 Γ A ′ N4→νe + e − /L and the bounds get stronger.Alternatively, if N 4 is too short-lived, cτ LAB 4 ≪ 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 µ4 are set by T2K [179] and MicroBooNE [180] (see also PS191 limits [181] and the discussion in Refs.[182,183]), while U τ 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 electronpositron 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) µ anomaly.The particle ψ 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 τ HN L ≪ 1 m cannot be obtained in the standard HNL scenario and require light dark photons and large kinetic mixing.

C. 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 − → A ′ γ, or the production of HNFs through off-shell dark pho-tons, e + e − → (A ′ ) * → ψ i ψ 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 photondark-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 2 X = s − 2E γ √ 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 α D , and so can be large for models where α D ≫ α [37,193].In terms of the cosine of the angle of ψ i with respect to the beam in the COM, c θ , the differential cross section is given by, where ∆ ij = (m i − m j )/m j , E CM = √ s is the center-ofmass collision energy, and Γ A ′ 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 ′ , 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 ψ 3 → ψ 2 → ψ 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 ε < 10 −3 for dark photon precisely in the region of interest for ∆a µ , m A ′ ≳ 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 Then, in both ISR and s-channel events, the dilepton invariant mass of ψ i → ψ j ℓ + ℓ − decays would constrain the mass HNF splittings through the inequality m ℓℓ < ∆ ij m j .In this case, displaced vertices would help isolate the different HNF decay cascades.Finally, with displaced vertices in both ISR and s-channel, the experiment would be able to extract more information on the boosts of the HNFs.The boosts will be larger for s-channel HNFs than in ISR.
We now comment on a short-term possibility that can be pursued.Current published datasets from searches for visible A ′ at KLOE/KLOE-2 [144,145] and BES-III [146] can shed light on semi-visible A ′ by looking for a broad invariant mass spectrum of dileptons on top of their smooth backgrounds.While visible resonances are better constrained due to the smaller backgrounds in a bump hunt, the smooth but often narrow distributions of dilepton invariant masses m ℓℓ can still be searched for.We identify the following channels as promising datasets for a semi-visible analysis: [10,144,194,195], • φ → η(X → e + e − ) [71,196], and [197].[198], where X and Y are some fully visible resonances.We leave the evaluation of these constraints to future literature.
Couplings to the Z boson -In addition to the A ′ coupling to the EM current, one can also explore the SM Z boson coupling to the dark current, shown in Eq. (3).For the large values of kinetic mixing explored here, Z decays can produce HNFs with branching ratios of where we neglected the small mass of the HNFs.While this BR is too small to be constrained by invisible Z decays, it can be used to look for lepton jets, as done at LEP by the DELPHI [189] and L3 [190,191] collaborations.The signature considered at LEP was a single HNL decaying to leptons or quarks, produced alongside a neutrino, Z → νN .This search can, in principle, also be used to constrain semi-visible dark photons, using the channels Z → ψ 1 ψ 2 , where ψ 2 decays either promptly or displaced inside the detector and ψ 1 would constitute missing energy.The limits will be modified due to the small splitting between parent and daughter HNFs, which decreases the dilepton energy.Channels like Z → ψ 2,3,... ψ 2,3,... could be much more common than in the HNL scenario, where they are doubly suppressed by neutrino mixing.We leave a detailed study of this interesting probe for future literature.
Neutrino experiments -Neutrino experiments can test semi-visible dark photon models in two ways.Firstly, in a DM interpretation, the HNFs can be produced in neutral meson decays and bremsstrahlung at the target, and the DM could travel to the detector, where it can coherently interact with nuclei N to produce its coannihilation partners [52,140,199], The decays of the co-annihilators can then produce displaced dileptons.This displacement is especially interesting for multi-component detectors like MINERvA and the near detector of T2K, ND280 [200,201], and can be explored at future experiments like DUNE and Hyper-Kamiokande [202,203].At high energies, experiments like IceCube and KM3NET can search for the upscattering signature by using the atmospheric production of DM particles.The DM particle can then upscatter via deepinelastic scattering and the subsequent decay of ψ 2,3,... would be sufficiently displaced from the vertex to form a double-bang signature.This has been explored in the context of HNLs in Refs.[204,205], but can be adapted to DM particles.Co-annihilators can also be produced at the target alongside ψ 1 .For small mass splittings, they would be long-lived and can be constrained using high-energy beam dumps, like CHARM and NuCAL [52], or searched for at forward or surface detectors at the LHC [206] Their prompt decays, however, would contribute to the flux of ψ 1 .In Ref. [199], the authors study the sensitivity of the SBN program at Fermilab to the production of iDM.The three Liquid Argon detectors placed along the Booster Neutrino Beam (BNB) are sensitive to the production of DM in the BNB, as well as those produced in the NuMI beam, which is located at an off-axis location.In addition, the BNB has the ability to run in off-target mode, directing the proton straight into the beam dump.In this way, the flux of neutrinos going through the detector is minimized due to the absence of focus from the magnetic horns.The sensitivity of SBN to iDM models can improve on CHARM and NuCAL constraints [52], especially in the mass range of interest for our work.
If instead an HNL interpretation is assumed, the decays in flight of the lightest HNF, N 4 = ψ 1 can be searched for.The HNL N 4 can only decay via the small mixing with SM neutrinos, through CC, NC, or dark photon interactions.Typically, for the masses and large values of α D ε 2 considered here, the decays of N 4 will proceed predominantly through the dark photon, and will make the branching ratios into ψ i → ν 1,2,3 ℓ + ℓ − dominate.The leading limits on this type of decay come from the T2K experiment [179].The dark photon contribution increases the decay rate of N 4 , and enhances the decay-in-flight event rate, opening new parameter space between cosmological and laboratory-based limits on the mixing of N 4 with active neutrinos [182].
Also in a HNL interpretation, there is a second possibility to produce the HNLs.As discussed in Section VI B, the scattering of neutrinos in the beam with the dirt or in the detector can produce HNLs, which subsequently decay.Through the exchange of a dark photon, active neutrinos coherently upscatter on nuclei N , The event rate is proportional to the mixing between the HNLs and active neutrinos.The mixing with the muon flavor is the most relevant in this case as most of the flux at accelerator and atmospheric experiments is composed of ν µ and ν µ .While both DIF and upscattering signatures are proportional to parameters that have no impact on the semivisible signatures at collider and fixed-target experiments, they cannot be uniquely determined in the parameter space shown in Section V. Nevertheless, evidence for either of these would indicate that HNFs mix with active neutrinos, ψ i = N i , and confirm an HNL hypothesis.
Kaon factories -Kaon decays can further constrain the parameter of semi-visible dark photons in two ways: i) through the direct production of dark photons, or ii) through the direct production of the HNFs.The latter possibility is, in particular, a powerful probe mixing between neutrinos and the HNFs.As discussed in Section III, dark photons can be produced directly via kinetic mixing in K → πA ′ as well as K + → ℓ + ν ℓ A ′ .The subsequent semivisible decays of A ′ would then lead to multi-lepton final states [207], albeit with at least two invisible particles.Direct production of A ′ via kinetic mixing is, however, suppressed by m 2 A ′ /m 2 K , and has limited reach (cf.Fig. 2).A much more promising channel, however, is the direct production of HNLs through their mixing with the electron or muon flavor.Just as the upscattering signatures discussed in the paragraph above, direct production of HNFs in kaon decays would provide direct evidence for their heavy neutral lepton interpretation, ψ i = N i .As pointed out in Ref. [176], NA62 can use a three-track search to look for the production and the decay products of N i , , where α, β ∈ {e, µ}. ( The striking feature of this signature is the reconstruction of the dark particle masses via the reconstructed quantities, The event rate is proportional to |U αNi | 2 , and the subsequent primary as well as any secondary decays would provide the additional lepton tracks at no additional cost to the rate.Displaced vertices in NA62 can be identified for proper lifetimes of the HNFs as small as cτ 0 ∼ 10 ps thanks to the O(10) cm vertex resolution of NA62.Future fixed target experiments (LDMX) -The next-generation fixed-target experiment LDMX [208,209] provides a unique setup to search for semi-visible signatures.The proposed design is focused on searches for bremsstrahlung-production of dark sector particles, e − N → e − A ′ N , through the missing-momentum technique.Differently from NA64, LDMX aims to measure both the energy and transverse momentum of the recoil electron, having more access to the kinematics of A ′ production.The proposal considers a primary beam of electron of ∼ 4 − 16 GeV at SLAC impinging on a thin target inside a magnetic field [210].The beam would be tracked with low-mass trackers up and downstream of the target and then stopped by a large detector with ECAL and HCAL components, where the total energy of the recoil electron can be measured.
Because the primary electrons will not shower until reaching the ECAL, the production of additional e + e − , µ + µ − , and π + π − pairs in prompt semi-visible A ′ decays would provide a striking signature in the experiment.In contrast to NA64, LDMX offers a lower-energy beam, enlarging its reach in HNF lifetime, and tracking capabilities, allowing the additional tracks to be seen in association with the recoiling electron.While a detailed background study is needed, we note that QED processes like trident production, e − N → e − e + e − N , and hard-photon conversions, e − N → e − (γ brem → e + e − )N , would not carry the missing momentum of the semi-visible signal.In this regard, semi-visible events have an advantage over fully invisible ones because of the high multiplicity of tracks.Finally, we note that the ∆a µ region of interest overlaps with many vector meson V resonances.The mass mixing of V and A ′ can provide an additional and powerful production mechanism of semi-visible HNFs [211].

VII. CONCLUSIONS
Semi-visible dark photons typically arise in rich U (1) D dark sectors with a non-minimal particle content.If a symmetry distinguishes the DS fields from the SM fields, the lightest DS particle provides an ideal candidate for dark matter below the GeV-scale.If no such symmetry is present, the HNFs can mix with the SM light neutrinos and are identified with heavy neutral leptons.In this case, the lightest HNL decays to SM particles with long lifetimes.
In this paper, we have systematically studied a range of models with increasing fermionic content in the dark sector.In particular, we have discussed the role of a chargeconjugation symmetry C in the dark sector, which ensures that the A ′ iψ i γ µ ψ j interactions are predominantly offdiagonal in the i and j indices, generalizing the idea behind the popular iDM model.This is necessary to suppress dark photon decays into two ψ 1 particles, as such decays contribute to the invisible branching ratio of the dark photon, which is severely constrained.As an addition to iDM, we propose the three fermion mixed-iDM model, where a mostly-sterile Majorana DM ψ 1 particle co-annihilates with a mostly-dark and heavier Dirac fermion, Ψ 2 .Due to the C symmetry, the DM can only couple to Ψ 2 through a small mixing α, while Ψ 2 can have O(1) self-couplings.This possibility has not been explored before in the context of DM.Within the three fermion scenario, we also consider more general models of three Majorana particles, with and without enforcing C symmetry.These models typically favor more pronounced mass hierarchies and a HNL interpretation over that of DM, since coannihilations are strongly suppressed.
We follow this with a discussion of the exact Dirac limit of a four Weyl-fermion model, recovering the inelastic Dirac dark matter (i2DM) scenario [45].In this case, a U (1) D -charged Dirac fermion mixes with a sterile Dirac particle.In contrast to the mixed-iDM case, the Dirac DM particle now has self-interactions, albeit suppressed by a small mixing angle β.
If A ′ is heavier than the HNFs, its decay to pairs of HNFs can be followed by their subsequent decays, ψ i → ψ j ℓ + ℓ − or ψ i → ψ j π + π − .As the decays of the A ′ do not lead to any visible resonances, resonance searches in the invariant mass spectra of dilepton pairs are not constraining.For large couplings and kinetic mixing, these HNF decay lengths can be much smaller than the size of a typical particle detector.The presence of visible particles and missing energy within the detector vetoes these semi-visible decays, modifying existing constraints on invisible dark photons.
In order to quantitatively assess these effects, we develop our own fast MC simulation of dark photon production and decay at BaBar and NA64, recasting the bounds from these experiments on the parameter space of semi-visible A ′ models.
We find a significant modification of the bounds on kinetic mixing in the region of m A ′ ∼ 0.3 − 1.3 GeV.This opens up new parameter space at large values of ε ∼ O(10 −3 − 10 −2 ), which has been fully excluded for visible and invisible A ′ .This region is of particular interest as dark photons with masses in this range can explain the discrepancy between dispersive calculations and the measurement of the anomalous magnetic moment of the muon, ∆a disp µ .In both iDM (cf., Fig. 6) and mixed-iDM scenarios (cf.. Fig. 10), we find that the lightest HNF can also constitute a thermal DM candidate.In i2DM, however, the mixing-suppressed self-annihilations of DM are still significant at late times, so CMB constraints exclude the entire ∆a µ region (cf., Fig. 13).
We point out that in the newly-open parameter space, the bremsstrahlung production of A ′ in the fixed target of NA64 can still probe the semi-visible dark photon without the need for displaced decays, the method employed in Ref. [24].By aggregating the additional energy of e + e − pairs produced by short-lived HNFs to the energy of the primary electron beam in the electromagnetic calorimeter, NA64 would be sensitive to the missing energy carried away by stable or long-lived HNFs.Our sensitivity curves show that this new signal definition could cover most of the open ∆a disp µ parameter space in an iDM model.A companion paper derives new NA64 limits and future sensitivity curves for iDM and i2DM models based on this method [38].Our projections are in good agreement with the experimental results given in Ref. [38], where a sophisticated detector simulation is used.
In addition, the newly-open regions also provide a realistic target for e + e − colliders.Re-analyses of existing BaBar, KLOE/KLOE-II and BESS-III data can target HNF production with multiple leptons associated, with or without initial state radiation.This also includes LEP, where the Z boson coupling to the dark current can be constrained using Z → ψ 1 ψ 2,3,... events.In the near future, monophoton searches at Belle-II can improve on the BaBar limits on invisible A ′ .The veto on additional leptons will be even more important in this case due to the improved hermeticity of the detector.Dedicated searches, such as those discussed here and in Refs.[30,31,193], will be essential in constraining a semi-visible A ′ .
Following a possible detection of HNFs in the semivisible decays of A ′ in fixed-target or collider experiments, a key question would arise on whether they are a DM or HNL particles, revealing the presence or absence of additional symmetries in the theory.Direct detection of non-relativistic DM particles would be challenging due to the large mass splittings and inelastic interactions.Boosted DM, produced via the DM upscattering by cosmic rays, could provide an interesting detection avenue to be further investigated, in particular, the potential to exploit large neutrino detectors with low energy thresholds.On the other hand, signatures associated with the HNF mixing with neutrinos would provide decisive evidence for the HNL interpretation and a possible connection to the origin of neutrino masses.The three most promising experimental strategies for this scenario are decay-in-flight searches for ψ 1 → νℓ + ℓ − , neutrino upscattering to ψ 2,3,... , and direct production of ψ i particles in leptonic kaon decays.
A semi-visible option for GeV-scale dark photons keeps the door open for large kinetic mixings that can be directly probed at low-energy experiments.This scenario provides a last chance for the kinetically-mixed dark photon interpretation of the muon ∆a disp µ , which has already been ruled out in the visible and invisible options.The class of semi-visible dark photon models is certainly within present experimental reach and may give us a clue as to whether nature prefers a rich and complex dark sector.

FIG. 1 .
FIG.1.A summary of the benchmark models we consider, including the HNF spectra and their interaction vertices with the dark photon, A ′ .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 (ψ1, ψ2, . . .), where the mixing angles are small α, β ≪ 1, and the entries denoted by × are arbitrary.Mixed-color states have no definite C property.

FIG. 2 .
FIG.2.Model independent limits on kinetic mixing ε 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, π 0 → γA ′ and K + → π + A ′ .Both assume A ′ 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 ′ (see text).In green, we show the preferred region to explain ∆a disp µ = 116 592 040 (54) × 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 BNL µ = 116 592 920 (63) × 10 −11 , provides an experimental average a comb µ = 116 592 061(41) × 10 −11 .
FIG.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 ′ → (ψ2 → ψ1µ + µ − )(ψ3 → ψ2e + e − → ψ1e + e − e + e − ), where ψ3 decayed promptly.
t e x i t s h a 1 _ b a s e 6 4 = " / 8 p Q v h u F S w w o I w + K J D z + P 2 P O e X b 3 S u M 7 j K M I J n M I 5 e H A J D b i F J v j A g M M z v M K b I 5 0 X 5 9 3 5 W L Q W n H z m G P 7 A + f w B n v u N 4 g = = < / l a t e x i t > e < l a t e x i t s h a 1 _ b a s e 6 4 = " / 8 p Q v h u F S w w o I w + K J D z + P 2 P O y R k = " > A A A B 7 H i c b V A 9 T w J B E J 3 D L 8 Q v 1 N J m I 5 j Y S O 5 o t C T a W G L i I Q m c Z G + Z g w 1 7 e 5 f d P R N C + A 0 2 F h p j 6 w + y 8 9 + 4 w B U K v m S S l / d m e X b 3 S u M 7 j K M I J n M I 5 e H A J D b i F J v j A g M M z v M K b I 5 0 X 5 9 3 5 W L Q W n H z m G P 7 A + f w B n v u N 4 g = = < / l a t e x i t >

1 FIG. 7 .
FIG.6.The kinetic mixing ε as a function of the dark photon mass m A ′ for BP1a (left) and BP1b (right) in the inelastic dark matter (iDM) model.We show the ∆aµ-preferred 1, 2, and 3σ 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.

FIG. 8 .
FIG. 8. Parameter space of αD versus m A ′ for the iDM model, fixing to the value needed to explain the ∆aµ, ε = εg−2.The mass splitting is fixed at ∆21 = 0.5 (left) and ∆21 = 0.4 (right).The constraints are shown with the same style used in Fig. 6.

FIG. 10 .
FIG. 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 α mixing angle.

FIG. 13 .
FIG. 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, β.

5 FIG. 14 .
FIG.14.Same as Fig.6but for BP3, corresponding to the inelastic Dirac dark matter (i2DM) model.The parameter space is shown as a function of ∆21 and the mixing angle β.The region dotted in black is excluded by CMB limits, providing a full exclusion of these slices of parameter space.

FIG. 15 .
FIG.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 ε versus dark photon mass m A ′ for fixed values of αD = 0.5, r = m1/m A ′ = 1/3, and ∆21.Each curve corresponds to a different value of the mixing angle β.The CMB limits exclude large values of ε.The ∆aµ 3σ preferred region is shown as a green band.

TABLE I .
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 = m1/m A ′ and ∆ij = (mi − mj)/mj.The dark photon coupling vertices Vij are defined in Eq. (4).