capbtabboxtable[][]

aainstitutetext: Center for AI and Natural Sciences, KIAS, Seoul 02455,
Republic of Koreabbinstitutetext: Department of Physics, Chungbuk National
University, Cheongju, Chungbuk 28644, Koreaccinstitutetext: Institut für
Theoretische Physik und Astrophysik, Uni Würzburg, Emil-Hilb-Weg 22, D-97074
Würzburg, Germanyddinstitutetext: School of Physics, KIAS, Seoul 02455,
Republic of Koreaeeinstitutetext: Center for Theoretical Physics of the
Universe, Institute for Basic Science, Daejeon 34126, Republic of Korea

# Uncovering doubly charged scalars with dominant three-body decays using
machine learning

Thomas Flacke<EMAIL_ADDRESS>b,d,e Jeong Han Kim<EMAIL_ADDRESS>c
Manuel Kunkel<EMAIL_ADDRESS>a,d Pyungwon Ko
<EMAIL_ADDRESS>b Jun Seung Pi<EMAIL_ADDRESS>c Werner Porod
<EMAIL_ADDRESS>c Leonard Schwarze  leonard.schwarze@stud-
mail.uni-wuerzburg.de

###### Abstract

We propose a deep learning-based search strategy for pair production of doubly
charged scalars undergoing three-body decays to $W^{+}t\bar{b}$ in the same-
sign lepton plus multi-jet final state. This process is motivated by composite
Higgs models with an underlying fermionic UV theory. We demonstrate that for
such busy final states, jet image classification with convolutional neural
networks outperforms standard fully connected networks acting on reconstructed
kinematic variables. We derive the expected discovery reach and exclusion
limit at the high-luminosity LHC.

††preprint: KIAS-A23009

##  1 Introduction

The first two runs of the Large Hadron Collider (LHC) have shown that the
Standard Model (SM) of particle physics continues to be an excellent
description of Nature. Many of its features have been tested to an extremely
high accuracy by collider experiments over the last four decades. The
discovery of the Higgs boson, a crucial prediction of the SM, by the ATLAS and
CMS collaborations at the LHC in 2012 ATLAS:2012yve ; CMS:2012qbp  has
strengthened confidence in the ability of the SM to accurately describe
subatomic physics. Despite its enormous success, several important features of
Nature are unexplained, thus motivating the search for physics beyond the SM.
Among them are: (i) The origin of the electroweak sector: why does the Higgs
sector have the SM-form, in particular why is the squared mass parameter
$\mu^{2}$ negative? (ii) Are there any additional scalars besides the observed
Higgs boson? (iii) Is the observed Higgs boson a fundamental scalar or does it
have an internal structure?

An interesting path for explaining the origin of electroweak symmetry breaking
is given in composite Higgs models Kaplan:1983fs ; Kaplan:1983sm ;
Contino:2003ve ; Agashe:2004rs . Composite extensions of the SM are very well
motivated due to their similarity with QCD as an example of a strongly coupled
theory yielding bound states at low scales which is realised in nature. In
many composite Higgs models, the Higgs boson is a pseudo Nambu-Goldstone boson
(pNGB), which explains its relatively small mass compared to other states
predicted in these models. A new gauge symmetry “hypercolour” is postulated
with new “hyperquarks” charged under this symmetry. The corresponding coupling
constant, although small at high energies, becomes large at a scale
$\Lambda_{\mathrm{HC}}\sim O($TeV$)$ where a condensate is formed. It is
further assumed that a global symmetry group $G$ is present above
$\Lambda_{\mathrm{HC}}$ which gets broken to a subgroup $H$ by the condensate,
yielding at least four Nambu-Goldstone bosons which form – within the SM
description – the corresponding Higgs doublet. In a second step the couplings
of the SM gauge fields, quarks, and leptons are taken into account which
amount to an explicit breaking of $G$. In this way a potential for the pNGBs
is radiatively generated, which then leads to electroweak symmetry breaking.
The hierarchical Yukawa couplings can be explained via a mechanism known as
partial compositeness, in which the SM fermions mix with composite states
Kaplan:1991dc .

The so-called “minimal” composite Higgs model Agashe:2004rs ; Agashe:2005dk
is based on the symmetry breaking pattern SO(5)/SO(4), where the only four
pNGBs precisely match the four degrees of freedom of the SM Higgs doublet
field. The phenomenology of this model has been widely explored, and we refer
the interested reader to the excellent reviews Contino:2010rs ;
Bellazzini:2014yua ; Panico:2015jxa . The minimal model, and variations
thereof, can be considered as a useful template for understanding the
phenomenology of a composite pNGB Higgs. However, models in this class are not
UV complete. From the point of view of an underlying fermionic gauge theory, à
la QCD, the symmetry breaking pattern SO(5)/SO(4) cannot be obtained as the
global symmetry group of the underlying fermions is unitary
Cacciapaglia:2014uja ; Cacciapaglia:2020kgq . The minimal breaking patterns
for a composite Higgs which arises from models with underlying fermions are
SU(4)/Sp(4), SU(5)/SO(5) or SU(4)$\times$SU(4)/SU(4) Dugan:1984hq ;
Galloway:2010bp ; Ferretti:2013kya ; Cacciapaglia:2014uja ; Ma:2015gra . All
of these models contain additional pNGBs besides the Higgs doublet. Several
aspects of their LHC phenomenology have been investigated in the literature,
see for example Belyaev:2016ftv ; Agugliaro:2018vsu ; Banerjee:2022izw ;
Banerjee:2022xmu ; Cacciapaglia:2022bax  and refs. therein. In particular, one
can obtain strong bounds on the masses of these additional pNGBs in final
states containing photons and multi-lepton final states, where the latter stem
from decays of vector bosons Banerjee:2022xmu ; Cacciapaglia:2022bax . While
these final states are clean at a hadron collider like the LHC, the
corresponding branching ratios can be rather small. In contrast, the branching
ratios for decays into quarks, in particular top and bottom quarks, are
usually significantly larger but the corresponding signals suffer from large
backgrounds. As a consequence, the decays into quarks are less strongly
constrained. Busy final states comprised of multiple jets and leptons face
numerous challenges such as multiple jet reconstructions at the cost of signal
efficiency, and the combinatorial problem to correctly pair them to
resonances. This motivates the development of new methods, and in this article
we address this challenge by using machine learning (ML) techniques.

Deep neural networks (DNN) have proven to be powerful for image recognition
LeCun:2012 ; Krizhevsky:2012 , and have become a ubiquitous tool in high
energy physics to search for new particles Bhimji:2017qvb ; Andrews:2019wng ;
Huang:2022rne ; Kim:2019wns ; Kim:2022miv ; Diefenbacher:2019ezd ; Kim:2023wuk
. Hadronic activity can be visualised in discretised calorimeter grids – so-
called jet images deOliveira:2015xxd ; Cogan:2014oua ; Kasieczka:2019dbj  –
which contain information about the direction and the energy of particles. Jet
images also contain information on the difference in colour flow between the
signal and background events. We note for completeness, that ML techniques
have also been proposed for generic BSM searches, see e.g. Hajer:2018kqm ;
Farina:2018fyg ; Heimel:2018mkt ; Blance:2019ibf ; CrispimRomao:2020ucc ;
Bortolato:2021zic ; Hallin:2021wme ; Jawahar:2021vyu ; Ngairangbam:2021yma ;
Karagiorgi:2021ngt , but they all come with certain theory assumptions and
limitations as has been critically accessed in Finke:2021sdf .

In this article we focus on pair production of doubly charged scalars which
subsequently decay via the three-body decay $S^{++}\rightarrow W^{+}t\bar{b}$
and the conjugate process. This signature is theoretically well motivated in
underlying models of a composite Higgs and has been shown to be only weakly
constrained by current LHC searches Cacciapaglia:2022bax . We focus on the
leptonic decay of either the $W^{+}t$ or the $W^{-}\bar{t}$ system which
yields a same-sign di-lepton signature. Moreover, the signal involves a lot of
hadronic activity from the decays of the other particles in this channel. We
target this hadronic activity by a combination of jet images and kinematic
variables of reconstructed objects and we show that deep neural networks can
provide powerful discrimination between signal and background.

This paper is organised as follows: In Sec. 2 we summarise the theoretical
background underlying the process we study, followed by details on the event
generation and preselection cuts in Sec. 3. Section 4 then describes the data
representation and the machine learning setup. We present our results in Sec.
5, and draw our conclusions and present an outlook in Sec. 6. We provide
additional information in two technical appendices. Appendix A gives details
about the neural network architectures and comparisons with networks not used
in the main text. Finally, in Appendix B we provide additional information
about the evaluation techniques and comparisons between different training
data sets.

##  2 Extended scalar sectors in composite Higgs models

The cosets of composite Higgs models with a fermionic UV completion contain
additional pNGBs besides the Higgs doublet. In case of SU(4)/Sp(4)
Galloway:2010bp ; Cacciapaglia:2014uja ; Ferretti:2016upr  the additional pNGB
is just a gauge singlet. For SU(5)/SO(5) Dugan:1984hq ; Ferretti:2016upr
there is a gauge singlet and three SU(2)L triplets with hypercharges $Y=0,\pm
1$. For SU(4)$\times$SU(4)/SU(4) Ma:2015gra ; Ferretti:2016upr  there are a
SU(2)L triplet with $Y=0$, an additional SU(2)L doublet with $Y=-1$ and three
SU(2)L singlets with $Y=0,\pm 1$. The decay modes of these additional pNGBs
depend on the details of the underlying model but fall into three broad
classes Ferretti:2016upr ; Agugliaro:2018vsu :

  * •

decays into electroweak gauge bosons via Wess-Zumino-Witten (WZW) terms
Wess:1971yu ; Witten:1983tw ,

  * •

decays into an electroweak gauge boson and another pNGB. These dominate over
the decays via WZW terms if the mass splitting between the pNGBs is larger
than a few GeV Banerjee:2022xmu ; Cacciapaglia:2022bax , or

  * •

decays into two quarks in case of models with partial compositeness. Here, one
expects the decays into top and bottom quarks as these have the largest
masses. These usually dominate for mass splitting between the pNGBs of up to
50 GeV Banerjee:2022xmu ; Cacciapaglia:2022bax .

