#  MS-TP-22-50 Combined constraints on dark photons and discovery prospects at
the LHC and the Forward Physics Facility

Amin<EMAIL_ADDRESS>, Mohammad Mahdi
<EMAIL_ADDRESS>, Michael Klasena,c333michael.klasen@uni-
muenster.de ,
Pran<EMAIL_ADDRESS>and Zhu-Yao
<EMAIL_ADDRESS>

aInstitut für Theoretische Physik, Westfälische Wilhelms-Universität Münster,
Wilhelm-Klemm-Straße 9, 48149 Münster, Germany
bLaboratoire de Physique Subatomique et de Cosmologie, CNRS/IN2P3,
53 Avenue des Martyrs, 38026 Grenoble, France
cSchool of Physics, The University of New South Wales, Sydney NSW 2052,
Australia
dDepartment of Physics, Northeastern University,
111 Forsyth Street, Boston, MA 02115-5000, USA

###### Abstract

Hidden sectors are ubiquitous in supergravity theories, in strings and in
branes. Well motivated models such as the Stueckelberg hidden sector model
could provide a candidate for dark matter. In such models, the hidden sector
communicates with the visible sector via the exchange of a dark photon (dark
$Z^{\prime}$) while dark matter is constituted of Dirac fermions in the hidden
sector. Using data from collider searches and precision measurements of SM
processes as well as the most recent limits from dark matter direct and
indirect detection experiments, we perform a comprehensive scan over a wide
range of the $Z^{\prime}$ mass and set exclusion bounds on the parameter space
from sub-GeV to several TeV. We then discuss the discovery potential of an
$\mathcal{O}$(TeV) scale $Z^{\prime}$ at HL-LHC and the ability of future
forward detectors to probe very weakly interacting sub-GeV $Z^{\prime}$
bosons. Our analysis shows that the parameter space in which a $Z^{\prime}$
can decay to hidden sector dark matter is severely constrained whereas limits
become much weaker for a $Z^{\prime}$ with no dark decays. The analysis also
favors a self-thermalized dark sector which is necessary to satisfy the dark
matter relic density.

###### Contents

  1. 1 Introduction
  2. 2 Dark photons in the Stueckelberg extension of the Standard Model
  3. 3 Experimental constraints
    1. 3.1 Fiducial measurements of SM processes using Contur
    2. 3.2 LHC searches: dijet, dilepton and monojet limits
    3. 3.3 BaBar, electron bremsstrahlung and beam dump experiments
    4. 3.4 Relic density constraints
    5. 3.5 Dark matter direct detection
    6. 3.6 Dark matter indirect detection
  4. 4 Exclusion limits
    1. 4.1 Heavy $Z^{\prime}$ bosons
    2. 4.2 Light $Z^{\prime}$ bosons
  5. 5 Discovery potential of the Stueckelberg $Z^{\prime}$ boson at the HL-LHC
    1. 5.1 Cut-and-count analysis
    2. 5.2 Boosted decision tree analysis
  6. 6 Detection of a Stueckelberg dark photon at the Forward Physics Facility
  7. 7 Conclusion
  8. A Rotation angles and dark matter couplings
  9. B More on the exclusion plots from Contur

##  1 Introduction

The identification of dark matter (DM) is one of the most important current
problems in astroparticle physics, given the overwhelming observational
evidence for it on many different length scales. Aside from visible sector
candidates for DM such as the weakly interacting massive particle (WIMP), one
may have dark matter arising from the hidden sector. Here we discuss
experimental limits on a model containing such a candidate that arises from
the extension of the electroweak sector of the Standard Model (SM) along with
an additional $U(1)_{X}$ gauge group of the hidden sector which kinetically
mixes with the gauge field of the hypercharge $U(1)_{Y}$ [2, 3]. For this, we
consider a well motivated extension with the Stueckelberg mechanism as the
source of mass generation for the extra gauge boson. As noted, the
communication between the hidden and the visible sectors arises due to kinetic
mixing of $U(1)_{X}$ and $U(1)_{Y}$. In this case, transition to the canonical
basis requires diagonalization of a $3\times 3$ mass matrix for the gauge
fields involving $U(1)_{X},U(1)_{Y}$ and the gauge field for the neutral
component of the group $SU(2)_{L}$. In the canonical basis one finds one
massless field which is the photon and two massive fields which can be
identified as the $Z$ boson and a massive dark field. This massive dark field
can be viewed as a dark photon, $Z^{\prime}$, since it is associated with a
$U(1)_{X}$ factor and mixes with the SM photon. In this work, $Z^{\prime}$
will take on a wide range of masses, from sub-GeV to multi-TeV values. In the
canonical basis, $Z^{\prime}$ will have interactions with the visible sector
quarks and leptons while the $Z$ boson will also have interactions with the
dark fermions which carry $U(1)_{X}$ quantum numbers. Thus communication
exists between the hidden and the visible sectors due to the exchange of
$Z^{\prime}$ and the exchange of $Z$. Several recent works in this framework
can be found in [4, 5, 6, 7]. For an overview of heavy $Z^{\prime}$ models see
ref. [8] and of dark photons see ref. [9].

We will discuss the Stueckelberg model in more detail in section 2. However,
here we give a brief review of the mechanism for easy reference. Thus, let us
consider the Lagrangian with an abelian vector boson $C_{\mu}$ coupled to a
pseudo-scalar $\sigma$ so that

$\displaystyle{\cal
L}=-\frac{1}{4}{C}_{\mu\nu}{C}^{\mu\nu}-\frac{1}{2}(mC_{\mu}+\partial_{\mu}\sigma)(mC^{\mu}+\partial^{\mu}\sigma)\
.$ (1.1)

This Lagrangian is gauge invariant under the transformations $\delta
C_{\mu}=\partial_{\mu}\epsilon$ and $\delta\sigma=-m\epsilon.$ Using the gauge
fixing term ${\cal L}_{\rm gf}=-\frac{1}{2\xi}\left(\partial_{\mu}A^{\mu}+\xi
m\sigma\right)^{2}$, the resulting Lagrangian has the form

$\displaystyle{\cal L}_{\rm
tot}=-\frac{1}{4}{C}_{\mu\nu}{C}^{\mu\nu}-\frac{m^{2}}{2}C_{\mu}C^{\mu}-\frac{1}{2\xi}(\partial_{\mu}C^{\mu})^{2}-\frac{1}{2}\partial_{\mu}\sigma\partial^{\mu}\sigma-\xi\frac{m^{2}}{2}\sigma^{2}.$
(1.2)

Here one finds that $C_{\mu}$ is massive and decoupled from $\sigma$. We can
include couplings with matter in the usual way using the interaction term
${\cal L}_{\rm int}=gA_{\mu}J^{\mu}$ along with the constraint
$\partial_{\mu}J^{\mu}=0$.

As noted earlier, in this work we will consider a Stueckelberg extension of
the SM. In the analysis, we use the coupling that $Z^{\prime}$ has with the SM
fermions to derive constraints on the model parameters based on searches at
the LHC as well as precision measurements of SM processes using Contur, a new
tool to set exclusion limits. Further, DM searches from direct and indirect
detection experiments add more constraints which also depend on the coupling
between $Z^{\prime}$ and the dark fermions. The analysis covers a wide range
of the $Z^{\prime}$ mass, from the sub-GeV to multi-TeV mass range, and we
discuss the parameter space remaining after the multitude of constraints from
a variety of experiments have been imposed. We also perform a detailed
analysis for a potential discovery of $Z^{\prime}$ at the high luminosity LHC
(HL-LHC) and discuss the mass reach of future forward detectors in the sub-GeV
mass region of the $Z^{\prime}$ boson. Several precision calculations have
been performed for LHC production of a heavy $Z^{\prime}$ [10, 11] and for
$t\bar{t}$ production with $Z^{\prime}$ and $W^{\prime}$ [12, 13, 14]. In this
work, our calculation is done at the NLO level without including the
resummation effects.

The outline of the remainder of the paper is as follows: In section 2 we
define the hidden sector model and its communication with the visible sector
and specifically of the $Z^{\prime}$ with the Standard Model particles and
with the dark fermion in the hidden sector. In section 3 we list the different
collider and dark matter direct and indirect detection experiments whose
limits are recasted and used to constrain our model which we show in section
4. In section 5 we give a detailed LHC analysis for a potential discovery of a
TeV mass $Z^{\prime}$. The sensitivity reach at forward detectors for a sub-
GeV dark photon is discussed in section 6. Conclusions are given in section 7.
Further details related to the model are given in Appendix A while exclusion
limits from Contur are shown in Appendix B.

##  2 Dark photons in the Stueckelberg extension of the Standard Model

We give now a brief account of the extension of the electroweak sector of the
Standard Model with an extra $U(1)_{X}$. The gauge content and matter content
of this sector consist of the gauge field $C_{\mu}$ and the dark Dirac fermion
$D$. We will assume a kinetic mixing of the $U(1)_{X}$ gauge field with the
hypercharge gauge field of the SM as well as a Stueckelberg mass growth of the
$U(1)_{X}$ gauge field and for generality we allow a mass mixing of the
$U(1)_{X}$ gauge field with the $U(1)_{Y}$ field. Thus, the extended
electroweak sector has the following Lagrangian [15]

$\displaystyle\mathcal{L}=\mathcal{L}_{\rm SM}+\Delta\mathcal{L},$ (2.1)

where $\Delta\mathcal{L}$ is the extended part of the Lagrangian given by

$\displaystyle\Delta\mathcal{L}=$
$\displaystyle-\frac{1}{4}C_{\mu\nu}C^{\mu\nu}+i\bar{D}\gamma^{\mu}\partial_{\mu}D-m_{D}\bar{D}D$
$\displaystyle-\frac{\delta}{2}C_{\mu\nu}B^{\mu\nu}-\frac{1}{2}(\partial_{\mu}\sigma+M_{1}C_{\mu}+M_{2}B_{\mu})^{2}$
$\displaystyle+g_{X}Q_{X}\bar{D}\gamma^{\mu}DC_{\mu}.$ (2.2)

Here the first line gives the free part of the Lagrangian for the $U(1)_{X}$
gauge boson $C^{\mu}$ and the hidden sector Dirac fermion $D$, the second line
gives the kinetic mixing and the Stueckelberg mass mixing for $C^{\mu}$ with
the hypercharge gauge field $B^{\mu}$, and the last line gives the interaction
of $C^{\mu}$ with the Dirac fermion $D$. To obtain the mass eigenstates, we
diagonalize $\Delta\mathcal{L}$ along with the Standard Model mass matrix for
the gauge fields $B^{\mu}$ and $A^{\mu}_{3}$, where the latter is the third
component of the $SU(2)_{L}$ gauge field $A^{\mu}_{a}$ ($a=1,2,3$). This
requires diagonalizing a $3\times 3$ matrix involving the fields
$C^{\mu},B^{\mu},A^{\mu}_{3}$ and in the eigen-frame where both the kinetic
and the mass squared matrices are diagonalized, one has the gauge bosons
$A^{\mu}_{\gamma},Z^{\mu},Z^{\prime\mu}$ corresponding to the particles
$\gamma,Z,Z^{\prime}$ where $\gamma$ is the photon, $Z$ is the $Z$-boson, and
$Z^{\prime}$ is the dark photon. In the canonically diagonalized frame the new
interactions are given by

$\displaystyle\Delta\mathcal{L}_{\rm int}$
$\displaystyle=\bar{D}\gamma^{\mu}(g_{Z^{\prime}}Z^{\prime}_{\mu}+g_{Z}Z_{\mu}+g_{\gamma}A_{\mu}^{\gamma})D+\frac{g_{2}}{2\cos\theta}\bar{\psi}_{f}\gamma^{\mu}(v^{\prime}_{f}-\gamma_{5}a^{\prime}_{f})Z^{\prime}_{\mu}\psi_{f}.$
(2.3)

Here $f$ stands for SM quarks and leptons. For the case when the kinetic
mixing or Stueckelberg mass mixing is small, one has $g_{Z^{\prime}}\simeq
g_{X}Q_{X}$, so $g_{Z^{\prime}}$ is of normal size. It is given together with
the couplings $g_{Z}$ and $g_{\gamma}$ in Appendix A. The dark photon can have
vector and axial vector couplings with the fermions $f$ in the visible sector
(quarks and leptons)

$\displaystyle v^{\prime}_{f}$ $\displaystyle=-\cos\psi[(\tan\psi-
s_{\delta}\sin\theta)T_{3f}-2\sin^{2}\theta(-s_{\delta}\csc\theta+\tan\psi)Q_{f}],$
(2.4) $\displaystyle a^{\prime}_{f}$ $\displaystyle=-\cos\psi(\tan\psi-
s_{\delta}\sin\theta)T_{3f}.$

Here $s_{\delta}=\sinh\delta$, $T_{3f}$ is the third component of isospin,
$Q_{f}$ is the electric charge for the fermion $f$ and $\theta$ and $\psi$ are
angles defined in Appendix A. Here we also note that the couplings of
$Z_{\mu}$ and $A^{\gamma}_{\mu}$ in the canonically diagonalized basis are
also modified and are given by [3]

