# Exploring Exotic Decays of the Higgs Boson to Multi-Photons at the LHC via
Multimodal Learning Approaches

A. Hammada, P. Kob, Chih-Ting Luc and Myeonghun Parkb,d

(  a Theory Center, IPNS, KEK, 1-1 Oho, Tsukuba, Ibaraki 305-0801, Japan.
b School of Physics, KIAS, Seoul 02455, Korea.
c Department of Physics and Institute of Theoretical Physics, Nanjing Normal
University, Nanjing, 210023, China.
d School of Natural Sciences, Seoultech, Seoul 01811, Korea. )

###### Abstract

The Standard Model (SM) Higgs boson, the most recently discovered elementary
particle, may still serve as a mediator between the SM sector and a new
physics sector related to dark matter (DM). The Large Hadron Collider (LHC)
has not yet fully constrained the physics associated with the Higgs boson,
leaving room for such possibilities. Among the various potential mass scales
of the dark sector, the sub-GeV mass range is particularly intriguing. This
parameter space presents significant challenges for DM direct detection
experiments that rely on nuclear recoils. Various innovative experimental
methods are currently under investigation to explore this sub-GeV dark sector.
The LHC, functioning as a Higgs factory, could explore this sector once the
challenge of identifying DM signals is resolved. Due to the significantly
lower mass of particles in the dark sector compared to the Higgs boson, these
particles are expected to be highly boosted following the Higgs boson’s decay.
However, detecting and identifying these highly boosted particles remains a
considerable challenge at the LHC, despite their eventual decay into SM
particles. We employ a well-motivated leptophobic $Z^{\prime}_{B}$ model as a
prototype to analyze the distinctive signatures from Higgs boson exotic decays
into multi-photons. These signatures consist of collimated photons that fail
to meet the photon isolation criteria, forming jet-like objects. Conventional
analyses relying solely on the purity of energy deposits in the
electromagnetic calorimeter would fail to detect these signatures, as they
would be overwhelmed by background events from Quantum Chromodynamics. To
effectively distinguish between such novel signal signatures and SM background
events, we leverage advanced machine learning techniques, specifically the
transformer encoder in a multimodal network structure. This neural network
successfully separates dark sector signals from SM backgrounds, outperforming
traditional event selection analyses.

###### Contents

  1. 1 Introduction
  2. 2 The simplified leptophobic $Z^{\prime}_{B}$ model
    1. 2.1 The model
    2. 2.2 Exotic Higgs boson decays
  3. 3 Constraining the model free parameters
  4. 4 Search strategy
    1. 4.1 Events generation and detector simulation
    2. 4.2 Preselection and event reconstruction
  5. 5 Transformer network
    1. 5.1 Self-attention
    2. 5.2 Networks structure
    3. 5.3 Data structure
  6. 6 Results
    1. 6.1 Projected limits on the model parameters
    2. 6.2 Decoding networks output
  7. 7 Conclusion
  8. A Analytical formulae of light $h_{1}$ partial decay widths
  9. B Analytical formulae of light $Z^{\prime}_{B}$ partial decay widths

##  1 Introduction

The discovery of the Higgs boson [1, 2] provides strong evidence supporting
the theory of electroweak symmetry breaking in the Standard Model of Particle
Physics [3, 4, 5]. However, unresolved puzzles remain beyond the Standard
Model, including the origin of tiny neutrino masses, dark matter (DM), and
matter-antimatter asymmetry in the universe. Dark sector models are a popular
choice to address these issues, and the Higgs boson can act as a portal
connecting the SM sector and the dark sector [6]111With the exception of the
Higgs portal, the dark photon/dark Z portal presents an attractive approach to
connect the visible sector with the dark sector [7, 8, 9, 10, 11].. In
particular, if new particles are lighter than the Higgs boson, the exotic
decays of the Higgs boson offer powerful ways to detect new physics beyond the
SM [12].

Another mystery in the SM is the origin of accidental global $U(1)$ symmetries
of baryon and lepton numbers which may suggest the presence of new physics
beyond the SM. A widely discussed proposal involves elevating these accidental
global $U(1)$ symmetries to gauged ones [13, 14, 15]. Furthermore, the
$U(1)_{B}$ gauge symmetry finds application in the $E_{6}$ grand unified
theory [16, 17, 18]. Consequently, the leptophobic $Z^{\prime}_{B}$ model
emerges as one of the well-motivated $U(1)^{\prime}$ models beyond the SM,
particularly those incorporating DM candidates. Examples of UV completions for
such models, with or without DM candidates, are detailed in Refs. [19, 20,
21]. The scenario featuring light (baryonic) DM, in particular, stands out as
one of the promising facets for DM models involving the leptophobic
$Z^{\prime}_{B}$ [22, 23]. In this study, our focus is on a simplified model
featuring $U(1)_{B}$ gauge symmetry, assuming that both the leptophobic
$Z^{\prime}_{B}$ and the new scalar boson are much lighter than the Higgs
boson. Specifically, the new scalar boson primarily decays into a pair of sub-
GeV $Z^{\prime}_{B}$ particles, with each $Z^{\prime}_{B}$ predominantly
decaying into a neutral pion and a photon, if kinematically allowed.

The decay of the Higgs boson into a pair of $Z^{\prime}_{B}$ or new scalar
boson results in both of these light particles being highly boosted.
Furthermore, the decay products from the highly boosted $Z^{\prime}_{B}$ or
new scalar boson tend to form a group of highly collimated photons. Notably,
depending on the detector resolution and photon isolation criteria, this group
of highly collimated photons can be identified as either a reconstructed
photon ($\gamma_{\text{rec}}$) or a photon-jet ($J_{\gamma}$) [24, 25]. If the
opening angle distribution ($\Delta R_{\gamma\gamma}$) between a pair of
photons within this group of collimated photons is smaller than the size of a
single cell (i.e., $\Delta R_{\gamma\gamma}<0.04$), the ATLAS detector cannot
distinguish whether the signal originates from a single photon or a group of
extremely collimated photons [26]. We refer to this signature as
$\gamma_{\text{rec}}$. On the other hand, if some $\Delta R_{\gamma\gamma}$
distributions within the group of collimated photons fall within the range of
$0.04<\Delta R_{\gamma\gamma}<0.4$, it is large enough to be distinguished
from a single photon but too collimated to satisfy the photon isolation
criteria. In such cases, we cannot identify them as multiple photons but
rather as a jet-like structure of photons, denoted as $J_{\gamma}$, in the
final state. The proposal for this kind of novel photon-jet signature was
presented in Refs. [27, 28]. Investigating the structure of photon-jets
involved the application of jet substructure techniques in Refs. [71, 30, 24,
31], and machine learning techniques were employed in Refs. [32, 33, 34].
Further studies on photon-jet signatures are available in Refs. [35, 36, 24,
37, 25, 38]. Signal events of the Higgs boson exotic decays to multi-photons
can classified into three categories, (1) two reconstructed photons [39, 40,
35, 25], (2) a reconstructed photon plus a photon-jet as the newly introduced
signature here, and (3) two photon-jets which have been studied with jet
substructures as above.

As the signal with photon-jet faces huge QCD-related backgrounds, one needs to
have additional handle more than an object tagging as previously studied. The
phase space of the signal events can be different from the one of background
events, either due to the different structure of event-topology or a different
weight from the amplitude. Thus, to maximize the performance in suppressing
backgrounds, we need to optimize the process of combining information from
kinematics and QCD [41, 42]. Accordingly, for object tagging, we utilize an
attention-based transformer encoder [43] to distinguish the signal from the SM
backgrounds. The proposed transformer model analyzes the events with the
structure of particle clouds, which is inherited by the points cloud
representation. Such a representation has the ability to share all the
advantages of particle based representations, especially the flexibility to
include arbitrary features for each particle. We endow this transformer model
with extra tokens in a multimodal structure to facilitate the model ability to
learn the full phase space of an event. Thus this network exploits the
different kinematics of the signal from backgrounds.

The paper is organized as follows. In Sec. 2, we introduce the simplified
leptophobic $Z^{\prime}_{B}$ model, while Sec. 3 summarizes the relevant
constraints for this model. Sec. 4.1 discusses the classification of signal
signatures and presents benchmark points. Sec. 4.2 showcases the signal-to-
background analysis for Higgs boson exotic decays into a photon-jet and a
reconstructed photon as well as two photon-jets at the LHC. All relevant
information for the transformer network used in this work is collected in Sec.
5. We present our results in Sec. 6. Finally, our conclusions are presented in
Sec. 7.

##  2 The simplified leptophobic $Z^{\prime}_{B}$ model

In this section, we initially introduce the simplified leptophobic
$Z^{\prime}_{B}$ model in Sec. 2.1. Subsequently, we delineate the exotic
partial decay widths of the SM-like Higgs boson $h_{2}$ in Sec. 2.2. The
partial decay widths of the new light scalar boson $h_{1}$ and
$Z^{\prime}_{B}$ can be found in the Appendix A and Appendix B, respectively.

###  2.1 The model

The model with $U(1)_{B}$ gauge symmetry is one of the well-motivated
$U(1)^{\prime}$ models beyond the SM [19, 20, 21]. In this simplified model,
we assume that all exotic new particles (exotic fermions for anomaly
cancellations and extra scalar bosons that provide masses to the exotic
fermions and/or make exotic charged/colored particles decay fast enough) are
at a relatively high energy scale, except for the light leptophobic
$Z^{\prime}_{B}$ and new scalar boson $h_{1}$. The associated $U(1)_{B}$
charges $Q^{i}_{B}$ for fields in our simplified model are listed in Table
1222Note that the Dirac (Majorana) fermion DM candidate $\chi_{D}$
($\chi_{M}$) with $Q^{i}_{B}(\chi_{D,M})=1$ can be involved in this simplified
leptophobic $Z^{\prime}_{B}$ model. Since the DM part is not the focus of this
study, we will address this issue in the future..

Field | quarks | leptons | $H$ | $\Phi$
---|---|---|---|---
$Q^{i}_{B}$ | $1/3$ | $0$ | $0$ | $-2$
Table 1: The associated $U(1)_{B}$ charges $Q^{i}_{B}$ for fields in our
simplified leptophobic $Z^{\prime}_{B}$ model.

Here, $H$ represents the SM Higgs doublet, and $\Phi$ is a new complex SM
singlet scalar responsible for generating the mass term for $Z^{\prime}_{B}$
via the Higgs mechanism.

The Lagrangian in this simplified model can be represented as

${\cal L}={\cal L}_{\text{SM}}+{\cal L}^{\prime}_{\text{scalar}}+{\cal
L}^{\prime}_{\text{gauge}},$ (1)

where

$\displaystyle{\cal
L}^{\prime}_{\text{scalar}}=D_{\mu}\Phi^{\dagger}D^{\mu}\Phi-
m_{\Phi}^{2}\Phi^{\dagger}\Phi-\lambda_{\Phi}(\Phi^{\dagger}\Phi)^{2}-\lambda_{H\Phi}H^{\dagger}H\Phi^{\dagger}\Phi,$
$\displaystyle{\cal
L}^{\prime}_{\text{gauge}}=-\frac{1}{4}X_{\mu\nu}X^{\mu\nu}-\frac{\epsilon}{2}B_{\mu\nu}X^{\mu\nu}+g_{B}\sum_{f}Q^{f}_{B}\overline{f}\hbox
to0.0pt{\,/\hss}Xf.$

The ${\cal L}_{\text{SM}}$ is the full SM Lagrangian and $f$ stands for SM
fermions. Covariant derivative can be represented as

$D_{\mu}=\partial_{\mu}-ig_{B}Q^{i}_{B}X_{\mu}.$ (2)

$\epsilon$ is the kinetic mixing parameter betweenthe two $U(1)$s, $g_{B}$ is
the $U(1)_{B}$ gauge coupling constant and $Q^{i}_{B}$ is shown in Table 1.

Expanding the scalar fields around the vacuum expectation values (VEVs) with
the unitary gauge as

$H(x)=\frac{1}{\sqrt{2}}\left(\begin{tabular}[]{c}0\\\
$v+h(x)$\end{tabular}\right)\;\;\;,\;\;\;\Phi(x)=\frac{1}{\sqrt{2}}\left(v_{B}+h_{B}(x)\right).$
(3)

After applying the condition for potential minimization, the VEVs $v$ and
$v_{B}$ can be written as

$v^{2}=\frac{4\lambda_{\Phi}\mu^{2}-2\lambda_{H\Phi}\mu^{2}_{\Phi}}{4\lambda_{\Phi}\lambda-\lambda^{2}_{H\Phi}},\quad
v^{2}_{B}=\frac{4\lambda\mu^{2}_{\Phi}-2\lambda_{H\Phi}\mu^{2}}{4\lambda_{\Phi}\lambda-\lambda^{2}_{H\Phi}},$
(4)

where $\lambda$ is the quartic self-coupling. Finally, these two neutral
scalar bosons $h(x)$ and $h_{B}(x)$ will mix to each other after the symmetry
breaking, and they can be diagonalized to the mass eigenstates $h_{1},h_{2}$
via

