Search for dark matter produced in association with a single top quark or a top quark pair in proton-proton collisions at $\sqrt{s} =$ 13 TeV

A search for dark matter produced in association with top quarks in proton-proton collisions at a center-of-mass energy of 13 TeV is presented. The data set used corresponds to an integrated luminosity of 35.9 fb$^{-1}$ recorded with the CMS detector at the LHC. Whereas previous searches for neutral scalar or pseudoscalar mediators considered dark matter production in association with a top quark pair only, this analysis also includes production modes with a single top quark. The results are derived from the combination of multiple selection categories that are defined to target either the single top quark or the top quark pair signature. No significant deviations with respect to the standard model predictions are observed. The results are interpreted in the context of a simplified model in which a scalar or pseudoscalar mediator particle couples to a top quark and subsequently decays into dark matter particles. Scalar and pseudoscalar mediator particles with masses below 290 and 300 GeV, respectively, are excluded at 95% confidence level, assuming a dark matter particle mass of 1 GeV and mediator couplings to fermions and dark matter particles equal to unity.


Introduction
Astrophysical observations provide evidence of the existence of nonluminous matter that can be inferred from gravitational effects on galaxies and other large scale objects in the Universe. While the nature of this dark matter (DM) is still unknown, a compelling candidate is the socalled weakly interacting massive particle [1]. This new particle is predicted to have weak interactions with standard model (SM) particles, allowing for direct-and indirect-detection experiments, as well as for searches at collider experiments.
Among all the possible interactions between the SM and DM sectors, it is of particular interest to investigate interactions mediated by a new neutral scalar or pseudoscalar particle that decays into DM particles, as these can be easily accommodated in models containing extended Higgs boson sectors [2][3][4][5]. Assuming that this DM scenario respects the principle of minimal flavor violation [6,7], the interactions of this new spin-0 mediator particle follow the same Yukawa coupling structure as in the SM. Therefore, the mediator would couple preferentially to heavy third-generation quarks. Assuming the DM particles to be Dirac fermions, the interaction Lagrangian terms for the production of a scalar (φ) or pseudoscalar (a) mediator particle can be expressed as: where the sum runs over the SM fermions f, y f = √ 2m f /v represents the Yukawa couplings, v = 246 GeV is the Higgs field vacuum expectation value, g χ is the DM-mediator coupling, and g q is the fermion-mediator coupling. The mediator particle subsequently decays into DM particles, which escape detection and leave an imbalance of momentum in the transverse plane, referred to as p miss T . Several theoretical studies of these types of models have been performed, in which the third-generation quark is either a top or bottom quark, leading to the production of DM in association with a pair of top (tt+DM) or bottom (bb+DM) quarks, respectively [8][9][10][11]. The main production diagram for tt+DM processes is shown in Fig. 1

(upper left).
Previous searches in these final states have been carried out by the ATLAS and CMS Collaborations at center-of-mass energies of 8 TeV [12, 13] and 13 TeV [14][15][16]. While the former results are based on an effective field theory (EFT) approach, the latter ones are interpreted in the context of simplified DM scenarios, where the mediator particle is explicitly modeled in the interaction. These interpretations have so far neglected the contribution from DM production in association with a single top quark (t/t+DM) in which the interaction is mediated by a neutral spin-0 particle, as pointed out in Ref. [17]. As in the SM, the single top quark is produced through processes mediated by a virtual t channel (Fig. 1, upper right) or through associated production with a W boson (Fig. 1, lower left and right) [17]. While the s channel production of a W boson is also possible, this process is found to have a negligible contribution for this search. The neutral DM mediator particle is then produced either by radiation from the top quark or via top quark fusion, as described in Ref. [18] for the associated production of DM with a top quark pair.
In this search, t/t+DM processes mediated by a neutral spin-0 particle are investigated for the first time. This additional production mechanism is predicted by the same interactions described in Eqs. (1) and (2) that also predict tt+DM events. For this reason, in the presented search t/t+DM and tt+DM processes are both considered. Searches for similar final states referred to as "monotop", which involve the production of a top quark and DM particles but . The underlying simplified models explored in these results, unlike the one presented in Eqs. (1) and (2), assume either the resonant production of a +2/3 charged and colored spin-0 boson that decays into a right-handed top quark and one DM particle, or a spin-1 mediator with flavor changing neutral current interactions. Considering these models, in addition to the DM particle, only one top quark is assumed to be produced in the final state, unlike the t/t+DM processes considered in this search where the top quark is produced through SM-like diagrams alongside a light quark or a W boson (Fig. 1).
In this paper we present a search for an excess of events above the SM background in the p miss T spectrum, as expected for the DM scenarios discussed earlier, for events that contain exactly one lepton (electron or muon) or zero leptons, henceforth assigned to the "single-lepton" (SL) region or to the "all-hadronic" (AH) region, respectively. The sensitivity of this analysis is improved beyond that of previous analyses by introducing a categorization of these signatures and new discriminating variables, as discussed in more detail in Section 4.