It has been shown in Banerjee:2022xmu ; Cacciapaglia:2022bax  that existing
LHC data put strong bounds on the masses of the electroweak pNGBs of up to 700
GeV if the decays into electroweak gauge bosons dominate. The bounds are
considerably weaker in the case that the decays into quarks dominate. In the
following, we investigate to which extent these bounds can be improved by
employing ML techniques. For concreteness, we focus on the SU(5)/SO(5) coset.
After electroweak symmetry breaking the triplets can be placed into multiplets
of the remaining custodial SU(2): a singlet $\eta_{1}$, a triplet $\eta_{3}$,
and a quintuplet $\eta_{5}$. The latter contains a doubly charged member
$S^{++}$ and a singly charged $S^{+}$, which will serve as reference states in
the following analysis. One expects that the members of a multiplet have
approximately the same mass Ferretti:2016upr ; Agugliaro:2018vsu , so we take
$m_{S^{++}}=m_{S^{+}}$.

Figure 1: The pair production of doubly charged scalars (left), and mixed pNGB
productions (right) via Drell-Yan processes.

At the LHC the electroweak pNGBs are produced in Drell-Yan processes as shown
in Fig. 1. One usually assumes that the $S^{++}$ either decays into two
$W^{+}$ or into two leptons which is motivated by seesaw type II models
Schechter:1980gr . While these decay channels have been searched for
ATLAS:2021jol ; Leban:2022pdz ; CMS:2017pet , in composite Higgs models they
can be subdominant to decays into quarks, which are the focus of this work.
The singly charged scalar decays according to $S^{+}\to t\bar{b}$. Although
the $S^{++}$ does not have a direct coupling to quarks, it decays via an off-
shell $S^{+}$ (see also Fig. 1)

$\displaystyle S^{++}\to W^{+}S^{+*}\to W^{+}t\bar{b}\,.$ (1)

In case of rather large mass splitting of 40 or more GeV among the pNGBs, the
decay into an off-shell $W$-boson and an on-shell $S^{+}$ becomes important
Cacciapaglia:2022bax . At the LHC, these scalars can either be produced as
decay products of heavier bound states like top partners, see e.g.
Banerjee:2022xmu ; Kim:2018mks ; Kim:2019oyh , or directly via Drell-Yan
processes. In what follows, we concentrate on the second option as it might be
that bound states other than the pNGBs are too heavy to be produced at the
LHC. In the model-specific examples below, the two scalars $S^{++}$ and
$S^{+}$ are assumed to have the same quantum numbers as in the quintuplet in
the SU(5)/SO(5) coset Agugliaro:2018vsu ; Cacciapaglia:2022bax .

##  3 Event generation and preselection cuts

The electroweak pNGBs described in Sec. 2 open up new and under-explored final
states at the LHC. Typical signatures comprise of multiple electroweak bosons
alongside with top and bottom quarks, giving rise to very busy final states.
However, what remains to be addressed is how we will identify such a
considerable amount of jets and leptons at the cost of signal efficiency, and
whether we can distinguish these signatures from SM backgrounds. To account
for these challenges, we need to develop new deep learning-based strategies
beyond traditional cut-and-count methods to increase the potential discovery
reach of these new particles.

As a first step towards this path, we consider the Drell-Yan pair production
of doubly charged scalars

$q\bar{q}\rightarrow
S^{++}S^{--}\rightarrow(W^{+}t\bar{b})\;(W^{-}\bar{t}b)\rightarrow(W^{+}W^{+}b\bar{b})\;(W^{-}W^{-}b\bar{b})\;,$
(2)

where the three-body decays are mediated by an off-shell $S^{\pm*}$ as shown
in Fig. 1 (left). In the scenario where the pNGBs are allowed to couple to SM
fermions, the pNGBs dominantly decay into third generation quarks because the
couplings with electroweak gauge bosons are one-loop suppressed. Once the top
quarks decay, therefore, the final state consists of four $W$ bosons and four
bottom quarks, where the same-sign $W$ bosons arise from the same decay chain.

In the presence of several other electroweak pNGBs in the same multiplet,
however, this final state has similar characteristics to the mixed pNGB
production shown in Fig. 1 (right)

$q\bar{q}^{\prime}\rightarrow
S^{\pm\pm}S^{\mp}\rightarrow(Wtb)\;(tb)\rightarrow(WWbb)\;(Wbb)\;,$ (3)

which also yields same-sign $W$ bosons, but in total one $W$ boson fewer than
the process in Eq. 2. Figure 2 shows that the production cross sections of
both processes are of a similar order of magnitude, and their kinematic
distributions are similar because they both result from Drell-Yan production.
Thus, after a series of selection cuts, the mixed pNGB production will likely
contribute to the signal as well.

The relative weight of the production cross sections between the two processes
depends on the specific realisation of the model. Here, to incorporate this
dependence, we explore two scenarios: i) when only the doubly charged scalar
production is present; ii) when both the doubly charged scalar and mixed pNGB
productions are present, where their relative cross sections are taken from
the $\eta_{5}$ multiplet described in Sec. 2.

Figure 2:  Drell-Yan production cross sections of $S^{++}S^{--}$ (black) and
$S^{\pm\pm}S^{\mp}$ (red) at $\sqrt{s}=14$ TeV.

Given multiple $W$ bosons in the final states, there are a number of possible
decay channels to explore. To focus on the channel with smaller SM
backgrounds, we consider the case where a pair of same-sign $W$ bosons decay
leptonically, while the remaining $W$ bosons decay hadronically. This gives
rise to the final state with two same-sign leptons, missing transverse
momentum, four $b$-tagged jets, and several light-flavour jets.

For the numerical study of the signal processes in Eqs. 2 and 3, we use the
public FeynRules Alloul:2013bka  implementation of the eVLQ model
Banerjee:2022xmu ; Flacke:eVLQ  which provides generic interactions of the
pNGBs and vector-like quarks. We use this implementation to generate a UFO
library Degrande:2011ua  which allows for event generation within the
MadGraph5_aMC@NLO Alwall:2014hca ; Hirschi:2015iia  framework. We generate the
signal events for a centre-of-mass energy of $\sqrt{s}=14$ TeV at the leading
order (LO) in QCD with NNPDF2.3QED parton distribution functions Ball:2013hta
using dynamical renormalisation and factorisation scales. We apply a flat
$K$-factor of 1.15 to the cross sections to account for NLO effects
Fuks:2013vua . In Fig. 2, we show the production cross sections as a function
of the mass of the pNGBs. In this analysis, we assume that all pNGBs are mass
degenerate.

We generate SM backgrounds for a centre-of-mass energy of $\sqrt{s}=14$ TeV at
the next-to-leading-order (NLO) in QCD (unless mentioned otherwise) using the
same PDF and scale choices. The irreducible $t\bar{t}t\bar{t}$ background is
generated at LO accuracy of which the cross section is normalised to the NLO
cross section of $1.17\times 10^{-2}$ pb. The next most important backgrounds
are $t\bar{t}h$ and $t\bar{t}V$, where $V$ denotes an electroweak gauge boson
$W$ or $Z$, whose NLO cross sections are $5.60\times 10^{-1}$ pb and $1.60$
pb, respectively. A subdominant background is given by $t\bar{t}VV$ with NLO
cross section of $1.89\times 10^{-2}$ pb. Other backgrounds, such as $VV$ and
$VVV$, turn out to be negligibly small (after preselection cuts and after
taking branching ratios into account), and we therefore we do not consider
them further in the following analysis.

After generating hard scattering events, the events are passed to Pythia8
Sjostrand:2014zea  for showering and hadronisation. We use Delphes 3.4.1
deFavereau:2013fsa  to include detector resolution effects based on modified
ATLAS configurations Kim:2019wns . For jet reconstruction, we utilise the
Fastjet 3.3.1 Cacciari:2011ma  implementation with an anti-$k_{T}$ algorithm
and a cone radius of $r=0.4$. We benchmark a HL-LHC $b$-tagging efficiency
from ATLAS report CERN-LHCC-2017-021 , where we use a flat $b$-tag rate of
$\epsilon_{b\rightarrow b}=0.8$, and a mistag rate for a $c$-jet (light-flavor
jet) being misidentified as a $b$-jet of $\epsilon_{c\rightarrow b}=0.2$
($\epsilon_{j\rightarrow b}=0.01$).

We require isolated leptons to satisfy
$p_{T}(\ell)/(p_{T}(\ell)+\sum_{i}p_{T,i})>0.7$ where $\sum_{i}p_{T,i}$ is the
sum of the transverse momenta of nearby particles with $p_{T,i}>0.5$ GeV and
$\Delta R_{i\ell}<0.3$. Isolated leptons are required to pass minimum cuts
$p_{T}(\ell)>20$ GeV and $|\eta(\ell)|<2.5$. The charge misidentification of
leptons is taken into account by Delphes.

Although it can be useful to fully identify final state objects to reconstruct
heavy resonances, it is typically accompanied by a poor signal acceptance due
to the preselection criteria described above. Besides, the combinatorial
problem makes it difficult to correctly associate all jets to the particles
they originate from, given the very busy final state. To cope with this, we
propose a more inclusive search by requiring at least three $b$-tagged jets
($N(b)\geq 3$) and at least three non-$b$-tagged jets ($N(j)\geq 3$). Finally,
events are required to have exactly two same-sign leptons, and to pass the
minimum missing transverse momentum (defined as in Ref. Kim:2019wns )
${\;/\\!\\!\\!\\!{P}_{T}}=|{\;/\\!\\!\\!\\!\vec{P}_{T}}|>20$ GeV, and $S_{T}>$
400 GeV where $S_{T}$ is the scalar sum of the transverse momenta of the
reconstructed jets and the two same-sign leptons.

The signal and background cross sections after this baseline selection are
summarised in Tab. 1. In this table, we take the reference pNGB mass for
signal events to be $m_{S}=400$ GeV. The first and second columns show the
list of signal and backgrounds and their corresponding efficiencies of the
preselection. The third and fourth columns show the cross sections and number
of events for the target integrated luminosity ${\cal L}$ = 3 ab-1 of the HL-
LHC.