$\tan
2\theta_{h}=\frac{\lambda_{H\Phi}vv_{B}}{\lambda_{H\Phi}v^{2}_{B}-\lambda
v^{2}},$ (5)

and

$M^{2}_{h_{1},h_{2}}=(\lambda v^{2}+\lambda_{H\Phi}v^{2}_{B})\mp\sqrt{(\lambda
v^{2}-\lambda_{H\Phi}v^{2}_{B})^{2}+\lambda_{H\Phi}v^{2}v^{2}_{B}}.$ (6)

In our convention, we assign $h_{2}$ as the SM-like Higgs boson with
$M_{h_{2}}=125$ GeV, and $h_{1}$ as the light scalar boson with $M_{h_{1}}\ll
M_{h_{2}}$.

On the other hand, the three neutral gauge bosons
$B_{\mu},W^{3}_{\mu},X_{\mu}$ are diagonalized to the mass eigenstates
$A_{\mu},Z_{\mu},Z^{\prime}_{B\mu}$333 The details of this digonalization can
be found in Ref. [44, 45, 46] and we will not repeat them again. ,

$M^{2}_{\gamma}=0,\quad M^{2}_{Z,Z^{\prime}_{B}}=(\Delta_{1}+\Delta_{2})/2,$
(7)

where

$\displaystyle\Delta_{1}=g^{2}_{B}(Q^{\Phi}_{B})^{2}v^{2}_{B}\eta^{2}+\frac{1}{4}(g^{2}+g^{\prime
2}\eta^{2})v^{2},\hskip
11.38109pt\Delta_{2}=\sqrt{\Delta_{1}^{2}-g^{2}_{B}(Q^{\Phi}_{B})^{2}v^{2}v^{2}_{B}(g^{2}+g^{\prime
2})\eta^{2}}.$

We assume $M_{Z^{\prime}_{B}}\ll M_{Z}$ and
$\eta\equiv\epsilon/\sqrt{1-\epsilon^{2}}$. Considering small kinematic mixing
$\epsilon$, as in this study, $M_{Z}$ and $M_{Z^{\prime}_{B}}$ can be
approximated as

$M_{Z}\simeq\sqrt{g^{2}+g^{\prime 2}}v/2,\quad M_{Z^{\prime}_{B}}\simeq
g_{B}|Q^{\Phi}_{B}|v_{B}.$ (8)

Although SM charged leptons do not carry the $U(1)_{B}$ charge, the
$Z^{\prime}_{B}$ cannot be completely decoupled from them. The kinematic
mixing between $Z^{\prime}_{B}$ and the photon, resulting from one-loop
radiative corrections involving SM quarks and heavy charged fermions, can
generate a parameter $\epsilon\sim\frac{eg_{B}}{(4\pi)}$ [47, 48]. We will
consider this effect in our calculations below.

###  2.2 Exotic Higgs boson decays

The light scalar boson $h_{1}$ and the leptophobic $Z^{\prime}_{B}$ can lead
to interesting exotic decays of the SM-like Higgs boson, $h_{2}$, via the
mixing between $H$ and $\Phi$. Depending on the mass spectrum of $h_{1}$ and
$Z^{\prime}_{B}$, there are at least four possibilities for the exotic decays
of $h_{2}$:

$h_{2}\rightarrow Z^{\prime}_{B}Z^{\prime}_{B},\quad h_{2}\rightarrow
h_{1}h_{1},\quad h_{2}\rightarrow h_{1}h^{\ast}_{1}\rightarrow
h_{1}Z^{\prime}_{B}Z^{\prime}_{B},\quad h_{2}\rightarrow h_{1}h_{1}h_{1}.$ (9)

We are aware that the three-body decays in the third and fourth cases suffer
from phase space suppression. Although these two decay channels can be
numerically calculated, as shown in Ref. [46]. As their contributions are much
smaller compared to the first two decay channels, we will only consider the
partial decay widths for the first two cases below.

The decay width of $h_{2}$ can be represented as

$\Gamma_{h_{2}}=\cos^{2}\theta_{h}\widehat{\Gamma}_{h_{2}}+\Gamma^{\prime}_{h_{2}},$
(10)

where $\widehat{\Gamma}_{h_{2}}$ is the SM Higgs boson width, $4.03$ MeV [49],
and

$\Gamma^{\prime}_{h_{2}}\approx\sin^{2}\theta_{h}\widehat{\Gamma}(h_{2}\rightarrow
Z^{\prime}_{B}Z^{\prime}_{B})+\Gamma(h_{2}\rightarrow h_{1}h_{1}),$ (11)

where

$\displaystyle\widehat{\Gamma}(h_{2}\rightarrow
Z^{\prime}_{B}Z^{\prime}_{B})=\frac{g^{2}_{B}M^{2}_{Z^{\prime}_{B}}}{8\pi
M_{h_{2}}}\sqrt{1-\frac{4M^{2}_{Z^{\prime}_{B}}}{M^{2}_{h_{2}}}}(3-\frac{M^{2}_{h_{2}}}{M^{2}_{Z^{\prime}_{B}}}+\frac{M^{4}_{h_{2}}}{4M^{4}_{Z^{\prime}_{B}}}),$
(12) $\displaystyle\Gamma(h_{2}\rightarrow
h_{1}h_{1})=\frac{\mu^{2}_{h_{2}h_{1}h_{1}}}{32\pi
M_{h_{2}}}\sqrt{1-\frac{4M^{2}_{h_{1}}}{M^{2}_{h_{2}}}},$ (13)

with

$\displaystyle\mu_{h_{2}h_{1}h_{1}}=$ $\displaystyle
2\cos^{3}\theta_{h}\lambda_{H\Phi}v+2\sqrt{2}\cos^{2}\theta_{h}\sin\theta_{h}(\lambda_{H\Phi}-6\lambda_{\Phi})v_{B}$
$\displaystyle+\cos\theta_{h}\sin^{2}\theta_{h}(6\lambda-4\lambda_{H\Phi})v-\sqrt{2}\sin^{3}\theta_{h}\lambda_{H\Phi}v_{B}.$

Assuming that $2M_{Z^{\prime}_{B}}<M_{h_{1}}\ll M_{h_{2}}$ and fix
$\alpha_{B}=10^{-5}$, $M_{Z^{\prime}_{B}}=0.4$ GeV and
$\sin\theta_{h}=1.6\times 10^{-4}$, the branching ratios of $h_{2}\rightarrow
Z^{\prime}_{B}Z^{\prime}_{B}$ and $h_{2}\rightarrow h_{1}h_{1}$ can be as
large as $4.35\times 10^{-4}$ and $8.70\times 10^{-4}$, respectively.

Figure 1:  The branching ratios of $h_{1}$ when $h_{1}\rightarrow
Z^{\prime}_{B}Z^{\prime}_{B}$ is kinematically allowed with fixed
$\alpha_{B}=10^{-7}$, $\sin\theta_{h}=10^{-3}$, and
$M_{Z^{\prime}_{B}}=0.4M_{h_{1}}$ (left panel) and the proper decay length
($c\tau$) [mm] of $h_{1}$ with fixed $\sin\theta_{h}=10^{-3}$,
$M_{Z^{\prime}_{B}}=0.4M_{h_{1}}$ and varying $\alpha_{B}=10^{-5}$, $10^{-6}$
and $10^{-7}$ (right panel).

The lighter scalar, $h_{1}$, can couple either directly to the
$Z^{\prime}_{B}$ pair or to SM particles via the mixing between $H$ and
$\Phi$. Analytical formulas for the partial decay widths of $h_{1}$ can be
found in Appendix A.

In this study, we focus on the case where$h_{1}\rightarrow
Z^{\prime}_{B}Z^{\prime}_{B}$ is kinematically allowed. We fix
$\alpha_{B}=10^{-7}$, $\sin\theta_{h}=10^{-3}$, and
$M_{Z^{\prime}_{B}}=0.4M_{h_{1}}$ to illustrate the branching ratios of
$h_{1}$ in the left panel of Fig. 1. The right panel of Fig. 1 displays the
proper decay lengths ($c\tau$) [mm] of $h_{1}$ with varying
$\alpha_{B}=10^{-5}$, $10^{-6}$, and $10^{-7}$ while keeping
$\sin\theta_{h}=10^{-3}$ and $M_{Z^{\prime}_{B}}=0.4M_{h_{1}}$ fixed. We
observe that once $h_{1}\rightarrow Z^{\prime}_{B}Z^{\prime}_{B}$ is
kinematically allowed, it becomes the dominant decay mode, and $h_{1}$
promptly decays to the final state for the entire parameter space.

Figure 2:  The decay branching ratios of sub-GeV $Z^{\prime}_{B}$ (left panel)
and its proper decay length ($c\tau$) [mm] with $\alpha_{B}=10^{-4}$,
$10^{-5}$, $10^{-6}$ and $10^{-7}$ (right panel). Here we assume
$\epsilon=0.1\times\frac{eg_{B}}{(4\pi)}$.

The important feature of the light $Z^{\prime}_{B}$ is its decay branching
ratios, which differ significantly from the usual kinematic mixing dark photon
scenario. For $M_{Z^{\prime}_{B}}\lesssim 1$ GeV, the dominant decay modes of
$Z^{\prime}_{B}$ are $Z^{\prime}_{B}\rightarrow\pi^{0}\gamma$ and
$Z^{\prime}_{B}\rightarrow\pi^{+}\pi^{-}\pi^{0}$. However,
$Z^{\prime}_{B}\rightarrow e^{+}e^{-}$ becomes the dominant decay mode when
pion decays are kinematically forbidden. The partial decay widths for
$Z^{\prime}_{B}$ to light mesons are calculated using the vector meson
dominance (VMD) model [50], and the detailed partial decay widths of
$Z^{\prime}_{B}$ are shown in Appendix B, following closely Ref. [51]. In the
left panel of Fig. 2, we present the decay branching ratios for sub-GeV
$Z^{\prime}_{B}$. The right panel of Fig. 2 displays the proper decay length
($c\tau$) [mm] of $Z^{\prime}_{B}$ with different values of
$\alpha_{B}=10^{-4}$, $10^{-5}$, $10^{-6}$, and $10^{-7}$. In order to avoid
constraints from KLOE and LHCb, which depend on the couplings of
$Z^{\prime}_{B}$ to lepton pairs, we assume
$\epsilon=0.1\times\frac{eg_{B}}{(4\pi)}$ in Fig. 2. Importantly, the light
$Z^{\prime}_{B}$ can decay promptly when $\pi^{0}\gamma$ is kinematically
allowed444For the case of a light $Z^{\prime}_{B}$ as a long-lived particle
(LLP) with a small enough $\alpha_{B}$ at the LHC, the charged track
information from the displaced vertex in the tracker system can record the
decay mode of $Z^{\prime}_{B}\rightarrow\pi^{+}\pi^{-}\pi^{0}$. However,
effectively finding a trigger for the decay mode of
$Z^{\prime}_{B}\rightarrow\pi^{0}\gamma$ when the $Z^{\prime}_{B}$ becomes the
LLP poses a challenge..

##  3 Constraining the model free parameters

In this section, we consider various constraints in the simplified leptophobic
$Z^{\prime}_{B}$ model. First, we address the constraints on $h_{2}$ and
$h_{1}$. The global fits of Higgs boson data from ATLAS in Ref. [52] suggest
that the allowed branching ratio of $h_{2}$ exotic decays is below $19\%$ at
the $95\%$ confidence level. We apply this constraint to our scenario as
follows:

$\frac{\Gamma^{\prime}_{h_{2}}}{\Gamma_{h_{2}}}\lesssim 0.19,$ (14)

where $\Gamma_{h_{2}}$ and $\Gamma^{\prime}_{h_{2}}$ are given in Eqs. (10)
and (11), respectively. Additionally, $h_{1}$ with $M_{h_{1}}\lesssim{\cal
O}(1)$ GeV is subject to severe constraints from fixed-target experiments and
$B$ meson decay processes if $h_{1}\rightarrow Z^{\prime}_{B}Z^{\prime}_{B}$
is kinematically forbidden. Regions of $M_{h_{1}}<280$ MeV with
$\sin\theta_{h}\geq 10^{-5}$ have already been ruled out. However, there are
still some allowed parameter space with sizable $\sin\theta_{h}$ for
$M_{h_{1}}>280$ MeV. Since we focus on $M_{h_{1}}\geqslant 1$ GeV and
$h_{1}\rightarrow Z^{\prime}_{B}Z^{\prime}_{B}$ becomes the dominant decay
channel in this work, the above constraints are much weaker in our scenario.
Therefore, we conservatively require $\sin\theta_{h}\lesssim 10^{-3}$ to the
parameter space of interest. For more details on constraints from a light
scalar mixed with the SM Higgs boson, refer to Refs. [53, 54, 55].

