o m #1 #2#2 #1 o [#1]Pythia o [#1]Whizard o [#1]Circe2 o [#1]HELAC-Onia

Snowmass Preprint

Prospects for searches for Higgs boson decays to dark photons at the ILC

S. Snyder1, C. Weber1, and D. Zhang2

1Brookhaven National Laboratory, Upton, NY, 11973
2Institute of High Energy Physics, Chinese Academy of Sciences, Shijingshan
District, Beijing, 100049

> ABSTRACT
>
>
>
>
> An interesting model of dark matter involves a hidden sector decoupled from
> Standard Model (SM) fields except for some portal interaction. A concrete
> realization of this is the Hidden Abelian Higgs Model, which gives rise to
> decays of the SM Higgs boson into a pair of new bosons, called
> $Z_{\text{d}}$ or dark photons. This note explores prospects for the search
> for such dark photons at the ILC with $\sqrt{s}=$250\text{\,}\mathrm{GeV}$$,
> where the dark photons decay promptly. For the $H\rightarrow
> Z_{\text{d}}Z_{\text{d}}\rightarrow 4\ell$ ($\ell=e,\mu$) final state, it
> follows closely recent similar searches at the LHC, while for the $2\ell 2j$
> and $4j$ final states a multivariate analysis approach is used. This study
> has not been approved by the SiD consortium.

Submitted to the Proceedings of the US Community Study

on the Future of Particle Physics (Snowmass 2021)

##  1 Introduction

A major open question in particle physics is the nature of the
astrophysically-motivated dark matter. An attractive strategy for
incorporating dark matter into the Standard Model (SM) is through a hidden
sector, decoupled from known SM fields except for some ‘portal’ interaction
[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11]. A concrete realization of this is the
Hidden Abelian Higgs Model (HAHM) [5, 6, 7, 8, 9, 10], in which a new
$U(1)_{d}$ dark gauge field kinetically mixes with the SM $U(1)_{Y}$
hypercharge gauge field with some strength $\epsilon$ [12, 13, 14]. This gives
rise to a new Higgs-like dark scalar $S$ along with the gauge boson of the new
field, $Z_{\text{d}}$, or ‘dark photon’. The scalar $S$ mixes with the SM
Higgs boson with strength $\kappa$, allowing decays of the SM Higgs boson into
pairs of $Z_{\text{d}}$ bosons via mixing with the $S$ scalar. For
$\epsilon\ll 1$, the decays of the $Z_{\text{d}}$ boson are largely determined
by the gauge couplings. Over the range
$$1\text{\,}\mathrm{GeV}$<m_{Z_{\text{d}}}<$60\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$$,
the branching fraction of the $Z_{\text{d}}$ boson into pairs of electrons or
muons could be 10%–15% [5]. These decays would be prompt for $\epsilon\gtrsim
10^{-5}$.

This note explores the prospects for a search for $H\rightarrow
Z_{\text{d}}Z_{\text{d}}$ at the ILC with
$\sqrt{s}=$250\text{\,}\mathrm{GeV}$$, with the $Z_{\text{d}}$ bosons decaying
promptly to $4\ell$ ($\ell\equiv e,\mu$), as illustrated in Fig. 1, as well as
$2\ell 2j$ and $4j$. For the $4\ell$ final state, it follows closely the
analysis of ATLAS at the LHC with $\sqrt{s}=$13\text{\,}\mathrm{TeV}$$ and an
integrated luminosity of $139\text{\,}{\rm fb}^{-1}$ [15]. Other similar
searches, including searches for pairs of light bosons decaying into muons,
$\tau$-leptons, photons, and/or jets, as well as searches for a single light
boson decaying into a pair of muons, using both
$\sqrt{s}=$8\text{\,}\mathrm{TeV}$$ and $13\text{\,}\mathrm{TeV}$ data, have
been performed by ATLAS [16, 17, 18, 19, 20], CMS [21, 22, 23, 24], and LHCb
[25]. Searches for long-lived signatures at ATLAS and CMS are reported in
Refs. [26, 27, 28, 29, 30, 31, 32, 33, 34], while further searches for a SM
Higgs boson decaying into undetected particles are reported in Refs. [35, 36].

Figure 1: Exotic decay of the Higgs boson into four leptons induced by
intermediate dark vector bosons via the Higgs portal, where $s$ is a dark
Higgs boson [5]. The $Z_{\text{d}}$ gauge boson decays into SM particles
through kinetic mixing with the hypercharge field (with branching ratios that
are nearly independent of $\epsilon$). The $HZ_{\text{d}}Z_{\text{d}}$ vertex
factor is proportional to $\kappa$.

The following section describes how the event samples used for these analyses
were simulated and reconstructed. This is followed by a description of the
$H\rightarrow Z_{\text{d}}Z_{\text{d}}\rightarrow 4\ell$, $2\ell 2j$, and $4j$
analyses, including the event selection, a discussion of quarkonia
backgrounds, and expected limits. The note ends with a summary and discussion
of possible future work.

##  2 Event simulation and reconstruction

The $ee\rightarrow ZH\rightarrow ffZ_{\text{d}}Z_{\text{d}}$ signal was
generated according to the HAHM [9, 10, 6, 5] implementation in
MadGraph5_aMC@NLO version 2.8.2 [37], with the Higgs boson mass set to
$m_{H}=$125\text{\,}\mathrm{GeV}$$ and $\epsilon$ and $\kappa$ both set to
$10^{-4}$. The $e^{-}$ beam polarization was set to $-80\%$ and $e^{+}$ to
$+30\%$. Leptonically decaying $Z_{\text{d}}$ bosons were forced to decay to
either a $ee$ or $\mu\mu$ pair. Final states with $\tau$-leptons were not
included. In the similar LHC analysis [15], the change in signal region yield
due to the omission of these decays was below 1% and thus neglected. Decays of
$Z_{\text{d}}$ bosons to jets were allowed to be inclusive. Other Higgs boson
production mechanisms, which are much smaller at the ILC with
$\sqrt{s}=$250\text{\,}\mathrm{GeV}$$, were also omitted. Showering was
performed via MadGraph5_aMC@NLO’s built-in interface to 8.244 [38]. Events in
the $4\ell$ final state were generated at
$m_{Z_{\text{d}}}=$1\text{\,}\mathrm{GeV}$$, every $2\text{\,}\mathrm{GeV}$ in
the range $$2\text{\,}\mathrm{GeV}$\leq
m_{Z_{\text{d}}}\leq$12\text{\,}\mathrm{GeV}$$, and every
$5\text{\,}\mathrm{GeV}$ in the range $$15\text{\,}\mathrm{GeV}$\leq
m_{Z_{\text{d}}}\leq$60\text{\,}\mathrm{GeV}$$. Signal samples with $2\ell 2j$
and $4j$ final states were produced at $mZd=$20\text{\,}\mathrm{GeV}$$,
$40\text{\,}\mathrm{GeV}$, and $60\text{\,}\mathrm{GeV}$. Each signal sample
was generated with 20,000 events per $m_{Z_{\text{d}}}$ value and per final
state.

The $ee\rightarrow X+H\rightarrow X+4\ell$ background was generated using
2.8.5 [39, 40] along with its internal version of 6.4 [41], with parameters
set corresponding to the OPAL tune provided in [42]. The $e^{-}$ beam
polarization was set to $-80\%$ and $e^{+}$ to $+30\%$, and the provided
`ilc250` parameterization was used. The Higgs boson was forced to decay into
either a $4e$, $2e2\mu$, or $4\mu$ final state. For this background, 20,000
events were generated.