Process | $\epsilon_{\mathrm{Preselection}}$ | Cross section [fb] | Events at 3 ab-1
---|---|---|---
$S^{++}S^{--}$ | $9.87\times 10^{-3}$ | $4.90\times 10^{-2}$ | 147
$S^{\pm\pm}S^{\mp}$ | $4.81\times 10^{-3}$ | $2.87\times 10^{-3}$ | 86
$t\bar{t}V$ | $1.70\times 10^{-4}$ | $2.72\times 10^{-1}$ | 816
$t\bar{t}h$ | $3.75\times 10^{-4}$ | $2.10\times 10^{-1}$ | 629
$t\bar{t}t\bar{t}$ | $1.63\times 10^{-2}$ | $1.91\times 10^{-1}$ | 572
$t\bar{t}VV$ | $1.74\times 10^{-3}$ | $3.29\times 10^{-2}$ | 98
$VVV$ | $2.08\times 10^{-6}$ | $1.05\times 10^{-3}$ | 3
Table 1: Signal and background efficiencies and cross sections after the
preselection cuts. For signal processes, we take the reference pNGB mass to be
$m_{S}=400$ GeV.

##  4 Input data and network architectures

The preselection cuts described in Sec. 3 provide a preparation for a more
sophisticated machine learning analysis. A deep learning-based analysis can
efficiently exploit differences in differential distributions, which, as
compared to traditional cut-and-count methods, increases the discrimination
power in searches for new particles. However, before coming to the description
of the network architecture, we need to decide on the data representation and
data pre-processing suitable for the individual architectures. Possibilities
range from training a DNN on the lowest level detector data to using neural
networks with human-engineered kinematic variables which are targeted to
specific features of the signal event (like for example reconstructing
invariant masses of intermediate states). For the analysis we propose here, we
use as inputs kinematic information of conventional reconstructed objects
(electrons, muons, jets, $b$-tagged jets) as well as a number of simple
derived kinematic quantities (invariant masses of two object pairs, radial
distance of object pairs, missing transverse energy, and total transverse
energy) as is described in more detail in Sec. 4.2. Additionally, we make use
of image recognition techniques for which events are translated to “jet
images” which are described in Sec. 4.1.

To distinguish the signal from the background events based on this data
representation, we have tested several network architectures. For details on
the various architectures and performance comparisons, we refer the interested
reader to the extensive Appendix A. In Sec. 4.3 we summarise the network
architecture which we found to perform best. This turned out to be a marriage
of two smaller networks: a convolutional neural network (CNN) used for
discrimination of jet images, and a fully connected (FC) deep neural network
which only works with kinematic variables. The outputs of both networks are
then interfaced for a combined classification.

###  4.1 Jet images

Jet images are based on the particle flow information in each event. Following
a similar procedure as described in Ref. Huang:2022rne , we divide the
particle flow into charged and neutral particles. The charged particles
include charged hadrons, while the neutral particles consist of neutral
hadrons as well as non-isolated photons. Leptons are removed from both
samples. We define the origin of the ($\eta$, $\phi$) coordinates to be the
centre of the two same-sign isolated leptons. All particles are mapped
accordingly in the ($\eta$, $\phi$) plane. The jet images are discretised into
$50\times 50$ calorimeter grids within a region of $-2.5\leq\eta\leq 2.5$ and
$-\pi\leq\phi\leq\pi$. The intensity of each pixel is given by the total
transverse momentum of particles passing through the pixel.111This information
is obtained after showering and hadronising at the Delphes level. When summing
transverse momenta of hadrons in each pixel, we only consider the hadrons with
$p_{T}>0.7$ GeV. The final jet images have a dimension of $(2\times 50\times
50)$ where 2 denotes the number of channels, charged and neutral particle
images. In order for the neural networks to fully take into account the
spatial correlation between images of final state particles, we project the
two same-sign leptons into the discretised $50\times 50$ calorimeter grids as
well. Combined with the jet images, we have a set of images for visible
particles whose data structure is represented by

$\displaystyle\mathcal{I}_{CN\ell}=\big{(}3\times 50\times 50\big{)}\;,$ (4)

where 3 denotes charged ($C$), neutral ($N$) particle images, and lepton
($\ell$) images.

In Fig. 3 we show an overlay of the jet images of all simulated events in the
background categories (left three columns) as well as the signal categories
(right three columns), where the respective scalar mass is indicated in
parentheses. When comparing the jet images of signal and background processes,
several differences are apparent: The signal images are generally “brighter”
than those of the backgrounds, indicating a larger total $p_{T}$. The
brightness increases with higher mass of the scalars. Furthermore, for the
signal events the leptons are mostly distributed near the origin, while they
are more broadly distributed for the backgrounds. The reason is that in case
of the backgrounds, the same-sign leptons originate from the decays of
different particles. This leads to a larger angular separation compared to the
signals, where both same-sign leptons come from the same $S^{++}$ or $S^{--}$.

Figure 3: Distribution of jet and lepton images of background (left) and
signal (right) processes, where the scalar mass is given in parentheses. The
final state signatures are discretised into $50\times 50$ calorimeter grids
within a region of $-2.5\leq\eta\leq 2.5$ and $-\pi\leq\phi\leq\pi$. The
origin of the coordinate system is the centre of the two same-sign leptons.
The columns show charged and neutral components of the jets, as well as
isolated leptons, respectively. The colourbar indicates the mean $p_{T}$ per
pixel.

One potential issue with using $\mathcal{I}_{CN\ell}$ is that pileup effects
could jeopardise the analysis, as the expected average number of pileup
interactions $\langle\mu\rangle$ at the HL-LHC is $\mathcal{O}(200)$
collisions per bunch crossing ATLAS:2019xli . We refer to Ref. Kim:2019wns
for the semi-realistic examination of pileup effects, and several methods to
mitigate the contamination. Although charged particles can be potentially
cleaned up from pileup by scrutinising the longitudinal vertex information
Bertolini:2014bba , neutral particles cannot be treated the same way. There
are, however, on-going experimental works to mitigate the impact of neutral
pileup collisions CMS:2020ebo ; Bertolini:2014bba . In addition to our
analysis based on the $\mathcal{I}_{CN\ell}$, we therefore also perform a more
conservative analysis based on the image data set excluding neutral particles
whose data structure is represented by

$\displaystyle\mathcal{I}_{C\ell}=\big{(}2\times 50\times 50\big{)}\;.$ (5)

The impact of excluding the neutral particle image will be discussed in Sec.
5.

###  4.2 Kinematic variables

After the basic selections described in Sec. 3, we consider two same-sign
leptons, three leading $b$-tagged jets, and three leading non-$b$-tagged jets
to construct kinematic variables. We do not utilise other inclusive jets in
the analysis. The most basic kinematic variables are four-momenta of the
reconstructed objects. From these, we can construct higher level kinematic
variables which often serve as efficient learnable features for neural
networks. We consider the following set of kinematic variables

$\mathcal{K}=\bigcup_{i\neq j}M_{ij}\cup\bigcup_{i\neq j}\Delta
R_{ij}\cup\bigcup_{i}p_{Ti}\cup\\{{\>/\\!\\!\\!\\!E_{T}},S_{T}\\}\;,$ (6)

where $M_{ij}$ and $\Delta R_{ij}$ denote the invariant mass and the angular
distance between two reconstructed objects $i$ and $j$, $p_{Ti}$ denotes the
transverse momentum of an object $i$, ${\>/\\!\\!\\!\\!E_{T}}$ denotes the
missing transverse energy, and $S_{T}$ is the scalar sum of transverse momenta
of the reconstructed jets and the two same-sign leptons.

###  4.3 Machine learning aspects

Figure 4: A schematic architecture of the neural networks used in this paper.
The separate DNN chain in the upper panel is used only when the kinematic
variables are included.

Using the jet images and kinematic data defined in the previous sections, we
train neural networks to differentiate the signal from its SM backgrounds. We
have implemented a number of different networks which are described in detail
and compared in Appendix A. For some of the networks we only use part of the
available data to assess which data sets yield the highest discriminatory
power:

  * •

$\mathcal{K}$: Fully connected layer using only the kinematic data.

  * •

$\mathcal{I}_{C\ell}/\mathcal{I}_{CN\ell}$: Convolutional neural network
LeCun:2012  using charged (and neutral) hadron images and lepton images.

  * •

$\mathcal{I}_{C(N)\ell}+\mathcal{K}$: CNN using charged (and neutral) hadron
images and lepton images, combined with a FC using the kinematic data. The
flattened part of the CNN and the kinematic FC are interfaced to produce a
combined output.

The schematic structure of the best-performing network
$\mathcal{I}_{CN\ell}+\mathcal{K}$ is illustrated in Fig. 4. The networks are
implemented in the PyTorch paszke2019pytorch  framework. Further details on
the network architectures are provided in Appendix A.

The data set consists of 130,000 events, which are split into training (64%),
validation (16%), and test (20%) sets. When training the network, the entire
training data is split into a smaller number of samples (referred to as mini-
batches) so that, in each iteration, the network is trained with a randomly
chosen mini-batch. It is iteratively trained with the Adam optimiser
Kingma:2014vow  to update the hyper-parameters (referred to as weights) of the
network in a way to minimise the cross entropy loss function.

The validation set is used to monitor the generalisation error and choose the
best model parameters as those with the lowest validation loss, and the test
set is exclusively used for final evaluation. All three data sets are
comprised to equal parts of signal and background events. The number of events
of the different backgrounds is proportional to the respective cross sections
after the preselection cuts. Unless stated otherwise, we train and validate
the network with a signal data set that contains several mass hypotheses:
$m_{S}=300$ GeV to 800 GeV in steps of 50 GeV, each appearing with equal
weight. For the test sets, only a single mass $m_{S}$ is used, $m_{S}=300$ GeV
to 1000 GeV in steps of 50 GeV.

While a full treatment of the systematic uncertainties is beyond the scope of
this article, we have control over the systematics that are due to the machine
learning: The limited number of Monte Carlo events makes the results
susceptible to statistical fluctuations. Another contribution is the
stochastic nature of the initialisation of the network parameters and the
gradient descent. In order to take these uncertainties into account, we
perform 20 independent training runs for each network and evaluate them
separately. In the following, the uncertainty bands indicate the $1\sigma$
fluctuations between the different runs.

##  5 Results

(a) $m_{S}=300$ GeV

(b) $m_{S}=550$ GeV

(c) $m_{S}=800$ GeV

Figure 5: Comparison of network performances with ROC curves. The markers
indicate the working points used in the following analysis.