The CMS detector and event reconstruction
The central feature of the CMS apparatus is a superconducting solenoid of 6 m internal diameter, providing a magnetic field of 3.8 T. Within the solenoid volume are a silicon pixel and strip tracker, a lead tungstate crystal electromagnetic calorimeter (ECAL), and a brass and scintillator hadron calorimeter (HCAL), each composed of a barrel and two endcap sections. Forward calorimeters extend the pseudorapidity (η) coverage provided by the barrel and endcap detectors. Muons are detected in gas-ionization chambers embedded in the steel flux-return yoke outside the solenoid. A more detailed description of the CMS detector, together with a definition of the coordinate system used and the relevant kinematic variables, can be found in Ref. [23].
Events of interest are selected using a two-tiered trigger system [24]. The first level, composed of custom hardware processors, uses information from the calorimeters and muon detectors to select events at a rate of around 100 kHz within a time interval of less than 4 µs. The second level, known as the high-level trigger, consists of a farm of processors running a version of the full event reconstruction software optimized for fast processing, and reduces the event rate to around 1 kHz before data storage.
The particle-flow (PF) algorithm [25] aims to reconstruct and identify each individual particle in an event, with an optimized combination of information from the various elements of the CMS detector. The energy of photons is obtained directly from the ECAL measurement and corrected for zero-suppression effects. The energy of electrons is obtained from a combination of the electron momentum at the primary interaction vertex as determined by the tracker, the energy of the corresponding ECAL cluster, and the energy sum of all bremsstrahlung photons spatially compatible with originating from the electron track. The muon track is obtained from the combination of central tracker and muon system information, and its curvature provides an estimate of the momentum. The energy of charged hadrons is determined from a combination of their momentum measured in the tracker and the matching ECAL and HCAL energy deposits, corrected for zero-suppression effects and for the response function of the calorimeters to hadronic showers. Finally, the energy of neutral hadrons is obtained from the corresponding corrected ECAL and HCAL energy.
The reconstructed vertex with the largest value of summed physics-object p 2 T , where p T is the transverse momentum, is taken to be the primary proton-proton (pp) interaction vertex. The physics objects are the jets and the associated p miss T , taken as the negative vector p T sum of those jets. For each event, hadronic jets are clustered from the particles reconstructed with PF (PF candidates) using the infrared-and collinear-safe anti-k T algorithm [26,27] with a distance parameter of 0.4. The jet momentum is determined as the vectorial sum of all particle momenta in the jet, and is found from simulation to be within 5 to 10% of the parton's generated momentum over the whole p T spectrum and detector acceptance. Additional pp interactions within the same or nearby bunch crossings (pileup) can contribute additional tracks and calorimetric energy depositions to the jet momentum. To mitigate this effect, tracks identified as originating from pileup vertices are discarded and an offset correction is applied to correct for remaining contributions [28]. Jet energy corrections are derived from simulation and applied to calibrate the jet momentum. In situ measurements of the momentum balance in dijet, pho-ton+jet, Z+jets, and multijet events are used to account for any residual differences in jet energy scale in data and simulation [29]. Additional selection criteria are applied to each jet to remove jets potentially dominated by anomalous contributions from various subdetector components or reconstruction failures [29].
The combined secondary vertex b tagging algorithm (CSVv2) is used to identify jets originating from the hadronization of bottom quarks [30], denoted in the following as "b-tagged jets". At the operating point of the tagging algorithm chosen for this analysis, the efficiency of identifying b quark jets in simulated tt events is about 80%, integrated over p T , and the misidentification rate for light-flavor jets is about 1%. Scale factors are applied to the simulated samples in order to reproduce the b tagging performance measured in data.
The missing transverse momentum vector p miss T is defined as the negative vector p T sum of all PF particles originating from the primary vertex; its magnitude is defined as p miss T . Jet energy scale and resolution corrections are also propagated to the p miss T calculation.