Figure 3:  Various constraints for the light leptophobic $Z^{\prime}_{B}$ on
the $(M_{Z^{\prime}_{B}},\alpha_{B})$ plane. Each constraint is explained in
the main text.

The sub-GeV leptophobic $Z^{\prime}_{B}$ is subject to severe constraints from
various meson decay processes. We will begin by summarizing potential
constraints from the $Z^{\prime}_{B}\rightarrow\pi^{0}\gamma$ decay channel:

  * •

$\eta\rightarrow\pi^{0}\gamma\gamma$: An updated measurement of
$\eta\rightarrow\pi^{0}\gamma\gamma$ was published by Crystal-Ball@MAMI [56],
yielding a more precise result for the partial decay width:
$\Gamma(\eta\rightarrow\pi^{0}\gamma\gamma)=(0.33\pm 0.03_{\text{tot}})$ eV. A
previous reanalysis of data from CrystalBall@AGS [57] also provided similar
results with $\Gamma(\eta\rightarrow\pi^{0}\gamma\gamma)=(0.285\pm
0.068_{\text{tot}})$ eV. Here we applied the value reported by the PDG,
BR$(\eta\rightarrow\pi^{0}\gamma\gamma)=(2.56\pm 0.22)\times 10^{-4}$ [58,
59]. This measurement can be applied to constrain the branching ratio,
BR$(\eta\rightarrow Z^{\prime}_{B}\gamma\rightarrow\pi^{0}\gamma\gamma)$, for
the $Z^{\prime}_{B}$, as indicated by the purple line in Fig. 3.

  * •

$\eta^{\prime}\rightarrow\pi^{0}\gamma\gamma$: A new measurement of
$\eta^{\prime}\rightarrow\pi^{0}\gamma\gamma$ was reported by BESIII [60],
yielding a more precise result for the nonresonant branching fraction:
BR$(\eta^{\prime}\rightarrow\pi^{0}\gamma\gamma)=(6.16\pm 0.64(\text{stat})\pm
0.67(\text{sys}))\times 10^{-4}$. This measurement can be used to constrain
the branching ratio BR$(\eta^{\prime}\rightarrow
Z^{\prime}_{B}\gamma\rightarrow\pi^{0}\gamma\gamma)$ for the $Z^{\prime}_{B}$,
as illustrated by the magenta line in Fig. 3.

  * •

$\phi\rightarrow\eta\pi^{0}\gamma$: The measurement of
$\phi\rightarrow\eta\pi^{0}\gamma$ by KLOE [61] provides a constraint on the
branching fraction, BR$(\phi\rightarrow\eta\pi^{0}\gamma)=(7.06\pm 0.22)\times
10^{-5}$. This constraint can be applied to the branching ratio
BR$(\phi\rightarrow\eta Z^{\prime}_{B}\rightarrow\eta\pi^{0}\gamma)$ for the
$Z^{\prime}_{B}$, as illustrated by the blue line in Fig. 3.

  * •

$\omega\rightarrow\eta\pi^{0}\gamma$: The upper limit on the branching
fraction, Br$(\omega\rightarrow\eta\pi^{0}\gamma)<3.3\times 10^{-5}$ ($90\%$
C.L.), was obtained from CMD-2@VEPP-2M [62]. This constraint can be applied to
the branching ratio BR$(\omega\rightarrow\eta
Z^{\prime}_{B}\rightarrow\eta\pi^{0}\gamma)$ for the $Z^{\prime}_{B}$, as
shown by the black line in Fig. 3.

Next, we consider the constraint from the
$Z^{\prime}_{B}\rightarrow\pi^{+}\pi^{-}\pi^{0}$ decay mode with
$M_{Z^{\prime}_{B}}>620$ MeV within the
$\eta^{\prime}\rightarrow\pi^{+}\pi^{-}\pi^{0}\gamma$ channel. CLEO conducted
a search for
$\eta^{\prime}\rightarrow\omega\gamma\rightarrow\pi^{+}\pi^{-}\pi^{0}\gamma$
from charmonium decays with BR$(\eta^{\prime}\rightarrow\omega\gamma)=(2.34\pm
0.30\pm 0.04)\%$ [63]. However, BESIII updated the measurement to
BR$(\eta^{\prime}\rightarrow\omega\gamma)=(2.489\pm 0.018\pm 0.074)\%$ [64].
We will apply the BESIII measurement to constrain BR$(\eta^{\prime}\rightarrow
Z^{\prime}_{B}\gamma\rightarrow\pi^{+}\pi^{-}\pi^{0}\gamma)$ for the
$Z^{\prime}_{B}$, as indicated by the red line in Fig. 3.

Notice some constraints from KLOE [65] and WASA [66] are dependent on
couplings of $Z^{\prime}_{B}$ to leptons which are proportional to $\epsilon$.
Here we assume $\epsilon$ is much smaller
($\epsilon<0.1\times\frac{eg_{B}}{(4\pi)^{2}}$) in this study, so those
constraints become much weaker [51, 67] as shown for the pink and green bulks
in Fig. 3. Finally, the hadronic decay of $\Upsilon(1S)$ state mediated by an
s-channel $Z^{\prime}_{B}$ provide $\alpha_{B}<0.014$ which is independent of
$M_{Z^{\prime}_{B}}$ if $\sqrt{s}\gg M_{Z^{\prime}_{B}}$ [48] as shown for the
yellow bulk in Fig. 3. Therefore, $\alpha_{B}\lesssim 10^{-5}$ is still
allowed for $M_{Z^{\prime}_{B}}\gtrsim 200$ MeV.

##  4 Search strategy

Among the various processes offered by the rich parameter space of our
benchmark leptophobic $Z^{\prime}_{B}$ model, we explore two intriguing
topologies arising from Higgs boson exotic decays into multiple photons in the
Higgs production through gluon fusion:

$\displaystyle\text{TP1}:gg\rightarrow h_{2}\rightarrow
Z^{\prime}_{B}Z^{\prime}_{B}\rightarrow(\pi^{0}\gamma)(\pi^{0}\gamma)\rightarrow
6\gamma,$ $\displaystyle\text{TP2}:gg\rightarrow h_{2}\rightarrow
h_{1}h_{1}\rightarrow(Z^{\prime}_{B}Z^{\prime}_{B})(Z^{\prime}_{B}Z^{\prime}_{B})\rightarrow(\pi^{0}\gamma)(\pi^{0}\gamma)(\pi^{0}\gamma)(\pi^{0}\gamma)\rightarrow
12\gamma.$ (15)

Since we consider $M_{h_{1}}\sim{\cal O}(1)$ GeV and a sub-GeV
$Z^{\prime}_{B}$, both $h_{1}$ and $Z^{\prime}_{B}$ are highly boosted.
Photons emerging from highly boosted $h_{1}$ and $Z^{\prime}_{B}$ are
collimated enough and can be identified as a single photon-like energy
cluster, $\gamma_{\text{rec}}$, or a jet-like structure, $J_{\gamma}$, in the
electromagnetic calorimeter depending on the opening angles among these
photons.

Figure 4:  Signal signatures at the detector: (a) Two reconstructed photons,
$\gamma_{\text{rec}}\gamma_{\text{rec}}$ (left panel); (b) A reconstructed
photon and a photon-jet, $\gamma_{\text{rec}}J_{\gamma}$ (middle panel); (c)
Two photon-jets, $J_{\gamma}J_{\gamma}$ (right panel).

Based on the EM calorimeter structure of ATLAS detector [26], the angular
separation between two photons in the final state $\Delta
R_{\gamma\gamma}\lesssim 0.04$ is approximately the same size as a standard
single photon energy cluster. In this case, existing triggers cannot
distinguish a calorimeter energy deposit resulting from highly collimated
photons from that of a single isolated photon. Therefore, we use the
Cambridge/Aachen (C/A) jet clustering algorithm [68, 69] with $R=0.04$ to
group these collimated photons. For those photons inside $R=0.04$, they form a
single reconstructed photon object, $\gamma_{\text{rec}}$. However, for
$0.04<R<0.4$, we use the anti-$k_{T}$ jet clustering algorithm [70] with the
jet radius of $R=0.4$ to group collimated photons as is done in ordinary QCD
jet analysis at the LHC. These photons are unable to pass the photon isolation
criteria and thus are defined to form a photon jet, $J_{\gamma}$ [71].
Accordingly, the signal signatures at the detector can be classified into

  1. (a)

Two reconstructed photons, $\gamma_{\text{rec}}\gamma_{\text{rec}}$,

  2. (b)

A reconstructed photon and a photon-jet, $\gamma_{\text{rec}}J_{\gamma}$,

  3. (c)

Two photon-jets, $J_{\gamma}J_{\gamma}$,

The three signal typologies are shown in Fig. 4 in which $\gamma_{\text{rec}}$
or $J_{\gamma}$ are defined according to the group of photons and the value of
the radius parameter. Notice that each $J_{\gamma}$ can include up to 3
collinear photons for the TP1 process and up to 6 collinear photons for the
TP2 process, respectively.

###  4.1 Events generation and detector simulation

To study these signal signatures, we establish a benchmark point (BP) for TP1
and four BPs for TP2, all with fixed values of $\alpha_{B}=10^{-5}$,
$M_{Z^{\prime}_{B}}=0.4$ GeV, and $\sin\theta_{h}=1.6\times 10^{-4}$. It’s
worth noting that for these BPs, $BR(h_{2}\rightarrow h_{1}h_{1})$ is
approximately $8.70\times 10^{-4}$, and $BR(h_{2}\rightarrow
Z^{\prime}_{B}Z^{\prime}_{B})$ is roughly $4.35\times 10^{-4}$. UFO model file
is generated from FeynRules [72] for the simplified leptophobic
$Z^{\prime}_{B}$. Madgraph5 aMC@NLO [73, 74] is used to compute the total
cross section and event generation. Generated events are showered by Pythia8
[75] to model the parton showering and hadronization. Delphes3 [76] with
modifed ATLAS template is used for the fast detecor simulation. FastJets [77]
is used to cluster collimated photons in the final state and apply the basic
photon isolation criteria in Delphes3, allowing to distinguish between
$\gamma_{\text{rec}}$ and $J_{\gamma}$. The results are summarized in Table 2
and the left panel of Fig. 5, with the average number of reconstructed photons
denoted as $\langle N_{\gamma}\rangle$. For TP1, the dominant classification
is $\gamma_{\text{rec}}\gamma_{\text{rec}}$ while TP2–BP1 and TP2–BP2 exhibit
sizable probabilities for $\gamma_{\text{rec}}J_{\gamma}$. In contrast,
TP2–BP3 and TP2–BP4 are primarily classified as $J_{\gamma}J_{\gamma}$.

BPs |  $M_{h_{1}}$ (GeV) | $\gamma_{\text{rec}}\gamma_{\text{rec}}$ | $\gamma_{\text{rec}}J_{\gamma}$ | $J_{\gamma}J_{\gamma}$ | $\langle N_{\gamma}\rangle$
---|---|---|---|---|---
TP1 | $-$ | $75.56\%$ | $21.70\%$ | $2.75\%$ | $1.65$
TP2–BP1 | $1.0$ | $5.55\%$ | $31.18\%$ | $63.27\%$ | $0.40$
TP2–BP2 | $1.2$ | $1.65\%$ | $18.71\%$ | $79.64\%$ | $0.21$
TP2–BP3 | $1.4$ | $4.97\times 10^{-3}$ | $9.49\%$ | $90.02\%$ | $0.10$
TP2–BP4 | $1.6$ | $2.58\times 10^{-3}$ | $4.78\%$ | $94.96\%$ | $0.05$
Table 2: Branching fractions for $\gamma_{\text{rec}}\gamma_{\text{rec}}$,
$\gamma_{\text{rec}}J_{\gamma}$, $J_{\gamma}J_{\gamma}$ and the number of
reconstructed photons, $\langle N_{\gamma}\rangle$ for the five BPs with fixed
$\alpha_{B}=10^{-5}$, $M_{Z^{\prime}_{B}}=0.4$ GeV and
$\sin\theta_{h}=1.6\times 10^{-4}$ .

Branching fractions for the number of reconstructed photons corresponding to
these BPs in the right panel of Fig. 5, for further confirmation of this
observation. , for these BPs, is also provided in Table 2. In the following a
detailed analysis, particularly regarding the $\gamma_{\text{rec}}J_{\gamma}$
and $J_{\gamma}J_{\gamma}$ signatures, are considered.

Figure 5:  The branching fractions for
$\gamma_{\text{rec}}\gamma_{\text{rec}}$, $\gamma_{\text{rec}}J_{\gamma}$, and
$J_{\gamma}J_{\gamma}$ (left panel), as well as the number of reconstructed
photons, $\langle N_{\gamma}\rangle$, (right panel) for the five BPs shown in
Table 2.

###  4.2 Preselection and event reconstruction

In this subsection, we begin by examining the impact of the
$\gamma_{\text{rec}}\gamma_{\text{rec}}$ signature on the measurement of the
SM-like Higgs boson to diphoton. Subsequently, we provide a detailed signal-
to-background analysis and the corresponding significance for TP1 and TP2, as
presented in Table 2, with a focus on the $\gamma_{\text{rec}}J_{\gamma}$ and
$J_{\gamma}J_{\gamma}$ signatures.

