# Universal New Physics Latent Space

Anna Hallin  Institut für Experimentalphysik, Universität Hamburg, Luruper
Chaussee 149, 22761 Hamburg, Germany  Gregor Kasieczka  Institut für
Experimentalphysik, Universität Hamburg, Luruper Chaussee 149, 22761 Hamburg,
Germany  Sabine Kraml  Laboratoire de Physique Subatomique et de Cosmologie
(LPSC), Université Grenoble-Alpes, CNRS/IN2P3, F-38026 Grenoble, France  André
Lessa  Centro de Ciências Naturais e Humanas, Universidade Federal do ABC,
Santo André, 09210-580 SP, Brazil  Louis Moureaux  Institut für
Experimentalphysik, Universität Hamburg, Luruper Chaussee 149, 22761 Hamburg,
Germany  Tore von Schwartz  Institut für Experimentalphysik, Universität
Hamburg, Luruper Chaussee 149, 22761 Hamburg, Germany  David Shih  NHETC,
Department of Physics and Astronomy, Rutgers University, Piscataway, NJ 08854,
USA

###### Abstract

We develop a machine learning method for mapping data originating from both
Standard Model processes and various theories beyond the Standard Model into a
unified representation (latent) space while conserving information about the
relationship between the underlying theories. We apply our method to three
examples of new physics at the LHC of increasing complexity, showing that
models can be clustered according to their LHC phenomenology: different models
are mapped to distinct regions in latent space, while indistinguishable models
are mapped to the same region. This opens interesting new avenues on several
fronts, such as model discrimination, selection of representative benchmark
scenarios, and identifying gaps in the coverage of model space.

##  1 Introduction

With Run 3 of Large Hadron Collider (LHC) underway, the roadmap for the High-
Luminosity LHC established, and discussions of future colliders ongoing,
searches for new physics are at a crossroads. The simplest and/or most popular
extensions of the Standard Model (SM) are severely challenged by the null
results of previous LHC runs.111For comprehensive reviews of Run 2 searches,
see [1, 2, 3]. Looking beyond these models quickly opens up a vast theory
space of all kinds of possible Beyond the Standard Model (BSM) scenarios [4,
5]. This large theory space makes it difficult to decide which signatures
should be probed next with high priority. Optimizing analyses for specific
target models needs to be balanced against a broad coverage of signature space
and convenient re-interpretability in terms of alternative theoretical
scenarios [6].

The large theory space also makes mapping the implications of theoretical
models to the LHC data a non-trivial and computationally intensive effort. It
is unavoidable that different models will yield predictions that are very
similar within the current experimental capabilities. Conversely, a single
experimental result will have implications for the feasibility of multiple
models, not all of which are explicitly considered by the experiment. Modern
anomaly detection based search strategies [7, 8, 9, 10] complicate the matter
further by aiming for model-agnostic sensitivity to various outliers and
overdensities in the data. While this wider sensitivity might be useful in
accelerating the discovery of BSM physics, it makes the assessment of which
models were potentially covered by a given search more difficult.

In this work, we propose a novel approach using machine learning (ML) to
address the issues of re-interpretation and experimental coverage of possible
BSM signatures. Using a fully connected encoder architecture, we build an
abstract low-dimensional representation of BSM models (a so-called latent
space) so that experimentally indistinguishable models are mapped to the same
region in this space, while models which are far apart in latent space have
distinct phenomenologies. Once realized, such a universal new physics latent
space would simplify the re-casting of results by projecting into the latent
space as well as identifying insufficiently explored regions. Models with
similar signal properties would not have to be individually searched for but
could be replaced by a single representative scenario, provided appropriate
reinterpretation is made possible [11, 6] and that the choice of
representative scenarios does not leave significant holes in latent space. In
addition, since distances between two models in latent space give an
approximate measure of how experimentally distinct these models are, the
latent space can also be useful for tackling the inverse problem, i.e.
inferring from data which kind of BSM scenario is realised.

Several ML techniques have been suggested for the task of determining what BSM
parameters are allowed within the constraints of the data. For a review of
these, see [12]. More specifically related to our work, the connection between
the latent space and the underlying physics has previously been studied in the
context of strings [13], $W$-jets [14], supersymmetry (SUSY) [15] and QCD vs.
BSM jets [16] where the aims and strategies of the latter are the most similar
to ours. But while [16] deals with single jets and constructs a loss function
that aims to preserve the distance between two events (defined via metrics on
the data space and the latent space), our model operates on event level rather
than jet level and uses a contrastive loss.

The rest of this paper is organised as follows: section 2 will introduce the
machine learning approach; sections 3–5 then apply this technique to three
different datasets with increasing complexity: simulated gluino decays from a
gluino-neutralino simplified model, simplified Dark Matter (DM) models leading
to one hard jet and missing energy, and the public Dark Machines Anomaly
Challenge dataset; section 6 offers a summary and conclusions. Additional
plots for the three datasets are collected in the Appendix.

##  2 Method

In this work, we learn a mapping of the high dimensional theory space of BSM
theories to a lower dimensional latent space using an encoder network trained
on physical (detector-level) observables such as missing transverse energy
(MET) and jet variables.