Data sample and simulation
The data used in this search were recorded with the CMS detector in 2016 and correspond to an integrated luminosity of 35.9 fb −1 . Several trigger criteria were used to collect the data, either requiring large amounts of p miss T or the presence of at least one high-p T lepton (electron or muon). Simulated samples are corrected to reproduce the observed trigger efficiencies in data.
Specifically, events that do not contain leptons are selected if they have p miss T and missing hadronic activity H miss T [24] above 120 GeV. This trigger is nearly 100% efficient for events with p miss T of at least 250 GeV. The second set of triggers requires the presence of at least one isolated electron (muon) with p T > 27 (25) GeV. The corresponding trigger efficiencies are above 90% for leptons with p T > 30 GeV. Trigger efficiencies are measured in data.
Monte Carlo (MC) simulated samples of the main SM backgrounds and of the DM signal processes are used to optimize the event selection, assess our sensitivity to the new-physics scenarios, and form the basis of our background estimation strategy. While the detailed background composition depends on the specific channel, the main sources arise from tt+jets, W+jets, and Z+jets production. Simulated events of tt+jets production and single top quark processes are generated at next-to-leading order (NLO) in quantum chromodynamics (QCD) using POWHEG v2 and POWHEG v1 [31][32][33], respectively. For tt+jets processes, the top quark p T distribution is reweighted to reproduce the differential cross section obtained from CMS measurements [34]. Samples of Z+jets, W+jets, and QCD multijet events are generated at leading order (LO) using MADGRAPH5 aMC@NLO [35] with the MLM prescription [36] for matching jets from the matrix element (ME) calculation to the parton shower description. Dedicated electroweak [37][38][39][40][41][42] and QCD (calculated with MADGRAPH5 aMC@NLO) NLO/LO K factors, parametrized as functions of the generated boson p T , are applied to Z+jets and W+jets events. Other SM backgrounds include rare processes, such as tt+W and tt+Z, which are simulated based on the NLO ME calculations implemented in MADGRAPH5 aMC@NLO and the FxFx [43] prescription to merge multileg processes. Diboson processes (WW, WZ, ZZ, WH, ZH) are generated at NLO using either MADGRAPH5 aMC@NLO or POWHEG v2. All background samples are normalized using the most accurate cross section calculations available, which generally incorporate NLO or next-to-NLO (NNLO) precision.
The signal process is simulated at LO with the MADGRAPH5 aMC@NLO v2.4.2 event generator using a simplified model investigated within the LHC Dark Matter Forum [44]. In this model, the DM particles χ are assumed to be Dirac fermions and the mediators are spin-0 particles φ (a) that couple preferentially to third-generation SM quarks through scalar (pseudoscalar) couplings whose strengths are parametrized by the factor g q . The coupling strength between the mediator and the DM particles is in turn given by the factor g χ . This simplified model has a minimal set of four free parameters: (m χ , m φ/a , g χ , g q ), and the benchmark scenarios assume g χ = g q = 1 as per recommendations of the LHC Dark Matter Working Group [45]. In addition, in this search we focus on the m χ = 1 GeV benchmark, which is a convenient signal reference as the production cross section is almost independent of m χ for on-shell mediators [44]. This simplified spin-0 model does not account for mixing between the φ scalar mediator and the SM Higgs boson, as discussed in Ref. [46]. Under these assumptions two distinct DM scenarios are possible: the associated production with a top quark pair (tt+DM) and the associated production with a single top quark (t/t+DM). Cross sections for both signal processes are calculated at LO with MADGRAPH5 aMC@NLO v2.4.2, with one (zero) additional partons for tt+DM (t/t+DM) events.
For all simulated samples, the initial-state partons are modeled with the NNPDF 3.0 [47] par-ton distribution function (PDF) sets at LO or NLO in QCD to match the ME calculation. Generated events are interfaced with PYTHIA 8.205 [48] for parton showering and hadronization using the CUETP8M1 tune [49], except for simulated tt+jets events where the CUETP8M2 tune customized by CMS with an updated strong coupling α S for initial-state radiation is employed [50]. All signal and background samples are processed using GEANT4 [51] to provide a full simulation of the CMS detector, including a simulation of the previously mentioned triggers. Correction factors are derived and applied to the simulated samples to match the trigger efficiencies measured in data. Additional corrections are applied to cover remaining residual differences between data and simulation that arise from the lepton identification and reconstruction efficiencies, as well as from b-tagged jet identification efficiencies.