We start by comparing the performance of the trained networks as classifiers.
To this end, the receiver operating characteristic (ROC) curves in Fig. 5 show
how many background events can be rejected for a given signal efficiency. As
expected, the performance increases from $\mathcal{K}$ to
$\mathcal{I}_{CN\ell}+\mathcal{K}$ as more data are used. All networks perform
better at higher signal masses $m_{S}$, but this effect is largest for the
kinematic data: $\mathcal{K}$ has the largest boost in performance, and only
for higher $m_{S}$ does the inclusion of kinematic data significantly improve
the image-based classification. This can be explained by the increased $p_{T}$
of the final state particles.

When calculating the $2\sigma$ exclusion bound Cowan:2010js , we require:

$Z_{\mathrm{exc}}\equiv\sqrt{-2\,\ln\bigg{(}\frac{L(S\\!+\\!B|B)}{L(B|B)}\bigg{)}}\geq
1.64,\;\;\;\;\;\text{with}\;\;\;L(x|n)=\frac{x^{n}}{n!}e^{-x}\,,$ (7)

where $L(x|n)$ is the likelihood of observing $n$ events when $x$ events were
expected, and $S$ and $B$ are the number of signal and background events,
respectively. For achieving a $5\sigma$ expected discovery reach, we require

$Z_{\mathrm{dis}}\equiv\sqrt{-2\,\ln\bigg{(}\frac{L(B|S\\!+\\!B)}{L(S\\!+\\!B|S\\!+\\!B)}\bigg{)}}\geq
5\,.$ (8)

When putting a cut on the neural network score (NN score), it is common
practice to choose the cut that maximizes Eq. 8 Huang:2022rne . However, our
networks are such strong classifiers, that this method regularly leads to cuts
with few or even less than 1 background event remaining, for which a
statistical treatment becomes challenging. We therefore conservatively choose
our NN score cut by demanding a fixed amount of background events. Details
which justify this approach are provided in Sec. B.2. Once the NN score cut is
fixed, we determine the $5\sigma/2\sigma$ cross sections by rescaling the
cross section iteratively until Eqs. 8 and 7 are satisfied.

(a) Discovery reach of $S^{++}S^{--}$

(b) Expected upper limit of $S^{++}S^{--}$

(c) Discovery reach of $S^{++}S^{--}$ and $S^{\pm\pm}S^{\mp}$

(d) Expected upper limit of $S^{++}S^{--}$ & $S^{\pm\pm}S^{\mp}$

Figure 6: Expected discovery reach (left) and exclusion limit (right) of
$S^{++}$ production at the HL-LHC (3 ab-1) for different network
architectures. In the first row simulated events of only $S^{++}S^{--}$ pair
production are taken into account. In the second row simulated events of
$S^{++}S^{--}$ and $S^{\pm\pm}S^{\mp}$ production with $m_{S^{++}}=m_{S^{+}}$
and the corresponding cross sections from Fig. 2 are taken into account. The
dashed lines indicate the reference cross sections in the SU(5)/SO(5) model of
$S^{++}S^{--}$ (top row) and the sum of the cross sections of $S^{++}S^{--}$
and $S^{\pm\pm}S^{\mp}$ (bottom row) at 14 TeV. For the solid lines, the
networks were trained on multiple masses simultaneously. The points marked by
stars were trained using only the respective masses.

The expected discovery reach and upper limits for the networks $\mathcal{K}$,
$\mathcal{I}_{C\ell}+\mathcal{K}$, and $\mathcal{I}_{CN\ell}+\mathcal{K}$ are
shown in Figs. 6(a) and 6(b), where the NN score cuts were placed to obtain
50, 5, and 5 background events after the cut, respectively, and we assume an
integrated luminosity of 3 ab-1. The corresponding signal efficiencies can be
read off from the markers in Fig. 5. The black dashed line indicates the 14
TeV pair production cross section in the SU(5)/SO(5) model. Using only
kinematic data, the signal process can just barely be discovered up to $440$
GeV and excluded up to $650$ GeV. If we combine them with charged jet and
lepton images, however, the signal can be discovered up to $550$ GeV and
excluded to $750$ GeV. If advances in pileup suppression also allow the usage
of neutral jet images, the discovery reach can be increased by another $90$
GeV while the projected exclusion is increased by $70$ GeV.

Discovery reach and expected bounds on $\sigma$ become significantly weaker
for low masses. This can be partially remedied by choosing a different
training data set. The solid lines were obtained by training on a combined
data set containing equal amounts of signal events with masses between 300 and
800 GeV in steps of 50 GeV. However, the stars in Figs. 6(a) and 6(b) show the
bound that can be reached if the networks are trained only on signal events
with the mass they are evaluated at. This has little impact for
$\mathcal{I}_{CN\ell}+\mathcal{K}$, but noticeably improves the bounds for
$\mathcal{K}$ and to a lesser extent also $\mathcal{I}_{C\ell}+\mathcal{K}$ at
$m_{S}=300$ GeV. At higher masses, training on a single mass is less
beneficial or even slightly detrimental compared to the combined data set. We
explore the impact of the training data set further in Sec. B.1.

In the second row of Fig. 6, also the singly charged scalar is taken into
account. The solid lines are now discovery reach (left) and exclusion limit
(right) for the combined $S^{++}S^{--}$ and $S^{\pm\pm}S^{\mp}$ production
with the relative cross sections shown in Fig. 2. The dashed line is the sum
of the two cross sections, so that an intersection indicates
discovery/exclusion of both processes. By including the second process, both
the discovery reach and the exclusion limit are increased by about 50 GeV for
all three networks.

Figure 7: Expected exclusion limit of $S^{++}S^{--}$ pair production at the
LHC with $\mathcal{L}_{\mathrm{int}}=139$ fb-1 for different network
architectures. The recast bounds are taken from Ref. Cacciapaglia:2022bax .
The black line indicates the 13 TeV reference cross sections in the
SU(5)/SO(5) model.

Having now seen the bounds from the neural network approach, it is natural to
ask how this compares to conventional cut-and-count methods. No direct search
for the present final state has been performed yet, but in Ref.
Cacciapaglia:2022bax  an upper limit on the cross section is derived from the
recast of an ATLAS search for $R$-parity violating supersymmetry using the
full Run-2 data set ATLAS:2021fbt . The recast results are compared with the
bounds derived from the neural networks with 139 fb-1 in Fig. 7. Naturally, a
dedicated search would provide stronger bounds than the recast, so we cannot
draw quantitative conclusions from Fig. 7. However, it is evident that the
performance of the network based on purely kinematic information is comparable
to the recast of ATLAS:2021fbt  for low $m_{S}$. The signal regions which
yield the dominant recast bound demand same-sign leptons and multiple
$b$-tagged jets and are thus similar to our preselection criteria, and the
selection criteria of ATLAS:2021fbt  are based purely on kinematic
information. Inclusion of jet image data of charged and potentially neutral
jet images processed by a CNN yield a substantial improvement in
discrimination power in our study.

##  6 Conclusion and outlook

The search for new physics is and remains one of the main aims of the LHC.
Composite Higgs models with extended scalar sectors are theoretically well
motivated, their signatures are far from being ruled out, and the LHC can
explore untested parameter space at its high luminosity run if dedicated
searches are performed. In this article, we proposed a search strategy for a
doubly charged scalar $S^{++}$ which is pair produced and which decays into
$W^{+}t\bar{b}$. An $S^{++}$ with this decay mode is predicted – for example –
by composite Higgs models with an underlying fermionic description based on
the global symmetry breaking pattern SU(5)$\rightarrow$SO(5).

In our proposed search, we focus on the process $pp\rightarrow
S^{++}S^{--}\rightarrow W^{+}t\bar{b}W^{-}\bar{t}b$. This process is not
constrained by current LHC searches Cacciapaglia:2022bax  and our proposal is
aiming to close this gap. We target the same-sign lepton final state of this
channel. We implement deep neural networks which aim to optimise the
discrimination power between signal and SM backgrounds based on
($\mathcal{K}$) kinematic information of reconstructed jets and leptons;
($\mathcal{K}+\mathcal{I}_{C\ell}$) additional information from jet images of
charged jets and leptons; ($\mathcal{K}+\mathcal{I}_{CN\ell}$) yet additional
information from jet images of neutral jets, as defined in detail in Sec. 4.
As a main result (shown in Fig. 6) we obtain projected discovery and exclusion
bounds at the high luminosity run of the LHC for (A) doubly charged scalars in
the process $pp\rightarrow S^{++}S^{--}\rightarrow W^{+}t\bar{b}W^{-}\bar{t}b$
and (B) models with doubly and singly charged scalars, where the latter is
contributing to the signal region via the process $pp\rightarrow
S^{++}S^{-}\rightarrow W^{+}t\bar{b}\bar{t}b$ and its conjugate process. We
find that the use of convolutional neural networks on charged jet and lepton
images substantially increases the discovery and exclusion potential, with
projected discovery reach of $m_{S}\leq 550$ GeV and $m_{S}\leq 600$ GeV for
the scenarios (A) and (B) and a projected exclusion potential up to $m_{S}\leq
750$ GeV and $m_{S}\leq 800$ GeV, respectively. A further improvement can be
achieved when including neutral jet images into the data set, which can raise
projected discovery reaches to $m_{S}\leq 640$ GeV and $m_{S}\leq 690$ GeV and
projected exclusion potential to $m_{S}\leq 820$ GeV and $m_{S}\leq 880$ GeV.
However, this latter improvement relies on sufficiently controlling pileup
effects. We do not study pileup effects in this article and instead present
results for $\mathcal{K}+\mathcal{I}_{C\ell}$ data as conservative and results
for $\mathcal{K}+\mathcal{I}_{CN\ell}$ data as a potentially realisable but
optimistic scenario.

For this study, we implemented a variety of neural networks which are
discussed and compared in Appendix A to which we refer especially readers
interested in the machine learning aspects of this proposed BSM search.
Appendix B contains additional details related to the choice of the cut on the
NN score and the composition of the training data set. We consider the
information therein methodologically relevant not only for the study presented
here, but more generally for BSM searches with unknown mass of the BSM
particle and with potentially very low SM backgrounds after NN score cuts.

This study demonstrates that the use of machine learning techniques in
general, and the use of jet images and CNNs in particular, provides a large
potential to improve the sensitivity of LHC searches for BSM final states with
a lot of hadronic activity. The channel $pp\rightarrow S^{++}S^{--}\rightarrow
W^{+}t\bar{b}W^{-}\bar{t}b$ in the same-sign lepton final state provides one
example. We expect that similar improvements are achievable for other final
states like Drell-Yan pair produced BSM scalars, pair produced vector-like
quarks with exotic decays, or signals from R-parity violating supersymmetry to
just name a few examples which deserve further investigation.