Since a set of events better represents the signal features than a single
event, a single input data point ($\vec{X}_{i}$) for the network corresponds
to the concatenation of observables (or input features) for a set of 10 Monte
Carlo (MC) events generated by the same BSM model222By BSM model we mean a
specific BSM scenario with fixed model parameters. For instance, for the
simplified gluino-neutralino model in the next section, two distinct sets of
gluino and neutralino masses are considered as two different BSM models.,
e.g.:

$\vec{X}_{i}=\left(\left\\{{\rm MET},p_{T}(j_{1}),...\right\\}_{\rm
event1},\;...,\;\left\\{{\rm MET},p_{T}(j_{1}),...\right\\}_{\rm
event10}\right)$ (1)

The events for each set were chosen at random, and the sets were recreated for
each epoch.

The encoder network architecture is a fully-connected MLP with hidden layer
dimensions (384, 512, 256, 128, …, 4) and ReLU activations. The final hidden
layer outputs to the latent space. We set the latent space dimension to two
for ease of visualization. Increasing the latent space to higher dimension was
not found to provide any further benefit in the different scenarios studied in
this work.

We use a contrastive loss function which is designed to learn the essential
differences between the various models. Contrastive learning [17] relies on
determining whether two data points are similar or not. Based on this
similarity measure, the data points either get separated or pulled together in
latent space. The resulting mapping should encode the pairwise relationships
between the data points.

Concretely, each pair of event sets $\vec{X}_{i}$ and $\vec{X}_{j}$ is
assigned a label $Y_{ij}$ that is set to zero for points from the same model
(similar points) or to 1 for points from distinct models (dissimilar points).
In other applications, like image recognition [18] or natural language
processing [19], similar samples are often generated using data augmentation
[20], which makes it possible to use this method even without labels. In this
work, we instead use the underlying model to define the similarity. The
contrastive loss function has the form

$\mathcal{L}(\theta)=\sum_{i\neq
j}L(\theta,(Y_{ij},\vec{X}_{i},\vec{X}_{j}))\,,$ (2)

where

$L(\theta,(Y_{ij},\vec{X}_{i},\vec{X}_{j}))=(1-Y_{ij})\frac{1}{2}D^{2}_{\theta}+Y_{ij}\frac{1}{2}\left(\max\left(0,d-D_{\theta}\right)\right)^{2}\,.$
(3)

Here, $D_{\theta}\equiv D_{\theta}(\vec{X}_{i},\vec{X}_{j})$ is the Euclidian
distance in the latent space, $\theta$ are the weights of the encoder network
being trained and $d$ is an additional hyperparameter that is used to deal
with unbounded representation spaces. Note that for bounded representation
spaces it is not necessary to introduce a maximum distance. In our case,
however, we do not impose any bounds on the latent space and $d$ is needed to
ensure a converging training. In what follows, we set $d=1$.

With the above definition, the loss function is minimized if $D_{\theta}^{2}$
is small for similar points and large for dissimilar points. If different
models occupy disjoint (i.e. non-overlapping) regions of phase space, the loss
can be trivially set to zero by collapsing each model to a distinct point in
latent space, separated from all the other points by distances $D_{\theta}\geq
1$. The interesting, non-trivial action of the loss arises when different
models occupy overlapping regions of phase space: when inputs from these
regions of phase space are fed to the network, the loss tries to minimize
$D_{\theta}$ when the inputs come from the same model, but tries to maximize
$D_{\theta}$ when the inputs come from different models. This tension results
in a nonzero value of $D_{\theta}$, i.e. the models cannot be collapsed to a
single point by the network. Instead, they are mapped to non-trivial regions
in the latent space which characterize and represent the phase space regions
occupied by each model.

Three datasets, described in more detail below, were tested with this
contrastive learning model. They all feature jets plus MET signatures but
represent different levels of complexity for model discrimination. The input
data points $\vec{X}_{i}$ for the training (see Eq. (1)) were thus made of the
MET plus various kinematic variables of the jets ($p_{T}$, $\eta$,
$\Delta\phi$, etc.); the specific inputs for each dataset will be described
below. Every input feature was normalized to have mean zero and standard
deviation one. The training was done using Adam [21] as optimizer with a
learning rate of $10^{-5}$. Once trained, the model state with the lowest
validation loss was selected.

##  3 MSSM gluino simplified model

###  3.1 Dataset

The first dataset to be studied is a supersymmetric scenario, with pair-
production of gluinos in the Minimal Supersymmetric Standard Model (MSSM).
Each gluino is assumed to decay with 100% branching ratio to a neutralino
$\tilde{\chi}^{0}_{1}$ and a quark-antiquark pair, $\tilde{g}\to
q\bar{q}\tilde{\chi}^{0}_{1}$, as illustrated in Fig. 1. This leads to events
with four hard jets and large MET. The dataset was generated with
MadGraph5_aMC@NLO v3.3.0 [22], Pythia 8.306 [23, 24], Delphes 3 [25], and
NNPDF 2.3 [26] at LO. Nine distinct BSM models were simulated, with gluino
masses of $m_{\tilde{g}}=1.1$, 1.6 and 2.1 TeV and neutralino masses of
$m_{\tilde{\chi}^{0}_{1}}=0.1$, 0.5 and 0.9 TeV. The events from all mass
configurations were combined into one single dataset.

Figure 1: Process simulated for the MSSM gluino dataset.

The features used in training from this dataset are the $p_{T},\eta,m$ and
$\Delta\phi$ of the first four jets in the event, where $p_{T}$ is the
transverse momentum, $\eta$ the pseudorapidity, and $m$ the invariant mass of
the jet; $\Delta\phi$ is the difference in azimuthal angle between the jet and
the MET, defined by