$\Delta\mathcal{L}^{\prime}_{\rm
SM}=\frac{g_{2}}{2\cos\theta}\bar{\psi}_{f}\gamma^{\mu}\Big{[}(v_{f}-\gamma_{5}a_{f})Z_{\mu}\Big{]}\psi_{f}+e\bar{\psi}_{f}\gamma^{\mu}Q_{f}A^{\gamma}_{\mu}\psi_{f}\,,$
(2.5)

where the modifications appear in the definition of the vector coupling
$v_{f}$ and the axial-vector coupling $a_{f}$ which are given by

$\displaystyle v_{f}$
$\displaystyle=\cos\psi[(1+s_{\delta}\tan\psi\sin\theta)T_{3f}-2\sin^{2}\theta(1+s_{\delta}\csc\theta\tan\psi)Q_{f}],$
(2.6) $\displaystyle a_{f}$
$\displaystyle=\cos\psi(1+s_{\delta}\tan\psi\sin\theta)T_{3f}.$

Here Eqs. (2.4) and (2.6) are written with the assumption of zero mass mixing,
i.e. $M_{2}=0$, which we take to be the case throughout this work. We note in
passing that in the literature the dark photon refers to a vector boson from
an extra $U(1)$ gauge field that has kinetic mixing with a massless SM photon.
In this case the dark photon does not mix with the $Z$ boson and so has no
coupling to neutrinos. The Stueckelberg analysis is different in that here,
the extra $U(1)$ gauge field mixes with the SM hypercharge $U(1)_{Y}$ gauge
field and since $U(1)_{Y}$ gauge field has mixing with the gauge field for the
neutral component of the $SU(2)_{L}$ triplet due to the SM Higgs mechanism one
has in general a mass matrix which mixes the three neutral fields. The
resulting massive $Z^{\prime}$ has couplings to all SM fermions just as the
$Z$ boson does. This is mentioned in order to draw a distinction between the
two approaches. For related works and extensions see [16, 17, 18, 19, 20].

##  3 Experimental constraints

In this Section we list the different experiments and methods used in our
analysis to constrain the parameter space of the Stueckelberg dark photon
model. The experimental limits we use pertain to the $Z^{\prime}$ as well as
to the dark fermion DM which we turn into constraints on the model parameter
space in the kinetic mixing-dark photon mass plane. The constraints are also
sensitive to the $Z^{\prime}$ decay channels, i.e., visible decays to SM
fermions and invisible decays to the dark fermion which we will refer to as
dark decays.

The scan of the model parameter space and calculation of the relevant
observables requires a tool chain which we describe next. The Stueckelberg
extension of the SM with a dark sector is implemented in SARAH [21, 22] which
produces the necessary SPheno [23, 24] model files for spectrum generation.
SARAH also automatically writes the CalcHep/CompHep [25, 26] files used by
micrOMEGAs [27] and the UFO [28] files needed by MadGraph5_aMC@NLO [29]. There
are more necessary tools that have been used which we will refer to later in
the paper.

###  3.1 Fiducial measurements of SM processes using Contur

Hundreds of differential cross section measurements were performed during the
first two runs of the LHC. These measurements, even though intended to test
the SM of particle physics, can still be used to investigate physics beyond
the SM (BSM) due to their model-independent nature. A software called Contur
or “Constraints On New Theories Using Rivet” [30, 31] scans the SM analyses
implemented in the Rivet (Robust Independent Validation of Experiment and
Theory) toolkit [32] to check whether a BSM signal is already excluded, and if
so, at which significance. Contur performs a $\chi^{2}$ test statistic to
evaluate the likelihood of the BSM model taking the experimental uncertainties
into account. The considered hypotheses are SM-only and SM+BSM. The CLs
technique [33, 34] is then used to derive the confidence level exclusion on
the BSM theory at a given point of the parameter space. This procedure is then
repeated for each set of parameter values resulting in a map of CL exclusions
(more information about Contur and its statistical method used to perform the
exclusion procedure can be found in [35]). Many studies were performed using
the Contur toolkit to check different BSM scenarios [36, 37, 38, 39] and it
was shown that limits from precision measurements of SM processes can be more
constraining than BSM searches at the LHC in some parts of the parameter space
[40]. In this work, we use Contur to set limits on the parameter space of our
model. The signal cross section was calculated at leading order using Herwig
[41] for a center of mass energy of $7$, $8$, and $13$ TeV. Details about the
LHC analyses that contributed to the exclusions are given in Appendix B. We
consider four cases in which the $Z^{\prime}$ is either heavy or light, and
whether it decays only to visible final states or to dark fermions as well.
Furthermore, since $Z^{\prime}$ mixes with $Z$, contributions to the $Z$ boson
mass and width are expected, so additional LEP constraints are included [42,
43, 44].

###  3.2 LHC searches: dijet, dilepton and monojet limits

We call $Z^{\prime}$ heavy if its mass is greater than the $Z$ boson mass and
light if it’s smaller. Models with an extra neutral gauge boson $Z^{\prime}$
are extensively tested and constrained by a myriad of collider searches. ATLAS
and CMS experiments have searched for a heavy $Z^{\prime}$ resonance and set
stringent bounds on the ratio $m_{Z^{\prime}}/g_{X}$, with $g_{X}$ being the
gauge coupling for a particular $U(1)_{X}$ extension of the SM. Searches for
light $Z^{\prime}$ in the mass range of 1 to $\sim 80$ GeV have been carried
out by CMS, LHCb and BaBar. Lighter masses have also been investigated at beam
dump experiments.

For heavy $Z^{\prime}$, ATLAS and CMS collaborations have looked for an excess
of events in the dijet invariant mass corresponding to the decay of a heavy
vector resonance [45, 46, 47, 48, 49, 50, 51, 52, 53] with up to 139 fb-1 of
data at 13 TeV. Since no significant excess has been found, limits were set on
the coupling and $Z^{\prime}$ mass based on a simplified model with a
Lagrangian

$\displaystyle{\cal L}_{\rm simp}=$
$\displaystyle-\frac{1}{4}F_{\mu\nu}F^{\mu\nu}-\frac{1}{2}M^{2}A_{\mu}A^{\mu}-\bar{\psi}(\frac{\hbar}{i}\gamma^{\mu}\partial_{\mu}+m_{D})\psi$
$\displaystyle+A_{\mu}\bar{\psi}(g_{V}+\gamma_{5}g_{A})\gamma^{\mu}\psi+g_{q}A_{\mu}\bar{q}\gamma^{\mu}q,$
(3.1)

where $A_{\mu}$ is the new massive vector boson field with mass $M$ and
$F_{\mu\nu}$ is its field strength, $\psi$ is the Dirac fermion with mass
$m_{D}$, $g_{V}(g_{A})$ are the vector (axial vector) coupling of $A_{\mu}$
with the Dirac fermion and $g_{q}$ is the coupling of the vector boson with
the SM quarks. The relevant part of the Lagrangian here is the term
$g_{q}A_{\mu}\bar{q}\gamma^{\mu}q$. A comparison between the theoretical and
observed cross sections using the simplified model of Eq. (3.1) is translated
to constraints on $g_{q}$. We recast the obtained limits to our model
parameters as constraints on the kinetic mixing $\delta$ in the case where
$Z^{\prime}$ does not decay to dark fermions, i.e., $m_{Z^{\prime}}<2m_{D}$.
To do so, we follow the procedure in refs. [54, 55] and construct the log-
likelihood

$\ln\mathcal{L}=-2\left[\frac{(v_{f}^{\prime 2}+a_{f}^{\prime
2})^{2}\times\text{BR}(Z^{\prime}\to q\bar{q})^{2}}{g_{q}^{4}}\right],$ (3.2)

where we have ignored the term in the vertex proportional to
$m_{f}^{2}/m_{Z^{\prime}}^{2}$. The terms in the numerator of Eq. (3.2) have a
non-trivial dependence on the kinetic mixing and is determined as the one
minimizing the log-likelihood function.

Unlike the simplified model of Eq. (3.1), our model is not leptophobic and we
have an important decay channel to leptons. So our $Z^{\prime}$ can have
dilepton decays which introduce stringent constraints on the model parameter
space. The dilepton channel is cleaner than the dijet since the latter is
contaminated by large QCD multijet background. Therefore, the dilepton
constraints are much more severe. In constraining the Stueckelberg model, we
use the most recent ATLAS dilepton search with Drell-Yan processes [56]. In
recasting limits from dilepton searches, we include possible interference
effects with SM processes involving $Z$ and $\gamma$ mediators using a
modified version of the code ZPEED [57].

For $m_{Z^{\prime}}\geq 2m_{D}$, the $Z^{\prime}\to D\bar{D}$ channel opens up
which means that the branching ratios to leptons and quarks become smaller. So
we expect in this case weaker dilepton and dijet limits. However, monojet
searches, $pp\to D\bar{D}+$jet, become relevant, i.e., missing energy
recoiling against a hard jet. For this, we use the most recent ATLAS and CMS
monojet searches [58, 59] to constrain the model parameter space. We updated
the monojet module in micrOMEGAs [60] with the most recent ATLAS and CMS data
and used it as our recasting tool.

For light $Z^{\prime}$, searches in the dimuon channel was carried out by CMS
[61] as well as LHCb [62, 63], where the latter investigated prompt and long-
lived $Z^{\prime}$. The null results from these experiments are translated
into constraints on the kinetic mixing coefficient for a dark photon model,
i.e., a $Z^{\prime}$ which kinetically mixes only with the SM photon. To
recast those limits to our model, we implemented the Stueckelberg model in
DarkCast [64, 65] which has a large repository of the most up-to-date searches
on dark photons.

###  3.3 BaBar, electron bremsstrahlung and beam dump experiments

Along with the LHC constraints, many other experiments have investigated a
light $Z^{\prime}$ and set constraints on the kinetic mixing in a simple dark
photon model. BaBar analyzed $Z^{\prime}$ production and decay to visible
final states from $e^{+}e^{-}$ annihilation, $e^{+}e^{-}\to Z^{\prime}\to
e^{+}e^{-}~{}(\mu^{+}\mu^{-})$ [66] and to invisible final states [67].
Furthermore, electron bremsstrahlung experiments such as APEX [68] and A1 [69]
studied the production and decay of $Z^{\prime}$ to $e^{+}e^{-}$ while NA64
[70] studied invisible decays of $Z^{\prime}$. Beam dump experiments such as
E137 [71], E141 [72], E774 [73], KEK [74], and Orsay [75] have studied long-
lived dark photons. All these limits are part of DarkCast which we use as our
recasting tool in this mass range.

###  3.4 Relic density constraints

Not only is our scan of the parameter space over a wide range of $Z^{\prime}$
mass but also over a wide range of the kinetic mixing coefficient. This adds a
complexity related to the fact that for small $\delta$, the dark sector,
comprised of dark fermions and $Z^{\prime}$, may not be in thermal equilibrium
with the SM. This means that one can assume the standard treatment of the
freeze-out scenario adopted in codes like micrOMEGAs and darkSUSY [76] only in
the case where $\delta$ is large enough to maintain thermal equilibrium. The
Boltzmann equation for the number density of $D$ is given by

$\frac{dn_{D}}{dt}+3Hn_{D}=C[f_{D}],$ (3.3)

where $C[f_{D}]$ is the collision term containing DM number-changing processes
such as $D\bar{D}\to f\bar{f}$ and $D\bar{D}\to Z^{\prime}Z^{\prime}$. For DM
lighter than the $Z^{\prime}$ mediator, the freeze-out DM relic density is set
by the annihilation processes into SM fermions. This processes is proportional
to $\sim(\delta g_{X})^{2}$. However, for DM heavier than $Z^{\prime}$, the
process $D\bar{D}\to Z^{\prime}Z^{\prime}$, which is proportional to
$g_{X}^{4}$, becomes kinematically accessible and will set the final DM relic
density. In the freeze-out scenario and considering a Maxwell-Boltzmann phase
space distribution, the collision term is given by

$C[f_{D}]=-\langle\sigma v\rangle(n_{D}^{2}-n_{D\,\text{eq}}^{2}),$ (3.4)

with the thermally averaged cross section given by

$\langle\sigma
v\rangle=\frac{1}{K_{2}(x)^{2}}\int_{1}^{\infty}d\tilde{s}\,4x\sqrt{\tilde{s}}\,(\tilde{s}-1)K_{1}(2\sqrt{\tilde{s}}x)\,\sigma_{D\bar{D}\to
XX},$ (3.5)

where the dimensionless parameters are $x=m_{D}/T$ and
$\tilde{s}=s/(4m_{D}^{2})$. We use micrOMEGAs [27, 60] to determine the DM
relic density in the region where the pure freeze-out mechanism is valid. In
our scan, we accept points whose relic density is less than or equal to that
measured by the Planck collaboration [77]

$(\Omega h^{2})_{\rm Planck}=0.120\pm 0.001,$ (3.6)

i.e., $f_{\rm DM}=(\Omega h^{2})_{D}/(\Omega h^{2})_{\rm Planck}\leq 1$, which
keeps the door open for multi-component DM scenarios.

In the very small kinetic mixing regime, the dark species will never reach
thermal equilibrium with the SM sector. Despite this, annihilation processes
of the type $f\bar{f}\to D\bar{D}$ and $f\bar{f}\to Z^{\prime}$ can gradually
populate the dark sector and set the DM relic density via the freeze-in
mechanism [78]. The situation becomes more involved if the coupling $g_{X}$
among the dark species becomes large enough so that the dark sector reaches
thermal equilibrium, i.e., $D$ and $Z^{\prime}$ enter thermal equilibrium. In
this case, the processes $D\bar{D}\leftrightarrow Z^{\prime}Z^{\prime}$ become
important and one needs to track the number density of $Z^{\prime}$ as well.
Therefore, the calculation of the relic density now requires solving the
coupled Boltzmann equations

