Artificial Proto-Modelling: Building Precursors of a Next Standard Model from Simplified Model Results

Unlike “full” models, the BSM particle content of a proto-model is not fixed. Since the number of degrees of freedom (including spin and multiplicity of states) for each particle is mostly relevant for the production cross section, which is treated as a free parameter, we do not specify the spin or multiplicity of the BSM particles. In the current version, we consider the following pool of 20 (SUSY inspired) particles:

• •

  Light quark partners $X_{q}$ ($q=u,d,c,s$): we allow for a single partner for each light-flavor quark. Unlike SUSY models, we do not consider two independent particles (left and right-handed squarks) for each flavor. As mentioned above, a possible multiplicity of (degenerate) states is accounted for by rescaling the production cross section.

• •

  Heavy quark partners $X_{b}^{i}$, $X_{t}^{i}$ ($i=1,2$): unlike the light quark partners, we consider two independent particles for each flavor. Since these particles have so extensively been searched for at the LHC, we include two states in order for the proto-models to have enough degrees of freedom to accommodate the data.

• •

  Gluon partner $X_{g}$: we introduce one new color-octet particle, analogous to a gluino in SUSY.

• •

  Electroweak partners $X_{W}^{i}$, $X_{Z}^{j}$ ($i=1,2$; $j=1,2,3$): we allow for two electrically charged and three neutral states. These might correspond to charginos and neutralinos in the MSSM (with the neutral higgsinos being exactly mass-degenerate), or to the scalars of an extended Higgs sector with a new conserved parity. The lightest neutral state ($X_{Z}^{1}$) is assumed to be the LBP.

• •

  Charged lepton partners $X_{\ell}$ ($\ell=e,\mu,\tau$): like for light-flavor quarks, we consider a single partner for each lepton flavor.

• •

  Neutrino partners $X_{\nu_{\ell}}$: again we consider one partner for each neutrino flavor ($\nu_{e},\nu_{\mu},\nu_{\tau}$).

As mentioned above, the lightest BSM particle is required to be the $X_{Z}^{1}$, hence all masses must satisfy $m(X)\geq m(X_{Z}^{1})$. In a given proto-model only particle masses below 2.4 TeV are considered, so they are within the LHC reach. States which appear in two- or three-fold multiplicities are mass ordered, e.g., $m({X_{t}^{2}})>m({X_{t}^{1}})$.

Some additional requirements are necessary for the masses of the colored new particles, to avoid that our machinery “discovers” light states in regions which are poorly covered in the database. Concretely, for the $X_{g}$ and $X_{q}$ we disallow masses below 310 GeV, since most of the simplified-model results we employ do not consider masses below this value. Also, for $X_{t}$ masses below 280 GeV, we disallow the region $150\;\mathrm{GeV}<m(X_{t})-m(X_{Z}^{1})<200\;\mathrm{GeV}$. The reason is that many CMS analyses do not make any statements for this “corridor” region, as events in there become too similar to $t\bar{t}$ events, see e.g. Ref. [15].

The decay modes of each new particle must of course be consistent with its quantum numbers. At present, we restrict the BSM particle decays to the channels shown in Table 1. Note that not all possibilities are considered in this work. In particular $X_{W}$ and $X_{Z}$ decays into $X_{\ell}+\ell/\nu$, $X_{\nu}+\nu/\ell$, and $X_{W,Z}+\gamma$ are not yet taken into account.

Specific decay modes are turned off if one of the daughter particles is not present in the proto-model, or if the decay is kinematically forbidden. The branching ratios (BRs) of the allowed channels are taken as free parameters, which however must add up to unity.

In addition to the masses and BRs of the BSM particles, their production cross sections are also allowed to vary freely. However, in order to have sensible starting values, the cross sections are first computed assuming the BSM particles to be MSSM-like, and then allowed to be rescaled freely by signal strength multipliers $\kappa$. For instance, the pair production of $X_{g}$ is given by:

with the mass of the gluino $\tilde{g}$ set to the $X_{g}$ mass. The rescaling factors $\kappa_{X_{i}X_{j}}$ are taken as free parameters of the proto-model.