$\Delta\phi=\min\left(\left|\phi_{j}-\phi_{\rm
MET}\right|,2\pi-\left|\phi_{j}-\phi_{\rm MET}\right|\right).$ (4)

In the set of features we also include the MET and the invariant mass of all
possible jet pairs:

$m^{2}_{j_{a},j_{b}}=2\,p_{T_{a}}\,p_{T_{b}}\left[\cosh(\eta_{a}-\eta_{b})-\cos(\phi_{a}-\phi_{b})\right].$
(5)

The distributions of the leading jet $p_{T}$, the MET and $m^{2}_{1,2}$ are
shown in Fig. 2, while the full set of features for all four jets can be found
in Appendix A.1. Since the amount of energy going into jets depends mostly on
the gluino-neutralino mass difference ($\Delta m$), we expect harder jets and
a harder MET spectrum for larger $\Delta m$. This is confirmed in Fig. 2,
where the hardest MET spectrum corresponds to $\Delta m=2$ TeV and the softest
to $\Delta m=0.2$ TeV. A similar behavior can be observed for the $p_{T}$ and
the invariant mass $m_{1,2}^{2}$. Furthermore, we see that models with similar
values of $\Delta m$ have similar features and are harder to distinguish
experimentally.333Note that we are not making use of the total rate of events,
which depends on the gluino cross-section and could in principle be used to
distinguish between models with different gluino masses.

Figure 2: Distributions of the leading jet $p_{T}$, plus MET and
$m^{2}_{1,2}$, for each of the mass configurations in the MSSM gluino dataset.

###  3.2 Results

Training the encoder with data from the nine models and using the signal
features discussed in the previous section, we obtain the latent space mapping
shown in Fig. 3. The plot was created by first mapping the data to the
embedding space and then fitting a probability density in this space using a
kernel density estimation with a Gaussian kernel. The contours in the plot
represent equal-density contours that enclose 50% of the probability density,
corresponding to a cumulative distribution function (CDF) value of 0.5. This
means that half of the events fall within these contours. The colors of the
different contours correspond to the mass difference $\Delta m$ between the
gluino and the neutralino. We see that the network has not only managed to
separate the different datasets into different clusters, but it has also
ordered the latent space according to $\Delta m$ along the diagonal axis. This
behavior agrees with the strong correlation that the signal features have with
$\Delta m$, as shown in Fig. 2. Furthermore, somewhat orthogonal to this main
diagonal, the network seems to sort the models according to $m_{\tilde{g}}$,
with smaller gluino masses corresponding to larger values for the $x$-axis.
This separation is not trivial, given that the signal features are weakly
dependent on the gluino mass.

Figure 3: Resulting latent space for the MSSM gluino dataset, showing contours
corresponding to a CDF value of 0.5 for each model. The different contours are
colored based on the mass difference $\Delta
m=m_{\tilde{g}}-m_{\tilde{\chi}^{0}_{1}}$, which is also given as the contour
label (in TeV); dark red: $\Delta m=2$ TeV; red: $\Delta m=1.6$ (dashed) and
$1.5$ (solid) TeV; green: $\Delta m=1.2$ (dash-dot), $1.1$ (dashed) and $1.0$
(solid) TeV; cyan: $\Delta m=0.7$ (dash-dot) and $0.6$ (dashed) TeV; blue:
$\Delta m=0.2$ TeV. The linestyles are chosen based on the gluino mass: solid
for $m_{\tilde{g}}=2.1$ TeV, dashed for $m_{\tilde{g}}=1.6$ TeV, and dotted
for $m_{\tilde{g}}=1.1$ TeV.

##  4 Dark Matter simplified model

###  4.1 Dataset

In the gluino-neutralino scenario discussed above, the theory space of BSM
models was two dimensional ($m_{\tilde{g}},m_{\tilde{\chi}_{1}^{0}}$) and was
successfully mapped to a two dimensional latent space. As the next step in
complexity, we analyze distinct classes of two dimensional models, which
corresponds to adding a third (discrete) dimension to the theory space. To
this end, we consider DM production at the LHC, which is typically searched
for in events with a single hard jet plus missing energy (additional soft jets
are also present from initial and final state radiation, ISR and FSR). The
following three models are used:

  * •

Vector mediator: in this model, the DM candidate is pair-produced through an
s-channel spin-1 mediator, as shown in Fig. 4a.

  * •

Pseudoscalar mediator: this model is similar to the above, but with the
production mediated by a pseudoscalar, which is produced through the gluon
fusion diagram shown in Fig. 4b.

  * •

Squark mediator: unlike the previous cases, this model presents a DM-squark
associated production, where the squark further decays to DM and a hard jet,
as shown in Fig. 4c.

In all these models, the DM candidate is a Majorana fermion and a hard jet is
present either from ISR in the vector and pseudoscalar mediator models, or
from the squark decay. We also point out that in the three scenarios the DM
production is dominated by distinct initial-state partons: the vector mediator
is produced from $q\bar{q}$ collision, the pseudoscalar from $gg$ collision,
and the squark-DM production requires a $qg$ initial state. Monte Carlo events
were generated for the following ranges of DM and mediator masses:

$\displaystyle m_{\rm DM}$ $\displaystyle=100-900\mbox{ GeV in steps of 100
GeV}\,,$ $\displaystyle M_{\rm med}$ $\displaystyle=600-2000\mbox{ GeV in
steps of 200 GeV}\,,$

where $M_{\rm med}=M_{Z^{\prime}},M_{\phi},M_{\tilde{q}}$ for the vector,
pseudoscalar and squark models, respectively and $m_{\rm DM}\leq 2M_{\rm med}$
was imposed. For each mass combination, 50k MC events were generated at LO
using MadGraph5_aMC@NLO v.3.4.1[27, 28], Pythia 8 v8.307 [23] and Delphes
v3.5.0 [25] with the PDF set NNPDF31_nnlo_as_0118 and the default Delphes card
for the CMS detector. In the case of the vector and pseudoscalar mediator
models, the DM production was computed with one additional extra jet and the
MLM matching scheme was used with xqcut $=M_{Z^{\prime}}/15$ or $M_{\phi}/12$
in order to properly model ISR jets.

$\bar{q}$$q$$Z^{\prime}$ (a)

$g$$g$$\phi$$t$$t$$t$ (b)

$q$$g$$q$$\tilde{q}$$\tilde{q}$ (c)

Figure 4: Examples of diagrams for DM production in the: a) vector mediator,
b) pseudoscalar mediator and c) squark mediator models.

The $p_{T}$, $\eta$, $\Delta\phi$ and mass of the leading two jets and the MET
of the event are used as input features, with $\Delta\phi$ as defined in Eq.
(4). In order to satisfy minimum trigger requirements, a cut on the transverse
momentum is applied by demanding $p_{T,1}>100$ GeV and $p_{T,2}>50$ GeV. Since
MET is expected to be the leading feature for these models, we show in Fig. 5
the MET distribution for the highest and lowest mediator masses for the
scenario where $m_{\rm DM}=100$ GeV. The full feature distributions of the
three mediator scenarios are shown in Appendix A.2 for the different mediator
masses. We see that the vector and pseudoscalar models have a MET distribution
which shifts slightly with the increase of the mediator mass. On the other
hand, the squark model presents a distribution which is more strongly
correlated to the mediator mass. This is expected since the $p_{T}$ of the
hard jet produced by the squark decay directly depends on the squark-
neutralino mass difference. It is also important to notice that, since the
scalar and vector mediators are produced on-shell, the signal features have a
weak dependence on $m_{\rm DM}$.

Figure 5: MET distribution of the three different mediator types with $m_{\rm
DM}=100$ GeV and the highest and lowest of the mediator masses. It is clearly
visible that while the peak of the distribution stays in the same place for
the vector and pseudoscalar mediators, for the squark mediator it gets shifted
with higher mediator mass.

###  4.2 Results

To get an idea of the event characteristics due to different types of
mediators, first separate trainings were performed for the datasets of the
individual scenarios. The arrangement of the mass configurations was then
investigated with the goal of tracing the latent space structure back to
physical properties. The contours corresponding to a CDF value of 0.5, created
with the method described in section 3.2, are shown in Fig. 6. To guide the
eye, the coloured contours show the distributions for $m_{\rm med}=600$–2000
GeV (in steps of 200 GeV) with the dark matter mass fixed to $m_{\rm DM}=100$
GeV; the contour lines for all other mass configurations are drawn in light
gray.

As can be seen, for the vector mediator data, the distributions of the
individual mass configurations are heavily overlapping. With increasing
mediator mass, the distributions get somewhat shifted and narrowed in
particular in latent dimension 1, but they remain overlapping even for very
different mediator masses, e.g. $m_{\rm med}=600$ GeV compared to $m_{\rm
med}=2000$ GeV. A similar behaviour is noticeable for the pseudoscalar
mediator data. All in all, for both these scenarios, different mass
configurations lead to overlapping representations. This is somewhat expected,
since the event features are strongly dependent on the ISR jets, which have a
weak dependence on the mediator mass, as illustrated by Fig. 5. A different
behaviour is observed for the squark mediator dataset, where the hard jet
spectrum is directly related to the squark-neutralino mass difference,
resulting in a clear (albeit still partly overlapping) separation between
different mediator masses in the latent space.

Figure 6: Latent space for the separate embeddings of the different mediator
types. The contours correspond to a CDF value of 0.5 for varying mediator
masses, with $m_{\rm DM}=100$ GeV. The contour lines for all other mass
configurations are drawn in light gray. Figure 7: Latent space for the
embeddings of pairs of different mediator types (left and center), and for all
mediator types (right). The contours correspond to a CDF value of 0.5 for the
varying mediator masses, with $m_{\rm DM}=100$ GeV. The contour lines for all
other mass configurations are drawn in light gray.

Having seen the internal structure of the embeddings for the theories with
different mediator particles, we next investigate how the distinct scenarios
behave when embedded together in the same latent space. To this end, a network
was trained on the combined data (including all mass configurations) from
pairs of mediator types. The resulting latent spaces for a selection of mass
parameters are shown in Fig. 7, left and center panels. For the combination of
vector and pseudoscalar mediator data, no clustering with respect to the type
of mediator is visible. Instead, distributions for the same mass parameters
for the different mediators are approximately mapped to the same region. As
seen in Fig. 5, in order to distinguish these models it would be required to
select events at the tails of the MET distribution, which corresponds to a
harder MET cut. Although possible, this cut would introduce a model bias which
we try to avoid here.