$\displaystyle\frac{dn_{D}}{dt}+3Hn_{D}$
$\displaystyle=-\frac{1}{2}\langle\sigma v\rangle_{D\bar{D}\to
f\bar{f}}(n_{D}^{2}-n_{D\,\text{eq}}^{2})-\frac{1}{2}\langle\sigma
v\rangle_{D\bar{D}\to
Z^{\prime}Z^{\prime}}\left(n_{D}^{2}-n_{D\,\text{eq}}^{2}\frac{n_{Z^{\prime}}^{2}}{n_{Z^{\prime}\,\text{eq}}^{2}}\right)$
$\displaystyle~{}~{}~{}-\frac{1}{2}\langle\sigma v\rangle_{D\bar{D}\to
Z^{\prime}}n_{D}^{2}+\langle\Gamma\rangle_{Z^{\prime}\to
D\bar{D}}n_{Z^{\prime}},$ (3.7)
$\displaystyle\frac{dn_{Z^{\prime}}}{dt}+3Hn_{Z^{\prime}}$
$\displaystyle=-\langle\sigma v\rangle_{Z^{\prime}Z^{\prime}\to
D\bar{D}}\left(n_{Z^{\prime}}^{2}-n_{Z^{\prime}\,\text{eq}}^{2}\frac{n_{D}^{2}}{n_{D\,\text{eq}}^{2}}\right)+\frac{1}{2}\langle\sigma
v\rangle_{f\bar{f}\to Z^{\prime}}n_{f\text{eq}}^{2}$
$\displaystyle~{}~{}~{}-\langle\Gamma\rangle_{Z^{\prime}\to
f\bar{f}}n_{Z^{\prime}}-\langle\Gamma\rangle_{Z^{\prime}\to
D\bar{D}}\left(n_{Z^{\prime}}-n_{Z^{\prime}\,\text{eq}}\frac{n_{D}^{2}}{n_{D\,\text{eq}}^{2}}\right).$
(3.8)

This treatment is not part of micrOMEGAs freeze-in routine [79] and so we use
our own numerical calculations with the help of MATLAB ode15s to determine the
DM relic density in the case when thermal equilibrium cannot be guaranteed.
The criteria we use to make this judgment is based on comparing the DM
annihilation rate, $n_{D}\langle\sigma v\rangle_{D\bar{D}\to f\bar{f}}$,
and/or the $Z^{\prime}$ decay rate, $\langle\Gamma_{Z^{\prime}\leftrightarrow
f\bar{f}}\rangle$, to the Hubble parameter $H(T)$. If the rates of both or any
of these processes are larger than $H(T)$ then a thermal equilibrium is
established between the two sectors and micrOMEGAs’s freeze-out routine is
able to handle this scenario.

One final comment regarding the validity of Eqs. (3.4) and (3.8) is in order.
In writing the Boltzmann equations in terms of the number density, one assumes
a well-defined phase space distribution (a Maxwell-Boltzmann distribution in
this case). This is justified if kinetic equilibrium can be maintained till
after chemical decoupling of DM species. In the standard freeze-out scenario,
kinetic equilibrium is maintained through efficient elastic scattering of DM
with SM particles. However, when the coupling between DM and the visible
sector becomes small, elastic scattering may become inefficient causing early
kinetic decoupling. In our model, kinetic equilibrium can still hold and this
is attributed to the dark sector itself. Once produced, dark matter self-
interactions can bring their momentum distribution to a thermal distribution
thus allowing one to use the Boltzmann equations for number density (see Fig.
2 in ref. [80]). Furthermore, DM-$Z^{\prime}$ elastic scattering is strong
enough (owing to the large gauge coupling $g_{X}$ in the dark sector) to also
keep $Z^{\prime}$ in kinetic equilibrium. One can then avoid solving the full
phase space Boltzmann equations and instead consider the number density
equations.

###  3.5 Dark matter direct detection

Experiments on DM direct detection involve the scattering of a DM particle off
the nucleus of a heavy material such as xenon. The voluminous experimental
apparatus operates for a period of time looking for an excess of nuclear
recoil events as a result of DM-nucleon spin-independent (SI) or spin-
dependent (SD) scattering. The differential event rate is given by

$\frac{dR}{dE_{R}}(E_{R})=\frac{\rho_{\rm DM}}{2m_{D}\mu^{2}_{r}}f_{\rm
DM}\sigma_{\rm Dn}F^{2}(E_{R})\eta(E_{R},t),$ (3.9)

where $\rho_{\rm DM}$ is the DM density at the Sun’s location, $\mu_{r}$ is
the DM-nucleon reduced mass, $F(E_{R})$ is the nuclear form factor,
$\eta(E_{R},t)$ contains all astrophysical information

$\eta(E_{R},t)=\int_{|\mathbf{v}|>v_{\rm
min}}\,d^{3}\mathbf{v}\frac{1}{v}f(\mathbf{v}),$ (3.10)

with $f(\mathbf{v})$ the DM velocity distribution and

$v_{\rm min}=\sqrt{\frac{m_{N}E_{R}}{2\mu^{2}_{r}}}.$ (3.11)

The DM-nucleon cross section, $f_{\rm DM}\sigma_{\rm Dn}$, refers to either
the SI or SD cross section. The SD cross section in our model is suppressed
compared to the SI one. Furthermore, experimental constraints on SD are much
weaker than those on SI and therefore we do not discuss them any further. We
calculate $f_{\rm DM}\sigma_{\rm SI}$ and the total event rate using
micrOMEGAs [81] which we also use to recast limits from several experiments
such as CDMSlite [82], CRESST-II [83], CRESST-III [84], DarkSide 50 [85], LUX
2016 [86], PICO-60 [87, 88], PandaX [89, 90] and Xenon1T [91]. We also take
into account the most recent limits from LUX-ZEPLIN (LZ) [92]. Note that there
are two diagrams contributing to the SI cross section: the $Z$ and
$Z^{\prime}$ exchange diagrams. Since setting $M_{2}=0$ would prevent DM from
acquiring a millicharge, diagrams with photon exchange are absent as well as
the ones with a Higgs since, unlike a Higgs portal model, our DM has no
coupling to the Higgs.

###  3.6 Dark matter indirect detection

Even though the DM particles have achieved a constant comoving number density,
annihilation of DM particles into the SM can still happen today especially in
regions with large density. The annihilation processes can result in $\gamma$
ray emissions, charged particles (such as electrons, positrons or even
composite particles such as antiprotons and antideuterons) as well as
neutrinos. The detection of charged particles cannot be easily attributed to
DM annihilation since it is difficult to trace back the origin of those
particles as they are deflected by magnetic fields. However, $\gamma$ ray
photons are not affected by magnetic fields and can be an important tool in DM
indirect detection. The well known gamma ray excess at the center of our
galaxy [93] is an example but it remains debatable as the galaxy center is a
rich source of gamma rays from other astrophysical sources.

As DM annihilate into SM particles, charged final states can radiate off
photons which are the source of prompt $\gamma$ ray emission. Photons can also
come from the decay of pions after final state quarks have hadronized. The
differential photon flux due to DM annihilation for an observation region
$d\Omega$ is given by

$\frac{d\phi}{dE\,d\Omega}(E_{\gamma})=\frac{r_{\odot}}{4\pi}\frac{1}{4}\left(\frac{\rho_{\rm
DM}}{m_{D}}\right)^{2}J\sum_{i}f^{2}_{\rm DM}\langle\sigma
v\rangle_{i}\frac{dN_{i}}{dE_{\gamma}},$ (3.12)

where $r_{\odot}$ is the location of the Sun in the galactic plane, the $J$
factor contains astrophysical information and the photon spectrum
$dN_{i}/dE_{\gamma}$ due to annihilation to some final state $i$ is determined
by PYTHIA [94]. The annihilation cross section as well as the photon flux is
calculated using micrOMEGAs where tabulated results of $dN_{i}/dE_{\gamma}$
from PYTHIA can be found. The results are compared to the 6 years of data from
the Fermi-LAT collaboration [95]. Notice that the photon flux is proportional
to $f^{2}_{\rm DM}$ which can help evade those constraints in the case of
multicomponent DM, i.e., for $f_{\rm DM}\ll 1$.

To determine the constraints on the DM thermally averaged annihilation cross
section from the Fermi-LAT measurements, we use the published data of 15 Milky
Way dwarf spheroidal galaxies (dSphs) from the Fermi-LAT collaboration [95].
The published six-year Pass 8 data pertain to the measured photon flux and the
bin-by-bin test statistic for each of the dSphs considered in the analysis.
Taking the LAT likelihood for target $i$ as
$\mathcal{L}_{i}(\bm{\mu},\bm{\theta}_{i}|\mathcal{D}_{i})$, where $\bm{\mu}$
contains the DM model parameters, $\bm{\theta}_{i}$ are the nuisance
parameters and $\mathcal{D}_{i}$ are the gamma ray data, we construct the
combined likelihood of 15 dSphs as

$\mathcal{L}(\bm{\mu},\bm{\theta}|\mathcal{D})=\prod_{i}\mathcal{L}_{i}(\bm{\mu},\bm{\theta}_{i}|\mathcal{D}_{i}).$
(3.13)

Then we define the test statistic

$q_{s}=-2\ln\left[\frac{\mathcal{L}(\bm{\mu_{0}},\bm{\hat{\theta}}|\mathcal{D})}{\mathcal{L}(\bm{\hat{\mu}},\bm{\hat{\theta}}|\mathcal{D})}\right],$
(3.14)

which determines the significance of the DM hypothesis. In our definition,
$\bm{\mu_{0}}$ represents the theory parameters under the null hypothesis
while the hatted variables are the best fit parameters under the DM
hypothesis. The upper limit on $\langle\sigma v\rangle$ is determined for
$q_{s}=2.71$ which represents a 90% quantile of a $\chi^{2}$ distribution. We
present in Fig. 1 the obtained upper limits for two cases: varying $g_{X}$
(left panel) and varying $m_{Z^{\prime}}$ (right panel). For the $D\bar{D}\to
Z^{\prime}Z^{\prime}$ channel, we notice a strong dependence on $g_{X}$ and
milder dependence on $m_{Z^{\prime}}$.

Figure 1: The 90% upper limits on the thermally averaged cross section using
the Fermi-LAT data [95] for $Z^{\prime}Z^{\prime}$ final state. We show the
limits for different $g_{X}$ values (left panel) and different
$m_{Z^{\prime}}$ (right panel) with the $1\sigma$ uncertainty band.

The right panel of Fig. 1 shows two sharp dips in the upper limits for
$m_{Z^{\prime}}=400$ GeV and 800 GeV. Notice the dips occur at
$m_{D}=m_{Z^{\prime}}/2$ which corresponds to the resonance region for the
process $D\bar{D}\to f\bar{f}$. In this case, the annihilation cross section
becomes very large which results in a large photon flux. Therefore, the
analysis of Fermi-LAT data produces very stringent bounds in this region and
so the upper limit on $\langle\sigma v\rangle$ becomes very small, i.e., more
stringent as can be seen from the dips. The limits we present in Fig. 1 will
be translated into constraints on the kinetic mixing and dark photon mass.

##  4 Exclusion limits

In this Section we give the exclusion limits on the Stueckelberg $Z^{\prime}$
model from the various experimental constraints discussed in the previous
section. The results pertain to the heavy and light $Z^{\prime}$ for a wide
range of kinetic mixing. Before we discuss the results, note that the
$m_{Z^{\prime}}$ in our model is not a free parameter as it depends on the
kinetic mixing $\delta$ and the mass parameter $M_{1}$. In the case of
vanishing mass mixing, recall that the $Z^{\prime}$ mass is
$m_{Z^{\prime}}^{2}=(q\pm p)/2$, where [3]

$\displaystyle p$
$\displaystyle=\sqrt{\left[M_{1}^{2}c^{2}_{\delta}+\frac{v^{2}}{4}\left(g_{Y}^{2}c^{2}_{\delta}+g_{2}^{2}\right)\right]^{2}-M_{1}^{2}v^{2}c^{2}_{\delta}(g_{Y}^{2}+g_{2}^{2})},$
(4.1) $\displaystyle q$
$\displaystyle=M_{1}^{2}c^{2}_{\delta}+\frac{v^{2}}{4}\left(g_{Y}^{2}c^{2}_{\delta}+g_{2}^{2}\right).$
(4.2)

The positive sign in $m^{2}_{Z^{\prime}}$ corresponds to a heavy $Z^{\prime}$
while the negative is for a light $Z^{\prime}$.

###  4.1 Heavy $Z^{\prime}$ bosons

We first consider the case of a heavy $Z^{\prime}$ with
$m_{D}>2m_{Z^{\prime}}$, whose decay width to SM fermions is given by