In Table 2, when the final state with a pair of reconstructed photons cannot
be distinguished from two isolated single photons in the ATLAS EM calorimeter,
the signal can mimic the SM-like Higgs boson decay to diphoton [40, 39].
According to Ref. [78], the observed (SM expected) production rate for
inclusive SM-like Higgs boson to diphoton is $(\sigma\times
B_{\gamma\gamma})_{\text{obs}}=127\pm 10$ fb ($(\sigma\times
B_{\gamma\gamma})_{\text{exp}}=116\pm 5$ fb). We estimate the contributions
from TP1 and TP2 to the SM-like Higgs boson to diphoton production rate as
follows: First, we adopt the total decay width of $h_{2}$ in Eq. (10) to
compute $BR(h_{2}\rightarrow\gamma\gamma)$. Then, we rescale $(\sigma\times
B_{\gamma\gamma})_{\text{exp}}$ in our simplified leptophobic $Z^{\prime}_{B}$
model accordingly. The total cross section, $(\sigma\times
B_{\gamma\gamma})_{\text{BSM}}$, can be estimated as

$(\sigma\times
B_{\gamma\gamma})_{\text{BSM}}\equiv\sigma\cos^{2}\theta_{h}\times\left(BR(h_{2}\rightarrow
Z^{\prime}Z^{\prime})\epsilon^{\text{iso}}_{\text{TP1}}+BR(h_{2}\rightarrow
h_{1}h_{1})\epsilon^{\text{iso}}_{\text{TP2}}\right)\times\left(\frac{\epsilon^{\text{ID}}_{\text{BSM}}}{\epsilon^{\text{ID}}_{\text{SM}}}\right)^{2}\,,$
(16)

where $\epsilon^{\text{iso}}$ represents the efficiency for collimated photons
to satisfy the photon isolation criteria by ATLAS [26]; more details about
this criteria discussed later. The $\epsilon^{\text{ID}}$ follows Eq. (1) in
Ref. [40], where "ID" refers to photon identification based on various shower
shape variables in different calorimeter layers. According to Fig. 7 in Ref.
[40], we estimate
$\epsilon^{\text{ID}}_{\text{BSM}}/\epsilon^{\text{ID}}_{\text{SM}}\sim 0.7$
for this study. Although, $\epsilon^{\text{ID}}_{\text{BSM}}$ for TP2 should
be smaller than the estimated value, the exact value can still computed using
the GEANT4 detector simulation [79], which is beyond the scope of this work.

We first define the combinations as follows: "TP1+TP2-BP1" as Scenario I,
"TP1+TP2-BP2" as Scenario II, "TP1+TP2-BP3" as Scenario III, and "TP1+TP2-BP4"
as Scenario IV. Next, we present $\epsilon^{\text{iso}}_{\text{TP1}}$,
$\epsilon^{\text{iso}}_{\text{TP2}}$, and $(\sigma\times
B_{\gamma\gamma})_{\text{BSM}}$ based on five BPs in Table 2 for these four
scenarios:

  * •

Scenario I: $\epsilon^{\text{iso}}_{\text{TP1}}=45.34\%$,
$\epsilon^{\text{iso}}_{\text{TP2-BP1}}=3.13\%$, $(\sigma\times
B_{\gamma\gamma})_{\text{BSM}}=5.62$ fb;

  * •

Scenario II: $\epsilon^{\text{iso}}_{\text{TP1}}=45.34\%$,
$\epsilon^{\text{iso}}_{\text{TP2-BP2}}=7.39\times 10^{-3}$, $(\sigma\times
B_{\gamma\gamma})_{\text{BSM}}=5.10$ fb;

  * •

Scenario III: $\epsilon^{\text{iso}}_{\text{TP1}}=45.34\%$,
$\epsilon^{\text{iso}}_{\text{TP2-BP3}}=1.09\times 10^{-3}$, $(\sigma\times
B_{\gamma\gamma})_{\text{BSM}}=4.97$ fb;

  * •

Scenario IV: $\epsilon^{\text{iso}}_{\text{TP1}}=45.34\%$,
$\epsilon^{\text{iso}}_{\text{TP2-BP4}}=6.60\times 10^{-5}$, $(\sigma\times
B_{\gamma\gamma})_{\text{BSM}}=4.97$ fb.

We proceed with a detailed analysis of the $\gamma_{\text{rec}}J_{\gamma}$ and
$J_{\gamma}J_{\gamma}$ signatures at the LHC. The primary sources of SM
background for these two signal signatures originate from

$\displaystyle pp\rightarrow\gamma\gamma,\quad pp\rightarrow\gamma j,\quad
pp\rightarrow jj.$ (17)

Anti-$k_{T}$ algorithm is used to cluster calorimeter cells into jets with
$R=0.4$. We adhere closely to the photon isolation criteria outlined in Ref.
[26]. Reconstructed photons must be isolated from other calorimeter energy
deposits and nearby charged tracks not associated with the photon. The
calorimeter isolation variable, denoted as $E^{\text{iso}}_{T}$, is defined as
the sum of energy deposits in the calorimeter within a cone of size $\Delta
R=0.4$ centered around the barycenter of the photon cluster (excluding the
energy contributed by the photon cluster itself), minus $0.022\times E_{T}$.
We require the calorimeter isolation variable to satisfy
$E^{\text{iso}}_{T}<2.45$ GeV. The track isolation variable, denoted as
$p^{\text{iso}}_{T}$, is the scalar sum of the transverse momenta of charged
tracks not associated with the photon within a cone of size $\Delta R=0.2$
centered around the barycenter of the photon cluster. It must satisfy
$p^{\text{iso}}_{T}<0.05\times E_{T}$. A photon candidate is defined as a
reconstructed photon, $\gamma_{\text{rec}}$, if it satisfies the following
criteria: $P^{\gamma}_{T}>20$ GeV, $|\eta_{\gamma}|<1.37$ or
$1.52<|\eta_{\gamma}|<2.37$, and meets the aforementioned photon isolation
criteria. After particle reconstruction, we define a photon-jet, $J_{\gamma}$,
which must fulfill the following three conditions:

$P^{J_{\gamma}}_{T}>20\text{ GeV}\quad\text{and}\quad|\eta_{J_{\gamma}}|<2.5,$
(18)
$\theta_{J_{\gamma}}=\frac{E_{J_{\gamma},\text{HCAL}}}{E_{J_{\gamma}}}\leqslant
0.1\,,$ (19) $\nu_{J_{\gamma}}=0\quad\text{with}\quad p_{T}>2\text{ GeV},$
(20)

where $\theta_{J_{\gamma}}$ represents the hadronic energy fraction, and
$\nu_{J_{\gamma}}$ denotes the number of charged particles within
$J_{\gamma}$.

We further apply the following baseline event selections to identify signal
events with the $\gamma_{\text{rec}}J_{\gamma}$ signature:

$\displaystyle P_{T}(\gamma_{\text{rec1}}),P_{T}(J_{\gamma 1})>0.35\
m_{\gamma_{\text{rec1}}J_{\gamma 1}},$ $\displaystyle 120\text{
GeV}<m_{\gamma_{\text{rec1}}J_{\gamma 1}}<130\text{ GeV},$
$\displaystyle|\eta_{\gamma_{\text{rec1}}}-\eta_{J_{\gamma 1}}|<2.5,$ (21)

where the subscript label corresponds to the leading momentum photon.
Reconstructed photon and the photon-jet are required to be highly boosted,
with their invariant mass distributions mainly falling within the Higgs boson
mass window. Similarly, we apply the following baseline event selections to
identify signal events with the $J_{\gamma}J_{\gamma}$ signature:

$\displaystyle P_{T}(J_{\gamma 1})>0.4m_{J_{\gamma 1}J_{\gamma
2}}\quad\text{and}\quad P_{T}(J_{\gamma 2})>0.3m_{J_{\gamma 1}J_{\gamma 2}},$
$\displaystyle 110\text{ GeV}<m_{J_{\gamma 1}J_{\gamma 2}}<140\text{ GeV},$
$\displaystyle|\eta_{J_{\gamma 1}}-\eta_{J_{\gamma 2}}|<2.5.$ (22)

Furthermore, both photon-jets need to be highly boosted, enabling the
reconstruction of the Higgs boson resonance peak from their invariant mass.
The pileup effect is not considered in this study because it can be mitigated
using traditional strategies, such as area subtraction [80] and charged hadron
subtraction [81]. Additionally, modern machine learning algorithms like PUMML
[82] and others [83, 84, 85] can effectively eliminate the pileup effect.

##  5 Transformer network

Transformer models, introduced by Ref. [43], have revolutionized the field of
deep learning and introduced in HEP in Refs. [86, 87, 42, 88, 89] . Unlike
traditional sequence-based models, transformers utilize multi-heads self-
attention mechanism to process input data in parallel, making them highly
efficient for capturing long-range dependencies. This inherent capability
makes transformers promising candidates for analyzing complex and
interconnected systems, such as particle clouds. Adapting transformer models
to analyze particle clouds involves encoding the spatial and temporal dynamics
inherent in event collisions. Representing each particle as an input token,
transformers can effectively capture interactions and correlations between
particles. Moreover, self-attention mechanism allows the model to weigh the
significance of each particle’s contribution to the overall cloud behavior. In
this manner, transformer model can capture the important features shared by
each individual particle to the whole cloud allowing for better classification
performance between different signal and background events.

In this paper, we employ two distinct networks utilizing attention-based
transformer layers for TP1 and TP2 analysis. Both networks are designed to
comprehensively analyze the event information encoded by the substructure of
the photon-jet and reconstructed high-level kinematics. Since the information
extracted from each benchmark point exhibits different structures, which
cannot be effectively processed by a single network, we adapt multimodal
networks comprising transformer and MLP layers for efficient modality
learning.

###  5.1 Self-attention

Self-attention mechanism leverages the relationships between particles within
a cloud to compute attention scores, which determine the importance of each
particle’s features. By aggregating information from neighboring particles,
self-attention enables the model to capture complex patterns and dependencies
within the particle cloud, facilitating various tasks such as event
classification and jet tagging. In this work we use the transformer encoder to
extract the substructure information of the photon-jet. Accordingly, we
construct the input to the transformer layers as unordered sequence of
particles inside the jet with a representation of

$X=x_{1}^{i},x_{2}^{i}\ \cdots\ x_{n}^{i}\,,$ (23)

where $i$ is the number of the corresponding features of each particle inside
the jet and $n$ is the maximum number of included particles. To compute the
attention score, input representation is divided into three sets of vectors,
Query (Q), Key (K) and Value (V). These vectors are obtained by linear
projection of a fully connected layer to assign a learnable weight matrices as

$Q=XW^{Q},\hskip 11.38109ptK=XW^{K},\hskip 11.38109ptV=XW^{V}\,,$ (24)

where $W$ are the learnable weight matrices assigned to each vector via the
linear projection.

Attention score can be obtained as scaled dot product of the three vectors as

$A=\frac{\text{exp}(Q.K^{T}/\sqrt{d})}{\sum\text{exp}(Q.K^{T}/\sqrt{d})}\,,$
(25)

where the sum is over the $K$ index and $d$ is the length of the input
data555Normalizing the attention score is important to avoid the divergence of
the dot product.. Attention score is a matrix with dimensions equal to the
number of the input particle tokens and provide a flexible mechanism for
modeling dependencies and relationships between particles in a sequence. By
assigning different weights to each particle based on its relevance to others,
attention mechanisms can adaptively focus on the most salient information for
the classification task. Moreover, it provides a degree of interpretability by
indicating which particles in the input sequence are most influential for
making predictions. The output of a single attention head is obtained as the
dot product between the attention score and the Value vector,
$\mathcal{Z}=A\cdot V$. The output of the attention head ensures the
permutation invariance over the particles in the dataset by summing all the
particle token via the matrix multiplication process. In order to speed up the
computation and increase the rich of the network to capture more relevant
information multi-heads can be used in parallel while their outputs can be
stacked as

$\widetilde{X}=\text{concatenate}\left(\mathcal{Z}_{1},\mathcal{Z}_{2}\
\cdots\ \mathcal{Z}_{k}\right)\times W\,,$ (26)

where $k$ is the number of the parallel attention heads and $W$ is a
normalizing matrix to keep the output dimensions as same as the inputs. The
output is then added to the input, via an addition wise layer, to weigh each
of the particle tokens according to its computed attention. Exploiting the
fact that transformer layer preserve the dimensions of the input dataset
allows to stack many layers for better network expressivity.

###  5.2 Networks structure

To fully leverage the capabilities of the two discussed channels
$\left(\gamma_{\rm rec}J_{\gamma},\,J_{\gamma}J_{\gamma}\right)$, we implement
two distinct configurations of the transformer encoder with self-attention
heads, as presented in Fig. 6. For the analysis of TP1, where we encounter a
single photon-jet and one reconstructed photon, we tailor a transformer model
with multi-head self-attention to analyze the constituents of the photon-jet,
augmented with additional distributions representing the reconstructed
isolated photon kinematics.