Once we combine the pseudoscalar and squark models, the mapping to the latent
space displays a more interesting arrangement: the network learns to embed
theories with a different mediator particle in different latent space regions
while maintaining their respective internal structure from Fig. 6. The fact
that this internal arrangement of the theories remains unchanged demonstrates
the stability of the embedding. Furthermore, we see that the squark model with
$m_{\rm med}=600$ GeV has the largest overlap with the pseudoscalar models, as
is expected from Fig. 5.

Finally, a training on the data from all theories was performed. The resulting
latent space is shown in the right panel of Fig. 7. The arrangement of the
theories in the latent space could be expected from the previous results.
Again the vector and pseudoscalar data are mapped to the same region, close to
the distributions of low-mass squark mediators. All the theories remain in the
same similarity relation as observed in the pairwise embeddings. It is also
interesting to notice that all the models are mostly distributed along a
vertical line, while each model presents some spread around the horizontal
dimension. This suggests that in this case a one dimensional latent space
might be sufficient to capture the main distinguishing features of the theory
space considered. We also note that using a single vector (or pseudoscalar)
benchmark model in addition to squark models with several $\Delta m$ values is
sufficient to provide a good coverage of the latent space.

Figure 8: Binned latent space after training on a combined dataset containing
all mediator types and masses. The bins are colored according to the average
value of the displayed input feature.

To provide a better understanding of the latent space structure, Fig. 8 shows
the correlation of the latent space position with the MET of the event, the
$p_{T}$ and $\Delta\phi$ of the first jet, and the $\Delta\phi$ of the second
jet (see Appendix A.2 for the same plots for the full selection of features).
Here, a binning is applied to the latent space from Fig. 7 (right) resulting
from the training on all mediator types. For each bin, we compute the mean of
the kinematic features over all data points ($\vec{X}_{i}$). However, since
each data point carries information about a set of 10 events, the mean over
the set is used to represent each data point. We can see that the input
features $p_{T,1}$ and MET both define the latent space structure along the
vertical latent dimension. The fact that they have the same orientation is
expected, since MET is strongly correlated with $p_{T,1}$. The other
(horizontal) latent dimension shows a strong correlation with
$\Delta\phi(j_{2})$. We note, however, that this feature is not very helpful
when distinguishing between different models, as seen in Fig. 7 (right). Plots
like these allow a better understanding of the learned embedding and make the
structure learned by the network easier to interpret.

##  5 Dark Machines Anomaly Challenge

###  5.1 Dataset

As a final example, we consider a subset of the Dark Machines Anomaly
Challenge [29, 8] models. The Dark Machines data consist of simulated events
for 12 different BSM models as well as a dataset with 26 different SM
background processes. From these we have selected the subset of BSM models
giving jets+MET signatures, shown in Table 1. These include the gluino-
neutralino scenario discussed in section 3, the vector mediator model from
section 4 and two additional SUSY simplified models corresponding to stop and
squark pair production with direct decays to neutralinos. We have also
included the SM background dataset, which is dominated by QCD (multijets),
$Z$+jets, $W$+jets and $t\bar{t}$.

The Dark Machines data is divided into four channels, corresponding to
different signal regions with different trigger requirements. In our analysis,
we focus only on channel 1, which requires:

  1. 1.

at least 4 (b-)jets with $p_{T}>50$ GeV,

  2. 2.

one (b-)jet with $p_{T}>200$ GeV,

  3. 3.

$H_{T}\geq 600$ GeV,

  4. 4.

MET$\geq 200$ GeV and ${\rm MET}/H_{T}\geq 0.2$,

where $H_{T}$ is the scalar sum of the transverse momenta of all (b-)jets. The
input features used for this analysis are the kinematic variables $E$,
$p_{T}$, $\eta$ and $\phi$ of the first four jets, and the MET of the event.
The distributions for all the features for the selected models are shown in
Appendix A.3.

BSM scenario |  Physical process and model parameters
---|---
DM Vector Mediator |  $pp\to Z^{\prime}\to\chi\chi$• $m_{Z^{\prime}}=2$ TeV, $m_{\rm DM}=50$ GeV
Gluino Simplified Models |  $pp\to\tilde{g}\tilde{g}$, $\tilde{g}\to qq+\tilde{\chi}_{1}^{0}$ • $m_{\tilde{g}}=1.4$ TeV, $m_{\chi^{0}}=1.1$ TeV • $m_{\tilde{g}}=1.6$ TeV, $m_{\chi^{0}}=0.8$ TeV
Stop Simplified Model |  $pp\to\tilde{t}\tilde{t}$, $\tilde{t}\to t+\tilde{\chi}_{1}^{0}$ • $m_{\tilde{t}}=1$ TeV, $m_{\chi^{0}}=0.3$ TeV
Squark Simplified Model |  $pp\to\tilde{q}\tilde{q}$, $\tilde{q}\to q+\tilde{\chi}_{1}^{0}$ • $m_{\tilde{q}}=1.8$ TeV, $m_{\chi^{0}}=0.8$ TeV
Table 1: Dark Machines BSM models considered in our analysis. In [29, 8] they
are called, from top to bottom: “$Z^{\prime}+$monojet”, “SUSY gluino-gluino
production”, “SUSY stop-stop production”, and “SUSY squark-squark production”;
we have renamed them to make the correspondence with the models in sections 3
and 4 apparent.