$\Gamma_{Z^{\prime}\rightarrow
f\bar{f}}=\frac{N_{c}\,g_{2}^{2}}{48\pi\,\cos^{2}\theta}m_{Z^{\prime}}\sqrt{1-\frac{4m_{f}^{2}}{m_{Z^{\prime}}^{2}}}\left[v_{f}^{\prime
2}+a_{f}^{\prime 2}+\frac{2m_{f}^{2}}{m_{Z^{\prime}}^{2}}(v_{f}^{\prime
2}-2a_{f}^{\prime 2})\right],$ (4.3)

where $N_{c}=3$ for quarks, 1 for leptons and $1/2$ for neutrinos. In Fig. 2
we show the relevant constraints in the kinetic mixing-$Z^{\prime}$ mass plane
which include dijet and dilepton searches from ATLAS and CMS as well as the
LEP constraint. We also show constraints from direct detection experiments
(Xenon1T and LZ), Fermi-LAT and the DM relic density. The latter only appears
in the left panel for reasons we discuss thereafter. The figure also exhibits
the 95% and 68% CL regions obtained from Contur using precision measurements
of SM processes. We identify in Fig. 12 of Appendix B the different analyses
pools giving the highest sensitivity for exclusion. In this figure we plot
$M_{1}$ rather than $m_{Z^{\prime}}$.

Figure 2: Exclusion limits from LHC analyses, LEP, DM relic density, direct
and indirect detection experiments in the kinetic mixing-$Z^{\prime}$ mass
plane for $m_{D}=2m_{Z^{\prime}}$ (left panel) and $m_{D}=1.5m_{Z^{\prime}}$
(right panel). The dashed black contour shows the projected reach of HL-LHC in
the dilepton channel. Note that the region bordered by a red solid (dashed)
line represents the 95% (68%) CL from Contur.

In plotting Fig. 2, we chose $m_{D}=2m_{Z^{\prime}}$ (left panel) and
$m_{D}=1.5m_{Z^{\prime}}$ (right panel) which means that $Z^{\prime}$ has no
dark decays and the process $D\bar{D}\to Z^{\prime}Z^{\prime}$ is
kinematically allowed. As a result, the DM relic density is set by the latter
annihilation process and so we expect the relic density to have mild to no
dependence on the kinetic mixing $\delta$. In fact, for a fixed
$m_{Z^{\prime}}$, the relic density is solely determined by $g_{X}$ which we
take here to be 1.0 in the left panel. The gold-colored area in the left panel
of Fig. 2 shows the region excluded by the relic density. One can see that for
a fixed $m_{Z^{\prime}}$ changing $\delta$ has no effect on the relic density.
But this begins to change for $\delta\gtrsim 0.3$ as the boundary starts
curving rightward for larger $m_{Z^{\prime}}$. The reason is that
$m_{Z^{\prime}}$ is not an independent parameter. According to Eqs. (4.1) and
(4.2), $m_{Z^{\prime}}$ depends on the scanning parameters $M_{1}$ and
$\delta$. For small $\delta$, $m_{Z^{\prime}}\approx M_{1}$, but for larger
values, $m_{Z^{\prime}}\gg M_{1}$. In this case, the kinetic mixing starts
affecting the relic density because it causes the $Z^{\prime}$ mass to change.
Larger $m_{Z^{\prime}}$ means that the process $D\bar{D}\to
Z^{\prime}Z^{\prime}$ becomes less efficient and so DM does not readily
depletes causing the relic density to shoot up. The relic density constraint
disappears from the right panel for the choice $m_{D}=1.5m_{Z^{\prime}}$ and
$g_{X}=2.0$ since in this case all points have a relic density smaller than
Eq. (3.6). As one can clearly see from Fig. 2, DM (in)direct detection and
dilepton searches as well as the DM relic density (for the left panel) are the
most constraining limits on the parameter space. However, there remains parts
of the parameter space that can still be explored as seen from the right
panel. We draw the projected reach in the kinetic mixing-mass plane at HL-LHC
in the dilepton channel (black dashed curve). We will explore in section 5 the
discovery potential of HL-LHC as a validity of this region drawn here. One
more comment is in order regarding the Contur limits from precision
measurements of SM processes. As seen from Fig. 2, those limits are very
competitive and are more stringent than the dijet limits and come close to the
reach of the dilepton limits near 1 TeV.

Figure 3: Limits from DM direct and indirect detection experiments, relic
density and LHC constraints from ATLAS dilepton, dijet and monojet searches at
139 fb-1 for $m_{D}=m_{Z^{\prime}}/4$ and $g_{X}=0.1$ (left panel) and
$g_{X}=1.0$ (right panel). The 95% CL and 68% CL contours from Contur and LEP
constraint are also shown.

We now allow $Z^{\prime}$ to decay to the dark fermions by taking
$m_{D}=m_{Z^{\prime}}/4$, thus opening a new decay channel with a decay width

$\Gamma_{Z^{\prime}\rightarrow
D\bar{D}}=\frac{m_{Z^{\prime}}}{12\pi}g_{X}^{2}(\mathcal{R}_{11}-s_{\delta}\mathcal{R}_{21})^{2}\sqrt{1-\frac{4m_{D}^{2}}{m_{Z^{\prime}}^{2}}}\left(1+\frac{2m_{D}^{2}}{m_{Z^{\prime}}^{2}}\right).$
(4.4)

This will weaken the limits from dilepton and dijet searches as one can
clearly see from Fig. 3, with the left panel corresponding to $g_{X}=0.1$ and
the right one to $g_{X}=1.0$. The LEP constraint as well as limits from
precision calculations of SM processes obtained from Contur are added along
with limits from direct and indirect detection experiments. Assuming thermal
production of DM, the final DM relic density for this setup is set by
annihilation to SM fermions since now $D\bar{D}\to Z^{\prime}Z^{\prime}$ is
not accessible. The processes $D\bar{D}\to f\bar{f}$ depend on
$\delta^{2}g_{X}^{2}$ and so the relic density can only be satisfied for large
enough $\delta$ as shown in Fig. 3. This region, however, is already excluded
by LEP and dilepton searches. Even for larger $g_{X}$, the yellow region opens
up but is still not enough to evade LHC constraints. One can thus see that
this parameter space is completely ruled out. But this only corresponds to the
DM-$Z^{\prime}$ mass relation used here. So it is important to check other
values by scanning over $m_{D}$ and $m_{Z^{\prime}}$ instead.

Figure 4: Constraints on the Stueckelberg model in the DM-$Z^{\prime}$ mass
plane for two choices of $\delta$ and $g_{X}$. Limits from direct detection
experiments, monojet and dilepton searches as well as the relic density
constraint eliminate most of the parameter space except near the resonance
region for smaller $\delta$ and $g_{X}$ (left panel).

Fig. 4 shows the different experimental constraints in the DM-$Z^{\prime}$
mass plane for a specific choice of $\delta$ and $g_{X}$. In addition to the
previous limits, monojet searches become relevant here and are also included.
In the left panel of Fig. 4 one can see that most of the parameter space is
excluded except for a narrow region along the diagonal which corresponds to
the resonant decay $Z^{\prime}\to D\bar{D}$. The relic density constraint
misses this region because the annihilation channel $D\bar{D}\to f\bar{f}$
proceeds on-shell which efficiently depletes the DM density. In the right
panel of Fig. 4 we exhibit the constraints for higher values of $\delta$ and
$g_{X}$. In this case the monojet and dilepton constraints become very strong
that they exclude all the parameter space including the resonance region.

The important takeaway here is that the Stueckelberg parameter space is
severely constrained for the case of dark decays of $Z^{\prime}$. This is due
to DM-related constraints, i.e., the relic density as well as direct and
indirect detection experiments. The only available window is near the
resonance region for smaller couplings. For heavier DM mass, no $Z^{\prime}$
dark decays occur and the process $D\bar{D}\to Z^{\prime}Z^{\prime}$ become
available which constitutes an important channel for depletion of the DM
abundance. For a large choice of $g_{X}$, the relic density constraint can be
severely weakened but it does not come with out a cost. The Fermi-LAT
constraint becomes stronger but would still allow a considerable part of the
parameter space to still be open. One can see here that the combination of all
these constraints points to a thermalized dark sector, i.e., a sector in which
its constituent species are in thermal equilibrium. In other words, the DM
relic density is controlled by annihilation within the dark sector itself
rather than annihilation into SM fermions.

###  4.2 Light $Z^{\prime}$ bosons

In this Section we present the experimental limits on our model for the light
$Z^{\prime}$ case, i.e., $m_{Z^{\prime}}<m_{Z}$. For this mass range, limits
from BaBar and LHCb become important. Since most of the limits are in the
$\delta<0.01$ region, one can safely take $M_{1}\approx m_{Z^{\prime}}$ so
that the vector and axial-vector couplings reduce to the simple forms

$\displaystyle v_{f}^{\prime}$
$\displaystyle\simeq\frac{\sin\theta_{W}}{1-(m_{Z^{\prime}}/m_{Z})^{2}}\left[2Q_{f}(\sin^{2}\theta_{W}-1)+\frac{m_{Z^{\prime}}^{2}}{m_{Z}^{2}}(2Q_{f}-T_{3f})\right]\delta,$
(4.5) $\displaystyle a_{f}^{\prime}$
$\displaystyle\simeq\frac{m_{Z^{\prime}}^{2}}{m_{Z^{\prime}}^{2}-m_{Z}^{2}}\sin\theta_{W}T_{3f}\delta.$
(4.6)

Also in the small $\delta$ limit,
$(\mathcal{R}_{11}-s_{\delta}\mathcal{R}_{21})\sim 1$ so that
$g_{Z^{\prime}}=g_{X}Q_{X}(\mathcal{R}_{11}-s_{\delta}\mathcal{R}_{21})\sim
g_{X}Q_{X}$, where we take $Q_{X}=1$. This approximation is implemented in
DarkCast to derive the experimental upper limits for our model in the kinetic
mixing-mass plane as shown in Fig. 5.

Figure 5: The different recasted experimental constraints from LEP, BaBar,
CMS, LHCb and (in)direct detection experiments for the light $Z^{\prime}$
case. We also show the 95% CL (solid line) and the 68% CL (dashed line)
excluded regions from Contur as well as the relic density constraint. Here
$m_{D}=2m_{Z^{\prime}}$ and so the limits correspond to the case of
$Z^{\prime}$ decaying to SM fermions only and no dark decays of the
$Z^{\prime}$.

One can see in Fig. 5 that for $m_{Z^{\prime}}>1$ GeV the parameter space is
excluded assuming that $g_{X}$ is large enough to produce a thermal dark
sector. Again here the channel $D\bar{D}\to Z^{\prime}Z^{\prime}$ is
responsible for setting the DM relic density. Now we allow the dark photon to
decay to DM by setting $m_{D}=m_{Z^{\prime}}/4$. In this case, DM annihilation
via $D\bar{D}\to f\bar{f}$ becomes important and since we are considering
small values of $\delta$ as shown in Fig. 6, thermal equilibrium between the
dark sector and the visible sector cannot be guaranteed. One can easily check
this by comparing $n_{D}\langle\sigma v\rangle_{D\bar{D}\to f\bar{f}}$ to the
Hubble parameter $H$, where $n_{D}\langle\sigma v\rangle_{D\bar{D}\to
f\bar{f}}<H$ for all temperatures means the two sectors have not reached
thermal equilibrium666Note that one should also compare $n_{f}^{\rm
eq}\langle\sigma v\rangle_{f\bar{f}\to Z^{\prime}}$ to $H(T)$.. To determine
the DM relic density, we numerically solve the coupled Boltzmann equations,
Eqs. (3.4) and (3.8), assuming the freeze-in mechanism. In other words, owing
to the small value of the kinetic mixing, we assume that DM has a negligible
initial abundance in the early universe and that this abundance gradually
increases due to annihilation of SM particles, i.e., $f\bar{f}\to D\bar{D}$
and $f\bar{f}\to Z^{\prime}$. In Fig. 6, we show in light blue the relic
density constraint for two values of $g_{X}$ along with the numerous limits
from other experiments. We notice that the available parameter space grows
with increasing $g_{X}$ due to increased DM depletion.

Figure 6: Constraints from various visible and invisible limits for a
Stueckelberg dark photon in the mass range less than the $Z$ boson mass. The
left panel corresponds to $g_{X}=0.1$ and the right for $g_{X}=1.0$ and with
the common assumption that $m_{D}=m_{Z^{\prime}}/4$. The light blue region
corresponds to the relic density constraint..

Note how the limits change as we increase the value of $g_{X}$ going from the
left panel to the right. Since the main experiments look into the decays of
$Z^{\prime}$ (visible and invisible decays), it is important to properly model
such decays especially that in this mass region, hadronic decays of the dark
photon can no longer be simply the sum of the decays to $q\bar{q}$. Here we
have used DarkCast to accurately determine the dark photon branching ratios.
This is shown in Fig. 7 for the cases of no dark decays (left panel) and with
dark decays for $g_{X}=0.1$ (middle panel) and $g_{X}=1.0$ (right panel).

Figure 7: Branching fractions of the Stueckelberg dark photon calculated using
DarkCast. In the figure legends, $\nu\nu$ represents the sum of all three
generations of neutrinos, ‘Hadrons’ represents the sum of all the branching
ratios of the hadronic channels. The left panel shows the branching fractions
when there’s no dark decays while the middle and right panels show the case of
dark decays for $g_{X}=0.1$ and $g_{X}=1.0$, respectively. Since we take
$m_{D}=m_{Z^{\prime}}/4$, the hidden sector decays are present for any
$Z^{\prime}$ mass. The analysis shows the strong effect of $g_{X}$ on the
branching ratio in the hidden sector from minimal at $g_{X}=0.1$ to dominant
at $g_{X}=1.0$.