Event selection
This search, similarly to a previous search for tt+DM events [16], defines several orthogonal signal regions (SRs) that are statistically combined in a simultaneous global fit of the p miss T spectrum. At the same time, various improvements are incorporated into this search to enhance the sensitivity to the t/t+DM final state over that of previous analyses [16].
At the analysis level, jet candidates are required to have p T > 30 GeV and are categorized as "central" if they lie within |η| < 2.4 and as "forward" if they are within 2.4 < |η| < 4.0. The b-tagged jets identified by the CSVv2 algorithm are also required to have p T > 30 GeV and in addition to lie within |η| < 2.4. Electrons and muons are selected with p T > 30 GeV and |η| < 2.1. Events containing additional leptons with p T > 10 GeV and |η| < 2.1 are vetoed. To ensure that candidate leptons are well-measured, identification requirements, based on hit information in the tracker and muon systems and on energy deposits in the calorimeters, are imposed. Leptons are further required to be isolated from hadronic activity, to reject leptons within jets that could arise, for example, from the decay of b quarks. A relative isolation quantity is defined as the scalar p T sum of all PF candidates within a ∆R = √ (∆η) 2 + (∆φ) 2 cone of radius 0.3 (0.4) centered around the electron (muon) candidate, where φ is the azimuthal angle in radians, divided by the lepton p T [52, 53]. This relative isolation is required to be less than 0.059 (0.057) for electrons in the barrel (endcap) and less than 0.15 for muons.
Events are separated into orthogonal categories based on the number of b-tagged jets (n b ), with n b = 1 or n b ≥ 2, and additional requirements on the number of forward jets are placed (0 or ≥1 forward jets) for the n b = 1 category. The mentioned categorization in terms of forward jets allows a further enhancement of t/t+DM t channel events. In fact, as shown in Fig. 1, this production mode leads to final states with one top quark and an additional jet, which tends to be in the forward region of the detector, while the additionally produced b quark is typically low in p T and therefore is not reconstructed. The minimum requirements on the number of jets is also lowered, with respect to the previous searches, to enhance the sensitivity specifically to the t/t+DM model. Control regions (CRs) enriched in the major background processes are included in the fit in order to improve the estimates of the background contributions.
Events are classified into two "channels", based on the number of leptons in the final state from the top quark decay: the single-lepton SL channel, containing events with exactly one electron or muon with p T > 30 GeV, and the all-hadronic AH channel, containing events with exactly zero leptons with p T > 10 GeV. A set of discriminating variables is identified, as discussed in more detail in Sections 4.1 and 4.2 for the SL SRs and the AH SRs, respectively. The selection requirements on these variables are optimized simultaneously to increase the signal significance, using as a figure of merit the ratio between the expected number of signal and the square root of the expected SM background events. The considered signal events are either t/t+DM events for a region that contains exactly one b-tagged jet (n b = 1) or tt+DM events for a region that contains two or more b-tagged jets (n b ≥ 2). The region with exactly one b-tagged jet is further divided into exactly zero or ≥1 forward jets.

Single-lepton signal regions
Events in the SL channel are required to contain ≥1 identified b-tagged jet, at least 2 jets with p T > 30 GeV, and p miss T > 160 GeV. After this selection, the dominant backgrounds in the SL channel are from tt and W+jets processes. Other backgrounds include single top quark, Drell-Yan, and diboson production.
To further improve the sensitivity and to reduce the dominant background from single-lepton tt and W+jets processes, we impose a requirement on the transverse mass m T , calculated as: where p T is the transverse momentum of the lepton and ∆φ is the opening angle between the lepton direction and the p miss T vector in the transverse plane. The m T variable is constrained by kinematic properties to be less than the W boson mass for leptonic on-shell W decays in tt and W+jets events, while for signals, off-shell W decays, or for dileptonic decays of tt, the m T variable is expected to exceed the W mass because of the additional p miss T in the event. A requirement of m T > 160 GeV therefore reduces the background from single-lepton events significantly and enhances the analysis sensitivity to the DM models.
After the m T selection, the remaining tt background is primarily from events where both top quarks decay leptonically (tt(2 )) and one lepton is not identified. This background can be further reduced by making use of the m W T2 variable [54], which is defined as the minimal value of the mass of a particle assumed to be pair produced and to decay to a W boson and a b quark jet. The W bosons are assumed to be produced on-shell and to decay leptonically, where one of the two leptons is not detected. Based on the variable definition, in tt(2 ) events the m W T2 distribution has a kinematic end point at the top quark mass, assuming perfect detector response, while this is not the case for signal events where two additional DM particles are present. The calculation of m W T2 requires two b-tagged jets from the decay of the top quarks, where one of these b-tagged jets comes from the same decay chain as the reconstructed lepton. If only one b-tagged jet is identified in the event, each of the first three (or two in three-jet events) leading non-b-tagged jets is considered as the second b-tagged jet in the calculation. The m W T2 is then evaluated for all possible jet-lepton combinations and the minimum m W T2 value is considered to discriminate between signal and background events. If two or more b-tagged jets are identified in the events, all b-tagged jets are considered and similarly all possible jetlepton combinations are used to calculate m W T2 values. The smallest of all the m W T2 values is taken as the event discriminant.
In addition, jets and the p miss T vector tend to be more separated in the transverse plane in signal events than in tt background processes. To improve the search sensitivity, the minimum opening angle min∆φ(j 1,2 , p miss T ) in the transverse plane between the direction of each of the first two leading-p T jets with |η| < 2.4 and the p miss T vector is required to be greater than 1.2 radians.
The tt background is further reduced by requiring that the transverse mass m b T of the p miss T vector and of a b-tagged jet is greater than 180 GeV, where m b T is defined similarly to Eq. (3) but considering a b-tagged jet instead of a lepton. In fact, for the remaining tt background m b T tends to have values below or around the top quark mass if the b-tagged jet belongs to the top quark whose lepton is not identified. For the calculation we choose the b-tagged jet with the highest CSVv2 discriminant value, if there is more than one candidate.
A summary of the selection criteria for the SL SRs is shown in the first three columns of Table 1. Each region is identified by a unique name, where 0 denotes exactly zero leptons, 1(2) b-tag represents exactly 1 (at least 2) b-tagged jet, and 0 FJ or 1 FJ denotes exactly zero or at least one forward jet. Single-lepton SRs All-hadronic SRs