The non-resonant $ee\rightarrow 4e,2e2\mu,4\mu$ background was also generated
with , using the same settings as described previously. To remove divergences,
the final state leptons were restricted to be within
$$20\text{\,}\mathrm{\SIUnitSymbolDegree}$<\theta<$160\text{\,}\mathrm{\SIUnitSymbolDegree}$$;
further, each pair of final state leptons had to have an invariant mass
greater than $1\text{\,}\mathrm{GeV}$. For this background, 250,000 events
were generated.

The generated events were processed through the full iLCSoft simulation and
reconstruction chain [43] using the `o2_V03` version of the SiD [44] geometry.
The electrons and muons used in the analysis were identified by
`IsolatedLeptonTaggingProcessor` [45]. This, however, relies on a cone-based
isolation algorithm that rejects leptons from
$Z_{\text{d}}\rightarrow\ell\ell$ decay if
$m_{Z_{\text{d}}}\lesssim$4\text{\,}\mathrm{GeV}$$. In order extend the
analysis to smaller values of $m_{Z_{\text{d}}}$, the isolation algorithm is
modified. When summing the particle flow objects inside an isolation cone
around an electron or muon candidate, the highest-energy same-flavor opposite-
sign lepton candidate is ignored. This recovers efficiency for
$Z_{\text{d}}\rightarrow\ell\ell$ at low $m_{Z_{\text{d}}}$ with no
significant increase in background for this analysis (as estimated from a
$b{\bar{b}}$ sample generated with ).

The $2\ell 2j$ and $4j$ final states rely on inclusive $ee\rightarrow ZH$ and
$ee\rightarrow ZZ$ samples from the SiD collaboration [46]. These were
generated using 2.6.4 using nominal ILC Technical Design Report polarization
fractions of 80% polarized electrons and 30% polarized positrons. However the
study is limited to samples with $-80\%$ and $e^{+}$ to $+30\%$ polarization
to match electron beam polarization of the signal samples, and has not been
approved by the SiD consortium.

##  3 $H\rightarrow Z_{\text{d}}Z_{\text{d}}\rightarrow 4\ell$ event selection

The event selection closely follows that of the ATLAS $H\rightarrow
Z_{\text{d}}Z_{\text{d}}\rightarrow 4\ell$ analysis [15], with a few
modifications for the ILC environment, such as removing detector-specific
particle-identification requirements and changing requirements on
$p_{\text{T}}$ and $\eta$ to $E$ and $\theta$. Although each signal event
contains a $Z$ boson decay in addition to the Higgs boson decay, no explicit
requirements are made for the $Z$ boson decay. Two sets of selections are
used. The first, called the high-mass (HM) selection is designed for
$$15\text{\,}\mathrm{GeV}$<m_{Z_{\text{d}}}<$60\text{\,}\mathrm{GeV}$$, and
the second, low-mass (LM), selection is designed for
$$1\text{\,}\mathrm{GeV}$<m_{Z_{\text{d}}}<$15\text{\,}\mathrm{GeV}$$. In this
latter region, the angular separation between the two leptons from
$Z_{\text{d}}\rightarrow\ell\ell$ decay becomes small. For the ATLAS analysis
[15], the LM analysis used only the $4\mu$ final state, due to a decreased
efficiency for identifying closely-spaced electrons. Simulated ILC events do
not show such a drop in efficiency, so here all final states are used for the
LM analysis.

Electrons and muons must satisfy $E>$7\text{\,}\mathrm{GeV}$$ and
$0.35<\theta<\pi-0.35$. Such leptons are formed into quadruplets consisting of
two same-flavor opposite-sign (SFOS) lepton pairs, giving $4e$, $2e2\mu$, and
$4\mu$ final states. If there are more than two such pairs, multiple
quadruplets are formed from all possible SFOS combinations. The invariant
masses of the two pairs are denoted by $m_{12}$ and $m_{34}$, with $m_{12}$
taken to be the one closest to the mass of the $Z$ boson:
$|m_{12}<m_{Z}|<|m_{34}-m_{Z}|$.

If all four leptons in a quadruplet have the same flavor, then one can also
define alternate pairings. The invariant mass $m_{14}$ is defined from the
positively charged lepton of the $m_{12}$ pair and the negatively charged
lepton of the $m_{34}$ pair. The other alternative pairing $m_{23}$ is defined
similarly.

For all quadruplets, the three highest-energy leptons must satisfy
$E_{1}>$20\text{\,}\mathrm{GeV}$$, $E_{2}>$15\text{\,}\mathrm{GeV}$$, and
$E_{3}>$10\text{\,}\mathrm{GeV}$$. For the HM analysis only, the angular
separation between all same-flavor leptons must satisfy $\Delta
R(\ell,\ell^{\prime})>0.1$ and for different-flavor leptons $\Delta
R(\ell,\ell^{\prime})>0.2$, where $(\Delta
R)^{2}=(\Delta\eta)^{2}+(\Delta\phi)^{2}$ and $\eta$ is the
pseudorapidity***Although an analysis in an $e^{-}e^{+}$ environment would
more naturally use $\theta$ than $\eta$, the use of $\Delta R$ is retained
here for consistency with the ATLAS analysis.. Events are required to have at
least one such quadruplet. If there is more than one, the quadruplet with the
smallest $\Delta m_{\ell\ell}=|m_{12}-m_{34}|$ is used.

For the HM event selection, the invariant mass of the four leptons must be
consistent with that of the SM Higgs boson:
$$115\text{\,}\mathrm{GeV}$<m_{4\ell}<$130\text{\,}\mathrm{GeV}$$. The
quadruplet must also not be consistent with the decay of $Z$ bosons
($Z$-veto): $$10\text{\,}\mathrm{GeV}$<m_{12,34}<$64\text{\,}\mathrm{GeV}$$.
For the $4e$ and $4\mu$ channels, it is possible that the leptons are
mispaired, so for these channels there is also a requirement on the
alternative lepton pairings:
$$5\text{\,}\mathrm{GeV}$<m_{14,23}<$75\text{\,}\mathrm{GeV}$$. Events with
lepton pairs consistent with $J/\psi$ or $\Upsilon$ decay are also rejected.
Events are rejected if any of $m_{12,34,14,23}$ are in the ranges
$(m_{J/\psi}-$0.25\text{\,}\mathrm{GeV}$)$ to
$(m_{\psi(2S)}+$0.30\text{\,}\mathrm{GeV}$)$ or
$(m_{\Upsilon(2S)}-$0.70\text{\,}\mathrm{GeV}$)$ to
$(m_{\Upsilon(3S)}+$0.75\text{\,}\mathrm{GeV}$)$, where the quarkonia masses
are taken to be $m_{J/\psi}=$3.096\text{\,}\mathrm{GeV}$$,
$m_{\psi(2S)}=$3.686\text{\,}\mathrm{GeV}$$,
$m_{\Upsilon(1S)}=$9.461\text{\,}\mathrm{GeV}$$, and
$m_{\Upsilon(3S)}=$10.355\text{\,}\mathrm{GeV}$$ [47]. Finally, a requirement
$m_{34}/m_{12}>0.85$ ensures that the two pairs have similar invariant masses.