Overall, we arrive here at the same conclusion we drew in the heavy
$Z^{\prime}$ case. For a kinetically mixed dark photon with dark decays, the
parameter space is severely constrained while more available parameter space
remains for the case of DM heavier than the mediator $Z^{\prime}$.

##  5 Discovery potential of the Stueckelberg $Z^{\prime}$ boson at the HL-LHC

The $Z^{\prime}$ mass reach projected at HL-LHC and exhibited in Fig. 2 shows
that a 1 TeV $Z^{\prime}$ with $\delta\sim\mathcal{O}(10^{-2})$ and heavier
can be probed. In this section, we perform a detailed analysis for potential
discovery of a Stueckelberg $Z^{\prime}$ with a TeV scale mass at HL-LHC. Here
we focus on the dilepton channel which, despite having a lower cross section
than the dijet channel, is clean and one can reconstruct the dilepton
invariant mass with much less SM background. However, the challenge here is
the naturally small production cross section due to a TeV mass $Z^{\prime}$
and a small kinetic mixing.

The production cross section of a Stueckelberg $Z^{\prime}$ is evaluated at
NLO. To do so, we implement the model in FeynRules [96] interfaced with NLOCT
[97] and FeynArts [98]. The obtained UFO files are used in MadGraph5_aMC@NLO
to determine the LO and NLO cross section of $pp\to
Z^{\prime}\to\ell^{+}\ell^{-}$ (see ref. [99] for a model-independent analysis
of $W^{\prime}$ and $Z^{\prime}$ production at the LHC). The LO and NLO cross
section of the process as a function of $m_{Z^{\prime}}$ is shown in Fig. 8
for three values of $\delta$. The $K$ factor defined as $K=\sigma_{\rm
NLO}/\sigma_{\rm LO}$ is shown in the bottom panel. We notice a factor of
$\sim 1.5$ increase from the LO prediction for $m_{Z^{\prime}}\sim 2$ TeV.

Figure 8: The LHC production cross section of a Stueckelberg $Z^{\prime}$ as a
function of $m_{Z^{\prime}}$ in the dilepton channel at $\sqrt{s}=14$ TeV. The
cross section is calculated at LO and NLO for three values of the kinetic
mixing $\delta$. The $K$ factor is shown in the bottom plot.

The signal Monte Carlo events are generated at LO using MadGraph5 and the
cross section is scaled accordingly using the obtained $K$ factors. The
dilepton final state has several sources of SM backgrounds. The dominant ones
are: diboson (mainly $WW$), $Z/\gamma^{*}$+jets, dilepton from off-shell
vector boson decay, $t\bar{t}$, single top and top associated production with
a vector boson. Background and signal events are generated with MadGraph5 and
showered with PYTHIA8 [100, 101] (adding ISR and FSR jets and ignoring
multiparticle interactions). Detector effects are included using Delphes [102]
which implements FastJet [103] for jet clustering with the anti-$k_{t}$ [104]
algorithm and jet radius $R=0.4$.

The event preselection is based on a recent analysis by the ATLAS
collaboration [56]. Electrons with a transverse energy of $E_{T}>30$ GeV and
located within $|\eta|<2.47$ are selected, while muons with a transverse
momentum $p_{T}>30$ GeV and $|\eta|<2.5$ are kept. Events are required to
contain at least two same flavor leptons. Candidate events with two muons are
required to have oppositely charged muons while such a requirement is not
forced on electrons because high $E_{T}$ electrons suffer from a higher
probability of charge misidentification. If an event contains more than two
leptons, then the electrons (muons) with the highest $E_{T}$ ($p_{T}$) are
kept. If an event is found to contain two lepton pairs, then the electron pair
is retained because the ATLAS detector has a better resolution and higher
efficiency for electrons. Next, the dilepton invariant mass is reconstructed
and a minimum cut of 220 GeV is applied as means to reject the overwhelming SM
events near the $Z$ pole mass.

###  5.1 Cut-and-count analysis

We select two benchmarks $(m_{Z^{\prime}},\delta)$ which lie within the region
of reach for HL-LHC. The benchmarks $(1$ TeV$,10^{-2})$ and $(2$ TeV$,3\times
10^{-2})$ have NLO cross sections of 0.241 fb and 0.097 fb, respectively. For
the cut-and-count analysis, we employ the kinematic variables

$m_{\ell\ell},~{}~{}~{}E_{T1},~{}~{}~{}E_{T2},~{}~{}~{}p_{T1},~{}~{}~{}p_{T2},$
(5.1)

where the subscripts ‘1’ and ‘2’ indicate leading and subleading leptons,
respectively. We perform a cut-and-count analysis where different cuts based
on the above kinematic variables are implemented with the aim to maximize the
$S/\sqrt{S+B}$ figure of merit. Using 3000 fb-1 as the maximum integrated
luminosity projected at HL-LHC, the figure of merit never reaches the
$5\sigma$ limit required for discovery. The main culprit here is the
irreducible SM dilepton background from off-shell decays of a vector boson.
This is clear from Fig. 9 where we show the signal and background distribution
in the invariant dilepton mass.

Figure 9: The dilepton invariant mass distribution for the signal and
background samples. The analysis is for the case of two benchmarks discussed
in the text, i.e., for $(m_{Z^{\prime}},\delta)$ cases $(1$ TeV, $10^{-2}$)
and $(2$ TeV, $3\times 10^{-2})$ using the cut-and-count method with cuts on
the variables of Eq. (5.1). This method is found not efficient for discovery
with a projected luminosity of 3000 fb-1 at HL-LHC. A better technique is
discussed in the next section.

###  5.2 Boosted decision tree analysis

To have a better discrimination between the SM background and the signal, we
use a boosted decision tree (BDT) which is part of the TMVA (Toolkit for
Multivariate Analysis) [105, 106] framework embedded in ROOT [107, 108]. We
train a BDT on the signal and background events using the above kinematic
variables. The training phase is following by a testing phase carried out on
statistically independent Monte Carlo samples of the signal and background
events, where the algorithm determines a new kinematic variable called the
‘BDT response’. This variable is a powerful discriminant necessary to enhance
$S/\sqrt{S+B}$. To include the effect of uncertainties, we actually use

$\frac{S}{\sqrt{S+B+(\delta_{S}S)^{2}+(\delta_{B}B)^{2}}}$ (5.2)

as the figure of merit for a $5\sigma$ discovery. Here, $\delta_{S}$
($\delta_{B}$) represents the systematic uncertainty in the signal
(background) which we take to be 10% (20%).

We show in Fig. 10 the distribution of signal and background events in the new
BDT variable for the two benchmarks of choice. The lower panels indicate the
effect of cuts on the significance defined in Eq. (5.2). One can see from the
left panel that a cut on the BDT response $>0.3$ produces a $5\sigma$
significance for an integrated luminosity of 2200 fb-1, while from the right
panel a cut $>0.4$ is required for discovery at an integrated luminosity of
2050 fb-1. Note here that additional cuts are required to arrive at the
desired results. Along with the cut on the BDT response, we require
$m_{\ell\ell}>500$ GeV and $E_{T1}>350$ GeV (left panel) and
$m_{\ell\ell}>1400$ GeV and $p_{T2}>350$ GeV (right panel).

Figure 10: Signal and background event distribution in the ‘BDT response’
variable obtained after training and testing the BDT on the background and
signal samples. The model points are the same as in Fig. 5.1. The bottom plots
show the significance given by Eq. (5.2) as a result of applying cuts on the
‘BDT response’. Here, both model points are discoverable at the HL-LHC.

##  6 Detection of a Stueckelberg dark photon at the Forward Physics Facility

As the kinetic mixing coefficient takes on smaller values, a dark photon
becomes a long-lived particle (LLP) which can still decay to the SM inside a
detector after having traveled a certain distance away from its production
vertex. In the sub-GeV regime and aside from direct production (for e.g.
Bremsstrahlung [109] and Drell-Yan production [110]), a dark photon can be
produced from the decay of a SM particle such as a meson. Light mesons are
copiously produced in the forward region at the LHC and therefore constitute
an important tool to study BSM physics with dark photons as leading
candidates. The Forward Physics Facility (FPF) [111, 112] is intended to host
several experiments which are able to capture the multitude of particles near
the beam line which are otherwise missed by the current LHC experiments. So
the FPF will be suited for studying and possibly detecting such particles. In
this Section we discuss the sensitivity reach of forward detectors at HL-LHC
and future colliders in terms of discovering a Stueckelberg dark photon.

For our analysis, we consider the two production modes of a dark photon:
direct production and production via meson decays. The latter requires a good
understanding of the meson spectra, which has been studied and greatly
improved over the years [113]. A new numerical package called FORESEE [114]
(FORward Experiment SEnsitivity Estimator) allows users to implement their
model and derive predictions on the sensitivity reach at future forward
detectors. The package also provides the meson spectra which is necessary to
determine the LLP flux generated from the decay of mesons. We implement our
Stueckelberg dark photon model in FORESEE and provide the long-lived dark
photon lifetime, its production rates and its decay branching ratios with the
latter estimated using DarkCast (see Fig. 7). After taking into account the
detector geometry and acceptance cuts, the number of surviving signal events
are counted and used in FORESEE to draw the contours reflecting the
sensitivity reach at forward detectors. We consider in this analysis the mass
reach at FASER [115, 116, 117], which is already installed at the LHC beam
line, FASER 2 which is planned for HL-LHC [118], as well as possible future
detectors at HE-LHC [119] and at the Future Circular Collider (FCC) [120]. The
corresponding limits are shown in Fig. 11.

Figure 11: The predicted sensitivity reach as determined by FORESEE for future
forward detectors: FASER (orange line), FASER2 (red line), HE-LHC (brown line)
and FCC (purple). Other experimental constraints are shown including CHARM
[121], $\nu$Cal [121, 122, 123], E137 [124], E141 [72], NA64 [125], NA48
[126], BaBar [66, 67], HPS [127], LHCb [63], Belle-2 [128], SHiP [129],
SeaQuest [130, 121] and NA62 [121]. The pink band represents the relic density
constraint consistent with Eq. (3.6). To ensure the dark photon will not decay
to dark fermions, we set $m_{D}=0.6m_{Z^{\prime}}$.

The projected limits are mainly derived from the dark photon decay channel
$Z^{\prime}\to e^{+}e^{-}$ whose branching ratio becomes progressively smaller
for larger $Z^{\prime}$ mass as hadronic decays become more favorable. This
means one can no longer reach a dark photon mass larger than $\sim 0.5$ GeV
when investigating dilepton final states. As one can see that there is a major
gain in the reach along the kinetic mixing axis when going from FASER to the
other future detectors. However, the gain is modest if one compares FASER2,
HE-LHC and FCC. The same is true along the mass axis, where the future
detectors can probe larger masses and major improvement from FASER is
observed. The pink region is the part of the parameter space where the DM
relic density is satisfied, i.e., consistent with Eq. (3.6). Several
phenomenological work has been put forth regarding light $Z^{\prime}$ at FASER
[131, 132, 133] and heavy $Z^{\prime}$ at future colliders [134].

##  7 Conclusion

In this work we have investigated the current constraints from collider
experiments as well as from DM direct and indirect detection experiments on a
well motivated extension of the electroweak sector of the Standard Model with
a $U(1)_{X}$ gauge group. The $U(1)_{X}$ belongs to the hidden sector which
also contains matter, that we assume to be a Dirac fermion field. Although the
hidden sector is neutral relative to the Standard Model gauge group, it can
still communicate with the visible sector via kinetic mixing and via
Stueckelberg mass mixing. We discussed the kinetic energy in the mass diagonal
basis, where one finds a massive dark photon or dark $Z^{\prime}$ which can
interact with the quarks and leptons in the visible sector with mixings
characterized by the kinetic and mass mixing parameters. However, the dark
photon has regular size couplings with the dark fermion allowing for a dark
freeze-out to happen in the hidden sector, which generates the desired relic
density for the dark fermions. In the analysis we have recast all the relevant
constraints for a wide range of $Z^{\prime}$ masses: above and below the $Z$
boson pole mass. For each case we considered both visible and dark decays of
$Z^{\prime}$, where our analysis shows that the model parameter space is
severely constrained for the case when a $Z^{\prime}$ decays to DM. The reason
is that the region with a large kinetic mixing is already excluded by the LHC
and LEP and this region is important to deplete the DM abundance via
$D\bar{D}\to f\bar{f}$ when the latter process is the only active one near
freeze-out. However, for heavier DM masses, i.e., for $m_{Z^{\prime}}>2m_{D}$,
the dark photon can no longer decay to DM, and the process $D\bar{D}\to
Z^{\prime}Z^{\prime}$ becomes the dominant channel for DM annihilation. This
process is proportional to $g_{X}^{4}$ and so a large enough $g_{X}$ can
weaken the relic density constraint thus opening up the available parameter
space. Note that an increase in $g_{X}$ leads to more severe constraints from
Xenon1T, LZ and Fermi-LAT. However, this is compensated by a small $f_{\rm
DM}$ and a competition between $f_{\rm DM}$ and $g_{X}$ can go either way. In
our analysis we see a slight increase in these constraints, but they are well
tamed in such a way that parts of the model parameter space remain viable. We
have also shown that unexplored regions of the parameter space can be
accessible at HL-LHC and in forward detectors at the LHC and at future
colliders. For the HL-LHC, we have carried out a detailed analysis for a
potential discovery of a TeV mass scale $Z^{\prime}$ and showed that one can
observe a $5\sigma$ excess for a 1 TeV or 2 TeV $Z^{\prime}$ using BDTs.
Furthermore, we demonstrated the sensitivity reach of forward detectors for a
sub-GeV dark photon at the LHC (FASER) and at future colliders (HE-LHC and
FCC).