In practice, we use a template SLHA input file for the MSSM to define masses and BRs of the proto-model (with $\tilde{q}_{R}$, $\tilde{\chi}^{0}_{4}$ and heavy Higgses decoupled); the reference SUSY cross sections are then computed with Pythia 8.2 [16] and NLLFast 3.1 [17, 18, 19, 20, 21, 22, 23]. We note that this way cross sections which are naturally suppressed or theoretically not allowed (such as $pp\to X_{e}^{-}X_{e}^{-}$) are automatically neglected. While this adds a certain bias, as processes which do not occur in the SUSY case (like production of two different top partners, $pp\to X_{t}^{1}\bar{X}_{t}^{2}$) will also be absent in the proto-model, this is not a problem, because there are currently no simplified-model results available for such cases.

The extent to which likelihoods can be computed crucially depends on the information available from the experimental collaboration.

1. 1.

   If only observed ULs are available, the likelihood becomes a constraint in the form of a step function at the observed 95% CL exclusion limit. This is in fact not useful for constructing $\mathrm{L}_{\mathbf{BSM}}$ per se, but will be used to determine the maximal allowed signal strength $\mu_{max}$ (see Section 3.2 below).

2. 2.

   If the expected ULs are available in addition to the observed ones, following [30] we approximate the likelihood as a truncated Gaussian:

   $$\mathrm{L}(\mu|D)=\frac{c}{\sqrt{2\pi}}\frac{\sigma_{\mathrm{ref}}}{{\bm{\sigma}}_{\mathrm{obs}}}\,e^{-(\mu\sigma_{\mathrm{ref}}-\sigma_{\mathrm{max}})^{2}/2{\bm{\sigma}}^{2}_{\mathrm{obs}}},\penalty\ \mathrm{for}\penalty\ \mu\geq 0\,.$$

   (3)

   Here, the likelihood is a function of the signal strength multiplier $\mu$ for any reference cross section $\sigma_{\mathrm{ref}}$ (in our case the cross section predicted for the given proto-model), while ${\bm{\sigma}}_{\mathrm{obs}}$ is approximated with the standard deviation of the expected Gaussian likelihood, ${\bm{\sigma}}_{\mathrm{obs}}\approx{\bm{\sigma}}_{\mathrm{exp}}=\frac{\sigma_{\mathrm{exp}}^{\mathrm{UL}}}{1.96}$, with $\sigma_{\mathrm{exp}}^{\mathrm{UL}}$ the expected 95% CL upper limit on the signal cross section. In addition, $\sigma_{\mathrm{max}}$ is chosen such that the approximate truncated Gaussian likelihood correctly reproduces the 95% CL observed limit on the production cross section, $\sigma_{\mathrm{obs}}^{\mathrm{UL}}$ . Finally, $c$ is a normalization constant, $c\geq 1$, which ensures that the integral over the likelihood from $\mu=0$ to infinity equals unity. While Eq. (3) is certainly a crude approximation in some cases, it does contain valuable information on possible excesses in the data.

   However, we cannot expect Eq. (3) to hold for large excesses. Therefore, if the observed UL differs from the expected one by more than two standard deviations, we “cap” the likelihood by artificially setting the observed UL to $\sigma_{\mathrm{obs}}^{\mathrm{UL}}=\sigma_{\mathrm{exp}}^{\mathrm{UL}}+2\bm{\sigma}_{\mathrm{exp}}$. With this procedure we avoid overly optimistic interpretations of excesses. While this affects only a small fraction of results, this is a crude approximation, and it could be avoided by EM results, see below.