For the LM event selection, the pair invariant masses are required to be in
the range $$0.8\text{\,}\mathrm{GeV}$<m_{12,34}<$20\text{\,}\mathrm{GeV}$$.
The requirement on the overall invariant mass is tightened to
$$120\text{\,}\mathrm{GeV}$<m_{4\ell}<$130\text{\,}\mathrm{GeV}$$ due to
smaller radiative tails in this regime. The $Z$-veto requirement is not
applied, and and only the alternate lepton pairings are used for the quarkonia
vetoes (see Section 4). The final requirement $m_{34}/m_{12}>0.85$ is the same
as for the HM analysis.

Both event selections are summarized in Section 3.

Table 1: Summary of event selection requirements for the HM and LM analyses.

|  |  High-mass (HM) analysis [t] |
---|---|---|---
$H\rightarrow Z_{\text{d}}Z_{\text{d}}\rightarrow 4\ell$ $(\ell=e,\mu)$ [b] | Low-mass (LM) analysis |  |
$H\rightarrow Z_{\text{d}}Z_{\text{d}}\rightarrow 4\ell$ $(\ell=e,\mu)$ |  |  |
Mass range |  $$15\text{\,}\mathrm{GeV}$<m_{X}<$60\text{\,}\mathrm{GeV}$$ |  $$1\text{\,}\mathrm{GeV}$<m_{X}<$15\text{\,}\mathrm{GeV}$$
Leptons |  Four isolated electrons or muons[t]
with $E>$7\text{\,}\mathrm{GeV}$$ and $0.35<\theta<\pi-0.35$[b] |  |  |
Quadruplet selection |  $e^{+}e^{-}e^{+}e^{-}$, $e^{+}e^{-}\mu^{+}\mu^{-}$, or $\mu^{+}\mu^{-}\mu^{+}\mu^{-}$;[t]
Three leading-$E$ leptons satisfying $E>$20\text{\,}\mathrm{GeV}$$, $15\text{\,}\mathrm{GeV}$, $10\text{\,}\mathrm{GeV}$ |  |  |
Define pairs $m_{12}$ and $m_{34}$ such that $|m_{12}-m_{Z}|<|m_{34}-m_{Z}|$ [b] |  |  |
|  |  [t]$\Delta R(\ell,\ell^{\prime})>0.10\,(0.20)$ for same-flavor (different-flavor) $\ell,\ell^{\prime}$ |  —
Quadruplet ranking  |  Select quadruplet with smallest $\Delta m_{\ell\ell}=|m_{12}-m_{34}|$
Event selection |  $m_{4\ell}$ |  $$115\text{\,}\mathrm{GeV}$<m_{4\ell}<$130\text{\,}\mathrm{GeV}$$ |  $$120\text{\,}\mathrm{GeV}$<m_{4\ell}<$130\text{\,}\mathrm{GeV}$$
$Z$-veto |  $$10\text{\,}\mathrm{GeV}$<m_{12,34}<$64\text{\,}\mathrm{GeV}$$[t] |
For $4e$ and $4\mu$ channels: $$5\text{\,}\mathrm{GeV}$<m_{14,23}<$75\text{\,}\mathrm{GeV}$$ [b] | — |  |
| Heavy-flavor veto |  Reject event if $m_{12,34,14,23}$ in:[t] |
$(m_{J/\psi}-$0.25\text{\,}\mathrm{GeV}$)$ to $(m_{\psi(2S)}+$0.30\text{\,}\mathrm{GeV}$)$, or |  |  |
$(m_{\Upsilon(1S)}-$0.70\text{\,}\mathrm{GeV}$)$ to $(m_{\Upsilon(3S)}+$0.75\text{\,}\mathrm{GeV}$)$ [b] |  Reject event if $m_{14,23}$ in:[t] |  |
$(m_{J/\psi}-$0.25\text{\,}\mathrm{GeV}$)$ to $(m_{\psi(2S)}+$0.30\text{\,}\mathrm{GeV}$)$, or |  |  |
$(m_{\Upsilon(1S)}-$0.70\text{\,}\mathrm{GeV}$)$ to $(m_{\Upsilon(3S)}+$0.75\text{\,}\mathrm{GeV}$)$ [b] |  |  |
|  Signal region |  $m_{34}/m_{12}>0.85$ |  $$0.8\text{\,}\mathrm{GeV}$<m_{12,34}<$20\text{\,}\mathrm{GeV}$$[t]
$m_{34}/m_{12}>0.85$ [b] |  |  |

The main backgrounds are from SM $H\rightarrow ZZ^{*}\rightarrow 4\ell$ decay
and also from nonresonant $4\ell$ production, with the latter source
dominating. These are estimated using the simulated samples described in
Section 2. Backgrounds in which jets are misidentified as leptons are assumed
to be negligible. Estimated backgrounds for a data sample of
$2000\text{\,}{\rm fb}^{-1}$ are shown in Table 2.

Table 2: Estimated backgrounds for the $H\rightarrow Z_{\text{d}}Z_{\text{d}}\rightarrow 4\ell$ analyses for a data sample of $2000\text{\,}{\rm fb}^{-1}$. No simulated events for the non-resonant background pass the LM selection, so a limit is given based on the background represented by one simulated event. HM selection | $4e$ | $2e2\mu$ | $4\mu$ | All
---|---|---|---|---
[t] $H\rightarrow ZZ^{*}\rightarrow 4\ell$ | $0.37\pm 0.04$ | $0.39\pm 0.04$ | $0.46\pm 0.04$ | $1.22\pm 0.07$
Non-resonant | $3.05\pm 1.15$ | $2.62\pm 1.07$ | $0.44\pm 0.44$ | $6.10\pm 1.63$
[b] Total | $3.42\pm 1.15$ | $3.00\pm 1.07$ | $0.90\pm 0.44$ | $7.32\pm 1.63$
LM selection | $4e$ | $2e2\mu$ | $4\mu$ | All
[t] $H\rightarrow ZZ^{*}\rightarrow 4\ell$ | $0.04\pm 0.01$ | $<0.01$ | $0.03\pm 0.01$ | $0.07\pm 0.02$
Non-resonant | $<0.44$ | $<0.44$ | $<0.44$ | $<0.44$
[b] Total | $<0.47$ | $<0.44$ | $<0.47$ | $<0.51$

##  4 Quarkonia backgrounds to $H\rightarrow
Z_{\text{d}}Z_{\text{d}}\rightarrow 4\ell$ final states

The ATLAS $H\rightarrow Z_{\text{d}}Z_{\text{d}}\rightarrow 4\ell$ analysis
[15] has no sensitivity in the regions
$$2\text{\,}\mathrm{GeV}$<m_{Z_{\text{d}}}<$4.4\text{\,}\mathrm{GeV}$$ and
$$8\text{\,}\mathrm{GeV}$<m_{Z_{\text{d}}}<$12\text{\,}\mathrm{GeV}$$ due to
the presence of large backgrounds from quarkonia production. However, at the
ILC, hadronic backgrounds such as this will be much smaller, and the excellent
lepton momentum resolution of the ILC detectors may reduce the size of the
$m_{Z_{\text{d}}}$ ranges affected by these backgrounds. Unfortunately, there
is no reliable, general-purpose simulation of quarkonia production at the ILC.
But one can still estimate these backgrounds, as described below.