In Fig. 9, we show the MET and the $p_{T}$ distributions for the first 3 jets
for the selected BSM models and the SM background. The squark model has the
largest mass difference ($\Delta m=1$ TeV) and each squark decays to a single
jet plus a neutralino. As a result, the MET distribution for this model is the
hardest of all the models considered, as seen in the first panel of Fig. 9.
All the other models have similar MET spectra, except for the light gluino
scenario, which has a small mass gap ($\Delta m=0.3$ TeV) and displays a
slightly softer distribution. A similar behavior is also seen for the first
jet $p_{T}$. Note that, despite having a large mass difference, all the jets
in the vector mediator model come from initial state radiation and are not
dependent on $\Delta m$. Once we consider the third jet $p_{T}$ we see that
the models start to display a different behavior. In particular, the heavy
gluino model shows a peak at higher $p_{T}$ values when compared to all the
other models. This is expected, since each produced gluino decays to two jets
and a neutralino, thus resulting in events with four hard jets. The same would
also be true for the light gluino model, except for its small $\Delta m$,
which leads to a softer spectrum. The other models, however, tend to display
softer third and fourth jets as they originate from ISR. As a result, the
stop, $Z^{\prime}$ and light gluino models all display a similar $p_{T,3}$
spectrum. Finally we point out that while the SM background has very distinct
(softer) MET and $p_{T,1}$ distributions, the transverse momentum
distributions of the second and third jets are similar to those of the light
gluino and the $Z^{\prime}$ models.

Figure 9: MET and $p_{T}$ distributions for the first 3 jets. The histograms
for the five BSM models from Table 1 are shown in color, while the SM
background corresponds to the filled gray histogram.

###  5.2 Results

The latent space resulting from training over all the five BSM models plus the
SM background is shown in Fig. 10. The left column shows the latent space
contour plot (created with the same method as described in section 3.2) for
the different signals, overlaid on the binned feature distributions of MET and
$p_{T,3}$ in latent space, while the right column shows the distributions of
these two quantities in physical space. These two features are highlighted
here because they seem to be the main sorting axes of the features in latent
space. The full binned feature distribution in latent space for all features
are shown in Appendix A.3.

Figure 10: Left: contours of the different models in the resulting latent
space overlaid on the mean of the MET and the $p_{T}$ of the third jet of the
events placed in each particular bin in latent space. The contours correspond
to a CDF value of 0.5 for each model. Right: the distributions of the MET and
the $p_{T}$ of the third jet for the different models.

Analyzing the contour plots, we see that the light gluino ($m_{\tilde{g}}=1.4$
TeV), the $Z^{\prime}$ and the stop models all fall into a similar region of
latent space. This is expected from the discussion in the previous section,
where it was shown that the MET and $p_{T}$ distributions for these models are
similar. In addition, the light gluino model displays a softer MET
distribution, which puts this model at slightly higher values for the second
(vertical) latent space dimension. The background has the softest MET
distribution, thus being located at the very top of the latent space plot. The
squark and the heavy gluino models, on the other hand, populate disjoint
regions. While the squark model is located at the bottom of the plot, which
corresponds to the highest MET values, the heavy gluino model occupies the
center-rightmost region. This region corresponds to intermediate MET values
and large $p_{T,3}$ values, as expected from the previous discussion.

In summary, the network was able to explore the signal features of each model
and cluster them accordingly. We also point out that the $Z^{\prime}$, light
gluino and stop models have similar features, which means that for
benchmarking purposes just one of them is sufficient to represent this region
in latent space. In turn, an additional benchmark model with a similar MET but
a harder $p_{T}(j_{3})$ spectrum would result in a better coverage of the
signature space.

##  6 Conclusions

Motivated by the large phase space of BSM models, we proposed a strategy to
build a joint embedding space of theories. The main property of this abstract
space is that models are clustered according to their phenomenology (i.e.,
their experimental signatures), so different models are mapped to different
regions in this space while indistinguishable models are mapped to the same
region.

We demonstrated the feasibility of this technique in three examples of
increasing complexity. First, we considered a simplified MSSM scenario with
different gluino and neutralino masses and observed that the embedding without
additional guidance correctly identifies the gluino to neutralino mass
difference as the phenomenologically most relevant parameter and arranges the
models by this quantity in latent space accordingly.

Next, we investigated three different models of dark matter production in
final states with a hard jet and missing energy: a vector mediator, a
pseudoscalar mediator, and a squark mediator model. We observed that the
learned representations for the vector mediator and pseudoscalar mediator do
not vary significantly with the mediator mass, mostly because the jets of
these events come from ISR. For the squark mediator, on the other hand, where
the mass of the jet is directly related to the mass difference between the
squark and neutralino, different masses are embedded in different regions of
latent space. When learning a joint embedding of all three models, the same
picture remains: vector and pseudoscalar mediator are mapped to one region,
which is phenomenologically indistinguishable from a low mediator mass squark
model. Higher mediator mass squark models extend away from this point, ordered
by their mass.