The network is fed with two input datasets: the photon-jet contents, presented
as a sequence of particles in a $(50,12)$ dimensional space, analyzed by
transformer layers, and a supplementary dataset comprising five distributions,
processed by a Fully Connected (FC) layer. The transformer layer is iterated
three times, each with seven attention heads. The output from the attention
heads is then integrated with the original input data via a skip connection
layer. The output of the final transformer layer is flattened and combined
with the second dataset, then fed into an averaging layer. Subsequently, the
output of the averaging layer is forwarded to an MLP with two FC layers of
$128$ and $64$ neurons, respectively, utilizing Rectified Linear (ReLU)
activation functions. Dropout layers, with dropout rate of $20\%$, follow each
FC layer of the MLP, and the final output is directed to an FC layer with two
neurons and softmax activation.

Figure 6:  Left: Network architecture used for TP1 analysis. Middle: Network
architecture used for TP2 analysis. Right: structure of the attention-based
transformer layer used in the two networks. N indicates the number of repeated
transformer layer.

The second network used for TP2 analysis, consists of twins transformer
encoders utilized to analyze the leading and subleading photon-jet
constituents. The network takes two separate data sets for the leading and
subleading jet constituents as input with dimensions of $(50,12)$, processed
through two distinct transformer layers. Each transformer layer is repeated
three times and has the same structure as the transformer layer in TP1
network. In order to efficiently merge the output of the two transformers, we
adopt an Additive attention layer [90]. Additive attention dynamically
computes attention weights based on the compatibility between the output of
the two transformers, allowing the model to focus on relevant parts of both of
them. Compatibility score between the output of the first transformer,
$H_{i}$, and the output of the second transformer, $E_{j}$, is computed as

$A_{i,j}=\upsilon^{T}\cdot\text{tanh}\left(W_{1}\cdot H_{i}+W_{2}\cdot
E_{j}\right)\,,$ (27)

where $\upsilon$ is a learnable parameter vector, $i.e,$ attention vector. The
output of the additive attention layer is obtained as

$\mathcal{O}_{i}=\sum_{j}\text{softmax}\left(A_{i,j}\right)\cdot E_{j}\,.$
(28)

The output vector capture the weighted combination of elements from the second
data set for each element of the first data set, allowing the model to
leverage information from both data sets to perform classification task. This
output then passes through a MLP comprising two fully connected layers with
dimensions $128$ and $64$, employing the GELU activation function. Following
each fully connected layer, a dropout layer with a dropout rate of $20\%$ is
applied. The output is then passed to the output layer for classification.

Output of the two networks is a probability with dimension
($\mathcal{P}_{sig},\mathcal{P}_{bkg}$), with $\mathcal{P}$ ranges from
$[0,1]$. For event classification if the $\mathcal{P}_{sig}>0.5$, the event is
classified as signal like event and if $\mathcal{P}_{sig}<0.5$ the event is
classified as most likely background event.

###  5.3 Data structure

Figure 7: Features of the leading photon-jet used as input to the transformer
layers. Each distribution is obtained by averaging over the jet constituents
in each event.

Configuring diverse data inputs for the network is achieved through particle
clouds, emphasizing the permutation symmetry of inputs to yield a promising
representation of jets. Initially, we pre-process the data sets for the
leading photon-jet contents by imposing a constraint of 50 constituents in
each. The particles are arranged based on their transverse momentum in
descending order. In instances where events comprise more than 50
constituents, only the top 50 are retained. For events with fewer
constituents, the remaining positions are padded with zeros, ensuring
conformity with the stipulated counts. An attention mask is used to notify the
network with the padded events to not be further considered. In this case the
results are independent on the prior hypothesized number of constituents. We
follow the variables input to the transformer layer from Ref. [86], with the
following:

  1. 1.

$\Delta\eta_{c}=\eta_{c}-\eta_{jet}$, where $\eta_{c}(\eta_{jet})$ is the
pseudorapidity of the jet constituents (reconstructed jet).

  2. 2.

$\Delta\phi_{c}=\phi_{c}-\phi_{jet}$, where $\phi_{c}(\phi_{jet})$ is the
azimuth angle of the jet constituents (reconstructed jet).

  3. 3.

$p_{T_{c}}$, transverse momentum of the jet constituents in GeV.

  4. 4.

$E_{c}$, energy of the jet constituents in GeV.

  5. 5.

$\log\left(\frac{p_{T_{c}}}{p_{T_{Jet}}}\right)$, transverse momentum of the
jet constituents relative to the jet transverse momentum.

  6. 6.

$\log\left(\frac{E_{c}}{E_{{Jet}}}\right)$, energy of the jet constituents
relative to the jet energy.

  7. 7.

$\Delta R_{c}=\sqrt{\Delta\eta^{2}_{c}+\Delta\phi^{2}_{c}}$, the angular
distance of the jet constituents from the jet axis.

  8. 8.

Particle ID for the jet constituents.

The particle ID features are fed to the transformer model as a one-hot-
encoding resulting extra 5 columns ending up with input features of dimensions
$(50,12)$. To optimize the network discriminative accuracy, it is imperative
to pre-process the jet contents, ensuring the manifestation of a multi-prong
structure specific to signal events. For this purpose, the following
transformations are applied before inputting the data into the network.
Features of the leading photon-jets which are used as input to the transformer
layers are shown in Fig.7.

Figure 8: Normalized events distributions for the extra tokens for the signal
and backgrounds used for TP1 process. With $E\left(\frac{Hadronic}{EM}\right)$
is the hadronic over electromagnetic energy of the leading jet and
$N\left(\frac{ChargedH}{NeutralH}\right)$ is the numebr of the charged over
the neutral particles inside the leading jet.

For TP1 process we consider full event information by investigating the
kinematics of the reconstructed photon. Consequently, the second data set
comprises five distributions representing the kinematics of the reconstructed
photon and the photon-jet, as shown in Fig.8. For this purpose we consider the
transverse momentum, pseudorapidity and azimuthal angle of the reconstructed
photon with the relative energy deposit of the photon-jet in the different sub
detectors.

After preparing the data sets, we proceed to train the networks to recognize
the intricate connections between the input data and their corresponding
labels, which denote signal and background events. Signal events are labeled
by Y = 1, while background events are labeled by Y = 0. To ensure the
network’s independence from the arrangement of signal and background events,
we merge them into a single data set and randomize the order along with their
assigned labels. During each epoch of network training (defined as a complete
pass through the entire data sets), the network adjusts the weights of its
neurons for each event through backward error propagation. The objective is to
minimize the disparity between its predictions and the actual labels by
converging to a global minimum of a defined loss function. This iterative
process continues until the desired accuracy is achieved.

##  6 Results

In this section, we present the signal significance obtained by setting a cut
on the networks output to maximize the signal to the background yield. Having
the considered analysis for the two BPs, TP1 and TP2–BP2, at the $14$ TeV LHC
with integrated luminosity of $3000$ fb-1 (HL-LHC), the analysis results are
presented in Fig. 9. The discriminating power of each of the networks will be
a measure of how well the signal and background characterised through their
different features. Measuring the classification efficiency of each network
Area Under the ROC (AUC) metric is reported with $0.986$ and $0.998$ for TP1
and TP2–BP2, respectively.

Figure 9: Normalized output score of the transformer networks for TP1 (left)
and TP2–BP2 (right). Blue and orange distributions present signal and
background-like events, respectively. Red dashed line indicates the signal
significance by setting a cut on the output score. Signal significance is
obtained at the HL-LHC.

In Fig. 9, the left panel illustrates the normalized output scores of the
networks for TP1 analysis, while the right panel depicts TP2–BP2 analysis. The
network output, represented as
$(\mathcal{P}{\text{sig}},1-\mathcal{P}{\text{bkg}})$, ranges from 0 to 1.
Events with a discriminant value close to 1 are classified as signal-like
events, depicted by the blue distribution, whereas those nearing 0 are
categorized as background-like events, represented by the orange distribution.
The intersection area between these two distributions measures the
misclassification rate of each network. Signal significance is assessed by
reweighting the network output for signal and background-like events based on
the number of events remaining after preselection cuts. For TP1, there are
$6375$ signal events and $5.6\times 10^{9}$ combined background events, while
for TP2–BP2, there are $8700$ signal events and $2.5\times 10^{10}$ background
events post the preselection cuts. Enhancing the signal to background yield is
done by maximizing the cut on the network output, which represented by the red
dashed line. The signal significance is approximated using the formula
$\frac{S}{\sqrt{S+B}}$, where $S$ and $B$ represent the number of signal and
background events, respectively. It shows a maximum value of $16\sigma$ for
TP1 and $11\sigma$ for TP2–BP2.

###  6.1 Projected limits on the model parameters

Based on the baseline event selections and the transformer encoder analysis
described above, we can extrapolate our results to establish future bounds on
model parameters at the HL-LHC. Using $\frac{S}{\sqrt{S+B}}=2$ for the $95\%$
Confidence Level (CL), we can determine the relevant cross-sections for signal
and background processes from Fig. 9. For TP1, the corresponding signal cross-
section is $0.27$ fb, while the total SM background cross-section is $53.5$
fb. For TP2-BP2, the corresponding signal cross-section is $0.49$ fb, and the
total SM background cross-section is $179$ fb. Although our analysis has
significantly reduced the size of background cross-sections, they still exceed
the signal cross-sections by two orders of magnitude.

Figure 10: Left: Projected limits (blue dashed line) at $95\%$ CL on the plane
$(M_{Z^{\prime}_{B}},\alpha_{B})$ for TP1 at the HL-LHC. We fix
$\sin\theta_{H}=10^{-3}$ and label the red star as our BP. For comparison, the
existing constraints from Fig. 3 are also shown. Right: The projection of
limits (blue dashed line) at $95\%$ CL on the plane
$(M_{h_{1}},|\mu_{h_{2}h_{1}h_{1}}|)$ for TP2 at the HL-LHC. Here,
$\sin\theta_{H}=10^{-3}$, $\alpha_{B}=10^{-5}$, and $M_{Z^{\prime}_{B}}=400$
MeV are fixed. Again, the red star is marked as our BP. The existing
constraint for Higgs boson exotic decays from Eq. 14 is shown for comparison.
The gray dashed lines correspond to including a $20\%$ systematic uncertainty
of the SM background events in Eq. 29.

For TP1, $h_{1}$ is not involved. To simplify the analysis, we can exclude the
contribution from $h_{1}$ and focus solely on the process $h_{2}\rightarrow
Z^{\prime}_{B}Z^{\prime}_{B}$. According to Eqs. (10), (11), and (12), there
are three free model parameters: $\sin\theta_{H}$, $\alpha_{B}$, and
$M_{Z^{\prime}_{B}}$. Assuming $\sin\theta_{H}=10^{-3}$, future bounds on the
plane $(M_{Z^{\prime}_{B}},\alpha_{B})$ can be predicted, as indicated by the
blue dashed line in the left panel of Fig. 10. The red star represents our BP,
which yielded $\alpha_{B}<4.41\times 10^{-7}$ for $M_{Z^{\prime}_{B}}=400$ MeV
at $95\%$ CL. It’s important to note that only the BP (red star) is from our
actual analysis; we assume that other $Z^{\prime}_{B}$ masses (blue dashed
line) have the same efficiency as the BP. Since the light $Z^{\prime}_{B}$ is
highly boosted from the Higgs boson decay, we expect the efficiency for other
$Z^{\prime}_{B}$ masses to be only mildly different from that of the BP.
However, the dashed line rises rapidly as $M_{Z^{\prime}_{B}}$ increases,
mainly due to the shape of the branching ratio of
$Z^{\prime}_{B}\rightarrow\pi^{0}\gamma$ in Fig. 3. Consequently, future
bounds on $\alpha_{B}$ may weaken as the size of $\sin\theta_{H}$ decreases.

For TP2, both $h_{1}$ and $Z^{\prime}_{B}$ are involved, leading to the
simultaneous appearance of Higgs boson exotic decays, $h_{2}\rightarrow
h_{1}h_{1}$ and $h_{2}\rightarrow Z^{\prime}_{B}Z^{\prime}_{B}$. In this
scenario, there are five free model parameters: $\sin\theta_{H}$,
$\alpha_{B}$, $M_{Z^{\prime}_{B}}$, $M_{h_{1}}$, and
$|\mu_{h_{2}h_{1}h_{1}}|$. By fixing $\sin\theta_{H}=10^{-3}$,
$\alpha_{B}=10^{-5}$, and $M_{Z^{\prime}_{B}}=400$ MeV, future bounds on the
plane $(M_{h_{1}},|\mu_{h_{2}h_{1}h_{1}}|)$ can be predicted, as indicated by
the blue dashed line in the right panel of Fig. 10. Here, the red star
represents our BP, which yielded $|\mu_{h_{2}h_{1}h_{1}}|<0.12$ GeV for
$M_{h_{1}}=1.2$ GeV at $95\%$ CL. Similarly, except for the BP (red star),
other $M_{h_{1}}$ values (blue dashed line) are assumed to have the same
efficiency as the BP. The constraint for Higgs boson exotic decays from Eq. 14
is also added in Fig. 10 for comparison. Moreover, future bounds on
$|\mu_{h_{2}h_{1}h_{1}}|$ can be weakened when one or both of $\sin\theta_{H}$
and $\alpha_{B}$ decrease.