First, consider direct, non-resonant production of $J/\psi$ and $\Upsilon$
pairs. These processes were estimated using version 2.0.1 [48, 49]. This can
calculate processes such as $e^{-}e^{+}\rightarrow J/\psi J/\psi$ and
$e^{-}e^{+}\rightarrow J/\psi J/\psi+q{\bar{q}}$, and similarly for
$\Upsilon$. The quarkonia can be produced as either color singlets or color
octets. For the purpose of this study, excited quarkonia states are not
considered. The initial beams are set to $e^{-}e^{+}$ with beam energies of
$125\text{\,}\mathrm{GeV}$ each, with initial-state radiation disabled. (does
not implement beam polarization effects, and enabling initial-state radiation
caused the subsequent showering step to fail.) For the processes with the
largest cross sections, generated events were then showered and hadronized
with  8.244, with the quarkonia forced to decay to either $e^{-}e^{+}$ or
$\mu^{-}\mu^{+}$. The quarkonia decay branching ratios were taken to be
$\mbox{\rm BR}(J/\psi\rightarrow\ell\ell)=0.119$ and $\mbox{\rm
BR}(\Upsilon\rightarrow\ell\ell)=0.049$ [47]. Events were then passed through
the detector simulation and analysis, yielding estimates of these quarkonia
backgrounds after the LM selection, shown in Table 3. These estimates are all
much smaller than other backgrounds for this selection.

Table 3: Quarkonia production processes calculated with . The 1 and 8 superscripts denote color-singlet and color-octet states, respectively. For the processes with larger cross sections, the number of expected events passing the LM selection for a data sample of $2000\text{\,}{\rm fb}^{-1}$ is shown. This is shown as a limit for cases where no simulated events pass the selection. Process | Cross section (ab) | Expected background
---|---|---
[t] $J/\psi^{{\bf 1}}J/\psi^{{\bf 1}}$ | $8.3$ | $<5.7\text{\times}{10}^{-6}$
$J/\psi^{{\bf 8}}J/\psi^{{\bf 8}}$ | $2.9\text{\times}{10}^{-6}$ |
$J/\psi^{{\bf 1}}J/\psi^{{\bf 1}}+d{\bar{d}}$ | $0.0033$ |
$J/\psi^{{\bf 8}}J/\psi^{{\bf 8}}+d{\bar{d}}$ | $11$ | $<0.0018$
$J/\psi^{{\bf 1}}J/\psi^{{\bf 1}}+u{\bar{u}}$ | $0.0017$ |
$J/\psi^{{\bf 8}}J/\psi^{{\bf 8}}+u{\bar{u}}$ | $45$ | $<0.0034$
$J/\psi^{{\bf 1}}J/\psi^{{\bf 1}}+s{\bar{s}}$ | $0.0033$ |
$J/\psi^{{\bf 8}}J/\psi^{{\bf 8}}+s{\bar{s}}$ | $11$ | $<0.0018$
$J/\psi^{{\bf 1}}J/\psi^{{\bf 1}}+c{\bar{c}}$ | $0.057$ |
$J/\psi^{{\bf 8}}J/\psi^{{\bf 8}}+c{\bar{c}}$ | $41$ | $6.3\text{\times}{10}^{-4}$
$J/\psi^{{\bf 1}}J/\psi^{{\bf 1}}+b{\bar{b}}$ | $0.0022$ |
$J/\psi^{{\bf 8}}J/\psi^{{\bf 8}}+b{\bar{b}}$ | $8.3$ | $1.0\text{\times}{10}^{-4}$
$J/\psi^{{\bf 1}}J/\psi^{{\bf 1}}+gg$ | $0.012$ |
$J/\psi^{{\bf 8}}J/\psi^{{\bf 8}}+gg$ | $1.6$ | †††failed to generate events for this process.
[b] $J/\psi^{{\bf 8}}J/\psi^{{\bf 8}}+g$ | $0.060$ | $1.1\text{\times}{10}^{-6}$
[t] $\Upsilon^{{\bf 1}}\Upsilon^{{\bf 1}}$ | $0.18$ | $<1.3\text{\times}{10}^{-7}$
$\Upsilon^{{\bf 8}}\Upsilon^{{\bf 8}}$ | $6.5\text{\times}{10}^{-11}$ |
$\Upsilon^{{\bf 1}}\Upsilon^{{\bf 1}}+d{\bar{d}}$ | $8.1\text{\times}{10}^{-4}$ |
$\Upsilon^{{\bf 8}}\Upsilon^{{\bf 8}}+d{\bar{d}}$ | $1.0\text{\times}{10}^{-5}$ |
$\Upsilon^{{\bf 1}}\Upsilon^{{\bf 1}}+u{\bar{u}}$ | $0.0034$ |
$\Upsilon^{{\bf 8}}\Upsilon^{{\bf 8}}+u{\bar{u}}$ | $4.0\text{\times}{10}^{-5}$ |
$\Upsilon^{{\bf 1}}\Upsilon^{{\bf 1}}+s{\bar{s}}$ | $8.1\text{\times}{10}^{-4}$ |
$\Upsilon^{{\bf 8}}\Upsilon^{{\bf 8}}+s{\bar{s}}$ | $1.0\text{\times}{10}^{-5}$ |
$\Upsilon^{{\bf 1}}\Upsilon^{{\bf 1}}+c{\bar{c}}$ | $0.0031$ |
$\Upsilon^{{\bf 8}}\Upsilon^{{\bf 8}}+c{\bar{c}}$ | $3.9\text{\times}{10}^{-5}$ |
$\Upsilon^{{\bf 1}}\Upsilon^{{\bf 1}}+b{\bar{b}}$ | $0.0022$ |
$\Upsilon^{{\bf 8}}\Upsilon^{{\bf 8}}+b{\bar{b}}$ | $8.4\text{\times}{10}^{-6}$ |
$\Upsilon^{{\bf 1}}\Upsilon^{{\bf 1}}+gg$ | $4.7\text{\times}{10}^{-5}$ |
$\Upsilon^{{\bf 8}}\Upsilon^{{\bf 8}}+gg$ | $8.0\text{\times}{10}^{-5}$ |
[b] $\Upsilon^{{\bf 1}}\Upsilon^{{\bf 1}}+g$ | $2.9\text{\times}{10}^{-7}$ |

Another possibility is the decay of a $Z$ boson into a quarkonia pair, or two
$Z$ bosons decaying to the same quarkonium state. Backgrounds from these
processes were estimated using current experimental limits/measurements for
such decays [47]:

$\mbox{\rm BR}(Z\rightarrow J/\psi+X)$ | $5.1\text{\times}{10}^{-3}$
---|---
$\mbox{\rm BR}(Z\rightarrow J/\psi J/\psi)$ | $<2.2\text{\times}{10}^{-6}$
$\mbox{\rm BR}(Z\rightarrow\Upsilon+X)$ | $1.0\text{\times}{10}^{-4}$
$\mbox{\rm BR}(Z\rightarrow\Upsilon\Upsilon)$ | $<1.5\text{\times}{10}^{-6}$

was used to generate a sample of $e^{-}e^{+}\rightarrow f{\bar{f}}Z$, where
$f$ is any fermion, configured as described in Section 2, except that the
showering was set to force the $Z$ boson to decay as either $Z\rightarrow
J/\psi J/\psi\rightarrow 4\ell$ or $Z\rightarrow\Upsilon\Upsilon\rightarrow
4\ell$. Events were then passed through the detector simulation and LM
selection. The total cross section calculated by for $f{\bar{f}}Z$ was
$7.1\text{\times}{10}^{3}\text{\,}\mathrm{f}\mathrm{b}$. Taking into account
the branching ratios above and the efficiency of the LM selection, the
expected background for a $2000\text{\,}{\rm fb}^{-1}$ data sample is
$<$2.2\text{\times}{10}^{-5}\text{\,}$$ for $Z\rightarrow J/\psi J/\psi$ and
$1.2\text{\times}{10}^{-5}\text{\,}$ for $Z\rightarrow\Upsilon\Upsilon$.