## Acknowledgements

We thank Raimund Ströhmer for discussions. This work has been supported by the
international cooperation program “GEnKO” managed by the National Research
Foundation of Korea (No. 2022K2A9A2A15000153, FY2022) and DAAD, P33 - projekt-
id 57608518. This work is supported by the Center for Advanced Computation at
Korea Institute for Advanced Study. T.F. is supported by a KIAS Individual
Grant (AP083701) via the Center for AI and Natural Sciences at the Korea
Institute for Advanced Study. P.K. is supported in part by a KIAS Individual
Grant No. PG021403, and by National Research Foundation of Korea (NRF)
Research Grant NRF- 2019R1A2C3005009. M.K. is supported by the
“Studienstiftung des deutschen Volkes”. J.K. is supported by the National
Research Foundation of Korea (NRF) grant funded by the Korea government (MSIT)
(No. 2021R1C1C1005076). J.K. acknowledges the hospitality at APCTP where part
of this work was done.

##  Appendix A Neural network architectures

Machine learning can bring significant improvements to results, but it is
important to understand how to effectively utilise them. For this work, we
tested various types of data representations and neural network architectures
and found that a combination of kinematic variables together with low level
calorimeter information represented in terms of jet images yields the best
results. In this section, we describe several neural network architectures we
have attempted and evaluate the performance of each network by analysing the
ROC curves. This analysis uses only a data set with the mass of the signal
$m_{S}=400$ GeV, and we define two additional kinematic variables.

$\displaystyle\mathcal{K}_{p_{\mu}}^{\rm(vis)}=\\{p_{\mu}(\ell_{1}),p_{\mu}(\ell_{2}),p_{\mu}(b_{1}),p_{\mu}(b_{2}),p_{\mu}(b_{3}),p_{\mu}(j_{1}),p_{\mu}(j_{2}),p_{\mu}(j_{3})\\}\,,$
(9) $\displaystyle\mathcal{K}_{2}=\bigcup_{i\neq j}M_{ij}\cup\bigcup_{i\neq
j}\Delta R_{ij}\;,$ (10)

where $\ell_{1}$, $\ell_{2}$, $b_{i}$, and $j_{i}$ denote two same-sign
leptons, $i$-th leading $b$-tagged jets, and $i$-th leading non-$b$-tagged
jets respectively.

###  A.1 Fully connected deep neural networks

Figure 8: A schematic architecture of the fully-connected NN (FC) used in this
paper.

Fully connected layer (FC) neural networks are simple structures that can
easily be used when an input layer is composed of the four-momenta of
reconstructed particles or combinations of kinematic variables. As shown in
Fig. 8, FCs are characterised by all neurons in each layer being connected to
all neurons in the next layer. In the FC we implemented, there are 8 hidden
layers, with the number of neurons decreasing from 2400 to 300, and a last
hidden layer is connected to an output layer. The output layer consists of two
neurons representing background or signal, respectively. The values associated
with these two neurons are denoted as $\hat{p}_{k}$, where $k$ is 0 and 1 for
a background and a signal, respectively. We normalise the values $\hat{p}_{k}$
as $p_{k}=e^{\hat{p}_{k}}/\sum e^{\hat{p}_{k}}$ using a softmax function,
where $p_{0}$ represents the “probability to be a background” that the network
assigns to this event, and analogous with $p_{1}$ and signal. We use the
Rectified Linear Unit (ReLU) pmlr-v15-glorot11a  activation function between
hidden layers. When training the network, we use a mini-batch size of 20 and a
learning rate of $5\times 10^{-7}$ of the Adam optimiser to minimise the cross
entropy loss function

$\displaystyle L=-y_{k}\log p_{k}-(1-y_{k})\log(1-p_{k})\;,$ (11)

where a background and a signal event are labeled as $y_{0}=0$ and $y_{1}=1$,
respectively. We include the $L_{2}$ regularisation222 The $L_{2}$
regularisation shifts the loss function by $L\rightarrow
L+\frac{1}{2}\lambda\|\textbf{W}\|^{2}$ where W represents all weights, and
$\lambda$ denotes weight_decay. In this paper, we have applied the $L_{2}$
regularisation to all neural networks. terms to the loss function with
weight_decay=$2\times 10^{-9}$ to alleviate the problem of data overfitting.
The learning rate and weight_decay are the same for the cases of
$\mathcal{K}_{p_{\mu}}^{\rm vis}$, $K_{p_{\mu}}^{\rm vis}+\mathcal{K}_{2}$,
and $\mathcal{K}$. In addition, we use a validation set to make sure that the
neural network is not overfitted. We train the model for 100 epochs. If the
validation loss does not reach a new minimum within 50 epochs, the training is
stopped and the parameters of the epoch with the minimal validation loss are
saved333 Unless otherwise stated, the activation function, the type of
optimiser, the loss function, the configuration of an output layer, and the
use of a validation set are the same for all other neural networks. .

###  A.2 Convolutional neural networks

Figure 9: A schematic CNN architecture used in this article. The separate FC
chain in the right-upper corner is used only when kinematic variables are
included.

CNNs are a class of machine learning algorithms commonly used for image
recognition tasks, and they are useful when a final state is represented as a
set of images. An input is a 3D image of $\mathcal{I}_{CN\ell}$
($\mathcal{I}_{C\ell}$) whose dimension is given by $3\times 50\times 50$
($2\times 50\times 50$) where 3 (2) denotes the number of image layers in Eq.
4 (Eq. 5).

To fully exploit spatial correlations among different layers, it is necessary
to utilise a 3D convolution. For this operation, we need to add a singleton
dimension to reshape the image dimension as $1\times 3\times 50\times 50$
($1\times 2\times 50\times 50$), where each entry is identified as a channel,
a depth, a height, and a width respectively. Then, we apply the 3D convolution
with kernel size of $3\times 3\times 3$ ($2\times 3\times 3$), stride of 1,
padding size of $(0,1,1)$, and 4 feature maps. After the 3D convolution, we
apply the batch normalisation and the ReLU activation function, and the
resulting image dimension becomes $4\times 1\times 50\times 50$. To remove the
singleton ($i.e.$ redundant) dimension in the second entry, we reshape the
image size again to make it as $4\times 50\times 50$.

Next, we apply a 2D convolution twice using the kernel size of $3\times 3$,
the stride of 1, the padding size of 1, and 4 feature maps. Throughout this
paper, for all 2D convolutions, padding size is fixed to 1, and after each 2D
convolution we apply the batch normalisation and the ReLU activation function.

The pooling method used in CNNs takes the most representative pixel in each
region of a feature map. For example, a max (average) pooling operation keeps
the maximum (average) value within a feature map, while discarding the other
remaining information. This is how most CNNs reduce a dimension of an input
image data through a series of combined convolutional and pooling operations
to learn important characteristics of a given image in various scales. As
shown in Fig. 9, we design a two-pronged architecture to include both the max
and the average pooling operations with a kernel size of $5\times 5$ and a
stride of 5, which reduces the image dimension down to $4\times 10\times 10$.
Each branch is followed by the 2D convolution and the pooling again, until its
size is reduced to $32\times 2\times 2$. We flatten out the final image into a
one-dimensional strip to combine it with 3 hidden layers with 600 neurons
each.

A neural network that combines kinematic variables with a CNN architecture has
been demonstrated in Ref. Huang:2022rne . In a similar way, we construct a
separate FC with 6 hidden layers, with the number of neurons decreasing from
1200 to 600 (as shown in the upper-right corner of the Fig. 9). The last FC
layer with 600 neurons is combined with the last hidden layer of the CNN. The
combined layer is followed by 4 hidden layers with the number of neurons
decreasing from 600 to 300. The last hidden layer is then connected to the
output layer. The network is trained using a mini-batch size of 20. The
learning rate is set to $5\times 10^{-4}$ for $\mathcal{I}_{C\ell}$ and
$\mathcal{I}_{CN\ell}+\mathcal{K}$, and $3\times 10^{-4}$ for
$\mathcal{I}_{CN\ell}$ and $\mathcal{I}_{C\ell}+\mathcal{K}$. We use
weight_decay=$2\times 10^{-4}$ for $\mathcal{I}_{C\ell}+\mathcal{K}$ and
$3\times 10^{-4}$ for the remaining networks.

###  A.3 Capsule neural networks

Figure 10: A schematic architecture of the CapsNet used in this paper.

As explained in the previous subsection, the pooling operation in CNNs is the
mechanism to reduce the image dimension to learn important features in various
scales. However, a major problem associated with this pooling operation is
that it loses information by throwing away other pixels. Eventually, it
accompanies the loss of spatial information about where things are. Therefore,
CNNs are agnostic to the geometric correlations among the pixels at higher and
lower levels. A capsule neural network (CapsNet) NIPS2017_2cad8fa4 ;
Diefenbacher:2019ezd  was proposed to address this problem by introducing a
concept of dynamic routing, which compensates for the shortcomings of the
pooling operations. For dedicated readers, we refer to Ref. Huang:2022rne  for
more detailed descriptions of the CapsNet and its variants. Although it was
originally proposed to improve CNNs, the previous study Huang:2022rne
demonstrated that the performance of the CapsNet was worse than CNNs in
general. We suppose that it would have been attributed to sparse jet images in
their example of the double Higgs production $(h\to b{\bar{b}})(h\to
W^{\pm}W^{*\mp}\to\ell^{+}\nu_{\ell}{\ell^{\prime}}^{-}\bar{\nu}_{\ell^{\prime}})$.
On the other hand, in our signal process, the busier final state delivers
denser jet images with more activated pixels. Ample features would appear in
jet images that the CapsNet might take advantage of. This motivates us to
apply the CapsNet to our image data as a comparison with CNNs to see how much
we gain from using a different type of image-based neural network.

Given a 3D image data set $\mathcal{I}_{CN\ell}$ ($\mathcal{I}_{C\ell}$), we
first apply a 3D convolution with a kernel size of $3\times 3\times 3$
($2\times 3\times 3$), a stride of 1, a padding of 0, and out channels of 256.
Then, we apply a 2D convolution with a kernel size of $3\times 3$, a stride of
2, and a padding of 1 until the image dimension is reduced down to $256\times
6\times 6$. The output neurons of the 2D convolution are reshaped to form
$6\times 6\times 32=1152$ primary capsule vectors of length 8, which contain
lower-level information of the input image. These primary capsules are
connected to two 16-dimensional higher-level capsules (digital capsule)
through a dynamic routing process NIPS2017_2cad8fa4 . The routed digital
capsule vectors $\textbf{v}_{j}$ are used to define the margin loss function,
where $j$ denotes a class label, 0 or 1