All-hadronic signal regions
Events categorized into the AH channel must contain at least 1 identified b-tagged jet and at least 3 jets with p T > 30 GeV, p miss T > 250 GeV, and min∆φ(j 1,2 , p miss T ) greater than 0.4 radians.
The dominant backgrounds after this selection arise from tt, W+jets, and Z → νν processes.
Other backgrounds include QCD multijet events, single top quark, Drell-Yan, and diboson production.
Semileptonic tt events populate this channel if the lepton in the final state is not identified. This tt(1 ) background is reduced by applying the same m b T selection as introduced in the SL channel. To further reduce the tt(1 ) background, together with that from Z → νν events, we make use of the p T (j 1 )/H T variable, which is defined as the ratio of the leading p T jet in the event divided by the total hadronic transverse energy in the event, H T , which is the scalar p T sum of the jets with p T > 30 GeV within |η| < 2.4. In the case of background, the distribution peaks at higher values with respect to tt+DM signal events. The t/t+DM events, however, tend to exhibit a distribution similar to that of the background. Events in the n b ≥ 2 category are required to have p T (j 1 )/H T < 0.5.
For QCD multijet events no intrinsic p miss T is expected. Therefore, events that pass our minimum p miss T selection contain mostly p miss T which arises from jet mismeasurements. For these events, the p miss T is often aligned with one of the leading jets. As a result, selecting events with min∆φ(j 1,2 , p miss T ) values greater than 1 radian reduces the background from QCD multijet production. This contribution to the SR, estimated through simulated samples, is negligible. The description of the QCD multijet background basic kinematic distributions is verified in a dedicated region enriched in multijet events, obtained by reversing the min∆φ(j 1,2 , p miss T ) selection, and the simulation is found to model the data well.
A summary of the selection criteria for the AH SRs is shown in the last three columns of Table 1. Each region is identified by a unique name, where 1 denotes exactly one muon or one electron, 1(2) b-tag represents exactly 1 (at least 2) b-tagged jet, and 0 FJ or 1 FJ denotes exactly zero or at least one forward jet.

Control regions
After events are categorized according to the selection presented in Table 1, the expected SM backgrounds in these different regions must be evaluated. In the SL SRs, the main backgrounds are dileptonic tt events, where one lepton is not identified, and W+jets events. For the AH regions the main backgrounds arise instead from single-lepton tt and W+jets events, where the lepton is not identified, and Z boson production, where the Z boson decays into two neutrinos and leads to a background with genuine p miss T . In order to improve the estimation of these main backgrounds, methods based on control samples in data are used. In particular, CRs enhanced in the different background sources are used to derive correction factors as a function of the p miss T from the comparison of the p miss T distribution between the data and the simulation. These corrections are extracted and simultaneously propagated across the CRs and SRs for a given channel in the context of a global fit, as explained in more detail in Section 6. The residual backgrounds processes are modeled with simulation.
The background CRs for the SL and AH channels are designed to be statistically independent from the corresponding SRs.