Similarly, was also used to generate a sample of $e^{-}e^{+}\rightarrow ZZ$,
where here each $Z$ boson was forced to decay as either $Z\rightarrow J/\psi
c{\bar{c}}$ or $Z\rightarrow\Upsilon b{\bar{b}}$. The total cross section for
this process was $1.8\text{\times}{10}^{3}\text{\,}\mathrm{f}\mathrm{b}$.
Again, taking into account branching ratios and efficiencies, the expected
background for a $2000\text{\,}{\rm fb}^{-1}$ data sample is 0.0029 for
$ZZ\rightarrow J/\psi J/\psi+c{\bar{c}}c{\bar{c}}$ and
$4.3\text{\times}{10}^{-7}\text{\,}$ for
$ZZ\rightarrow\Upsilon\Upsilon+b{\bar{b}}b{\bar{b}}$.

A final possibility is $H\rightarrow J/\psi J/\psi$ or
$H\rightarrow\Upsilon\Upsilon$. Such a background would be particularly
concerning since it could not be removed by the requirement that the overall
invariant mass be consistent with that of the SM Higgs boson. However,
calculations give $BR(H\rightarrow J/\psi
J/\psi)=$5.9\text{\times}{10}^{-10}\text{\,}$$ and
$BR(H\rightarrow\Upsilon\Upsilon)=$4.3\text{\times}{10}^{-10}\text{\,}$$ [50,
51], which are much smaller than the $H\rightarrow Z_{\text{d}}Z_{\text{d}}$
branching ratio to which the LM analysis is sensitive.

Although this does not exhaust all possibilities for quarkonia background
processes, it should be a representative sample. All processes examined result
in backgrounds that are much smaller than the other (already-small)
backgrounds to the LM selection. Therefore, quarkonia production is unlikely
to be a significant background to this analysis. However, evaluating this with
more confidence would likely require improved codes for calculating quarkonia
processes.

Further, the invariant mass distributions of dilepton decays of quarkonia are
shown in Section 3. These events were generated by +and processed with the
full detector simulation. (sets the quarkonia masses to be exactly the sum of
the masses of the constituent quarks, so the positions of the peaks are
shifted from the true quarkonia masses.) These peaks are very narrow,
especially for the $\mu^{-}\mu^{+}$ decays. Therefore, even if quarkonia
backgrounds were to be significant, they could be effectively suppressed by
rejecting a much smaller range in $m_{Z_{\text{d}}}$ than was done in the
ATLAS analysis.

(a)  (c)  (d)

Figure 2: Invariant mass distributions for quarkonia decays to lepton pairs as
generated by +: (a) $J/\psi\rightarrow e^{-}e^{+}$; (b)
$J/\psi\rightarrow\mu^{-}\mu^{+}$; (c) $\Upsilon\rightarrow e^{-}e^{+}$; (d)
$\Upsilon\rightarrow\mu^{-}\mu^{+}$. Note that sets the quarkonia masses to be
exactly the sum of the constituent quark masses; hence, the peaks are shifted
from the true quarkonia masses.

##  5 $H\rightarrow Z_{\text{d}}Z_{\text{d}}\rightarrow 4\ell$ expected limits

Expected limits are set based on the distribution of the average of the
invariant masses of the two lepton pairs in a quadruplet,
$\left<m_{\ell\ell}\right>=\frac{1}{2}\left(m_{12}+m_{34}\right)$. The
likelihood function describing the data for a channel $j$ consists of a
Poisson factor for each histogram bin $i$:

${\cal L}(N)=\prod_{i}\textrm{Pois}\left(N_{ij};\mu
S_{ij}(m_{Z_{\text{d}}})+B_{ij}\right),$ (1)

where $S$ and $B$ are the predicted numbers of signal and background events
for each bin and channel and $\mu$ is the signal strength. For this study,
systematic uncertainties are assumed to be negligible. The signal shape as a
function of $m_{Z_{\text{d}}}$ is found by fitting a Gaussian to the simulated
signal at each generated mass point and then interpolating in the fit mean and
width as a function of $m_{Z_{\text{d}}}$. Background histograms are smoothed
using the `RooKeysPdf` class of RooFit [52, 53], except that if there are less
then ten simulated events surviving in a channel, the background is taken to
be flat with respect to $m_{Z_{\text{d}}}$.

Following [15], a set of generator-level fiducial requirements, described in
LABEL:tab:4lfiducial, are used to factorize the event selection into a largely
model-independent ‘efficiency’ and a model-dependent ‘acceptance’. For the
purpose of these selections, the four-momenta of photons close to a lepton
($\Delta R<0.1$) are added to that of the lepton. This accounts for the
effects of quasi-collinear electromagnetic radiation from the leptons [54].
The efficiency for a channel is defined as the fraction of events passing the
generator-level fiducial selection that also passes the full event selection,
while the acceptance is defined as the fraction of generator-level events that
pass the fiducial selection. The efficiency and acceptance for the analyses
described here are shown in LABEL:fig:4l-effacc. For the HM selection, the
acceptance falls for low $m_{Z_{\text{d}}}$ for the $4e$ and $4\mu$ channels
due to the alternate pair requirement of the $Z$-veto. (Similar behavior was
seen in the ATLAS analysis [15].)

Table 4: Summary of the fiducial phase-space definitions for the HM and LM
analyses. Objects are considered at generator-level, with photons nearby
leptons summed with those leptons.  (c)

Figure 4: Per-channel expected upper limits at 95% CL on the fiducial cross
sections for the $H\rightarrow Z_{\text{d}}Z_{\text{d}}\rightarrow 4\ell$
process, for the (a) $4e$, (b) $4\mu$, and (c) $2e2\mu$ final states. The
discontinuities at $m_{Z_{\text{d}}}=$15\text{\,}\mathrm{GeV}$$ are due to the
change from the LM to HM selection.

(a)

(b)

Figure 5: (a) Expected upper limit at 95% CL for the cross section of the
$e^{+}e^{-}\rightarrow H+X\rightarrow Z_{\text{d}}Z_{\text{d}}+X\rightarrow
4\ell+X$ process, assuming SM Higgs boson production. (b) Expected upper limit
at 95% CL for the cross section times the model-dependent branching ratio
divided by the SM Higgs boson production cross section
($319\text{\,}\mathrm{f}\mathrm{b}$) for the $H\rightarrow
Z_{\text{d}}Z_{\text{d}}$ process for the benchmark HAHM. In both cases, all
final states are combined. The discontinuities at
$m_{Z_{\text{d}}}=$15\text{\,}\mathrm{GeV}$$ are due to the change from the LM
to HM selection.

Compared to limits from the similar ATLAS analysis with $139\text{\,}{\rm
fb}^{-1}$ of data and $\sqrt{s}=$13\text{\,}\mathrm{TeV}$$, the expected
branching ratio limits here are a factor of 5–10 higher. This is not
unexpected: since the background is quite small, even at the LHC, the
sensitivity is driven mainly by the total number of Higgs bosons produced,
which was about ten times larger at the LHC than would be expected in
$2000\text{\,}{\rm fb}^{-1}$ of ILC data at
$\sqrt{s}=$250\text{\,}\mathrm{GeV}$$. The exception is in the mass ranges
$$2\text{\,}\mathrm{GeV}$<m_{Z_{\text{d}}}<$4.4\text{\,}\mathrm{GeV}$$ and
$$8\text{\,}\mathrm{GeV}$<m_{Z_{\text{d}}}<$12\text{\,}\mathrm{GeV}$$, where
the ATLAS analysis has no sensitivity due to quarkonia backgrounds.