3. 3.

   EM results contain more information to construct a proper likelihood. In the absence of a full statistical model (see point 4 below), we assume $p(\theta)$ to follow a Gaussian distribution centered around zero and with a variance of $\delta^{2}$, whereas $P(D)$ corresponds to a counting variable and is thus properly described by a Poissonian. The likelihood is then given by

   $$\mathrm{L}(\mu|D)=\frac{(\mu s+b+\theta)^{n_{obs}}e^{-(\mu s+b+\theta)}}{n_{obs}!}{\rm exp}\left(-\frac{\theta^{2}}{2\delta^{2}}\right)\,,$$

   (4)

   where $n_{obs}$ is the number of observed events in the signal region under consideration, $b$ is the number of expected background events, and $\delta^{2}=\delta_{b}^{2}+\delta_{s}^{2}$ the signal+background uncertainty. The nuisances $\theta$ can be profiled or marginalized over; our default is profiling. This is often referred to as a simplified likelihood [31, 32, 33]. While $n_{obs}$, $b$ and $\delta_{b}$ are directly taken from the experimental publications, we assume a default 20% systematical uncertainty on the signal. See [24] for more detail.

   Generally, in each analysis we use only the signal region with the highest $\sigma_{\mathrm{obs}}^{\mathrm{UL}}/\sigma_{\mathrm{exp}}^{\mathrm{UL}}$ ratio. A large ratio indicates an excess has been observed and allows us to identify potential dispersed signals. Although it is not a given that this is the best criterion when combining results from distinct analyses, it allows us to drastically reduce the number of possible signal region combinations across analyses.

4. 4.

   The best case is to have EM results together with a statistical model, which describes the correlations of uncertainties across signal regions [11, 10]. In this context, CMS sometimes provides covariance matrices, which allow for the combination of signal regions as discussed in [32]. (N.b., this is still a simplified likelihood assuming Gaussian uncertainties.) So far, however, a covariance matrix together with simplified model EMs is available for only one analysis [34].

   ATLAS, on the other hand, has recently started to publish full likelihoods using a plain-text (JSON) serialization [35] of the likelihood, which describes statistical and systematic uncertainties, and their correlations across signal regions, at the same fidelity as used in the experiment.³³ 3 This JSON format describes the HistFactory family of statistical models [36], used by the majority of ATLAS searches; it can conveniently be used via the pyhf package [37] in a scientific Python environment. At present, EM results with JSON likelihoods are available in the SModelS database for three ATLAS SUSY analyses [38, 39, 40] for full Run 2 luminosity (139 fb^-1), leading to clear improvements in the statistical evaluation for these analyses [29].⁴⁴ 4 While, as we will see, [38, 39, 40] do not impact our results, we stress that full likelihoods are a boon for reinterpretation studies. Once available on a broader basis, they will be extremely beneficial in particular for global analyses, including the one presented in this paper.

If we had perfect knowledge of the LHC results and their mutual correlations, we would simply compute a single joint likelihood for all results. However, cross-analysis correlations are very difficult to quantify and require a dedicated effort by the experimental collaborations. Examples are the combinations of Higgs measurements [41, 42], which do provide channel-by-channel correlations as auxiliary public material. For searches, statistical combinations of several analyses (considering different final states) are rare⁵⁵ 5 One example is the CMS combination of electroweakino searches [43]. and if done, the cross-analysis correlations have not been made public so far.

In this work, we therefore treat such correlations in a binary way—a given pair of analyses is either considered to be approximately uncorrelated, in which case it may be combined (by multiplying the respective likelihoods), or it is not at all considered for combination. We always assume results from different LHC runs, and from distinct experiments (ATLAS or CMS) to be approximately uncorrelated. Furthermore, we also treat results with clearly different final states in their signal regions (e.g., fully hadronic final states vs. final states with leptons) to be uncorrelated.

Figure 2 summarizes our treatment of correlations in the form of four binary analysis correlation matrices. We understand that, although well motivated, this is but an educated guess from our side. Future runs of our algorithm will be based on more reliable estimates of the results’ mutual statistical independence. Indeed work on systematically constructing such analyses correlation matrices is underway, see [44] and contribution 16 in [45]. However, while a good number of analyses can be combined in our approach, we also see from the white bins in Fig. 2, that about half (one third) of the 13 (8) TeV results are observed ULs only, and thus no proper likelihood can be computed for them. Nonetheless, under the above assumptions it is possible to construct an approximate combined likelihood for subsets of LHC results. We will refer to each of these subsets as a combination of results. Each combination must satisfy:

• •

  any pair of results in the subset must be considered as uncorrelated and

• •

  any result which is allowed to be added to the combination, is in fact added.