Regarding the systematic uncertainties arising from QCD di-jet measurements,
which may affect the reported classification performance, we consider a $20\%$
systematic uncertainty to estimate this effect. The formula for signal
significance is modified accordingly:

$Z=\sqrt{2\cdot\left[(S+B)\cdot
ln\left(\frac{(S+B)(B+\sigma^{2}_{B})}{B^{2}+(S+B)\sigma^{2}_{B}}\right)-\frac{B^{2}}{\sigma^{2}_{B}}\cdot
ln\left(1+\frac{\sigma^{2}_{B}S}{B(B+\sigma^{2}_{B})}\right)\right]}\,,$ (29)

where $\sigma_{B}$ represents the systematic uncertainty of the SM background
events $B$. In our analysis results, the background events $B$ are much larger
than the signal events $S$, so the projected limits accounting for a $20\%$
systematic uncertainty on the model parameters only change slightly at $95\%$
CL, as shown by the gray lines in Fig. 10.

###  6.2 Decoding networks output

Interpretable ML methods are crucial because they provide transparency and
insights into how models make decisions, which is essential for trust,
accountability, and debugging. Analyzing the network’s latent space, in
particular, is crucial because it reveals how the model internally represents
and processes data, shedding light on the underlying structure and patterns it
has learned. This understanding helps in diagnosing and correcting errors,
enhancing model robustness, and ensuring that the model’s decision-making
process aligns with the physics of the considered process.

After determining the signal significance of the two benchmark points by
optimizing the cut on the networks’ output, we proceed to delve into the
network decision making process by conducting a numerical analysis of each
network’s flow. In Fig. 11, we present the latent space distribution of both
the transformer and MLP, utilized for additional tokens analysis, weighted by
the output values of the element-wise average layer.

Figure 11:  Latent space distribution of the jet constituents and extra tokens
for $10^{5}$ test events (left) and background events (right). Color bar
indicates the output of the element-wise average layer.

The latent space representation of the MLP indicates a degree of similarity
between signal and background events, while the latent distribution of the
transformer layers showcases a broader variation among background events
compared to signal events. This disparity stems from the inherent structure of
the photon-jet in TP1, characterized by a prong structure where the particles
with the highest momenta are concentrated around the jet axis. In contrast,
the QCD jet lacks such a prong structure, resulting in varied momenta of the
particles within the jet distributed around the QCD jet axis.

Figure 12: Attention weights of the particle tokens for the leading and sub-
leading jet for the signal events (left) and background events (right). The
attention weights are normalized for the $10^{5}$ test images. Color bar
indicates the value of the attention scores as in Eq. (27).

For TP2, we investigate the network decision making by visualizing the
attention score, $A_{i,j}$ in Eq. (27). These scores have dimensions
corresponding to the number of particle tokens in the leading and subleading
jets, i.e., $(50,50)$, offering insight into the relationships and
dependencies between different particles in the respective datasets.
Essentially, attention scores highlight which particle tokens from the two
datasets, the leading and subleading jet constituents, play pivotal roles in
the network’s decision-making process.

Figure 12 visualizes the attention scores, $A_{i,j}$, generated by the
additive attention layer for $10^{5}$ signal events (left panel) and
background events (right panel). Intriguingly, the network predominantly
focuses on the first five highest momentum particles inside the leading and
subleading jet to identify signal events, whereas for background events, the
network distributes high attention across a wider range of particle contents.
We interpret this behavior as the network learning the distinct substructures
of photon and QCD jets.

##  7 Conclusion

Although the SM has been so successful in a wide range of energy scales, there
are a number of reasons to consider its extensions in various directions. One
possible extension would be $U(1)_{B}$, which may help to understand better
the longevity of proton or the origin of global $U(1)_{B}$ symmetry of the SM.
In this paper we considered a simplified model for $U(1)_{B}$ extensions of
the SM, assuming all the exotic particles are heavy except for the $U(1)_{B}$
gauge boson $Z_{B}^{{}^{\prime}}$ and the $U(1)_{B}$-breaking Higgs boson
$h_{B}\simeq h_{1}$. Here $h_{B}\simeq h_{1}$ can be pair-produced in the SM
Higgs decays, where $h_{\rm SM}(\simeq h_{2})$ decays into $h_{1}h_{1}$ if
kinematically allowed. Additionally, if this channel is open, $h_{1}$ will
mainly decay into a pair of $Z_{B}^{{}^{\prime}}$. For a certain mass region
of $Z_{B}^{{}^{\prime}}$ ($m_{\pi^{0}}\lesssim m_{Z_{B}^{{}^{\prime}}}\lesssim
620$ MeV), the main decay channel would be
$Z_{B}^{{}^{\prime}}\rightarrow\pi^{0}\gamma\rightarrow(\gamma\gamma)\gamma$,
and not $\pi^{+}\pi^{-}$ as shown in Fig. 2. This is because
$Z_{B}^{{}^{\prime}}$ is an iso-singlet vector meson, just like
$\omega^{0}(783)$ in ordinary QCD 666$\rho^{0}(770)$ is in the similar mass
region as $Z_{B}^{{}^{\prime}}$, and its main decay channel is
$\rho^{0}(770)\rightarrow\pi^{+}\pi^{-}$, since it is an isovector ($I=1$)
vector ($J=1$) meson. Note that $\rho^{0}\rightarrow\pi^{0}\pi^{0}$ is
strictly forbidden by the total angular momentum conservation and spin-
statistics theorem.. When $Z_{B}^{{}^{\prime}}$ boson is boosted, its decay
products will make a reconstructed photon or a photon-jet, unique signatures
of our simplified $U(1)_{B}$ model with light $Z_{B}^{{}^{\prime}}$ and
$h_{1}$ bosons as shown in Fig. 4. One has to pay extra care for its study at
high energy colliders, especially at the HL-LHC. Compared with previous
photon-jet studies at the LHC, such signatures emerging from $Z^{\prime}_{B}$
decays have not been considered before, as far as we know. This motivates us
to focus on these types of signatures.

For this purpose, we adopted advanced machine learning techniques, in
particular the transformer encoder to identify the novel photon-jet signature
and distinguish it from the SM background events. Unraveling the distinctive
characteristics of these typologies at the LHC presents a formidable challenge
due to significant contamination from the QCD di-jet process within the signal
events. To address this, and to enhance the signal-to-background ratio, we
employ two ML methodologies centered on self-attention transformer encoders.
Our investigation focuses on two specific processes: $h_{2}\to
Z^{\prime}_{B}Z^{\prime}_{B}$ (TP1) and $h_{2}\to h_{1}h_{1}$ (TP2), yielding
$6$ and $12$ photons, respectively. Within these final states, we examine
processes involving $\gamma_{reco}+J_{\gamma}$ and $J_{\gamma}+J_{\gamma}$.
Recognizing that each final state encapsulates distinct event features, we
employ two slightly different ML networks. The first network, analyze the
$\gamma_{reco}+J_{\gamma}$ final state, comprises a transformer layer for
analyzing photon-jet substructure, complemented by a MLP encoding the
kinematics of reconstructed photons and photon-jets. The second network,
focusing on the $J_{\gamma}+J_{\gamma}$ final state, incorporates two
transformer layers to analyze the substructure of both photon-jets. These two
network architectures are displayed in Fig. 6.

The structured architecture of these networks enables robust discrimination
between signal and background events, achieving classification accuracies
exceeding $98\%$, as shown in Fig. 9. Using our baseline event selections and
the transformer encoder analysis, we can forecast the future bounds of model
parameters in this simplified leptophobic $Z^{\prime}_{B}$ model at the HL-
LHC. For TP1, fixing the mixing angle between $h_{B}$ and $h_{\rm SM}$ at
$\sin\theta_{H}=10^{-3}$, we found that the baryonic force fine structure
constraint, $\alpha_{B}$, can be as low as $5\times 10^{-7}$ for
$M_{Z^{\prime}_{B}}\sim 400$ MeV at $95\%$ CL. For TP2, fixing both
$\sin\theta_{H}=10^{-3}$ and $\alpha_{B}=10^{-5}$, the three-point interaction
$h_{2}$-$h_{1}$-$h_{1}$ can be constrained to $|\mu_{h_{2}h_{1}h_{1}}|<0.12$
GeV for $M_{h_{1}}\sim 1$ GeV at $95\%$ CL. Ultimately, our results suggest
that probing photon-jet signatures remains feasible at the HL-LHC, even when
accounting for $20\%$ of systematic uncertainties. Therefore, searching for
Higgs boson exotic decays in energy frontier experiments can link the physics
at the electroweak scale to the dark sector at the low energy scale and serve
as a complementary probe to search for light $Z^{\prime}_{B}$ and $h_{1}$ in
intensity frontier experiments.

## Acknowledgment

This work is supported in part by the KIAS Individual Grants under Grant No.
PG021403 at Korea Institute for Advanced Study (P.K.), National Research
Foundation of Korea NRF-2021R1A2C4002551 (M. P. ), the National Natural
Science Foundation of China (NNSFC) under grant No. 12335005 (C.-T. Lu). AH is
funded by grant number 22H05113, “Foundation of Machine Learning Physics”,
Grant in Aid for Transformative Research Areas and 22K03626, Grant-in-Aid for
Scientific Research (C).

##  Appendix A Analytical formulae of light $h_{1}$ partial decay widths

In this appendix, we show analytical formulae of light $h_{1}$ partial decay
widths used in our work [46, 91]. For $M_{h1}\sim{\cal O}(1)$ GeV, the $h_{1}$
total width, $\Gamma_{h_{1}}$, can be written as

$\Gamma_{h_{1}}=\cos^{2}\theta_{h}\widehat{\Gamma}(h_{1}\rightarrow
Z^{\prime}_{B}Z^{\prime}_{B})+\sin^{2}\theta_{h}\left(\widehat{\Gamma}(h_{1}\rightarrow
f\overline{f})+\widehat{\Gamma}(h_{1}\rightarrow\gamma\gamma)\right)$ (30)

where $f\overline{f}$ can be charged lepton pairs or light quarks to form
meson pairs. These $h_{1}$ partial decay widths are

$\displaystyle\widehat{\Gamma}(h_{1}\rightarrow
Z^{\prime}_{B}Z^{\prime}_{B})=\frac{g^{2}_{B}M^{2}_{Z^{\prime}_{B}}}{8\pi
M_{h_{1}}}\sqrt{1-\frac{4M^{2}_{Z^{\prime}_{B}}}{M^{2}_{h_{1}}}}(3-\frac{M^{2}_{h_{1}}}{M^{2}_{Z^{\prime}_{B}}}+\frac{M^{4}_{h_{1}}}{4M^{4}_{Z^{\prime}_{B}}}),$
(31)
$\displaystyle\widehat{\Gamma}(h_{1}\rightarrow\ell^{+}\ell^{-})=\frac{m^{2}_{\ell}M_{h_{1}}}{8\pi
v^{2}_{B}}(1-\frac{4m^{2}_{\ell}}{M^{2}_{h_{1}}})^{\frac{3}{2}},$ (32)
$\displaystyle\widehat{\Gamma}(h_{1}\rightarrow\pi^{+}\pi^{-})=\frac{M^{3}_{h_{1}}}{324\pi
v^{2}_{B}}\sqrt{1-\frac{4m^{2}_{\pi}}{M^{2}_{h_{1}}}}(1+\frac{11m^{2}_{\pi}}{2M^{2}_{h_{1}}})^{2},$
(33)
$\displaystyle\widehat{\Gamma}(h_{1}\rightarrow\pi^{0}\pi^{0})=\frac{1}{2}\times\widehat{\Gamma}(h_{1}\rightarrow\pi^{+}\pi^{-}),$
(34) $\displaystyle\widehat{\Gamma}(h_{1}\rightarrow
K^{+}K^{-})=\widehat{\Gamma}(h_{1}\rightarrow
K^{0}\overline{K^{0}})=\frac{M^{3}_{h_{1}}}{324\pi
v^{2}_{B}}\sqrt{1-\frac{4m^{2}_{K}}{M^{2}_{h_{1}}}}(1+\frac{11m^{2}_{K}}{2M^{2}_{h_{1}}})^{2},$
(35)
$\displaystyle\widehat{\Gamma}(h_{1}\rightarrow\eta\eta)=\frac{M^{3}_{h_{1}}}{648\pi
v^{2}_{B}}\sqrt{1-\frac{4m^{2}_{\eta}}{M^{2}_{h_{1}}}}(1+\frac{11m^{2}_{\eta}}{2M^{2}_{h_{1}}})^{2},$
(36)
$\displaystyle\widehat{\Gamma}(h_{1}\rightarrow\gamma\gamma)=\frac{\alpha^{2}_{\text{EM}}M^{3}_{h_{1}}}{256\pi^{3}v^{2}_{B}}\lvert\sum_{i}N^{i}_{c}(Q^{i})^{2}F_{\frac{1}{2}}(\frac{4m^{2}_{i}}{M^{2}_{h_{1}}})+F_{1}(\frac{4m^{2}_{W}}{M^{2}_{h_{1}}})\rvert^{2},$
(37)