Single-lepton control regions
The first set of CRs is defined to isolate dileptonic tt events by requiring exactly two leptons (1 electron and 1 muon, 2 electrons, or 2 muons), n jet ≥ 2, n b ≥ 1, and p miss T > 160 GeV. In order to statistically enhance these CRs the m T , m W T2 , and forward jet selections are removed. The second set of CRs is designed to isolate W+jets events by requiring exactly one lepton (electron or muon), n jet ≥ 2, n b = 0, p miss T > 160 GeV, and m T > 160 GeV. The n b = 0 requirement makes this CR orthogonal to the SL SR and allows the events in the m T tail to be modeled without extrapolation from a lower-m T region.
Both of these selections are summarized in the first two columns of Table 2. Table 2: Control regions defined for the main backgrounds of the SL SRs (first two columns, tt(2 ) and W+jets) and the AH SRs (last 3 columns, tt(1 ), W+jets, and Z → ). Some selections applied in the SRs are removed in the corresponding CRs to increase the available statistics and are therefore not listed. The p miss T selection for the Z → CR refers to the hadronic recoil.

Single-lepton CRs
All-hadronic CRs

All-hadronic control regions
For the AH SRs, three independent sets of CRs are defined. The first set of CRs is enhanced in single-lepton tt events selecting events with exactly one lepton (electron or muon), n jet ≥ 3, n b ≥ 1, p miss T > 250 GeV, and, in order to avoid overlap with the SL SRs, m T < 160 GeV.
The second set of CRs is defined to enhance single-lepton W+jets events. Events are selected with exactly one lepton (electron or muon), n jet ≥ 3, n b = 0, p miss T > 250 GeV, and in order to avoid overlap with the SL W+jets CR, m T < 160 GeV.
The third and last set of CRs are designed to model the background due to Z+jets production, where the Z boson decays into a pair of neutrinos (Z → νν). Here we use the Z boson decays to an opposite-sign, same-flavor dilepton pair (Z → ), as proxy events to emulate the kinematic properties of the Z+jets process. Events are selected requiring 2 leptons, which have the same flavor (i.e., ee or µµ), and opposite charge, and that satisfy a requirement on their invariant mass of 60 < m < 120 GeV. Additionally, events must contain at least 3 jets, but events with b-tagged jets are vetoed (n b = 0). In order to reproduce the p T spectrum of Z → νν events, the two leptons are added to the p miss T , referred to as hadronic recoil.
A summary of the different AH CRs can be found in the last three columns of Table 2.

Systematic uncertainties
Several sources of uncertainty are considered that affect either the simulation of the background processes or the underlying theoretical modeling. We distinguish between two types of uncertainties, ones that only affect the normalization of a process and others that additionally affect the shape of the p miss T distribution. All of these uncertainties are included in the global simultaneous fit, described in detail later. The largest impacts on the final results stem from the uncertainties in the b tagging scale factors and the limited statistical precision of the dilepton tt CR, where the latter is the main determining factor for the contribution of tt events in the SL SRs.
The following sources of uncertainty correspond to constrained normalization nuisance parameters in the fit (unless specified, the source of uncertainty applies to all search channels): • Lepton reconstruction, selection, and trigger. Scale factors are applied to the simulation in order to mimic the measured lepton reconstruction and selection efficiencies in data. The measured uncertainties in these scale factors are of the order of 2.2% per electron and 1% per muon, and are p T and η dependent [52, 53]. The effect of these uncertainties is found to be independent of the p miss T spectrum.
• p miss T trigger. At values of p miss T > 250 GeV the applied triggers are almost fully efficient; a normalization uncertainty of 2% is assigned. This uncertainty is only applied in the AH channel.
• b tagging efficiency scale factors. The b tagging and light-flavor mistag efficiencies scale factors and the respective uncertainties are measured in independent control samples [30], and propagated to the analysis. In the range of p miss T considered, these scale factor uncertainties do not alter the shape of the p miss T distribution.
• Forward jets. Inclusive CRs in terms of forward jet multiplicity are considered to constrain the major background in the 0 and ≥1 forward jets SRs. The impact of this extrapolation in forward jets multiplicity on the background estimation is evaluated and assigned as additional systematic uncertainty. The extrapolation effect is evaluated by splitting each CR into a 0f and 1f category, and a systematic uncertainty is assigned based on the ratio of the correction factors, where each correction factor is the ratio of the data to the simulation in its category. This uncertainty ranges from approximately 2% (W+jets AH) to about 7% (tt SL). • QCD multijet background normalization. An uncertainty of 100% in the normalization is considered for QCD processes to cover effects in the kinematic tails that may not be well-modeled by the simulation. This has little overall impact on the final result, since the contribution from QCD multijet events is reduced to a negligible amount in this analysis.
• Single top quark background normalization. An uncertainty of 20% in the normalization is considered for single top quark processes, accounting for the uncertainty in the PDF and the effects from varying the factorization and renormalization scale parameters.
• Uncertainty related to ECAL mistiming. Partial mistiming of signals in the forward regions of the ECAL endcaps led to a minor reduction in trigger efficiency. To cover this effect, an additional uncertainty is applied on the signal acceptances of up to 10% in the forward jet categories. A potential effect on the background extrapolation into regions with forward jets is already taken into account by a dedicated systematic uncertainty.
The following sources of uncertainty affect the shape of the p miss T distribution, as well as the normalization of the various backgrounds and the signal, and are applied to all search channels: • PDF uncertainties. Uncertainties due to the choice of PDF are estimated by reweighting the samples with the NNPDF3.0 [47] replicas [58] and are applied to all backgrounds except for the single top quark, as these uncertainties are covered by the associated background normalization uncertainty.
• W/Z+heavy-flavor fraction. The uncertainty in the fraction of W/Z+heavy-flavor (HF) jets in W+jets and Z+jets event is taken into account. The relative contribution of W+HF and Z+HF are allowed to vary within 20% [59-62].
• Electroweak and QCD K factors. Uncertainties in the NLO/LO K factors calculated for W+jets and Z+jets processes are considered. These uncertainties account for missing higher-order corrections. For QCD, this comes from variations due to factorization and renormalization scales. For electroweak processes, an estimate of the size of the missing higher-order corrections is obtained by taking the difference between applying and not applying the NLO/LO electroweak K factors.
• Top quark p T reweighting. Differential measurements of the top quark p T spectrum in top quark pair production events [34] show that the measured p T spectrum is softer than in simulation. In order to improve the description of top quark pair events, simulated samples are reweighted to match the measurements. An associated systematic uncertainty is estimated by taking the difference between applying and not applying the reweighting.
• Factorization and renormalization scales. The uncertainties in the choice of the factorization and renormalization scale parameters are taken into account for the tt, tt+V, and diboson processes by applying a set of weights that represent a change of these scales by a factor of 2 or 0.5.
• Simulation sample size. Uncertainties due to the limited size of the simulated signal and background samples are included by allowing each bin of the distributions used in the signal extraction to fluctuate independently according to the statistical uncertainties in simulation, following Ref. [63].