Likelihoods from uncorrelated analyses can simply be multiplied: $\mathrm{L}_{\mathbf{BSM}}(\mu)=\prod_{i=1}^{n}\mathrm{L}_{i}(\mu)$, where the product is over all $n$ uncorrelated analyses and $\mu$ is the global signal strength. Information from all the other analyses, which are not included in the combination, is accounted for as a constraint on the global signal strength $\mu$. Given the upper limit on the signal cross section obtained from the most sensitive analysis/signal region, we compute an upper limit on $\mu$:

where $\sigma$ is the relevant signal cross section and the 1.3 factor allows for a 30% violation of the 95% CL observed limit. We accept such mild transgressions in order to account for the fact that when simultaneously checking limits from a large number of analyses, a few are statistically allowed to be violated. Imposing the above limit on the global signal strength $\mu$ corresponds to truncating the likelihood at $\mu_{max}$ and ensures constraints from correlated analyses, as well as analyses for which no proper likelihood can be computed, are approximately taken into account.

The above procedure is illustrated in Fig. 3. The left panel in this figure shows how likelihoods are combined for a proto-model where it is possible to construct a likelihood for all the relevant results, which can all be considered as uncorrelated. In the example shown, “Result A” represents an analysis which favors $\mu>0$, while “Result B” corresponds to a “constraining analysis”. The latter observed an under fluctuation in the data with respect to the expected SM background, thus favoring $\mu\leq 0$. The grey line shows the maximum of the combined likelihood ($\hat{\mu}$), while the red line corresponds to the limit on $\mu$ ($\mu_{max}$) obtained from “Result B” alone.

The right panel in Fig. 3 displays our approximation for the case where it is not possible to construct a proper likelihood for the most constraining analysis (“Result B” in this example) or it cannot be considered as uncorrelated with the other results. In this case we truncate the likelihood for “Result A” at $\mu_{max}$ obtained from “Result B”. Hence, the value of $\hat{\mu}$ for the truncated likelihood shifts slightly with respect to the full result. In this work we apply this approximation whenever it is not possible to construct a full likelihood.

As mentioned above, during the MCMC-type walk, the algorithm aims to maximize the proto-model test statistic $K$. Naturally, $K$ must be designed such that it increases for models which better satisfy all the constraints (which includes better fitting potential dispersed signals). Furthermore, it is desirable to reduce the test statistic of models with too many degrees of freedom in order to enforce the law of parsimony, that is, Occam’s razor [46]. Given a proto-model, we define for each combination of results $c\in C$ the auxiliary quantity

Some explanations are in order regarding the choice of the prior. Since $\pi(\mathbf{SM})$ is a common factor for all combinations and does not affect the comparison between distinct proto-models, we wish to define the prior such that

Moreover, as mentioned above, the proto-model prior $\pi(\mathbf{BSM})$ should penalize the test statistic for newly introduced particles, branching ratios, or signal strength multipliers. We therefore heuristically choose the prior to be of the form

where $n_{\mathrm{particles}}$ is the number of new particles present in the proto-model, $n_{\mathrm{BRs}}$ is the number of non-trivial branching ratios⁶⁶ 6 This means, for only one decay mode with 100% BR, $n_{\mathrm{BRs}}=0$., and $n_{\mathrm{SSMs}}$ the number of signal strength multipliers.⁷⁷ 7 Note that particle-particle and particle-antiparticle production modes can have different signal strength multipliers. In order to favor democratic decays, if two decay channels for a given particle have BRs differing by less than 5%, they are counted only once. The parameters $a_{1}$, $a_{2}$, and $a_{3}$ are chosen to be 2, 4, and 8, respectively. This way, one particle with one non-trivial decay and two production modes is equivalent to one free parameter in the Akaike Information Criterion (AIC) [47]. For the SM, $n_{\mathrm{particles}}=n_{\mathrm{BRs}}=n_{\mathrm{ssm}}=0$ and the prior in Eq. (9) satisfies the normalization set by Eq. (8). The test statistic thus roughly corresponds to a $\Delta\chi^{2}$ of the proto-model with respect to the SM, with a penalty for the new degrees of freedom. Note that the prior is not normalized in the space of all proto-models.