Our analysis shows that the presence of a thermal hidden sector can weaken the
current constraints on dark $Z^{\prime}$ models while the available regions of
the model parameter space remain within reach of standard LHC searches as well
as at forward detectors.

Acknowledgments: The research of AA and MK was supported by the BMBF under
contract 05P21PMCAA and by the DFG through the Research Training Network 2149
“Strong and Weak Interactions - from Hadrons to Dark matter” and grant KL
1266/10-1, while the research of PN and ZYW was supported in part by the NSF
Grant PHY-2209903. The research of MMA was supported by the French Agence
Nationale de la Recherche (ANR) grant no. ANR-21-CE31-0023 (PRCI SLDNP) and by
the National Science Centre, Poland, under research grant 2017/26/E/ST2/00135.
MK thanks the School of Physics at the University of New South Wales in
Sydney, Australia for its hospitality and financial support through the Gordon
Godfrey visitors program.

##  Appendix A Rotation angles and dark matter couplings

Further details of the analysis presented in the main body of the paper are
given in this Appendix. As noted in section 2, we have mixing of three vector
bosons: $C^{\mu},B^{\mu},A^{\mu}_{3}$ which leads to diagonalization of both
the kinetic energy matrix and the mass matrix involving the fields. In general
this results in a diagonalization of a $3\times 3$ vector boson mass-square
matrix ${\cal M}^{2}$ which, however, is symmetric and can be diagonalized by
an orthogonal transformation ${\cal R}$ so that

$\displaystyle\mathcal{R}^{T}\mathcal{M}^{2}\mathcal{R}=\text{diag}(m^{2}_{Z^{\prime}},m_{Z}^{2},0),$
(A.1)

where $\mathcal{M}^{2}$ is defined by Eq. (21) of [3]. The three Euler angles
$\theta,\phi,\psi$ are given by

$\tan\phi=-\sinh{\delta},~{}~{}~{}\tan\theta=\frac{g_{Y}}{g_{2}}\cosh{\delta}\cos\phi,~{}\tan
2\psi=\frac{2m^{2}_{Z}\sin\theta\tan\phi}{m^{2}_{Z^{\prime}}-m^{2}_{Z}+(m^{2}_{Z^{\prime}}+m^{2}_{Z}-m^{2}_{W})\tan^{2}\phi},$
(A.2)

with the assumption of zero mass mixing. Thus, the couplings
$g_{Z},g_{\gamma}$ and $g_{Z^{\prime}}$ that appear in $\Delta\mathcal{L}_{\rm
int}$ (see Eq. (2.3)) are given by

$g_{Z}=g_{X}Q_{X}(\mathcal{R}_{12}-s_{\delta}\mathcal{R}_{22}),~{}g_{\gamma}=g_{X}Q_{X}(\mathcal{R}_{13}-s_{\delta}\mathcal{R}_{23}),~{}g_{Z^{\prime}}=g_{X}Q_{X}(\mathcal{R}_{11}-s_{\delta}\mathcal{R}_{21}).$
(A.3)

##  Appendix B More on the exclusion plots from Contur

In this appendix, we give the exclusion plots from precision measurements of
the SM obtained using Contur. The plots are drawn for the parameters $\delta$
and $M_{1}$ over which the scan is made. The limits in these plots are then
converted to the kinetic mixing-$m_{Z^{\prime}}$ plane using
$m_{Z^{\prime}}^{2}=(q\pm p)/2$ and Eqs. (4.1) and (4.2). Let us begin by
explaining the different data pools used by Contur.

An event that passes the cuts of a specific measurement can also be accepted
in measurements that share similar final states. In order to avoid multiple
counting of such events, and due to the lack of information about the
correlations between different measurements, the analyses in Rivet are grouped
into orthogonal pools based on three criteria: the experiment that conducted
the measurement, the center of mass energy, and the considered final state.
For each pool of analyses, a likelihood is built for every distribution taking
the correlation between its bins777Note that this can only be done if the
correlation information is provided by the experiment. If not, Contur will
only consider the most sensitive bin of the histogram. into account. The
likelihoods of the orthogonal histograms within the pool are then combined.
Finally, Contur constructs the total likelihood by combining the likelihoods
of the different pools.

The plots in this Appendix show the most sensitive pool at each point of the
parameter space for the different scenarios that we consider in this paper,
i.e. heavy $Z^{\prime}$ bosons decaying to SM fermions only (Fig. 12) or also
DM (Fig. 13) as well as light $Z^{\prime}$ bosons decaying to SM fermions only
(Fig. 14) or also DM (Fig. 15).

Figure 12: The breakdown of Contur’s most sensitive analysis pool for each
scan point in the heavy $Z^{\prime}$ case. The solid (dashed) white line
corresponds to the 95% (68%) CL exclusion on the kinetic mixing $\delta$
versus the $M_{1}$ parameter. Here $m_{D}=2m_{Z^{\prime}}$ and so the limits
correspond to the case of $Z^{\prime}$ decaying to SM fermions only.

ATLAS $E_{T}^{\rm miss}$+jet |  CMS high-mass Drell-Yan $\ell\ell$ |  ATLAS $ee$+jet
---|---|---
ATLAS $\ell\ell$+jet |  ATLAS high-mass Drell-Yan $\ell\ell$ |  ATLAS $\ell\ell\gamma$
CMS $\mu\mu$+jet |  ATLAS $\mu\mu$+jet |  ATLAS jets

Figure 13: The 95% CL (solid) and the 68% CL (dashed line) exclusions on the
kinetic mixing $\delta$ versus the $M_{1}$ parameter for the case of heavy
$Z^{\prime}$ invisible decay to DM fermions and for different values of
$g_{X}$. The left panel corresponds to $g_{X}=0.1$ while the right panel
corresponds to $g_{X}=1.0$. We take $m_{D}=m_{Z^{\prime}}/4$. The colored bins
correspond to the pool of analyses giving the dominant exclusion.

ATLAS $E_{T}^{\rm miss}$+jet |  ATLAS high-mass Drell-Yan $\ell\ell$ |  ATLAS $\ell\ell$+jet
---|---|---
ATLAS $\mu\mu$+jet |  ATLAS jets |  ATLAS $\ell\ell\gamma$
CMS $\mu\mu$+jet |  CMS jets |  CMS high-mass Drell-Yan $\ell\ell$
CMS $\ell$+$E_{T}^{\rm miss}$+jet |  |

Figure 14: The 95% CL (solid) and the 68% CL (dashed line) exclusions on the
kinetic mixing $\delta$ versus the $M_{1}$ parameter . Here
$m_{D}=2m_{Z^{\prime}}$ and so the limits correspond to the case of light
$Z^{\prime}$ decaying to SM fermions only. The scan on $\delta$ is done in the
range $10^{-5}$ to $1$ and the colored bins correspond to the analyses pool
giving the dominant contribution.

ATLAS $\gamma$ |  ATLAS low-mass Drell-Yan $\ell\ell$ |  ATLAS $ee$+jet
---|---|---
ATLAS $\ell_{1}\ell_{2}$+$E_{T}^{\rm miss}$ |  CMS high-mass Drell-Yan $\ell\ell$ |  ATLAS $e$+$E_{T}^{\rm miss}$+jet
ATLAS $E_{T}^{\rm miss}$+jet |  CMS $\ell$+$E_{T}^{\rm miss}$+jet |  CMS jets
ATLAS jets |  ATLAS $\ell\ell$+jet |  ATLAS $\mu\mu$+jet
ATLAS $\ell_{1}\ell_{2}$+$E_{T}^{\rm miss}$+jet |  LHCb $\ell$+jet |  ATLAS $\mu$+$E_{T}^{\rm miss}$+jet

Figure 15: The 95% CL (solid) and the 68% CL (dashed line) exclusions on the
kinetic mixing $\delta$ versus the $M_{1}$ parameter for the case of light
$Z^{\prime}$ invisible decay to DM fermions with $m_{D}=m_{Z^{\prime}}/4$ and
for $g_{X}=0.1$ (left panel) and $g_{X}=1.0$ (right panel). The color-shading
scheme specifies which SM measurement has the dominant exclusion.

ATLAS $\mu\mu$+jet |  CMS $\ell$+$E_{T}^{\rm miss}$+jet |  ATLAS $\gamma$
---|---|---
ATLAS $ee$+jet |  ATLAS $E_{T}^{\rm miss}$+jet |  ATLAS jets
ATLAS $\ell_{1}\ell_{2}$+$E_{T}^{\rm miss}$ |  ATLAS high-mass Drell-Yan $\ell\ell$ |  CMS $\ell\ell$+jet
ATLAS $\ell_{1}\ell_{2}$+$E_{T}^{\rm miss}$+jet |  ATLAS $\mu$+$E_{T}^{\rm miss}$+jet |  ATLAS $\ell\ell$+jet
ATLAS low-mass Drell-Yan $\ell\ell$ |  LHCb $\ell$+jet |  CMS jets
CMS high-mass Drell-Yan $\ell\ell$ |  ATLAS $e$+$E_{T}^{\rm miss}$+jet |  ATLAS $\ell\ell\gamma$