where $\alpha_{\text{EM}}$ is the fine structure constant, $N_{c}$ is the
color factor and $Q$ is the electric charge. The functions
$F_{1}(x),F_{\frac{1}{2}}(x)$ are

$\displaystyle F_{1}(x)=2+3x+3x(2-x)f(x),$ (38) $\displaystyle
F_{\frac{1}{2}}(x)=-2x[1+(1-x)f(x)],$ (39)

which

$\displaystyle
f(x)=\left\\{\begin{array}[]{lr}\left(\sin^{-1}\sqrt{1/x}\right)^{2},&\quad
x\geq 1,\\\
-\frac{1}{4}\left[\ln\left(\frac{{1+\sqrt{1-x}}}{{1-\sqrt{1-x}}}\right)-i\pi\right]^{2},&\quad
x<1.\end{array}\right.$ (42)

Notice there is no reliable description of non-perturbative QCD effects can be
presented in the mass region $M_{h_{1}}\simeq 1.5-2.5$ GeV [92]. Therefore, we
only focus on $M_{h_{1}}\lesssim 1.5$ GeV in this appendix.

##  Appendix B Analytical formulae of light $Z^{\prime}_{B}$ partial decay
widths

In this appendix, we collect analytical formulae of $Z^{\prime}_{B}$ partial
decay widths used in our work which closely follows Ref. [51]. The light
$Z^{\prime}_{B}$ total decay width, $\Gamma_{Z^{\prime}_{B}}$, can be written
as

$\displaystyle\Gamma_{Z^{\prime}_{B}}=$
$\displaystyle\Gamma(Z^{\prime}_{B}\rightarrow\pi^{0}\gamma)+\Gamma(Z^{\prime}_{B}\rightarrow\pi^{+}\pi^{-}\pi^{0})+\Gamma(Z^{\prime}_{B}\rightarrow\pi^{+}\pi^{-})+\Gamma(Z^{\prime}_{B}\rightarrow\eta\gamma)$
$\displaystyle+\Gamma(Z^{\prime}_{B}\rightarrow l^{+}l^{-})$ (43)

where

  * •

$\Gamma(Z^{\prime}_{B}\rightarrow\pi^{0}\gamma)$:

$\Gamma(Z^{\prime}_{B}\rightarrow\pi^{0}\gamma)=\frac{\alpha_{B}\alpha_{\text{EM}}M^{3}_{Z^{\prime}_{B}}}{96\pi^{3}f_{\pi}^{2}}(1-\frac{m^{2}_{\pi}}{M^{2}_{Z^{\prime}_{B}}})^{3}\lvert
F_{\omega}(M^{2}_{Z^{\prime}_{B}})\rvert^{2}$ (44)

where $f_{\pi}\simeq 93$ MeV is the pion decay constant and
$\alpha_{B}\equiv\frac{g^{2}_{B}}{4\pi}$. The form factor $F_{\omega}(s)$ is
given by $F_{\omega}(s)\approx
1-s/m^{2}_{\omega}-i\Gamma_{\omega}/m_{\omega}$.

  * •

$\Gamma(Z^{\prime}_{B}\rightarrow\pi^{+}\pi^{-}\pi^{0})$:

$\Gamma(Z^{\prime}_{B}\rightarrow\pi^{+}\pi^{-}\pi^{0})=\frac{g^{4}_{\rho\pi\pi}\alpha_{B}M_{Z^{\prime}_{B}}}{192\pi^{6}f_{\pi}^{2}}{\cal
I}(M^{2}_{Z^{\prime}_{B}})\lvert F_{\omega}(M^{2}_{Z^{\prime}_{B}})\rvert^{2}$
(45)

where $g_{\rho\pi\pi}$ is the $\rho\pi\pi$ coupling and fixed to be
$g^{2}_{\rho\pi\pi}/(4\pi)\simeq 3.0$ from the observed
$\rho\rightarrow\pi\pi$ decay rate. The details of the integral over phase
space ${\cal I}(M^{2}_{Z^{\prime}_{B}})$ can be found in the Appendix of Ref.
[51].

  * •

$\Gamma(Z^{\prime}_{B}\rightarrow\pi^{+}\pi^{-})$:

$\Gamma(Z^{\prime}_{B}\rightarrow\pi^{+}\pi^{-})=\frac{\alpha_{\text{EM}}\epsilon^{2}M_{Z^{\prime}_{B}}}{12}(1-\frac{4m^{2}_{\pi}}{M^{2}_{Z^{\prime}_{B}}})^{\frac{3}{2}}\lvert
F_{\pi}(M^{2}_{Z^{\prime}_{B}})\rvert^{2}$ (46)

where

$F_{\pi}(s)=F_{\rho}(s)[1+\frac{1+\delta}{3}\frac{\overline{\Pi}(s)}{s-m^{2}_{\omega}+im_{\omega}\Gamma_{\omega}}].$
(47)

If $\delta=0$, then $\overline{\Pi}(m^{2}_{\omega})=-3500\pm 300$ MeV2. Due to
the direct mixing between $Z^{\prime}_{B}$ and $\omega$ meson,
$\delta=2g_{B}/(\epsilon e)$ is considered. For lack of a better understanding
of $\overline{\Pi}(s)$, we assume a constant value
$\overline{\Pi}(s)\approx\overline{\Pi}(m^{2}_{\omega})$ here, although this
assumption is not well justified. On the other hand, we take a simple Breit-
Wigner form $F_{\rho}(s)\approx 1-s/m^{2}_{\rho}-i\Gamma_{\rho}/m_{\rho}$ as
an approximation.

  * •

$Z^{\prime}_{B}\rightarrow\eta\gamma$:

$\Gamma(Z^{\prime}_{B}\rightarrow\eta\gamma)=\frac{\alpha_{B}\alpha_{\text{EM}}M^{3}_{Z^{\prime}_{B}}}{96\pi^{3}f_{\eta}^{2}}(1-\frac{m^{2}_{\eta}}{M^{2}_{Z^{\prime}_{B}}})^{3}\lvert
F_{\phi}(M^{2}_{Z^{\prime}_{B}})\rvert^{2}$ (48)

where the form factor $F_{\phi}(s)$ is given by $F_{\phi}(s)\approx
1-s/m^{2}_{\phi}-i\Gamma_{\phi}/m_{\phi}$.

  * •

$Z^{\prime}_{B}\rightarrow\ell^{+}\ell^{-}$:

$\Gamma(Z^{\prime}_{B}\rightarrow\ell^{+}\ell^{-})=\frac{\alpha_{\text{EM}}\epsilon^{2}M_{Z^{\prime}_{B}}}{3}(1+\frac{2m^{2}_{\ell}}{M^{2}_{Z^{\prime}_{B}}})\sqrt{1-\frac{4m^{2}_{\ell}}{M^{2}_{Z^{\prime}_{B}}}}.$
(49)

## References

  * [1] G. Aad et al. [ATLAS], Phys. Lett. B 716, 1-29 (2012) doi:10.1016/j.physletb.2012.08.020 [arXiv:1207.7214 [hep-ex]].
  * [2] S. Chatrchyan et al. [CMS], Phys. Lett. B 716, 30-61 (2012) doi:10.1016/j.physletb.2012.08.021 [arXiv:1207.7235 [hep-ex]].
  * [3] P. W. Higgs, Phys. Rev. Lett. 13, 508-509 (1964) doi:10.1103/PhysRevLett.13.508
  * [4] F. Englert and R. Brout, Phys. Rev. Lett. 13, 321-323 (1964) doi:10.1103/PhysRevLett.13.321
  * [5] G. S. Guralnik, C. R. Hagen and T. W. B. Kibble, Phys. Rev. Lett. 13, 585-587 (1964) doi:10.1103/PhysRevLett.13.585
  * [6] B. Patt and F. Wilczek, [arXiv:hep-ph/0605188 [hep-ph]].
  * [7] M. Fabbrichesi, E. Gabrielli and G. Lanfranchi, doi:10.1007/978-3-030-62519-1 [arXiv:2005.01515 [hep-ph]].
  * [8] M. Graham, C. Hearty and M. Williams, Ann. Rev. Nucl. Part. Sci. 71, 37-58 (2021) doi:10.1146/annurev-nucl-110320-051823 [arXiv:2104.10280 [hep-ph]].
  * [9] S. Baek, J. Kim and P. Ko, [arXiv:2204.04889 [hep-ph]].
  * [10] A. Datta, A. Hammad, D. Marfatia, L. Mukherjee and A. Rashed, JHEP 03, 108 (2023) doi:10.1007/JHEP03(2023)108 [arXiv:2210.15662 [hep-ph]].
  * [11] S. Ghosh and P. Ko, [arXiv:2311.14099 [hep-ph]].
  * [12] D. Curtin, R. Essig, S. Gori, P. Jaiswal, A. Katz, T. Liu, Z. Liu, D. McKeen, J. Shelton and M. Strassler, et al. Phys. Rev. D 90, no.7, 075004 (2014) doi:10.1103/PhysRevD.90.075004 [arXiv:1312.4992 [hep-ph]].
  * [13] R. Foot, G. C. Joshi and H. Lew, Phys. Rev. D 40, 2487-2489 (1989) doi:10.1103/PhysRevD.40.2487
  * [14] P. Fileviez Perez and M. B. Wise, Phys. Rev. D 82, 011901 (2010) [erratum: Phys. Rev. D 82, 079901 (2010)] doi:10.1103/PhysRevD.82.079901 [arXiv:1002.1754 [hep-ph]].
  * [15] M. Duerr, P. Fileviez Perez and M. B. Wise, Phys. Rev. Lett. 110, 231801 (2013) doi:10.1103/PhysRevLett.110.231801 [arXiv:1304.0576 [hep-ph]].
  * [16] K. S. Babu, C. F. Kolda and J. March-Russell, Phys. Rev. D 54, 4635-4647 (1996) doi:10.1103/PhysRevD.54.4635 [arXiv:hep-ph/9603212 [hep-ph]].
  * [17] K. S. Babu, C. F. Kolda and J. March-Russell, Phys. Rev. D 57, 6788-6792 (1998) doi:10.1103/PhysRevD.57.6788 [arXiv:hep-ph/9710441 [hep-ph]].
  * [18] T. G. Rizzo, Phys. Rev. D 59, 015020 (1998) doi:10.1103/PhysRevD.59.015020 [arXiv:hep-ph/9806397 [hep-ph]].
  * [19] P. Ko and Y. Omura, Phys. Lett. B 701, 363-366 (2011) doi:10.1016/j.physletb.2011.06.009 [arXiv:1012.4679 [hep-ph]].
  * [20] P. Ko, Y. Omura and C. Yu, JHEP 01, 147 (2012) doi:10.1007/JHEP01(2012)147 [arXiv:1108.4005 [hep-ph]].
  * [21] P. Ko, Y. Omura and C. Yu, JHEP 04, 098 (2016) doi:10.1007/JHEP04(2016)098 [arXiv:1601.00586 [hep-ph]].
  * [22] P. Ko, Y. Omura and C. Yu, Phys. Lett. B 710, 197-200 (2012) doi:10.1016/j.physletb.2012.02.074 [arXiv:1104.4066 [hep-ph]].
  * [23] P. Gondolo, P. Ko and Y. Omura, Phys. Rev. D 85, 035022 (2012) doi:10.1103/PhysRevD.85.035022 [arXiv:1106.0885 [hep-ph]].
  * [24] B. Sheff, N. Steinberg and J. D. Wells, Phys. Rev. D 104, no.3, 036009 (2021) doi:10.1103/PhysRevD.104.036009 [arXiv:2008.10568 [hep-ph]].
  * [25] S. D. Lane, H. S. Lee and I. M. Lewis, [arXiv:2305.00013 [hep-ph]].
  * [26] M. Aaboud et al. [ATLAS], Phys. Rev. D 99, no.1, 012008 (2019) doi:10.1103/PhysRevD.99.012008 [arXiv:1808.10515 [hep-ex]].
  * [27] B. A. Dobrescu, G. L. Landsberg and K. T. Matchev, Phys. Rev. D 63, 075003 (2001) doi:10.1103/PhysRevD.63.075003 [arXiv:hep-ph/0005308 [hep-ph]].
  * [28] N. Toro and I. Yavin, Phys. Rev. D 86, 055005 (2012) doi:10.1103/PhysRevD.86.055005 [arXiv:1202.6377 [hep-ph]].
  * [29] S. D. Ellis, T. S. Roy and J. Scholtz, Phys. Rev. Lett. 110, no.12, 122003 (2013) doi:10.1103/PhysRevLett.110.122003 [arXiv:1210.1855 [hep-ph]].
  * [30] S. D. Ellis, T. S. Roy and J. Scholtz, Phys. Rev. D 87, no.1, 014015 (2013) doi:10.1103/PhysRevD.87.014015 [arXiv:1210.3657 [hep-ph]].
  * [31] D. Wang, L. Wu, J. M. Yang and M. Zhang, Phys. Rev. D 104, no.9, 095016 (2021) doi:10.1103/PhysRevD.104.095016 [arXiv:2102.01532 [hep-ph]].
  * [32] J. Ren, D. Wang, L. Wu, J. M. Yang and M. Zhang, JHEP 11, 138 (2021) doi:10.1007/JHEP11(2021)138 [arXiv:2106.07018 [hep-ph]].
  * [33] D. Wang, J. H. Cho, J. Kim, S. Lee, P. Sanyal and J. Song, Phys. Rev. D 109, no.1, 015017 (2024) doi:10.1103/PhysRevD.109.015017 [arXiv:2310.17741 [hep-ph]].
  * [34] X. Ai, W. Y. Feng, S. C. Hsu, K. Li and C. T. Lu, [arXiv:2401.15690 [hep-ph]].
  * [35] M. Bauer, M. Neubert and A. Thamm, Phys. Rev. Lett. 119, no.3, 031802 (2017) doi:10.1103/PhysRevLett.119.031802 [arXiv:1704.08207 [hep-ph]].
  * [36] A. Chakraborty, A. M. Iyer and T. S. Roy, Nucl. Phys. B 932, 439-470 (2018) doi:10.1016/j.nuclphysb.2018.05.019 [arXiv:1707.07084 [hep-ph]].
  * [37] C. T. Lu, Phys. Rev. D 108, no.11, 115029 (2023) doi:10.1103/PhysRevD.108.115029 [arXiv:2210.15648 [hep-ph]].
  * [38] G. Alonso-Álvarez, J. Jaeckel and D. D. Lopes, [arXiv:2302.12262 [hep-ph]].
  * [39] P. Draper and D. McKeen, Phys. Rev. D 85, 115023 (2012) doi:10.1103/PhysRevD.85.115023 [arXiv:1204.1061 [hep-ph]].
  * [40] [ATLAS], ATLAS-CONF-2012-079.
  * [41] K. Ban, K. Kong, M. Park and S. C. Park, [arXiv:2311.16674 [hep-ph]].
  * [42] A. Hammad, S. Moretti and M. Nojiri, JHEP 03, 144 (2024) doi:10.1007/JHEP03(2024)144 [arXiv:2401.00452 [hep-ph]].
  * [43] Attention is all you need, Vaswani, Ashish and Shazeer, Noam and Parmar, Niki and Uszkoreit, Jakob and Jones, Llion and Gomez, Aidan N and Kaiser, Łukasz and Polosukhin, Illia Advances in neural information processing systems -2017
  * [44] D. Feldman, Z. Liu and P. Nath, Phys. Rev. D 75, 115001 (2007) doi:10.1103/PhysRevD.75.115001 [arXiv:hep-ph/0702123 [hep-ph]].
  * [45] J. D. Wells, [arXiv:0803.1243 [hep-ph]].
  * [46] C. F. Chang, E. Ma and T. C. Yuan, JHEP 03, 054 (2014) doi:10.1007/JHEP03(2014)054 [arXiv:1308.6071 [hep-ph]].
  * [47] C. D. Carone and H. Murayama, Phys. Rev. D 52, 484-493 (1995) doi:10.1103/PhysRevD.52.484 [arXiv:hep-ph/9501220 [hep-ph]].
  * [48] A. Aranda and C. D. Carone, Phys. Lett. B 443, 352-358 (1998) doi:10.1016/S0370-2693(98)01309-4 [arXiv:hep-ph/9809522 [hep-ph]].
  * [49] S. Dittmaier et al. [LHC Higgs Cross Section Working Group], doi:10.5170/CERN-2011-002 [arXiv:1101.0593 [hep-ph]].
  * [50] J. J. Sakurai, Annals Phys. 11, 1-48 (1960) doi:10.1016/0003-4916(60)90126-3
  * [51] S. Tulin, Phys. Rev. D 89, no.11, 114008 (2014) doi:10.1103/PhysRevD.89.114008 [arXiv:1404.4370 [hep-ph]].
  * [52] [ATLAS], ATLAS-CONF-2020-027.
  * [53] K. Schmidt-Hoberg, F. Staub and M. W. Winkler, Phys. Lett. B 727, 506-510 (2013) doi:10.1016/j.physletb.2013.11.015 [arXiv:1310.6752 [hep-ph]].
  * [54] J. D. Clarke, R. Foot and R. R. Volkas, JHEP 02, 123 (2014) doi:10.1007/JHEP02(2014)123 [arXiv:1310.8042 [hep-ph]].
  * [55] J. Chang, K. Cheung, S. C. Hsu and C. T. Lu, Phys. Rev. D 95, no.3, 035012 (2017) doi:10.1103/PhysRevD.95.035012 [arXiv:1607.07550 [hep-ph]].
  * [56] B. M. K. Nefkens et al. [A2 at MAMI], Phys. Rev. C 90, no.2, 025206 (2014) doi:10.1103/PhysRevC.90.025206 [arXiv:1405.4904 [hep-ex]].
  * [57] S. Prakhov, B. M. K. Nefkens, C. E. Allgower, V. Bekrenev, W. J. Briscoe, J. R. Comfort, K. Craig, D. Grosnick, D. Isenhower and N. Knecht, et al. Phys. Rev. C 78, 015206 (2008) doi:10.1103/PhysRevC.78.015206
  * [58] P. A. Zyla et al. [Particle Data Group], PTEP 2020, no.8, 083C01 (2020) doi:10.1093/ptep/ptaa104
  * [59] R. Escribano, S. Gonzàlez-Solís and E. Royo, Phys. Rev. D 106, no.11, 114007 (2022) doi:10.1103/PhysRevD.106.114007 [arXiv:2207.14263 [hep-ph]].
  * [60] M. Ablikim et al. [BESIII], Phys. Rev. D 96, no.1, 012005 (2017) doi:10.1103/PhysRevD.96.012005 [arXiv:1612.05721 [hep-ex]].
  * [61] F. Ambrosino et al. [KLOE], Phys. Lett. B 681, 5-13 (2009) doi:10.1016/j.physletb.2009.09.022 [arXiv:0904.2539 [hep-ex]].
  * [62] R. R. Akhmetshin et al. [CMD-2], Phys. Lett. B 580, 119-128 (2004) doi:10.1016/j.physletb.2003.11.046 [arXiv:hep-ex/0310012 [hep-ex]].
  * [63] T. K. Pedlar et al. [CLEO], Phys. Rev. D 79, 111101 (2009) doi:10.1103/PhysRevD.79.111101 [arXiv:0904.1394 [hep-ex]].
  * [64] M. Ablikim et al. [BESIII], Phys. Rev. Lett. 122, no.14, 142002 (2019) doi:10.1103/PhysRevLett.122.142002 [arXiv:1902.03823 [hep-ex]].
  * [65] D. Babusci et al. [KLOE-2], Phys. Lett. B 720, 111-115 (2013) doi:10.1016/j.physletb.2013.01.067 [arXiv:1210.3927 [hep-ex]].
  * [66] P. Adlarson et al. [WASA-at-COSY], Phys. Lett. B 726, 187-193 (2013) doi:10.1016/j.physletb.2013.08.055 [arXiv:1304.0671 [hep-ex]].
  * [67] P. Ilten, Y. Soreq, M. Williams and W. Xue, JHEP 06, 004 (2018) doi:10.1007/JHEP06(2018)004 [arXiv:1801.04847 [hep-ph]].
  * [68] Y. L. Dokshitzer, G. D. Leder, S. Moretti and B. R. Webber, JHEP 08, 001 (1997) doi:10.1088/1126-6708/1997/08/001 [arXiv:hep-ph/9707323 [hep-ph]].
  * [69] M. Wobisch and T. Wengler, [arXiv:hep-ph/9907280 [hep-ph]].
  * [70] M. Cacciari, G. P. Salam and G. Soyez, JHEP 04, 063 (2008) doi:10.1088/1126-6708/2008/04/063 [arXiv:0802.1189 [hep-ph]].
  * [71] S. D. Ellis, T. S. Roy and J. Scholtz, Phys. Rev. Lett. 110, no.12, 122003 (2013) doi:10.1103/PhysRevLett.110.122003 [arXiv:1210.1855 [hep-ph]].
  * [72] A. Alloul, N. D. Christensen, C. Degrande, C. Duhr and B. Fuks, Comput. Phys. Commun. 185, 2250-2300 (2014) doi:10.1016/j.cpc.2014.04.012 [arXiv:1310.1921 [hep-ph]].
  * [73] J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, O. Mattelaer, H. S. Shao, T. Stelzer, P. Torrielli and M. Zaro, JHEP 07, 079 (2014) doi:10.1007/JHEP07(2014)079 [arXiv:1405.0301 [hep-ph]].
  * [74] R. Frederix, S. Frixione, V. Hirschi, D. Pagani, H. S. Shao and M. Zaro, JHEP 07, 185 (2018) [erratum: JHEP 11, 085 (2021)] doi:10.1007/JHEP11(2021)085 [arXiv:1804.10017 [hep-ph]].
  * [75] T. Sjostrand, S. Mrenna and P. Z. Skands, Comput. Phys. Commun. 178, 852-867 (2008) doi:10.1016/j.cpc.2008.01.036 [arXiv:0710.3820 [hep-ph]].
  * [76] J. de Favereau et al. [DELPHES 3], JHEP 02, 057 (2014) doi:10.1007/JHEP02(2014)057 [arXiv:1307.6346 [hep-ex]].
  * [77] M. Cacciari, G. P. Salam and G. Soyez, Eur. Phys. J. C 72, 1896 (2012) doi:10.1140/epjc/s10052-012-1896-2 [arXiv:1111.6097 [hep-ph]].
  * [78] [ATLAS], ATLAS-CONF-2020-026.
  * [79] S. Agostinelli et al. [GEANT4], Nucl. Instrum. Meth. A 506, 250-303 (2003) doi:10.1016/S0168-9002(03)01368-8
  * [80] M. Cacciari, G. P. Salam and G. Soyez, JHEP 04, 005 (2008) doi:10.1088/1126-6708/2008/04/005 [arXiv:0802.1188 [hep-ph]].
  * [81] [CMS], CMS-PAS-JME-14-001.
  * [82] P. T. Komiske, E. M. Metodiev, B. Nachman and M. D. Schwartz, JHEP 12, 051 (2017) doi:10.1007/JHEP12(2017)051 [arXiv:1707.08600 [hep-ph]].
  * [83] J. Arjona Martínez, O. Cerri, M. Pierini, M. Spiropulu and J. R. Vlimant, Eur. Phys. J. Plus 134, no.7, 333 (2019) doi:10.1140/epjp/i2019-12710-3 [arXiv:1810.07988 [hep-ph]].
  * [84] P. Berta, L. Masetti, D. W. Miller and M. Spousta, JHEP 08, 175 (2019) doi:10.1007/JHEP08(2019)175 [arXiv:1905.03470 [hep-ph]].
  * [85] B. Maier, S. M. Narayanan, G. de Castro, M. Goncharov, C. Paus and M. Schott, Mach. Learn. Sci. Tech. 3, no.2, 025012 (2022) doi:10.1088/2632-2153/ac7198 [arXiv:2107.02779 [physics.ins-det]].
  * [86] H. Qu, C. Li and S. Qian, [arXiv:2202.03772 [hep-ph]].
  * [87] T. Finke, M. Krämer, A. Mück and J. Tönshoff, JHEP 06, 184 (2023) doi:10.1007/JHEP06(2023)184 [arXiv:2303.07364 [hep-ph]].
  * [88] M. He and D. Wang, Eur. Phys. J. C 83, no.12, 1116 (2023) doi:10.1140/epjc/s10052-023-12293-1 [arXiv:2307.04723 [hep-ph]].
  * [89] A. Hammad and M. M. Nojiri, [arXiv:2404.14677 [hep-ph]].
  * [90] D. Bahdanau, K. Cho and Y. Bengio, [arXiv:1409.0473 [cs.CL]].
  * [91] J. Liu, N. Weiner and W. Xue, JHEP 08, 050 (2015) doi:10.1007/JHEP08(2015)050 [arXiv:1412.1485 [hep-ph]].
  * [92] F. Bezrukov and D. Gorbunov, JHEP 05, 010 (2010) doi:10.1007/JHEP05(2010)010 [arXiv:0912.0390 [hep-ph]].