In Table 2 we illustrate how the prior affects the model test statistic $K$ for a few choices of parameters. It should be noted that, with our choice of prior, $K$ ceases to have a clear, probabilistic interpretation. $K$ can even become negative. We interpret such an outcome as a preference of the SM over the BSM model.

In order to illustrate the walker algorithm, we first apply it to an extremely reduced version of the database with only three 13 TeV analyses: the ATLAS search for stops in the 1 lepton channel (ATLAS-SUSY-2016-16) [48], the CMS search for hadronically decaying tops (CMS-SUS-16-050) [49], and the generic CMS search in the multi-jet plus $E_{T}^{\rm miss}$  channel (CMS-SUS-19-006) [15]. All three present results for simplified models with 2 and 4 tops or b-jets in the final state.⁸⁸ 8 ATLAS-SUSY-2016-16 has EM results in addition to UL maps; the results from the two CMS analyses consist of expected and observed ULs.

As discussed previously, the aim of the algorithm is to find proto-models which can fit dispersed signals and at the same time satisfy the other constraints. In this example, the potential dispersed signals appear in the ATLAS and CMS stop searches, where $\sim 1$–$2\,\sigma$ excesses have been observed in some signal regions. The third analysis (CMS-SUS-19-006), on the other hand, has seen an under-fluctuation in data and will play the role of the constraining result, or critic. Indeed, since the two CMS searches included cannot be considered as uncorrelated (see Fig. 2), they are not allowed to be combined into a single likelihood. The ATLAS analysis, however, can be combined with either one of the CMS searches, and since the small excesses are found in ATLAS-SUSY-2016-16 and CMS-SUS-16-050, these will correspond to the best combination discussed at the beginning of Section 4.

In Fig. 4 we display the evolution of the walker during 200 steps, by showing the masses of the particles present in the proto-model every 10 steps. Each point is scaled according to the proto-model test statistic $K$, so bigger points correspond to higher $K$ values. Since the database considered only contains results for simplified models with tops and bottoms in the final state, proto-models without top and bottom partners find no corresponding data and no test statistic can be computed; hence they are set to zero. Around step 20, a proto-model with a top partner ($X_{t}^{1}$) is created, thus resulting in an increase in the test statistic. The following steps modify the $X_{Z}^{1}+X_{t}^{1}$ model by randomly changing its parameters and adding or removing particles. Since the addition of new states only reduces the model test statistic due to the prior (see Eq. (9)), the walker usually reverts back to the minimal scenario with one $X_{Z}$ and one $X_{t}$. We also see that for a few steps, such as step 70, the top partner is replaced by another state, resulting in a drastic reduction in $K$. After 200 steps, the highest scoring proto-model is that of step 94 and only contains the $X_{Z}^{1}$ and $X_{t}^{1}$ with masses of roughly 160 GeV and 1.1 TeV, respectively.

The behavior of the likelihoods at a few relevant steps is illustrated in Fig. 5, where the $\mathrm{L}_{\mathbf{BSM}}(\mu)$ curves, normalized by $\mathrm{L}_{\mathbf{SM}}$, are shown for the best combination of analyses, ATLAS-SUSY-2016-16 and CMS-SUS-16-050. The maximum signal strength $\mu_{\rm max}$ allowed by the constraining analysis (CMS-SUS-19-006) is indicated by the dash-dotted vertical lines. The first two panels show the steps immediately before the highest-score model. The first (leftmost) panel corresponds to a proto-model containing a $X_{\mu}$, $X_{t}^{1}$ and $X_{Z}^{1}$ with masses around 1.8 TeV, 1 TeV and 160 GeV, respectively. Since only the top partner contributes for fitting the excesses in the best combination, $X_{\mu}$ adds irrelevant degrees of freedom, thus reducing the proto-model score due to the penalty from the prior. In the following step, shown in the center panel, $X_{\mu}$ is removed, without impact on the likelihood ratio, as expected. Nonetheless, due to the change in prior, the test statistic $K$ increases by about one unit. The next step, shown in the right panel, modifies the $X_{t}$ mass, which gives a better fit to the excess and an increase of the likelihood ratio. With $K=3.65$ this turns out to be the highest-score model of the toy walk.

