Theory, phenomenology, and experimental avenues for dark showers: a Snowmass 2021 report

Abstract

In this work, we consider the case of a strongly coupled dark/hidden sector, which extends the Standard Model (SM) by adding an additional non-Abelian gauge group. These extensions generally contain matter fields, much like the SM quarks, and gauge fields similar to the SM gluons. We focus on the exploration of such sectors where the dark particles are produced at the LHC through a portal and undergo rapid hadronization within the dark sector before decaying back, at least in part and potentially with sizeable lifetimes, to SM particles, giving a range of possibly spectacular signatures such as emerging or semi-visible jets. Other, non-QCD-like scenarios leading to soft unclustered energy patterns or glueballs are also discussed. After a review of the theory, existing benchmarks and constraints, this work addresses how to build consistent benchmarks from the underlying physical parameters and present new developments for the pythia Hidden Valley module, along with jet substructure studies. Finally, a series of improved search strategies is presented in order to pave the way for a better exploration of the dark showers at the LHC.

1 Introduction

As the experimental program of the LHC searches for physics beyond the Standard Model (SM) is maturing, the community has started devoting significant effort to investigating alternative models and their associated phenomenology, especially those providing exotic signatures which would not have been directly addressed yet in the existing searches. Of interest here is the case of a strongly coupled dark sector or hidden sector, which extends the SM with an additional non-Abelian gauge group. Considering the non-trivial structure of the SM QCD, we should be open to the idea of a potentially complicated dark sector via non-Abelian gauge groups. These extensions generally contain matter fields, much like the SM quarks and gauge fields similar to the SM gluons. There are no a priory expectations on the gauge group dimension (number of colors), or that of matter fields (number of flavors) that the theory may have.

When the dark sector confines below some confinement scale ($\Lambda _\textrm{D}$) dark hadrons are formed, and depending on the symmetries of the theory, some of them could be stable leading to dark matter candidates. To allow for the production of dark states at the LHC, which could be either dark quarks or hadrons, the dark sector is coupled to the SM via a portal. The realisation of associated LHC phenomenology of such dark/hidden sector is however very much dependent on the details of the model. Nevertheless, some generic expectations can be set. For example, at the LHC, cases where dark quark masses ($m_{\textrm{q}_\textrm{D}}$) and corresponding confinement scale are much smaller than the collider centre-of-mass energy ($m_{\textrm{q}_\textrm{D}} \lesssim \Lambda _\textrm{D} \ll \sqrt{s}$) lead to spectacular signatures in terms of emerging or semi-visible jets [1,2,3]. Increasing $\Lambda _\textrm{D}$ implies heavier bound states, which in turn decreases the final state multiplicities for a given $\sqrt{s}$, as the allowed phase space decreases. This means that as the limit $\Lambda _\textrm{D} \sim \sqrt{s}$ is approached, depending on the relevant production mechanisms, $2\rightarrow 2$ SM initial state to dark meson final state processes become prevalent and resonance-like searches for dark bound states may prove useful [4,5,6,7]. Finally, cases where $m_{\textrm{q}_\textrm{D}} \gg \Lambda _\textrm{D}, m_{\textrm{q}_\textrm{D}} \lesssim \sqrt{s}$ lead to unusual signals known as quirks [8,9,10]. If the strongly-interacting sector is non-QCD like, other signatures such as Soft Unclustered Energy Patterns are also possible [9, 10]. For a review pertaining to this discussion see [6].

In this work, we focus on LHC exploration of such dark/hidden sectors where $\Lambda _\textrm{D} \ll \sqrt{s}$. In such cases, dark quarks could be produced at the LHC through a portal and undergo rapid hadronization within the dark sector before decaying back, at least in part and potentially with sizeable lifetimes, to SM particles. Models with such hidden dark sectors have been discussed e.g. in the context of twin Higgs models [11,12,13,14], composite and/or asymmetric dark matter scenarios [5, 15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32], and string theory [33, 34]. Some example studies focusing on hidden valley phenomenology can be found in [6, 10, 35,36,37,38,39,40,41,42,43,44,45,46,47,48] with additional examples referred to elsewhere in this write-up. The above references show a rather large activity throughout the last 4 decades.

These hidden valley models differ from most other Beyond the Standard Model scenarios because the infrared (IR) parameters of the theory can not be computed from ultraviolet (UV) definitions using perturbative techniques. This is in contrast to many other models e.g. MSSM, which have a well-defined relationships between UV and IR, even if their parameter space is high dimensional. The other well known approach to characterise new physics is the use of simplified models which do not rely on such top-down priors but are defined by their minimality; these may not be effective for hidden valley scenarios either, due to the inherent dependence on UV parameters in strongly interacting theories. Therefore, neither principles used otherwise to analyse new physics scenarios apply to hidden valley models.

While on the one hand these considerations motivate avenues for model-building with applications to shortcomings of the SM such as e.g. dark matter, LHC searches for hidden valleys are primarily motivated by the exotic phenomenology as stated above. Concretely, this means that we do not have a strong theory prior on e.g. the number of colors and flavors, the mass hierarchies amongst the matter fields, or the possible patterns of flavor breaking. Out of the multiple choices at hand, however, it is nevertheless possible to try and build internally coherent models and develop tools to predict their phenomenology and guide the searches.

Despite the complexity and the challenges in analysing such non-Abelian new sectors, there has been an increased activity in the recent years that has focused on understanding the signature parameter space of such models, both on the theory and the experimental sides. This report presents some of these developments, particularly concentrating on jet-like signatures, and puts in perspective the efforts necessary to make systematic progress in understanding, classifying and searching for such non-Abelian scenarios. Throughout this report we will consider dark sector scenarios where the dark quarks are uncharged under any SM group and the dark sector communicates to the SM via an additional mediator.

The report is organised as follows. QCD-like dark-sector scenarios are first reviewed in Sect. 2, addressing the theories in the s- and t-channels, and the existing benchmarks and limits for such models. Section 3 will address two possibilities of dark sector beyond the QCD-like scenarios: soft unclustered energy patterns (SUEP) and glueballs. Section 4 will be devoted to simulation tool limitations and how to build consistent benchmarks from the underlying physical parameters for semi-visible jets. After a discussion of consistent parameter setting, some improvements to the pythia Hidden Valley module and their validation will be presented, followed by some phenomenological studies on the jet substructure effects of varying the physical parameters. Finally, a series of improved search strategies will be discussed in Sect. 5, based on event-level variables, deep neural networks, autoencoder-based anomaly detection or better triggering algorithms.

2 QCD-like scenarios of dark sector

In SM QCD, the strong coupling constant $\alpha _s$ becomes weaker as the energy increases. This is known as asymptotic freedom. New non-Abelian sectors which display such asymptotic freedom fall into the category of QCD-like scenarios.

This section outlines possible exotic signatures that QCD-like dark sector scenarios may exhibit, including a discussion of benchmark models that have been or are currently being employed by the community. We also discuss existing limits in the context of these benchmarks. These signatures and results in turn provide a strong motivation behind detailed studies of such scenarios, motivating further theory effort, as will be discussed in Sect. 4.

2.1 Theories of dark QCD

Contributors: Timothy Cohen and Christiane Scherb

As an organizing principle, we will assume that the dark sector communicates with the Standard Model via a so-called portal. The Standard Model admits three renormalizable portals, in that there are three Standard Model gauge singlet operators with mass dimension less than four: the dark sector can couple to the field strength of the Hypercharge gauge boson $B_{\mu \nu }$ [49], to the Higgs bi-linear $|H|^2$ [50], or to the neutrino via HL [51]. In this paper, we will focus on introducing a new mediator particle that serves as the portal. This could be a $\textrm{Z}^{\prime }$ which would mediate s-channel production [52,53,54,55,56,57,58,59,60,61,62,63], as was proposed in the original hidden valley paper [36], or it could be a new scalar bi-fundamental which would mediate t-channel production [20, 21, 64,65,66,67]. In both of these cases, in the limit that the mediator mass is large, it may also be appropriate to integrate it out which would induce a contact operator. There are of course many options beyond these two examples, but to keep our scope finite we will only discuss these s- and t-channel production models.

Dark sector particles can then be produced at hadron colliders via the portal. Similar to SM quarks, dark quarks shower and hadronize and form dark jets. The properties of dark jets are determined by the dynamics of the dark sectors, namely the coupling strength, the ratio of unstable to stable dark hadrons inside the dark jets, and the mass scale of the dark hadrons. Production of dark sector particles at hadron colliders lead to a broad class of exotic signatures: Depending on the lifetime of the dark hadrons final states can contain semi-visible jets, lepton jets, emerging jets, soft bombs, quirks, etc. [1, 2, 9, 63, 68,69,70,71,72,73,74,75].

Our focus will be on characterizing exotic signatures that could result at the LHC. The space of possible dark sector models is vast, and furthermore many models can yield essentially the same LHC phenomenology. For this reason, we work with a simplified model-like parameterization of the dark sector. The phenomenology can be largely determined by specifying the dynamics of the dark sector shower (the number of dark colors, dark quark flavors, and the dark confinement scale), the mass spectrum, and decay patters of the dark mesons. We will largely frame the phenomenological implications of having a strongly coupled dark sector in terms of these variables.

2.1.1 s-channel

As discussed above, we take as an organizing principle that the dark sector communicates with the visible sector via a portal. If the portal is heavy, then one can describe it by integrating out the mediator to obtain a so called contact operator [76,77,78], for example

$$\begin{aligned} {\mathcal {L}} \supset \frac{c_{ij\alpha \beta }}{\Lambda }\left( \overline{q}_i\gamma ^\mu q_j\right) \left( \overline{\textrm{q}}_{\textrm{D} \alpha } \gamma _{\mu }\textrm{q}_{\textrm{D} \beta } \right) , \end{aligned}$$

(1)

where q are SM fermions, $\textrm{q}_\textrm{D}$ are dark sector quarks, $c_{ij\alpha \beta }$ are ${\mathcal {O}}(1)$ couplings encoding a possible flavor structure, and $\Lambda$ is the scale of the operator. Generally we use Roman indices as SM flavor indices and Greek indices for the dark sector flavor indices.

We are assuming that the portal couples the dark sector to the Standard Model quarks. Therefore, the observables of interest will be jets (which are expected to have non-QCD-like features) and missing energy that is likely to be aligned with the jets. One way to organize thinking about the possible signature space is in terms of the average fraction of invisible particles that are contained within a final state jet