##  6 $H\rightarrow Z_{\text{d}}Z_{\text{d}}\rightarrow 2\ell 2j$, $4j$ event
selection

The event selection for the $2\ell 2j$ and $4j$ final states proceeds in two
steps, beginning with a cut-based preselection followed by a selection based
on boosted decision trees (BDT).

The $2\ell 2j$ final state requires at least one pair of opposite-sign
electrons or muons along with four jets formed from particle flow objects
(PFO) not associated with isolated leptons or photons. If there are multiple
eligible lepton pairs, the one with invariant mass closest to the $Z$-boson
mass $m_{Z}$ is selected as the $Z$ or $Z_{d}$ candidate. The $4j$ final state
has no requirement on the number of isolated leptons, but requires six jets
from non-isolated PFO objects. Either case uses the standard jet
reconstruction to four or six jets, respectively.

For both final states, all possible jet pairs are constructed. The invariant
masses of the jet pairs, along with the dilepton’s invariant mass in the
$2\ell 2j$ final state, are compared against $m_{Z}$ and the pair with the
invariant mass closest to $m_{Z}$ is selected as the $Z$-boson candidate. If
in the $2\ell 2j$ final state preselection a jet pair is selected as the
$Z$-candidate, the remaining jet and lepton pairs are selected as the $Z_{d}$
candidates. In all other cases all possible jet pairs from the remaining four
jets are constructed and the two jet pairs with invariant mass closest to each
other, i.e. minimizing $|m_{j_{1}j_{2}}-m_{j_{3}j_{4}}|$, are selected as the
two $Z_{d}$ candidates.

The jets, as well as leptons in the $2\ell 2j$ case, must satisfy the
requirement on their polar angle $|\cos(\theta)|<0.9$.

The final requirement for the preselection is that the four fermions
constituting the two $Z_{d}$ candidates have a total invariant mass broadly
consistent with that of the Higgs boson:
$$90\text{\,}\mathrm{GeV}$<m_{ffff}<$160\text{\,}\mathrm{GeV}$$.

The preselection efficiency $\times$ acceptance for the signal samples is
around 31% for the $4j$ final state and 24% for the $2\ell 2j$ one. For the
background samples the efficiencies $\times$ acceptances are 17% and 2% for
the $4j$ and $2\ell 2j$ final states, respectively.

The signal regions are defined by boosted decision trees individually trained
for each final state. Half of the generated events for both the background and
signal samples are randomly assigned for training, while the remainder are
used for evaluation. The input variables for the BDT are:

  * •

The transverse momentum, total energy, invariant mass, and $\cos(\theta)$ for
each $Z$-boson and $Z_{d}$-boson candidate;

  * •

the $\Delta R$ between each possible boson candidate pair;

  * •

and the transverse momentum, total energy, invariant mass of the Higgs boson
candidate.

The efficiencies after the preselection and the BDT selection for the signals
and background are shown in Table 5, and Table 6 shows the signal region event
yields assuming an integrated luminosity of $2000\text{\,}{\rm fb}^{-1}$ and
an $H\rightarrow Z_{\text{d}}Z_{\text{d}}\rightarrow 2\ell 2j$, $4j$ cross
section of $1\text{\,}\mathrm{fb}$.

Table 5: Selection efficiencies $\times$ acceptance for the $4j$ and $2\ell 2j$ signal regions after preselection and boosted decision tree selection. |  signal efficiency $\times$ acceptance |  background efficiency $\times$ acceptance
---|---|---
[t] $m_{Z_{d}}$ | $20\text{\,}\mathrm{GeV}$ | $40\text{\,}\mathrm{GeV}$ | $60\text{\,}\mathrm{GeV}$ |
$4j$ final state | 9.4% | 5.2% | 4.7 % | 0.02%
$2\ell 2j$ final state | 24% | 22% | 24% | 0.6%
Table 6: Expected signal region event yields for the $2\ell 2j$ and $4j$ final states. Both are scaled to an integrated luminosity of $2000\text{\,}{\rm fb}^{-1}$. The signal yields assume a $H\rightarrow Z_{\text{d}}Z_{\text{d}}\rightarrow 2\ell 2j$, $4j$ cross section of $1\text{\,}\mathrm{fb}$. | signal yields | background yield
---|---|---
[t] $m_{Z_{d}}$ | $20\text{\,}\mathrm{GeV}$ | $40\text{\,}\mathrm{GeV}$ | $60\text{\,}\mathrm{GeV}$ |
$4j$ final state | 187 | 103 | 93 | 8400
$2\ell 2j$ final state | 484 | 448 | 487 | 131

##  7 $H\rightarrow Z_{\text{d}}Z_{\text{d}}\rightarrow 2\ell 2j$, $4j$
expected limits

The likelihood function describing the data for the $4j$ and $2\ell 2j$ final
states follows Eq. 1. Both final state channels are fitted concurrently so the
model in this case is given by:

${\cal L}(N)=\prod_{j=4j\text{,}2l2j}\prod_{i}\textrm{Pois}\left(N_{ij};\mu
S_{ij}(m_{Z_{\text{d}}})+B_{ij}\right).$ (2)

The distributions used to evaluate the likelihood are the average
$Z_{\text{d}}$ mass:
$\left<m_{Z_{d}}\right>=\frac{1}{2}\left(m_{Z_{\text{d1}}}+m_{Z_{\text{d2}}}\right)$.
The limits on the total cross section for the combined $4j$ and $2\ell 2j$
branching ratio are shown in Fig. 5a.

Fig. 5b shows that the $2\ell 2j$ and $4j$ final states do not yield a
stricter expected limit on $\mbox{\rm BR}(H\to Z_{\text{d}}Z_{\text{d}})$ than
the $4\ell$ final state, despite the more favorable $\mbox{\rm
BR}(Z_{\text{d}}\to 2j)$ branching ratio [6].

(a)

(b)

Figure 6: (a) Expected upper limit at 95% CL for the cross section of the
$e^{+}e^{-}\rightarrow H+Z\rightarrow Z_{\text{d}}Z_{\text{d}}+Z\rightarrow
4j,2\ell 2j+Z$ process, assuming SM Higgs boson production. (b) Expected upper
limit at 95% CL on the $H\rightarrow Z_{\text{d}}Z_{\text{d}}$ branching
ratio, derived from the $e^{+}e^{-}\rightarrow H+Z\rightarrow
Z_{\text{d}}Z_{\text{d}}+Z\rightarrow 4j,2\ell 2j+Z$ cross section.

##  8 Summary and future work

Expected limits have been presented for a search for dark photons in the
$H\rightarrow Z_{\text{d}}Z_{\text{d}}\rightarrow 4\ell$, $2\ell 2j$, and $4j$
final states. It is seen that compared to the LHC, searches for these channels
are not competitive at the ILC, except for $Z_{\text{d}}$ masses close to the
$J/\psi$ and $\Upsilon$ quarkonia resonances. The LHC analysis does not have
sensitivity in those region due to large hadronic backgrounds, but it should
be possible to derive limits in those regions at the ILC where such
backgrounds are much smaller. Doing this reliably, however, will likely
require progress in theoretical calculations of decays involving quarkonia
states.