## References

  * [1]
  * [2] B. Holdom, Phys. Lett. B 166, 196-198 (1986) doi:10.1016/0370-2693(86)91377-8
  * [3] D. Feldman, Z. Liu and P. Nath, Phys. Rev. D 75, 115001 (2007) doi:10.1103/PhysRevD.75.115001 [arXiv:hep-ph/0702123 [hep-ph]].
  * [4] A. Aboubrahim, W. Z. Feng, P. Nath and Z. Y. Wang, Phys. Rev. D 103, no.7, 075014 (2021) doi:10.1103/PhysRevD.103.075014 [arXiv:2008.00529 [hep-ph]].
  * [5] A. Aboubrahim, T. Ibrahim, M. Klasen and P. Nath, Eur. Phys. J. C 81, no.8, 680 (2021) doi:10.1140/epjc/s10052-021-09483-0 [arXiv:2012.10795 [hep-ph]].
  * [6] A. Aboubrahim, W. Z. Feng, P. Nath and Z. Y. Wang, JHEP 06, 086 (2021) doi:10.1007/JHEP06(2021)086 [arXiv:2103.15769 [hep-ph]].
  * [7] A. Aboubrahim, P. Nath and Z. Y. Wang, JHEP 12, 148 (2021) doi:10.1007/JHEP12(2021)148 [arXiv:2108.05819 [hep-ph]].
  * [8] P. Langacker, Rev. Mod. Phys. 81, 1199-1228 (2009) doi:10.1103/RevModPhys.81.1199 [arXiv:0801.1345 [hep-ph]].
  * [9] M. Fabbrichesi, E. Gabrielli and G. Lanfranchi, doi:10.1007/978-3-030-62519-1 [arXiv:2005.01515 [hep-ph]].
  * [10] B. Fuks, M. Klasen, F. Ledroit, Q. Li and J. Morel, Nucl. Phys. B 797, 322-339 (2008) doi:10.1016/j.nuclphysb.2008.01.017 [arXiv:0711.0749 [hep-ph]].
  * [11] M. Klasen, F. Lyonnet and F. S. Queiroz, Eur. Phys. J. C 77, no.5, 348 (2017) doi:10.1140/epjc/s10052-017-4904-8 [arXiv:1607.06468 [hep-ph]].
  * [12] R. Bonciani, T. Jezo, M. Klasen, F. Lyonnet and I. Schienbein, JHEP 02, 141 (2016) doi:10.1007/JHEP02(2016)141 [arXiv:1511.08185 [hep-ph]].
  * [13] M. M. Altakach, T. Ježo, M. Klasen, J. N. Lang and I. Schienbein, Phys. Rev. D 103, no.11, 115026 (2021) doi:10.1103/PhysRevD.103.115026 [arXiv:2012.14855 [hep-ph]].
  * [14] T. Jezo, M. Klasen, D. R. Lamprea, F. Lyonnet and I. Schienbein, JHEP 12, 092 (2014) doi:10.1007/JHEP12(2014)092 [arXiv:1410.4692 [hep-ph]].
  * [15] B. Kors and P. Nath, Phys. Lett. B 586, 366-372 (2004) doi:10.1016/j.physletb.2004.02.051 [arXiv:hep-ph/0402047 [hep-ph]].
  * [16] K. Cheung and T. C. Yuan, JHEP 03, 120 (2007) doi:10.1088/1126-6708/2007/03/120 [arXiv:hep-ph/0701107 [hep-ph]].
  * [17] D. Feldman, B. Kors and P. Nath, Phys. Rev. D 75, 023503 (2007) doi:10.1103/PhysRevD.75.023503 [arXiv:hep-ph/0610133 [hep-ph]].
  * [18] A. Aboubrahim, W. Z. Feng and P. Nath, JHEP 02, 118 (2020) doi:10.1007/JHEP02(2020)118 [arXiv:1910.14092 [hep-ph]].
  * [19] A. Aboubrahim and P. Nath, JHEP 09, 084 (2022) doi:10.1007/JHEP09(2022)084 [arXiv:2205.07316 [hep-ph]].
  * [20] M. Du, Z. Liu and P. Nath, Phys. Lett. B 834, 137454 (2022) doi:10.1016/j.physletb.2022.137454 [arXiv:2204.09024 [hep-ph]].
  * [21] F. Staub, Comput. Phys. Commun. 185, 1773-1790 (2014) doi:10.1016/j.cpc.2014.02.018 [arXiv:1309.7223 [hep-ph]].
  * [22] F. Staub, Adv. High Energy Phys. 2015, 840780 (2015) doi:10.1155/2015/840780 [arXiv:1503.04200 [hep-ph]].
  * [23] W. Porod, Comput. Phys. Commun. 153, 275-315 (2003) doi:10.1016/S0010-4655(03)00222-4 [arXiv:hep-ph/0301101 [hep-ph]].
  * [24] W. Porod and F. Staub, Comput. Phys. Commun. 183, 2458-2469 (2012) doi:10.1016/j.cpc.2012.05.021 [arXiv:1104.1573 [hep-ph]].
  * [25] A. Pukhov, [arXiv:hep-ph/0412191 [hep-ph]].
  * [26] E. E. Boos, M. N. Dubinin, V. A. Ilyin, A. E. Pukhov and V. I. Savrin, [arXiv:hep-ph/9503280 [hep-ph]].
  * [27] G. Bélanger, F. Boudjema, A. Pukhov and A. Semenov, Comput. Phys. Commun. 192, 322-329 (2015) doi:10.1016/j.cpc.2015.03.003 [arXiv:1407.6129 [hep-ph]].
  * [28] C. Degrande, C. Duhr, B. Fuks, D. Grellscheid, O. Mattelaer and T. Reiter, Comput. Phys. Commun. 183, 1201-1214 (2012) doi:10.1016/j.cpc.2012.01.022 [arXiv:1108.2040 [hep-ph]].
  * [29] J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, O. Mattelaer, H. S. Shao, T. Stelzer, P. Torrielli and M. Zaro, JHEP 07, 079 (2014) doi:10.1007/JHEP07(2014)079 [arXiv:1405.0301 [hep-ph]].
  * [30] J. M. Butterworth, D. Grellscheid, M. Krämer, B. Sarrazin and D. Yallup, JHEP 03, 078 (2017) doi:10.1007/JHEP03(2017)078 [arXiv:1606.05296 [hep-ph]].
  * [31] M. M. Altakach, [arXiv:2104.10608 [hep-ph]].
  * [32] C. Bierlich, A. Buckley, J. Butterworth, C. H. Christensen, L. Corpe, D. Grellscheid, J. F. Grosse-Oetringhaus, C. Gutschow, P. Karczmarczyk and J. Klein, et al. SciPost Phys. 8, 026 (2020) doi:10.21468/SciPostPhys.8.2.026 [arXiv:1912.05451 [hep-ph]].
  * [33] T. Junk, Nucl. Instrum. Meth. A 434, 435-443 (1999) doi:10.1016/S0168-9002(99)00498-2 [arXiv:hep-ex/9902006 [hep-ex]].
  * [34] A. L. Read, J. Phys. G 28, 2693-2704 (2002) doi:10.1088/0954-3899/28/10/313
  * [35] A. Buckley, J. M. Butterworth, L. Corpe, M. Habedank, D. Huang, D. Yallup, M. Altakach, G. Bassman, I. Lagwankar and J. Rocamonde, et al. SciPost Phys. Core 4, 013 (2021) doi:10.21468/SciPostPhysCore.4.2.013 [arXiv:2102.04377 [hep-ph]].
  * [36] J. Butterworth, J. Heeck, S. H. Jeon, O. Mattelaer and R. Ruiz, [arXiv:2210.13496 [hep-ph]].
  * [37] M. M. Altakach, J. M. Butterworth, T. Ježo, M. Klasen and I. Schienbein, [arXiv:2111.15406 [hep-ph]].
  * [38] M. M. Altakach, J. M. Butterworth, T. Ježo, M. Klasen and I. Schienbein, [arXiv:2204.10577 [hep-ph]].
  * [39] J. M. Butterworth, M. Habedank, P. Pani and A. Vaitkus, SciPost Phys. Core 4, 003 (2021) doi:10.21468/SciPostPhysCore.4.1.003 [arXiv:2009.02220 [hep-ph]].
  * [40] A. Buckley, J. M. Butterworth, L. Corpe, D. Huang and P. Sun, SciPost Phys. 9, no.5, 069 (2020) doi:10.21468/SciPostPhys.9.5.069 [arXiv:2006.07172 [hep-ph]].
  * [41] J. Bellm, G. Bewick, S. Ferrario Ravasio, S. Gieseke, D. Grellscheid, P. Kirchgaeßer, M. R. Masouminia, G. Nail, A. Papaefstathiou and S. Platzer, et al. Eur. Phys. J. C 80, no.5, 452 (2020) doi:10.1140/epjc/s10052-020-8011-x [arXiv:1912.06509 [hep-ph]].
  * [42] S. Schael et al. [ALEPH, DELPHI, L3, OPAL, SLD, LEP Electroweak Working Group, SLD Electroweak Group and SLD Heavy Flavour Group], Phys. Rept. 427, 257-454 (2006) doi:10.1016/j.physrep.2005.12.006 [arXiv:hep-ex/0509008 [hep-ex]].
  * [43] t. Electroweak [LEP, ALEPH, DELPHI, L3, OPAL, LEP Electroweak Working Group, SLD Electroweak Group and SLD Heavy Flavor Group], [arXiv:hep-ex/0312023 [hep-ex]].
  * [44] M. Tanabashi et al. [Particle Data Group], Phys. Rev. D 98, no.3, 030001 (2018) doi:10.1103/PhysRevD.98.030001
  * [45] A. M. Sirunyan et al. [CMS], JHEP 05, 033 (2020) doi:10.1007/JHEP05(2020)033 [arXiv:1911.03947 [hep-ex]].
  * [46] G. Aad et al. [ATLAS], JHEP 03, 145 (2020) doi:10.1007/JHEP03(2020)145 [arXiv:1910.08447 [hep-ex]].
  * [47] M. Aaboud et al. [ATLAS], Phys. Rev. Lett. 121, no.8, 081801 (2018) doi:10.1103/PhysRevLett.121.081801 [arXiv:1804.03496 [hep-ex]].
  * [48] T. Aaltonen et al. [CDF], Phys. Rev. D 79, 112002 (2009) doi:10.1103/PhysRevD.79.112002 [arXiv:0812.4036 [hep-ex]].
  * [49] M. Aaboud et al. [ATLAS], Phys. Lett. B 788, 316-335 (2019) doi:10.1016/j.physletb.2018.09.062 [arXiv:1801.08769 [hep-ex]].
  * [50] A. M. Sirunyan et al. [CMS], Phys. Rev. D 100, no.11, 112007 (2019) doi:10.1103/PhysRevD.100.112007 [arXiv:1909.04114 [hep-ex]].
  * [51] M. Aaboud et al. [ATLAS], Phys. Lett. B 795, 56-75 (2019) doi:10.1016/j.physletb.2019.03.067 [arXiv:1901.10917 [hep-ex]].
  * [52] A. M. Sirunyan et al. [CMS], Phys. Rev. Lett. 123, no.23, 231803 (2019) doi:10.1103/PhysRevLett.123.231803 [arXiv:1905.10331 [hep-ex]].
  * [53] G. Aad et al. [ATLAS], ATLAS-CONF-2018-052.
  * [54] C. Chang, P. Scott, T. E. Gonzalo, F. Kahlhoefer, A. Kvellestad and M. White, [arXiv:2209.13266 [hep-ph]].
  * [55] E. Bagnaschi, J. C. Costa, K. Sakurai, M. Borsato, O. Buchmueller, A. Roeck, M. J. Dolan, J. R. Ellis, H. Flächer and K. Hahn, et al. Eur. Phys. J. C 79, no.11, 895 (2019) doi:10.1140/epjc/s10052-019-7382-3 [arXiv:1905.00892 [hep-ph]].
  * [56] G. Aad et al. [ATLAS], Phys. Lett. B 796, 68-87 (2019) doi:10.1016/j.physletb.2019.07.016 [arXiv:1903.06248 [hep-ex]].
  * [57] F. Kahlhoefer, A. Mück, S. Schulte and P. Tunney, JHEP 03, 104 (2020) doi:10.1007/JHEP03(2020)104 [arXiv:1912.06374 [hep-ph]].
  * [58] G. Aad et al. [ATLAS], Phys. Rev. D 103, no.11, 112006 (2021) doi:10.1103/PhysRevD.103.112006 [arXiv:2102.10874 [hep-ex]].
  * [59] A. M. Sirunyan et al. [CMS], CMS-PAS-EXO-20-004.
  * [60] D. Barducci, G. Belanger, J. Bernon, F. Boudjema, J. Da Silva, S. Kraml, U. Laa and A. Pukhov, Comput. Phys. Commun. 222, 327-338 (2018) doi:10.1016/j.cpc.2017.08.028 [arXiv:1606.03834 [hep-ph]].
  * [61] A. M. Sirunyan et al. [CMS], [CMS], CMS-PAS-EXO-19-018.
  * [62] R. Aaij et al. [LHCb], Phys. Rev. Lett. 120, no.6, 061801 (2018) doi:10.1103/PhysRevLett.120.061801 [arXiv:1710.02867 [hep-ex]].
  * [63] R. Aaij et al. [LHCb], Phys. Rev. Lett. 124, no.4, 041801 (2020) doi:10.1103/PhysRevLett.124.041801 [arXiv:1910.06926 [hep-ex]].
  * [64] P. Ilten, Y. Soreq, M. Williams and W. Xue, JHEP 06, 004 (2018) doi:10.1007/JHEP06(2018)004 [arXiv:1801.04847 [hep-ph]].
  * [65] C. Baruch, P. Ilten, Y. Soreq and M. Williams, JHEP 11, 124 (2022) doi:10.1007/JHEP11(2022)124 [arXiv:2206.08563 [hep-ph]].
  * [66] J. P. Lees et al. [BaBar], Phys. Rev. Lett. 113, no.20, 201801 (2014) doi:10.1103/PhysRevLett.113.201801 [arXiv:1406.2980 [hep-ex]].
  * [67] J. P. Lees et al. [BaBar], Phys. Rev. Lett. 119, no.13, 131804 (2017) doi:10.1103/PhysRevLett.119.131804 [arXiv:1702.03327 [hep-ex]].
  * [68] S. Abrahamyan et al. [APEX], Phys. Rev. Lett. 107, 191804 (2011) doi:10.1103/PhysRevLett.107.191804 [arXiv:1108.2750 [hep-ex]].
  * [69] H. Merkel, P. Achenbach, C. Ayerbe Gayoso, T. Beranek, J. Bericic, J. C. Bernauer, R. Böhm, D. Bosnar, L. Correa and L. Debenjak, et al. Phys. Rev. Lett. 112, no.22, 221802 (2014) doi:10.1103/PhysRevLett.112.221802 [arXiv:1404.5502 [hep-ex]].
  * [70] D. Banerjee, V. E. Burtsev, A. G. Chumakov, D. Cooke, P. Crivelli, E. Depero, A. V. Dermenev, S. V. Donskov, R. R. Dusaev and T. Enik, et al. Phys. Rev. Lett. 123, no.12, 121801 (2019) doi:10.1103/PhysRevLett.123.121801 [arXiv:1906.00176 [hep-ex]].
  * [71] A. Konaka, K. Imai, H. Kobayashi, A. Masaike, K. Miyake, T. Nakamura, N. Nagamine, N. Sasao, A. Enomoto and Y. Fukushima, et al. Phys. Rev. Lett. 57, 659 (1986) doi:10.1103/PhysRevLett.57.659
  * [72] E. M. Riordan, M. W. Krasny, K. Lang, P. De Barbaro, A. Bodek, S. Dasu, N. Varelas, X. Wang, R. G. Arnold and D. Benton, et al. Phys. Rev. Lett. 59, 755 (1987) doi:10.1103/PhysRevLett.59.755
  * [73] J. D. Bjorken, S. Ecklund, W. R. Nelson, A. Abashian, C. Church, B. Lu, L. W. Mo, T. A. Nunamaker and P. Rassmann, Phys. Rev. D 38, 3375 (1988) doi:10.1103/PhysRevD.38.3375
  * [74] A. Bross, M. Crisler, S. H. Pordes, J. Volk, S. Errede and J. Wrbanek, Phys. Rev. Lett. 67, 2942-2945 (1991) doi:10.1103/PhysRevLett.67.2942
  * [75] M. Davier and H. Nguyen Ngoc, Phys. Lett. B 229, 150-155 (1989) doi:10.1016/0370-2693(89)90174-3
  * [76] T. Bringmann, J. Edsjö, P. Gondolo, P. Ullio and L. Bergström, JCAP 07, 033 (2018) doi:10.1088/1475-7516/2018/07/033 [arXiv:1802.03399 [hep-ph]].
  * [77] N. Aghanim et al. [Planck], Astron. Astrophys. 641, A6 (2020) [erratum: Astron. Astrophys. 652, C4 (2021)] doi:10.1051/0004-6361/201833910 [arXiv:1807.06209 [astro-ph.CO]].
  * [78] L. J. Hall, K. Jedamzik, J. March-Russell and S. M. West, JHEP 03, 080 (2010) doi:10.1007/JHEP03(2010)080 [arXiv:0911.1120 [hep-ph]].
  * [79] G. Bélanger, F. Boudjema, A. Goudelis, A. Pukhov and B. Zaldivar, Comput. Phys. Commun. 231, 173-186 (2018) doi:10.1016/j.cpc.2018.04.027 [arXiv:1801.03509 [hep-ph]].
  * [80] A. Hryczuk and M. Laletin, Phys. Rev. D 106, no.2, 2 (2022) doi:10.1103/PhysRevD.106.023007 [arXiv:2204.07078 [hep-ph]].
  * [81] G. Belanger, A. Mjallal and A. Pukhov, Eur. Phys. J. C 81, no.3, 239 (2021) doi:10.1140/epjc/s10052-021-09012-z [arXiv:2003.08621 [hep-ph]].
  * [82] R. Agnese et al. [SuperCDMS], Phys. Rev. Lett. 116, no.7, 071301 (2016) doi:10.1103/PhysRevLett.116.071301 [arXiv:1509.02448 [astro-ph.CO]].
  * [83] G. Angloher et al. [CRESST], Eur. Phys. J. C 76, no.1, 25 (2016) doi:10.1140/epjc/s10052-016-3877-3 [arXiv:1509.01515 [astro-ph.CO]].
  * [84] A. H. Abdelhameed et al. [CRESST], Phys. Rev. D 100, no.10, 102002 (2019) doi:10.1103/PhysRevD.100.102002 [arXiv:1904.00498 [astro-ph.CO]].
  * [85] P. Agnes et al. [DarkSide], Phys. Rev. D 98, no.10, 102006 (2018) doi:10.1103/PhysRevD.98.102006 [arXiv:1802.07198 [astro-ph.CO]].
  * [86] D. S. Akerib et al. [LUX], Phys. Rev. Lett. 118, no.2, 021303 (2017) doi:10.1103/PhysRevLett.118.021303 [arXiv:1608.07648 [astro-ph.CO]].
  * [87] C. Amole et al. [PICO], Phys. Rev. Lett. 118, no.25, 251301 (2017) doi:10.1103/PhysRevLett.118.251301 [arXiv:1702.07666 [astro-ph.CO]].
  * [88] C. Amole et al. [PICO], Phys. Rev. D 100, no.2, 022001 (2019) doi:10.1103/PhysRevD.100.022001 [arXiv:1902.04031 [astro-ph.CO]].
  * [89] A. Tan et al. [PandaX-II], Phys. Rev. Lett. 117, no.12, 121303 (2016) doi:10.1103/PhysRevLett.117.121303 [arXiv:1607.07400 [hep-ex]].
  * [90] X. Cui et al. [PandaX-II], Phys. Rev. Lett. 119, no.18, 181302 (2017) doi:10.1103/PhysRevLett.119.181302 [arXiv:1708.06917 [astro-ph.CO]].
  * [91] E. Aprile et al. [XENON], Phys. Rev. Lett. 121, no.11, 111302 (2018) doi:10.1103/PhysRevLett.121.111302 [arXiv:1805.12562 [astro-ph.CO]].
  * [92] J. Aalbers et al. [LZ], [arXiv:2207.03764 [hep-ex]].
  * [93] D. Hooper, D. P. Finkbeiner and G. Dobler, Phys. Rev. D 76, 083012 (2007) doi:10.1103/PhysRevD.76.083012 [arXiv:0705.3655 [astro-ph]].
  * [94] T. Sjostrand, S. Mrenna and P. Z. Skands, Comput. Phys. Commun. 178, 852-867 (2008) doi:10.1016/j.cpc.2008.01.036 [arXiv:0710.3820 [hep-ph]].
  * [95] M. Ackermann et al. [Fermi-LAT], Phys. Rev. Lett. 115, no.23, 231301 (2015) doi:10.1103/PhysRevLett.115.231301 [arXiv:1503.02641 [astro-ph.HE]].
  * [96] A. Alloul, N. D. Christensen, C. Degrande, C. Duhr and B. Fuks, Comput. Phys. Commun. 185, 2250-2300 (2014) doi:10.1016/j.cpc.2014.04.012 [arXiv:1310.1921 [hep-ph]].
  * [97] C. Degrande, Comput. Phys. Commun. 197, 239-262 (2015) doi:10.1016/j.cpc.2015.08.015 [arXiv:1406.3030 [hep-ph]].
  * [98] T. Hahn, Comput. Phys. Commun. 140, 418-431 (2001) doi:10.1016/S0010-4655(01)00290-9 [arXiv:hep-ph/0012260 [hep-ph]].
  * [99] B. Fuks and R. Ruiz, JHEP 05, 032 (2017) doi:10.1007/JHEP05(2017)032 [arXiv:1701.05263 [hep-ph]].
  * [100] T. Sjöstrand, S. Ask, J. R. Christiansen, R. Corke, N. Desai, P. Ilten, S. Mrenna, S. Prestel, C. O. Rasmussen and P. Z. Skands, Comput. Phys. Commun. 191, 159-177 (2015) doi:10.1016/j.cpc.2015.01.024 [arXiv:1410.3012 [hep-ph]].
  * [101] C. Bierlich, S. Chakraborty, N. Desai, L. Gellersen, I. Helenius, P. Ilten, L. Lönnblad, S. Mrenna, S. Prestel and C. T. Preuss, et al. [arXiv:2203.11601 [hep-ph]].
  * [102] J. de Favereau et al. [DELPHES 3], JHEP 02, 057 (2014) doi:10.1007/JHEP02(2014)057 [arXiv:1307.6346 [hep-ex]].
  * [103] M. Cacciari, G. P. Salam and G. Soyez, Eur. Phys. J. C 72, 1896 (2012) doi:10.1140/epjc/s10052-012-1896-2 [arXiv:1111.6097 [hep-ph]].
  * [104] M. Cacciari, G. P. Salam and G. Soyez, JHEP 04, 063 (2008) doi:10.1088/1126-6708/2008/04/063 [arXiv:0802.1189 [hep-ph]].
  * [105] A. Hocker, P. Speckmayer, J. Stelzer, J. Therhaag, E. von Toerne, H. Voss, M. Backes, T. Carli, O. Cohen and A. Christov, et al. [arXiv:physics/0703039 [physics.data-an]].
  * [106] P. Speckmayer, A. Hocker, J. Stelzer and H. Voss, J. Phys. Conf. Ser. 219, 032057 (2010) doi:10.1088/1742-6596/219/3/032057
  * [107] I. Antcheva, M. Ballintijn, B. Bellenot, M. Biskup, R. Brun, N. Buncic, P. Canal, D. Casadei, O. Couet and V. Fine, et al. Comput. Phys. Commun. 180, 2499-2512 (2009) doi:10.1016/j.cpc.2009.08.005 [arXiv:1508.07749 [physics.data-an]].
  * [108] I. Antcheva, M. Ballintijn, B. Bellenot, M. Biskup, R. Brun, N. Buncic, P. Canal, D. Casadei, O. Couet and V. Fine, et al. Comput. Phys. Commun. 182, 1384-1385 (2011) doi:10.1016/j.cpc.2011.02.008
  * [109] J. L. Feng, I. Galon, F. Kling and S. Trojanowski, Phys. Rev. D 97, no.3, 035001 (2018) doi:10.1103/PhysRevD.97.035001 [arXiv:1708.09389 [hep-ph]].
  * [110] A. Berlin and F. Kling, Phys. Rev. D 99, no.1, 015021 (2019) doi:10.1103/PhysRevD.99.015021 [arXiv:1810.01879 [hep-ph]].
  * [111] L. A. Anchordoqui, A. Ariga, T. Ariga, W. Bai, K. Balazs, B. Batell, J. Boyd, J. Bramante, M. Campanelli and A. Carmona, et al. Phys. Rept. 968, 1-50 (2022) doi:10.1016/j.physrep.2022.04.004 [arXiv:2109.10905 [hep-ph]].
  * [112] J. L. Feng, F. Kling, M. H. Reno, J. Rojo, D. Soldin, L. A. Anchordoqui, J. Boyd, A. Ismail, L. Harland-Lang and K. J. Kelly, et al. [arXiv:2203.05090 [hep-ex]].
  * [113] K. Akiba et al. [LHC Forward Physics Working Group], J. Phys. G 43, 110201 (2016) doi:10.1088/0954-3899/43/11/110201 [arXiv:1611.05079 [hep-ph]].
  * [114] F. Kling and S. Trojanowski, Phys. Rev. D 104, no.3, 035012 (2021) doi:10.1103/PhysRevD.104.035012 [arXiv:2105.07077 [hep-ph]].
  * [115] A. Ariga et al. [FASER], [arXiv:1811.10243 [physics.ins-det]].
  * [116] A. Ariga et al. [FASER], Phys. Rev. D 99, no.9, 095011 (2019) doi:10.1103/PhysRevD.99.095011 [arXiv:1811.12522 [hep-ph]].
  * [117] A. Ariga et al. [FASER], [arXiv:1812.09139 [physics.ins-det]].
  * [118] X. Cid Vidal, M. D’Onofrio, P. J. Fox, R. Torre, K. A. Ulmer, A. Aboubrahim, A. Albert, J. Alimena, B. C. Allanach and C. Alpigiani, et al. CERN Yellow Rep. Monogr. 7, 585-865 (2019) doi:10.23731/CYRM-2019-007.585 [arXiv:1812.07831 [hep-ph]].
  * [119] E. Todesco and F. Zimmermann, doi:10.5170/CERN-2011-003 [arXiv:1111.7188 [physics.acc-ph]].
  * [120] M. Mangano, doi:10.23731/CYRM-2017-003 [arXiv:1710.06353 [hep-ph]].
  * [121] Y. D. Tsai, P. deNiverville and M. X. Liu, Phys. Rev. Lett. 126, no.18, 181801 (2021) doi:10.1103/PhysRevLett.126.181801 [arXiv:1908.07525 [hep-ph]].
  * [122] J. Blumlein and J. Brunner, Phys. Lett. B 701, 155-159 (2011) doi:10.1016/j.physletb.2011.05.046 [arXiv:1104.2747 [hep-ex]].
  * [123] J. Blümlein and J. Brunner, Phys. Lett. B 731, 320-326 (2014) doi:10.1016/j.physletb.2014.02.029 [arXiv:1311.3870 [hep-ph]].
  * [124] S. Andreas, C. Niebuhr and A. Ringwald, Phys. Rev. D 86, 095019 (2012) doi:10.1103/PhysRevD.86.095019 [arXiv:1209.6083 [hep-ph]].
  * [125] D. Banerjee et al. [NA64], Phys. Rev. D 101, no.7, 071101 (2020) doi:10.1103/PhysRevD.101.071101 [arXiv:1912.11389 [hep-ex]].
  * [126] J. R. Batley et al. [NA48/2], Phys. Lett. B 746, 178-185 (2015) doi:10.1016/j.physletb.2015.04.068 [arXiv:1504.00607 [hep-ex]].
  * [127] N. Baltzell et al. [HPS], Nucl. Instrum. Meth. A 859, 69-75 (2017) doi:10.1016/j.nima.2017.03.061 [arXiv:1612.07821 [physics.ins-det]].
  * [128] E. Kou et al. [Belle-II], PTEP 2019, no.12, 123C01 (2019) [erratum: PTEP 2020, no.2, 029201 (2020)] doi:10.1093/ptep/ptz106 [arXiv:1808.10567 [hep-ex]].
  * [129] C. Ahdida et al. [SHiP], Eur. Phys. J. C 81, no.5, 451 (2021) doi:10.1140/epjc/s10052-021-09224-3 [arXiv:2011.05115 [hep-ex]].
  * [130] A. Berlin, S. Gori, P. Schuster and N. Toro, Phys. Rev. D 98, no.3, 035011 (2018) doi:10.1103/PhysRevD.98.035011 [arXiv:1804.00661 [hep-ph]].
  * [131] K. Asai, A. Das, J. Li, T. Nomura and O. Seto, Phys. Rev. D 106, no.9, 095033 (2022) doi:10.1103/PhysRevD.106.095033 [arXiv:2206.12676 [hep-ph]].
  * [132] K. Cheung, C. J. Ouseph and T. Wang, JHEP 12, 209 (2021) doi:10.1007/JHEP12(2021)209 [arXiv:2111.08375 [hep-ph]].
  * [133] K. Cheung and C. J. Ouseph, JHEP 10, 196 (2022) doi:10.1007/JHEP10(2022)196 [arXiv:2208.04523 [hep-ph]].
  * [134] A. Das, P. S. B. Dev, Y. Hosotani and S. Mandal, Phys. Rev. D 105, no.11, 115030 (2022) doi:10.1103/PhysRevD.105.115030 [arXiv:2104.10902 [hep-ph]].