$\displaystyle L_{j}=$ $\displaystyle
T_{j}\,\text{max}\Big{(}0,m^{+}-||\textbf{v}_{j}||\Big{)}^{2}+\lambda(1-T_{j})\,\text{max}\Big{(}0,||\textbf{v}_{j}||-m^{-}\Big{)}^{2}\;,$
(12)

where a background and a signal event are labeled as $T_{0}=0$ and $T_{1}=1$,
$m^{+}=0.9$, $m^{-}=0.1$, and $\lambda=0.5$. In addition, the digital capsule
vectors are connected to 3 hidden layers with increasing number of neurons
from 512 to 7500 (5000) via a parallel routine to reconstruct the input image.
A reconstruction loss function is defined by the sum of squared differences in
pixel intensities between the reconstructed
($\mathcal{I}^{\text{(reco)}}_{i}$) and the input
($\mathcal{I}^{\text{(input)}}_{i}$)

$\displaystyle
L_{\text{reco}}=\frac{1}{N}\sum_{i=1}(\mathcal{I}^{(\text{reco})}_{i}-\mathcal{I}^{\text{(input)}}_{i})^{2}\;,$
(13)

where the index $i$ runs from 1 to the total number of pixels in the image,
and $N$ is a normalisation factor defined by a total number of pixels in the
image data times the total number of training events. The total loss function
is given by

$\displaystyle L=L_{j}+\alpha L_{\text{reco}}\,,$ (14)

where we multiply the reconstructed loss function by a scaling factor of
$\alpha=5.0\times 10^{-4}$ to prevent it from significantly exceeding the
margin loss. The network is trained using a mini-batch size of 20. The
learning rate is set to $5\times 10^{-3}$ for $\mathcal{I}_{C\ell}$ and
$6\times 10^{-3}$ for $\mathcal{I}_{CN\ell}$. We use weight_decay=$2.0\times
10^{-7}$ in both networks.

###  A.4 Graph neural networks

Figure 11:  A schematic architecture of the GNN used in this paper.

Graph neural networks (GNN) 1555942 ; 4700287  are specialised to learn graph
structured information from the data that is comprised of nodes and edges
connecting them. In this framework, each reconstructed object ($i.e.$ jets and
leptons) can be represented as a node, and each node $i$ has a feature vector
$\textbf{x}_{i}$. Similarly, an edge between two nodes $i$ and $j$ is
represented by a vector $\textbf{e}_{ij}$ which contains an angular distance
and an invariant mass of two particles.

Among various realisations of GNN models, we employ the Message Passing Neural
Network (MPNN) DBLP:journals/corr/GilmerSRVD17 . In this implementation, each
node $i$ has a 4-dimensional feature vector $\textbf{x}^{0}_{i}=(p_{x},\
p_{y},\ p_{z},\ E)_{i}$. Each edge vector $\textbf{e}_{ij}$ contains an
angular distance and an invariant mass of two particles $i$ and $j$. A
prerequisite step in the MPNN algorithm is to multiply each feature vector
$\textbf{x}^{0}_{i}$ by a weight matrix $\textbf{W}^{0}$

$\displaystyle\textbf{m}^{0}_{i}=\textbf{W}^{0}\textbf{x}^{0}_{i}\;.$ (15)

The feature vector evolves with a series of multiplications with nearby edge
vectors and weight matrices, and its recurrent relation between layers $n$ and
$n-1$ is given by

$\displaystyle\textbf{m}_{i}^{n}=\sum_{j\in\mathcal{E}}\textbf{W}^{n}\Big{(}\textbf{m}_{i}^{n-1}\oplus\big{(}\textbf{W}^{{}^{\prime}n}(\textbf{m}_{j}^{n-1}\oplus\textbf{e}_{ij})\big{)}\Big{)},$
(16)

where $\textbf{W}^{n}$ and $\textbf{W}^{{}^{\prime}n}$ denote weight matrices,
and $\mathcal{E}$ denotes a set of nodes that are connected with a node $i$.
This process, referred to as a graph convolution, is repeated three times.
After the last graph convolution, each feature vector is concatenated in one
dimension and then multiplied by a weight matrix to become a vector of length
2

$\displaystyle\hat{\textbf{p}}_{k}=\textbf{W}^{N}\textbf{m}^{N}\;,$ (17)

where $k$ denotes the class label, 0 or 1 for a background or a signal,
respectively. As previously mentioned, the configuration of the activation
function and the loss function is the same as other networks except CapsNet.
The network is trained using a mini-batch size of 20, the learning rate is set
to $7.0\times 10^{-7}$, and weight_decay=$2.0\times 10^{-4}$.

###  A.5 Comparison of different networks

Figure 12: ROC curves for selected NN architectures evaluated on the same test
sample.

We present ROC curves in Fig. 12 for several architectures described above
where $\epsilon_{S}$ and $1/\epsilon_{B}$ denote the signal efficiency and the
background rejection, respectively. The area under the ROC curve (AUC) is a
useful measure for evaluating the performance of each architecture, where the
closer the AUC is to 1, the better the network is able to classify signal and
backgrounds.

We find that the CNN with $\mathcal{I}_{CN\ell}+\mathcal{K}$ has the highest
AUC (0.964), followed by the CNN with $\mathcal{I}_{CN\ell}$ (0.951), the
CapsNet with $\mathcal{I}_{CN\ell}$ (0.934), the FC with $\mathcal{K}$
(0.866), the FC with $\mathcal{K}_{p_{\mu}}^{\rm vis}+\mathcal{K}_{2}$
(0.852), and the GNN with $\mathcal{K}_{p_{\mu}}^{\rm vis}+\mathcal{K}_{2}$
(0.849). This indicates that the neural networks cannot efficiently learn
important features of the signal and backgrounds when they are trained only
with four-momenta information. Even when the engineered kinematic variables
$\mathcal{K}$ are used, the results can be only mildly improved. On the other
hand, the image-based architectures generally outperform the networks that are
based on the kinematic variables only. This implies that the use of image
recognition techniques in the busy final states comprised of multiple jets and
leptons provides a powerful tool in the search for new particles. As also can
be seen, combining the image-based CNN with the FC acting on kinematic
variables provides a further gain in performance.

Due to this comparison, we choose $\mathcal{I}_{CN\ell}+\mathcal{K}$ as data
representation and a combination of a CNN acting on $\mathcal{I}_{CN\ell}$ and
an FC acting on $\mathcal{K}$ in the main article. In the main article, the
additional data set $\mathcal{I}_{C\ell}+\mathcal{K}$ is used and
discussed.444Comparison between CNN and CapsNet acting on
$\mathcal{I}_{C\ell}$ (not shown in Fig. 12) show that also for this data set,
CNN outperforms CapsNet. This data set is not introduced for technical or
performance reasons but because data including neutral jet images is
potentially subject to corrections and uncertainties resulting from pileup
effects, i.e. from the problem of separating the detector response from
several particle collisions which happen simultaneously in the same bunch
crossing. For charged jet images, the tracks of charged particles can be very
efficiently used in order to subtract the pileup such that our projections
based on $\mathcal{I}_{C\ell}+\mathcal{K}$ are conservative. This method has
limitations for removing neutral particles, but there are on-going
experimental works to cope with the problem CMS:2020ebo ; Bertolini:2014bba ,
which are beyond the scope of this article. Thus, by assuming that the pileup
can be controlled, the analysis based on $\mathcal{I}_{CN\ell}+\mathcal{K}$ is
an optimistic (best case) scenario.

##  Appendix B Neural networks: Training and evaluation details

###  B.1 Choice of training data sets and comparisons

As described in Sec. 4.3 the data set is split into training (64%), validation
(16%), and test (20%) sets. For training and validation data sets, we use data
comprised of background events (50%) in which each background is added
according to its respective cross section after the preselection cuts, and of
signal events (50%) for which we add equal numbers of signal events simulated
with $m_{S}=300$ GeV to 800 GeV in steps of 50 GeV. The network is thus not
trained on a specific BSM scalar mass, but instead the training aims to
identify common features in signal events which distinguish the BSM events
from background and apply to signal events from a broad mass window. We choose
this approach because of several reasons: First, the signal contains missing
energy (because we demand same-sign leptons and thus leptonic decays of a top
and a $W$ boson) and a mass reconstruction is not as straight forward as in a
2-body resonance search. Second, training is expensive on computing resources,
and as we show here, results from this collective training yield comparable
results to re-training for a series of mass hypotheses. Third and last,
training with a mixed set of different masses allows to straight forwardly
apply efficiencies and bounds for all masses within the interval of
$m_{S}=300$ GeV to 800 GeV and it gives confidence that extrapolations above
800 GeV are feasible.555For a series of trainings at different benchmark
masses, one could determine efficiencies by evaluating each trained network on
each test data set with different mass signal events which can then be
combined, but training costs and evaluation complexity are much higher in such
an approach.

The network is trained and validated on this “mixed” data set. The 15 test
data sets are each specific to a single mass $m_{S}$ of $m_{S}=300$ GeV to
1000 GeV in steps of 50 GeV and comprise of (50%) background events and (50%)
signal events simulated at mass $m_{S}$. To obtain projections of discovery
reach and exclusion bounds at a given mass $m_{S}$, the fully trained network
is applied to the test data set of the respective mass. The NN score is chosen
as detailed in the next section, yielding a certain number of signal and
background events, $S$ and $B$. We then replace $S\to\mu S$ to iteratively
rescale the signal until the target significance is reached, see Eqs. 8 and 7.

(a) Network $\mathcal{K}$

(b) Network $\mathcal{I}_{C\ell}+\mathcal{K}$

(c) Network $\mathcal{I}_{CN\ell}+\mathcal{K}$

Figure 13: Comparison of networks trained on a single fixed mass.

For the main results in Figs. 5, 6 and 7, this training and evaluation
procedure is independently performed 20 times for each network and each mass
point, and the bands indicate the 1$\sigma$ fluctuations between the 20
different runs. For the supplementary results in Fig. 13, the bands are based
on 10 independent runs.

(a) Network $\mathcal{K}$