Signal extraction
As previously discussed, the potential DM signal is expected to have the signature of tt or single top quark events with additional p miss T , therefore leading to an excess of events above the SM prediction in the p miss T spectrum. The DM signal is extracted from a simultaneous fit to the binned p miss T distribution in the various SRs and CRs, including all previously mentioned uncertainties. This global fit is performed as a binned maximum likelihood fit employing the ROOSTATS statistical package [64]. The main SM backgrounds were discussed previously in Section 4, and are dileptonic tt+jets and W+jets events for the SL SRs, and Z → νν, singlelepton tt+jets, and W+jets events for the AH SRs.
The effect of the systematic uncertainties in the shape and normalization of the p miss T spectrum, as discussed in the previous section, is taken into account by introducing nuisance parameters, which are constrained by the magnitude of the corresponding source of uncertainty. Uncertainties that affect normalization only are modeled using nuisances with log-normal probability densities. These parameters are treated as correlated between p miss T bins and between the different CRs and SRs within each channel. The sources common between SL and AH SRs and CRs are correlated across channels.
To improve the estimation of the main backgrounds, an unconstrained multiplicative parameter is assigned separately to each background for each bin of the p miss T spectrum. These multiplicative parameters scale the normalization of the associated background process simultaneously in the SRs and CRs for a given channel. For example, in a given p miss T bin of the SL selection, there is one multiplicative parameter for tt that links the tt background in the tt enhanced 2 CR, the W+jets enhanced 1 CR, and the SR. Therefore, the effect of contributions of the same background process in the different CRs is also taken into account. Additionally, potential contributions from the DM signals are included for all CRs and SRs, and scaled by a signal strength modifier µ = σ/σ th , i.e., the ratio between the measured and theoretical cross sections. Regions containing leptons (electrons and muons) are separated by lepton flavor.
The simultaneous fit to the binned p miss T distribution is performed combining the SL and AH regions. The values for the background multiplicative factors extracted from the fit are on average close to one, with a root-mean-square deviation that ranges from 5% to 21%, depending on the background processes and on the category considered (SL or AH). The post-fit distributions assuming the absence from the DM signal (i.e., the background only fit) are shown in Figs. 2 and 3 (4 and 5) for the SL and AH CRs (SRs), respectively. No significant excess at high p miss T in the SRs is observed. The SRs, both for the SL and AH channels, are divided into: 1 b-tagged jet and 0 forward jets, 1 b-tagged jet and ≥1 forward jets, and ≥2 b-tagged jets. The plots also contain the pre-fit distributions, represented by the dashed magenta line. The statistical and systematic uncertainties in the prediction are represented by hatched uncertainty bands, while the lower panels show the ratio of data and the post-fit prediction, and the bottom panels show the difference between the observed data events and the post-fit total background, divided by the full statistical and systematic uncertainty.