Although better results are likely to be found with a larger number of steps, this example illustrates how the algorithm walks in the proto-model parameter space to find models consistent with dispersed signals and the constraints. In the following section, we will present more meaningful results, found using the same algorithm but for the full SModelS database and with a much larger number of steps.

In order to compute the significance of the test statistic and, ultimately, its global $p$-value, we produce “MC toy” versions of the database, in which we replace the observed data by values sampled from the SM background models. More specifically,

• •

  for EM-type results we sample a normal distribution with a mean of the background estimate, and a variance of the squared error of the background estimate. The sampled value is then entered as the lambda parameter of a Poissonian distribution. Finally, the value drawn from the Poissonian is used as the “fake” observation, i.e. the fake event yield;

• •

  for UL-type results, for which expected upper limits are known to SModelS, we estimate the SM uncertainty on the expected limit assuming a large number of observed events and small systematical uncertainties [30]. The observed upper limits are then sampled from a Gaussian distribution:

  $$\sigma^{\mathrm{UL}}_{\mathrm{obs}}\sim\mathrm{Norm}\left(\sigma^{\mathrm{UL}}_{\mathrm{exp}},\left(\delta\sigma^{\mathrm{UL}}_{\mathrm{exp}}\right)^{2}\right),$$

  (10)

  where $\delta\sigma^{\mathrm{UL}}_{\mathrm{exp}}=\sigma^{\mathrm{UL}}_{\mathrm{exp}}/1.96$ is the estimated width of the distribution;

• •

  UL-type results, for which only observed upper limits are available, enter the fake database as they are, since it is not possible to estimate the SM uncertainty in this case.

Using the sampled values described above, we produce 50 fake databases and apply the walker algorithm to all of them. The proto-models with the highest $K$ values from each of these runs are then used to estimate the density of the test statistic $K$ under the SM hypothesis via a Kernel density estimator:

Here, $N=50$ is the number of runs with fake databases, and $K^{i}_{\mathrm{fake}}$ the test statistic for the $i$-th run; $\mathrm{kern}_{w}$ denotes the choice of Kernel in the Kernel density estimation. We choose a Gaussian with a width $w$ determined by Scott’s rule [50]. Figure 6 shows the $K^{i}_{\mathrm{fake}}$ values generated for each run as well as the result obtained for $\rho(K)$. As we can see, under the SM hypothesis, values up to $K\approx 11$ are expected, with the density peaking at $K\approx 4$.

As discussed in Section 4, the walker algorithm looks for dispersed signals, with $K$ computed as the maximum over all possible combinations of results. Therefore, large $K$ values are expected to be found even when considering the SM hypothesis, where upward background fluctuations can mimic dispersed signals. In addition, it has been suggested [51] that the experimental results tend to be conservative and overestimate the background uncertainties.

We can strengthen this claim by considering the $p$-values computed for all signal regions of the EM results contained in the database. If the observed data is well described by the estimated backgrounds and the uncertainties thereon, one would expect a flat distribution of $p$-values. In fact, for the fake databases generated under the SM hypothesis, this is the case, as shown by the red histogram in Fig. 7. When we consider the real data, however, the distribution is not flat, as seen from the blue histogram in the same figure. This is expected if the background uncertainties are overestimated. As a consequence, the density $\rho(K)$ of the SM hypothesis in Fig. 6 is conservative; we can safely assume that the actual density would be shifted towards lower values of $K$ and that our results generally lie on the conservative side.

We are now ready to apply the walker algorithm to the actual SModelS database with the real LHC results. We perform 10 runs, each employing $50$ walkers and $1,000$ steps/walker. The results are summarized in Fig. 8, where we display the proto-models with the highest $K$ value from each run. Besides the $X_{Z}^{1}$ LBP, all models include one top-partner, $X_{t}^{1}$, and one light-flavor quark partner, $X_{d,c}$, and their test statistics are at $K=6.76\pm 0.08$ thus showing the stability of the algorithm. The $X_{\mu}$ particle introduced in run 5 is due to small $\approx 1\sigma$ excesses in the CMS search for sleptons, CMS-SUS-17-009 [52] and the ATLAS search for electroweakinos, ATLAS-SUSY-2016-24 [53]; the prior’s penalty for introducing the additional particle, $\Delta K=-1.25$, is overruled by the increase in the likelihood ratio, $\Delta K=1.94$.