(b) Network $\mathcal{I}_{C\ell}+\mathcal{K}$

(c) Network $\mathcal{I}_{CN\ell}+\mathcal{K}$

Figure 14: Discovery reach as a function of number of background events, i.e.
the NN score cut.

To test how well the training on a mixed data set is performing, we choose
benchmark points at $m_{S}=300$ GeV, 550 GeV and 800 GeV for which we train
additional networks on a single mass. When each of these networks are
evaluated at the mass they were trained for, one would expect better
performance at this mass as compared to the performance of the network trained
on the mixed training set. In Figs. 6(a) and 6(b), the stars indicate the
results of the mass-specific trainings. At $m_{S}=300$ GeV, results from the
mass specific training are indeed significantly better. However, for 550 GeV
and 800 GeV, this effect is reduced or even lost, which shows that at higher
masses, the training with the mixed training data set performs comparably
while saving computation resources. To present some more detailed results on
the mass specific training, Fig. 13 shows the discovery reach of the networks
trained at a fixed mass $m_{S}=300$ GeV, 550 GeV, and 800 GeV as compared to
the discovery reach of the networks trained on the mixed mass data set. At
$m_{S}=300$ GeV, the network trained on 300 GeV events performs best, but it
gets quickly overtaken by the network trained on 550 GeV in the medium mass
range. Only above 800 GeV is there an advantage to training on high mass
events. Even then, the improvement over the $550$ GeV network is only
marginal.

Overall, the mixed data set performs well when compared to networks trained on
a single mass. Only in the low mass region is there a significant difference.
Therefore, a hybrid approach may be best to maximise the sensitivity of the
search. Two separate networks could be trained: one for the region $m_{S}\leq
400$ GeV, and one for $m_{S}\geq 400$ GeV. The latter can either be trained on
a mixed to data set or a single mass from the region of interest, e.g. 550
GeV.

###  B.2 Cut on the neural network score

When evaluating a trained network on the test set, we obtain a NN score
distribution, i.e. each event is assigned a score between 0 (background-like)
and 1 (signal-like). In order to derive physical results from this
distribution, we need to place a cut on the NN score, below which events are
discarded. To maximise the sensitivity of our search, the cut aught to be
placed such as to maximise the signal to background ratio. This is typically
enforced by choosing the cut to maximise the discovery significance Eq. 8.
However, we found that this method often places the cut very close to a NN
score of 1, resulting in very few background events. Designing a search with
an expected number of background events of $<1$ raises questions about the
statistical soundness of the result.

To address this challenge, we plot the discovery reach as a function of the
number of background events passing the NN score cut, see Fig. 14.
Fortunately, this reveals that the 5$\sigma$ cross section reach is not very
sensitive to the number of background events. In other words: For the specific
signature and the networks we study, decreasing the NN cut increases the
number of signal and background events in such a way that there is only a mild
change in discovery reach. We therefore fix the NN score by demanding a
certain number of background events to pass the cut. For the networks
$\mathcal{K}$, $\mathcal{I}_{C\ell}+\mathcal{K}$, and
$\mathcal{I}_{CN\ell}+\mathcal{K}$, we require 50, 5, and 5, background
events, respectively. We fully admit that these numbers are chosen ad hoc.
They give sufficient statistics and at the same time do not grossly degrade
the discovery potential.