## Acknowledgments

This work is supported in part by the U.S. Department of Energy under contract
DE-AC02-98CH10886 with Brookhaven National Laboratory.

## References

  * [1] Pierre Fayet “Light spin-1/2 or spin-0 dark matter particles” In _Phys. Rev. D_ 70, 2004, pp. 023514 DOI: 10.1103/PhysRevD.70.023514
  * [2] Douglas P. Finkbeiner and Neal Weiner “Exciting dark matter and the INTEGRAL/SPI 511 keV signal” In _Phys. Rev. D_ 76, 2007, pp. 083519 DOI: 10.1103/PhysRevD.76.083519
  * [3] Nima Arkani-Hamed, Douglas P. Finkbeiner, Tracy R. Slatyer and Neal Weiner “A theory of dark matter” In _Phys. Rev. D_ 79, 2009, pp. 015014 DOI: 10.1103/PhysRevD.79.015014
  * [4] Emilian Dudas, Yann Mambrini, Stefan Pokorski and Alberto Romagnoni “Extra U(1) as natural source of a monochromatic gamma ray line” In _JHEP_ 10, 2012, pp. 123 DOI: 10.1007/JHEP10(2012)123
  * [5] David Curtin, Rouven Essig, Stefania Gori and Jessie Shelton “Illuminating dark photons with high-energy colliders” In _JHEP_ 02, 2015, pp. 157 DOI: 10.1007/JHEP02(2015)157
  * [6] David Curtin et al. “Exotic decays of the 125 GeV Higgs boson” In _Phys. Rev. D_ 90.7, 2014, pp. 075004 DOI: 10.1103/PhysRevD.90.075004
  * [7] Hooman Davoudiasl, Hye-Sung Lee, Ian Lewis and William J. Marciano “Higgs decays as a window into the dark sector” In _Phys. Rev. D_ 88.1 APS, 2013, pp. 015022 DOI: 10.1103/PhysRevD.88.015022
  * [8] Hooman Davoudiasl, Hye-Sung Lee and William J. Marciano ““Dark” $Z$ implications for parity violation, rare meson decays, and Higgs physics” In _Phys. Rev. D_ 85, 2012, pp. 115019 DOI: 10.1103/PhysRevD.85.115019
  * [9] James D. Wells “How to Find a Hidden World at the Large Hadron Collider”, 2008 arXiv:0803.1243 [hep-ph]
  * [10] Shrihari Gopalakrishna, Sunghoon Jung and James D. Wells “Higgs boson decays to four fermions through an abelian hidden sector” In _Phys. Rev. D_ 78, 2008, pp. 055002 DOI: 10.1103/PhysRevD.78.055002
  * [11] Jim Alexander “Dark Sectors 2016 Workshop: Community Report”, 2016 arXiv:1608.08632 [hep-ph]
  * [12] Peter Galison and Aneesh Manohar “Two Z’s or not two Z’s?” In _Phys. Lett. B_ 136, 1984, pp. 279 DOI: 10.1016/0370-2693(84)91161-4
  * [13] Bob Holdom “Two U(1)’s and $\epsilon$ charge shifts” In _Phys. Lett. B_ 166, 1986, pp. 196 DOI: 10.1016/0370-2693(86)91377-8
  * [14] Keith R. Dienes, Christopher F. Kolda and John March-Russell “Kinetic mixing and the supersymmetric gauge hierarchy” In _Nucl. Phys. B_ 492, 1997, pp. 104–118 DOI: 10.1016/S0550-3213(97)00173-9
  * [15] ATLAS Collaboration “Search for Higgs bosons decaying into new spin-0 or spin-1 particles in four-lepton final states with the ATLAS detector with $139\,\text{fb}^{-1}$ of $pp$ collision data at $\sqrt{s}=13\,\text{TeV}$”, 2021 arXiv:2110.13673 [hep-ex]
  * [16] ATLAS Collaboration “Search for Higgs boson decays into two new low-mass spin-0 particles in the $4b$ channel with the ATLAS detector using $pp$ collisions at $\sqrt{s}=13\,\text{TeV}$” In _Phys. Rev. D_ 102, 2020, pp. 112006 DOI: 10.1103/PhysRevD.102.112006
  * [17] ATLAS Collaboration “Search for Higgs boson decays into pairs of light (pseudo)scalar particles in the $\gamma\gamma jj$ final state in $pp$ collisions at $\sqrt{s}=13\,\text{TeV}$ with the ATLAS detector” In _Phys. Lett. B_ 782, 2018, pp. 750 DOI: 10.1016/j.physletb.2018.06.011
  * [18] ATLAS Collaboration “Search for the Higgs boson produced in association with a vector boson and decaying into two spin-zero particles in the $H\rightarrow aa\rightarrow 4b$ channel in $pp$ collisions at $\sqrt{s}=13\,\text{TeV}$ with the ATLAS detector” In _JHEP_ 10, 2018, pp. 031 DOI: 10.1007/JHEP10(2018)031
  * [19] ATLAS Collaboration “Search for Higgs boson decays into a pair of light bosons in the $bb\mu\mu$ final state in $pp$ collision at $\sqrt{s}=13\,\text{TeV}$ with the ATLAS detector” In _Phys. Lett. B_ 790, 2019, pp. 1 DOI: 10.1016/j.physletb.2018.10.073
  * [20] ATLAS Collaboration “Search for Higgs bosons decaying to $aa$ in the $\mu\mu\tau\tau$ final state in $pp$ collisions at $\sqrt{s}=8\,\text{TeV}$ with the ATLAS experiment” In _Phys. Rev. D_ 92, 2015, pp. 052002 DOI: 10.1103/PhysRevD.92.052002
  * [21] CMS Collaboration “A search for pair production of new light bosons decaying into muons” In _Phys. Lett. B_ 752, 2016, pp. 146 DOI: 10.1016/j.physletb.2015.10.067
  * [22] CMS Collaboration “Search for light bosons in decays of the $125\,\text{GeV}$ Higgs boson in proton–proton collisions at $\sqrt{s}=8\,\text{TeV}$” In _JHEP_ 10, 2017, pp. 076 DOI: 10.1007/JHEP10(2017)076
  * [23] CMS Collaboration “Search for a light pseudoscalar Higgs boson in the boosted $\mu\mu\tau\tau$ final state in proton–proton collisions at $\sqrt{s}=13\,\text{TeV}$” In _JHEP_ 08, 2020, pp. 139 DOI: 10.1007/JHEP08(2020)139
  * [24] CMS Collaboration “Search for long-lived particles decaying into muon pairs in proton–proton collisions at $\sqrt{s}\ =13\,\text{TeV}$ collected with dedicated high-rate data stream”, 2021 arXiv:2112.13769 [hep-ex]
  * [25] LHCb Collaboration “Search for Dark Photons Produced in 13 TeV $pp$ Collisions” In _Phys. Rev. Lett._ 120.6, 2018, pp. 061801 DOI: 10.1103/PhysRevLett.120.061801
  * [26] ATLAS Collaboration “Search for long-lived neutral particles decaying into lepton jets in proton–proton collisions at $\sqrt{s}=8\,\text{TeV}$ with the ATLAS detector” In _JHEP_ 11, 2014, pp. 088 DOI: 10.1007/JHEP11(2014)088
  * [27] ATLAS Collaboration “A search for prompt lepton-jets in $pp$ collisions at $\sqrt{s}=8\,\text{TeV}$ with the ATLAS detector” In _JHEP_ 02, 2016, pp. 062 DOI: 10.1007/JHEP02(2016)062
  * [28] ATLAS Collaboration “Search for massive, long-lived particles using multitrack displaced vertices or displaced lepton pairs in $pp$ collisions at $\sqrt{s}=8\,\text{TeV}$ with the ATLAS detector” In _Phys. Rev. D_ 92, 2015, pp. 072004 DOI: 10.1103/PhysRevD.92.072004
  * [29] CMS Collaboration “Search for long-lived particles that decay into final states containing two electrons or two muons in proton–proton collisions at $\sqrt{s}=8\,\text{TeV}$” In _Phys. Rev. D_ 91, 2015, pp. 052012 DOI: 10.1103/PhysRevD.91.052012
  * [30] ATLAS Collaboration “Search for long-lived particles in final states with displaced dimuon vertices in $pp$ collisions at $\sqrt{s}=13\,\text{TeV}$ with the ATLAS detector” In _Phys. Rev. D_ 99, 2019, pp. 012001 DOI: 10.1103/PhysRevD.99.012001
  * [31] ATLAS Collaboration “Search for displaced vertices of oppositely charged leptons from decays of long-lived particles in $pp$ collisions at $\sqrt{s}=13\,\text{TeV}$ with the ATLAS detector” In _Phys. Lett. B_ 801, 2020, pp. 135114 DOI: 10.1016/j.physletb.2019.135114
  * [32] ATLAS Collaboration “Search for light long-lived neutral particles produced in $pp$ collisions at $\sqrt{s}=13\,\text{TeV}$ and decaying into collimated leptons or light hadrons with the ATLAS detector” In _Eur. Phys. J. C_ 80, 2020, pp. 450 DOI: 10.1140/epjc/s10052-020-7997-4
  * [33] ATLAS Collaboration “Constraints on mediator-based dark matter and scalar dark energy models using $\sqrt{s}=13\,\text{TeV}$ $pp$ collision data collected by the ATLAS detector” In _JHEP_ 05, 2019, pp. 142 DOI: 10.1007/JHEP05(2019)142
  * [34] CMS Collaboration “A search for pair production of new light bosons decaying into muons in proton–proton collisions at $13\,\text{TeV}$” In _Phys. Lett. B_ 796, 2019, pp. 131 DOI: 10.1016/j.physletb.2019.07.013
  * [35] ATLAS Collaboration “Combination of Searches for Invisible Higgs Boson Decays with the ATLAS Experiment” In _Phys. Rev. Lett._ 122, 2019, pp. 231801 DOI: 10.1103/PhysRevLett.122.231801
  * [36] CMS Collaboration “Search for invisible decays of a Higgs boson produced through vector boson fusion in proton–proton collisions at $\sqrt{s}=13\,\text{TeV}$” In _Phys. Lett. B_ 793, 2019, pp. 520 DOI: 10.1016/j.physletb.2019.04.025
  * [37] J. Alwall et al. “The automated computation of tree-level and next-to-leading order differential cross sections, and their matching to parton shower simulations” In _JHEP_ 07, 2014, pp. 079 DOI: 10.1007/JHEP07(2014)079
  * [38] Torbjörn Sjöstrand et al. “An introduction to PYTHIA 8.2” In _Comput. Phys. Commun._ 191, 2015, pp. 159 DOI: 10.1016/j.cpc.2015.01.024
  * [39] Wolfgang Kilian, Thorsten Ohl and Jurgen Reuter “WHIZARD: Simulating Multi-Particle Processes at LHC and ILC” In _Eur. Phys. J. C_ 71, 2011, pp. 1742 DOI: 10.1140/epjc/s10052-011-1742-y
  * [40] Mauro Moretti, Thorsten Ohl and Jurgen Reuter “O’Mega: An Optimizing matrix element generator”, 2001, pp. 1981–2009 arXiv:hep-ph/0102195
  * [41] Torbjorn Sjöstrand, Stephen Mrenna and Peter Z. Skands “PYTHIA 6.4 physics and manual” In _JHEP_ 05, 2006, pp. 026 DOI: 10.1088/1126-6708/2006/05/026
  * [42] Daniel Jeans and Chris Potter “Snowmass Energy Frontier ILC Analysis Walkthrough” Talk given at the Oct. 2020 Snowmass tutorial URL: https://indico.fnal.gov/event/45721/contributions/198054/attachments/136223/169374/snowmass-walkthrough-v3.pdf
  * [43] authors “iLCSoft Project Page” URL: https://github.com/iLCSoft
  * [44] Halina Abramowicz “The International Linear Collider Technical Design Report - Volume 4: Detectors”, 2013 arXiv:1306.6329 [physics.ins-det]
  * [45] Claude Duerig and Junping Tian “A new MVA based Isolated Lepton Tagging Processor” URL: https://github.com/iLCSoft/MarlinReco/tree/master/Analysis/IsolatedLeptonTagging
  * [46] Chris Potter “SiD Simulation & Analysis for ILC Snowmass Physics LoIs” In _arXiv preprint arXiv:2106.00819_ , 2021
  * [47] P.A. Zyla “Review of Particle Physics” In _Prog. Theor. Exp. Phys._ 2020, 2020, pp. 083C01 DOI: 10.1093/ptep/ptaa104
  * [48] Hua-Sheng Shao “HELAC-Onia 2.0: an upgraded matrix-element and event generator for heavy quarkonium physics” In _Comput. Phys. Commun._ 198, 2016, pp. 238–259 DOI: 10.1016/j.cpc.2015.09.011
  * [49] Hua-Sheng Shao “HELAC-Onia: An automatic matrix element generator for heavy quarkonium physics” In _Comput. Phys. Commun._ 184, 2013, pp. 2562–2570 DOI: 10.1016/j.cpc.2013.05.023
  * [50] Dao-Neng Gao and Xi Gong “Higgs boson decays into a pair of heavy vector quarkonia”, 2022 arXiv:2203.00514 [hep-ph]
  * [51] V. Kartvelishvili, A.. Luchinsky and A.. Novoselov “Double vector quarkonia production in exclusive Higgs boson decays” In _Phys. Rev. D_ 79, 2009, pp. 114015 DOI: 10.1103/PhysRevD.79.114015
  * [52] Wouter Verkerke and David Kirkby “The RooFit toolkit for data modeling”, 2003 arXiv:physics/0306116 [physics.data-an]
  * [53] Kyle S. Cranmer “Kernel estimation in high-energy physics” In _Comput. Phys. Commun._ 136, 2001, pp. 198–207 DOI: 10.1016/S0010-4655(00)00243-5
  * [54] ATLAS Collaboration “Proposal for particle-level object and observable definitions for use in physics measurements at the LHC”, ATL-PHYS-PUB-2015-013, 2015 URL: https://cds.cern.ch/record/2022743
  * [55] Alexander L. Read “Presentation of search results: the $CL_{S}$ technique” In _J. Phys. G_ 28, 2002, pp. 2693 DOI: 10.1088/0954-3899/28/10/313
  * [56] Glen Cowan, Kyle Cranmer, Eilam Gross and Ofer Vitells “Asymptotic formulae for likelihood-based tests of new physics” In _Eur. Phys. J. C_ 71, 2011, pp. 1554 DOI: 10.1140/epjc/s10052-011-1554-0