Although the scenarios from Fig. 8 are statistically all equally plausible, it is instructive to discuss the absolute “winning” one, i.e. the proto-model with $K=6.90$ generated in step 582 of the 29th walker in run #9, in some more detail. It has $X_{t}^{1}$, $X_{d}$ and $X^{1}_{Z}$ masses of 1166, 735 and 163 GeV, respectively, and produces signals in the $t\bar{t}+E_{T}^{\rm miss}$ and jets$+E_{T}^{\rm miss}$ final states with SUSY-like cross sections. Concretely, the effective signal strength multipliers are $\hat{\mu}\times\kappa_{X_{t}^{1}\bar{X}_{t}^{1}}\approx 1.2$ and $\hat{\mu}\times\kappa_{X_{d}\bar{X}_{d}}\approx 0.5$, corresponding to $\sigma(pp\to X_{t}^{1}\bar{X}_{t}^{1})\approx 2.6$ fb and $\sigma(pp\to X_{d}\bar{X}_{d})\approx 24$ fb at $\sqrt{s}=13$ TeV, and both the $X_{t}^{1}$ and $X_{d}$ directly decay to the lightest state ($X^{1}_{Z}$) with 100% BR.

The analyses, which contain possible dispersed signals and drive this model (as well as the other ones in Fig. 8), are listed in Table 3. The light quark partner with mass around 700 GeV is introduced in order to fit the observed excesses in the multijet (i.e., $0\ell$ + jets) +$E_{T}^{\rm miss}$ analyses from ATLAS at $\sqrt{s}=8$ and 13 TeV [54, 55] and from CMS at $\sqrt{s}=8$ TeV [56]. The $K$ value obtained is, however, mostly driven by the $\approx 1.5\,\sigma$ and $2\,\sigma$ excesses observed in the 13 TeV CMS and ATLAS stop searches [49, 48], which lead to the introduction of a top partner ($X_{t}^{1}$) with a mass around 1.2 TeV. Despite corresponding to small excesses, identifying the presence of such potential dispersed signals is one of the main goals of the algorithm presented here. Furthermore, the fact that these excesses appear in distinct ATLAS and CMS analyses and can be explained by the introduction of a single top partner is another interesting outcome of the whole procedure.

It is worth noting that the signal injected by the winning proto-model is typically smaller than the one favored by the excesses, as seen when comparing the Obs and the Signal columns in Table 3. This is due to a tension with the constraints imposed on the signal strength by all the other analyses in the database (the critic). The model discussed above is mostly constrained by the 13 TeV CMS multijet analyses for 36 fb^-1 [57, 58] and 137 fb^-1 [15], which all observed under-fluctuations with respect to the expected background, see Table 4. The most serious challenges comes from the 2 jets+$E_{T}^{\rm miss}$ and $t\bar{t}$+$E_{T}^{\rm miss}$ limits of [15], for which the proto-model under consideration exceeds both the observed and expected UL. Concretely, the analysis [15] gives $r_{\rm obs}$ $(r_{\rm exp})=1.3$ ($1.11$) for the $X_{d}$, and $r_{\rm obs}$ $(r_{\rm exp})=1.28$ ($1.26$) for the $X_{t}^{1}$. ¹⁰¹⁰ 10 Recall that we accept mild transgressions up to $r\lesssim 1.3$ to account for the fact that when simultaneously checking limits from a large number of analyses, a few are statistically allowed to be violated.

Finally, we point out that very similar conclusions are drawn for all the high score proto-models shown in Fig. 8. Of course, all these conclusions may be challenged (or strengthened) by new results for full Run 2 luminosity, which are not yet included in the SModelS database, e.g. the stop searches [60, 61].

Some comments are in order on the significance of the excesses and $K$ values found for the models in Fig. 8. In Fig. 9 we show the $K$ distribution found under the SM-only hypothesis (see the discussion in Section 5.1) and the values for the 10 highest-scoring proto-models obtained with the “real” database (green stars). As we can see, despite being at the tail of the distribution, the values still seem compatible with the SM-only hypothesis. In order to quantify this compatibility, we compute the $p$-value given by the relative frequency of the $K$ values under the SM hypothesis which were found above the observed $K$ values ($K_{\mathrm{obs}}$):