The third showcase was taken from the Dark Machines Anomaly Challenge dataset.
Here, we observed a sorting of the latent space along the direction of missing
transverse energy and the $p_{T}$ of the third jet. The network’s focus on MET
was to be expected, given how distinguished the MET distributions of the
background and the squark model are compared to the remaining models. In
addition, the network was able to separate models with similar MET spectra by
exploiting the $p_{T}$ of the third jet, which puts the heavy gluino model far
from the stop and DM vector mediator models. Therefore, the latent space can
not only cluster together models with similar features, but also indicate
possible discriminating observables.

The three examples investigated here have illustrated that the proposed
approach can achieve the clustering of models with similar features, thus
providing a way to select a minimal set of benchmark models with a maximal
coverage of the feature space. In addition, the correlation between the latent
space dimensions and the physical features can be used to identify suitable
observables for experimentally distinguishing BSM models. Finally, regions in
latent space not covered by any of the chosen models can indicate gaps in
signature space, which could be filled by different benchmark choices, thus
maximizing the physical signatures explored at the LHC and future colliders.

Together, these results demonstrate the potential usefulness of a learned
representation to combine and understand data from different theoretical
models. However, to fully realize the potential of this approach, more work is
needed. Most important is an inclusion of the cross-section of a given process
in addition to the distributions of features, as this is an important property
that would be useful for separating different models. Next a more expressive
basis of features and higher capacity learning algorithms need to be explored.
It will also be useful to develop an algorithm that is capable of sampling the
latent space and reconstruct the physical features corresponding to these
sampled points. This would allow us to explore unpopulated regions of latent
space, which could potentially represent signal features which are not present
in any of the BSM models considered. Together this would motivate the building
of a universal embedding of BSM theories to guide the exploration of model
space at the LHC and beyond.

## Acknowledgements

This research was supported in part through the Maxwell computational
resources operated at Deutsches Elektronen-Synchrotron DESY, Hamburg, Germany.
The work of AH, GK, LM and TS was supported by the Deutsche
Forschungsgemeinschaft under Germany’s Excellence Strategy – EXC 2121 Quantum
Universe – 390833306, and under PUNCH4NFDI – project number 460248186. The
work of SK was supported in part by the IN2P3 theory master project
“DataMATTER”. AL is supported by FAPESP grants no. 2018/25225-9 and
2021/01089-1. DS is supported by the U.S. Department of Energy (DOE), Office
of Science grant DOE-SC0010008

We thank the organizers of the 2022 MITP workshop A Deep-Learning Era for
Particle Theory, the Les Houches PhysTeV 2023 BSM workshop, and the Aspen
Center for Physics (supported by National Science Foundation grant
PHY-2210452) for opportunities to discuss this work.

##  Appendix A Supplementary material

###  A.1 MSSM gluino dataset

Figure 11 shows the distributions of the features used in the training on the
MSSM gluino dataset.

Figure 11: Feature distributions of the MSSM gluino dataset.

###  A.2 Dark Matter dataset

Figures 12, 13 and 14 show the feature distributions for the three different
mediator scenarios, for the different mediator masses and $m_{\rm DM}=100$
GeV. Figure 15 shows the binned latent space for the different features after
training on a combined dataset of all mediator types of all mass combinations.

Figure 12: Features of the vector mediator scenario for the different mediator
masses with $m_{\rm DM}=100$ GeV. Figure 13: Features of the pseudoscalar
mediator scenario for the different mediator masses with $m_{\rm DM}=100$ GeV.
Figure 14: Features of the squark mediator scenario for the different mediator
masses with $m_{\rm DM}=100$ GeV. Figure 15: Binned latent space after
training on a combined dataset containing all mediator types and all mass
configurations in the Dark Matter dataset.

###  A.3 Dark Machines dataset

Figure 16 shows the distributions of the different features used for training
in the Dark Machines dataset. Figure 17 shows the binned feature distributions
in the latent space of the Dark Machines study.

Figure 16: Feature distributions for the selected signals from the Dark
Machines dataset.  Figure 17: Binned feature distributions in the latent space
of the Dark Machines study.