## References

  * (1) ATLAS collaboration, G. Aad et al., _Observation of a new particle in the search for the Standard Model Higgs boson with the ATLAS detector at the LHC_ , _Phys. Lett. B_ 716 (2012) 1–29, [1207.7214].
  * (2) CMS collaboration, S. Chatrchyan et al., _Observation of a New Boson at a Mass of 125 GeV with the CMS Experiment at the LHC_ , _Phys. Lett. B_ 716 (2012) 30–61, [1207.7235].
  * (3) D. B. Kaplan and H. Georgi, _SU(2) x U(1) Breaking by Vacuum Misalignment_ , _Phys. Lett. B_ 136 (1984) 183–186.
  * (4) D. B. Kaplan, H. Georgi and S. Dimopoulos, _Composite Higgs Scalars_ , _Phys. Lett. B_ 136 (1984) 187–190.
  * (5) R. Contino, Y. Nomura and A. Pomarol, _Higgs as a holographic pseudoGoldstone boson_ , _Nucl. Phys. B_ 671 (2003) 148–174, [hep-ph/0306259].
  * (6) K. Agashe, R. Contino and A. Pomarol, _The Minimal composite Higgs model_ , _Nucl. Phys. B_ 719 (2005) 165–187, [hep-ph/0412089].
  * (7) D. B. Kaplan, _Flavor at SSC energies: A New mechanism for dynamically generated fermion masses_ , _Nucl. Phys. B_ 365 (1991) 259–278.
  * (8) K. Agashe and R. Contino, _The Minimal composite Higgs model and electroweak precision tests_ , _Nucl. Phys. B_ 742 (2006) 59–85, [hep-ph/0510164].
  * (9) R. Contino, _The Higgs as a Composite Nambu-Goldstone Boson_ , in _Theoretical Advanced Study Institute in Elementary Particle Physics: Physics of the Large and the Small_ , pp. 235–306, 2011, 1005.4269, DOI.
  * (10) B. Bellazzini, C. Csáki and J. Serra, _Composite Higgses_ , _Eur. Phys. J. C_ 74 (2014) 2766, [1401.2457].
  * (11) G. Panico and A. Wulzer, _The Composite Nambu-Goldstone Higgs_ , vol. 913.  Springer, 2016, 10.1007/978-3-319-22617-0.
  * (12) G. Cacciapaglia and F. Sannino, _Fundamental Composite (Goldstone) Higgs Dynamics_ , _JHEP_ 04 (2014) 111, [1402.0233].
  * (13) G. Cacciapaglia, C. Pica and F. Sannino, _Fundamental Composite Dynamics: A Review_ , _Phys. Rept._ 877 (2020) 1–70, [2002.04914].
  * (14) M. J. Dugan, H. Georgi and D. B. Kaplan, _Anatomy of a Composite Higgs Model_ , _Nucl. Phys. B_ 254 (1985) 299–326.
  * (15) J. Galloway, J. A. Evans, M. A. Luty and R. A. Tacchi, _Minimal Conformal Technicolor and Precision Electroweak Tests_ , _JHEP_ 10 (2010) 086, [1001.1361].
  * (16) G. Ferretti and D. Karateev, _Fermionic UV completions of Composite Higgs models_ , _JHEP_ 03 (2014) 077, [1312.5330].
  * (17) T. Ma and G. Cacciapaglia, _Fundamental Composite 2HDM: SU(N) with 4 flavours_ , _JHEP_ 03 (2016) 211, [1508.07014].
  * (18) A. Belyaev, G. Cacciapaglia, H. Cai, G. Ferretti, T. Flacke, A. Parolini et al., _Di-boson signatures as Standard Candles for Partial Compositeness_ , _JHEP_ 01 (2017) 094, [1610.06591].
  * (19) A. Agugliaro, G. Cacciapaglia, A. Deandrea and S. De Curtis, _Vacuum misalignment and pattern of scalar masses in the SU(5)/SO(5) composite Higgs model_ , _JHEP_ 02 (2019) 089, [1808.10175].
  * (20) A. Banerjee, D. B. Franzosi and G. Ferretti, _Modelling vector-like quarks in partial compositeness framework_ , _JHEP_ 03 (2022) 200, [2202.00037].
  * (21) A. Banerjee et al., _Phenomenological aspects of composite Higgs scenarios: exotic scalars and vector-like quarks_ , 2203.07270.
  * (22) G. Cacciapaglia, T. Flacke, M. Kunkel, W. Porod and L. Schwarze, _Exploring extended Higgs sectors via pair production at the LHC_ , _JHEP_ 12 (2022) 087, [2210.01826].
  * (23) Y. Lecun, L. Bottou, Y. Bengio and P. Haffner, _Gradient-based learning applied to document recognition_ , _Proceedings of the IEEE_ 86 (1998) 2278–2324.
  * (24) A. Krizhevsky, I. Sutskever and G. E. Hinton, _Imagenet classification with deep convolutional neural networks_ , in _Advances in Neural Information Processing Systems_ (F. Pereira, C. Burges, L. Bottou and K. Weinberger, eds.), vol. 25, Curran Associates, Inc., 2012, https://proceedings.neurips.cc/paper_files/paper/2012/file/c399862d3b9d6b76c8436e924a68c45b-Paper.pdf.
  * (25) W. Bhimji, S. A. Farrell, T. Kurth, M. Paganini, Prabhat and E. Racah, _Deep Neural Networks for Physics Analysis on low-level whole-detector data at the LHC_ , _J. Phys. Conf. Ser._ 1085 (2018) 042034, [1711.03573].
  * (26) CMS collaboration, M. Andrews, M. Paulini, S. Gleyzer and B. Poczos, _Exploring End-to-end Deep Learning Applications for Event Classification at CMS_ , _EPJ Web Conf._ 214 (2019) 06031.
  * (27) L. Huang, S.-b. Kang, J. H. Kim, K. Kong and J. S. Pi, _Portraying double Higgs at the Large Hadron Collider II_ , _JHEP_ 08 (2022) 114, [2203.11951].
  * (28) J. H. Kim, M. Kim, K. Kong, K. T. Matchev and M. Park, _Portraying Double Higgs at the Large Hadron Collider_ , _JHEP_ 09 (2019) 047, [1904.08549].
  * (29) J. Kim et al., _Large-Scale Deep Learning for Multi-Jet Event Classification_ , 2207.11710.
  * (30) S. Diefenbacher, H. Frost, G. Kasieczka, T. Plehn and J. M. Thompson, _CapsNets Continuing the Convolutional Quest_ , _SciPost Phys._ 8 (2020) 023, [1906.11265].
  * (31) T. Kim, J. H. Kim, S. Kumar, A. Martin, M. Münchmeyer and Y. Tsai, _Probing Cosmological Particle Production and Pairwise Hotspots with Deep Neural Networks_ , 2303.08869.
  * (32) L. de Oliveira, M. Kagan, L. Mackey, B. Nachman and A. Schwartzman, _Jet-images — deep learning edition_ , _JHEP_ 07 (2016) 069, [1511.05190].
  * (33) J. Cogan, M. Kagan, E. Strauss and A. Schwarztman, _Jet-Images: Computer Vision Inspired Techniques for Jet Tagging_ , _JHEP_ 02 (2015) 118, [1407.5675].
  * (34) A. Butter et al., _The Machine Learning landscape of top taggers_ , _SciPost Phys._ 7 (2019) 014, [1902.09914].
  * (35) J. Hajer, Y.-Y. Li, T. Liu and H. Wang, _Novelty Detection Meets Collider Physics_ , _Phys. Rev. D_ 101 (2020) 076015, [1807.10261].
  * (36) M. Farina, Y. Nakai and D. Shih, _Searching for New Physics with Deep Autoencoders_ , _Phys. Rev. D_ 101 (2020) 075021, [1808.08992].
  * (37) T. Heimel, G. Kasieczka, T. Plehn and J. M. Thompson, _QCD or What?_ , _SciPost Phys._ 6 (2019) 030, [1808.08979].
  * (38) A. Blance, M. Spannowsky and P. Waite, _Adversarially-trained autoencoders for robust unsupervised new physics searches_ , _JHEP_ 10 (2019) 047, [1905.10384].
  * (39) M. Crispim Romão, N. F. Castro and R. Pedro, _Finding New Physics without learning about it: Anomaly Detection as a tool for Searches at Colliders_ , _Eur. Phys. J. C_ 81 (2021) 27, [2006.05432].
  * (40) B. Bortolato, A. Smolkovič, B. M. Dillon and J. F. Kamenik, _Bump hunting in latent space_ , _Phys. Rev. D_ 105 (2022) 115009, [2103.06595].
  * (41) A. Hallin, J. Isaacson, G. Kasieczka, C. Krause, B. Nachman, T. Quadfasel et al., _Classifying anomalies through outer density estimation_ , _Phys. Rev. D_ 106 (2022) 055006, [2109.00546].
  * (42) P. Jawahar, T. Aarrestad, N. Chernyavskaya, M. Pierini, K. A. Wozniak, J. Ngadiuba et al., _Improving Variational Autoencoders for New Physics Detection at the LHC With Normalizing Flows_ , _Front. Big Data_ 5 (2022) 803685, [2110.08508].
  * (43) V. S. Ngairangbam, M. Spannowsky and M. Takeuchi, _Anomaly detection in high-energy physics using a quantum autoencoder_ , _Phys. Rev. D_ 105 (2022) 095004, [2112.04958].
  * (44) G. Karagiorgi, G. Kasieczka, S. Kravitz, B. Nachman and D. Shih, _Machine Learning in the Search for New Fundamental Physics_ , 2112.03769.
  * (45) T. Finke, M. Krämer, A. Morandini, A. Mück and I. Oleksiyuk, _Autoencoders for unsupervised anomaly detection in high energy physics_ , _JHEP_ 06 (2021) 161, [2104.09051].
  * (46) G. Ferretti, _Gauge theories of Partial Compositeness: Scenarios for Run-II of the LHC_ , _JHEP_ 06 (2016) 107, [1604.06467].
  * (47) J. Wess and B. Zumino, _Consequences of anomalous Ward identities_ , _Phys. Lett. B_ 37 (1971) 95–97.
  * (48) E. Witten, _Global Aspects of Current Algebra_ , _Nucl. Phys. B_ 223 (1983) 422–432.
  * (49) J. Schechter and J. W. F. Valle, _Neutrino Masses in SU(2) x U(1) Theories_ , _Phys. Rev. D_ 22 (1980) 2227.
  * (50) ATLAS collaboration, G. Aad et al., _Search for doubly and singly charged Higgs bosons decaying into vector bosons in multi-lepton final states with the ATLAS detector using proton-proton collisions at $\sqrt{\mathrm{s}}$ = 13 TeV_, _JHEP_ 06 (2021) 146, [2101.11961].
  * (51) ATLAS collaboration, B. Leban, _Search for doubly charged Higgs boson production in multi-lepton final states using 139 fb -1 of proton-proton collisions at $\sqrt{s}$ = 13 TeV with the ATLAS detector_, _PoS_ ICHEP2022 (2022) 1081.
  * (52) CMS collaboration, _A search for doubly-charged Higgs boson production in three and four lepton final states at $\sqrt{s}=13~{}\mathrm{TeV}$_, Tech. Rep. CMS-PAS-HIG-16-036, 2017.
  * (53) J. H. Kim and I. M. Lewis, _Loop Induced Single Top Partner Production and Decay at the LHC_ , _JHEP_ 05 (2018) 095, [1803.06351].
  * (54) J. H. Kim, S. D. Lane, H.-S. Lee, I. M. Lewis and M. Sullivan, _Searching for Dark Photons with Maverick Top Partners_ , _Phys. Rev. D_ 101 (2020) 035041, [1904.05893].
  * (55) A. Alloul, N. D. Christensen, C. Degrande, C. Duhr and B. Fuks, _FeynRules 2.0 - A complete toolbox for tree-level phenomenology_ , _Comput. Phys. Commun._ 185 (2014) 2250–2300, [1310.1921].
  * (56) T. Flacke, B. Fuks, M. Kunkel and L. Panizzi, “eVLQ S012.” http://feynrules.irmp.ucl.ac.be/wiki/NLOModels.
  * (57) C. Degrande, C. Duhr, B. Fuks, D. Grellscheid, O. Mattelaer and T. Reiter, _UFO - The Universal FeynRules Output_ , _Comput. Phys. Commun._ 183 (2012) 1201–1214, [1108.2040].
  * (58) J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, O. Mattelaer et al., _The automated computation of tree-level and next-to-leading order differential cross sections, and their matching to parton shower simulations_ , _JHEP_ 07 (2014) 079, [1405.0301].
  * (59) V. Hirschi and O. Mattelaer, _Automated event generation for loop-induced processes_ , _JHEP_ 10 (2015) 146, [1507.00020].
  * (60) NNPDF collaboration, R. D. Ball, V. Bertone, S. Carrazza, L. Del Debbio, S. Forte, A. Guffanti et al., _Parton distributions with QED corrections_ , _Nucl. Phys. B_ 877 (2013) 290–320, [1308.0598].
  * (61) B. Fuks, M. Klasen, D. R. Lamprea and M. Rothering, _Precision predictions for electroweak superpartner production at hadron colliders with Resummino_ , _Eur. Phys. J. C_ 73 (2013) 2480, [1304.0790].
  * (62) T. Sjöstrand, S. Ask, J. R. Christiansen, R. Corke, N. Desai, P. Ilten et al., _An introduction to PYTHIA 8.2_ , _Comput. Phys. Commun._ 191 (2015) 159–177, [1410.3012].
  * (63) DELPHES 3 collaboration, J. de Favereau, C. Delaere, P. Demin, A. Giammanco, V. Lemaître, A. Mertens et al., _DELPHES 3, A modular framework for fast simulation of a generic collider experiment_ , _JHEP_ 02 (2014) 057, [1307.6346].
  * (64) M. Cacciari, G. P. Salam and G. Soyez, _FastJet User Manual_ , _Eur. Phys. J._ C72 (2012) 1896, [1111.6097].
  * (65) ATLAS collaboration, _Technical Design Report for the ATLAS Inner Tracker Pixel Detector_ , Tech. Rep. ATLAS-TDR-030, 2017.
  * (66) ATLAS collaboration, _Expected performance of the ATLAS detector at the High-Luminosity LHC_ , Tech. Rep. ATL-PHYS-PUB-2019-005, 2019.
  * (67) D. Bertolini, P. Harris, M. Low and N. Tran, _Pileup Per Particle Identification_ , _JHEP_ 10 (2014) 059, [1407.6013].
  * (68) CMS collaboration, A. M. Sirunyan et al., _Pileup mitigation at CMS in 13 TeV data_ , _JINST_ 15 (2020) P09018, [2003.00503].
  * (69) A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan et al., _Pytorch: An imperative style, high-performance deep learning library_ , in _Advances in Neural Information Processing Systems_ (H. Wallach, H. Larochelle, A. Beygelzimer, F. d’ Alché-Buc, E. Fox and R. Garnett, eds.), vol. 32, Curran Associates, Inc., 2019, https://proceedings.neurips.cc/paper_files/paper/2019/file/bdbca288fee7f92f2bfa9f7012727740-Paper.pdf.
  * (70) D. P. Kingma and J. Ba, _Adam: A Method for Stochastic Optimization_ , 1412.6980.
  * (71) G. Cowan, K. Cranmer, E. Gross and O. Vitells, _Asymptotic formulae for likelihood-based tests of new physics_ , _Eur. Phys. J. C_ 71 (2011) 1554, [1007.1727].
  * (72) ATLAS collaboration, G. Aad et al., _Search for R-parity-violating supersymmetry in a final state containing leptons and many jets with the ATLAS experiment using $\sqrt{s}=13$ TeV proton–proton collision data_, _Eur. Phys. J. C_ 81 (2021) 1023, [2106.09609].
  * (73) X. Glorot, A. Bordes and Y. Bengio, _Deep sparse rectifier neural networks_ , in _Proceedings of the Fourteenth International Conference on Artificial Intelligence and Statistics_ (G. Gordon, D. Dunson and M. Dudík, eds.), vol. 15 of _Proceedings of Machine Learning Research_ , (Fort Lauderdale, FL, USA), pp. 315–323, PMLR, 11–13 Apr, 2011, https://proceedings.mlr.press/v15/glorot11a.html.
  * (74) S. Sabour, N. Frosst and G. E. Hinton, _Dynamic routing between capsules_ , in _Advances in Neural Information Processing Systems_ , vol. 30, Curran Associates, Inc., 2017, https://proceedings.neurips.cc/paper/2017/file/2cad8fa47bbef282badbb8de5374b894-Paper.pdf.
  * (75) M. Gori, G. Monfardini and F. Scarselli, _A new model for learning in graph domains_ , in _Proceedings. 2005 IEEE International Joint Conference on Neural Networks, 2005._ , vol. 2, pp. 729–734 vol. 2, 2005, DOI.
  * (76) F. Scarselli, M. Gori, A. C. Tsoi, M. Hagenbuchner and G. Monfardini, _The graph neural network model_ , _IEEE Transactions on Neural Networks_ 20 (2009) 61–80.
  * (77) J. Gilmer, S. S. Schoenholz, P. F. Riley, O. Vinyals and G. E. Dahl, _Neural message passing for quantum chemistry_ , _CoRR_ abs/1704.01212 (2017) , [1704.01212].