Here $N=50$ is the number of fake SM-only runs, $\mathds{1}_{X}(x)$ denotes the indicator function that equals unity for all $x\in X$ and zero otherwise, $\bar{K}_{\mathrm{obs}}$ is the average of $K_{\mathrm{obs}}$ of the ten runs over the “real” database, and $\rho(K)$ refers to the estimated density of the test statistic as defined in Eq. (11). The resulting value corresponds to the shaded area under the curve shown in Fig. 9. Using the $K$ values in Fig. 8 we obtain the global $p$-value for the SM hypothesis:

Our results thus indicate a very mild disagreement with the SM hypothesis. As already pointed out in Section 5.1, we expect this to be conservative due to potentially overestimated background uncertainties. Furthermore, this result requires no look-elsewhere correction, since the SM density $\rho(K)$ was derived through the exact same procedure as applied for the real database.

Complementary information on the distributions of the test statistic $K$ and posterior densities can be found in appendix A.2.

As discussed above, the $K$ values obtained for the real database are in a mild disagreement with the SM hypothesis. It is relevant, however, to investigate whether the procedure proposed in this work would be able to reproduce the underlying model. For this kind of closure test, we employ a scheme analogous to the one used to produce fake databases under the SM-only hypothesis in Section 5.1. This time, however, using the highest-score model from Fig. 8, we create fake observed data under the “BSM hypothesis”. Concretely, we sample from the signal plus background model, assuming that the corresponding uncertainties are dominated by the errors on the background. However, when producing the fake “BSM hypothesis” databases, we scale the background uncertainties by a “fudge factor” of 0.65, aiming to produce fake data closely resembling those observed in the real database. This is motivated by the fact that, as discussed in Section 5.1, the observed data in the real database point to an overestimation of the background uncertainties.

The theory predictions from the injected proto-model are included as follows: for EM-type results the signal + background model has central values given by the SM expected value plus $\sigma\times\mathrm{BR}\times\epsilon$ for the injected signal. The number of observed events is then determined by sampling the signal + background model. For UL-type results, on the other hand, we simply shift the observed limit by the predicted cross section in the region of the UL map that has masses within 100 GeV of the assumed proto-model values. This procedure is a good approximation within the Gaussian limit if the signal efficiencies do not change significantly within the shifted region.

We generate 10 such “fake signal” databases and perform a run with each of them. The resulting distribution of the test statistics $K$ is shown in Fig. 9. As can be seen, under the BSM hypothesis (the fake signals) $K$ varies between $K\approx 5.8$ and $K\approx 12.3$. Figure 10 shows the particle spectra of the high-score models of these fake signal runs and compares them to the injected signal. We see that $X_{t}^{1}$ was reconstructed nicely at around the right mass scale, in 8 out of 10 runs. The quark partner, $X_{q}$, was found in 9 out of 10 runs, though once at a significantly lower mass. Finally, runs 5, 6, and 9 introduced spurious lepton partners, all due to various background fluctuations created in the statistical sampling.

All in all, the differences between the injected and reconstructed signals shown in Fig. 10 are statistical in nature and in fact to be expected in our procedure, given that the injected signal is of the size of the background uncertainties. To verify this, we have also performed a purely technical closure test, where background and signal sampling were turned off, so the number of observed events for EM results is simply set to the corresponding central values of the signal plus background model. In this case, the particle content of the injected signal is correctly reproduced in all runs, including the masses. The result from injecting the high-score model from Fig. 8 is shown in Fig. 11. It nicely confirms (i) the robustness of our procedure, and (ii) the conclusion that the small deviations in Fig. 10 are due to fluctuations generated by the statistical sampling.

The computational results presented were obtained using the CLIP cluster [65]. We thank Nishita Desai, Sylvain Fichet, Jan Heisig, Suchita Kulkarni, Harrison Prosper and Helga Wagner for discussions related to this work. The work of S.K. was supported in part by the IN2P3 master project “Théorie – BSMGA”.