## References

  * [1] ATLAS Collaboration, G. Aad et al., “The quest to discover supersymmetry at the ATLAS experiment,” arXiv:2403.02455 [hep-ex].
  * [2] CMS Collaboration, A. Hayrapetyan et al., “Dark sector searches with the CMS experiment,” arXiv:2405.13778 [hep-ex].
  * [3] CMS Collaboration, A. Hayrapetyan et al., “Review of searches for vector-like quarks, vector-like leptons, and heavy neutral leptons in proton-proton collisions at $\sqrt{s}$ = 13 TeV at the CMS experiment,” arXiv:2405.17605 [hep-ex].
  * [4] T. Bose et al., “Report of the Topical Group on Physics Beyond the Standard Model at Energy Frontier for Snowmass 2021,” arXiv:2209.13128 [hep-ph].
  * [5] P. J. Fox et al., “TF08 Snowmass Report: BSM Model Building,” arXiv:2210.03075 [hep-ph].
  * [6] S. Bailey et al., “Data and Analysis Preservation, Recasting, and Reinterpretation,” arXiv:2203.10057 [hep-ph]. Snowmass 2021 white paper.
  * [7] G. Kasieczka, B. Nachman, D. Shih, et al., “The LHC Olympics 2020 a community challenge for anomaly detection in high energy physics,” Rept. Prog. Phys. 84 no. 12, (2021) 124201, arXiv:2101.08320 [hep-ph].
  * [8] T. Aarrestad et al., “The Dark Machines Anomaly Score Challenge: Benchmark Data and Model Independent Event Classification for the Large Hadron Collider,” SciPost Phys. 12 no. 1, (2022) 043, arXiv:2105.14027 [hep-ph].
  * [9] G. Karagiorgi, G. Kasieczka, S. Kravitz, B. Nachman, and D. Shih, “Machine learning in the search for new fundamental physics,” Nature Rev. Phys. 4 no. 6, (2022) 399–412.
  * [10] V. Belis, P. Odagiu, and T. K. Aarrestad, “Machine learning for anomaly detection in particle physics,” Rev. Phys. 12 (2024) 100091, arXiv:2312.14190 [physics.data-an].
  * [11] LHC Reinterpretation Forum Collaboration, W. Abdallah et al., “Reinterpretation of LHC Results for New Physics: Status and Recommendations after Run 2,” SciPost Phys. 9 no. 2, (2020) 022, arXiv:2003.07868 [hep-ph].
  * [12] R. Baruah, S. Mondal, S. K. Patra, and S. Roy, “Probing intractable beyond-standard-model parameter spaces armed with Machine Learning,” arXiv:2404.02698 [hep-ph].
  * [13] A. Mütter, E. Parr, and P. K. S. Vaudrevange, “Deep learning in the heterotic orbifold landscape,” Nucl. Phys. B 940 (2019) 113–129, arXiv:1811.05993 [hep-th].
  * [14] J. H. Collins, “An Exploration of Learnt Representations of W Jets,”  9, 2021.  arXiv:2109.10919 [hep-ph].
  * [15] J. Baretz, N. Carrara, J. Hollingsworth, and D. Whiteson, “Visualization and efficient generation of constrained high-dimensional theoretical parameter spaces,” JHEP 11 (2023) 062, arXiv:2305.12225 [hep-ph].
  * [16] S. E. Park, P. Harris, and B. Ostdiek, “Neural embedding: learning the embedding of the manifold of physics data,” JHEP 07 (2023) 108, arXiv:2208.05484 [hep-ph].
  * [17] R. Hadsell, S. Chopra, and Y. LeCun, “Dimensionality reduction by learning an invariant mapping,” in 2006 IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR’06), vol. 2, pp. 1735–1742.  2006\.
  * [18] T. Chen, S. Kornblith, M. Norouzi, and G. Hinton, “A Simple Framework for Contrastive Learning of Visual Representations,” arXiv:2002.05709 [cs.LG].
  * [19] Q. Cheng, X. Yang, T. Sun, L. Li, and X. Qiu, “Improving Contrastive Learning of Sentence Embeddings from AI Feedback,” arXiv:2305.01918 [cs.CL].
  * [20] P. Khosla, P. Teterwak, C. Wang, A. Sarna, Y. Tian, P. Isola, A. Maschinot, C. Liu, and D. Krishnan, “Supervised contrastive learning,” arXiv:2004.11362 [cs.LG].
  * [21] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” arXiv:1412.6980 [cs.LG].
  * [22] J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, O. Mattelaer, H. S. Shao, T. Stelzer, P. Torrielli, and M. Zaro, “The automated computation of tree-level and next-to-leading order differential cross sections, and their matching to parton shower simulations,” JHEP 07 (2014) 079, arXiv:1405.0301 [hep-ph].
  * [23] C. Bierlich, S. Chakraborty, N. Desai, L. Gellersen, I. Helenius, P. Ilten, L. Lönnblad, S. Mrenna, S. Prestel, C. T. Preuss, T. Sjöstrand, P. Skands, M. Utheim, and R. Verheyen, “A comprehensive guide to the physics and usage of PYTHIA 8.3,” SciPost Phys. Codeb. 2022 (2022) 8, arXiv:2203.11601 [hep-ph].
  * [24] C. Bierlich, S. Chakraborty, N. Desai, L. Gellersen, I. Helenius, P. Ilten, L. Lönnblad, S. Mrenna, S. Prestel, C. T. Preuss, T. Sjöstrand, P. Skands, M. Utheim, and R. Verheyen, “Codebase release 8.3 for PYTHIA,” SciPost Phys. Codebases (2022) 8–r8.3.
  * [25] DELPHES 3 Collaboration, J. de Favereau, C. Delaere, P. Demin, A. Giammanco, V. Lemaître, A. Mertens, and M. Selvaggi, “DELPHES 3, A modular framework for fast simulation of a generic collider experiment,” JHEP 02 (2014) 057, arXiv:1307.6346 [hep-ex].
  * [26] R. D. Ball et al., “Parton distributions with LHC data,” Nucl. Phys. B 867 (2013) 244–289, arXiv:1207.1303 [hep-ph].
  * [27] J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, O. Mattelaer, H. S. Shao, T. Stelzer, P. Torrielli, and M. Zaro, “The automated computation of tree-level and next-to-leading order differential cross sections, and their matching to parton shower simulations,” JHEP 07 (2014) 079, arXiv:1405.0301 [hep-ph].
  * [28] R. Frederix, S. Frixione, V. Hirschi, D. Pagani, H. S. Shao, and M. Zaro, “The automation of next-to-leading order electroweak calculations,” JHEP 07 (2018) 185, arXiv:1804.10017 [hep-ph]. [Erratum: JHEP 11, 085 (2021)].
  * [29] DarkMachines community, “Unsupervised-Hackathon.” https://doi.org/10.5281/zenodo.3961917, July, 2020.