1 Introduction↩︎

Astrophysical evidence points to a non-baryonic dark component of the Universe, but the microscopic identity of dark matter (DM) remains unknown. Among the experimental approaches, high-energy colliders can produce DM in controlled environments and probe its couplings to Standard Model (SM) particles through missing transverse momentum (\(E_{\mathrm{T}}^{\mathrm{miss}}\)) recoiling against visible final-state objects [3]. The mono-\(Z\) topology—a leptonically decaying \(Z\) boson balanced by \(E_{\mathrm{T}}^{\mathrm{miss}}\)—offers a relatively clean signature: the dilepton system constrains the \(Z\) mass, electroweak backgrounds are smaller than in hadronic mono-jet searches, and the channel is sensitive to simplified models with \(s\)-channel mediators that also radiate an on-shell \(Z\) boson [4]–[6].

Published LHC searches in this final state have progressed from early reinterpretations and multivariate projections to full Run 2 results from CMS and ATLAS [7]–[9]. These analyses typically bin events in \(E_{\mathrm{T}}^{\mathrm{miss}}\) or related kinematic variables, fit parametric or template background models, and compare data to signal hypotheses defined within the ATLAS/CMS Dark Matter Forum benchmark framework [6]. Complementary mono-\(Z\) and mono-\(Z'\) phenomenology has also been studied through dilepton mass spectra, angular distributions, and resonance-plus-\(E_{\mathrm{T}}^{\mathrm{miss}}\) topologies [10]–[13].

We investigate an alternative strategy that keeps the mono-\(Z\) event selection but replaces hand-crafted discriminant variables with a likelihood-ratio score built from learned event densities. Neural Spline Flows (NSFs) [1], [2] model the SM background and three DM mediator hypotheses as continuous distributions over a fixed set of cleaned kinematic features; the per-event score is the log-density difference between the selected DM hypothesis and the channel-specific SM hypothesis. Normalizing flows supply exact density evaluation for correlated inputs [2], [14], and surjective or simulation-aware extensions have been explored for collider event modeling and anomaly detection [15], [16]. Our implementation uses CMS Run 2015D open data [17], [18] for SM-dominated control samples and publicly released MonoZToLL MC samples with vector, axial-vector, and scalar mediators [19]–[22] for signal studies, with one DM flow trained per mediator model and identical physics-feature definitions in the \(\mu\mu\) and \(ee\) channels.

Section 2 reviews the mono-\(Z\) literature and related machine-learning methods in more detail. Section 1.1 states the analysis regions and score definition used here; Section 1.2 summarizes our contributions and the paper structure. All region definitions, training splits, and numerical results reported below follow our analysis pipeline unless explicitly attributed to external work.

1.1 Analysis objective↩︎

The physics goal is to test whether the number and kinematics of events in a signal region are consistent with SM production of \(Z+E_{\mathrm{T}}^{\mathrm{miss}}\), or require a DM contribution. The working signal hypotheses are DM pair production with vector, axial-vector, and scalar mediators, simulated with the same lepton selection as data.

Operationally, we partition events into a control region (CR) and a signal region (SR). The CR (\(E_{\mathrm{T}}^{\mathrm{miss}}<50\,\text{GeV}\)) is dominated by Drell–Yan \(Z\)+jets and is used only to train SM density models; it is disjoint from the SR by construction. The SR applies \(E_{\mathrm{T}}^{\mathrm{miss}}\ge 50\,\text{GeV}\), \(|\Delta\phi(E_{\mathrm{T}}^{\mathrm{miss}},Z)|>2.5\), and \(n_{\mathrm{jets}}\le 1\). The SR imposes no explicit upper \(E_{\mathrm{T}}^{\mathrm{miss}}\)requirement; in practice the cleaned dataset caps \(E_{\mathrm{T}}^{\mathrm{miss}}\)at approximately 200  (Section 3), so events enter the fit through \(E_{\mathrm{T}}^{\mathrm{miss}}\) values up to this bound. The validation interval (\(50\leE_{\mathrm{T}}^{\mathrm{miss}}<100\,\text{GeV}\)) additionally serves as a background-normalisation sideband in the simultaneous SR+VR fit: it constrains the per-channel normalisation nuisances using VR data independently of the SR, providing a genuine sideband-based background prediction.

Two SM flows are trained, one on CR double-muon events and one on CR double-electron events (70%/30% train/validation split in each channel). Three DM flows are trained, one each on the full vector, axial-vector, and scalar MC samples, and are shared across lepton channels. The branch lepton_flavor (\(11\) for electrons, \(13\) for muons) selects the channel but is not included in the NSF input vector.

For an event \(\mathbf{x}\) in the SR and mediator hypothesis \(h\in\{\mathrm{vector},\mathrm{axial},\mathrm{scalar}\}\), \[\mathcal{S}_h(\mathbf{x}) = \log p(\mathbf{x}\mid\mathrm{DM}_h) - \log p(\mathbf{x}\mid\mathrm{SM}_{\mu\mu\;\mathrm{or}\;ee}). \label{eq:score}\tag{1}\] Large positive values indicate kinematics more characteristic of the trained DM sample than of the corresponding SM CR model.

Before interpreting the SR, we check the SM-only flow in a 50–100 \(E_{\mathrm{T}}^{\mathrm{miss}}\) validation band to verify stable channel-specific SM behavior. The primary search constructs six SR score arrays, \(\mathcal{S}_h^{\mu\mu}\) and \(\mathcal{S}_h^{ee}\) for the three mediator hypotheses, and fits three combined binned \(\mu\mu\)+ee likelihoods, one per signal hypothesis [23].

1.2 Contributions and paper outline↩︎

This paper documents:

1. extraction of 40 event-level features from CMS open data and MonoZToLL MC;

2. a cleaning stage with fixed bounds, quantile-based tail rules, and reference-aware caps for MC samples;

3. CR/SR definitions, NSF training protocol, and likelihood-ratio search strategy;

4. final SR profile-likelihood fits and CL\(_s\) limits for the Run 2015D dataset.

Section 3 lists data and MC inputs. Section 4 describes feature extraction and cleaning. Section 5 defines CR, validation, and SR event categories. Section 6 summarizes NSF architecture and training. Section 7 details the score, fit, and cross-checks. Section 8 presents distributions and limits. Section 9 concludes.

2 Literature Review↩︎

2.1 Collider constraints and mono-\(Z\) phenomenology↩︎

Early collider studies showed that effective-field-theory (EFT) descriptions of DM coupling to quarks and gluons can yield limits on interaction scales that complement, and in some kinematic regimes surpass, direct-detection bounds—particularly for light DM and spin-dependent scattering [3]. The mono-\(Z\) channel was proposed as an electroweak counterpart to mono-jet and mono-photon searches: a \(Z\) boson recoils against a pair of stable DM particles, with the \(Z\) emitted from initial-state quarks or an internal line in a simplified model [4]. Ref. [5] extended mono-\(X\) reinterpretations to events with a reconstructed \(Z\to\ell^+\ell^-\) system, deriving limits on EFT operators with quark bilinear and electroweak-boson couplings from ATLAS \(Z\)+\(E_{\mathrm{T}}^{\mathrm{miss}}\) measurements.

Related signatures broaden the mono-\(Z\) program beyond a single on-shell \(Z\) recoiling against DM alone. Ref. [10] studied \(E_{\mathrm{T}}^{\mathrm{miss}}\) in association with a dilepton or dijet resonance from a \(Z'\) boson, highlighting sensitivity in mass regions where standard resonance searches lose efficiency. Ref. [11] pointed out that DM with scalar mediators can induce loop-level box corrections to Drell–Yan production, producing a “monocline” feature near \(m_{\ell\ell}\sim 2m_\chi\) rather than a narrow resonance. Ref. [24] discussed pseudoscalar-portal fermion DM in two-Higgs-doublet setups, where chiral couplings suppress direct detection while preserving thermal relic targets. Ref. [25] compared mono-\(X\) searches with direct searches for the mediating particle across several simplified models, finding that mono-\(X\) is most competitive when the spectrum is compressed or when electroweak bosons participate directly in the mediator–DM coupling. Ref. [12] showed that lepton angular distributions in the Collins–Soper frame can distinguish spin-0, spin-1, and spin-2 mediator scenarios in mono-\(Z\) production and that shape information beyond inclusive \(E_{\mathrm{T}}^{\mathrm{miss}}\) can strengthen coupling limits.

2.2 Simplified models, benchmarks, and experimental searches↩︎

The ATLAS/CMS Dark Matter Forum report [6] codified simplified models —\(s\)-channel vector, axial-vector, scalar, and pseudoscalar mediators, \(t\)-channel colored mediators, and EFT benchmarks—together with recommended parameter scans and implementation guidance for Run 2 \(E_{\mathrm{T}}^{\mathrm{miss}}+X\) analyses. These models underpin the MonoZToLL open-data samples used in our signal studies; vector, axial-vector, and scalar mediator samples are each extracted and assigned a separate DM flow for hypothesis-specific scoring [6], [19]–[22].

On the experimental side, CMS published mono-\(Z\) searches at \(\sqrt{s}=13\) TeV interpreting spin-0 and spin-1 mediators, invisible Higgs decays, and other \(E_{\mathrm{T}}^{\mathrm{miss}}\)+\(\ell^+\ell^-\) scenarios [8], [9]. The most recent CMS result uses the full Run 2 dataset (\(137~\mathrm{fb}^{-1}\)) and sets limits on mediator and DM masses for vector, axial-vector, scalar, and pseudoscalar mediators, with comparisons to direct-detection cross sections [9]. Ref. [7] projected the mono-\(Z\) reach at \(13\) TeV with a multivariate likelihood discriminant, reporting roughly a factor-of-two improvement over cut-and-count methods once background normalization systematics exceed a few percent. High-mass dilepton searches from ATLAS [26], [27] bound \(Z'\) resonances and four-fermion contact interactions; while not mono-\(Z\) analyses, they constrain overlapping \(Z\)-associated new-physics scenarios probed in Refs. [10], [13]. Ref. [13] applied CMS open data from Run 1 to dimuon angular distributions at high \(m_{\mu\mu}\), interpreting the data in a mono-\(Z'\) simplified model and finding agreement with SM Drell–Yan expectations.

2.3 Density estimation and learning-based search methods↩︎

Normalizing flows map a simple base density to a target distribution through invertible transformations, enabling both sampling and exact likelihood evaluation [2]. Neural Spline Flows replace elementwise affine couplings with monotonic rational-quadratic splines, increasing expressivity while retaining analytic inversion [1]. Ref. [14] examined how several flow architectures scale with input dimensionality on toy HEP-like datasets, a relevant consideration for multivariate collider features. Ref. [15] introduced surjective and stochastic flow layers to handle permutation symmetry, variable particle multiplicities, and discrete quantum numbers in matrix-element-level and detector-level event generation, as well as anomaly-detection studies. Ref. [16] proposed simulation-assisted likelihood-free anomaly detection (SALAD), reweighting background simulation to data in sidebands before supervised learning in a signal-sensitive region—a hybrid of data-driven and simulation-based approaches distinct from our explicit two-hypothesis density ratio. Closely related anomaly-detection methods that use sideband density estimates include ANODE [28] and CATHODE [29]; both construct signal-region classifiers from control-region density ratios, a strategy adjacent to our direct NSF log-ratio score. Ref. [7] remains the closest prior mono-\(Z\) reference for multivariate discrimination, though it used engineered kinematic inputs and a likelihood function rather than normalizing-flow densities trained separately on control-region data and signal MC.

Our analysis sits at the intersection of these threads: we adopt the mono-\(Z\) simplified-model context and CR/SR philosophy of the experimental and Forum literature [6], [9], but estimate SM and DM event densities with NSFs and use their log-ratio as the test statistic, following asymptotic inference prescriptions [23] in the signal region. To our knowledge, a published mono-\(Z\) DM search combining NSF-based likelihood ratios with CMS open Run 2015D data in both \(\mu\mu\) and \(ee\) channels has not been reported; the present work documents that pipeline through the final SR likelihood fit.

3 Datasets and Signal Model↩︎

SM events come from CMS Run 2015D MINIAOD open data: the DoubleMuon primary dataset (record 24127) [17] and the DoubleEG primary dataset (record 24132) [18], both at \(\sqrt{s}=13\) TeV. These samples provide opposite-sign dilepton events after analyzer-level selection; they are used as SM-dominated inputs for CR training and for the SR search. The Run 2015D dataset corresponds to an integrated luminosity of approximately \(2.32~\mathrm{fb}^{-1}\). The extracted local TTrees contain 1,738,984 selected DoubleMuon events split across four output files and 1,418,467 selected DoubleEG events split across eight output files. The stored SM data schema contains the 40 physics features plus audit-only provenance and quality-control branches used to validate trigger and \(b\)-tag behavior; the row counts are summarized in Appendix 10.1, and the complete feature schema with formulas is given in Appendix 11.1.

The signal inputs are MonoZToLL MINIAODSIM samples generated in the same \(\sqrt{s}=13\) TeV era for vector, axial-vector, and scalar mediators [6]. The vector sample is record 16630 [19]; the axial-vector inputs combine records 16575 and 16597 [20], [21]; the scalar sample is record 16601 [22]. The extracted selected samples contain 27,602 vector-mediator events, 45,710 axial-vector events, and 32,792 scalar-mediator events. The raw MC inputs are organized as ROOT source files; the axial-vector training set combines four source files (one for \(M_\chi=10\,\text{GeV}\), \(M_V=20\,\text{GeV}\) and three for \(M_\chi=50\,\text{GeV}\), \(M_V=200\,\text{GeV}\)), totaling about 100,000 raw entries before the final mono-\(Z\) selection. The same 40 physics branches listed in Appendix 11.1 are written for data and MC; the MC output adds only the auxiliary gen_weight branch. Although no explicit upper \(E_{\mathrm{T}}^{\mathrm{miss}}\)bound is imposed in the selection, the data-driven tail clipping applied during feature cleaning caps \(E_{\mathrm{T}}^{\mathrm{miss}}\)at approximately 200 in the extracted dataset; events above this value are absent from the SR score distribution. The axial-vector signal sample combines two physically distinct benchmark points — \(M_\chi=10\,\text{GeV}\), \(M_V=20\,\text{GeV}\) (\(\sigma=1.856\,\mathrm{pb}\)) and \(M_\chi=50\,\text{GeV}\), \(M_V=200\,\text{GeV}\) (\(\sigma=0.158\,\mathrm{pb}\)) — stored across four ROOT files. Because the two points differ in cross-section by a factor of \(\sim\)​12, a single effective cross-section would not correspond to any real benchmark; cross-section limits are therefore reported separately for each axial-vector point using the same fitted \(\mu\). Generator-level cross-sections were extracted directly from the GenRunInfoProduct branch of each MINIAODSIM ROOT file.

4 Feature Extraction and Cleaning↩︎

Raw MINIAOD and MINIAODSIM files are read directly with uproot, with jagged collections handled in awkward and converted to fixed per-event NumPy arrays after object selection. The analysis products are ROOT Events TTrees written with uproot.recreate, mktree, and extend; CSV files are used only for cut-flow, file-audit, and EDA summary outputs. Per-event quantities include \(E_{\mathrm{T}}^{\mathrm{miss}}\), dilepton and \(Z\) kinematics, jet and \(b\)-tag summaries, hadronic recoil variables, and angular correlations (40 physics features in the extracted schema). The complete list of stored branches and formulas for derived observables is given in Appendix 11.1; Table 1 lists the selection cuts implemented in the extraction notebooks.

The DoubleMuon and DoubleEG extractors use a custom trigger-object decoder because the 2015 MiniAOD trigger information is not exposed as simple one-branch-per-HLT-path booleans. The decoder reads the selectedPatTrigger path-name, pathLastFilterAccepted, and filter-label collections. When path names are available, events are matched to the channel-specific HLT prefixes; when path names are empty, the extractor falls back to accepted final-filter labels such as the double-muon or double-electron filter tokens. Trigger decisions are accumulated in TRIGGER_AUDIT and the retained event mask is applied before physics-feature writing. For the data \(b\) veto, the notebooks also decode the MiniAOD pat::Jet::pairDiscriVector_ payload using uproot’s AsBinary interpretation. A slow decoder discovers the discriminator labels, a cached fixed-layout fast path extracts the CSVv2 scores, and a validation step checks the fast decoder against the slow decoder on sampled events. The medium 2015/76X CSVv2 working point (\(0.800\)) is then used to compute n_bjets, max_btag_csvv2, and the \(b\)-veto audit branches.

Cleaning is applied only to SM data TTrees. The SM cleaner clips fixed-domain variables, estimates high-quantile caps for heavy-tailed variables, and uses the opposite SM lepton channel as the reference sample when deriving channel-consistent caps. The DoubleMuon cleaner therefore references DoubleEG parts, and the DoubleEG cleaner references DoubleMuon parts. The merged SM bounds are exported to Output/sm_feature_bounds.json; DM MC ROOT files are not cleaned by cleaner.py. Instead, the exported SM bounds are applied later at DataLoader time, so genuine high-\(E_{\mathrm{T}}^{\mathrm{miss}}\)or high-\(p_{\mathrm{T}}\) DM tails are preserved in the extracted MC files until the training pipeline applies the common SM feature domain. The fixed bounds and tail-rule families used by the cleaning script are summarized in Appendix 12.1.

The cleaned NSF feature manifest retains 35 extracted physics features and excludes lepton_flavor, n_bjets, gen_weight, rho, n_vertices, and raw HT; lepton_flavor is retained only for channel routing and is not passed to any flow. Two binary jet-presence indicators computed at DataLoader time (has_jet1\(=\mathbb{1}[n_\mathrm{jets}\ge 1]\), has_jet2\(=\mathbb{1}[n_\mathrm{jets}\ge 2]\)) are appended to form the final 37-dimensional NSF input.

5 Control and Signal Regions↩︎

Region masks are built from \(E_{\mathrm{T}}^{\mathrm{miss}}\), \(|\Delta\phi(E_{\mathrm{T}}^{\mathrm{miss}},Z)|\), and \(n_{\mathrm{jets}}\) after cleaning:

• CR: \(E_{\mathrm{T}}^{\mathrm{miss}}< 50\,\text{GeV}\).

• Validation: \(50 \le E_{\mathrm{T}}^{\mathrm{miss}}< 100\,\text{GeV}\), with the same \(|\Delta\phi(E_{\mathrm{T}}^{\mathrm{miss}},Z)|\) and jet cuts as the SR.

• SR: \(E_{\mathrm{T}}^{\mathrm{miss}}\ge 50\,\text{GeV}\), \(|\Delta\phi(E_{\mathrm{T}}^{\mathrm{miss}},Z)| > 2.5\), and \(n_{\mathrm{jets}}\le 1\).

The CR (\(E_{\mathrm{T}}^{\mathrm{miss}}<50\,\text{GeV}\)) and SR (\(E_{\mathrm{T}}^{\mathrm{miss}}\ge 50\,\text{GeV}\)) are disjoint by construction at the \(E_{\mathrm{T}}^{\mathrm{miss}}=50\,\text{GeV}\) boundary. The validation interval lies inside the SR \(E_{\mathrm{T}}^{\mathrm{miss}}\)range and is reserved for SM closure checks before interpreting the full SR scores in data. The SR imposes no explicit upper \(E_{\mathrm{T}}^{\mathrm{miss}}\)requirement; in practice the cleaned dataset caps \(E_{\mathrm{T}}^{\mathrm{miss}}\)at approximately 200  (Section 3).

6 Neural Spline Flow Models↩︎

Each NSF maps a feature vector \(\mathbf{x}\in\mathbb{R}^d\) to a base density through a composition of coupling layers with rational-quadratic spline transforms [1]. We implement flows in PyTorch using zuko , with hyperparameters (layer count, hidden width, bin count) fixed before training and recorded in run manifests.

Training protocol:

• NSF_SM_mumu: CR double-muon data, 70% train / 30% validation.

• NSF_SM_ee: CR double-electron data, 70% train / 30% validation.

• NSF_DM_vector: full vector-mediator MC, 70% train / 30% validation, all flavors.

• NSF_DM_axial: full axial-vector MC, 70% train / 30% validation, all flavors.

• NSF_DM_scalar: full scalar-mediator MC, 70% train / 30% validation, all flavors.

\(\dagger\) The significance is capped at 8.0 due to floating-point underflow of the chi-square survival function at large \(q_0\). The underlying \(q_0\) values (233–327) indicate a background-modelling residual (Section 7.2), not a genuine \(8\sigma\) discovery. Early stopping uses validation negative log-likelihood.

The five NSFs share the following hyperparameters: 8 masked autoregressive coupling transforms with rational-quadratic splines (8 bins per spline, randomised feature permutation per transform), conditioner networks with two hidden layers of 256 units each, and a standard-normal base distribution. Training uses the Adam optimiser (\(\eta=3\times10^{-4}\), \(\lambda=10^{-5}\) weight decay) with a 500-step linear warm-up followed by cosine annealing to \(\eta_\mathrm{min}=10^{-6}\) over 200 epochs, a mini-batch size of 2048, and gradient clipping at \(\|\nabla\|_2\le 1\). Early stopping monitors validation NLL with a patience of 20 epochs and restores the best-epoch weights. In the present training runs, validation NLL continued to decrease at epoch 200 for all five flows; early stopping was not triggered for any flow before reaching MAX_EPOCHS. Every _best.pt checkpoint therefore corresponds to the final epoch rather than a converged minimum, and we flag potential residual undertraining as a limitation affecting the learned densities and reported limits. The discrete feature n_jets is dequantized during training by adding uniform noise \(\mathcal{U}(-0.5,0.5)\); raw integer values are used at inference.

All five flows operate in the same standardized feature space. SM flows fit a per-feature \(z\)-score scaler on their respective training splits and save the parameters to a JSON manifest. DM flows reuse the NSF_SM_mumu scaler rather than fitting new scale parameters, ensuring that DM event densities are evaluated in the SM feature domain. Before standardization, DM MC features are clipped to the SM bounds exported in Output/sm_feature_bounds.json; this preserves genuine high-\(E_{\mathrm{T}}^{\mathrm{miss}}\)tails while preventing out-of-domain extrapolation of the learned SM density.

The total training campaign therefore consists of five independent NSF runs: two channel-specific SM flows and three mediator-specific DM flows. Closure checks use the 50–100 validation band to verify the channel-specific SM models before SR score arrays are interpreted. The axial-vector \(\mu\mu\) score distribution shows a localised bin-wobble pattern near \(\mathcal{S}_h\approx0\)–50, consistent with residual NSF spline instability from non-convergence; this feature appears in all three fit variants and is therefore a property of the trained flow rather than an artifact of any single background model.

7 Likelihood-Ratio Search↩︎

The search uses \(\mathcal{S}_h\) as the per-event test statistic for each mediator hypothesis (Eq. 1 ). For \(\mu\mu\) SR events, vector, axial-vector, and scalar DM log densities are compared with NSF_SM_mumu; for \(ee\) SR events the same three DM log densities are compared with NSF_SM_ee. This yields six score arrays in total. We perform three simultaneous SR+VR binned profile-likelihood fits for signal strength \(\mu\), one for each DM hypothesis. Each fit combines the \(\mu\mu\) and \(ee\) score histograms from both the SR and VR simultaneously: \[\ln\mathcal{L} = \ln\mathcal{L}_{\mathrm{SR}} + \ln\mathcal{L}_{\mathrm{VR}} - \tfrac{1}{2}\theta_{\mu\mu}^2 - \tfrac{1}{2}\theta_{ee}^2,\] where \(\mathcal{L}_{\mathrm{SR}}\) and \(\mathcal{L}_{\mathrm{VR}}\) are binned Poisson likelihoods over the score histograms in each region, the signal strength \(\mu\) and the per-channel normalisation nuisances \(\theta_{\mu\mu}\), \(\theta_{ee}\) are shared between regions, and the background templates in both regions are derived from a single-step VR\(\to\)SR shape transfer: the SM VR score histogram is renormalised to the respective region yield and used as the nominal background prediction for that region. The observed data histograms are the raw SR and VR score counts, which differ from the renormalised VR background template by real Poisson fluctuations and by genuine shape differences between the VR and SR score distributions. The VR component constrains the background normalisation using \(\sim\)​6,970 (\(\mu\mu\)) and \(\sim\)​6,644 (\(ee\)) VR events independently of the SR, providing a genuine sideband-based background prediction. Signal MC contributes to both SR and VR components with the same \(\mu\). Asymptotic formulae [23] are used for the discovery test statistic and CL\(_s\) upper limits.

7.1 Systematic uncertainties↩︎

The profile-likelihood fit includes independent Gaussian-constrained normalisation nuisance parameters for the \(\mu\mu\) and \(ee\) channels with a prior width of 5% (\(\sigma=0.05\)). In the simultaneous SR+VR fit the normalisation nuisances are constrained additionally by the VR data; the posterior uncertainty on \(\theta_{\mathrm{norm}}\) is approximately 0.17 (compared to the prior of 1.0 in sigma units), corresponding to an effective normalisation constraint of \(\sim\)​0.85% from the VR statistics alone.

7.1.0.1 Integrated luminosity (propagated).

The Run 2015D dataset corresponds to \(\mathcal{L}=2.32\pm0.037\,\mathrm{fb}^{-1}\) (uncertainty 1.6%), taken from the CMS precision luminosity measurement [30]. A correlated log-normal luminosity nuisance is applied across both channels.

7.1.0.2 Lepton efficiency (propagated).

Combined identification, reconstruction, and trigger efficiency uncertainties of approximately 3% per channel are applied as independent per-channel log-normal nuisances, following the treatment in the CMS Run 2 mono-\(Z\) search at the same centre-of-mass energy [8]. Lepton momentum-scale uncertainties of 1% (\(\mu\)) and 2%/5% (electrons, barrel/endcap) are included as additional shape nuisances on the lepton \(p_{\mathrm{T}}\) distributions.

7.1.0.3 MET resolution and pileup (not propagated).

Evaluating \(E_{\mathrm{T}}^{\mathrm{miss}}\)resolution and pileup-reweighting uncertainties requires re-scoring SR, VR, and CR events with shifted inputs, which constitutes new analysis work not performed in this study. These sources are acknowledged as a gap in the systematic budget.

7.1.0.4 Combined systematic budget.

Treating luminosity as fully correlated across channels and lepton efficiency as per-channel independent, the combined propagated external systematic is approximately \(\sqrt{1.6^2+3^2+3^2}\approx4.5\%\) per channel, on top of the 0.85% data-driven normalisation constraint. The dominant uncertainty is the unresolved high-\(E_{\mathrm{T}}^{\mathrm{miss}}\)tail modelling residual (Section 7.2).

7.2 High-\(E_{\mathrm{T}}^{\mathrm{miss}}\)tail modelling residual↩︎

Events with \(E_{\mathrm{T}}^{\mathrm{miss}}\ge100\,\text{GeV}\) constitute 2.4–2.5% of the SR (\(\sim\)​160–181 events across both channels) but carry mean NSF scores 140–195 units higher than the VR bulk — a shift exceeding the score distribution’s own standard deviation of 40–60 units. The single-step VR background has negligible support at these scores, so the tail population appears as a systematic upward fluctuation in the fit.

To characterise this effect, we tested linear and quadratic \(E_{\mathrm{T}}^{\mathrm{miss}}\)-dependent models of the CR/VR score trend, validating each against held-out VR bins before extrapolating to the SR tail (Section 7.3). Neither functional form reliably predicts the observed tail behaviour: the linear model underpredicts the true SR-tail mean score by 90–115 units, while the quadratic model overpredicts by 300–330 units. The limited tail statistics (160–181 events across 100 ) do not support more flexible extrapolation.

The VR-shape background is therefore used as nominal, and the resulting inflation of \(\hat{\mu}\), \(q_0\), and the observed-to-expected limit ratios reflects this unresolved modelling residual rather than evidence of DM production.

7.3 NSF CR\(\to\)SR extrapolation validation↩︎

To validate that the NSF SM density, trained only on the control region (\(E_{\mathrm{T}}^{\mathrm{miss}}<50\,\text{GeV}\)), extrapolates in a characterised way into the signal region, we measured the mean NSF score as a function of \(E_{\mathrm{T}}^{\mathrm{miss}}\)across CR and VR bins for all three mediator hypotheses and both lepton channels (six combinations), then tested two extrapolation models against held-out data.

7.3.0.1 Procedure.

Events in the CR and VR were binned by \(E_{\mathrm{T}}^{\mathrm{miss}}\)in 10 GeV intervals. For each bin the mean NSF score \(\langle\mathcal{S}_h\rangle\) was computed. The two highest-\(E_{\mathrm{T}}^{\mathrm{miss}}\)VR bins (80–100 GeV) were held out as a validation set. Linear and quadratic functions of \(E_{\mathrm{T}}^{\mathrm{miss}}\)were fit to the remaining bins and extrapolated to the SR tail (\(E_{\mathrm{T}}^{\mathrm{miss}}\ge100\,\text{GeV}\)), where the true mean score was measured independently.

7.3.0.2 Results.

The linear fit undershot the held-out VR bins by 0.3–0.6\(\sigma\) (systematic across all six combinations). When extrapolated to \(E_{\mathrm{T}}^{\mathrm{miss}}\sim150\,\text{GeV}\) it overshot the true SR-tail mean score by 90–115 units. The quadratic fit overshot the held-out bins by 1.1–3.3\(\sigma\) (in the opposite direction from the linear case) and overshot the true SR-tail mean by 300–330 units; the positive curvature in all six fits indicates the quadratic form was dominated by the steep low-\(E_{\mathrm{T}}^{\mathrm{miss}}\)rise rather than the true behaviour at higher \(E_{\mathrm{T}}^{\mathrm{miss}}\).

7.3.0.3 Conclusion.

Neither functional form reliably extrapolates the CR-trained NSF density into the high-\(E_{\mathrm{T}}^{\mathrm{miss}}\)SR tail. Combined with the limited tail statistics (160–181 events across 100 ), this demonstrates that the NSF’s CR-trained density does not extrapolate smoothly or predictably into the sparse high-\(E_{\mathrm{T}}^{\mathrm{miss}}\)signal region. This validates and motivates the treatment of the resulting background modelling residual as an explicit, quantified limitation (Section 7.2) rather than an uncorrected prediction.

8 Results↩︎

The normalizing-flow likelihood-ratio discriminant was evaluated using disjoint control-region (CR), validation-region (VR), and signal-region (SR) samples. The Standard Model density flows were trained only on CR events with \(E_{\mathrm{T}}^{\mathrm{miss}}<50\,\text{GeV}\). The final statistical interpretation uses the full SR selection, \(E_{\mathrm{T}}^{\mathrm{miss}}\ge50\,\text{GeV}\), \(|\Delta\phi(E_{\mathrm{T}}^{\mathrm{miss}},Z)|>2.5\), and \(n_{\mathrm{jets}}\le1\). Although no explicit upper \(E_{\mathrm{T}}^{\mathrm{miss}}\)bound is imposed in the selection, the cleaned dataset caps \(E_{\mathrm{T}}^{\mathrm{miss}}\)at approximately 200 , and the VR interval (\(50\leE_{\mathrm{T}}^{\mathrm{miss}}<100\,\text{GeV}\)) is retained as a validation and stress-test sample rather than as the nominal background template for the final limits.

For each mediator hypothesis the NSF score \[S_h(\mathbf{x}) = \log p(\mathbf{x}\mid\mathrm{DM}_h) - \log p(\mathbf{x}\mid\mathrm{SM}_{\ell\ell})\]

was constructed separately for the vector, axial-vector, and scalar mediator hypotheses. Large positive values correspond to signal-like events, while negative values indicate SM-like kinematics.

Control-region closure tests show stable channel-specific behaviour of the trained SM flows. The CR distributions of \(\log p(\mathbf{x}\mid\mathrm{SM})\) peak near \(\sim 125\) with similar widths in the \(\mu\mu\) and \(ee\) channels, consistent with a well-converged density model on Drell–Yan-dominated events.

Validation-region score distributions show clear separation between observed data and simulated dark matter events for all mediator hypotheses. The strongest discrimination is observed for the scalar mediator, whose signal distribution occupies a largely distinct region of score space. Consistent behaviour is observed in both lepton channels, demonstrating that the NSF score generalises across detector signatures and event topologies.

The SR score distributions are shown in Fig. 4. In all cases, the observed data are concentrated at lower score values, whereas simulated dark matter events populate the high-score tail. The scalar hypothesis produces the largest signal-background separation, followed by the vector and axial-vector models.

A binned profile-likelihood fit was performed using the NSF score distributions. Initial histograms were constructed using 30 score bins and subsequently merged to ensure a minimum nominal background yield of 20 events per bin. The nominal background shapes are derived from a single-step VR\(\to\)SR transfer, while the final statistical interpretation uses a simultaneous SR+VR profile-likelihood fit with independent 5% normalization nuisance parameters in the \(\mu\mu\) and \(ee\) channels. The validation region acts as a normalisation sideband and constrains the background prediction through nuisance-parameter profiling, providing an independent cross-check of the SR-only result shown in Appendix 14.5. A pure VR-extrapolated shape alternative was also investigated but failed closure in the high-score tail; the corresponding VR-only post-fit distributions are provided in Appendix 14.4 as a diagnostic study and are not used for the final statistical interpretation.

The fitted score distributions for all mediator hypotheses and both lepton channels are shown in Fig. 3. In every case the observed data show a pull relative to the VR-based background prediction, most pronounced in the high-score tail (Section 7.2). The discovery test statistic \(q_0\) is large (233–327 across the three hypotheses; Table 2), formally corresponding to a significance capped at \(Z=8.0\) by floating-point underflow; we attribute this to the characterized high-\(E_{\mathrm{T}}^{\mathrm{miss}}\)background-modelling residual rather than a dark matter signal (Sections 7.2–7.3).

Upper limits on the signal strength parameter \(\mu\) were derived using the asymptotic CL\(_s\) method. The resulting limits are summarized in Table 2, while the corresponding confidence-level scans are shown in Fig. 5. The scalar mediator provides the strongest exclusion, yielding an observed 95% confidence-level upper limit of \(\mu < 0.0177\). The vector and axial-vector hypotheses yield observed limits of \(\mu < 0.0362\) and \(\mu < 0.0498\), respectively.

For all three mediator hypotheses the observed limits exceed the median expected limits by factors of roughly 7–12 across the nominal and cross-check background constructions, reflecting the unresolved high-\(E_{\mathrm{T}}^{\mathrm{miss}}\)tail modelling residual described in Section 7.2. The superior performance of the scalar model is consistent with its significantly stronger separation in NSF score space. The axial-vector hypothesis exhibits the weakest sensitivity owing to the greater overlap between signal and background score distributions.

9 Summary and Outlook↩︎

We presented a mono-\(Z\) DM search that replaces hand-crafted cut flows with likelihood-ratio scores from Neural Spline Flows trained on CR data and DM MC. The workflow—open-data extraction, audited cleaning, disjoint CR/SR design, and three combined \(\mu\mu\)+ee fits from five trained flows—provides a reproducible template for density-based searches at the LHC. The final simultaneous SR+VR fits yield observed 95% CL limits of \(\mu<0.0177\) (scalar), \(\mu<0.0362\) (vector), and \(\mu<0.0498\) (axial-vector), corresponding to the cross-section limits reported in Table 2. The observed-to-expected gaps are driven by the high-\(E_{\mathrm{T}}^{\mathrm{miss}}\)tail modelling residual rather than by a dark matter signal. A pure VR-extrapolated shape alternative fails closure in the high-score tail (Appendix 14.4) and is not used for the final limits.

Appendix↩︎

10 Dataset and Extraction Summary↩︎

10.1 Rows read and extracted↩︎

Table 3 summarizes the row counts used after extraction. For Run 2015D data, raw rows read are aggregated from the per-file audit CSVs under Data/Extracted. For DM MC, the source files are ROOT files; the Rows read column reports the total raw entries processed before extraction, and the extracted rows are the selected event counts in the output TTrees and EDA summaries.

10.2 HLT trigger paths and b-tag configuration↩︎

Table 4 summarizes the trigger and \(b\)-tag settings recorded by the extraction notebooks. Trigger decisions are read from selectedPatTrigger path-name and final-filter collections when available; the filter tokens are used as the fallback audit path for 2015 MiniAOD files with empty path-name payloads. The data extractors apply the medium CSVv2 \(b\)-veto, while the DM MC extractor retains n_bjets only for schema compatibility.

10.3 Signal MC parameter points↩︎

Table 5 lists the signal samples used for the three DM density models. The parameter values are taken from the CERN Open Data sample names in the local bibliography.

11 Feature Definitions↩︎

11.1 Physics feature branches↩︎

The extraction notebooks write the 40 physics branches listed in Tables 6 and 7 to the output Events TTrees. The DoubleMuon, DoubleEG, and MonoZToLL DM extractors use the same branch set; data outputs add audit/provenance branches, while DM MC adds gen_weight. We use \(\Delta\phi(a,b)=|((a-b+\pi)\bmod 2\pi)-\pi|\), \(m_Z=91.1876\,\text{GeV}\), and \(\epsilon=10^{-6}\) for protected divisions.

11.2 NSF input feature order↩︎

Table 8 gives the 37-dimensional input order used when loading the trained NSF checkpoints and scoring events. This order is read from Output/results/scoring_audit.json; has_jet1 and has_jet2 are computed from n_jets before clipping.

12 Data Cleaning↩︎

12.1 Fixed domain bounds↩︎

Table 9 lists the fixed-domain bounds applied by the SM cleaner. Features not shown in the fixed-bound table are either pass-through variables or receive only the tail-rule treatment in Table 11.

12.2 Data-derived tail caps↩︎

Table 10 lists the numerical caps stored in Output/sm_feature_bounds.json for variables governed by data-derived tail rules. Positive rules clip to \([0,\mathrm{cap}]\), while symmetric rules clip to \([-\mathrm{cap},+\mathrm{cap}]\).

12.3 Tail-rule families↩︎

13 Region Definitions and Event Yields↩︎

13.1 Region selection criteria↩︎

Table 12 records the CR, VR, and SR masks from Output/results/scoring_audit.json. The VR is nested in the SR kinematic definition and is used as a 50–100 closure and stress-test region.

13.2 Event counts per region and channel↩︎

Table 13 summarizes the region counts stored in the scoring audit. DM MC is scored in both SM channel scaler spaces, so the mediator-specific VR and SR counts are identical for the \(\mu\mu\) and \(ee\) scoring passes.

14 Fit Configuration and Numerical Results↩︎

14.1 Fit configuration↩︎

Table 14 lists the final simultaneous SR+VR fit settings recorded in paper/json/sr+vr/fit_summary.json and associated output files. The fit combines the \(\mu\mu\) and \(ee\) score histograms for each mediator hypothesis.

14.2 Full fit results↩︎

Table 15 reports the numerical results from paper/json/sr+vr/fit_summary.json. All three signal-plus-background fits converged; the large \(q_0\) values are attributed to the high-\(E_{\mathrm{T}}^{\mathrm{miss}}\)background-modelling residual discussed in Section 7.2.

14.3 Background estimation method note↩︎

The simultaneous SR+VR fit uses a single-step VR\(\to\)SR shape transfer: the SM VR score histogram is renormalised to the region yield and used as the nominal background template for that region. The VR component of the likelihood constrains the per-channel normalisation nuisances \(\theta_{\rm norm}\) using \(\sim\)​6,970 (\(\mu\mu\)) and \(\sim\)​6,644 (\(ee\)) VR events; the posterior uncertainty on \(\theta_{\rm norm}\) is \(\sim\)​0.17 in sigma units, significantly tighter than the 5% Gaussian prior alone. Signal MC contributes to both SR and VR with the same signal strength \(\mu\), so the VR signal contamination is accounted for. This approach avoids a pure SR-direct circularity while keeping the SR score shape as the primary discriminant. Bins with zero nominal background are excluded from the likelihood.

14.4 VR-extrapolated validation fit↩︎

As a sideband stress test, an alternative fit used the SM VR (\(50\leE_{\mathrm{T}}^{\mathrm{miss}}<100\,\text{GeV}\)) score histograms, normalized to the SR yield, as the nominal background template. This procedure does not close in the high-score tail: the VR template underpredicts the SR tail, and the signal-plus-background fit absorbs the resulting shape difference as a positive signal strength. Table 16 therefore records the VR-extrapolated fit as a diagnostic rather than as the final limit result.

With a signal-strength scan extended to bracket \(\hat{\mu}\) for all three mediators, the VR-extrapolated fit yields observed limits 9.5–10.8\(\times\) weaker than expected, consistent in direction and magnitude with the nominal fit’s 7.2–9.8\(\times\) gap (Table 2). This corroborates that the high-\(E_{\mathrm{T}}^{\mathrm{miss}}\)tail residual, not a background-construction artifact specific to one method, drives the observed/expected discrepancy.

14.5 SR-only validation fit↩︎

As an additional robustness check, the likelihood fit was repeated using only the SR score distributions without the VR normalisation constraint. The fitted signal strengths remain positive with capped \(Z=8.0\), as expected from the same high-\(E_{\mathrm{T}}^{\mathrm{miss}}\)residual, while the resulting limits are compatible with those obtained from the nominal simultaneous SR+VR fit. The corresponding CL\(_s\) scans and post-fit distributions are shown in Figs. 8 and 9.

15 Additional Validation Figures↩︎

15.1 VR-extrapolated CL\(_s\) scans↩︎

15.2 VR-extrapolated post-fit distributions↩︎

15.3 SR-only CL\(_s\) scans↩︎

15.4 SR-only post-fit distributions↩︎

References↩︎