Results
Overall, data are found to be in agreement with the expected SM background in the SRs. Upper limits at 95% confidence level (CL) are computed on the ratio between the measured and theoretical cross sections µ, which is calculated with respect to the expected number of events for a scalar or pseudoscalar mediator and either the t/t+DM or tt+DM production modes separately, or summed together, where the results are referred to here as t, tt+DM. The theoretical cross sections for both signal models are obtained at LO. The limits are calculated using a modified frequentist approach with a test statistic based on the profile likelihood in the asymptotic approximation and the CL s criterion [65][66][67]. We test different mediator mass scenarios with m χ = 1 GeV and g q = g χ = 1 and the results are shown in Fig. 6 for scalar (left) and pseudoscalar (right) models. The expected limit for the t/t+DM signal alone is depicted by the blue dash-dotted line, while the expected tt+DM limit alone is given by the red dash-dotted line. The observed limit on the sum of both signals is represented by the black solid line, while its expected value is shown by the black dashed line with the 68 and 95% CL uncertainty bands in green and yellow, respectively.
For masses of the mediator particle below 200 GeV for the scalar model and below 300 GeV for the pseudoscalar model, the leading contribution to the sensitivity of the analysis stems from tt+DM. This behavior is mostly driven by the larger cross section for the tt+DM process when compared to the sum of the production processes for t/t+DM. However, the t/t+DM cross section drops less rapidly as a function of mediator particle mass in comparison to the tt+DM mode. Additionally, the p miss T spectrum for a given mediator mass leans towards higher values for the t/t+DM signal model when compared to the tt+DM model. These two features, together with the analysis specifically designed for both DM production modes and the statistical combination of the different SRs, lead up to a factor of two improvement at high mediator masses on the limits when compared to previous results [16]. In particular, the ≥1 forward jet category, which is specifically designed to enhance t/t+DM t channel events, improves the final results up to 14%.  Figure 6: The expected and observed 95% CL limits on the DM production cross sections, relative to the theory predictions, shown for the scalar (left) and pseudoscalar (right) models. The expected limit for the t/t+DM signal alone is depicted by the blue dash-dotted line, while the tt+DM limit alone is given by the red dash-dotted line. The observed limit on the sum of both signals is shown by the black solid line, while the expected value is shown by the black dashed line with the 68 and 95% CL uncertainty bands in green and yellow, respectively. The solid horizontal line corresponds to σ/σ th = 1. Table 3 represents the final combined limits (SL + AH) for the t/t+DM and tt+DM processes separately, and for the sum of the two processes.
Overall, we exclude mediator masses below 290 and 300 GeV for the scalar and pseudoscalar hypotheses, respectively. Table 3: Upper limits at 95% CL on the cross section ratio with respect to the expected DM signal for different scalar (φ) or pseudoscalar (a) mediator masses, m χ = 1 GeV, and g χ = g q = 1 for the combination of SL and AH signal regions. The median expected value and its 68 and 95% confidence intervals (CIs) are given.

Summary
The first search at the LHC for dark matter (DM) produced in association with a single top quark or a top quark pair in interactions mediated by a neutral scalar or pseudoscalar particle in proton-proton collisions at a center-of-mass energy of 13 TeV has been presented. The data correspond to an integrated luminosity of 35.9 fb −1 recorded by the CMS experiment in 2016. No significant deviations with respect to standard model predictions are observed and the results are interpreted in the context of a simplified model in which a scalar or pseudoscalar mediator particle couples to the top quark and subsequently decays into two DM particles.
Scalar and pseudoscalar mediator masses below 290 and 300 GeV are excluded at 95% confidence level assuming a DM particle mass of 1 GeV and mediator couplings to fermions and DM particles equal to unity. This analysis provides the most stringent limits derived at the LHC for these new spin-0 mediator particles.   [13] ATLAS Collaboration, "Search for dark matter in events with heavy quarks and missing transverse momentum in pp collisions with the ATLAS detector", Eur. Phys. J. C 75 (2015) 92, doi:10.1140/epjc/s10052-015-3306-z, arXiv:1410.4031.
[50] CMS Collaboration, "Investigations of the impact of the parton shower tuning in PYTHIA 8 in the modelling of tt at √ s = 8 and 13 TeV", CMS Physics Analysis Summary CMS-PAS-TOP-16-021, 2016.