$$\begin{aligned} r_{\text {inv}} \equiv \left\langle { \frac{\# \mathrm{{stable\, dark\, hadrons}}}{\# \mathrm{{dark\, hadrons}}}}\right\rangle . \end{aligned}$$

(2)

We will present results in terms of this variable in what follows.

One option for UV completing Eq. 1 is to introduce a so-called s-channel mediator. In such models, pairs of dark quarks can be produced via a heavy resonance $\textrm{Z}^{\prime }$, that also couples to SM quarks via [52,53,54,55,56,57,58,59,60,61,62,63]

$$\begin{aligned} {\mathcal {L}} \supset -Z_{\mu }^\prime \left( g_{\textrm{q}} \overline{q}_i\gamma ^\mu q_i + g_{\textrm{q}_\textrm{D}} \overline{\textrm{q}}_{\textrm{D} \alpha } \gamma ^\mu \textrm{q}_{\textrm{D} \alpha } \right) , \end{aligned}$$

(3)

where $g_{q,\textrm{q}_\textrm{D}}$ are the respective coupling constants. In general, $\textrm{Z}^{\prime }$ can also couple to other SM particles, which would lead to many possibilities in the final state. For concreteness here, we will focus on dark showers that result in SM jets + missing energy signatures. Therefore, we limit ourselves here to coupling the $\textrm{Z}^{\prime }$ to quarks, see Fig. 1. We will also simply give the $\textrm{Z}^{\prime }$ a mass, and will not worry about the associated Higgs mechanism or related effects. We will not discuss the additional particle content needed to cancel anomalies, nor the $Z-\textrm{Z}^{\prime }$ mixing structure needed for $g_{\textrm{q}_\textrm{D}} \ne g_{\textrm{q}}$ of models with a heavy $\textrm{Z}^{\prime }$ here (c.f. e.g. [79]), but will simply focus on the phenomenology.

Fig. 1

Diagram for the pair production of dark quarks through a $\textrm{Z}^{\prime }$ portal

Heavy resonances $\textrm{Z}^{\prime }$ are produced at a hadron collider in Drell Yan processes and will have a non-trivial branching ratio to decay to two dark quarks, which shower and hadronize in the dark sector. Then some of the dark sector hadrons are assumed to decay back to Standard Model quarks, which subsequently shower and hadronize as usual. Consequently, the phenomenology of s-channel models is governed by the following parameters: the $\textrm{Z}^{\prime }$ mass $m_{Z'}$, its couplings to visible and dark quarks $g_{\textrm{q}}$ and $g_{\textrm{q}_\textrm{D}}$, the dark sector shower (governed by the number of dark colors, dark flavors, and the scale of dark sector confinement $\Lambda _\textrm{D}$), the characteristic scale of the dark hadrons $m_\textrm{D}$, and the average fraction of stable hadrons that are aligned with the visible jet $r_{\text {inv}}$. While the coupling to SM quarks determines the $\textrm{Z}^{\prime }$ production cross section, the other parameters determine the final state. The details of the shower and mass scale of the dark hadrons determine how many dark sector particles are produced, and $r_{\text {inv}}$ determines the amount of missing energy in a dark jet. Depending on these parameters several interesting signatures and search methods can be defined, e.g. semi-visible jets and searches for dark matter in the jet substructure.

There are two dominant strategies to search for the signatures of these models, that largely depend on the choice of $r_{\text {inv}}$. When $r_{\text {inv}}$ is small, most of the final state associated with the resonance is visible, and so a normal bump hunt strategy can be employed. Then as $r_{\text {inv}}$ gets larger, it becomes advantageous to perform a bump hunt using the standard transverse mass variable . Then, once $r_{\text {inv}}$ approaches unity, the final state is essentially dominated by missing energy, and so a “mono-jet” style strategy is most sensitive. This is illustrated in Fig. 2 [63], where we show the projected limits on the s-channel model for these different strategies. This approach relies on very simple criteria, and so there is clearly much room for improvements that rely on additional characteristics of these models. Existing limits will be discussed in Sect. 2.3, while strategies for improvements will be addressed in Sect. 5.

Fig. 2

Figure taken from [63]

Estimates of the projected limits for the s-channel model. The ‘contact’ limit uses a mono-jet search strategy, where $\Delta \phi$ is the angle between the missing transverse energy and the closest jet.

2.1.2 t-channel

Another simple option for UV completing the contact operator is to introduce a so-called t-channel mediator. The t-channel UV completion is determined by the following interaction [20, 21, 64,65,66,67]

$$\begin{aligned} {\mathcal {L}}\supset -\left( \kappa _{\alpha i } \textrm{q}_{\textrm{D} \alpha } \Phi \overline{q}_{R i}\right) + h.c., \end{aligned}$$

(4)

Fig. 3

Diagram for the pair production of mediators and the subsequent decay to dark and SM quarks

which can be realized in either Minimal Flavor Violating models, where in addition to the SM flavor symmetry $U_q(3)\times U_u(3)\times U_d(3)$ the dark flavor symmetry $U_D(3)$ of the dark quarks $\textrm{q}_\textrm{D}$ is introduced [80,81,82], or by enhancing the SM gauge group by a dark flavor symmetry $SU_D(N_{c_\textrm{D}})$ and introducing $N_{f_\textrm{D}}$ dark quarks $\textrm{q}_\textrm{D}$ [1, 21, 83]. Then, the visible and dark sector communicate via a scalar bi-fundamental mediator $\Phi$ charged under both the SM and the dark flavor symmetry and $q_R$ represent right-handed up-type and down-type quarks. We consider the case where $\Phi$ is an SU(2)-singlet, and so it will not have any couplings to the left-handed quark doublets. (We note that generally SU(2)-doublet mediators coupling to left-handed SM quark doublets are also possible.) Depending on the hypercharge of $\Phi$, the dark sector communicates with either the up- ($Y_\Phi = 1/3$) or the down-type quarks ($Y_\Phi = -2/3$). In the following we will always use $N_{c_\textrm{D}} = N_{f_\textrm{D}} = 3$. In addition, we assume $m_{\textrm{q}_\textrm{D}} < \Lambda _\textrm{D}$, so that the pseudo-Nambu–Goldstone bosons (we will denote them dark pions in the following) of the spontaneously broken dark chiral symmetry are parametrically lighter than other dark hadrons. Consequently, heavier dark sector states, e.g. heavier dark hadrons or glueballs, will decay into dark pions and the dark pions will govern the phenomenology of such models. It is also worth pointing out that for $N_{f_\textrm{D}} >3$ an unbroken $SU(N_{f_\textrm{D}}-3)$ symmetry leads to one or more stable dark pion [83].

A particularly interesting feature of such models is the fact that the dark sector inherits the SM flavor structure via $\kappa _{\alpha \beta i j}$. The coupling $\kappa _{\alpha i}$ can generally be expressed as

$$\begin{aligned} \kappa = VDU \end{aligned}$$

(5)

with D a diagonal $3\times 3$ matrix of the form [80]

$$\begin{aligned} D=diag\left( \kappa _0+\kappa _1,\kappa _0+\kappa _2,\kappa _0-(\kappa _1+\kappa _2)\right) \end{aligned}$$

(6)

and V and U Hermitian $3\times 3$ matrices. For $m_{Q_{\alpha \beta }}=\delta _{\alpha \beta }m_{Q_{\alpha \beta }}$ the resulting dark flavor symmetry $U_d(3)$ can be used to rotate V away. Finally, U can be decomposed into

$$\begin{aligned} U = U_{23}U_{13}U_{12}, \end{aligned}$$

(7)

with $U_{ij}$ the rotational matrices for ij.

The phenomenology of t-channel dark sectors and dark mesons has been studied e.g. in [1, 4, 20, 21, 63,64,65,66,67, 80,81,82,83,84,85,86,87]. At colliders, the particle content of t-channel dark sector models can be produced in various channels such as from mediator pair production ($gg/qq \rightarrow \Phi \Phi$) or associated mediator production ($gq \rightarrow \Phi \textrm{q}_\textrm{D}$), as well as the direct production of dark quarks.

Here, we will focus on mediator pair production. The diagrams for this process are shown in Fig. 3. Both mediators decay subsequently to a visible and a dark quark, which undergo showering and hadronization, forming a SM and a dark jet. Heavier dark hadrons decay promptly into the lightest dark hadron, the dark pions. Dark pions decay back into SM particles. Depending on the lifetime of the dark pions three different final states are possible:

• The dark pions decay promptly and the final states consists of four prompt jets,

• The dark pions have intermediate lifetimes ($c\tau \sim 0.001{-}1$ m) and form emerging jets,

• The dark pions are stable on collider scales and are recorded as missing energy.

Fig. 4

Figure taken from [1]

Emerging jet signature for t-channel dark QCD models in comparison to the displaced dijet signature.

While in the first and last case the final states consists of typical SM objects, emerging jets provide a very distinct signature: Each pion of a dark jet will decay at a different length due to the boost of each dark pion depending on its individual momentum and the exponential distribution of the actual decay points for a given lifetime. Therefore, from a radial perspective, a dark jet deposits very little energy at the interaction point and then emerges with every dark pion decay into visible particles. The topology of an emerging jet is shown in Fig. 4 in comparison to the displaced dijet signature.

The emerging jet signature was first studied in [1] for a dark sector coupling to right-handed down-type quarks and a first search for this signature was performed with CMS [88]. In [87] this search has been combined with recasts of four jet and two jet plus missing energy. (For more details on the recast c.f. Sect. 2.3.5.) It was found that the dark pion mass does not change these bounds in a significant way. The obtained limits can be shown in the usual dark matter mass-mediator mass frame. To do so assumptions about the ratio of the dark pion and dark matter candidate, here taken as the dark proton ($p_D$), must been made.In Fig. 5 the exclusion limits, combined with the constraints from direct detection experiments are shown as a function of the dark matter mass and the mediator mass for $m_{p_D} = 10{m_{\pi _{\textrm{D}}}}$ for two different choices of couplings: The left panel corresponds to $\kappa = \text {diag}\left( 1,1,1\right)$, right to $\kappa = \text {diag}\left( 0.1,1,1\right)$. This shows clearly how the coupling structure will influence the best search method.

Fig. 5

Figure taken from [87]

Bounds from direct detection and collider searches for $m_{p_D} = 10{m_{\pi _{\textrm{D}}}}$. Left: $\kappa = \text {diag}(1,1,1)$, right: $\kappa = \text {diag}\left( 0.1,1,1\right)$. In the gray shaded region the dark matter is larger than the mediator mass.

Due to the connection to the SM flavor structure such models also contribute to flavor processes such as neutral meson mixing and flavor violating Kaon B and D decays. The impact on flavor physics for dark sectors coupled to the down-type quarks has been studied in [80, 82, 83], and for couplings to up-type quarks in [81, 82], as well as models with CP violation [89], and for a simplified model in [86]. In both cases the parameter space is largely unconstrained for dark pion masses above a few GeV. For the case of couplings to the up-type quarks this region could for example be probed via emerging jets from flavor violating top decays (c.f. [90]).

2.2 Existing benchmarks

Contributors: Elias Bernreuther, Florian Eble, Alison Elliot, Giuliano Gustavino, Simon Knapen, Benedikt Maier, Kevin Pedro, Jessie Shelton, Daniel Stolarski

Several attempts have been done in order to parametrize QCD-like dark sector theories in the literature. These are partly motivated by observations of parametric relationships between the confinement scale and rho and pion masses in the SM, and partly by inferred relationships between e.g. quark masses and meson masses within the pythia 8 HV module. We list below some of the efforts in this context. It is important to note that these do not necessarily imply consistent UV and IR parameters of the underlying non-Abelian dynamics itself, however they have been useful in providing first phenomenological insight in the behaviour of such theories to guide the experiments.

2.2.1 CMS emerging jet search

The CMS emerging jet search [88] follows the class of models introduced in Ref. [1]. The specific process investigated is ${{\textrm{pp}}} \rightarrow \Phi \overline{\Phi }$, $\Phi \rightarrow \textrm{q} \overline{\textrm{q}}_\textrm{D}$, depicted in Fig. 3, where $\Phi$ is a bifundamental scalar mediator with Yukawa couplings $\kappa _{\alpha i }$ between dark quarks $\textrm{q}_{\textrm{D} \alpha }$ and SM quarks $\textrm{q} _i$, as shown in Eq. 4.

The parameters of these models are briefly summarized here:

• $N_{c_\textrm{D}} = 3$

• $N_{f_\textrm{D}} = 7$

• $\Lambda _\textrm{D} = m_{\textrm{q}_\textrm{D}}$

• $m_{\textrm{q}_\textrm{D}} = 2{m_{\pi _{\textrm{D}}}}$

• $m_{\rho _\textrm{D}} = 4{m_{\pi _{\textrm{D}}}}$

• $m_{\Phi } = 400{-}2000\,\textrm{GeV}$

• ${m_{\pi _{\textrm{D}}}} = 1{-}10\,\textrm{GeV}$

• $c\tau _{\pi _{\textrm{D}}} = 1{-}1000~\text {mm}.$

The first two parameters are inspired by the dark matter model from [21]. These two along with the next two are pythia parameters that control the shower but are not directly observable. The mass of $\rho _\textrm{D}$ is set such that $\rho _\textrm{D} \rightarrow \pi _{\textrm{D}} \pi _{\textrm{D}}$ is kinematically allowed, and that decay will be dominant. The last three parameters in this list are treated as free parameters with their values varied as indicated.

The production cross section is controlled by $m_{\Phi }$. The distinct emerging jets phenomenology is controlled by the dark pion decay length $c\tau _{\pi _{\textrm{D}}}$, treated as a free parameter encapsulating variations of both the dark pion decay constant $f_{\pi _{\textrm{D}}}$ and the Yukawa coupling $\kappa _{\alpha i }$:

$$\begin{aligned}{} & {} c\tau _{\pi _{\textrm{D}}} \approx 80{\text { mm}} \left( \frac{1}{\kappa _{\alpha i }}\right) ^4 \left( \frac{2\,\textrm{GeV}}{f_{\pi _{\textrm{D}}}} \right) ^2 \left( \frac{100\,\textrm{MeV}}{m_{\textrm{q}}} \right) ^2\nonumber \\{} & {} \qquad \qquad \times \left( \frac{2\,\textrm{GeV}}{{m_{\pi _{\textrm{D}}}}} \right) \left( \frac{m_{\Phi }}{1\,\textrm{TeV}} \right) ^4, \end{aligned}$$

(8)

Because only pair production of $\Phi$ is considered, $\kappa _{\alpha i }$ does not influence the production cross section.

The event generation and hadronization are done with the pythia 8.212 Hidden Valley module, modified to allow running of the dark coupling constant $\alpha _\textrm{D}$. Only the Yukawa couplings to down quarks are non-zero. The dark rho mesons decay promptly to pairs of dark pions with branching fraction 0.999 or directly to pairs of down quarks with branching fraction 0.001, while the dark pions decay exclusively to pairs of down quarks with decay length $c\tau _{\pi _{\textrm{D}}}$. The parameter $P_{\rho _\textrm{D}}$, the probability of producing a vector rather than pseudoscalar meson during hadronization, is set to its default value of 0.75. Dark baryons are expected to be stable and can act as dark matter candidates as in the model of [21], but are expected to be produced rarely compared to dark mesons (${\sim }10\%$ for SM QCD [1]), so their presence is not simulated.

2.2.2 Flavored emerging jet model

The class of models described in Sect. 2.2.1 was extended to the case where the dark sector has a flavor structure related to SM QCD [83]. The result is multiple scenarios in which different dark mesons have different lifetimes, varying over a wide range of values. In particular, the following parameters are considered:

• $N_{c_\textrm{D}} = 3$

• $N_{f_\textrm{D}} = 3$

• $\Lambda _\textrm{D} > m_{\textrm{q}_\textrm{D}}.$

The last condition implies that heavier dark mesons decay promptly to lighter dark mesons, so only the latter influence the final state kinematic behavior. The spectrum of lighter dark mesons can be understood in analogy to SM pions and kaons: $\pi _{\textrm{D}} ^{\text {``}0\text {''}}$ ($\textrm{q}_{\textrm{D} \alpha } \overline{\textrm{q}}_{\textrm{D} \alpha }$), $\pi _{\textrm{D}} ^{\text {``}\pm \text {''}}$ ($\textrm{q}_\textrm{D 1} \overline{\textrm{q}}_\textrm{D 2}$, $\textrm{q}_\textrm{D 2} \overline{\textrm{q}}_\textrm{D 1}$), $\text {K}_\textrm{D} ^{\text {``}0\text {''}}$ ($\textrm{q}_\textrm{D 2} \overline{\textrm{q}}_\textrm{D 3}$, $\textrm{q}_\textrm{D 3} \overline{\textrm{q}}_\textrm{D 2}$), $\text {K}_\textrm{D} ^{\text {``}\pm \text {''}}$ ($\textrm{q}_\textrm{D 1} \overline{\textrm{q}}_\textrm{D 3}$, $\textrm{q}_\textrm{D 3} \overline{\textrm{q}}_\textrm{D 1}$); the quotation marks indicate that the superscripts do not represent actual charges. The mass splittings between these different dark meson species are taken to be negligible.

In the simplest version of this flavored model, called the “aligned” scenario, there is no mixing of neutral dark mesons. This scenario can be generated using a custom modification of pythia 8.230¹ that gives the correct proportions of the dark meson species listed above. The dark pion and kaon decay widths can be calculated as follows, for dark quark content $\overline{\textrm{q}}_{\textrm{D} \alpha } \textrm{q}_{\textrm{D} \beta }$ and SM quark products $\overline{\textrm{q}}_{i} \textrm{q}_{j}$:

$$\begin{aligned}{} & {} \Gamma _{\alpha \beta ij}=\frac{N_{c_\textrm{D}} {m_{\pi _{\textrm{D}}}} f_{\pi _{\textrm{D}}} ^{2}}{8\pi m_{\Phi } ^{4}}\left|\kappa _{\alpha i } \kappa ^{*}_{\beta j } \right|^2 \left( m_{\textrm{q}_{i}} ^{2}+m_{\textrm{q}_{j}} ^{2}\right) \nonumber \\{} & {} \quad \times \sqrt{\left( 1-\frac{(m_{\textrm{q}_{i}} +m_{\textrm{q}_{j}})^2}{{m_{\pi _{\textrm{D}}}}^{2}}\right) \left( 1-\frac{(m_{\textrm{q}_{i}}-m_{\textrm{q}_{j}})^2}{{m_{\pi _{\textrm{D}}}}^{2}}\right) }. \end{aligned}$$

(9)

2.2.3 CMS semi-visible jet search

The CMS search for semi-visible jets [91] is based on the class of models introduced in Refs. [2, 63]. The specific process investigated is ${\textrm{pp}} \rightarrow \textrm{Z}^{\prime } \rightarrow \textrm{q}_\textrm{D} \overline{\textrm{q}}_\textrm{D}$, where $\textrm{Z}^{\prime }$ is a leptophobic vector mediator with couplings to SM quarks $g_{\textrm{q}}$ and couplings to dark quarks $g_{\textrm{q}_\textrm{D}}$, as shown in Eq. 3.

The parameters of the models used in the CMS search are summarized below:

• $N_{c_\textrm{D}} = 2$

• $N_{f_\textrm{D}} = 2$

• $m_{\textrm{q}_\textrm{D}} = {m_{\pi _{\textrm{D}}}}/2$

• $m_{\textrm{Z}^{\prime }} = 1500{-}5100\,\textrm{GeV}$

• ${m_{\pi _{\textrm{D}}}} = m_{\rho _\textrm{D}} = 1{-}100\,\textrm{GeV}$

• $r_{\text {inv}} = 0.0{-}1.0$ (see Eq. 2)

• $\alpha _\textrm{D} = \alpha _\textrm{D} ^{\text {low}} {-}\alpha _\textrm{D} ^{\text {high}}.$

The last four parameters in this list are treated as free parameters with their values varied as indicated.

The unstable dark pions, as pseudoscalars, decay via a mass insertion preferentially to the most massive allowed SM quark species (bottom quarks unless ${{m_{\pi _{\textrm{D}}}} < 2m_{\textrm{b}}}$). The branching fractions for the mass insertion decays are calculated with quark mass running included. The unstable dark rho mesons, as vectors, decay democratically to pairs of any allowed SM quark. All decays to SM quarks are assumed to be prompt.

Fig. 6

Top left: the average number of dark hadrons, both stable and unstable, created per $\textrm{Z}^{\prime }$ event by pythia 8.230 for different $m_{\pi _{\textrm{D}}}$ and $\Lambda _\textrm{D}$ or $\alpha _\textrm{D}$ values. The behavior for ${\Lambda _\textrm{D} < m_{\textrm{q}_\textrm{D}}}$ may not be physically accurate. Top right: For each $m_{\pi _{\textrm{D}}}$ value, the corresponding $\Lambda _\textrm{D}$ value that maximizes the number of dark hadrons in the shower, $\Lambda _\textrm{D} ^{\text {peak}}$, is identified from the top left plot. The relationship between $m_{\pi _{\textrm{D}}}$ and $\Lambda _\textrm{D} ^{\text {peak}}$ is represented with a power law. Bottom: A comparison of different $\alpha _\textrm{D}$ variations for a representative dark hadron mass ${m_{\pi _{\textrm{D}}}} = 20~\textrm{GeV}$

The treatment of the dark coupling constant $\alpha _\textrm{D}$, or equivalently the dark scale $\Lambda _\textrm{D}$, follows a relationship that is derived based on the behavior shown in Fig. 6. The former can be expressed in terms of the latter: $\alpha _\textrm{D} (\Lambda _\textrm{D}) = \pi /\left( b_{0}\log \left( Q_\textrm{D}/\Lambda _\textrm{D} \right) \right)$, with the factor $b_{0} = (11N_{c_\textrm{D}}- 2N_{f_\textrm{D}})/6$ and $Q_\textrm{D} = 1\,\textrm{TeV}$. The benchmark value for the scale is chosen according to the empirical fit ${\Lambda _\textrm{D} ^{\text {peak}} = 3.2}{({m_{\pi _{\textrm{D}}}})^{0.8}}$. This maximizes the production of dark hadrons in the dark showers, because the effect of $\Lambda _\textrm{D}$ depends on the dark hadron mass $m_{\pi _{\textrm{D}}}$. The corresponding value $\alpha _\textrm{D} ^{\text {peak}}$ is then varied by ${\pm }50\%$ to give $\alpha _\textrm{D} ^{\text {low}}$, $\alpha _\textrm{D} ^{\text {high}}$. It is important to note that the results generated from pythia for ${\Lambda _\textrm{D} < m_{\textrm{q}_\textrm{D}}}$ may not be physically accurate. This is discussed further in Sect. 4.1.

The event generation and hadronization are done with the pythia 8.230 Hidden Valley module. The invisible fraction $r_{\text {inv}}$ is implemented by reducing the branching fractions of dark hadrons to SM quarks, so the “decay” of a dark hadron to invisible particles represents a stable dark hadron. A ${\mathbb {Z}}_{2}$ symmetry filter is enforced to reject events with an odd number of stable dark hadrons, as they would always be produced in pairs in a complete model. The exact pythia settings used in Ref. [88] can be found on HEPData [92].

2.2.4 Further semi-visible jet models under study by CMS

The published CMS searches described in Sects. 2.2.1 and 2.2.3 use pythia for both event generation and hadronization. The processes that the pythia Hidden Valley module can generate are those shown in Figs. 1 and 3.

Fig. 7

Processes available via FeynRules in MadGraph5: ${\textrm{pp}} \rightarrow \textrm{g} \textrm{Z}^{\prime } \rightarrow \textrm{g} \textrm{q}_\textrm{D} \overline{\textrm{q}}_\textrm{D}$ (boosted), ${\textrm{qg}} \rightarrow \textrm{q}_\textrm{D} \Phi \rightarrow \textrm{q}_\textrm{D} \textrm{q} \overline{\textrm{q}}_\textrm{D}$, $\textrm{q} \overline{\textrm{q}} \rightarrow \textrm{q}_\textrm{D} \overline{\textrm{q}}_\textrm{D}$ (non-resonant)

Studies are in progress to extend the CMS semi-visible jet program to low-mass boosted resonances and t-channel production via the bifundamental scalar mediator $\Phi$. To generate these processes, depicted in Fig. 7, MadGraph5 2.6.5 is used. (MadGraph5 can also generate the processes in Figs. 1 and 3.) The FeynRules definitions are obtained from Ref. [93], associated with Ref. [63]. The pythia Hidden Valley module is still used for hadronization. To obtain accurate results, several additional steps are needed:

• Ensure that mediator decay widths are properly computed (set param_card decay [id] auto)

• Increase the number of events when making gridpacks from 2000 to 10,000 (to overcome instability in t-channel phase space integration)

• Convert PDG IDs to pythia conventions in LHE output.

Both MadGraph5 and pythia have some limitations for handling dark shower generation and hadronization. For pythia, some useful future additions would be:

• Add processes like $\textrm{qq} \rightarrow \textrm{q}_\textrm{D} \Phi \rightarrow \textrm{q}_\textrm{D} \textrm{q} \overline{\textrm{q}}_\textrm{D}$, $\textrm{gg} \rightarrow \textrm{q}_\textrm{D} \overline{\textrm{q}}_\textrm{D}$ (non-resonant) that are currently only available via FeynRules, to facilitate generator-level studies.

• Include more theory uncertainties as event weights: PDF variations, renormalization and factorization scale variations, and Hidden Valley-specific parameters such as those that control the hadronization models.

• Allow user control over the dark hadron spectrum for studies of flavored models (Sect. 2.2.2, Ref. [83]).

• Add dark baryons for completeness (currently only dark mesons are explicitly produced).

Updates to the pythia Hidden Valley module made in the context of this Snowmass project and described in Sect. 4.2, are intended to address the last two points. For MadGraph5, some useful updates and improvements would include:

• Better fixes for the items mentioned in the previous paragraph.

• Central support for common processes such as t-channel production, to reduce the need for manual FeynRules implementations.

• Ability to add new SU(N) gauge groups in a complete and consistent way, in order to model Hidden Valley radiation explicitly. A Hidden Valley jet matching procedure would also be needed to avoid overlap with radiation from pythia during the hadronization process.

2.2.5 Semi-visible jet models under study by ATLAS

The semi-visible jet models under study are based on those introduced in Refs. [2, 63]; the fully visible jet models described in [72] are also considered. The $\textrm{Z}^{\prime }$ and $\Phi$ portals introduced in Sects. 2.1.1 and 2.1.2 are considered, with the following parameter settings:

• $N_{c_\textrm{D}} = 2$

• $N_{f_\textrm{D}} = 2$

• $m_{\textrm{q}_\textrm{D}} = {m_{\pi _{\textrm{D}}}}/2$

• ${m_{\pi _{\textrm{D}}}} = m_{\rho _\textrm{D}} = 20\,\textrm{GeV}$

• $m_{\textrm{Z}^{\prime }}/ m_{\Phi } = 750{-}5000\,\textrm{GeV}$

• $r_{\text {inv}} = 0.0{-}1.0.$

The $r_{\text {inv}}$ and $m_{\textrm{Z}^{\prime }}$ parameters are varied through the values as displayed in the above list. The $\alpha _\textrm{D}$ coupling is chosen to be a running constant as it is in SM QCD.

The generated models differ between s-channel and t-channel in terms of the jet matching applied. For the t-channel, an MLM matching scheme is used, while for the s-channel, a CKKM-L matching has been used.

The semi-visible jet models are generated in MadGraph5, and pythia  is used for showering. Fo the s-channel, MadGraph5  2.9.3 and pythia  8.245 are used. For the t-channel, MadGraph5  2.8.1 and pythia  8.244 are used. The fully visible jet models are entirely generated with pythia  8.

2.2.6 Aachen model

This model was introduced in Ref. [61] and subsequently used in Refs. [73, 94]. It is designed to satisfy cosmological constraints and reproduce the observed dark matter relic abundance. The dark sector is connected to the SM by a heavy vector mediator $\textrm{Z}^{\prime }$ arising from a new $U(1)^{\prime }$ gauge group.

The baseline parameters of this model are summarized below:

• $N_{c_\textrm{D}} = 3$

• $N_{f_\textrm{D}} = 2$

• $\Lambda _\textrm{D} = m_{\rho _\textrm{D}}$

• $m_{\textrm{q}_\textrm{D}} = 0.5\,\textrm{GeV}$

• $m_{\textrm{Z}^{\prime }} = 1\,\textrm{TeV}$

• ${m_{\pi _{\textrm{D}}}} = 4\,\textrm{GeV}$

• $m_{\rho _\textrm{D}} = 5\,\textrm{GeV}$

• $r_{\text {inv}} = 0.75.$

The particle masses $m_{\textrm{q}_\textrm{D}}$, $m_{\textrm{Z}^{\prime }}$, ${m_{\pi _{\textrm{D}}}}$, and $m_{\rho _\textrm{D}}$ are treated as free parameters that are set to the benchmark values indicated above. In the benchmark model, only the $\rho _\textrm{D} ^{0}$ mesons decay (where 0 indicates the $U(1)^{\prime }$ charge), while all other dark mesons are stable on the scale of the detector. Dark baryons and other dark bound states, such as dark eta and dark omega mesons, are assumed to be too heavy to be produced frequently. The study in Ref. [73] includes variations ${m_{\pi _{\textrm{D}}}} = m_{\rho _\textrm{D}} = 5{-}20\,\textrm{GeV}$ and $r_{\text {inv}} = 0.1{-}0.9$.

MadGraph5 2.6.4 with FeynRules 2.3.13 is used to generate the studied events, with pythia 8.240 used for hadronization. For the benchmark case with $r_{\text {inv}} = 0.75$, the pythia Hidden Valley default settings are modified to $P_{\rho _\textrm{D}} = 0.5$ in order to obtain the correct proportion of unstable dark mesons. When $r_{\text {inv}}$ is varied, it is implemented in the same manner described in Sect. 2.2.3. In this case, all unstable dark hadrons are assumed to decay democratically to pairs of $\textrm{u}$, $\textrm{d}$, $\textrm{s}$, or $\textrm{c}$ quarks. (The assumption that only $\rho _\textrm{D} ^{0}$ mesons decay is relaxed.)

2.2.7 Decay portals

Hidden valley models differ from most other models in two crucial aspects. The first is the lack of a very well-defined theory prior on what the model should look like, as already mentioned earlier in this document. The second complication is that, even if we somehow had a strongly preferred model, it would still be very challenging to extract accurate predictions for all observables due to the non-perturbative dynamics in the dark sector.

Both these challenges indicate that we may be best served by a suite of searches that is as model-independent as possible. This is easier said than done however, and some theory priors are always needed to design an experimental analysis, especially in the initial phases of this program. Reference [95] has advocated to inject these theory priors in the decay portals that allow the dark sector to decay back to the Standard Model. This approach has the following advantages:

• 1. The number of plausible options is relatively limited once one restricts to the set of decay portals that do not introduce dangerous flavor-changing neutral currents. A systematic survey is therefore very feasible.

• 2. We have good theoretical control over these decay portals, in terms of both the dark particle’s decay length and its allowed branching ratios. This is in contrast to the process of dark sector hadronization, where we must resort to parameterizing our ignorance as best we can.

• 3. If one moreover insists on a relatively minimal UV completion and/or no more than moderate fine tuning, one can moreover derive an approximate lower bound on the lifetime of the dark mesons. This gives the models a little bit more predictive power. The resulting lifetimes and branching ratios of the visibly-decaying dark particle are crucial for the design of experimental search strategies and a systematic survey of models featuring a minimal suite of decay portals can therefore be a good starting point for a comprehensive experimental search program.

Table 1 Overview of the decay portals considered in [95]. The $\eta _\textrm{D}$ and $\omega _\textrm{D}$ represent respectively the lightest spin-zero and spin-one meson in the dark sector, and $F'_{\mu \nu }$ is the field strength for an elementary dark photon $\textrm{A}^{\prime }$. The decay portal column indicates the operator(s) that allow the unstable dark meson to decay, defining the model

Fig. 8

Approximate lower bound on the proper lifetime of the visibly decaying particle (VDP) in the dark sector, for the decay portals in Table 1

The portals considered in Ref. [95] are summarized in Table 1. To maintain compatibility with the pythia Hidden Valley module at the time, a simplified dark sector was considered with a single (pseudo)scalar $\eta _\textrm{D}$ and vector $\omega _\textrm{D}$ meson. For the dark photon portal an additional elementary dark photon $\textrm{A}^{\prime }$ was added. For each of these portals the branching ratios and lower bound on the lifetime of the visibly decaying particle were computed (see Fig. 8).

For the gluon portal one assumes that the dark sector pseudoscalar meson ($\eta _\textrm{D}$) decays through a dimension 5 coupling to the Standard Model gluons, leading to very hadron-rich final states. This coupling requires a low-scale UV completion, especially because $\eta _\textrm{D}$ itself is a composite state and the $\eta _\textrm{D} G\tilde{G}$ coupling should be suppressed by dimensional transmutation in models with perturbative parton showers. Such a UV completion moreover requires the presence of new, colored particles, which could be produced at the LHC. The bound in Fig. 8 was obtained by assuming that any such states must have a mass $\gtrsim 2\,\textrm{TeV}$. This assumption can be relaxed by devising an elaborate extension of the model to hide the new colored particles from existing searches at ATLAS and CMS.

The photon portal works in the same way as the gluon portal, except that the $\eta _\textrm{D}$ couples to the Standard Model photons instead of the gluons. This portal leads to very photon-rich final states. Its lifetime bound is informed by collider constraints on new charged particles, which are required in the UV completion of this portal.

The vector portal assumes that the Standard Model photon mixes with the dark vector meson $\omega _\textrm{D}$, similar to the photon-$\rho$ mixing in the Standard Model. In this portal it was assumed that the $\eta _\textrm{D}$ was absolutely stable, leading to a semi-visible jet phenomenology. The UV completion of this portal requires the introduction of an additional, elementary $\textrm{Z}^{\prime }$ vector field. The lifetime bound in Fig. 8 was informed by the bounds on the mass and coupling of the $\textrm{Z}^{\prime }$ from direct searches and electroweak precision observables.

In the Higgs portal scenario, the $\eta _\textrm{D}$ decays by mixing with the Standard Model Higgs, leading to a heavy flavor-rich phenomenology. In this scenario, no new particles are needed in the UV completion. The $H^\dagger H$ operator does however contribute directly to the masses of the constituents of the $\eta _\textrm{D}$ mesons, which leads to the lower bound in Fig. 8. This bound can be evaded by fine-tuning the masses of the constituents of the $\eta _\textrm{D}$.

Finally, the dark photon portal is inspired by the Standard Model $\eta \rightarrow \gamma \gamma$ decay, as it assumes a light, elementary vector field ($\textrm{A}^{\prime }$) which couples to the $\eta _\textrm{D}$ through a chiral anomaly. The $\textrm{A}^{\prime }$ itself can then mix with the Standard Model photon, resulting fairly lepton-rich final states. With this portal, the lifetime constraint is very mild and comes exclusively from direct searches for the $\textrm{A}^{\prime }$ itself.

All assumptions and branching ratios are encoded in the public python script https://gitlab.com/simonknapen/dark_showers_tool which can generate pythia configuration cards for all five decay portals.

2.3 Existing constraints

Contributors: Giuliano Gustavino, Steven Lowette, Kevin Pedro, Pedro Schwaller, Andrii Usachov, Carlos Vázquez Sierra

The discussion so far has not only established that QCD-like dark sectors are theoretically interesting, but that they can also lead to exotic signatures for the experiments. In light of this observation, several searches at the LHC have already been carried out, with many other studies also underway. The search results are currently reported using one or more of the benchmarks discussed in Sect. 2.2. In this section, we illustrate some of the public search results, and existing constraints on the theory landscape. Given that the searches rely on generic signal characteristics, it may be possible to reinterpret the results in terms of other, UV/IR coherent models.

2.3.1 ATLAS search program

While no direct constraints have been published so far by the ATLAS Collaboration on these models, a broad set of semi-visible jet scenarios in both t- and s-channel production processes are being studied by the collaboration as mentioned in the previous section. Besides possible dedicated searches, the recasting of previously published results in other channels might prove useful in constraining the parameter space. Indeed, existing exclusion limits obtained in analyses looking at di-jet [96] and +jet [97] final states should already be able to constrain a phase-space predicted by dark QCD models in different $r_{\text {inv}}$ ranges. Furthermore, some studies also focus on scenarios with emerging jets, where dark hadrons have a non-negligible lifetime.

2.3.2 CMS search for emerging jets

The CMS emerging jet search [88] follows the class of models introduced in Ref. [1], as described in Sect. 2.2.1. The search considers the following parameter variations: $m_{\Phi } = 400{-}2000\,\textrm{GeV}$, ${m_{\pi _{\textrm{D}}}} = 1{-}10\,\textrm{GeV}$, $c\tau _{\pi _{\textrm{D}}} = 1{-}1000{\text { mm}}$.

The search requires four high-$p_{{\textrm{T}}}$ jets and triggers on such events using the scalar sum of their momenta, $H_{{\textrm{T}}}$. Per-jet quantities indicating displaced tracks – including the 2D impact parameter, the 3D impact parameter significance, and the fraction of the $p_{{\textrm{T}}}$ from prompt tracks associated with the primary vertex – are used to identify or “tag” jets as emerging. The signal region definition requires either two jets tagged as emerging, or one jet tagged as emerging along with substantial missing transverse momentum. The latter option increases sensitivity to models with larger $c\tau _{\pi _{\textrm{D}}}$ values where many dark hadrons decay outside of the detector. The misidentification rate for this tagging procedure is measured and used to estimate the QCD multijet background. Heavy flavor jets from B hadrons are found to be misidentified more frequently, as expected because of their non-negligible decay lengths. Multiple signal regions with different selection requirements are defined, and for each model, the signal region with the highest expected sensitivity is used.

Fig. 9

Reproduced from Ref. [88]

95% CL limits on the cross section for pair production of $\Phi$ leading to emerging jets, in the plane of $c\tau _{\pi _{\textrm{D}}}$ versus $m_{\Phi }$ (here called $\text {c}\tau _{\pi _\text {DK}}$ and ${\textrm{m}}_{{\textrm{X}}_{\textrm{DK}}}$) where $m_{\pi _{\textrm{D}}}$ (here called $m_{\pi _\text {DK}}$) is set to 5 GeV.

The search results are shown in Fig. 9. Using a $13\,\textrm{TeV}$ dataset with $16.1\,\textrm{fb}^{-1}$ of integrated luminosity, no significant excess above the SM prediction is observed. This search excludes models with $400< m_{\Phi } < 1250\,\textrm{GeV}$ for $5< c\tau _{\pi _{\textrm{D}}} < 255{\text { mm}}$ at 95% confidence level (CL). It is found that the limits are similar over the range of $m_{\pi _{\textrm{D}}}$ values explored. Work is ongoing to incorporate the remainder of the LHC Run 2 dataset, up to $138\,\textrm{fb}^{-1}$, and to improve the sensitivity to other models, such as Ref. [83], which is summarized in Sect. 2.2.2.

2.3.3 CMS search for semi-visible jets

The CMS search for semi-visible jets [91] is based on the class of models introduced in Refs. [2, 63], as described in Sect. 2.2.3.

The search requires two high-$p_{{\textrm{T}}}$ wide jets and triggers on the jet $p_{{\textrm{T}}}$ and the $H_{{\textrm{T}}}$. This dijet system is combined with the missing transverse momentum to compute the transverse mass $M_{{\textrm{T}}}$, which has a falling spectrum for the SM backgrounds, while the signal has a kinematic edge at the $\textrm{Z}^{\prime }$ mediator mass. The QCD multijet background is rejected by requiring high values of the transverse ratio , while the electroweak backgrounds ($\textrm{t}$ $\overline{\textrm{t}}$, ${\textrm{W}}{(\rightarrow \ell \nu )\text {+jets}}$, ${\text {Z}}{(\rightarrow \nu \overline{\nu })\text {+jets}}$) are rejected by vetoing identified and isolated leptons ($\textrm{e}$, $\mu$) and requiring a small minimum angle between the jets and the . Various sources of instrumental background and misreconstruction, which introduce artificial , are also rejected. Two signal regions are defined in terms of the $R_{{\textrm{T}}}$ variable to provide additional sensitivity.

Fig. 10

Reproduced from Ref. [91]

95% CL limits on the cross section for $\textrm{Z}^{\prime } \rightarrow \text {semi-visible jets}$ in a two-dimensional plane with variations of $m_{\textrm{Z}^{\prime }}$ and $r_{\text {inv}}$ from the model-independent (left) and model-dependent (right) searches, where ${m_{\pi _{\textrm{D}}}}=m_{\rho _\textrm{D}} =20~\textrm{GeV}$ (here called ${\textrm{m}}_{\hbox {dark}}$) and $\alpha _\textrm{D} =\alpha _\textrm{D} ^{\text {peak}}$ (here called $\alpha _\text {dark}$).

The full $13\,\textrm{TeV}$ dataset of $138\,\textrm{fb}^{-1}$ is analyzed and no indication of a resonance in the $M_{{\textrm{T}}}$ spectrum is observed. Both model-independent and model-dependent results are obtained, as shown in Fig. 10. The former excludes models with $1.5< m_{\textrm{Z}^{\prime }} < 4.0\,\textrm{TeV}$ and $0.07< r_{\text {inv}} < 0.53$ at 95% CL, depending on the other model parameters. The latter uses a boosted decision tree (BDT) that combines jet substructure variables to tag jets as semi-visible. It extends the exclusions to $1.5< m_{\textrm{Z}^{\prime }} < 5.1\,\textrm{TeV}$ and $0.01< r_{\text {inv}} < 0.77$ for the specific models described in Sect. 2.2.3 that were used to train the BDT.

2.3.4 CMS search for SIMPs as a link to signatures of trackless jets

Dark sector models could give rise to experimental signatures where jets are formed with visible regular hadrons arising from decays of hidden-sector particles that are long-lived, and thus make these hadrons appear displaced within the jet. With sufficient displacement, jets can arise in which none of the tracks of the constituent charged hadrons can be reconstructed, thus making the jets appear neutral. Such neutral jets are extremely rare among high-momentum jets from regular standard-model quarks or gluons, and can thus be a very sensitive probe of physics beyond the standard model from dark sectors.

The CMS collaboration has performed a search for a pair of such trackless jets [98] using $16.1\,\textrm{fb}^{-1}$ of integrated luminosity recorded in 2016. The signature probed consisted of a pair of back-to-back high-momentum trackless jets, where experimentally the trackless nature was sought by looking for the ratio of the jet energy carried by charged particles to the energy carried by neutral particles to be less than 5%. To illustrate the effectiveness of this requirement in suppressing standard model QCD background jets, a background rejection of over $10^5$ in data was reported for this 5% requirement.

This CMS search for trackless jets was inspired by and interpreted in a model proposing a new interaction through a low-mass mediator with a new dark matter fermion [3]. The interaction leads to very high interaction cross sections, which are not necessarily excluded, though many model assumptions need to be made to avoid the many cosmology, particle physics and astrophysical constraints. The model considered is not a dark sector model per se, as indeed there is no decay back from a dark sector leading to missing charged hadrons, but rather the jets constitute solely of SIMP particles interacting in the calorimeters, thus generating neutral jets. The aspect of a displaced decay is thus missing, though the similar experimental signature does potentially impact dark sector searches.

The interpretation of the CMS trackless jets search in the SIMPs model has been found to be difficult at large SIMP masses, above ${\sim }100\,\textrm{GeV}$. At such high masses, the modeling of the SIMP-nucleon interaction is complicated by the SIMP mass, and the approximation used in the CMS analysis, which involved treating the SIMP in the Geant simulation as a massive neutron-like object, becomes exploratory. As such, the strong exclusion limits on the SIMP pair production cross section as a function of the SIMP mass, obtained by CMS in absence of an excess of trackless jets in the analyzes data, shown in Fig. 11, are reported with this caveat above $100\,\textrm{GeV}$.

Fig. 11

Figure taken from [98]

Cross section upper limits for the production of a pair of SIMPs (represented by $\chi$) in the context of the CMS search for trackless jets. See text for details.

2.3.5 Collider constraints on t-channel models

A plethora of new physics searches are ongoing at the LHC and have pushed the limits on the masses of new particles above the TeV scale in many cases. While dark showers are a spectacular and unique signature, existing new physics searches still retain some sensitivity in regions of parameter space where the signal looks SM-like.

For t-channel dark sectors [1, 21, 83], often the leading production mode is pair production of the heavy mediator, which gives rise to a signature with two ordinary QCD jets and two dark showers, as shown in Fig. 3. It is clear that missing energy searches should become efficient in the limit where the dark sector particles become very long lived, while searches for prompt multijet signals can probe the regime of very short lifetimes. Due to the stochastic nature of particle decays, these searches also retain some sensitivity in the intermediate lifetime regime.

Fig. 12

Figure taken from [87]

Constraints on the mass of a t-channel mediator $\Phi$ (here called X) as a function of the dark sector particle lifetime. See text for details.

A recast of a di-jet plus MET search [99], a search for paired prompt di-jet resonances [100] and of the dedicated emerging jets search [88] was performed in Ref. [87]. As can be seen in Fig. 12, the dedicated emerging jets search performs best in the intermediate lifetime regime, while the recast searches can probe all the remaining range of lifetimes. As expected, mediator masses below the TeV scale are already strongly constrained. In the short lifetime regime, the constraints are weaker. This is partially because the published search uses only a limited amount of data, but also because fighting the QCD multi-jet background is hard. Another option to constrain the t-channel mediators in the regime of short dark sector lifetimes is through their contribution to angular correlations in dijet events [101], when the mediator is exchanged in the t-channel instead. This constraint however will depend on the magnitude of the Yukawa coupling of the mediator.

The t-channel mediators naturally allow couplings that break SM flavor symmetries. While most flavor violating couplings are constrained by low-energy flavor observables [83, 86], the LHC has the potential to constrain flavor changing couplings of the top quark [90]. Such scenarios could also be of interest for searches for dark showers produced in association with tops or even stemming from exotic top decays.

2.3.6 Existing constraints and projections from LHCb

The LHCb experiment, originally designed for heavy-flavor physics, has shown its potential as a general purpose detector in the recent years. An excellent secondary vertex resolution, low momentum thresholds and particle identification capabilities make LHCb a natural candidate to search for dark QCD signatures in the low-mass region. These searches are now becoming part of the LHCb physics program, where world-leading constraints have already been set for hidden valley scenarios, with reinterpretations in other dark sector models as well as a large number of very encouraging sensitivity projections, described in the following paragraphs.

LHCb has published a search for low-mass resonances decaying into pairs of muons, using 5.1 fb$^{-1}$ of data collected at 13 TeV [102]. In this article, a model-independent search of both prompt and displaced resonances X is performed. For the displaced case, the secondary $X~{\rightarrow }~\mu ^+\mu ^-$ decay vertex is required to be transversely displaced from the primary vertex in the range $12<\rho _T<30$ mm, allowing for the resonance to become a long-lived signature in the detector. Then, limits on the cross-section $\sigma (X~{\rightarrow }~\mu ^+\mu ^-)$ are placed, and interpreted in various production models. One of these models in the regard of dark QCD is that of a hidden valley scenario, where constraints are set on the kinetic mixing strength $\gamma -Z_{HV}$ between a heavy hidden valley boson $Z_{HV}$ with photon-like couplings, and a photon, fixing the average multiplicity of hidden valley hadrons to a value of 10. These are the most stringent constraints placed up to date, for masses of a composite hidden valley vector boson, X, up to 3 GeV, as presented in Fig. 13.

Fig. 13

Figure taken from Ref. [102]

Upper bounds at 90% CL on the kinetic mixing strength $\gamma -Z_{HV}$. The grey box shows a vetoed region due to the large doubly misidentified $K_S^0$ background.

These results have been also re-interpreted in the context of Z-initiated dark showers, assuming various benchmark scenarios [103]. In Fig. 14, projections for one of the scenarios are shown, showing the capabilities of LHCb to probe Z branching fractions down to 10$^{-7}$ during the high-luminosity phase.

Fig. 14

Projections at 90% CL of LHCb sensitivity to a muon rich dark showers initiated by the decay of a Z into dark quarks (here called $\psi '$). In this benchmark scenario, a dark pion mass (here called $\hat{\pi }$) of 650 MeV, an average multiplicity of 7 and a branching fraction of the dark pion into two muons of 96%, are assumed. More details can be found in Ref. [103]

In a more general sense, the capabilities of LHCb to probe other dark QCD models are summarized in Ref. [104], in the context of benchmark scenarios featuring a range of dark hadron and mediator masses, for different assumptions on the average dark hadron multiplicity in the dark sector. The projections described in this major report show the outstanding potential of the LHCb experiment to place very stringent constraints in the low-mass range, in complementarity with ATLAS and CMS.

3 Dark sector beyond QCD-like scenarios

3.1 SUEP

Contributors: Cari Cesarotti, Carlos Erice, Karri Folan DiPetrillo, Chad Freer, Luca Lavezzo, Christos Papageorgakis, Christoph Paus, Matt Strassler

Dark showers produced in hidden valley models need not result in collimated jets like in SM QCD. Events with spherically-symmetric, large multiplicities of low momentum charged particles, are also a possible phenomenology of strongly-coupled hidden valley models. This section discusses the motivation for these so called soft-unclustered-energy patterns (SUEPs), tools available for simulation, and typical phenomenology. Experimental challenges for SUEP searches at LHC general purpose detectors are also summarized.

3.1.1 Theoretical motivation

Although the production of quarks in a QCD-like confining theory leads inevitably to jets of hadrons, the details are not always the same. The width of the jets, the hadron multiplicity per jet, and the jet multiplicity all depend on the value of the running coupling. More specifically, jets arise from the partonic shower, which depends on the running ’t Hooft coupling $\lambda \equiv \alpha _\textrm{D} N_{c_\textrm{D}}$ ($\alpha _s N_c$ in QCD), evaluated at and somewhat below the energy scale at which the jets are produced.

In QCD-like theories, like in QCD, jets at lower energy, closer to the confinement scale and thus at larger $\lambda$, are broader, because gluon radiation at larger angles is more common with the large coupling. One might then imagine that if $\lambda$ could somehow be taken large without reducing the energy of the collisions, then the jets produced might become so broad and numerous as to blur together, creating a smooth distribution and a non-jetty final state. In QCD-like theories this question is almost academic, since $\lambda$ only becomes large near the confinement scale. In $e^+e^-$ collisions near the confinement scale, only a small handful of hadrons are produced, making it hard to define jets in the first place. Jets were only identified in $e^+e^-$ collisions at scales $\sim 10$ GeV, where $\lambda < 1$, and gluon jets only at 27 GeV.

However, in 1998 it was shown [105] that certain classes of supersymmetric conformal field theories (roughly speaking, these are theories with whose beta functions are all zero and are thus scale-invariant) are equivalent to string theories on certain curved spaces. These theories can have arbitrary values of $\alpha _\textrm{D}$ (now constant) and $N_{c_\textrm{D}}$. The duality can be used to compute many of properties of these theories when $\alpha _\textrm{D} N_{c_\textrm{D}} \gg 1\gg \alpha _\textrm{D}$. It was noted in [106] that rapid pdf evolution in this regime leads to an absence of hard partons inside hadrons. Since similar dynamics controls jet evolution, this naturally suggests that jets will be absent in this regime as well. Several groups [107,108,109] argued that there should be no jets in this limit, and it was proven convincingly in [109] that the correlation function of energy operators indicates that a partonic shower in this regime, allowed to proceed over a wide range of energies, will approach a spherically symmetric distribution on average. The distribution of momenta of these partons is not determined, however.

Note, importantly, that a spherically symmetric shower is not a consequence of conformal symmetry. A conformally invariant theory at small $\lambda$ will exhibit jets much like QCD itself; indeed, QCD at high energy is nearly conformally invariant, with scaling violated only by small logs (a fact which played an important role in the discovery of quarks.) Only at when $\lambda \gg 1 \gg \alpha _\textrm{D}$ are roughly spherical showers expected, and corrections to the spherical shape are of order $1/\sqrt{\lambda }$. Thus, for events that typically differ from spherical by $< 10\%$, one probably must have $N_{c_\textrm{D}} \gg 100$.

A conformally-invariant hidden sector is generally observable only as (except for small rare processes discussed in the unparticle literature) as the energy will be shared down to massless partons. Interesting hidden valley signatures arise only if the conformal invariance is broken at some scale $Q_c$ much lower than the production scale M. The shower that follows production of hidden partons is converted at $Q_c$ to a large number of hidden particles of small mass $m_\textrm{D} \lesssim Q_c$. If some of these are able to decay to the SM, then this can lead to a signature of many particles which are roughly spherically distributed in some frame of reference, not necessarily the lab frame [9, 107]. This signature is defined as a Soft Uncorrelated Energy Pattern, or SUEP.

We note that there have been disputes in the theoretical literature concerning whether the SUEP arising in this context is related directly to cascades of 5d KK states, with different points of view taken by [68, 110] versus [107, 111]. From the phenomenological point of view this may not matter very much in the near term, but the issue could arise if in future KK-cascade-based generators are created and assumed (perhaps incorrectly) to describe the same phenomenon as SUEP.

If the conformal symmetry breaking involves the complete Higgsing of a gauge group, it is clear that a near-spherical shower leads to a near-spherical SUEP. If the conformal symmetry breaking involves confinement, this is less clear and has not been proven; it may be that it is somewhat model-dependent. We will nevertheless assume it is the case, and that SUEPs can arise from the SM decay products of dark hadrons produced from a spherical shower of dark gluons. We will further assume that those dark hadrons can resemble those from a QCD-like spectrum, with low-lying spin-0 and spin-1 mesons and with limited roles for heavier mesons. Baryons should have mass of order $N_{c_\textrm{D}} \gg 1$ and will play no role.

Since techniques to calculate SUEPs, even in the few models where they are known to occur, do not yet exist, the simulations [9, 112] to study them are necessarily somewhat ad hoc. The assumption is made that particles are produced spherically according to a thermal distribution with an unknown temperature $T_\textrm{D}$ of order $Q_c$, the conformal breaking scale; see Sect. 3.1.2. (Note $m_\textrm{D}$, the typical hadron mass, may be of order $Q_c$ or may be much smaller, especially for pion-like states; in the latter case, multiplicities may be surprisingly low.) The value of $T_\textrm{D}/Q_c$ should be varied and treated as a quantifiable uncertainty.

Less clear is how to treat the uncertainties of the thermal approximation in the first place. A thermal model for hadronization in QCD works moderately well; perhaps this reflects the tendency for complex statistical systems to approach thermal ones. However, there is no proof that it would work for all or even most confining theories. One way to approach the uncertainty might be to vary the thermal spectrum in one or another way, using insights from studies of how systems equilibrate.

Finally, an uncertainty arises from the fact that in realistic theories the spherical approximation will suffer corrections, and at the current time those corrections have not been characterized theoretically nor incorporated into the simulation tools. Some effort to quantify this uncertainty ought to be undertaken.

Since SUEPs and the simulation packages used to simulate them represent a certain idealized situation, real signatures may differ from this idealized model. In this regard, the following should probably be kept in mind:

• Only a small number of theories of this class have been proven to exist, all of them supersymmetric and with $N_{f_\textrm{D}} \ll N_{c_\textrm{D}}$, where $N_{f_\textrm{D}}$ is the number of quarks in the fundamental representation.

• The statement that expected distributions are spherical receives substantial corrections, of order $1/\sqrt{\lambda }$ at the confinement scale. For an observably spherical SUEP, one may need $N_{c_\textrm{D}} \sim 100$. Because the particles of the hidden sector can appear in loops, coupling a mediator to a theory with $N_{c_\textrm{D}} \sim 100$ or more can have large consequences for the mediator or even the Standard Model sector; care must be taken in defining a consistent model.

• Event-to-event fluctuations can lead to large deviations from the SUEP idealization. When the multiplicity of visibly-decaying hadrons is low, Poisson fluctuations are large. For instance, if sixty dark hadrons are produced near-spherically but only ten decay visibly, the observed event will be far from spherical and potentially very asymmetric.

3.1.2 Simulation tools

In order to search for these novel signatures at high energy colliders, it is essential to generate events that capture the phenomenological characteristics of the energy pattern. Perturbation theory breaks down in the large coupling regime, and standard approaches to event simulation are unreliable. Novel methods are necessary to generate SUEP events.

Several simulation methods that can generate quasi-spherical energy patterns exist, such as black hole generators [113,114,115] or simplified 5d models [111]. However, these tools are not ideal for developing new analyses or triggering strategies. LHC experiments already have extensive search programs for black holes [116,117,118,119,120]. The latter tool does not have an obvious portal to connect to the SM. We will therefore discuss the utility of another tool in this section: the SUEP generator [112].

The simplified model used by this generator is described in Ref. [9]. In this framework, a hidden valley (HV) of new physics with confining dynamics is accessed via a heavy scalar mediator. A wide class of mediators can be used to connect the SM and HV, but a scalar portal from gluon fusion was chosen to both explore a triggering ‘nightmare scenario’ as well as study a potential rare Higgs decay mode. The scalar then decays into a high multiplicity of light HV mesons of a single flavor, of mass $m_\textrm{D}$, that follow a Maxwell–Boltzmann distribution

$$\begin{aligned} \frac{dN}{d^3{{{\textbf {p}}}}} \sim \exp \left( -\sqrt{{{\textbf {p}}}^2+m_\textrm{D} ^2}/T_\textrm{D} \right) . \end{aligned}$$

(10)

An illustration of this process is shown in Fig. 15. The user is free to select the temperature, $T_\textrm{D}$, and $m_\textrm{D}$, but for reasonable results the two parameters should satisfy $T_\textrm{D}/m_\textrm{D} > 1$ such that the final state is high multiplicity. Note that in a confining theory like QCD, the mass of the lightest meson can be much less than the confinement scale $\Lambda _\textrm{D} \sim T_\textrm{D}$, so $T_\textrm{D} \gg m_\textrm{D}$ is a motivated and physical choice. For a fixed scalar and meson mass, a higher temperature will correspond to fewer particles with a more significant boost. The tool is only intended to work in the high multiplicity limit, which means that the user must set $m_S\gg T_\textrm{D}, m_\textrm{D}$ for consistent results. A numerically small amount of momentum conservation violation may be observed in some events.

Fig. 15

Figure adapted from Ref. [9]

A visualization of a generic SUEP event. A heavy scalar S accesses a confining HV and showers into light mesons (e.g. $\pi _{\textrm{D}}$, here denoted $\phi$) before decaying to SM particles.

After production and showering, the dark mesons must decay back to SM for a visible signal. In this benchmark model, it is assumed that the dark mesons couple to a new U(1) gauge boson $A'_\gamma$ that kinetically mixes with SM hypercharge. The dark mesons decay into a $A'_{\gamma }A'_{\gamma }$ pair. Each $A'_\gamma$ then decays into SM particles, for example dilepton ($e^+e^-$ and $\mu ^+ \mu ^-$) or hadronized final states. The $A'_{\gamma }$ mass and branching ratio to Standard Model particles is configurable by the user.

This tool has the potential to inform and test diverse future analysis techniques. The phenomenology of the model depends on the choices of scalar mediator mass $m_S$, the dark meson mass $m_\textrm{D}$ and temperature $T_\textrm{D}$, as well as the decay branching ratios of the new particles. Depending on the signal of interest, the user can configure these values to achieve different multiplicity or species final states. Simple extensions of the package could include different possible mediators and decay portals.

3.1.3 Phenomenology

The classification of ‘SUEP’ is on the final state signature rather than a specific type of model. A SUEP (previously called ’soft bomb’ [9]) is usually an event with a high multiplicity of soft particles distributed quasi-isotropically in their rest frame. The underlying physics that produces such events can be varied. As discussed earlier, if the new particles interact strongly over a wide energy window, the shower develops by soft and isotropic emissions [36, 106, 109]. However, SUEPs can also develop from kinematics due to phase space arguments. If a new physics model includes many unstable particles with small mass splittings [110, 111], subsequent decays are unboosted and can approximate a spherical energy distribution for sufficiently large particle multiplicity. A common example of such event is R-parity violating SUSY [121, 122].

Since many different new physics scenarios can produce a SUEP signature, it is compelling to generically search for such events at colliders. As their common underlying feature is their global radiation pattern, event shape observables can serve as useful analysis tools. By studying the global event shape, it is both possible to quantify new physics and distinguish signal from background.

Event shape observables have been used to study QCD and measure $\alpha _s$ [123,124,125,126,127,128,129,130,131,132,133,134,135,136,137]. There have been several observables developed to quantify the degree to which a collider event is isotropic versus jet-like, including thrust [123, 138, 139], sphericity [140, 141], spherocity [142], and the C- and D-parameters [143,144,145]. While these observables have provided indispensable insights to QCD, they are most sensitive to deviations from dijet events rather than a robust probe of isotropy. To capture the phenomenology of a SUEP event it is necessary to define new observables which probe the opposite regime.

A new observable called event isotropy aims to study deviations from truly isotropic events [146]. This observable is defined using the Energy Mover’s Distance (EMD) [147, 148], the particle physics application of the Earth Mover’s Distance [149,150,151,152,153]. Event isotropy quantifies how ‘far’ an event is from isotropic, with smaller values indicating an event is more isotropic. Figure 16 shows event isotropy for SUEP benchmark models with different mediator masses and temperatures.

Fig. 16

Isotropy distributions are shown for different SUEP mediator masses on the left and for different temperatures on the right. The isotropy is calculated for ring geometry with segmentation 64 for generator level tracks. Tracks are defined as status $=1$ charged particles with $p_T > 0.1$ GeV and $|\eta | < 2.4$. Only events that pass the requirement $\sum p_T > 100$ GeV are kept. The production mode is Gluon Fusion (GF) with a dark photon branching fraction of $BR(A'_{\gamma } \rightarrow e\overline{e},\mu \overline{\mu },\pi \overline{\pi }) = (40,40,20)\%$

Figure 17 shows results from a related study that demonstrates how event isotropy is not strongly correlated with final state multiplicity, or reconstructed number of jets. Additionally, it correlates with canonical event shape observables much less than they correlate with each other in the quasi-isotropic regime. While there can be correlation with traditional event shapes, both types of observables can be used to better characterize the underlying physics.

Fig. 17

Figures adapted from [154]

Correlations between event isotropy and other canonical event shape observables calculated on quasi-isotropic events generated using the framework of Ref. [111], such as particle multiplicity (left top), jet multiplicity (right top), and the maximum eigenvalue of the sphericity tensor (left bottom), which is closely related to the C and D parameter in quasi-isotropic radiation patterns. For reference, we also include the correlation plot of the same samples in thrust and maximum eigenvalue (bottom right). The contours enclose 99% of the events for all of the samples.

3.1.4 Experimental aspects

The diffuse low momentum nature of SUEP events strongly resembles soft-QCD backgrounds at the LHC. Without high momentum final state particles, SUEP events with all-hadronic final states can easily be mistaken for pile-up collisions, and pose extreme challenges for general purpose detectors such as ATLAS and CMS.

Fig. 18

The $H_{T}$ distribution for different SUEP mediator masses is shown on the left. The $H_{T}$ is calculated by taking the scalar sum of $p_{T}$ for generator level jets with $p_{T}>20$ GeV and $|\eta |>4.7$. The number of tracks per event for different SUEP mediators is shown on the right. Tracks are defined as status $=1$ charged particles with $p_T > 0.7$ GeV and $|\eta | < 2.5$. Both plots show Gluon Fusion (GF) production mechanism with a dark photon branching fraction of $BR(A'_{\gamma } \rightarrow u\overline{u})=100\%$

The first and most challenging step for any SUEP analysis is to identify a trigger strategy. The trigger systems of ATLAS and CMS operate in a two-step process to determine which events are saved for analysis. The first stage, Level 1, makes a fast decision incorporating coarse calorimeter and muon information. The second stage, the High Level Trigger (HLT), makes use of refined calorimeter and muon information and adds limited tracking. SUEP events typically have low efficiency for traditional triggers, which are designed to reject pile-up. However, there remain several possible analysis strategies utilizing data already collected during Run 2 of the LHC.

In order to characterize the difficulty of observing lower mass mediators, several benchmarks points are used: mediator masses of 1000 GeV, 750 GeV, 400 GeV, and 125 GeV. Here, the choice of the 125 GeV benchmark is motivated by the observed Higgs boson, which may itself serve as a portal to the hidden sector. The mediator is assumed to be produced via gluon fusion (GF). For the Higgs portal benchmark, associated production (ZH) also considered, where the Z boson decays leptonically. All signals are produced assuming dark mesons have a mass of 2 GeV, and the system has a temperature of 2 GeV. These dark mesons subsequently decay into a pair of dark photons. Multiple dark photon decays are explored for GF production using dark photon branching ratios; $BR(A'_{\gamma } \rightarrow u\overline{u})=100\%$, with $m(A'_{\gamma })=1$ GeV, $BR(A'_{\gamma } \rightarrow e\overline{e},\mu \overline{\mu },\pi \overline{\pi }) = (40,40,20)\%$ with $m(A'_{\gamma })= 0.5$ GeV labelled as leptonic, and $BR(A'_{\gamma } \rightarrow e\overline{e},\mu \overline{\mu },\pi \overline{\pi }) = (15,15,70)\%$ with $m(A'_{\gamma })= 0.7$ GeV labelled as hadronic.

Trigger strategies based on the scalar sum of hadronic activity ($H_{T}$) in the event can be used to target mediators produced via GF. In order to prevent extremely high background rates from QCD, ATLAS and CMS typically required $H_{T} > 500$ GeV at Level 1, and $H_{T} > 1$ TeV at HLT throughout Run 2. The vast majority of signal events which pass these requirements involve a mediator recoiling against an initial state radiation jet with high $p_{T}$. Trigger efficiency significantly decreases as the mediator mass decreases and the total energy deposition decreases. Figure 18 shows the $H_{T}$ distribution as well as the number of tracks for different SUEP mediator masses.

pending on the production mechanism and branching fraction of the dark photon, SUEP events may contain multiple final state muons. In these cases, muon triggers designed to target vector boson or b-physics processes can provide much higher trigger efficiency than $H_{T}$ triggers, especially for lower mass mediators. Figure 19 shows the number of muons and the leading muon $p_{T}$ for several Higgs portal SUEP scenarios. Signals with larger leptonic branching ratios result in large multiplicities of moderate momentum muons, and can be targeted with tri-muon triggers. Leptons produced via associated production can be targeted with single and double-muon triggers.

Fig. 19

The number of muons at the generator level for different SUEP production mechanisms and dark photon decay branching fractions is shown on the left. Muons are selected with $p_T > 3$ GeV and $|\eta | < 2.5$. The muon with the highest $p_T$ is shown on the right for different SUEP production mechanisms and dark photon decay branching fractions

The upcoming Run 3 at the LHC offers an opportunity to design new triggers that specifically target SUEP signatures. One possibility is to use a standard $H_{T}$ or multi-jet trigger at Level 1, with a large multiplicity of tracks in the High Level Trigger. SUEP events are likely to be boosted after passing the Level 1 trigger, with SUEP decay products recoiling against Standard Model jets from initial state radiation. Figure 20 shows how events become less isotropic with increasing $\sum p_T$ requirements. As a result, using event shape observables would be suboptimal at HLT. In contrast, it is possible to choose an HLT track multiplicity requirement that is highly efficient for SUEP events, and also reduces backgrounds such that the $H_{T}$ requirement can be kept as low as Level 1 thresholds.

Fig. 20

Isotropy distributions are shown for different values of the $\sum p_T$ event selection. The $\sum p_T$ is calculated for generator level tracks that are defined as status $=1$ charged particles with $p_T > 0.1$ GeV and $|\eta | < 2.4$. The production mode is Gluon Fusion (GF) with a dark photon branching fraction of $BR(A'_{\gamma } \rightarrow e\overline{e},\mu \overline{\mu },\pi \overline{\pi }) = (40,40,20)\%$

A trigger which requires $H_T>500$ GeV and at least 150 tracks per event would yield a QCD efficiency one order of magnitude smaller than the currently employed $H_T>1050$ GeV trigger, while recovering nearly all the SUEP events that already pass the $H_T>500$ GeV selection. Note that this study assumes tracks can be reconstructed with nearly full eff

iciency and associated to the primary vertex.

If track reconstruction is too computationally expensive to run in the High Level trigger or has sub-optimal performance, it is possible to design a trigger which counts the number of hits in the innermost layers of the ATLAS and CMS tracking detectors. This approach was studied in Ref. [9]. In addition to being less computationally expensive, hit counting is sensitive to even the softest particles in SUEP events. Charged particles with $p_{T}>{\mathcal {O}}(10)$ MeV can reach the innermost layer of the tracker and produce a hit. One trade-off is that hits cannot be associated to a primary vertex, and pile-up collisions increase the number of hits per event for background. Hits from SUEP events are likely to be localized in z near the primary vertex, and this can be used to discriminate against pile-up.

There are also several potential strategies to trigger on SUEP events directly at Level 1. When SUEP events are not boosted, final state particles create a band which spans all $\phi$ and are centered around a definite value of $\eta$. It may be possible to reduce thresholds by designing a trigger which looks for high $H_T$ in one $\eta$-slice compared to the rest of the event.

Additional Level 1 strategies depend on the final state particles in the SUEP event. An all-electron or all-photon signal would create an abnormally high ratio of energy deposited in the Electromagnetic Calorimeter compared to the Hadronic Calorimeter. This method could be further refined by requiring Electromagnetic energy be centered around a single $\eta$-slice. Alternatively, SUEP events with dark matter particles in the final state, semi-visible SUEP could potentially be accessed via a missing-transverse-momentum trigger.

Both ATLAS and CMS will upgrade their detectors and trigger schemes at the High Luminosity LHC, offering further opportunities to design new SUEP triggers. The projected CMS trigger at the HL-LHC will reconstruct charged particles with $p_{T}>2$ GeV at Level 1. This new scheme would enable even more optimization of the trigger design at the L1 level. One could optimize a selection on the number of tracks at L1 to recover an even greater quantity of low $H_T$ events. This strategy would be most effective for higher temperature scenarios.

Once a trigger strategy has been determined, it is essential to reconstruct the low-momentum charged particles associated with SUEP signatures. Standard ATLAS and CMS track reconstruction is highly efficient for charged particles which traverse roughly 8 layers of the silicon trackers [155, 156]. In CMS reconstruction is roughly $90\%$ efficient for charged particles with $p_{T} \ge 1$ GeV, and roughly $60\%$ efficient for particles with $p_{T}\sim 300$ MeV. In ATLAS, nominal track reconstruction has a minimum requirement of $p_{T} \ge 500$ MeV. Standard reconstructed tracks typically have impact parameter resolutions on the order of $100~\upmu {\textrm{m}}$ [157], which enables tracks to be associated to the primary vertex (proton-proton collision of interest). This primary vertex requirement ensures that the computed track multiplicity is not biased by nearby pile-up collisions.

There are additional possibilities to reconstruct charged particles with even lower momenta. In ATLAS, tracks with $p_{T}>100$ MeV can be reconstructed for a small subset of events [158]. These events would be processed with an additional pass of tracking, using leftover hits, in a region of interest of $z \sim 1$ mm around the primary vertex. To access even lower momentum charged particles, $p_{T} \sim 10$ MeV, it is possible to count the number of hits in the inner most layers of the detector. While this strategy would increase the acceptance for extremely low-momentum SUEP particles, it is impossible to associate hits to a particular primary vertex. This strategy also requires accessing low level data which may not be possible due to storage constraints.

For SUEP signatures with higher temperatures, a significant number of final state Standard Model particles will have moderate momenta, $p_{T}>3$ GeV. For these scenarios, it is possible to identify final state particles as muons, electrons, or hadrons. Events with a large multiplicity of low momentum leptons would be a particularly striking signal, and this information could be used to further improve signal to background discrimination.

After reconstruction, the unique properties of SUEP events can then be used to separate against signal from QCD background. The two most powerful observables include the characteristic high track multiplicity, and the isotropic distribution of such tracks.

For analyses using data that were collected during Run 2, trigger strategies will likely rely on the presence of additional objects produced in association with the SUEP shower, either based on a high quantity of ISR QCD or a massive gauge boson decaying to triggering objects. In both cases, a simple procedure to recover the SUEP shower can be followed. First, the triggering object can be identified with standard reconstruction techniques and subtracted from the overall event representation. Second, the remaining tracks can be boosted against the triggering object. In case the later is produced back-to-back with the mediator, these would allow to recover both the track multiplicity and spherically symmetrical distribution of the SUEP shower.

Recently, the anomaly detection techniques as a generic way to search for new physics have been incorporated to SUEP searches. A proposal based on autoencoders [159] shows the strength of such techniques for low mediator masses and temperatures. Additional details on this approach are presented in Sect. 5.4.

3.2 Glueballs

Contributors: David Curtin, Caleb Gemmell, Christopher B. Verhaaren

In the $N_{f_\textrm{D}} = 0$ limit for $SU(N_{c_\textrm{D}})$ Yang–Mills theories, the only hadronic states that form below the confinement scale are glueballs, composite gluon states. This limit is unique because there are no light degrees of freedom below the confinement scale. In the absence of such light states color flux tubes cannot break via the creation of quark-antiquark pairs. This process is essential to present QCD hadronization models [160, 161], thus the usual understanding of hadronization does not directly apply to the $N_{f_\textrm{D}} = 0$ limit. This qualitative difference has hindered efforts to study dark glueball showers in scenarios with pure Yang–Mills dynamics.

Sectors with $N_{f_\textrm{D}} = 0$ $SU(N_{c_\textrm{D}})$ Yang–Mills descriptions commonly appear in neutral naturalness theories, such as Twin Higgs models [11, 12], Folded Supersymmetry [162] and many others [163,164,165,166]. The glueballs of a hidden confining sector have also been considered as dark matter candidates [23, 26, 167, 168]. Thus, $N_{f_\textrm{D}} = 0$ $SU(N_{c_\textrm{D}})$ Yang–Mills models are motivated by possible solutions to the Little Hierarchy problem and the unknown nature of dark matter. Current ignorance of the pure-glue hadronization process has left these models in an largely unstudied corner of motivated parameter space. The hadronization process determines the final state multiplicity and energy distribution of dark glueballs, which are essential for both collider and indirect detection studies.

The properties of the glueballs themselves are relatively well-known, having been studied in lattice gauge theory [169,170,171,172,173]. In the absence of external couplings, these studies have established a spectrum of 12 stable glueball states characterized by their $J^{PC}$ quantum numbers. When considered as part of a dark sector, these glueballs can be stable or decay through a variety of portals to the SM [43, 45], with possibly long lifetimes on collider or cosmological time scales. This spectrum can entirely be parameterised by the confinement scale of the theory, or equivalently lightest glueball mass, $m_0 \sim 6\Lambda _\textrm{D}$. These glueball properties are also known for several $N_{c_\textrm{D}} \ne 3$ which paves the way for studying exotic dark sectors outside the standard dark SU(3) case [23, 167, 174, 175].

Recently, efforts have been made to enable quantitative studies of pure Yang–Mills parton showers and hadronization [75]. This includes the creation of a new public python package, GlueShower.² This package allows users to simulate dark glueball showers produced from an initial pair of dark gluons. GlueShower combines a perturbative pure-glue parton shower with a self-consistent and physically motivated parameterization of our ignorance regarding the unknown glueball hadronization behaviour. Two qualitatively different hadronization possibilities are included: a more physically motivated jet-like assumption, and a more exotic plasma-like option that accounts for the possibility that color-singlet gluon-plasma-states are created by hypothetical non-perturbative effects far above the confinement scale. Each such plasma-ball then decays isotropically to glueballs in its restframe, somewhat akin to dark hadron production in SUEP scenarios [9, 107, 159]. Within each hadronization option, two nuisance parameters control the hadronization scale and the hadronization temperature, mostly controlling the glueball multiplicity and relative abundance of different glueball species respectively.

The study [75] defines a set of 4 benchmark points for these nuisance parameters in each option to represent the range of physically reasonable glueball hadronization possibilities. For phenomenological studies, the range of predictions spanned by these hadronization benchmarks can be interpreted as a theoretical uncertainty on predictions for glueball production. Despite the wide range of possibilities for hadronization that are considered, most glueball observables are predicted within an ${\mathcal {O}}(1)$ factor. This is in large part due to the modest hierarchy between the glueball mass and the confinement scale, which makes most inclusive glueball observables dominantly dependant on the physics of the perturbative gluon shower. However, exclusive production of certain glueball states can vary by up to a factor of 10 depending on hadronization assumptions, a range that accurately represents our current theoretical uncertainty.

In summary, to support the increasing interest in dark showers, efforts are being made to ensure the possibility space is covered. Until recently the $N_{f_\textrm{D}} = 0$ limit has been largely ignored due to difficulties in studying pure-glue hadronization. However, such hidden sectors are motivated by both naturalness concerns and as a possible dark matter particle. The GlueShower tool aims to effectively facilitate studying the $N_{f_\textrm{D}} = 0$ limit, even without a full understanding of the underlying non-perturbative hadronization physics. Despite the current theoretical limitations, this demonstrates that quantitative studies of and searches for glueball signatures can be reliably conducted if the underlying uncertainties are accurately accounted for.

4 Simulation tool limitations and how to build consistent benchmarks from the underlying physical parameters for semi-visible jets

4.1 Consistent parameter setting and roadmap for improving on the simulation of dark showers

Contributors: Suchita Kulkarni, Seán Mee, Matt Strassler

The non-perturbative nature of QCD-like strongly interacting scenarios makes it impossible to set consistent UV and IR parameters based purely on perturbative analysis. In this section, we address this problem, going beyond existing efforts in the literature. First, we sketch the importance of lattice calculations to set low energy bound state masses given UV parameters. Second, we illustrate the importance of portal phenomenology and associated symmetry breaking patterns, and use chiral Lagrangian techniques to set the interactions of the low energy bound states among themselves and to the SM final states. Finally, we comment on the hadronization parameters necessary for LHC phenomenology and make some observations for a subset of them. After this, we turn to the simulation of benchmark models consistent with the above observations. We describe the recent improvements to the pythia 8 Hidden Valley module, and present a few benchmarks which are used in later sections for studies.

4.1.1 UV scenarios: SM extension with non-Abelian gauge groups

We suppose the Standard Model is extended with a new sector, consisting of an additional non-Abelian gauge group $SU(N_{c_\textrm{D}})$ with $N_{f_\textrm{D}}$ degenerate Dirac fermions $\textrm{q}_{\textrm{D} \alpha }$ in the fundamental representation, with current mass $m_{\textrm{q}_\textrm{D}}$. We will refer to the new sector as the “dark sector,” though we note this is fully equivalent to a “hidden valley” as defined in [36]. This sector has $SU_L(N_{f_\textrm{D}}) \times SU_R(N_{f_\textrm{D}})\times U(1)_B$ global symmetry broken by the mass term to a diagonal $SU(N_{f_\textrm{D}})\times U(1)$.³ We will assume $N_{c_\textrm{D}},N_{f_\textrm{D}}$ are such that the theory confines and has a chiral $\textrm{q}_{\textrm{D} \alpha } \overline{\textrm{q}}_{\textrm{D} \alpha }$ condensate, which in the absence of $m_{\textrm{q}_\textrm{D}}$ would spontaneously break the chiral symmetries and lead to Nambu–Goldstone bosons. Instead, as in QCD itself, we have pseudo-Nambu–Goldstone bosons with masses that are proportional to $\sqrt{m_{\textrm{q}_\textrm{D}}}$.

If the only connection between this sector and the SM is through a mediator (or “portal”) which is either massive or ultra-weakly coupled, then typical confining Hidden-Valley-type phenomenology inevitably results. Specifically, production of the “dark quarks” leads to production of dark hadrons. These will be collimated in jets if the theory is QCD-like and the invariant mass of the produced $\textrm{q}_\textrm{D} \overline{\textrm{q}}_\textrm{D}$ pair is far above the dark confinement scale. We will specifically consider a mediator in the form of a heavy $U(1)^\prime$ leptophobic $\textrm{Z}^{\prime }$ mediator between SM quarks and the dark quarks. These are the same models introduced in Sect. 2.1.1, used widely in the semi-visible jet searches and are also considered in the dark matter working group [176]. The process $\textrm{q} \overline{\textrm{q}} \rightarrow \textrm{Z}^{\prime } \rightarrow \textrm{q}_\textrm{D} \overline{\textrm{q}}_\textrm{D}$ allows dark hadrons to be produced, while the $Z^\prime$ also allows some dark hadrons to decay to SM quark-antiquark pairs, leading to an all-hadronic signal. We will assume throughout that $m_{\textrm{Z}^{\prime }}$ is much larger than the confinement scale, by at least $\sim 30$, so that the physics in the dark sector actually leads to jets of dark hadrons.⁴ Because a fraction of the dark hadrons are typically stable, at least on LHC-detector time scales, the observable signal usually consists of at least two relatively fat jets with considerable substructure, and often a high multiplicity of SM hadrons, along with roughly collinear . These “semi-visible jets” (SVJ), introduced in Sect. 2.1.1, are the target of the searches in question here.

The mediator’s couplings to the dark sector will break the $SU(N_{f_\textrm{D}})\times U(1)$ flavor symmetry to a smaller subgroup $G_f$. Without this breaking, the majority of the dark hadrons would be charged in the adjoint of $SU(N_{f_\textrm{D}})$ and would be unable to decay to a SM final state, which would be a singlet under this symmetry. (Note that decays with both dark and SM particles in the final state would still be permitted.) If the couplings are vector-like, we assign the dark quarks charges $Q_{\alpha }$ under the $U(1)^\prime$, and define the charge matrix $\textbf{Q}$ as a diagonal matrix with eigenvalues $Q_{\alpha }$; the group $G_f$ is the subgroup of $SU(N_{f_\textrm{D}})\times U(1)$ which commutes with $\textbf{Q}$. (In chiral models the left- and right-handed quarks have different charges $Q_{\alpha }$ and $\tilde{Q}_{\alpha }$, and get their masses from a Higgs field with charge $Q_{\alpha }-\tilde{Q}_{\alpha }$. We will not consider such models in detail here.) The precise choice of $\textbf{Q}$ has a significant impact on the phenomenology.

The ultraviolet Lagrangian for the hidden sector is

(11)

where $G_{D,\mu \nu }$ is the gauge field strength tensor and $m_{\textrm{q}_\textrm{D}}$ is the current mass of the dark quarks. This part of the theory has two discrete parameters $N_{c_\textrm{D}}$ and $N_{f_\textrm{D}}$, and two continuous parameters, the running gauge coupling $\alpha _\textrm{D} (\mu )$, with $\mu$ a renormalization scale, and the “current” quark mass $m_{\textrm{q}_\textrm{D}}$. Since neither of these parameters has direct contact with the observable phenomena, we replace them with the confinement scale $\Lambda _\textrm{D}$, or some proxy for it, such as the one-loop dimensional transmutation scale, and the mass ${m_{\pi _{\textrm{D}}}}$ of the light pseudoscalar mesons $\pi _{\textrm{D}}$.

The interaction of this sector with the SM via a $\textrm{Z}^{\prime }$ mediator takes the form

$$\begin{aligned} {\mathcal {L}}_{\textrm{int}}\subset -e_D \textrm{Z}^{\prime } _{\mu } \sum _{\alpha }\overline{\textrm{q}}_{\textrm{D} \alpha } Q_{\alpha } \gamma ^{\mu } \textrm{q}_{\textrm{D} \alpha }- g_{\textrm{q}} \,\textrm{Z}^{\prime } _{\mu } \sum _{i} \overline{\textrm{q}} _i\gamma ^{\mu }\textrm{q} _i \end{aligned}$$

(12)

if the $\textrm{Z}^{\prime }$ couplings are vectorlike and the $\textrm{q}_{\textrm{D} \alpha }$ have charge $Q_{\alpha }$. (If the couplings to the hidden sector are chiral, separate charge $Q_{\alpha }$ and $\tilde{Q}_{\alpha }$ must be assigned for left- and right-handed hidden quarks.) We will assume all SM quarks have the same charge under $U(1)^\prime$ for simplicity, though in realistic models one must account for the differences, which can affect observables, such as the SM heavy flavor fractions in the SVJs.

One other issue of importance is whether the $\textrm{Z}^{\prime }$ and the $\textrm{q}_\textrm{D}$ obtain their masses from the same source, in such a way that the longitudinal polarization of the $\textrm{Z}^{\prime }$ mixes with the $\pi _{\textrm{D}}$ from chiral symmetry breaking. Similar mixing of the SM charged pion with the $\textrm{W}$ bosons allows the classic decay $\pi \rightarrow \mu \nu$. By analogy, the $\textrm{Z}^{\prime }$ mixing with the $\pi _{\textrm{D}}$ affects the decays of the latter to the SM.

4.1.2 From ultra-violet theories to infrared parameters

The $SU(N_{c_\textrm{D}})$ confines at around the scale $\Lambda _\textrm{D}$ and various dark hadronic bound states are produced. In the exact $SU(N_{f_\textrm{D}})$ limit, the lightest hadrons consist of the spin-0 flavor-adjoint $\pi _{\textrm{D}}$, which are pseudo-Nambu–Goldstone bosons (PNGBs), the spin-1 flavor-adjoint $\rho _\textrm{D}$, the spin-1 flavor-singlet $\omega _\textrm{D}$, and the spin-0 flavor-singlet $\eta '_\textrm{D}$. In general ${m_{\pi _{\textrm{D}}}}<m_{\rho _\textrm{D}} \lessapprox m_{\omega _\textrm{D}}$, while the $\eta '_\textrm{D}$ mass depends on the anomaly, which scales like $\sqrt{N_{f_\textrm{D}}/N_{c_\textrm{D}}}$ when it is dominant. Thus for $N_{f_\textrm{D}}$ flavors, the theory contains $N_{f_\textrm{D}} ^2 -1$ mass-degenerate pions and an equal number of degenerate rho mesons, along with an omega which will be slightly heavier than the rhos, and an eta-prime which may be near-degenerate with the pions for $N_{f_\textrm{D}} \ll N_{c_\textrm{D}}$ but is much heavier for $N_{f_\textrm{D}} \sim N_{c_\textrm{D}}$. In addition there are baryons and antibaryons with mass of order $N_{c_\textrm{D}} \Lambda _\textrm{D}$ (bosons for $N_{c_\textrm{D}}$ even, fermions otherwise), except for $N_{c_\textrm{D}} =2$ in which case they are exactly degenerate with the pions and are themselves PNGBs. We note from these remarks that $N_{f_\textrm{D}} =1$ and $N_{c_\textrm{D}} =2$ are special cases, which must be treated with care.

There are ambiguities in defining the scale $\Lambda _\textrm{D}$. One way to define it is via the running gauge coupling $\alpha _\textrm{D} (\mu )$ at one loop. As pointed out by ’t Hooft, physics in the large $N_{c_\textrm{D}}$ limit depends mainly on $\alpha _\textrm{D} (\mu )N_{c_\textrm{D}}$, up to $1/N_{c_\textrm{D}}$ corrections. Lattice results show that $N_{c_\textrm{D}} =3$ is already close to $N_{c_\textrm{D}} =\infty$, and moreover the physics of a QCD-like shower has very small $1/N_{c_\textrm{D}}$ corrections, a fact that the pythia showering routines take advantage of. With this in mind, the one-loop running coupling can be written in a form familiar from QCD:

$$\begin{aligned} \alpha _\textrm{D} (\mu ^2) N_{c_\textrm{D}} = \left[ {\frac{1}{2\pi }\left( \frac{11}{3}-\frac{2}{3} \,\frac{N_{f_\textrm{D}}}{N_{c_\textrm{D}}}\right) \log \left( \frac{\mu }{\Lambda _\textrm{D}}\right) }\right] ^{-1}. \end{aligned}$$

(13)

This form emphasizes that the physics of this theory is really a function of $N_{f_\textrm{D}}/N_{c_\textrm{D}}$, with $1/N_{c_\textrm{D}}$ corrections, for large $N_{c_\textrm{D}}$. To this end in Fig. 13(left), we show the running of $\alpha _D$ at one lop for several values of $N_{c_\textrm{D}}/N_{f_\textrm{D}}$ for a fixed $\Lambda _D = 1$ GeV.

Although the one-loop formula for $\Lambda _\textrm{D}$ is currently used in the pythia 8 HV module and by us throughout this document, this situation should be viewed as temporary. Inevitably, this method indexes non-perturbative hadronic masses to a scale which is perturbative, and this may pose challenges for interpreting results from lattice gauge theory in which $\Lambda _\textrm{D}$ is often defined non-perturbatively, for instance through the string tension. Moreover, the connection between this one-loop estimate of $\Lambda _\textrm{D}$ and the physical confinement scale becomes less and less accurate as $N_{f_\textrm{D}}/N_{c_\textrm{D}}$ increases; two-loop effects become important, with the effect that non-perturbative definitions of $\Lambda _\textrm{D}$ will be much smaller than the one-loop definition. It seems likely a two-loop perturbative definition would be significantly closer to non-perturbative ones. To illustrate for this effect in Fig. 13(right) we demonstrate the effect of two loop correction for running of $\alpha _D$ (solid lines) in comparison with the corresponding one loop correction (dashed lines), for two different values of $N_{f_\textrm{D}} = 3,6$ with fixed value of $N_{c_\textrm{D}} = 3$ keeping $\Lambda _D = 1$ GeV. In order to derive this result, we have used the procedure as described in [177] with beta function as defined in [178].

In fact for sufficiently large $N_{f_\textrm{D}}/N_{c_\textrm{D}}$ the theory will no longer confine because its running coupling reaches an infrared fixed point. The values for $N_{c_\textrm{D}},N_{f_\textrm{D}}$ where this occurs are not precisely known. For various $N_{c_\textrm{D}}$, estimates of the value of $N_{f_\textrm{D}}$ above which the theory is believed not to confine are computed in e.g. [179] and tabulated in Table 2. These results suggest that theories with $N_{f_\textrm{D}} < 3N_{c_\textrm{D}}$ likely confine, and we will refer to them as “QCD-like” theories.

Again, we have at this time only used the one-loop formula above to define $\Lambda _\textrm{D}$, and we use the same definition for the parameter Lambda in pythia 8, which is perhaps reasonable for $N_{f_\textrm{D}}/N_{c_\textrm{D}} \sim 1$ or below. But it is important to note that for $N_{f_\textrm{D}}/N_{c_\textrm{D}} \rightarrow 3$, the one-loop running is inaccurate and at a minimum a two-loop formula (within which the fixed points at large $N_{c_\textrm{D}}, N_{f_\textrm{D}}$ can be observed) ought to be used (Fig. 21).

Table 2 Values of $N_{c_\textrm{D}}, N_{f_\textrm{D}}$ which lead to asymptotically free $SU(N_{c_\textrm{D}})$ theories. It is important to stay well below these values in order to remain within the QCD-like scenarios of current interest to us

Fig. 21

Left: The ’t Hooft coupling $N_{c_\textrm{D}} \,\alpha _\textrm{D} (\mu )$, computed at one loop, as a function of $\mu /\Lambda _\textrm{D}$ for several values of $N_{f_\textrm{D}}/N_{c_\textrm{D}}$. Right: Two loop running of $\alpha _\textrm{D} (\mu )$ for $N_{c_\textrm{D}} = 3, N_{f_\textrm{D}} = 3, \,6 \text { and } \Lambda _\textrm{D} = 1\,\text {GeV}$

The masses and couplings of the low-lying bound states are a direct consequence of UV parameters, but are not calculable analytically. Rough estimates of these quantities can be obtained by combining lattice gauge theory calculations with the chiral Lagrangian for spin-0 mesons, extended by including spin-one mesons as though they were flavor-symmetry gauge bosons. It is convenient to set the overall scale of the dark hadrons using a non-perturbative definition of the confinement scale, which we will call $\tilde{\Lambda }_\textrm{D}$ and specify later, and to express all other dimensionful dark-hadronic quantities in terms of this parameter. Once we have fixed $\tilde{\Lambda }_\textrm{D}$ as the overall hadronic scale, and specified the ratio

$$\begin{aligned} \frac{{m_{\pi _{\textrm{D}}}}}{\tilde{\Lambda }_\textrm{D}}\propto \sqrt{\frac{m_{\textrm{q}_\textrm{D}}}{\tilde{\Lambda }_\textrm{D}}} \end{aligned}$$

(where the relation ${m_{\pi _{\textrm{D}}}}\sim \sqrt{m_{\textrm{q}_\textrm{D}}}$ follows from the chiral Lagrangian), everything else should be computable in principle.

Such lattice calculations for mass degenerate fermions in the fundamental representation of SU(N) gauge theories are available in abundance, albeit in the quenched approximation, see e.g. [180,181,182,183]. A particularly useful resource is [184], which summarises the spectrum of mesons in the large-N limit of QCD-like theories. These calculations can be used to determine the ratios of the dark hadron masses as a function of the hidden sector parameters.

Using lattice calculations and fits plotted in Fig. 19 of [185], we can relate the dark quark mass to the dark pion mass and the dark rho mass. We express these in terms of $\tilde{\Lambda }_\textrm{D}$ defined as the chiral limit ($m_\pi \rightarrow 0$) of the $\rho$ mass divided by 2.37. (In terms of the physical units used in Fig. 19 of [185], this puts the analogue of $\tilde{\Lambda }_\textrm{D}$ for physical QCD at 300 MeV.) These relations are concretely shown in Fig. 22. Our analytic fits to the curves shown are

$$\begin{aligned} \frac{{m_{\pi _{\textrm{D}}}}}{\tilde{\Lambda }_\textrm{D}} = 5.5\,\sqrt{ \frac{m_{\textrm{q}_\textrm{D}}}{\tilde{\Lambda }_\textrm{D}}} \quad \frac{m_{\rho _\textrm{D}}}{\tilde{\Lambda }_\textrm{D}} = \sqrt{5.76 + 1.5 \frac{{m_{\pi _{\textrm{D}}}}^2}{\tilde{\Lambda }_\textrm{D} ^2}}. \end{aligned}$$

(14)

The fit functions and coefficients shown in Eq. 14 are appropriate for small ${m_{\pi _{\textrm{D}}}}/\tilde{\Lambda }_\textrm{D}$, though they work far beyond this expectation, and begin to differ from lattice computations by $>10\%$ only for ${m_{\pi _{\textrm{D}}}}/\tilde{\Lambda }_\textrm{D} > 2.3$.

Fig. 22

Fits given in Eq. 14 for the $\rho _\textrm{D}$ mass (left) and $\pi _{\textrm{D}}$ mass (right) to results from lattice simulations [185]. The left panel also indicates the kinematic thresholds for $\rho _\textrm{D}$ to decay to $\pi _{\textrm{D}} \pi _{\textrm{D}}$ and $\pi _{\textrm{D}} \pi _{\textrm{D}} \pi _{\textrm{D}}$

The relation between the perturbative $\Lambda _\textrm{D}$ of Eq. 13 (or its higher-loop version) and the non-perturbative $\tilde{\Lambda }_\textrm{D}$, defined in terms of the chiral limit of $m_{\rho _\textrm{D}}$ (or some other similar definition), is not established. Although they are proportional, the proportionality depends on $N_{c_\textrm{D}}$ and $N_{f_\textrm{D}}$ (mostly on $N_{f_\textrm{D}}/N_{c_\textrm{D}}$.) In what follows, we will assume they are the same. The uncertainties and inaccuracies that result from this choice can only be reduced in future through more careful matching between perturbative and non-peturbative quantities.

The spin-1 singlet is expected to be nearly degenerate with the spin-1 adjoint hadrons, so there is little importance in giving it a different mass.⁵ By contrast, the spin-0 singlet is expected to be heavier than the spin-0 adjoint, possibly by a large amount if $N_{f_\textrm{D}} \sim N_{c_\textrm{D}}$, due to the axial anomaly. Some analysis of QCD hadrons using the chiral Lagrangian, which we will present elsewhere, suggests

$$\begin{aligned} m_{\eta '_\textrm{D}} ^2 \approx {m_{\pi _{\textrm{D}}}}^2 + \frac{N_{f_\textrm{D}}}{N_{c_\textrm{D}}}(3\tilde{\Lambda }_\textrm{D})^2. \end{aligned}$$

(15)

Thus for $N_{f_\textrm{D}} \sim N_{c_\textrm{D}}$, as in SM QCD, the splitting of the spin-0 singlet from the adjoint is large. Note that the factor of 3 in front of $\tilde{\Lambda }_\textrm{D}$ depends on the precise definition of the corresponding $\tilde{\Lambda }_{QCD}$, and will retain some uncertainties until this definition is handled more carefully. In any case, our current benchmark models presented below do not account for the anomaly term, and instead treat the $\eta '_\textrm{D}$ as degenerate with the $\pi _{\textrm{D}}$. However, the new version of pythia 8 includes a parameter HiddenValley:separateFlav, described in Sect. 4.1.4 below; when it is set “on,” the $\eta '_\textrm{D}$ mass can be set separately from that of the $\pi _{\textrm{D}}$ states.

For the scenarios we are interested in, a hard process leads to production of dark quarks which shower over a wide energy range; the shower then subsequently hadronizes. Hadronization in the dark sector is a far more challenging problem, because neither lattice calculations nor effective field theory methods are applicable to this process. All we know of it arises from studies of QCD data, and in particular through the use of phenomenological models (such as the Lund string model used in pythia, or the clustering model used in Herwig) whose parameters are tuned to fit experimental results. Currently there is no theoretical insight into how these models or their parameters should be adjusted for different values of $N_{c_\textrm{D}}$, $N_{f_\textrm{D}}$, or ${m_{\pi _{\textrm{D}}}}/\Lambda _\textrm{D}$. Consequently we take existing parameterizations from data as a starting point. One must vary these parameters within reason to obtain a sense of uncertainties.

In the present context, we are using pythia 8’s HV module, whose four main parameters HiddenValley:aLund, HiddenValley:bmqv2, HiddenValley:rFactqv, HiddenValley:sigmamqv parallel those of the corresponding QCD hadronization routine. The dimensionless HV parameters are set to exactly the same values as the dimensionless QCD parameters, while those with dimensions are scaled by the ratio of constituent quark mass parameter mqv = 4900101:m0 (which is not the current quark mass $m_{\textrm{q}_\textrm{D}}$ but rather a phenomenological parameter, of order $\Lambda _\textrm{D}$ for small $m_{\textrm{q}_\textrm{D}}$) to the constituent quark mass in QCD, which for u, d quarks is 330 MeV. There are other parameter tunes proposed by pythia 8 experts, see for example the Monash tune [186]. It is probably wise to try two or more tunes that are known to work in QCD as a means of estimating a minimum systematic error from this source. However, we have not studied this, and so further investigation is needed before an informed recommendation can be made.

There are three other hadronization parameters that are currently in use in the HV module. About the parameter HiddenValley:probVector, which gives the probability that a new meson formed in the hadronization should be assigned to spin-1 rather than spin-0, we have two pieces of information. Were spin-0 and spin-1 mesons mass-degenerate (appropriate for ${m_{\pi _{\textrm{D}}}}/\Lambda _\textrm{D} \gg 1$ and bordering on unphysical for the Lund model), we would expect probVector = 0.75 based on spin counting (three spin-1 states versus one spin-0 state.) Data from QCD, with $m_\pi /\Lambda _{QCD}\sim 0.5$, suggests use of probVector = 0.5, downweighting spin-1 presumably because of phase space.⁶ From this we learn that the appropriate probVector is a slowly increasing function of $m_\pi /\Lambda _{QCD}$. It would be reasonable to choose a phase-space-motivated functional form for this function, with a smooth ${m_{\pi _{\textrm{D}}}}\rightarrow 0$ limit, but we have not made an effort to do this. Little is known about the limit ${m_{\pi _{\textrm{D}}}}\rightarrow 0$; it is not even clear that the Lund model is accurate there.

When the parameter HiddenValley:separateFlav (included in the new version of pythia 8 and described in Sect. 4.1.4) is set “on”, the parameter HiddenValley:probKeepEta1 downweights the probability of producing a singlet $\eta '_\textrm{D}$ meson relative to other diagonal mesons. This should be set to 1 when $N_{c_\textrm{D}} \gg N_{f_\textrm{D}}$ since in the large-$N_{c_\textrm{D}}$ limit (with $N_{f_\textrm{D}}$ fixed) the axial anomaly is negligible and the $\eta '_\textrm{D}$ is like the adjoint-flavor bosons. Conversely it should be set to a small value when $N_{f_\textrm{D}}$ is of order or greater than $N_{c_\textrm{D}}$ and the $\eta '_\textrm{D}$ is heavy, as it is in QCD; in pythia 8, the corresponding QCD parameter StringFlav:etaPrimeSup is set by default to 0.12.

Finally, the option of allowing baryons in hadronization for $N_{c_\textrm{D}} =3$ can be controlled with the parameter HiddenValley:probDiquark; this determines the likelihood of pair-producing diquarks, which, for $N_{c_\textrm{D}} =3$ only, combine with a quark to form a baryon. We have not validated this parameter and recommend that for now baryons (at most a 10% effect, which is probably smaller than hadronization uncertainties) should not yet be used.

4.1.3 Decays of dark hadronic bound states

The dark hadronic bound states are either stable, undergo decays within the dark sector, or decay to final states that include SM particles. The decay patterns depend on the charge matrix $\textbf{Q}$ and the dark hadron mass hierarchy. In particular, when the mass $m_{\rho _\textrm{D}}$ is larger than twice ${m_{\pi _{\textrm{D}}}}$, the $\rho _\textrm{D}$ decay to $\pi _{\textrm{D}} \pi _{\textrm{D}}$. In the regime where such decays are not allowed some of the $\rho _\textrm{D}$ may decay back to the SM via mixing with the $\textrm{Z}^{\prime }$. The details of these decay modes however are determined by the group algebra and need careful treatment. We outline below the salient considerations in setting such decay modes.

Region 1: Dark sector decays $2{m_{\pi _{\textrm{D}}}} < m_{\rho _\textrm{D}}$

When the decay channel $\rho _\textrm{D} \rightarrow \pi _{\textrm{D}} \pi _{\textrm{D}}$ is open, it dominates all other decays since $g_{\rho _\textrm{D} \pi _{\textrm{D}} \pi _{\textrm{D}}}$ is large compared to any other coupling (Fig. 23 left). The width

$$\begin{aligned} \Gamma (\rho _\textrm{D} ^a\rightarrow \pi _{\textrm{D}} ^b\pi _{\textrm{D}} ^c ) \propto |f^{abc}|^2 \ \frac{g_{\rho _\textrm{D} \pi _{\textrm{D}} \pi _{\textrm{D}}}^2}{16\pi } m_{\rho _\textrm{D}}, \end{aligned}$$

(16)

where $f^{abc}$ are the structure constants of $SU(N_{f_\textrm{D}})$, is non-zero for all $\rho _\textrm{D}$ mesons, and large unless $N_{c_\textrm{D}}$ is enormous.⁷

Without any mixing between the $\textrm{Z}^{\prime }$ and the $\pi _{\textrm{D}}$, the latter is stable and invisible, so we assume that we are considering a model where such mixing occurs, allowing the decay

$$\begin{aligned} \pi _{\textrm{D}} \rightarrow (\textrm{Z}^{\prime })^* \rightarrow \textrm{q} \overline{\textrm{q}}. \end{aligned}$$

(17)

This mixing typically arises because a Higgs field (whose scalar is assumed too heavy to be of interest here) gives mass to both the $\textrm{Z}^{\prime }$ and the quarks $\textrm{q}_\textrm{D}$, typically along with some additional flavor violation. The exact details of the mixing and corresponding lifetimes depends on precise model-building.

Because the decay in Eq. 17 is helicity-suppressed, the width for this process is of order $|y_q|^2 {m_{\pi _{\textrm{D}}}}^5/m_{\textrm{Z}^{\prime }} ^4$ or smaller, where $y_q$ is the Yukawa coupling of the Standard Model quark; the heaviest kinematically-accessible quarks dominate. Note this width is parametrically small and low-mass $\pi _{\textrm{D}}$ will have displaced decays. To determine if a particular $\pi _{\textrm{D}}$ decays promptly, its lifetime needs to be calculated in a consistent leptophobic model, but to our knowledge the relevant model building has not been done.

Furthermore, we do not treat the flavor singlet $\eta '_\textrm{D}$ in detail here as our analysis is not yet complete. For $N_{f_\textrm{D}} \gtrsim N_{c_\textrm{D}}$ it is heavy, as in QCD, and rarely produced. Similarly, since baryons are only available for $N_{c_\textrm{D}} =3$, where they are a small effect, we do not discuss them here.

We conclude by noting that one should keep in mind that the $\textrm{Z}^{\prime }$ charge assignments also determines the decay branching fractions of the $\textrm{Z}^{\prime }$. Especially when dark hadron multiplicities are small in $\textrm{Z}^{\prime }$ decays, this introduces a small but significant correlation in the flavors of the dark hadrons, This will become more important in future studies with non-degenerate quark masses.

Fig. 23

Decay modes of diagonal and off-diagonal dark rho mesons for regime 1 (left hand side) and regime 2 (middle and right hand side)

Region 2: Dark sector decays $2{m_{\pi _{\textrm{D}}}} > m_{\rho _\textrm{D}}$

In this case we will assume that there is no mixing between the $\textrm{Z}^{\prime }$ and the $\pi _{\textrm{D}}$ – that the quarks and the $\textrm{Z}^{\prime }$ get their masses from separate sources. The hidden pions charged under $G_f$ are then stable and invisible. The precise fate of the other singlets needs further investigation. Some may be stable due to a discrete symmetry. Others are obviously unstable due to standard flavor anomalies, though the details depend on the matrix $\textbf{Q}$. Decays to one or two SM $\textrm{q} \overline{\textrm{q}}$ pairs have small widths because of powers of $1/m_{\textrm{Z}^{\prime }} ^2$ factors along with loop factors or phase space factors. In general these particles will be very long-lived on LHC detector time-scales, and thus a source only of . This statement may however be model-dependent and so one must be careful to compute the lifetimes for these states in a particular model. We assume here that all $\pi _{\textrm{D}}$ are LHC-detector stable.

The decays of the $\rho _\textrm{D}$, however, can be observed. First, there can be mixing between the $\textrm{Z}^{\prime }$ and the $\rho _\textrm{D}$ mesons which are singlets under the group $G_f$. These decays

$$\begin{aligned} \rho _\textrm{D} \rightarrow (\textrm{Z}^{\prime })^* \rightarrow \textrm{q} \overline{\textrm{q}} \end{aligned}$$

are not helicity suppressed and are thus faster than the corresponding $\pi _{\textrm{D}}$ decays that we discussed in the previous section (Fig. 23 middle).

For those $\rho _\textrm{D}$ that are non-singlet under $G_f$, flavor symmetry would not prevent the decay (Fig. 23 right)

$$\begin{aligned} \rho _\textrm{D} \rightarrow \pi _{\textrm{D}} + (\textrm{Z}^{\prime })^* \rightarrow \pi _{\textrm{D}}\ \textrm{q}\ \overline{\textrm{q}}. \end{aligned}$$

The $\rho _\textrm{D}$ and the $\pi _{\textrm{D}}$ in this decay have the same $G_f$ quantum numbers, while the $\textrm{q} \overline{\textrm{q}}$ are a flavor singlet. This decay would be prohibited by the naive symmetry $\pi _{\textrm{D}} \rightarrow -\pi _{\textrm{D}}$ in the chiral Lagrangian, but this symmetry is violated by the usual chiral anomaly that mediates $\pi _0\rightarrow \gamma \gamma$ decay in QCD, and allows a $\rho _\textrm{D} \rho _\textrm{D} \pi _{\textrm{D}}$ coupling in this context. Mixing between a $G_f$-singlet $\rho _\textrm{D}$ and a $\textrm{Z}^{\prime }$ then induces a $\rho _\textrm{D} \textrm{Z}^{\prime } \pi _{\textrm{D}}$ coupling, which permits this decay to proceed. Specifically

$$\begin{aligned} \Gamma (\rho _\textrm{D} ^a\rightarrow \pi _{\textrm{D}} ^b \textrm{q} \overline{\textrm{q}}) \propto |d^{abc}\textrm{Tr}(T^c\textbf{Q})|^2 \propto |\textrm{Tr}(\{T^a , T^b\} \,\textbf{Q})|^2 \end{aligned}$$

where $d^{abc}$ appears in the anti-commutator [187]

$$\begin{aligned} \{T^a, T^b\} = \frac{1}{N_{c_\textrm{D}}}\delta ^{ab} + d^{abc} T^c. \end{aligned}$$

Importantly, however, Tr$(\{T^a , T^b\} \,\textbf{Q})$ can vanish. If this occurs, then this decay channel is not available. For instance, if $T^a$ is the matrix whose $(\alpha ,\beta )$ entry is 1 and whose other entries are all zero, then the above trace is proportional to $Q_{\alpha }+Q_\beta$. Equal and opposite eigenvalues in $\textbf{Q}$ then assure that the corresponding $\rho _\textrm{D}$ does not decay via the anomaly. Although this does not guarantee that this particle is stable against decay via higher-order processes, it does mean that it has a very long lifetime and is likely LHC-stable.

Because this decay has a 3-body phase space and because by assumption $m_{\rho _\textrm{D}}- {m_{\pi _{\textrm{D}}}} < \frac{1}{2} m_{\rho _\textrm{D}}$, this decay is heavily suppressed, and will lead to displaced vertices if $m_{\rho _\textrm{D}}/m_{\textrm{Z}^{\prime }}$ is too small. In the limit $\Lambda _\textrm{D} \ll m_{\textrm{Z}^{\prime }}$ as we have assumed in this section,

$$\begin{aligned}{} & {} \Gamma (\rho _\textrm{D} ^a \rightarrow \pi _{\textrm{D}} ^b q \overline{q}) \nonumber \\{} & {} \quad =\frac{\left| \textrm{Tr}\left( \{T^a,T^b\}\textbf{Q}\right) \right| ^2 e^2_D g^2_D N_{c_\textrm{D}} ^2 m_{\rho _\textrm{D}} ^{11}}{5898240\,m_{\textrm{Z}^{\prime }} ^4 \,\pi ^7 f^6_{\pi _{\textrm{D}}}}F\left( \frac{{m_{\pi _{\textrm{D}}}}^2}{m_{\rho _\textrm{D}} ^2} \right) \nonumber \\{} & {} \quad \textrm{where} \ F(x)\equiv \ 1-15x-80 x^2+80 x^3+15 x^4-x^5\nonumber \\{} & {} \quad -60 (x+1) x^2 \log (x). \end{aligned}$$

(18)

In addition, we have used

$$\begin{aligned} |g_{\rho _\textrm{D} \rho _\textrm{D} \pi _{\textrm{D}}}| = \frac{N_{c_\textrm{D}} g^2_{\rho _\textrm{D} \pi _{\textrm{D}} \pi _{\textrm{D}}}}{8\pi ^2 f_{\pi _{\textrm{D}}}} \end{aligned}$$

(19)

and $g_{\rho _\textrm{D} \pi _{\textrm{D}} \pi _{\textrm{D}}} = m_{\rho _\textrm{D}}/(\sqrt{2}f_{\pi _{\textrm{D}}})$. The former relationship and in particular the factor of $N_{c_\textrm{D}}$ arises from $SU(N_{c_\textrm{D}})$ symmetry [188] while the latter is KSFR relationship [189, 190], and $e_D, g_{\textrm{q}}$ are defined in Eq. 12.⁸ For these decays to be prompt, the ratio $m_{\rho _\textrm{D}} ^{11}/m_{\textrm{Z}^{\prime }} ^4$ must not be too small. It should further be noted that $f_{\pi _{\textrm{D}}}$ also includes a mild $N_{c_\textrm{D}}$ dependence [181].

As above, we do not discuss the $\eta '_\textrm{D}$, the $\omega _\textrm{D}$ or baryons here.

4.1.4 Updates and inputs for PYTHIA 8 hidden valley module

The pythia 8 hidden valley module has received an update in version 8.307, after having been stable for some years. One update is substantive; the previous versions were overproducing very soft hidden hadrons (mainly pions) at low $p_{{\textrm{T}}}$. This bug fix slightly affects many plots, as we will see in Sect. 4.3; for example it affects the total multiplicity of hadrons. Fortunately, the methods used in previous SVJ analyses are not very sensitive to this effect, which leaves total visible energy in a jet, the aimed in its direction, and most substructure variables roughly unchanged. The possible exception is for variables which are not infrared safe, most notably $p_{{\text {T}}}D$; see 4.3.

The other main change has been to increase the flexibility of the module. Depending on a newly introduced flag separateFlav, the simulation in each regime may proceed in two ways. An imperfect but often sufficient simulation, which was already available in pythia 8.150, is available with separateFlav=off; in this case the full adjoint multiplets of spin-0 and spin-1 mesons are each simplified into two states, one flavor-diagonal and one flavor-off-diagonal. This division is not consistent with most choices of $\textbf{Q}$, as it requires that $G_f=U(1)^{N_{f_\textrm{D}}}$. However, as long as all dark hadrons are stable (on LHC timescales) or decay promptly, it is possible to mock up other choices of $\textbf{Q}$, where for instance only a fraction of the flavor-diagonal states would decay visibly, by assigning the flavor-diagonal meson a probability to decay to a visible SM state and a corresponding probability to decay invisibly.

Alternatively, the setting separateFlav=on allows full control over all the spin-0 and spin-1 states; separate lifetimes and decay modes can be assigned to each. The particle ID number for the spin-0 (spin-1) meson with quark i and anti-quark j is 4900ij1 (4900ij3), at least for $i\ne j$. For $i=j$ the situation is more complicated since the diagonal mesons are flavor mixtures; for example, with $N_{f_\textrm{D}} =3$, the pion is a $u\overline{u}-d\overline{d}$ state and the $\eta$ is a $u\overline{u}+d\overline{d}-2s\overline{s}$ state (ignoring normalizations). Typically one may order the diagonal flavor-adjoint mesons in a canonical way, through the increasing number of quark flavors appearing in their wavefunctions (or equivalently by the increasing number of non-zero entries in the corresponding diagonal $SU(N_{f_\textrm{D}})$ generator). The flavor-singlet state is always $(1/\sqrt{N_{f_\textrm{D}}})\sum _{\alpha } \textrm{q}_{\textrm{D} \alpha } \overline{\overline{\textrm{q}}_{\textrm{D} \alpha }}$ and is always assigned particle 4900FF1 (4900FF3) where $F\equiv N_{f_\textrm{D}}$. This is important because this state has special status, see below.

This setting then requires the user to create a full decay table for of order $N_{f_\textrm{D}} ^2$ dark hadrons. Although we will comment on the settings for this below, we have not yet automated this task, so at this time we have no benchmarks for separateFlav=on. It also permits the hidden quarks to have different masses, but we have not yet validated this capability and more studies are needed.

Other changes, not utilized below, are in the treatment of the flavor singlets, especially for spin-0, and baryons for $N_{c_\textrm{D}} =3$. Since the flavor singlets can be given different masses with separateFlav=on, this allows for a more accurate spectrum. The masses of the singlets should be assigned to the spin-0 and spin-1 particle ID codes 4900FF1 and 4900FF3. This is especially important for spin-0 because the singlet can have a much larger mass than the adjoint due to the axial anomaly, as for the $\eta ^{\prime }$ in QCD. As we mentioned above, an additional parameter probKeepEta1, which can be chosen between 0 and 1, has been added; this reduces the probability of producing of the $\eta '_\textrm{D}$ relative to other spin-0 mesons in the hadronization process. Meanwhile the routines for producing baryons in the SM sector have been activated for the HV sector as well, but only work for $N_{c_\textrm{D}} =3$. (For $N_{c_\textrm{D}} >3$ this is not a concern since baryon production would be highly suppressed. For $N_{c_\textrm{D}} =2$ a special routine must be written, because baryons, antibaryons and mesons are all degenerate; this is why the current HV module should not be used for $N_{c_\textrm{D}} =2$, at least not without careful consideration of how to reinterpret its results). For now, only one type of diquark is produced, that of $\textrm{q}_\textrm{D 1} \textrm{q}_\textrm{D 1}$. All the baryons produced are assumed to have spin 3/2 and to have one of the $N_{f_\textrm{D}}$ quark flavors i combined with a single flavor of diquark, with particle ID code 490i114. For separateFlav=off, all of these states are conflated into the state with $i=1$.

We have mentioned that pythia 8’s hadronization routine cannot simulate a theory with $N_{f_\textrm{D}} =1$ or $N_{c_\textrm{D}} =2$, but it may fail for other reasons. For any choices of $N_{f_\textrm{D}}$ and $N_{c_\textrm{D}}$, one should avoid overly small or large values of ${m_{\pi _{\textrm{D}}}}/\Lambda _\textrm{D}$. At small values approaching the chiral limit, theoretical understanding of hadronization is lacking, and the Lund string model used in pythia 8 may not function in any case; meanwhile at large values other hadrons (glueballs, in particular) will become as important as pions or rhos, but are not included in the Lund string model. To be conservative, we suggest limiting studies to $0.25< {m_{\pi _{\textrm{D}}}}/\Lambda _\textrm{D} < 2$ until there has been further theoretical work on this issue.

We now turn to the pythia 8 parameters that must be set to simulate the models discussed above. We begin with those that are independent of whether separateFlav=off or separateFlav=on.

• HiddenValley:Ngauge, HiddenValley:nFlav – These are $N_{c_\textrm{D}}$ and $N_{f_\textrm{D}}$; the former should always be set greater than 2 and the latter should always be set greater than 1. (For $N_{c_\textrm{D}} =2$ or $N_{f_\textrm{D}} =1$, pythia 8 is currently missing essential dark hadrons and gives an inaccurate simulation.)

• Constituent dark quark mass 4900101:m0 – The quark mass defined in pythia 8 is the constituent quark mass, not the current quark mass. This quantity has never been given a theoretical definition, but may be roughly defined by $m_{q_{const}} \approx m_{\textrm{q}_\textrm{D}} + {\mathcal {O}}(1)\times \Lambda _\textrm{D}$. For definiteness we will use this relation with the coefficient fixed to 1, namely $m_{q_{const}} \equiv m_{\textrm{q}_\textrm{D}} + \Lambda _\textrm{D}$.

• Confinement scale HiddenValley:Lambda – This can be defined in multiple ways, but we take it for now to be the scale at which the running gauge coupling constant diverges at 1-loop order, since currently the PYTHI8 HV module has implemented the running coupling at one loop. As we have mentioned above, further consideration of this definition is warranted. The associated behaviour is illustrated in Fig. 21.

• Lund model hadronization parameters must be set: HiddenValley:aLund, HiddenValley:bmqv2, HiddenValley:rFactqv, HiddenValley:sigmamqv; see Sect. 4.1.1. The effects of these parameters on the underlying phenomenology have not yet been investigated, so we make no specific recommendations for them beyond the existing default settings.

• Certain hadronization parameters must be set, such as HiddenValley:probVector; see Sect. 4.1.1.

Next, if separateFlav=off, only two additional parameters must be defined.

• Dark pion mass 4900111:m0, 4900211:m0. These are the masses of the bound state spin-0 multiplets; they should always be taken equal.⁹ Within the chiral regime, these may be related to the confinement scale $\Lambda _\textrm{D}$ and current quark mass $m_{\textrm{q}_\textrm{D}}$ via Eq. 14. However, we advise taking this observable as an input parameter and viewing $m_{\textrm{q}_\textrm{D}}$, which is scheme-dependent, as an output.

• Dark rho mass 4900113:m0, 4900213:m0. These are the masses of the bound state spin-1 multiplets, and should always be taken equal.¹⁰ Within the chiral regime, these are related to the confinement scale $\Lambda _\textrm{D}$ and the dark pion mass 4900111:m0 using Eq. 14.

In addition, decay channels and lifetimes for these four states must be defined by the user.

If instead separateFlav=on, then even for the mass-degenerate case, all spin-0 and all spin-1 mesons must have separately defined masses, 4900ij1:m0, 4900ij3:m0 for $N_{f_\textrm{D}} \ge i\ge j\ge 1$. Note the flavor singlets have particle ID codes 4900iis with $i=N_{f_\textrm{D}}$ and $s=1,3$; the user may wish to change probKeepEta1 which can be used to suppress the spin-0 singlet production. Again the user must define all lifetimes and decay channels, now for a much larger set of particles. Depending on the model, it may be very important to ensure that the flavor structure of the decays is precisely specified, as is emphasized in the earlier Eqs. 16 and 18.

4.1.5 Proposed benchmarks

We have created benchmarks for the purpose of the studies in Sect. 4.3. We are implicitly assuming dark hadron lifetimes are short enough to be considered prompt, as appropriate for the SVJ signatures. For low-mass dark hadrons, this is far from obvious. Lifetimes need to be calculated in the context of complete models, but constructing such models is no simple matter in the context of a leptophobic $\textrm{Z}^{\prime }$ because of potential $U(1)'$ gauge anomalies that would make the theory inconsistent. We are not aware of any complete calculations of dark hadron lifetimes in this context, so we must warn the user that some of these benchmarks, especially those with light $\pi _{\textrm{D}}$, may not be realizable theoretically.

Let us first note what all the benchmarks have in common. In each case

• $m_{\textrm{Z}^{\prime }} =1$ TeV;

• We take separateFlav=off.

Versions of the benchmarks with separateFlav=on would be more accurate in their treatment of flavor-singlets, but will have to be created at a later time.

We have several benchmarks with ${m_{\pi _{\textrm{D}}}}<\frac{1}{2} m_{\rho _\textrm{D}}$.

• All have $N_{c_\textrm{D}} =N_{f_\textrm{D}} =3$.

• All have ${m_{\pi _{\textrm{D}}}}=0.6 \Lambda _\textrm{D}$ (and thus $m_{\rho _\textrm{D}} =2.6\ \Lambda _\textrm{D}$ by Eq. 14).

• Because of this choice, the parameter probVec is taken to be 0.5, since ${m_{\pi _{\textrm{D}}}}/\Lambda _\textrm{D}$ is similar to its value used in real-world QCD.

• Three choices of $\Lambda _\textrm{D}$ are considered: 5 GeV, 10 GeV and 50 GeV.

• For each of these, the number of stable diagonal spin-0 mesons is $k=0$, 1 or 2, with $3-k$ decaying to the SM; since the six off-diagonal pions are stable in this model, this gives $r_{\text {inv}} =(6+k)/9.$

• The dark pions are assumed to decay promptly and only to $c\overline{c}$ (charm being the heaviest kinematically-allowed SM quark for the smaller values of $\Lambda _\textrm{D}$).

The choice of k depends on mixing among the singlet and diagonal adjoint pions and the $\textrm{Z}^{\prime }$. The details, especially the interplay between mixings and lifetimes, require careful model-building. We are not aware of any papers in which this has been done.

Note that the use of separateFlav=off means that we do not treat the $SU(N_{f_\textrm{D}})$ flavor singlets separately from the other mesons. For $N_{c_\textrm{D}} =3$ the $\eta '_\textrm{D}$ is much heavier than the other states, and this is not correctly modeled. In particular, it leads to a small correction to $r_{\text {inv}}$. For $N_{c_\textrm{D}} \gg N_{f_\textrm{D}}$, the splitting between the flavor adjoint and singlet states becomes small, so the use of separateFlav=off is less problematic there.

For ${m_{\pi _{\textrm{D}}}}<\frac{1}{2} m_{\rho _\textrm{D}}$, we have so far defined only one benchmark

• $N_{c_\textrm{D}} =3$, $N_{f_\textrm{D}} =4.$

• $\Lambda _\textrm{D} =10$ GeV, ${m_{\pi _{\textrm{D}}}}=17$ GeV, $m_{\rho _\textrm{D}} = 31.8$ GeV

• The parameter probVec is taken to be 0.58, in between the values of 0.5 (as used for QCD) and 0.75 (as appropriate for ${m_{\pi _{\textrm{D}}}}\approx m_{\rho _\textrm{D}}$).

• $\textbf{Q}= \{-1,2,3,-4\}$, a choice that ensures that no $\rho _\textrm{D}$ are stable.

• All spin-0 mesons are assumed to be stable on LHC-detector time-scales.

• All diagonal spin-1 mesons (including the singlet, which we do not treat carefully) decay to all available SM $\textrm{q} \overline{\textrm{q}}$ pairs.

• All off-diagonal spin-1 mesons decay to SM $\textrm{q} \overline{\textrm{q}}$ plus an invisible spin-0 meson.

Table 3 summarises our current benchmarks.

Table 3 Current benchmarks for ${m_{\pi _{\textrm{D}}}} > m_{\rho _\textrm{D}}/2$ and ${m_{\pi _{\textrm{D}}}} < m_{\rho _\textrm{D}}/2$ regimes. In the former case all $\pi _{\textrm{D}}$ are stable and source of , while for the later, the $\rho _\textrm{D}$ mesons decay to $\pi _{\textrm{D}}$ which further decay to $c\overline{c}$ final states at the LHC. The benchmarks assume that the decays of $\rho _\textrm{D}, \pi _{\textrm{D}}$ are prompt

To compose benchmarks with separateFlav=on, there are a number of additional steps needed. For ${m_{\pi _{\textrm{D}}}}<\frac{1}{2} m_{\rho _\textrm{D}}$, the decays $\rho _\textrm{D} ^a \rightarrow \pi _{\textrm{D}} ^b\pi _{\textrm{D}} ^c$ need to be correctly programmed. For example, for $N_{f_\textrm{D}} =3$, $\rho _\textrm{D} ^3$, the diagonal member of the rho isotriplet (particle ID 4900113), decays to spin-0 bosons $\pi _{\textrm{D}} ^{ij}$ (particle ID 4900ij1) in the following pattern:

$$\begin{aligned}{} & {} \rho _\textrm{D} ^3 (4900113) \rightarrow \pi _{\textrm{D}} ^{12}\pi _{\textrm{D}} ^{21} \ (66\%), \pi _{\textrm{D}} ^{13}\pi _{\textrm{D}} ^{31} \ (16\%),\\{} & {} \quad \pi _{\textrm{D}} ^{23}\pi _{\textrm{D}} ^{21} \ (16\%), \end{aligned}$$

the 4:1:1 branching ratios reflecting the relative isospins-squared of these spin-0 states. All of these details need to be correctly laid out in the pythia decay table in order that spin-0 mesons be produced in the right abundances. It is also important to decide how to treat the singlet states, especially the spin-0 singlet whose mass and production rate in hadronization may be quite different from the others. Finally, all the spin-0 meson decays and lifetimes must be separately entered into the pythia decay table.

For ${m_{\pi _{\textrm{D}}}}>\frac{1}{2} m_{\rho _\textrm{D}}$, similar efforts are required to ensure that the flavor structure of the diagonal and off-diagonal $\rho _\textrm{D}$ decays are correctly implemented in the decay table.

4.1.6 Final remarks

Before we proceed with phenomenological studies using the benchmarks proposed in this section, we would like to emphasize that we have only laid out an initial road for defining consistent phenomenology in the context of semi-visible jets. We have considered the leptophobic $\textrm{Z}^{\prime }$ SM–DS portal widely used in the semi-visible jets literature, and pointed out the crucial role of charge assignments in determining the phenomenology, though we have not worked out the details. Dark hadron masses may potentially be extracted from a combination of lattice simulation of the hidden sector and general theoretical considerations. But their decay channels and lifetimes are highly model-dependent, and calculating them involves careful consideration of the detailed charge assignments of the $\textrm{Z}^{\prime }$, its mixing with various dark hadrons, and the spectrum and interactions of the dark hadrons (including anomalies) as obtained from symmetry considerations and the chiral Lagrangian. These sometimes intricate calculations must be performed in each model, unless an over-arching theoretical treatment, covering all models in this class, can be given.

In the context of semi-visible jets, the lifetimes of the various states are particularly important. This signature is defined to be one in which all objects either are stable, producing , or decay promptly to SM-hadronic final states. Long-lived particles with lifetimes greater than a few centimeters and less than 10 meters (in the lab frame) would move the signature into a different regime, outside the semi-visible jet framework. It is therefore imperative to identify all unstable dark hadrons and calculate their lifetimes correctly. We have estimated lifetimes and have moderate confidence that all particles in our benchmarks decay promptly or are stable on LHC detector scales, but we have not by any means done a thorough analysis.

We would also like to note that there are still significant issues with hadronization that we have not begun to address. We have made a few observations about the hadronization parameters used in the pythia 8 HV module, but have neither attempted to explore the impact of their uncertainty on the underlying phenomenology, nor made concrete statements about what ranges of values they might take. These questions, and even deeper ones about how hadronization models perform in other regimes, such as the chiral limit (${m_{\pi _{\textrm{D}}}}\ll \Lambda _\textrm{D}$), must be left for future studies.

4.2 Improvements on the PYTHIA8 Hidden Valley Module and their validation

Contributors: Guillaume Albouy, Cesare Cazzaniga, Annapaola de Cosa, Florian Eble, Marie-Hélène Genest, Nicoline Hemme, Suchita Kulkarni, Stephen Mrenna, Ana Peixoto, Akanksha Singh, Torbjörn Sjöstrand, Matt Strassler

4.2.1 Sample generation

The signal process considered for the validation of new Hidden Valley (HV) Module of pythia 8 [192, 193] consists of semi-visible jets [2] produced in the s-channel via a heavy $\textrm{Z}^{\prime }$ mediator. A set of signal samples has been produced with different versions of pythia 8 for proton-proton collisions (and also electron-positron collisions for completeness) at the benchmark centre-of-mass energy of 13 TeV (1 TeV). Namely, in order to test the new implementation of the HV Module we have produced three main groups of samples as illustrated in Table 4 for different dark sector color charge $N_{c_\textrm{D}}$ and dark quark flavor $N_{f_\textrm{D}}$ choices.

Table 4 SVJ MC samples generated with pythia 8 for a $\textrm{Z}^{\prime }$ mass $m_{\textrm{Z}^{\prime }} = 1$ TeV. For all these categories, only the decay of the $\textrm{Z}^{\prime }$ to dark quarks is simulated

In particular, the first type of samples have been produced with an older pythia 8.245 release [194],¹¹ while the second one have been generated with the new pythia 8 version 8.307 [195].

In the new pythia 8 release, it is now possible to set the masses of all 8 dark quarks and associated 64 mesons for each pseudo-scalar and vector multiplet individually. Even if this allows to consider mass split scenarios, we consider only mass degenerate dark quarks since a consistent treatment for UV to IR settings in mass split scenario is not yet available. As an outcome of this choice, flavor symmetry leads to mass degenerate pseudo-scalar and vector multiplets. However, it is still crucial to have the possibility to set all dark mesons properties individually since the lifetimes of these different states can differ according to the model and mediator. Following these necessities, compared to pythia 8.245 (8.306) release, in the newest version a more detailed handling of dark hadrons is possible with the setting HiddenValley:SeparateFlav = on. As shown in Table 4, a third sample has been added in order to test this new option. In particular, using the flavor splitting option, each of the quark and meson flavors are shown explicitly. The quark names now are $\textrm{q}_{\textrm{D} i}$, with $i \in \{ 0, \ldots , N_{f_\textrm{D}} \}$. Similarly, meson names are $\pi _{\textrm{D} ij}$ and $\rho _{\textrm{D} ij}$, where $i = j$ are the flavor-diagonal mesons, and else $i > j$, with j representing the anti-quark. The identity codes then are 4900ij1 for pseudo-scalars and 4900ij3 for vectors. An anti-meson comes with an overall negative sign, and here i gives the anti-quark. The data tables by default contain identical properties for all diagonal mesons in a multiplet. All nondiagonal mesons of a multiplet are also assumed to be identical and stable by default.

An advantage of the SeparateFlav=on option, is the possibility of setting masses (as well as decay modes) of spin-0,1 flavour singlets differently than the corresponding multiplets. As discussed in Sect. 4.1, the exact computation of the flavor singlet mass with respect to the flavor multiplets, especially for spin-0 states, is an open question. There are indications that spin-0 singlets tend to be heavier than their multiplet counterparts, and therefore for these states a suppression of the production rate is also expected. For this reason, the option HiddenValley:probKeepEta1 can be set in pythia 8.307 in order to specify the suppression factor for the spin-0 flavor singlets production rates. This feature has been tested in the validation procedure, but we do not report plots related to this.

Fixing the color charge to $N_{c_\textrm{D}} =3$ using the option HiddenValley:Ngauge = 3, two configurations for the number of flavors $N_{f_\textrm{D}} =3,8$ have been considered in this study, using the setting HiddenValley:Nflav. $N_{f_\textrm{D}} = 3$ corresponds to the smallest possible configuration with more particles than the triplet representation used in pythia 8.245 , while $N_{f_\textrm{D}} = 8$ is the maximal number of flavors implemented in pythia HV module. Choosing these two values, we thus test the extremes of the flavor configurations. The hidden valley partners $F_D$ of the SM particles (charged both under both SM and hidden valley group) are assumed to be decoupled in our case, such that they will not produce interleaved showers between the hidden sector and the SM [192, 193]. While in the pythia 8.245 release only dark mesons originating from string fragmentation are implemented, in pythia 8.307 tested in this study an option to produce dark baryons has been added with the line HiddenValley:probDiquark = on. With this option, it is possible to set the probability that a string breaks by “diquark–antidiquark” production rather than quark–antiquark one. This then leads to an adjacent baryon–antibaryon pair in the flavor chain. Currently only one kind of diquark is implemented, implying at most eight different Delta baryons $\Delta _{\textrm{D} i}$ if HiddenValley:SeparateFlav = on. In the validation procedure of pythia 8.307 of the HV module we have considered decoupled Delta baryons.

A minimal number of input parameters have to be specified in pythia 8 when the Hidden Valley module is called with the option HiddenValley:fragment = on. In particular, the masses of the dark hadrons have to be fixed as well as the dark sector hadronization scale $\Lambda _\textrm{D}$ (set to 10 GeV in this study). Furthermore, the masses of pseudo-scalar states are set to 6 GeV, and the masses of the vector mesons are chosen to be 25 GeV. These settings correspond to $\pi _{\textrm{D}}/\Lambda _\textrm{D} = 0.6$ same as that considered in benchmarks in Sect. 4.1. The final states configuration that we chose for our study is simply a fully invisible signature where all the dark hadrons are considered to be stable. A further relevant setting which must be specified is the running of the dark sector coupling $\alpha _\textrm{D}$ which can be switched on with the option HiddenValley:alphaOrder = 1.

For the purposes of this study, for efficient MC generation, we consider a simplified scenario where the $\textrm{Z}^{\prime }$ mediator decays only to dark quarks, even if in a real physics case the non-vanishing coupling to SM-quarks contributes to the branching ratios. By default the $\textrm{Z}^{\prime }$ mediator nominal width $\Gamma _{\textrm{Z}^{\prime }}$ of the $\textrm{Z}^{\prime }$ boson is set to 20 GeV and the mass $m_{\textrm{Z}^{\prime }} = 1$ TeV. Figure 24 shows the invariant mass distribution for the $\textrm{Z}^{\prime }$ boson using the dark quarks before parton shower and hadronization in the hidden sector. The distribution deviates from the Breit–Wigner showing an excess of events in the low mass tail. This effect can be explained from the factorisation theorem considering that parton distribution functions blows up for low transferred momentum fractions for the SM incoming partons. Since we are only interested in typical events where a $\textrm{Z}^{\prime }$ boson is created to have a consistent comparison between the different samples, we choose to cut away the low mass tail requiring the generated invariant mass of the $\textrm{Z}^{\prime }$ to be within the range [800, 1200] GeV.

Fig. 24

Distribution of the $\textrm{Z}^{\prime }$ mediator mass from dark quarks (nominal value set in the simulation $m_{\textrm{Z}^{\prime }} = 1$ TeV and $\Gamma _{\textrm{Z}^{\prime }}=20$ GeV)

4.2.2 Validation plots

PYTHIA8 triplet implementation

In a two flavor theory there are 4 spin-0 states (and 4 spin-1 states); 1 diagonal and 2 off-diagonals, which make up the triplet, and an additional singlet. In the current pythia 8 release there are only 3 PIDs for dark pions (3 PIDs for $\rho _\textrm{D}$ mesons), which signify the positive and the negative off-diagonal and the diagonal dark pion. However, the singlet is still produced in pythia 8 and shares the same PID as the diagonal dark pion. With HiddenValley:SeparateFlav=off option, it is thus impossible to separate out the singlet: it is produced with the same probability as that of the diagonal dark pion. As the singlet is considered to be another diagonal dark pion in pythia 8, the ratio of diagonal to off-diagonal dark pions is 1:1 for $N_{f_\textrm{D}} =2$. In other words, pythia 8 will create an even amount of diagonal and off-diagonal dark mesons, and hence the PID for the diagonal dark pions ($\rho _\textrm{D}$ mesons), 111 (113), is equally as likely as the PIDs for off-diagonal dark pions ($\rho _\textrm{D}$ mesons) when considered together, 211 (213) and $-211$ $(-213).$ This is clearly illustrated in Fig. 25a. Similarly, in a theory with $N_{f_\textrm{D}} =3$ there is an octet and a singlet, or 3 diagonal and 6 off-diagonal dark mesons. In the current pythia 8 release there are also only 3 PIDs for these 9 dark pions (and 3 PIDs for 9 $\rho _\textrm{D}$ mesons), following the same logic as for $N_{f_\textrm{D}} =2$. The ratio of diagonal to off-diagonal dark mesons is now 1:2 so pythia creates twice as many off-diagonal dark mesons as diagonal ones. In this situation, 111 (113) is only half as likely as 211 (213) and $-211$ $(-213)$ together, see Fig. 25b.

Fig. 25

PdgID distributions for final-state dark particles without the 4900 prefix for a $N_{f_\textrm{D}} =2$ model and b $N_{f_\textrm{D}} =3$ model. FV splitting is turned off to simulate as the standard pythia 8 version. The parameter $\text {probVec}=0.75$ for both models

The pythia 8.307 includes individual PIDs for all the multiplets and singlets, as well as a new parameter called HiddenValley:probKeepEta1, which determines the probability to create the singlet state. This probability is set relative to the probability of producing spin-0 multiplets. The default setting is 1, but it can be set to 0 such that the singlet is not produced at all.

The handling of dark PIDs affects the expected value of $r_{\text {inv}}$. In pythia 8.245 release it is not possible to turn off the production of the singlet state and so this must be taken into account in the calculation of $r_{\text {inv}}$. Take as an example an $N_{f_\textrm{D}} =2$ model with diagonal $\rho _\textrm{D}$ mesons (113) promptly decaying to the SM through a vector portal. Firstly, the $\text {probVec}=0.75$ parameter dictates that 3/4 of the dark mesons will be $\rho _\textrm{D}$ mesons, of which half are diagonal. This means that 3/8 of the dark mesons will be unstable, while the remaining 3/8 off-diagonal $\rho _\textrm{D}$ mesons and 1/4 dark pions are stable, resulting in a ratio of stable dark mesons to all dark mesons of 5/8 or 0.625. The value of $r_{\text {inv}}$ can be calculated at the generator level by counting separately final-state, stable dark mesons and all dark mesons (including decayed $\rho _\textrm{D}$ mesons) in the event and taking the ratio of these two sums. The distribution of $r_{\text {inv}}$ for such a model with FV splitting turned off can be seen in Fig. 26.

PYTHIA full n-plate implementation

The validation of the new pythia HV module was performed through a phenomenological analysis of the distinct variables obtained for three different cases: HV module in pythia 8.245 and in pythia 8.307 with either HiddenValley:SeparateFlav = off or HiddenValley:SeparateFlav = on. All dark mesons were set to be stable. The distributions of angular and kinematic variables of the different final state particles were compared for those three cases. Although the most important variables are related to the dark pion and $\rho _\textrm{D}$ mesons, the missing transverse energy and the produced jets were also considered for this validation study. The reconstruction of the jets is done by clustering the generator level objects obtained after parton-shower and hadronization, using the radius parameter $\Delta R=1.4.$ As the dark mesons are set to be stable, the jets we study in this section are therefore not a result of hadronization in the dark sector. They originate e.g. from initial state radiation and subsequent hadronization in the SM sector. As mentioned before, the validation of the pythia 8.307 HV module is executed for two specific models: $N_{f_\textrm{D}} =3$ and $N_{f_\textrm{D}} =8$ (with $N_{c_\textrm{D}} =3$ for both). For simplicity, only the results from the pp analysis are shown, as the same conclusions were obtained for the $e^+e^-$ study.

$N_{c_\textrm{D}} = 3$, $N_{f_\textrm{D}} =3$ model

Fig. 26

Distribution of $r_{\text {inv}}$ for a model with $\text {probVec}=0.75$ and unstable diagonal $\rho _\textrm{D}$ mesons

Changing the $N_{f_\textrm{D}}$ value from 2 to 3 results in additional PIDs being produced, as can be seen by comparing the dark pions and $\rho _\textrm{D}$ mesons particle ID shown in Fig. 27a, b, respectively, to the ones shown in Fig. 25a, b. One can see that dark pions and $\rho _\textrm{D}$ with ID 311 (313) and $-311$ $(-313)$ are produced when setting $N_{f_\textrm{D}}$ to 3 and HiddenValley:SeparateFlav = on. For the case where HiddenValley:SeparateFlav = on in particular, a total of 9 pseudo-scalar and 9 vector particles are identified, with the same production rates for all states in a given multiplet. The distributions of the multiplicity and the transverse momentum of the dark hadrons are represented in Fig. 28a, b. A lower multiplicity of these dark hadrons and a softer transverse momentum can be seen for the HiddenValley:SeparateFlav = on scenario compared with pythia 8.245. These changes are expected with the new HV module due to the bug fix related to the newly implemented $p_{\text {T}}$ suppression for mini-string fragmentation, as discussed in Sect. 4.1.4. An overall agreement can be found for the HiddenValley:separateFlav = on and HiddenValley:separateFlav = off with the new HV module. Concerning the specific case of the dark pions, the corresponding multiplicity and transverse momentum can be found in Fig. 29a, b. With similar conclusions as for the dark pions, the distributions of the same variables corresponding to the $\rho _\textrm{D}$ mesons are shown in Fig. 30a, b. From Figs. 29a and 30a, it can be concluded that the pseudo-scalars have lower multiplicity with respect to vector mesons. The difference between the distributions for the diagonal and off-diagonal dark pions and $\rho _\textrm{D}$ mesons was studied. The multiplicity and the transverse momentum of the diagonal and off-diagonal dark hadrons were consistent with the previous conclusions, with an agreement between the new and old HV modules with the different HiddenValley:separateFlav options. For completeness, the distributions of the missing transverse energy and the minimum azimuthal angle between jets and missing transverse energy can also be found in Fig. 31a, b. The latter shows that the missing transverse energy is recoiling against jets, as expected in the fully invisible scenario investigated here. The use of the new HV module does not have any impact on the event kinematics, as expected.

Fig. 27

$N_{c_\textrm{D}} = 3$, $N_{f_\textrm{D}} =3$ model: a PdgId distribution for dark pions, b PdgId distribution for $\rho _\textrm{D}$ mesons

Fig. 28

$N_{c_\textrm{D}} = 3$, $N_{f_\textrm{D}} =3$ model: a distribution of the number of dark hadrons, b dark hadrons $p_T$ distribution. The bottom panels show the ratio of the distributions of the new HV module scenarios to the nominal one

Fig. 29

$N_{c_\textrm{D}} = 3$, $N_{f_\textrm{D}} =3$ model: a distribution of the number of dark pions, b dark pion $p_T$ distribution. The bottom panels show the ratio of the distributions of the new HV module scenarios to the nominal one

Fig. 30

$N_{c_\textrm{D}} = 3$, $N_{f_\textrm{D}} =3$ model: a distribution of the number of $\rho _\textrm{D}$ mesons, b $\rho _\textrm{D}$ meson $p_T$ distribution. The bottom panels show the ratio of the distributions of the new HV module scenarios to the nominal one

Fig. 31

$N_{c_\textrm{D}} = 3$, $N_{f_\textrm{D}} =3$ model: a Distribution of the generator level missing transverse energy , b distribution of the minimum azimuthal angle between jets and . The bottom panels show the ratio of the distributions of the new HV module scenarios to the nominal one

$N_{c_\textrm{D}} = 3$, $N_{f_\textrm{D}} =8$ model

Setting $N_{f_\textrm{D}} =8$ brings a whole new set of PIDs both for dark pions and $\rho _\textrm{D}$ mesons, as confirmed in Fig. 32a, b. For the case with HiddenValley:separateFlav = on, additional PIDs from 311 (313) and $-311$ $(-313)$ to 811 (813) and $-811$ $(-813)$ are produced with a total of 64 dark pions or $\rho _\textrm{D}$ mesons, with the same production rate for all states in each multiplet. The multiplicity and transverse momentum of the dark pions are shown in Fig. 33a, b and similarly, for the $\rho _\textrm{D}$ mesons in Fig. 34a, b. In agreement with the previous model analyzed, a lower multiplicity and a softer transverse momentum of the dark hadrons are observed with the new pythia 8 HV module. The same conclusions stand when looking at diagonal and off-diagonal dark pions and $\rho _\textrm{D}$ mesons separately. The missing transverse energy and the minimum azimuthal angle between jets and missing transverse energy can be found in Fig. 35a, b. Once again, these distributions agree for the three cases considered.

Fig. 32

$N_{c_\textrm{D}} = 3$, $N_{f_\textrm{D}} =8$ model: a PdgId distribution for dark pions, b PdgId distribution for $\rho _\textrm{D}$ mesons

Fig. 33

$N_{c_\textrm{D}} = 3$, $N_{f_\textrm{D}} =8$ model: a distribution of the number of dark pions, b dark pion $p_T$ distribution. The bottom panels show the ratio of the distributions of the new HV module scenarios to the nominal one

Fig. 34

$N_{c_\textrm{D}} = 3$, $N_{f_\textrm{D}} =8$ model: a distribution of the number of $\rho _\textrm{D}$ mesons, b $\rho _\textrm{D}$ meson $p_T$ distribution. The bottom panels show the ratio of the distributions of the new HV module scenarios to the nominal one

Fig. 35

$N_{c_\textrm{D}} = 3$, $N_{f_\textrm{D}} =8$ model: a distribution of the generator level missing transverse energy , b distribution of the minimum azimuthal angle between jets and . The bottom panels show the ratio of the distributions of the new HV module scenarios to the nominal one

Through this validation, we thus highlight some of the differences between pythia 8.245 and pythia 8.307 hidden valley module. We also demonstrate that switching HiddenValley:SeparateFlav = on or off does not lead to physics differences in the production rates or kinematics of the events, but allows to access additional meson PIDs whose masses and branching ratios can be manipulated according to theory predictions.

4.3 Phenomenological studies of jet substructure observables

Contributors: Cesare Cazzaniga, Florian Eble, Aran Garcia-Bellido, Nicoline Hemme, Nukulsinh Parmar

In this section we exemplify the kinematic distributions resulting from benchmarks proposed in Sect. 4.1, focusing on the benchmark with $\Lambda _D = 10$ GeV and $N_{f_\textrm{D}} = 3$, and belonging to the regime ${m_{\pi _{\textrm{D}}}} < m_{\rho _\textrm{D}}/2$ for which the $\rho _\textrm{D} \rightarrow \pi _{\textrm{D}} \pi _{\textrm{D}}$ decay mode is open. The mass of $\textrm{Z}^{\prime }$ boson is set to 1 TeV. We then consider either 1, 2 or 3 diagonal pions decaying to SM particles. The $\pi _{\textrm{D}}$ mesons decay to $c\overline{c}$ as this is the heaviest allowed fermion pair. We simulate this signal using ${\textsc {pythia}} 8.307$ with HiddenValley:separateFlav = off and pass it through DELHPES3 using the HL-LHC card. Jets are clustered using FastJet [196, 197] using anti-kt algorithm [198]. We produce 50k events for the distributions shown in Sect. 4.3.1, and 500k events for those shown in Sect. 4.3.2.

4.3.1 Basic kinematic distributions

As the $\rho _\textrm{D}$ mesons all decay within the dark shower in this benchmark, they are not included in the calculation of $r_{\text {inv}}$. Figure 36 shows the $r_{\text {inv}}$ parameter distribution. As expected, the 1-$\pi _{\textrm{D}}$ decay model has a an average $r_{\text {inv}}$ of $\simeq \frac{8}{9}$ as all $\rho _\textrm{D}$ mesons decay to $\pi _{\textrm{D}}$ and only 1 of the 9 $\pi _{\textrm{D}}$ is unstable. The 2-$\pi _{\textrm{D}}$ decays has a mean of $r_{\text {inv}} \simeq \frac{7}{9}$ and for the 3-$\pi _{\textrm{D}}$, a mean of $r_{\text {inv}} \simeq \frac{2}{3}$ is obtained.

Fig. 36

Comparison of $r_{\text {inv}}$ for 1-, 2- and 3-$\pi _{\textrm{D}}$ decay models

In Fig. 37, some basic kinematic variable distributions are compared for the 3 dark pion decay models. These generator-level distributions are computed with jets of radius $R=0.4$ and $p_{{\textrm{T}}} > 25$ GeV. The $p_{{\textrm{T}}}$ distribution in Fig. 37a shows that more dark pion decays result in a higher average lead jet $p_{{\textrm{T}}}$, as expected when more dark particles decay to SM particles that can be detected. The distribution shown in Fig. 37b reveals very similar values between the 3 different models, which may seem contrary to what one would expect, i.e. more SM-decaying pions might be expected to result in lower ; however, while more stable dark pions truly gives higher missing or invisible energy in the system, the additional invisible particles may be evenly distributed between the two back-to-back jets and therefore not appear in the detector as additional .

Figure 37c shows the distributions of the transverse mass, $M_{{\textrm{T}}}$, of the leading and sub-leading jet and the . As can be seen, having more SM-decaying dark pions generally yields a higher transverse mass. As the remains relatively stable but the jet $p_{{\textrm{T}}}$ increases with the number of unstable dark pions, this results in higher $M_{{\textrm{T}}}$ values.

Fig. 37

Comparison of kinematic variable distribution for 1, 2 and 3 dark pions decay: a $p_{{\textrm{T}}}$ distribution of the leading jet in the event, b distribution, c distribution of $M_{{\textrm{T}}}$ of the leading and sub-leading jets plus

4.3.2 Jet substructure consistency

Experimental searches and phenomenological studies for dark showers exploit jet substructure (JSS) observables to tag jets as dark jets [72, 74, 91, 199, 200]. Comparisons of jet suCohen:2020afvbstructure variables of interest, between the former and the new Hidden Valley pythia modules, between different dark vector meson production probabilities, and between different number of unstable dark pions $\pi _{\textrm{D}}$, are presented in this section. In the pythia 8 Hidden Valley module [201, 202] the probability to produce a dark vector meson can be changed by setting the parameter HiddenValley:probVector. There is no precise theoretical prediction for the fraction of dark vector mesons produced after string fragmentation in the hidden sector. Assuming a mass degeneracy between vector and pseudo-scalar states, it is reasonable to fix HiddenValley:probVector = 0.75 as pseudo-scalars have 1 degree of freedom while vector mesons have 3 degrees of freedom. However, generically the $\rho _\textrm{D}$ mesons and dark pions are not mass degenerate, hence the production rate of pseudo-scalars is enhanced compared to mass-degenerate scenarios due to the larger phase space available for lighter states. In this specific case, a reasonable value is HiddenValley:probVector = 0.5, very much like in QCD.

For this study, generator-level jets have been clustered with the inclusive anti-kt algorithm [198], choosing a cone size $R=0.8$ and a minimum $p_{{\textrm{T}}}$  of 200 GeV. Jets were clustered from all visible SM particles and jet constituents were used for computing the jet substructure. The JSS observables studied here are the generalized angularities $\lambda ^{\kappa }_{\beta }$, the N-subjettiness [203] $\tau _N$ and jet major and minor axes.

Fig. 38

Adapted from [205]

Visualization of the space of the generalized angularities $\lambda ^{k}_{\beta }$.

Generalized angularities are presented in Fig. 38 and are defined from the constituents $i \in \{ 1, \ldots , N \}$ carrying momentum fraction $z_i$ inside a jet of cone size R as:

$$\begin{aligned} \lambda ^{\kappa }_{\beta } = \sum _{i \in jet} z_i^{\kappa } \left( \frac{\Delta R_{i,jet}}{R} \right) ^{\beta }. \end{aligned}$$

(20)

N-subjettiness $\tau ^{\beta }_{N}$ are designed to count the number of subjets inside a jet. In specific, N-subjettiness is defined as:

$$\begin{aligned} \tau ^{\beta }_{N} = \sum _{i}p_{{{\text {T}}},i} {\text {min}}({R^{\beta }_{1,i}, R^{\beta }_{2,i}, R^{\beta }_{3,i}, \ldots , R^{\beta }_{N,i}}) \end{aligned}$$

(21)

where the sum is over the jet constituents, and $R^{\beta }_{N,i}$ is the distance between the Nth subjet and the ith constituent of the jet. $\tau ^{\beta }_{N}$ measures departure from N-parton energy flow: if a jet has N subjets, $\tau ^{\beta }_{N-1}$ should be much larger than $\tau ^{\beta }_{N}$. Originally, $\tau ^{\beta }_{N}$ have been introduced in order to identify hadronically-decaying boosted objects and reject QCD background. In those studies, the angular parameter $\beta$ has been fixed to 1 as done in previous studies for boosted objects discrimination [204].

The shape of the jet can be approximated by an ellipse in the $\eta - \phi$ plane. The major and minor axes are the two principal components of this ellipse and are defined from the following symmetric matrix M:

$$\begin{aligned} M = \begin{pmatrix} \sum _{i} p_{{\text {T}},i}^{2}\Delta \eta _{i}^{2} &{} - \sum _{i} p_{{\text {T}},i}^{2}\Delta \eta _{i}\Delta \phi _{i} \\ - \sum _{i} p_{{\text {T}},i}^{2}\Delta \eta _{i}\Delta \phi _{i} &{} \sum _{i} p_{{\text {T}},i}^{2}\Delta \phi _{i}^{2} \end{pmatrix} \end{aligned}$$

(22)

where the sum runs over all constituents of the jet and $\Delta \eta$, $\Delta \phi$ are the differences in $\eta$ and $\phi$ with respect to the jet axis. The major and minor axes are defined from the eigenvalues $\lambda _1$ and $\lambda _2$ of M as:

$$\begin{aligned} \sigma _{{\text {minor,major}}} = \sqrt{\frac{\lambda _{1,2}}{\sum _i p_{{\text {T}},i}^{2}}} \end{aligned}$$

(23)

Generalized angularities belong to the category of jet shape variables and they have been originally built to measure the quantity of radiation inside a jet in order to discriminate between jets initiated by quarks and those initiated by gluons [205, 206]. Indeed, for the gluon jets the values of the generalized angularities are usually expected to be larger since gluons are expected to radiate more due to the larger color factor. In the same way, these observables have been used in analyses to discriminate between SM jets and dark jets [98]. In particular, the dark jets are expected to be wider than SM jets due to the double hadronization process and the mass splitting between the dark bound states and the SM quarks.

We first start by comparing JSS observables between the old and new pythia Hidden Valley modules. Comparison of quark-gluon discriminant variables and N-subjettiness variables are shown in Figs. 39 and 40. Some systematic differences are observed for the jet transverse momentum dispersion $p_{{\text {T}}}D$, due to the different number and different $p_{{\textrm{T}}}$ spectrum of the dark mesons $\pi _{\textrm{D}}$ and $\rho _\textrm{D}$ in the new pythia Hidden Valley module. N-subjettiness are smaller with the new module when decreasing $r_{\text {inv}}$ and looking for high number of subjets. No large systematic difference is observed for the other substructure variables.

Fig. 39

Comparison of quark-gluon discriminant variables between old (grey solid line) and new (colored dashed line) Hidden Valley modules. The comparison is made for 1 (left), 2 (middle) and 3 (right) dark pions decay. The plotted ratio is the ratio of new to old pythia module

Fig. 40

Comparison of N-subjettiness variables between old (grey solid line) and new (colored dashed line) Hidden Valley modules. The comparison is made for 1 (left), 2 (middle) and 3 (right) dark pions decay. The plotted ratio is the ratio of new to old pythia module

Next, we studied the differences in the JSS observables between two dark vector meson production fractions: 50% and 75%. Comparison of quark-gluon discriminant variables, N-subjettiness and number of constituents are shown in Figs. 41, 42 and 43. Some systematic differences are observed for all variables. It is clear that the number of constituents in jets is lower for higher vector meson fraction. Jets with large number of soft constituents are characterized by low $p_{{\text {T}}}D$ while $p_{{\text {T}}}D$ is higher for jets where just a few constituents carry most of the momentum. The fact that $p_{{\text {T}}}D$ is higher for higher vector dark meson fractions is certainly an effect of the lower number of constituents. Jet girth, axes and N-subjettiness are all smaller in the case of probVector=0.75 compared to probVector=0.5. This indicates that jets are narrower since with larger values of vector mesons fraction we observe a harder $p_{{\textrm{T}}}$ spectrum for the dark hadrons decaying visibly.

Fig. 41

Comparison of quark-gluon discriminant variables between probVector=0.5 (grey solid line) and probVector=0.75 (colored dashed line). The comparison is made for 1 (left), 2 (middle) and 3 (right) dark pions decay. The plotted ratio is the ratio of probVector=0.75 to probVector=0.5

Fig. 42

Comparison of N-subjettiness between probVector=0.5 (grey solid line) and probVector=0.75 (colored dashed line). The comparison is made for 1 (left), 2 (middle) and 3 (right) dark pions decay. The plotted ratio is the ratio of probVector=0.75 to probVector=0.5

Fig. 43

Comparison of number of constituents in jets between probVector=0.5 (grey solid line) and probVector=0.75 (colored dashed line). The comparison is made for 1 (left), 2 (middle) and 3 (right) dark pions decay. The plotted ratio is the ratio of probVector=0.75 to probVector=0.5

We then studied how the number of unstable diagonal $\pi _{\textrm{D}}$ mesons affects the jet substructure. Plots of quark-gluon discriminant, number of constituents and photon energy fraction for different number of unstable diagonal dark pions are provided in Fig. 44. Multiplicity is higher for lower $r_{\text {inv}}$, which is expected as the multiplicity is directly related to the number of unstable dark pions. Major and minor axes as well as girth are higher for lower $r_{\text {inv}}$, suggesting that the jet is wider.

In conclusion, we have noticed that the variation of the hidden sector parameters such as probVector can impact JSS distributions at generator level leading to a harder spectrum for the dark hadrons and consequently narrower jets. Notably, only two benchmark points for the vector meson fraction have been investigated, and further studies are encouraged to understand better the impact of the parameters of the hidden sector on the observable JSS distributions.

Fig. 44

Comparison between different number of unstable diagonal dark pions

4.3.3 Infrared-collinear safety of JSS observables

Traditional calculations in perturbative quantum chromodynamics are based on an order-by-order expansion in the strong coupling $\alpha _s$. Observables that are calculable in this way are known as “safe” [207]. As it is well-known, divergences of different nature can appear in the perturbative series. For the ultraviolet divergences appearing in loop diagrams, since QCD is a renormalisable theory, such infinities can be consistently cured. Moreover, real-emission diagrams exhibit singularities in particular corners of the phase-space. More specifically, the singular contributions have to do with collinear splittings of massless partons and emissions of soft gluons off both massless and massive particles. Virtual diagrams also exhibit analogous infra-red and collinear (IRC) singularities and theorems [208,209,210] assure that such infinities cancel at each order of the perturbative series, when real and virtual corrections are added together, thus leading to physical transition probabilities that are free of IRC singularities. An observable ${\mathcal {O}}(\{p_i \})$ calculated from a system of particles with momenta $p_i$ is defined to be infrared safe if adding a soft particle with momentum $\epsilon$ the following relation holds:

$$\begin{aligned} {\mathcal {O}}(\{p_i \}) = \lim _{\epsilon \rightarrow 0} {\mathcal {O}}(\epsilon ,\{p_i \}). \end{aligned}$$

(24)

Instead, if we consider a particle $p_1$ splitting into 2 particles $p_1 \rightarrow p^{(a)}_1 + p^{(b)}_1$ with angle between them $\theta _{1,ab} \rightarrow 0$, the observable ${\mathcal {O}}(\{p_i \})$ is said to be collinear safe if:

$$\begin{aligned} {\mathcal {O}}(\{p_i \}) = \lim _{\theta _{1,ab} \rightarrow 0} {\mathcal {O}}(p^{(a)}_1, p^{(b)}_1 ,\{p_i \}). \end{aligned}$$

(25)

To check IRC safety of JSS observables, we computed them at different stages of the shower/hadronization going from the dark sector to the SM sector. For IRC unsafe observables, large fluctuations in the showering process are expected, while IRC safe observables should be more stable during the evolution. Therefore for collinear splittings or soft emissions happening during the parton shower, the IRC unsafe observables will tend to diverge from the original value calculated in previous stages of the showering. Due to this feature, the IRC unsafe observables if not validated on data can introduce important model dependence in analyses exploiting them in supervised classifiers . This is particularly relevant in the case of dark shower studies where the MC-data agreement for signal cannot be assessed, and therefore there is no real control on IRC unsafe observables due to the unknown details of the hidden sector (for example the dark hadronization scale $\Lambda _\textrm{D}$). Specifically, given an observable ${\mathcal {O}}(\{p_i \})$, changes in the Hidden sector parameters such as $\Lambda _\textrm{D}$ are expected to produce a power law scaling for IRC safe observables given the jet pt $p_{{\textrm{T}}j}$: $\langle \delta O_{safe} \rangle \sim (\Lambda _\textrm{D}/p_{{\textrm{T}}j})^{\alpha }$. On the other hand, the scaling is logarithmic in the case of IRC unsafe observables: $\langle \delta O_{unsafe} \rangle \sim \log (\Lambda _\textrm{D}/p_{{\textrm{T}}j})$. This means that depending on the hadronization scale of the dark sector, the IRC unsafe observables can undergo large fluctuations for $\Lambda _\textrm{D} \ll p_{{\textrm{T}}}$, which means that without knowing $\Lambda _\textrm{D}$ these observables are correctly described by the parton shower but are also dependent of the details of the hidden sector.

In this study we test IRC safety of JSS observables by calculating them at 3 levels in the evolution of the shower: dark sector hadrons, SM quarks and SM hadrons. As previously mentioned, we expect the IRC safe observables to fluctuate less in the evolution. For the test we consider two generalized angularities, namely $p_{{\text {T}}}D$ which is an IRC unsafe observable, and the jet girth, which is IRC safe. The collinear unsafety of $p_{{\text {T}}}D$ is due to its dependence on the squared of the transverse momenta of the jet constituents. Therefore, taking a particle with transverse momentum $p_{1,{\textrm{T}}}$, if the particle splits into 2 particles with transverse momenta $p^{(a)}_{1,{\textrm{T}}}$ and $p^{(b)}_{1,{\textrm{T}}}$, $p_{{\text {T}}}D$ becomes:

$$\begin{aligned} (p_{{\text {T}}}D)^2 \propto (p_{1,{\textrm{T}}})^2 = (p^{(a)}_{1,{\textrm{T}}})^2 + (p^{(b)}_{1,{\textrm{T}}})^2 + p^{(a)}_{1,{\textrm{T}}} p^{(b)}_{1,{\textrm{T}}} cos(\theta _{12}). \end{aligned}$$

(26)

Therefore, we expect $p_{{\text {T}}}D$ to fluctuate more during the showering compared to the jet girth. Our results for the test of IRC safety for JSS observables is presented in Fig. 45. The plots show the following ratios for the tested JSS observable: unstable dark hadrons vs SM quarks, SM quarks vs SM hadrons and unstable dark hadrons vs SM hadrons. We expect the distributions of the ratios for the collinear unsafe observable calculated at different steps of the shower to differ from unity. For a fair comparison between the same observable calculated at different stages of the showering we consider only jets with a multiplicity of SM quarks which is twice the dark hadrons one. Moreover, because the girth of jets with one constituent is a special case as girth is close to 0, we consider only jets with a number of unstable dark hadrons strictly larger than one. The main result of this study is that even if IRC unsafe observables are expected to be described quite well by the parton shower, the application of IRC unsafe observables in the context of dark shower searches should be carefully validated in control regions by comparing Monte-Carlo and data. Secondly, as the dark hadronization scale is unknown, the effect of changing $\Lambda _\textrm{D}$ on JSS observables must be evaluated. The usage of such variables especially in Hidden Valley searches can lead to important limitations in terms of interpretability of the results due to their strong dependence on the unknowns of the Hidden sector.

Fig. 45

Test of IRC safety. Ratio of jet substructure variables (top: girth, bottom: $p_{{\text {T}}}D$) computed at three different levels: unstable dark hadrons, SM quarks from dark hadrons and SM hadrons. Large variations are observed for $p_{{\text {T}}}D$, which is IRC-unsafe, while ratios of girth, which is IRC-safe, peak at 1

4.3.4 Study of JSS observables after jet reconstruction in Delphes

After checking how the different parameters of the model affect the generator-level jets, we perform a similar study at reconstructed level using Delphes output. This is important to understand the impact of detector effects on the JSS observables that can be used by the experiments to tag dark jets efficiently. Delphes was configured for a CMS-like detector at the HL-LHC, and in particular Particle Flow candidates have been clustered with four different distance parameters, $R = 0.4,$ 0.8, 1.0 and 1.2, using FastJet [196, 197]. Jets with larger radius help in containing more of the radiation of the dark jet. Jets are required to have at least two tracks, and $\vert \eta |< 2.5$ and a minimum $p_{T}$ for clustering of 25 GeV. Figures 46 and 47 show the difference between the samples with probVector=0.5 and 0.75 when the jets are clustered with $R=0.8$.

Fig. 46

Comparison of reco-level variables between probVector=0.5 and probVector=0.75 for different number of unstable diagonal dark pions. The first row shows the number of reconstructed jets, the second row shows the constituent multiplicity of all jets, and the third row shows the girth of all jets in the event. The plotted ratio is the ratio of probVector=0.75 to probVector=0.5

Fig. 47

Comparison of reco-level variables between probVector=0.5 and probVector=0.75 for different number of unstable diagonal dark pions. The first row shows the minor axis, the second row shows the major axis, the third row shows the n-subjettiness ratio $\tau _{21}$, and the fourth row shows $\tau _{32}$. The plotted ratio is the ratio of probVector=0.75 to probVector=0.5

Figures 48 and 49 show the effect of varying the number of unstable diagonal pions on the JSS observables, for different distance parameters and probVector=0.5. The variables $p_{{\text {T}}}D$ and the N-subjettiness ratios show the most discrimination between the different samples.

4.3.5 Conclusion

Setting the IR parameters in accordance with the UV physics in general leads to a more cohesive modelling of the signal. This modelling however necessarily suffers from uncertainties due to a lack of knowledge of the precise hadronization parameters. These parameters can be varied to understand their effect on the resulting kinematic observables. In this section we have considered several jet substructure variables. We illustrated that changes in probVector can lead to changes in the observed jet substructure variables. It should be noted that this study concentrates only on one specific benchmark point and two values of probVector settings. It nevertheless shows the importance of understanding the effects of hadronization uncertainties. We also discussed the importance of Infrared and Collinear (IRC) safety when using substructure variables and in particular demonstrated that $p_{{\text {T}}}D$ is not IRC safe. Our studies thus highlight the need of a more detailed analysis of widely used jet substructure techniques in the light of dark showers phenomenology.

Fig. 48

Comparison of reco-level variables between different number of unstable diagonal dark pions for probVector=0.5. Left column show the $R=0.4$ jets and right column $R=0.8$ jets. First row is the number of reconstructed jets, second row is the constituent multiplicity of all jets, third row is the $p_{{\text {T}}}D$ of all jets, and the fourth row is the girth of all jets in the event

Fig. 49

Comparison of reco-level variables between different number of unstable diagonal dark pions for probVector=0.5. Left column show the $R=0.4$ jets and right column $R=0.8$ jets. First row is the axis minor, second row is the axis major, third row is the n-subjettiness ratio $\tau _{21}$ and fourth row is $\tau _{32}$

5 Improved search strategies

The wide variety of signatures coming from the dark/hidden sector scenarios considered throughout this work also motivates advanced techniques which may enable us to distinguish between signal and background at the LHC. These techniques may involve new kinematic variables, jet substructure information (as briefly discussed in Sect. 4.3), machine learning, or advanced triggering strategies. In this section, we illustrate some the avenues which have been explored in the literature, using the dark/hidden sector parametrizations presented in Sect. 2.2. It would be of great interest to also perform such analyses in the light of the new developments presented in Sect. 4.1.

5.1 Event-level variables

Contributors: Hugues Beauchesne, Giovanni Grilli di Cortona

Semi-visible jets are a characteristic signature of many confining dark sectors and consist of jets of visible hadrons intermixed with invisible stable particles. Up to now, two main search strategies have been pursued: tagging semi-visible jets (see e.g. Refs. [72,73,74, 199, 200, 211,212,213]) and exploiting the special relation between the azimuthal direction of the semi-visible jets and the missing transverse momentum (see e.g. Refs. [2, 63, 214]). In Ref. [215], it was shown that these two approaches can be combined to define new event-level variables that considerably increase the sensitivity of semi-visible jet searches. The central idea is that semi-visible jets are responsible for most of in signals and that tagging specifies which jets are semi-visible. The tagging information then predicts the direction and magnitude of , which can be compared to its measurement. In this section, we present a summary of Ref. [215] and refer to it for technical details.

For illustration purposes, consider the following benchmark model. Assume a new confining group ${\mathcal {G}}$. Introduce a dark quark $\textrm{q}_\textrm{D}$ that is a fundamental of ${\mathcal {G}}$ and neutral under the Standard Model gauge groups. Introduce a scalar mediator S that is an antifundamental of ${\mathcal {G}}$ and has an hypercharge of $-1$. These fields allow the Lagrangian

$$\begin{aligned} {\mathcal {L}} = \lambda _i S^\dagger \overline{\textrm{q}}_\textrm{D} P_R E_i + \text {h.c.}, \end{aligned}$$

(27)

where $E_i$ are the Standard Model leptons. Assume for simplicity that the only non-negligible $\lambda _i$ is the one corresponding to the electron. If the mediators are pair-produced, they will each decay to an electron and a dark quark. The experimental signature will then be two electrons and two semi-visible jets. This is similar to the signature of leptoquark pair-production and as such preselection cuts are applied based on typical leptoquark cuts. The event is also required to contain two jets tagged as semi-visible. We focus on the $\textrm{t}$ $\overline{\textrm{t}}$ background. Events are generated using MadGraph5 [216], pythia  8 [217] and Delphes 3 [218]. The Hidden Valley module of pythia is used with the following parameters:

Setting	Value	Setting	Value
NGauge	3	Dark pion mass	10 GeV
nFlav	1	Dark rho mass	21 GeV
FSR	On	pTminFSR	11 GeV
alphaOrder	1	fragment	On
Lambda	10 GeV	probVec	0.75
Dark quark mass	10 GeV

All other parameters are left to their default value. Finally, we define $1 - r_{\text {inv}}$ as the average fraction of the dark pions that decay back to Standard Model particles.

Consider a signal event. Label the transverse momenta of the two dark quarks produced from the decay of the two S as ${\textbf{p}}_{\textrm{T}}^{\textrm{q}_{\textrm{D} i}}$, with $i \in \{1, 2\}$. Each $\textrm{q}_\textrm{D}$ leads to a visible jet of transverse momentum ${\textbf{p}}_{\textrm{T}}^{D_i}\sim (1 - r_{\text {inv}}) {\textbf{p}}_{\textrm{T}}^{\textrm{q}_{\textrm{D} i}}$ and a contribution to of $\sim r_{\text {inv}} {\textbf{p}}_{\textrm{T}}^{\textrm{q}_{\textrm{D} i}}$. This gives

(28)

Consider the decomposition . The coefficients $a_1$ and $a_2$ should then peak at $\sim r_{\text {inv}}/(1 - r_{\text {inv}})$ and can be combined in a single test statistics. This could be done in multiple ways, but a simple and powerful one is to train a fully supervised neural network on the $a_1$ and $a_2$ of both the signal and the background. Alternatively, one can encode much of the same reasoning in a single variable. Define

(29)

where $\phi _{{\textbf{p}}_{\textrm{T}}^{D}}$ () is the azimuthal angle of ${\textbf{p}}_{\textrm{T}}^{D_1} + {\textbf{p}}_{\textrm{T}}^{D_2}$ (). This quantity should peak at 0 for the signal, but unfortunately contains no information on the norm of . We introduce two comparisons. First, the standard procedure up to now has been to compute the minimal difference in azimuthal angle between and the leading jets [63]

(30)

where in this case four jets are considered. Second, we consider a supervised neural network using . This is only meant as a comparison, as fully supervised neural network are susceptible to simulation artefacts and sculpting.

Fig. 50

ROC curves for $\Delta \phi$ (blue), $\Delta \phi _{CLLM}$ (orange) and for the neural networks using the input variables $a_1$ and $a_2$ (green) or the set x (red) for a mediator mass (here called L) of 500 GeV. Taken from Ref. [215]

Receiver Operating Characteristic (ROC) curves are shown in Fig. 50 for different values of $r_{\text {inv}}$. As can be seen, the coefficients $a_1$ and $a_2$ typically provide the strongest results. They sometimes exceed the fully supervised neural network by exploiting information on the magnitude of the momenta which are not provided to the neural network. The coefficients outperform the standard approach of $\Delta \phi _{\text {CLLM}}$ by an order of magnitude for a signal rejection rate of 0.5. The variable $\Delta \phi$ also generally outperforms $\Delta \phi _{\text {CLLM}}$.

5.2 Casting a graph net to catch dark showers

Contributors: Elias Bernreuther

To increase the sensitivity to dark shower signals consisting of promptly decaying dark hadrons, it is crucial to reduce the large QCD background. While backgrounds from mismeasured QCD jets mimic the signal with regards to event-level observables, such as $\Delta \phi$, differences are expected at the level of jet substructure. These can arise from differences in the shower evolution between QCD and the dark sector, the presence of visibly decaying heavy dark mesons in the jets, or invisible dark hadrons that are interspersed with visible particles. See e.g. Refs. [74, 200] for recent studies of dark shower signals in terms of classic jet substructure variables. In contrast, advances in tagging jets with modern machine learning techniques make use of low-level properties of jet constituents. Here, we summarize the results of Ref. [73], which studies the potential of deep neural networks for identifying semi-visible jets from dark showers.

As a benchmark, dark showers of nearly mass-degenerate GeV-scale dark mesons which are produced at the LHC via a heavy $\textrm{Z}^{\prime }$ vector mediator with mass on the TeV scale were considered. The underlying dark sector is the Aachen model summarized in Sect. 2.2.6 and motivated by cosmological and experimental constraints [61]. The dark quark production process $p p \rightarrow \textrm{q}_\textrm{D} \overline{\textrm{q}}_\textrm{D}$ was simulated with MadGraph5 2.6.4 [216] using a UFO file generated with FeynRules [219] and performing MLM matching with up to one additional hard jet. Showering and hadronization, both in QCD and in the hidden sector, were carried out using pythia  8.240 [46, 47, 220]. The settings used in pythia ’s Hidden Valley module for a signal with dark meson mass $m_\textrm{D}$ are summarized in Table 5. The parameter probVector was set to 0.5 such that 25% of dark mesons are unstable, flavor-diagonal vector mesons as predicted by the benchmark model. Jet clustering is performed by FastJet [196] using the anti-$k_T$ algorithm with jet radius $R=0.8$.

Table 5 Settings of the pythia Hidden Valley module used for generating the dark shower signal in Ref. [73]

A priori, it is not clear what the optimal jet representation and neural network architecture are to optimally distinguish dark shower jets from QCD jets. In Ref. [73] it was shown that a dynamic graph convolutional neural network (DGCNN) [221, 222] operating on particle clouds outperforms convolutional neural networks (CNNs) based on jet images [223] and a network operating on ordered lists of Lorentz vectors [224]. While a standard CNN carries out convolutions over neighboring pixels in a jet image, a DGCNN performs convolutions over edges of a graph constructed from jet constituents that are neighbors in feature space. While graph networks also represent the state of the art in tagging boosted top jets [225], their advantage over a CNN or a Lorentz Layer network is considerably larger in identifying semi-visible jets. A comparison of ROC curves showing the QCD jet background rejection $1/\epsilon _B$ as a function of the dark shower signal efficiency $\epsilon _S$ for $m_\textrm{D} = 5$ GeV is shown in Fig. 51.

Fig. 51

Figures taken from Ref. [73]

Left: ROC curves showing the semi-visible jet tagging efficiency $\epsilon _S$ and QCD background rejection $1/\epsilon _B$ for a DGCNN compared to a CNN and a LoLa network operating on jet images and Lorentz vectors, respectively. Right: ROC curves for a DGCNN trained on a mixed sample containing a number of different dark meson masses $m_\textrm{D}$ (here called $m_{\textrm{meson}}$) and tested on dark showers with $m_\textrm{D}$ as stated in the legend (dashed lines), compared to a DGCNN trained and tested on dark showers with identical $m_\textrm{D}$ (solid lines).

Since the parameters of the dark sector are a priori unknown it is a crucial question how well the classification performance of the DGCNN generalizes to dark showers with different parameter values than were used for training. Varying $r_{\text {inv}}$ and $m_\textrm{D}$, the performance continuously degrades the further the parameters of the dark showers in the test sample are from those in the training sample. While the effect is modest for $r_{\text {inv}}$, it is much more substantial for the dark meson mass. For example, for a network trained with $m_\textrm{D} =5$ GeV, the background rejection rate for signal efficiencies between 0.1 and 0.3 is reduced by nearly an order of magnitude when tested on samples with $m_\textrm{D} =20$ GeV. This suggests that the network learns to reconstruct this mass from the jet constituents. Importantly, this behavior can be mitigated by training the network on mixed samples which contain jets with a range of different dark meson masses. This yields a much more general classifier as reflected in the ROC curves in Fig. 51.

Finally, it was investigated how much the sensitivity of an experimental search for dark showers can be improved by applying a DGCNN as a semi-visible jet tagger. As an example, an ATLAS search for mono-jet events with a luminosity of 36.1 fb$^{-1}$ [226] was considered, which is sensitive to signal events where one of the two dark showers remains invisible and, thus, $\Delta \phi \approx \pi$. For an event to be accepted, it had to fulfil the original selection criteria of the search and contain at least one fat jet that is classified as a semi-visible jet by the network. The training sample consisted of jets from a dark shower signal with the benchmark parameters stated in Sect. 2.2.6 and from the dominant $\text {Z}$+jets background. The expected number of background events with and without the DGCNN tagger is shown in Table 6 for the signal region EM4, which is the region most sensitive to the signal when $m_{\textrm{Z}^{\prime }} = 1$ TeV. In addition, the table compares the resulting expected 95% CL limit $S^{95}_{\textrm{exp}}$ on the number of signal events in the region with and without the tagger and shows the corresponding improvement of the projected limit on the dark quark production cross section. In the benchmark scenario shown in Table 6 a DGCNN for tagging semi-visible jets can improve the sensitivity of the search to dark showers by more than one order of magnitude.

Table 6 Number of background events B with systematic uncertainty in the signal region EM4 of the search with and without the dark shower tagger and corresponding expected 95% CL limit $S^{95}_{\textrm{exp}}$ on the number of signal events. In addition, the improvement in the limit on the dark quark production cross section for the benchmark scenario described in the main text is shown relative to the search without a tagger. Table adapted from Ref. [73]

5.3 Autoencoders for semi-visible jets

Contributors: Annapaola de Cosa, Jeremi Niedziela, Kevin Pedro

Semi-visible jets arise from Hidden Valley models of dark matter, which include strong interaction in the dark sector. They constitute a challenging experimental signature in which a fraction of jet constituents is invisible to the detector, leading to missing transverse energy being aligned with the jet.

The details of the kinematics are mainly affected by the following theory parameters: $m_{\textrm{Z}^{\prime }}$ (the mass of the mediator), $m_\textrm{D}$ (the mass of the dark hadrons) and $r_{\text {inv}}$ (the fraction of stable, invisible dark hadrons). However, a large total number of unknown theory parameters leads to a vast model space with a huge number of possible scenarios that can easily evade any constraints from e.g. cosmological measurements. Since it is impractical to perform dedicated searches for all possible model variations, we propose to use autoencoders (AE) as anomalous jets taggers instead [199].

The autoencoder-based anomaly detection strategy is robust against both detector effects and details of the model implementation. AEs are designed to detect objects significantly different from the training sample, without prior knowledge of signal characteristics. For reference, the AE introduced here is compared to a Boosted Decision Tree (BDT) trained on the QCD background and a mixture of different signals. For completeness, we have also studied alternative anomaly detection techniques, namely Variational Autoencoders (VAE) and Principal Component Analysis (PCA).

All architectures mentioned above were trained on high-level properties of jets: $\eta$ and $\phi$ coordinates and invariant mass $m_{j}$, as well as jet substructure variables: jet $p_{{\textrm{T}}}$ dispersion $p_{{\text {T}}}D$, jet ellipse minor and major axes, EFP₁, and ECF ratios: C₂ and D₂. We have also considered including four-momenta of jet constituents in the training, but they were ultimately discarded since no improvement was observed.

Fig. 52

This figure is reproduced from Ref. [199]

Left and middle panels: comparison of AUC values from the autoencoder and a BDT trained on a mixture of all signal models with $m_\textrm{D} = 20\,\textrm{GeV}$. Right panel: AUC values for a BDT trained on a signal with parameters different from the signal used for testing, as indicated by the arrows. For example, the AUC value presented in the top left corner of the table comes from a model trained on the sample from the lower right corner.

The performance of different approaches is quantified by comparing the area under the ROC (receiver operator characteristic) curve (AUC), shown in Fig. 52. It was demonstrated that an AE-based jet tagger can provide satisfactory performance, compared with the fully supervised BDT approach. The PCA proved to be less efficient than other approaches. The VAE was found give the best results when trained exclusively on reconstruction loss, leading its variance to collapse to zero and therefore becoming equivalent to a regular AE.

Fig. 53

This figure is reproduced from Ref. [199]

Comparison of AUC values of the autoencoder and BDT with varying $m_\textrm{D}$ values (here called $m_{\text {dark}}$). The BDT was trained on a mixture of all signals with $m_\textrm{D} = 20~\textrm{GeV}$.

Robustness against unknown model parameters was also assessed. As shown in the rightmost panel of Fig. 52 and in Fig. 53, in certain cases the AE can outperform the BDT when the latter was trained on an incorrect signal hypothesis. Another interesting observation that can be made in the right panel of Fig. 52 is that a BDT trained on $r_{\text {inv}}$ = 0.3 and tested on $r_{\text {inv}}$ = 0.7 performs better then the one trained on a mixture of different signals (left panel). This is caused by the fact that the low $r_{\text {inv}}$ signal is more similar to the background, and therefore the BDT has to learn how to distinguish between the two more precisely. This results in a performance boost when tested on an easier case of large $r_{\text {inv}}$.

5.4 Autoencoders for SUEP

Contributors: Jared Barron, David Curtin, Gregor Kasieczka, Tilman Plehn, Aris G.B. Spourdalakis

5.4.1 Searching for hadronic SUEP

The theoretical motivation and experimental phenomenology of SUEPs are described in Sect. 3.1. Strategies to overcome the experimental challenges of searches for SUEP at the LHC are still being developed. For the nightmare scenario of prompt, hadronically decaying SUEP, a search strategy was proposed in [159], employing an autoencoder neural network as an anomaly detector.

An autoencoder is an unsupervised neural network trained on background events, which attempts to minimize the difference between its output and input. Ideally the autoencoder learns to do this efficiently only for inputs that are similar to its training data, so that when evaluated on an event from outside the background distribution, a high reconstruction error flags the event as anomalous. In the case of a search for SUEP, the background events are soft, highly isotropic QCD events. The unsupervised nature of this analysis avoids the model dependence that comes from using signal simulation to develop a classifier.

5.4.2 Signal generation

While the use of unsupervised machine learning techniques removes the need for signal events in the training dataset, a simulated signal dataset is still necessary to evaluate the autoencoder’s performance as an anomaly detector. For this purpose, SUEP events were generated using a statistical toy model of the dark shower.

The highly isotropic hadronic SUEP toy model simulated events were generated beginning with the production of Higgs bosons in association with a $\textrm{W}$ or $\text {Z}$ boson, simulated at center-of-mass energy 14 TeV in pythia 8 [220]. The vector boson was then required to decay leptonically. The hard lepton(s) from the vector boson decay were used to sidestep the issue of how to trigger on SUEP for this analysis. The decay of the Higgs to a shower of dark mesons was performed with the SUEP_Generator plugin in pythia 8.243, which models the dark shower as being a completely isotropic cloud with Boltzmann-distributed momenta as was presented in Eq. 10, and for which the parameter $T_\textrm{D}$ controls the energy distribution of the dark mesons, and represents the Hagedorn temperature of the dark sector. Only one flavor of dark meson is assumed, with mass $m_\textrm{D}$. Each dark meson was then forced to decay hadronically to a $u\overline{u}$ quark pair. From this point the parton showering and hadronization were performed by pythia as normal. Signal simulation was generated for $m_\textrm{D}$ from 0.4 to 8 GeV, and for $T_\textrm{D}/m_\textrm{D}$ from 0.25 to 4. Detector simulation was performed with Delphes 3 with CMS detector settings [218]. Due to the difficulty of disentangling the highly diffuse energy depositions of SUEP from pile-up, only charged track information was used for the analysis.

5.4.3 Background generation

The simulated QCD background events necessary for the training and test datasets were created by generating di-jet plus lepton(s) events with a reduced jet $p_{{\textrm{T}}}$ threshold of 15 GeV in MadGraph5_aMC@NLO 2.6.6 with hadronization by pythia 8 and detector simulation by Delphes 3 [216, 227].

5.4.4 Analysis

A trigger-level selection was applied to all simulated events requiring at least one charged lepton with $p_{{\textrm{T}}} >40$ GeV, or two opposite-charged leptons with $p_{{\textrm{T}}} >30(20)$ GeV, as well as hadronic $H_{{\textrm{T}}} >30$ GeV. A further set of pre-selection cuts were then applied. Before feeding events into the autoencoder as training data, 98% of the initial simulated background events were discarded by cutting on three high-level observables that encode the essential features of SUEP.

First, the multiplicity of charged tracks was required to be $N_{charged}\ge 70$. Second, the event ring isotropy variable introduced in [146], measuring the Wasserstein distance between a given event and a uniformly isotropic distribution of energy, was required to be ${\mathcal {I}}<0.07$. Finally, the inter-particle $\Delta R_{ij}$ distance averaged over all pairs of tracks in the event was required to be $\overline{\Delta R}>3$. Signal efficiency of these cuts varied from $1{-}30\%$ with $m_\textrm{D}$ and $T_\textrm{D}$.

A fully connected autoencoder with five layers was trained using QCD background events that passed the pre-selection cuts as training data. Each event was represented using a modified inter-particle distance matrix $\Delta \tilde{R}_{ij}$ of the 70 highest-$p_{{\textrm{T}}}$ charged tracks in the event.

$$\begin{aligned} \Delta \tilde{R}_{ij} \equiv \left\{ \begin{array}{ll} \Delta R_{ij}=\sqrt{(\Delta \eta _{ij})^{2} + (\Delta \phi _{ij})^{2}} &{} \quad i > j\\ p_{T,i}/{\textrm{GeV}} &{}\quad \text {for} \ i = j\\ 0 &{}\quad i < j. \end{array} \right. \end{aligned}$$

(31)

A modified mean-squared-error loss quantified the reconstruction error for each event.

5.4.5 Results

After training on background events, the autoencoder was fed test data including both background events and signal events across the range of simulated $m_\textrm{D}$ and $T_\textrm{D}$ points. Using the reconstruction loss as an anomaly score, ROC curves were constructed for each parameter point. To estimate the physical sensitivity of the model, the minimal excludable branching ratio of Higgs to SUEP for which $S/\sqrt{B + u_{sys}B^{2}}>2$ was computed. Statistical uncertainties due to the limited size of the simulated background sample became dominant as the cut threshold was increased before the classifier’s performance began to deteriorate, indicating that the sensitivity of a real search using this method could be even higher than we report here.

Fig. 54

Minimum excludable ${\textrm{Br}}(h \rightarrow {\textrm{SUEP}})$ at the HL-LHC, assuming $1\%$ systematic uncertainty on QCD background for fully connected autoencoder

As Fig. 54 illustrates, this autoencoder-based analysis could exclude Higgs branching ratios to SUEP down to 1% for dark meson mass $m_\textrm{D} <~ 1$ GeV and $T_\textrm{D}/m_\textrm{D} <~ 1$. If the dark shower temperature $T_\textrm{D}$ is $<~0.5 m_\textrm{D}$, branching ratios down to 5% could be probed for $m_\textrm{D}$ up to $\approx 8$ GeV.

Using a neural network architecture and event representation tailored to the essential characteristics of the SUEP signature, but without relying on the details of any signal simulation model, this study demonstrates that even the maximally challenging scenario of entirely prompt and hadronic SUEP can be probed at the HL-LHC.

5.5 Triggering on emerging jets

Contributors: Daniel Stolarski

The original Emerging Jets (EJ) theory paper [1] as well as the CMS search [88] (see also Sect. 2.1.2 of this white paper) considers a model with a colored mediator $\Phi$ that is pair produced at the LHC, which leads to a final state with two QCD jets and two EJs. Those works also consider mediator masses in the regime of $m_{\Phi } \gtrsim 600$ GeV. Given those assumptions, the vast majority of EJ events have substantial $H_{{\textrm{T}}}$ and thus the trigger efficiency is very high. In this section we consider relaxing both of the above assumptions and explore how one can still trigger on Emerging Jets.

In [228], an s-channel mediator was considered (see also Sect. 2.1.1 of this report), focusing on a $\textrm{Z}^{\prime }$ that couples to the quark current in the SM and the dark quark current. Such a mediator produces events that typically do not include additional hard jets. That work also considered the possibility of relatively light $\textrm{Z}^{\prime }$ down to masses of $\sim 50$ GeV. The typical $H_{{\textrm{T}}}$ of such events, particularly in the light $\textrm{Z}^{\prime }$ regime, is considerably lower than typical trigger thresholds at the LHC experiments, and other techniques are needed to increase the trigger efficiency.

The events were generated using a modified spin-1 mediator model¹² [63] implemented using the FeynRules  [219] package. The hard process is generated with Madgraph5_aMC@NLO [216] using a centre of mass energy of 13 TeV. This output is interfaced to the Hidden Valley [46, 47] module of pythia 8 [220], which simulates showering and hadronization in the dark sector as well as decays of dark hadrons to either other dark hadrons or to SM states. The $\textrm{Z}^{\prime }$ mass is varied and a $\textrm{Z}^{\prime }$ width of $\Gamma _{\textrm{Z}^{\prime }} = m_{\textrm{Z}^{\prime }}/100$ is used. The remaining dark sector parameters are varied across a few benchmark models shown in Table I of [228].

Initial state radiation (ISR) in QCD or EW is included at leading order in the hard processes. The resulting hadrons are clustered into jets using the Anti-$k_{t}$ algorithm [229] implemented in FASTJET [196] with a jet angular parameter $R = 0.4$ and a maximum pseudorapidity of $|\eta | < 2.49$ to be compatible with the ATLAS inner tracker. MLM matching and merging procedure [230] is employed for extra QCD radiation with XQCut of $m_{\textrm{Z}^{\prime }}$/10. A crude detector volume cut is implemented at the pythia 8 stage for which particles that are outside of a cylinder of $(r = 3000$ mm, $z = 3000$ mm) are considered stable.

Two main strategies are explored to increase the trigger efficiency. The first is exploiting the possibility of SM radiation from the initial state. While electroweak ($\textrm{W}/\text {Z}/\gamma$) radiation was explored, the most effective strategy was to use additional QCD radiation. This radiation can increase the trigger efficiency in two complimentary ways. First, the additional hard jet(s) can be used to trigger on directly. Second, the emerging jets tend to be boosted and carry more energy. This in turn will increase the $H_{{\textrm{T}}}$ () if the dark pion states are short (long) lived.

Fig. 55

(Figure 4 from [228]) Cross section times efficiency of various processes (leading order, 1-jet ISR, 2-jet ISR, electro weak ISR) scaled by their respective leading order process. The left plot uses the (here called MET) trigger and Model A which has a dark pion lifetime of 150 mm. The right uses the $H_{{\textrm{T}}}$ trigger and Model B which has a dark pion lifetime of 5 mm. The dotted line is the leading order process. On the right plot, the blue region on the left has zero of events that we simulated pass the trigger and thus an efficiency $\epsilon \lesssim \epsilon _{\textrm{min}} = 1/ 400{,}000$

Using ATLAS trigger thresholds from [231], we estimate the improvement in rate achieved by including radiation and the results are shown in Fig. 55. In addition to increasing the trigger efficiency, events with extra radiation have reduced rates, therefore Fig. 55 shows the ratio of the cross section times trigger efficiency for events with radiation to those without. We see that the largest improvement is additional radiation of two extra jets (green line). The left panel is a benchmark with a dark pion lifetime of 150 mm (Model A) and uses the missing energy trigger. We see that for a light $\textrm{Z}^{\prime }$, more than an order of magnitude improvement in rate is possible. The right side is a model with a dark pion lifetime of 5 mm (Model B) and uses the $H_{{\textrm{T}}}$ trigger. In that benchmark, the efficiency of the leading order process is below what was simulated for $m_{\textrm{Z}^{\prime }} \lesssim 350$ GeV, and the improvement is potentially even larger.

The first method considered above uses existing triggers, but [228] also considers implementing new triggers using modern machine learning techniques. As ISR is no longer relevant, pythia 8’s hidden valley production process $f\overline{f} \rightarrow \textrm{Z}^{\prime }$ processes is employed to generate events. Regardless of the lifetime of the dark pions, the detector subsystem with the largest number of decays is the inner tracker. Therefore, the strategy employed (which is also similar to that of [232] proposed for b-tagging), is to use the tracker information but not reconstruct tracks. Rather [228] proposes to use hit patterns in different layers of the tracker as an input to a support vector machine¹³ (SVM) from the TMVA toolkit [233]. A proper detector simulation of the inner tracker is outside of the scope, but a crude detector simulation with code used in [9] which encompasses the ATLAS tracker from the Inner B-layer (IBL) to the Transition Radiation Tracker (TRT). This detector simulation assumes simple models of energy loss through each thin layer of the detector.

When proposing new triggers, backgrounds must also be considered, and the main background for this strategy is $b\overline{b}$ jets as they have a very large rate and also produce displaced hadrons. Simulations are performed using $gg \rightarrow b\overline{b}$ with pythia 8’s heavy flavor hard $b\overline{b}$ processes. The inclusive background cross section is taken from the pythia 8. Pileup is added to both signal and background events with pythia 8’s minimum bias events. For each signal or background event, a number of minimum bias events are added randomly sampled from a poisson distribution with mean of $\mu = 50$, mimicking the Run 2 conditions.

Fig. 56

(Figure 8 from [228]) On the left, discrimination of signal (blue) from $b\overline{b}$ background (green) using a support vector machine. The flat bars (points) correspond to the training (test) set. On the right, Receiver Operation Characteristic ROC for four different lifetimes of signal. At a given background efficiency, the expected signal efficiencies increase as the dark pion lifetimes lower. The required background rejection is estimated to lie between the horizontal dotted lines. Both figures use a mediator mass of 500 GeV

The left panel of Fig. 56 shows the ability of the SVM to distinguish signal in blue from the dominant $b\overline{b}$ background. On the right panel we show the ROC curve as a function of lifetime. A background rejection of $\sim 10^{-2}{-}10^{-3}$ is needed for a novel high level trigger, and we see that efficiencies of ${\mathcal {O}}(10\%)$ are achievable, with larger efficiencies at lower dark pions lifetimes. It is also found that using an SVM trained on one signal benchmark can also give good acceptance for other signal benchmarks, showing great promise for such a new trigger.

6 Summary and perspectives

In this report, we have summarised the work performed in the context of the Dark Shower Snowmass project: it is the first comprehensive effort to gather the large, pre-existing theoretical, phenomenological and experimental communities working in this field, following initial discussions in the LHC Long-lived Particles Working Group [234] and also some presentations in the LHC Dark Matter Working Group. This report also concretely describes pathways for a systematic exploration of strongly interacting theories. In this context, we mainly concentrated on QCD-like scenarios leading to jetty signatures at the LHC, but we also discussed signatures such as SUEPs and glueballs which are typically associated with non-QCD like theories.

QCD-like scenarios, which are the main focus of this report, are inherently non-trivial to analyse due their non-perturbative nature. In such theories, confinement in the IR leads to bound states whose masses and interactions are governed by the UV dynamics. While the SM QCD has been analysed in great detail in terms of UV versus IR parametrizations, little is known for arbitrary gauge groups and flavor contents. Nevertheless, due to the interesting new signatures the strongly interacting scenarios could produce at the LHC, their phenomenology is being actively explored.

In this context, we began this report (see Sect. 2) with a review of the existing efforts and phenomenological parametrizations of QCD-like scenarios. We qualitatively illustrated the phenomenological differences obtained for various mediator mechanisms, giving rise to exotic LHC signatures such as emerging or semi-visible jets. We also discussed some existing experimental results constraining these models and ongoing efforts to search for these signatures.

If the dark sector is instead non QCD-like, other classes of spectacular signatures can be obtained in terms of SUEPs or glueballs. These were discussed Sect. 3, in which original phenomenological SUEP studies were presented, along with recent preliminary simulation tools for these scenarios.

After this overview of existing efforts and of the signature landscape, the report also addressed in Sect. 4 possible pathways for consistent theory frameworks, especially concentrating on semi-visible jets. In that section, lattice calculations, chiral perturbation theory and an analysis of symmetry breaking due to SM-DS portals were combined to exemplify avenues in theoretical model building. Improvements to the pythia 8 Hidden Valley module, made in the context of this Dark Shower Snowmass project, were also presented along with their validation. Combining the theory developments with the new Hidden Valley module, we then illustrated their impact on the phenomenology of semi-visible jets.

In the final section of the report, Sect. 5, we discussed some proposed improvements to LHC search strategies. These include efforts using machine learning, trigger considerations and the definition of new event level variables.

Strongly-interacting dark sectors are an exciting class of scenarios in which a vibrant community of theorists, phenomenologists and experimentalists is being invested. They could lead to spectacular signatures which have not yet been systematically explored at the LHC. In view of the large phenomenological interest of such theories, a more concentrated effort in theoretical work is needed, covering model building and classification of associated LHC signatures, a deeper understanding of hadronization physics, as well as studies of cross correlation with open problems of the SM such as the nature of dark matter. It is clear from this report that such a work involves communication among experts in SM QCD, lattice, and collider physics as well as in dark matter. We hope that our report lays down the foundations for such a wider exchange, and that this may help devising better strategies that could ultimately lead to a breakthrough in finding signals of strongly-interacting theories.