Search for charginos and neutralinos in final states with two boosted hadronically decaying bosons and missing transverse momentum in pp collisions at √s = 13 TeV with the ATLAS detector

A search for charginos and neutralinos at the Large Hadron Collider using fully hadronic final states and missing transverse momentum is reported. Pair-produced charginos or neutralinos are explored, each decaying into a high-p T Standard Model weak boson. Fully hadronic final states are studied to exploit the advantage of the large branching ratio, and the efficient rejection of backgrounds by identifying the high-p T bosons using large-radius jets and jet substructure information. An integrated luminosity of 139 fb − 1 of proton-proton collision data collected by the ATLAS detector at a center-of-mass energy of 13 TeV is used. No significant excess is found beyond the Standard Model expectation. Exclusion limits at the 95% confidence level are set on wino or higgsino production with various assumptions about the decay branching ratios and the type of lightest supersymmetric particle. A wino (higgsino) mass up to 1060 (900) GeVis excluded when the lightest supersymmetry particle mass is below 400 (240) GeVand the mass splitting is larger than 400 (450) GeV. The sensitivity to high-mass winos and higgsinos is significantly extended relative to previous LHC searches using other final states.


Introduction
Supersymmetry (SUSY) [1][2][3][4][5][6] is a theoretical framework that extends the Standard Model (SM) by introducing new particles ("superpartners") that have the same quantum numbers as the SM particles except for their spins. In the Minimal Supersymmetric Standard Model (MSSM) [7,8], the bino ( ), wino ( ) and higgsino ( ) are the superpartners of the U(1) Y and SU(2) gauge fields and the Higgs field, respectively. These are collectively referred to as "electroweakinos" and form the chargino ( ± 1 , ± 2 ) and neutralino ( 0 1 , 0 2 , 0 3 , 0 4 ) mass eigenstates through mixing, with the subscripts indicating increasing mass. The lightest neutralino ( 0 1 ) is often considered to be the lightest SUSY particle (LSP), since it is then a viable candidate for weakly interacting massive particle (WIMP) dark matter [9,10] when -parity conservation is assumed [11]. Depending on the signal model, the superpartner of the graviton (gravitino, ) and the axion (axino,˜) are alternatively considered to be the LSP and dark matter candidate in this search.
Electroweakinos with masses of the order of 0.1-1 TeV are motivated by various phenomenological arguments: (1) the mass of the neutralino LSP dark matter candidate is constrained to be less than a few TeV by the observed relic density [12,13]; (2) the higgsino mass is also motivated to be of the same order as the boson mass by naturalness arguments [14][15][16][17]; (3) the MSSM parameter space explaining the discrepancy between the measured muon anomalous magnetic moment [18] and its SM predictions [19] tends to include electroweakinos with masses from 200 GeV to 1 TeV [20][21][22].
This search targets the pair production of electroweakinos ( heavy ), where each of them decays into a lighter one ( light ) and an on-shell , or SM Higgs boson (ℎ). A mass splitting Δ ( heavy , light ) greater than 400 GeV is considered in the search. The heavy can be either wino-or higgsino-like, and light can be a bino-, wino-, higgsino-like chargino/neutralino, gravitino or axino as discussed more in Section 2. Cases where both heavy and light are wino-or higgsino-like are not considered as they lead to very small Δ ( heavy , light ). The light is either the LSP or an electroweakino nearly degenerate with it, leading to missing transverse momentum (p miss T , with magnitude miss T ) in the decay signature. 1 The analysis focuses on the hadronic decay modes of the , and ℎ bosons, namely →¯, →¯/ā nd ℎ →¯, where (¯) represents light-flavor (anti)quarks , , , (¯,¯,¯,¯). 2 Two fully hadronic final states are considered: the final state which involves two / bosons decaying into two light-flavor (anti)quarks, and the final state where a / boson decays into two light-flavor (anti)quarks and a /ℎ boson decays into¯, as illustrated by the diagrams in Figure 1. While the ATLAS and CMS experiments have typically searched for pair production of electroweakinos using leptonic decay modes [23][24][25][26][27][28][29][30][31][32][33][34][35][36], the fully hadronic final states are used here to take advantage of the larger , , and ℎ hadronic branching ratios, and thereby provide sensitivity to the production of heavier electroweakinos, despite their smaller production cross-sections. Final states with four -quarks are not considered in this search, 3 however they are targeted by dedicated searches [37,38] complementing the sensitivity.
The multĳet background is highly suppressed by requiring large miss T , and the dominant backgrounds in the analysis are (→ ) + jets, (→ ℓ ) + jets, and diboson production. These can be effectively suppressed when targeting large mass splittings between the produced electroweakinos and the LSP by selecting high-T kinematics and explicitly reconstructing the two boosted SM bosons. The "boson tagging" technique developed by the ATLAS and CMS experiments [39][40][41][42] is employed, where a single large-radius jet is used to capture the two collimated energetic jets from each boosted SM boson decay. Analysis of the jet substructure helps to identify the hadronic decays of , , and ℎ. are interpreted in terms of various models with different electroweakino types for ( heavy , light ) and their branching ratio assumptions. While each model predicts different electroweakino branching ratios into , or ℎ, the search combines the dedicated event selections for each decay to achieve only a small model dependency.

Target physics scenarios and the signal models
Three physics scenarios are considered in the analysis: • a baseline MSSM scenario where bino, wino and higgsino are considered as heavy or light ; • a scenario with a gravitino LSP and light higgsinos inspired by the general gauge mediation (GGM) [43][44][45][46][47] models and naturalness; • a scenario with an axino LSP assuming the SM extension with a QCD axion and light higgsinos driven by naturalness.
The signal models considered in the analysis, derived from each physics scenario, are described in the following subsections. Each model is labeled as ( , ), where and represent the dominant component of heavy and light respectively. All the SUSY particles other than heavy and light are assumed to be decoupled in mass [48,49]. The production modes, final states, and the branching ratio assumptions for heavy are summarized in Table 1.
The wino-like states form a doublet consisting of one chargino and one neutralino, while the higgsino-like triplet includes an additional neutralino. The mass degeneracy within these wino/higgsino multiplets is dictated by the extent of mixing with the other electroweakino states, which is characterized by Δ ( heavy , light ). With Δ ( heavy , light ) > 400 GeV, the mass splittings within the wino/higgsino multiplets are typically small (< 10 GeV), such that the decay products are almost never reconstructed in this analysis. Therefore, these multiplets are treated as approximately degenerate; the masses of the wino-like chargino and wino-like neutralino are treated as identical; and the masses of the higgsino-like chargino and the heavier higgsino-like neutralino are set 1 GeV heavier than the lighter higgsino-like neutralino. 4 Particles originating from decays within the wino doublet or higgsino triplet are ignored so that only the decays of the form heavy → light + / /ℎ are taken into account in the signature.

Bino/wino/higgsino LSP models:
( , ), ( , ), ( , ), ( , ) Two of , , and are assumed to be light enough to be produced at the LHC while the others are decoupled. Four mass hierarchies are experimentally explorable under this regime: ( , ), ( , ), ( , ), ( , ). A bino is not considered as heavy here because the production cross-section is negligible when sfermions and the non-SM Higgs bosons are decoupled.
The ( , ), ( , ), and ( , ) hierarchies are typically predicted by the MSSM parameter space that explains the muon − 2 anomaly with loop contributions including a wino and/or higgsino. 5 The ( , ) hierarchy is additionally motivated when the mass of a bino-like LSP is half of the or ℎ boson mass, where the LSP dark matter can annihilate via the /ℎ resonance (" /ℎ-funnel" dark matter). This is a special case of well-tempered neutralino dark matter where a large higgsino-bino mass splitting is favored in order to realize the observed relic density [50][51][52][53][54][55][56][57][58].
A signal model is defined for each of the four hierarchies with a set of assumptions described below. The mass spectra, corresponding mass eigenstates, and the decays considered are illustrated for each model in Figure 2. A boson is generated when a chargino decays into a neutralino or vice versa; a or ℎ boson is emitted when a chargino decays into a chargino, or a neutralino decays into a neutralino.
The production modes considered in each model are shown in Table 1. Processes involving charginochargino or chargino-neutralino pair production are taken into account in the wino production models: ( , ) and ( , ); and neutralino-neutralino production is additionally included for the higgsino production models: ( , ) and ( , ). It should be noted that neutral wino pairs can only be produced via a -channel exchange of squarks, which is prohibited when squarks are assumed to be decoupled.
In ( , ) and ( , ) models, however, the branching ratios of the produced chargino and neutralino(s) are largely dictated by the three MSSM parameters: the wino mass parameter 2 , the higgsino mass parameter , and the ratio of vacuum expectation values of the two Higgs fields, tan . For a given set of ( 2 , , tan ) the branching ratios are coherently derived using SOFTSUSY 4.1.7 [60,61], with all the SUSY mass parameters except for 2 and being set as decoupled. The signal models are tested using various 5 In this case, smuons and muon sneutralinos need to be only mildly decoupled in mass. 6 Produced heavier higgsinos are assumed to decay 100% directly into bino rather than the lighter higgsino, which is reasonable given the much smaller higgsino mass splitting compared with the higgsino-bino mass splitting. 7 Particularly for ( , ), the branching ratio strongly depends on the higgsino mass as the decay has to rely on the small higgsino component, even though it is assumed to be decoupled. production, ( ± 1 / 0 2 ) < 640 GeV is excluded for ( 0 1 ) < 300 GeV when the 0 2 is assumed to decay into and 0 1 with 100% probability [24-26, 32, 33]. Alternatively, ( ± 1 / 0 2 ) < 740 GeV is excluded for ( 0 1 ) < 250 GeV [27-29, [34][35][36] when the 0 2 decays solely into ℎ and 0 1 . Figure 2: The electroweakino mass spectra and corresponding mass eigenstates in each model in the bino/wino/higgsino LSP scenario. The solid (dashed) arrows represent the decay modes emitting a ( or ℎ) boson. A boson is generated when a chargino decays into a neutralino or vice versa; a or ℎ boson is emitted when a chargino decays into a chargino or a neutralino decays into a neutralino.

GGM/naturalness-driven gravitino LSP model: ( , )
General gauge mediation (GGM), a class of SUSY breaking scenarios characterized by a messenger sector to which only SM gauge bosons can couple, typically predicts a nearly massless gravitino ( ) as the LSP. Motivated also by the naturalness argument, the production of a relatively light higgsino triplet ( ± 1 , 0 2 , 0 1 ) decaying into a gravitino LSP has been explored at ATLAS [30, 37] and CMS [33], as illustrated in Figure 3(a). All of the four production modes are considered together: ± 1 ∓ 1 , ± 1 0 1 , ± 1 0 2 , 0 1 0 2 . A moderately small higgsino-gravitino coupling is considered in this analysis, where the produced heavy higgsinos ( ± 1 / 0 2 ) always decay into a gravitino via the lightest neutral higgsino ( 0 1 ), while the 0 1 still has a short enough lifetime to be regarded as decaying promptly. In this model, the 0 1 decays into a gravitino and either a or ℎ boson, where the branching ratio is treated as a free parameter and scanned in the limit setting. The previous searches have excluded masses of 0 1 lighter than 650-880 GeV depending on the branching ratio [30, 33, 37].

Naturalness-driven axino LSP model: ( ,˜)
While the QCD Lagrangian generally allows for CP violation, the absence of such observation suggests a highly unnatural tuning of the parameters in the theory, referred to as the "strong CP problem". The Peccei-Quinn mechanism aims to solve this problem by introducing an additional chiral U(1) symmetry [62]. Through its spontaneous symmetry breaking, the CP-violating term vanishes dynamically, leaving a Nambu-Goldstone boson known as the axion [63,64]. In the SUSY extension, the axino is introduced as the superpartner of the axion. A model including a light higgsino triplet ( ± 1 , 0 2 , 0 1 ) decaying into an axino LSP [65] is proposed in the spirit of pursuing naturalness as well as axion/axino dark matter [66,67] under -parity conservation.
A diagram of the model is shown in Figure 3(b). Higgsinos are produced by each of the four modes: ± 1 ∓ 1 , The produced heavier higgsinos ( ± 1 / 0 2 ) are assumed to always decay into the axino via the lightest natural higgsino ( 0 1 ). This is typically valid when the wino and bino are reasonably decoupled so as to maintain approximate mass degeneracy of the higgsino triplet, and when the conventionally motivated range of the axion coupling constant is assumed [65]. A prompt 0 1 decay into an axino and a or ℎ boson is considered in the search. The value of the branching ratio B ( 0 1 →˜) (= 1 − B ( 0 1 → ℎ˜)) is scanned over 25%, 50%, 75%, and 100% in the interpretation. The model is similar to ( , ), except that the LSP can be massive.

ATLAS detector
The ATLAS experiment [69,70] is a multipurpose detector with a forward-backward symmetric cylindrical geometry and nearly 4 coverage in solid angle. 8 It consists of an inner tracking detector surrounded by a thin superconducting solenoid providing a 2 T axial magnetic field, electromagnetic and hadron calorimeters, and a muon spectrometer.
The inner detector (ID) consists of pixel and microstrip silicon detectors covering the pseudorapidity range | | < 2.5 and a transition radiation tracker covering | | < 2.0. Outside the ID, a lead/liquid-argon (LAr) electromagnetic calorimeter (ECAL) and a steel/scintillator-tile hadronic calorimeter cover the | | < 3.2 and | | < 1.7 ranges, respectively. In the forward regions, a copper/LAr endcap calorimeter extends the coverage of hadronic measurements to 1.7 < | | < 3.2, while copper/LAr and tungsten/LAr forward calorimeters are employed for electromagnetic and hadronic measurements in the 3.1 < | | < 4.9 region. The muon spectrometer (MS) surrounds the calorimeters and comprises three layers of trigger and high-precision tracking chambers spanning | | < 2.4 and | | < 2.7, respectively. A magnetic field is 8 ATLAS uses a right-handed coordinate system with its origin at the nominal interaction point (IP) in the center of the detector and the -axis along the beam pipe. The -axis points from the IP to the center of the LHC ring, and the -axis points upwards. Cylindrical coordinates ( , ) are used in the transverse plane, being the azimuthal angle around the -axis. The pseudorapidity is defined in terms of the polar angle as = − ln tan( /2). Angular distance is measured in units of Δ ≡ √︁ (Δ ) 2 + (Δ ) 2 . Table 1: Summary of the production modes, final states, and signal regions (SRs) used for the hypothesis tests, and the branching ratio assumptions for the signal models targeted in the search. The notation and definition of the SRs are described in Section 6.2. The ( , ) and ( , ) models are used to optimize the selection, and the rest are considered in the interpretation. The ( , ) simplified models (( , )-SIM) discussed in Section 4.2.2 are also interpreted in order to allow comparisons with the ATLAS electroweakino search results [23-25, 29, 68].

Model Production
Final states SRs simultaneously fitted Branching ratio provided by a system of three superconducting air-core toroidal magnets with eight coils each. The field integral of the toroids ranges between 2.0 and 6.0 T m across most of the detector.
Events of interest are selected and collected by the ATLAS trigger system [71], consisting of a hardwarebased first-level trigger (L1) and a software-based high-level trigger (HLT). The L1 trigger is designed to accept events from the 40 MHz bunch crossings at a rate below 100 kHz, and the HLT reduces this to about 1 kHz, the rate at which events are recorded to disk. An extensive software suite [72] is used for real and simulated data reconstruction and analysis, for operation and in the trigger and data acquisition systems of the experiment.

Data sample
The data events used in the analysis are from proton-proton collisions at √ = 13 TeV, recorded during stable beam conditions at the LHC during 2015-2018. The collected dataset corresponds to an integrated luminosity of 139 fb −1 after applying the data quality criteria [73]. The primary dataset was collected by triggers targeting large missing transverse momentum [74]. Events were accepted when the miss T calculated at trigger level was greater than 70-110 GeV, with the threshold rising with the increased instantaneous luminosity during the data-taking period. The efficiency reaches approximately 100% for events reconstructed offline with miss T > 200 GeV, which is generally required in the search. Auxiliary data samples used to validate the background estimation were selected using triggers requiring at least one isolated electron, muon or photon [75,76]. The thresholds were T = 24 (26) GeV for electrons, T = 20 (26) GeV for muons, and T = 120 (140) GeV for photons in data taken in 2015 (2016-2018).

Monte Carlo simulation
Monte Carlo (MC) simulations are used to estimate the SM backgrounds and the signals contributing to the analysis regions. All of the generated events were propagated through the ATLAS detector simulation [77] based on G 4 [78]. Multiple proton-proton collisions in the same and neighboring bunch crossings (pileup) were modeled by overlaying the hard-scatter events with minimum-bias events simulated by P 8.186 [79] with a set of tuned parameters called the A3 tune [80] and NNPDF2.3 parton distribution function (PDF) set [81].
The simulated events are processed with the same trigger and reconstruction algorithms as the data. The lepton and photon trigger efficiencies in the simulation are corrected to match those in data using scale factors that depend on the T and of the leptons and photons, as derived from control samples [75,76].

Standard Model backgrounds
Events with a leptonically decaying or boson associated with hadron jets were simulated using the S 2.2.1 [82] generator. The matrix elements were calculated for up to two partons at next-toleading-order (NLO) accuracy and up to four jets at leading-order (LO) accuracy using the Comix [83] and O L [84,85] generators. The NLO matrix elements for a given jet multiplicity were matched to the parton shower using a color-exact variant of the MC@NLO algorithm [86]. Different jet multiplicities were then merged into an inclusive sample using an improved Catani-Krauss-Kuhn-Webber (CKKW) matching procedure [87,88] which is extended to NLO accuracy using the MEPS@NLO prescription [89]. The NNPDF3.0 PDF sets [90] were used. Prompt single-photon production, denoted by +jets, was simulated using the same configuration except that the generator version is S 2.2.2. The photons must be isolated according to a smooth-cone isolation criterion [91].
Samples of¯and single-top-quark ( + ) events were generated with P B [92][93][94][95] v2 at NLO with the NNPDF3.0 PDF sets. The top quark mass was set to 172.5 GeV. For the¯generation, the ℎ damp parameter, which controls the T of the first additional emission beyond the Born configuration in P , was set to 1.5 times the top quark mass. This is to regulate the high-T emission recoiling against the¯system so as to reproduce the data [96]. The parton shower, fragmentation, and underlying event were simulated using P 8.230 with the NNPDF2.3 PDF set and the A14 tune [97]. The decays of bottom and charm hadrons were performed by E G 1.6.0 [98]. The diagram removal scheme [99] was employed to account for the interference between¯and single-top production.
Events containing¯with additional heavy particles, such as¯+ / /ℎ,¯+ ,¯and¯¯, are collectively referred to as¯+ in this paper. , where all the processes at orders of ( 4 , 2 s ) and ( 6 , 0 s ) are taken into account, including the off-shell contributions and those mediated by Higgs bosons. The NNPDF3.0 PDF set is used.
The production of ℎ and ℎ (collectively denoted by ℎ) was modeled by P B 2.2 interfaced with P 8.186. The ℎ mass was set to 125 GeV. The NNPDF3 PDF set and the AZNLO tune were used.
Theoretical cross-sections are used to normalize the generated background samples. The¯sample is normalized to the cross-section predicted at NNLO in QCD, including the resummation of next-to-next-toleading-logarithmic (NNLL) soft-gluon terms calculated using T ++ 2.0 [101][102][103][104][105][106][107]. The cross-sections of single-top-quark -and -channel production are calculated using the H 2.1 program [108,109], while the -channel calculation followed the prescriptions from Refs. [110,111]. The and¯+ / samples are normalized to the cross-sections calculated at NLO [112][113][114]. The cross-sections for the +jets and +jets samples are calculated at NNLO [115]. The remaining samples, including the production of ℎ, , ,¯+ , ,¯+ℎ,¯, and¯¯, are normalized to the cross-sections calculated at LO by the generators.

Signals
The ( , ), ( , ), ( , ), ( , ), and ( ,˜) models discussed in Section 2 were simulated by combining the "simplified model" signals in which a fixed production mode and decay chain are considered. Six sets of simplified model samples, derived as variants of the ( , ) and ( , ) models, were generated to cover the different final states: • wino pair production, with each decaying into a bino LSP (referred to as ( , ) simplified models, or ( , )-SIM):

N2N3-ZZ:
The symbols C1, N2 and N3 in the model names represent ± 1 , 0 2 and 0 3 respectively. These samples are used to model event kinematics for each final state; the ( , )-SIM samples are used for events containing , or ℎ, while ( , )-SIM samples are employed to model events with , ℎ or ℎℎ. The events are reweighted to account for the production cross-section and the heavy branching ratios considered in the model. The underlying assumption for the method is that the event kinematics depend only on ( light ) and ( heavy ), and not on the production mode, type of LSP, or other MSSM variables such as tan . This is validated in the phase space considered in the analysis by using generator-level samples.
The ( , )-SIM models are also considered in the limit interpretation in order to allow comparisons with previous ATLAS searches, as these models were the most commonly studied [23, [25][26][27][28][29][31][32][33][34][35][36]. For the ( , ) model, a dedicated sample was generated with all of the production modes ( ± Decays of the produced electroweakinos are simulated using P . As the performance of the boosted / boson tagging is known to be sensitive to the boson polarization [39], / bosons from the electroweakino decays were carefully modeled either by using M S 2.7.3 [117,118] or by reweighting the helicity angle distribution such that the overall cross-section remains unchanged. The / bosons are typically longitudinally polarized when Δ ( heavy , light ) > 400 GeV is considered. The decays of bottom and charm hadrons were performed by E G 1.2.0.

Event reconstruction
The primary reconstructed objects used in the analysis are large-radius (large-) jets, denoted by . These are reconstructed from locally calibrated topo-clusters [126] using the anti-algorithm [127] implemented in the F J package [128] with a radius parameter = 1.0. A trimming algorithm [129]  , is calculated according to the combined mass prescription [131] in order to achieve the best mass resolution, which is given by the weighted sum of masses computed using only the calorimeter information and with tracking information included. The T and mass scales are calibrated using simulation, followed by an in situ calibration [132] to correct for residual differences between data and MC simulation. Large-jets used in the analysis are selected with T > 200 GeV, | | < 2.0, and > 40 GeV.
Track-jets are used to identify large--jet subjets that contain -hadrons. They are reconstructed from ID tracks by using the anti-algorithm with a sliding radius parameter = 30 GeV/ T truncated at 0.02 and 0.4 [133]. The MV2c10 -tagging algorithm [134] is applied to track-jets satisfying T > 20 GeV and | | < 2.5. This algorithm is a multivariate discriminator that utilizes track impact parameters, the presence of secondary vertices, and the trajectories of -and -hadrons inside the jet. A working point is chosen such that -jets from simulated¯events are identified with 85% efficiency, with rejection factors of 3 against -jets and 33 against jets originating from other light-flavor quarks or gluons [134]. Efficiency correction factors are applied to simulated samples to account for the efficiency difference observed between MC and data events in dedicated measurement regions in which¯samples [134,135] are used for -jets and mistagged -jets, and +jets samples are used for the light-flavor mistagged jets [136].
The -jet multiplicity of each large-jet is defined by the number of -tagged track-jets it contains. A -tagged track-jet is contained in a large-jet if the two jet axes have an angular separation Δ < 1.0.
Two types of boosted boson tagging are employed for the preselected large-jets to identify the SM boson decays: ( )-tagging targeting ( ) → , and (ℎ )-tagging targeting (ℎ) → . The ( )-tagging utilizes cuts on , the energy correlation function 2 , and the track multiplicity trk [39,137]. In order to maintain orthogonality with the (ℎ )-tagging, the -jet multiplicity of the large-jet is required to be less than two. The 2 variable is defined as a ratio of three-point to two-point energy correlation functions [138, 139] based on the energies and pairwise angular separations of particles within a jet, and trk is the number of tracks matched to the large-jet by ghost association [140] before trimming is applied. While the upper bound on trk is fixed, the cut values applied to and 2 are shifted smoothly as a function of T to maximize the rejection for typical single-parton-initiated jets, maintaining a constant efficiency for signal jets that contain the decay products of ( ) → . The selections are optimized separately for targeting and , and -tagging is defined by a selection satisfying either the -or -tagging. The cut values applied to and 2 are the same as for the "50% efficiency / tagger" in Ref. [137]; however, a loosened trk cut, from ≤ 26 to ≤ 32 (34) for the ( )-tagging, is applied to achieve the optimum sensitivity for the analysis and better modeling in MC simulation. The efficiency correction factors are rederived according to this refinement, using a methodology similar to that described in Ref. [137]. A correction factor of 0.85-1.05 is typically obtained. The performance of the -and -tagging is summarized in Figures 4(a) and 4(b). For a sample of preselected largejets ( T > 200 GeV, | | < 2.0, > 40 GeV), the tagging efficiency is about 50% for the signal jets originating from electroweakino decays, 9 while the background rejection is typically about 10 (40) at T = 200 (1000) GeV per jet in (→ ) + jets events.

The
(ℎ )-tagging is applied to large-jets containing exactly two -jets (denoted by ) by applying a jet mass window cut that selects the peak consistent with (ℎ) bosons. The jet mass is corrected by adding the momentum of the highest-T muon identified inside the large-jet in order to improve the resolution of the mass peak. The > 40 GeV) in the simulated ( , )-SIM signal events with Δ ( heavy , light ) ≥ 400 GeV. The jets are matched with generator-level / ( /ℎ) bosons within Δ < 1.0 which decay into¯(¯). The background jet rejection factor is calculated using preselected large-jets in the sample of simulated (→ ) + jets events, dominated by initial-state radiation jets. As in the /ℎ -tagging, the rejection factor is shown as a function of the number of -or -quarks contained in the large-jet within Δ < 1.0. The efficiency correction factors are applied to the signal efficiency and background rejection for the / -tagging. The uncertainty is represented by the hashed bands, which includes the MC statistical uncertainty and the systematic uncertainties discussed in Section 8.1.
The analysis also uses reconstructed electrons and muons (collectively referred to as "leptons"), as well as photons and small-radius (small-) jets, for kinematic selection, validation of the background estimation, and the miss T computation. Electron candidates are reconstructed from energy clusters that are consistent with electromagnetic showers in the ECAL and are matched to tracks in the ID, which are calibrated in situ using → samples [141]. Muon candidates in the detector are typically reconstructed by matching tracks in the MS to tracks in the ID, and they are calibrated in situ using → and / → samples [142]. Small-jet candidates are reconstructed from particle-flow objects [143] calibrated at the electromagnetic scale using the anti-algorithm with a radius parameter of = 0.4. After subtracting the expected energy contribution from pileup using the jet area technique [144], the jet energy scale and resolution are corrected to particle level using MC simulation as well as by in situ calibration using +jets, +jets, and multĳet events [145]. Photon candidates are reconstructed either as electromagnetic clusters with no matching ID track or as + − pairs from photon conversions in the ID material [141].
Reconstructed electrons, muons, small-jets, and photons are subject to two sets of identification criteria: the looser "baseline" criteria and the tighter "signal" criteria. The baseline objects are used for the miss T computation, event cleaning, and the overlap removal procedure that resolves ambiguities between reconstructed objects as described below. Baseline electrons are required to have T > 4.5 GeV and | | < 2.47, and meet the Loose criteria of the likelihood-based identification [141]. Baseline muons are required to have T > 3 GeV and | | < 2.7, and to meet the Medium identification criteria defined in Ref. [142]. To suppress the contributions from pileup, baseline leptons are also required to have a trajectory consistent with the primary vertex, i.e. | 0 sin | < 0.5 mm. 10 Baseline small-jets must have T > 30 GeV and | | < 4.5. Baseline photons must meet the Tight identification criteria [141] in addition to satisfying T > 50 GeV and | | < 2.4.
To prevent the reconstruction of a single particle as multiple objects, an overlap removal procedure is applied to the baseline leptons, photons, and jets in the following order. Any electron sharing an ID track with a muon or other electrons is removed, such that only the highest-T electron is kept when multiple electrons share the same ID track. Photons around the remaining electrons and muons are removed if the photon-lepton separation is Δ < 0.4. Next, small-jets are removed if they are Δ < 0.2 from a remaining electron or photon, or are Δ < 0.4 from a muon and the jet has fewer than three associated tracks with T > 500 MeV. Leptons (photons) are removed if they are separated from a remaining smalljet by Δ < min(0.4, 0.04 + 10 GeV/ T,ℓ ) (Δ < 0.4). Finally, large-jets are removed if they are separated by Δ < 1.0 from any remaining electrons.
The missing transverse momentum, with magnitude miss T , is calculated as the negative vectorial sum of the transverse momenta of all baseline leptons, photons and small-jets calibrated to their respective energy scales, and an additional "soft term" constructed from tracks originating from the primary vertex but not associated with any of the baseline objects [146].
Signal objects are defined by applying additional stringent criteria in event selections to ensure a high selection purity. Signal leptons and photons are used only for validating the background estimation. Signal electrons and muons must satisfy T > 10 GeV. Additionally, the Tight identification criteria are imposed for signal electrons. To reduce the contribution from nonprompt decays of heavy-flavor hadrons, the significance of the transverse impact parameter is required to satisfy | 0 / ( 0 )| < 5 (3) for signal electrons (muons). An isolation selection is imposed to further suppress the residual misidentified leptons originating from jets. The Tight (HighPtCaloOnly) working point defined in Ref. [141] is used for signal electrons with T < 200 GeV ( T > 200 GeV), and the Tight working point defined in Ref.
[142] is applied for signal muons. Signal small-jets are selected within | | < 2.8, and those with T < 120 GeV and | | < 2.5 must satisfy the Tight quality criteria of the track-based jet vertex tagger [144,147] in order to suppress jets originating from pileup. Signal photons are defined as baseline photons with T > 200 GeV. 6 Event selection 6

.1 Common preselection
As described below, a common preselection is used when defining the signal regions (SRs), as well as the control regions (CRs) and validation regions (VRs) used for background estimation.
After the trigger requirement discussed in Section 4.1, both the data and MC events are required to have at least one reconstructed vertex that is associated with two or more tracks with T > 500 MeV. The primary vertex of each event is selected as the vertex with the largest 2 T of associated tracks [148]. A set of event cleaning criteria are applied to ensure the quality of the measurements as well as to veto noncollision backgrounds. Events with baseline muons consistent with cosmic muons (| 0 | > 1 mm or | 0 | > 0.2 mm) are removed. To avoid pathological miss T reconstruction, events are also vetoed if any baseline jet points to a module of the tile hadronic calorimeter that was not operational during the data taking, or if any baseline muon suffers from a poor momentum measurement ( ( / )/( / ) > 0.4 where / is the measured charge divided by momentum and ( / ) is its uncertainty). Track-jets can overlap due to their T -dependent variable radius. To avoid the ambiguous cases of concentric jets, events with a track-jet overlapping with another track-jet are removed. An overlap is defined by Δ < min where Δ is the angular distance between a given pair of track-jets and min is the smaller of their radii.
Beam-induced background is one of the major noncollision backgrounds. Particles (typically muons) generated by interactions between the beam and the upstream collimators may directly hit the detector material. Particular care is needed for this background since it can cause high-energy jet-like signatures and, correspondingly, large spurious miss T [149]. Events with baseline jets failing the Loose cleaning [150] are removed. Further cleaning is applied to events with no baseline leptons and photons by requiring consistency between the miss T and the alternative miss T,track , which is computed using only the good-quality ID tracks associated with the primary vertex. While those events satisfy miss T > 200 GeV due to the trigger requirements, miss T,track > 75 GeV and Δ ( miss T , miss T,track ) < 2.0 are required in the cleaning procedure. Finally, events are required to contain at least two large-jets. Events with three or more large-jet are not vetoed in the analysis, in order to be as inclusive as possible in model coverage; however, the applied event selection is always based on the two highest-T large-jets.

Signal region selection
After the preselection, events with no baseline leptons are selected. Two orthogonal signal region categories, 4Q and 2B2Q, based on the absence or presence of a large-jet containing exactly two -tagged track-jets ( ) are defined in order to target the and final states respectively.
Boson tagging is required for the two leading large-jets. The two leading large-jets must pass the -tagging ( ( ) = 2) in SR-4Q. On the other hand, SR-2B2Q requires to satisfy the -or ℎ -tagging, while the other jet (denoted by ) must satisfy the -tagging criteria. Multiple SRs are defined in each SR category to target the different final states of the signal processes. The WW (ZZ) region in SR-4Q requires both leading large-jets to pass the ( )-tagging, while the WZ region is defined to contain events with at least one -tagged jet and at least one -tagged jet. The inclusive bin VV, the logical union of WW, WZ and ZZ, is designed to cover the signals that can have various heavy branching ratios into relative to . Similarly, the 2B2Q category accommodates the WZ/ZZ/Wh/Zh region, which varies by -versus -tagging or -versus ℎ -tagging. The VZ (Vh) region serves as the inclusive SR defined by the logical union of WZ and ZZ (Wh and Zh). The SR segmentation is summarized in Table 2 and illustrated schematically in Figure 5. In total, four and six SRs are defined in SR-4Q and SR-2B2Q, respectively. SR-4QXX, SR-2B2QXZ and SR-2B2QXh (where X = W/Z/V) are mutually orthogonal, and are statistically combined according to the models considered in Table 1. Table 2: Definition of each SR in the 4Q and 2B2Q categories, where ( ), ( ), ( ), ( ), and (ℎ ) are respectively the number of large-jets passing the -, -, -, -, and ℎ -tagging of the two leading large-jets. SR-4Q-WZ requires ( ), ( ) ≥ 1 instead of ( ) = ( ) = 1 because the selections in -and -tagging are not exclusive. The overlap and the segmentation between the SRs are illustrated in Figure 5.
Finally, a series of background rejection cuts are applied to further suppress the main SM backgrounds at this stage: (→ ) + jets, (→ ℓ ) + jets, and¯. To reject backgrounds including top quarks, the -jets that do not originate from boosted boson candidates are vetoed. This is done by requiring unmatched where unmatched -jet denotes the number of -tagged track-jets that are not matched with any of the two leading large-jets by the Δ < 1.0 criterion. For the 4Q category, the total number of -tagged track-jets in the event ( -jet ) must be less than two in order to further suppress the¯background. The effective mass variable, eff , defined as the scalar sum of the T of the two leading large-jets and miss T , is used to select events with hard kinematics together with miss T . Selections of miss T > 300 (200) GeV and eff > 1300 (1000) GeV are applied in the 4Q (2B2Q) regions. Event shape information is also useful in distinguishing the signals from the backgrounds. Signal events have a relatively spherical shape, where the jets tend to be isolated from miss T , indicating that the heavy electroweakinos are produced nearly at rest. Background events, meanwhile, have a higher chance of containing a jet aligned with a / boson, due to boosted top decays in¯events or the emission of collinear radiation in +jets and +jets events. The minimum azimuthal angle separation between miss T and any signal small-jets must satisfy min Δ ( miss T , ) > 1.0 for both SR categories. The use of small-jets is motivated by the need to effectively identify low-T jets and it also provides better resolution in terms of jet alignment with miss T . In 2B2Q, the stransverse mass variable T2 [151, 152] is also used, constructed by assigning each of the two leading large-jets to the visible particle legs. A selection of T2 > 250 GeV is found to effectively suppress the SM backgrounds, particularly¯which exhibits a kinematic cutoff at T2 ∼ 200 GeV, driven by the top quark mass constraint. 11 11 The hypothetical missing-particle mass is set to 100 GeV and this offset is subtracted from the calculated T2 , although the ≈ 800 GeV. The discovery significance is used as the metric of sensitivity. The obtained cuts are also found to be nearly optimal for the other signal models.
The estimation strategy varies between the "reducible" and "irreducible" backgrounds. The irreducible backgrounds in this search are due to SM events including at least two hadronic / /ℎ decays and large miss T from high-T neutrinos. These consist of and¯+ , and are estimated using MC simulation. The contribution from fully-hadronic and¯are negligible due to the stringent miss T requirement.
The reducible backgrounds are all that remain, including the dominant (→ ) + jets production. These backgrounds are characterized by the presence of at least one large-jet that originates from a process dependency on the choice of missing-particle mass is very small.
other than a / → decay and is referred to as a "fake boson jet". The fake boson jets are typically caused by two collimated high-T initial-state radiation (ISR) jets that are clustered together as a single large-jet. A partly data-driven method is used to estimate the reducible backgrounds; a control region (CR0L) 12 is defined in the phase space adjacent to a SR, where the MC sample is normalized to the CR data. The SR expectation is obtained using the normalized MC sample, assuming the modeling of the SR/CR0L yield ratio ("0L transfer factor") is reliable. This assumption is tested in a number of validation regions (VR) in data. A CR0L is defined for each 4Q and 2B2Q category by reversing the -tagging requirement on one of the two leading large-jets in the SR, and is denoted by CR0L-4Q and CR0L-2B2Q, respectively. The multĳet and noncollision background contributions are found to be negligible using the data-driven methods or from the estimation in a similar phase space carried out in Ref. [153].
The following subsections discuss the methodology and results of the irreducible and reducible background estimations.

Irreducible background estimation
The (¯+ ) events account for at most 10% of the total background in SR-4Q (SR-2B2Q), and are negligible in SR-2B2Q (SR-4Q). Given their minor contribution, these backgrounds are estimated directly from the MC predictions and assigned conservative uncertainties.
The dominant¯+ component in SR-2B2Q is¯(→ ) + (→ ). To validate the MC modeling, a dedicated validation region VRTTX is defined in a 3-lepton region populated by¯(→ ℓ ) + (→ ℓℓ). The selections are summarized in Table 3. Exactly three baseline and signal leptons are required, with the leading lepton of T > 30 GeV firing the single-lepton trigger. At least one large-jet is required in the event, and it must contain exactly two -tagged track-jets. No further kinematic cuts are applied, in order to maintain a sufficient data sample size in the region. The¯+ purity in this region is about 70%. Sixty-eight data events are observed in VRTTX, while 46.5 events are predicted by MC simulation (31.5 from¯+ , 12.1 from , and 2.9 from the others). A 70% uncertainty is therefore assigned to the¯+ normalization to fully cover the observed discrepancy; this has only a small impact on the total background estimation, given the minor contribution in the SRs.
The contribution to SR-4Q is mainly from (→ ) + (→ ) processes. The process has only recently been observed at the LHC [154,155]. As the data sample size in that phase-space region is insufficient and loosening the selection leads to poor purity, no control or validation regions were designed and this background is estimated directly from the MC prediction. A 50% uncertainty is assigned for the normalization based on the precision of the production cross-section measurement performed by the CMS experiment [155].

Reducible background estimation
The CR0L-4Q is defined by the same selection as SR-4Q-VV except that one of the two leading large-jets must fail the -tagging, and the CR0L-2B2Q is constructed as the logical union of SR-2B2Q-VZ and SR-2B2Q-Vh with the failing the -tagging. For the reducible backgrounds, the extrapolation from a CR0L to the SR is mainly characterized by the -tagging response for a fake boson jet. In order to maintain a sufficient data sample size and to suppress signal contamination, the eff selection and the mass window cut for CR0L-2B2Q are loosened relative to the SR, from eff > 1000 GeV to eff > 900 GeV and from 70 < ( ) < 135 GeV to 70 < ( ) < 150 GeV. The signal contamination in CR0L-4Q and CR0L-2B2Q is evaluated using the ( , )-SIM and ( , ) samples. While a contribution of up to 15% (24%) of the expected backgrounds in CR0L-4Q (CR0L-2B2Q) can be caused by the non-excluded signals, the introduced bias in the estimate is still smaller than the total uncertainty in the SRs and VRs, and therefore has only a small impact on the final sensitivity. The selections are summarized in Table 3, and the relation with the SR is illustrated in Figure 6.  Table 3 for details.
One advantage of this CR0L definition is that all physics processes contributing to the background ( +jets, +jets,¯, etc.) have comparable 0L transfer factors. This is because the fake boson jets have similar origins and kinematics, and therefore the same boson tagging efficiency, confirmed by the simulation. Consequently, they can be treated as a single combined component with a common normalization factor assigned to correct their normalizations at once. The normalization is performed in 4Q and 2B2Q separately with an independent normalization factor, based on the "background-only fits" described in Section 9.1.
Good MC modeling of the 0L transfer factor is essential for the estimation. Data events with exactly one lepton or one photon are utilized for the validation. This is motivated by the fact that in these regions the main backgrounds, (→ ) + jets in SR/CR0L, (→ ℓ ) + jets in the 1-lepton region and +jets in the 1-photon region, have similar ISR jet kinematics and, relative to the minor backgrounds, contribute similarly to the respective regions when a compatible kinematic phase space is chosen. The 1-lepton (1-photon) regions corresponding to CR0L and the SR are constructed, denoted by CR1L(1Y) and VR1L(1Y) respectively. The level of data-vs-MC agreement in the ratio VR1L(1Y)/CR1L(1Y) ("1L(1Y) transfer factor") is validated as a proxy for the MC modeling of the 0L transfer factor.
The selections applied for VR(CR)1L and VR(CR)1Y are listed in Table 3. The VR(CR)1L is defined by requiring exactly one baseline lepton and a signal lepton with T > 30 GeV that fires the single-lepton trigger. In addition, miss T > 50 GeV is required to suppress contributions from QCD multĳet events with poorly measured miss T . The VR(CR)1Y requires exactly one baseline and signal photon to fire the single-photon trigger, and vetoes events with at least one baseline lepton. For both VR(CR)1L and VR(CR)1Y, at least two large-jets are required, and the events are separated into the 4Q and 2B2Q categories, based on the absence or presence of , respectively A -jet veto strategy similar to that for the SRs and the CR0L bins is applied: unmatched -jet = 0 in all the regions, -jet ≤ 1 in the 4Q regions. A stricter veto of -jet = 0 is applied in CR(VR)1L-4Q to suppress the large single-top contribution relative to the other 4Q regions.
The kinematic selection in the 1L (1Y) regions is applied with miss T replaced by T ( ) ( T ( )) in the variables. This is to ensure that compatible phase spaces are probed for (→ ) + jets in SR/CR0L, (→ ℓ ) + jets in VR1L/CR1L, and +jets in VR1Y/CR1Y, since miss T in (→ ) + jets events typically represents the T of the boson. The T ( ) variable is defined as the vector sum of the T of the lepton and miss T in the 1L regions. The miss T , T ( ) and T ( ) are collectively denoted by T ( ). The min Δ ( miss T , ) variable, as well as those analogously derived from T ( ) and T ( ), are collectively denoted by min Δ ( , ). The cuts in T ( ) and eff are moderately loosened relative to the SRs to ensure a data sample size sufficient for the validation. For 2B2Q regions, the mass window cut for the /ℎ -tagging is loosened to 70 < < 150 GeV to match CR0L-2B2Q, further increasing the data sample size. The contributions from events containing fake leptons and fake photons are estimated by data-driven methods [156,157] and are found to be negligible.
The same estimation procedure used to estimate the backgrounds in the SRs is employed to predict the backgrounds in the VRs. A background-only fit, described in Section 9.1, is performed, where the sum of the reducible backgrounds in MC samples is normalized to the data in each CR1L (CR1Y). The modeling of the transfer factor is examined using the level of agreement between the data and the post-fit background expectation in the corresponding VR1L (VR1Y). Four sets of fits are performed in the 1L-4Q, 1L-2B2Q, 1Y-4Q and 1Y-2B2Q categories separately. Including all the systematic uncertainties discussed in Section 8, the observed data event yields and the post-fit background expectations in the VRs are summarized in Table 4, and compared graphically in Figure 7. The distributions of eff in VR1L(1Y)-4Q and T2 in VR1L(1Y)-2B2Q are shown in Figure 8. The data agree reasonably well with the estimated backgrounds across the VRs. The largest disagreement is found in VR1L-4Q, corresponding to a statistical significance of 1.8 .

Systematic uncertainties
Uncertainties in the expected signal and background yields account for the statistical uncertainties of the MC samples, the experimental systematic uncertainties associated with the detector measurements and reconstruction, and the theoretical systematic uncertainties in the MC simulation modeling. For the signals and irreducible backgrounds, the uncertainties are assigned directly to the event yields. For the reducible backgrounds, however, they are assigned to the SR(VR)/CR ratio ("transfer factor", TF) as a consequence of the normalization performed in the CRs. The reducible backgrounds are also subject to uncertainties due to the limited data sample size in the CRs. Each systematic uncertainty is treated as fully correlated across the analysis regions but not across physics processes, unless explicitly stated otherwise. A summary of the background prediction uncertainties is shown in Figure 9. The post-fit values are quoted after a background-only fit described in Section 9. The MC statistical uncertainties give the largest contribution to the systematic uncertainty, mainly from the limited size of the (→ ) + jets background sample used for the extrapolation. However, this is not a limiting factor for the analysis sensitivity since the total uncertainty in the SRs is dominated by the statistical uncertainty due to the low number of data events in Table 3: Summary of selections for the SRs, CRs, and VRs. ( ) ( (! )) represents the number of large-jets passing (failing) the -tagging of the two highest-T large-jets. The same selection is applied to the SR (VR) and CR in the same category except for the -tagging and some kinematic selections that are explicitly indicated in parentheses. The trigger selection and event cleaning described in Section 5 are also applied. VRTTX is a validation region used to validate the¯+ modeling as described in Section 7.1. T ( ) is the vector sum of the T of the lepton and miss T in the 1L regions. In the 1L (1Y) regions, miss T is replaced by T ( ) ( T ( )) when calculating the kinematic variables eff , min Δ ( miss T , ) and T2 . Details are given in Section 7.2.
the SRs. Details of the experimental and theoretical systematic uncertainties are described in the following subsections.

Experimental uncertainties
The first class of experimental uncertainties is related to the reconstruction and identification efficiencies for large-jets, small-jets, leptons and photons considered in the analysis. These are assigned as uncertainties in the efficiency correction factors applied to the MC samples, which correct for discrepancies between the efficiency predicted by MC simulation and the efficiency in data, as measured using dedicated control samples.
The uncertainty in ( )-tagging contributes the most to this class of efficiency uncertainties. It originates from the MC modeling uncertainty in the jet substructure variable distributions ( , 2 and trk ), as well as the precision of the efficiency determination in data. The MC modeling uncertainty is evaluated by comparing samples from different MC generator configurations [137]. The -tagging efficiency in data is measured using¯for signal jets (large-jets that contain the → decays) and QCD multĳet/ +jets for background jets, following the prescription described in Ref. [137]. The total uncertainty of the -tagging efficiency correction factor ranges from 12% to 23% for signal jets, and from 7% to 15% for background jets, depending on the T . The same efficiency correction factors are  applied to the -tagging selection, with a simulation-based additional uncertainty of 4%-5% to account for differences between and bosons.
The uncertainty in the (ℎ )-tagging is obtained from the -tagging efficiency uncertainty and the mass distribution's shape uncertainty. For the -tagging uncertainty, 1%-10%, 15%-50%, and 50%-100% uncertainties are assigned to the MC correction factor for -jets, -jet mistagging, and light-flavor jet mistagging, respectively, driven by theoretical uncertainties in the MC-simulated efficiency and the precision of the efficiency measurement in data. The mass scale uncertainty is estimated using the trk method [132], and a relative scale uncertainty of 2%-8% is assigned depending on T, and . A 20% relative uncertainty is assigned for the mass resolution based on the variation in the simulation [159].
Other efficiency uncertainties related to triggering, identification, reconstruction, and isolation requirements of electrons [75,141], muons [142] and photons [75], and the jet vertex tagger selection for small-jets, are found to be negligible.
The second class of experimental uncertainties is related to the energy (or momentum) determination for the reconstructed objects, namely large-jets [132], small-jets [145], electrons [141], muons [160] and photons [141]. These typically come from the precision of simulation-based and in situ calibrations of the energy (or momentum) scale and resolution. These per-object uncertainties are propagated through the miss T calculation, with additional uncertainties accounting for the scale and resolution of the soft term [146].
Additionally, an uncertainty in the integrated luminosity used to normalize the MC samples is considered. A 1.7% uncertainty is quoted for the combined 2015-2018 integrated luminosity obtained primarily using the Table 4: Number of observed data events and the post-fit SM background prediction in the VR1L (1Y) bins and the corresponding CR1L (1Y) bins. Negligibly small contributions are indicated by "-". Within each CR the reducible backgrounds have the same relative uncertainty in their expected yield because a common normalization factor is assigned to all of them in the fit.

Theoretical uncertainties
Theoretical uncertainties in the main reducible backgrounds ( / +jets, +jets and production) are estimated with varied generator parameters. Uncertainties due to the choice of QCD renormalization and factorization scales are evaluated by varying them up and down by a factor of two relative to their nominal values [163]. For / +jets, the uncertainties related to the choice of CKKW merging scale are also considered. These are assessed by shifting the merging scale to 15 GeV or 30 GeV from the default scale of 20 GeV. For the¯background, the nominal P +P 8 sample is compared with two alternative samples: one from M G 5_ MC@NLO to estimate the hard-scatter modeling uncertainty, and the other from P B interfaced to H 7.0.4 [164] and H7UE set of tuned parameters [164] to assess the uncertainty due to the choice of parton shower scheme and hadronization model. Variations in the2 initial-and final-state radiation modeling, and renormalization and factorization scales, are also considered following the prescription described in Ref. [165]. Uncertainties related to the choice of NNPDF3.0 PDF sets are assigned to the / +jets, +jets, ,¯and single-top backgrounds. These are derived by taking the envelope of the eigenvector variations from 100 propagated uncertainties.
For the reducible background estimation, an additional uncertainty is assigned for the modeling of the relative background composition. This is because different physics processes ( +jets, +jets,¯, etc.) are considered as a single component in the fits, and therefore its composition is predicted solely by the simulation. While the estimation is insensitive to the composition at first order since the TFs of those physics processes are similar, the residual TF difference can cause a bias in the estimation when the composition is significantly mismodeled by the simulation. The impact of the potential composition mismodeling is evaluated by the variation in the combined TF when shifting the normalization of each physics process up and down by a factor of two. An uncertainty of about 2%-8% is assigned in the 4Q category, while 7%-10% is assigned for those in the 2B2Q category.
For the irreducible backgrounds, cross-section uncertainties are assigned to account for their normalization.
A 50% uncertainty is quoted for production based on the cross-section measurement by the CMS experiment [155], while a 70% normalization uncertainty is assigned to¯+ based on the data/MC discrepancy observed in VRTTX, as discussed in Section 7.1.  Figure 9: The total post-fit uncertainty and the breakdown in each of the SRs and VRs. "MC statistical" and "CR statistical" present the statistical uncertainty due to the limited sample size either in MC simulation (mainly due to the (→ ) + jets sample size) or in the CR data, respectively. "Boson tagging" indicates the uncertainty in the ( )-tagging efficiency, while "experimental" shows the contribution from the rest of the experimental uncertainties. "Reducible BG composition" refers to the uncertainty in the reducible background estimation due to the relative composition uncertainty. "Theory" represents the total uncertainty from all the theoretical uncertainties.
The uncertainty in the signal yields consists of the cross-section uncertainty and the shape uncertainties. The cross-section uncertainty ranges from 6% to 20% for the production of electroweakinos with masses between 400 GeV and 1 TeV, driven mainly by the PDF uncertainty [125]. The shape uncertainties comprise uncertainties in the choice of renormalization/factorization scales and parton shower modeling, affecting the signal acceptance by 5%-10%.

Statistical analysis
Final background estimates are obtained by performing a profile log-likelihood fit [166] simultaneously in all CRs and SRs relevant to a given interpretation. The H F [167] framework is employed. Systematic uncertainties are treated as Gaussian-distributed nuisance parameters in the likelihood, while the statistical uncertainties of the MC samples are treated as Poisson-distributed nuisance parameters.
Three types of fit configurations are used to derive the results.
• A "background-only fit" is performed considering only the CRs and assuming no contribution from signals. The normalization of the total reducible background is allowed to float and is constrained by the fit using the data in the CRs. The normalization factors and nuisance parameters are adjusted by maximizing the likelihood. Three independent sets of fits are performed in the 0L, 1L, and 1Y categories, respectively. The 4Q and 2B2Q regions in each category are fitted simultaneously, but with independent normalization factors assigned in 4Q and 2B2Q. The normalization factors obtained from the fits range from 0.7 to 1.3.
• A "discovery fit" performs the hypothesis test for a generic beyond-the-SM (BSM) signal, setting upper limits on the number of events and visible cross-section for the signal. The fit uses only a single SR and the associated CR0L bin(s), constraining the backgrounds following the same method as in the background-only fit. Any contribution from signals is allowed only in the SR, and the signal-strength parameter is defined to be strictly positive.
• An "exclusion fit" is performed to set the exclusion limit for a given signal model. The SRs and the corresponding CR0L bins are fit simultaneously to determine the reducible background normalization factors and constrain the systematic uncertainties. The signal contamination in CR0L is also taken into account according to the model predictions.
For each discovery or exclusion fit, the compatibility of the observed data with the background-only or signal-plus-background hypotheses is quantified by calculating a one-sided 0 -value with the profile likelihood ratio used as a test statistic [166]. The upper limits and exclusion are derived using the CL s prescription [168] where the 95% confidence level (CL) exclusion is defined by CL s < 0.05.

Signal region yields
The observed data yields in each SR and CR0L together with their SM background expectations are summarized in Table 5 and shown in Figure 10. No significant excess is found in any of the SRs. The distributions of eff in SR-4Q-VV and T2 in SR-2B2Q-VZ and SR-2B2Q-Vh are shown in Figure 11 with some representative signal samples overlaid to illustrate the sensitivity.     (

Model-independent upper limits
A discovery fit is performed for each SR to derive the expected and observed 95% CL upper limits on the number of BSM signal events ( 95 exp and 95 obs ) as well as the one-sided -value ( 0 ) of the background-only hypothesis. Pseudo-experiments with toy MC are used for the calculation. An upper limit on the cross-section, 95 obs where represents the efficiency times acceptance of the SR for the given signal, is obtained by dividing 95 obs by the integrated luminosity. The upper limits and the 0 value associated with each SR are summarized in Table 6. Two additional "discovery SRs" are defined in order to set model-independent upper limits in the inclusive phase space. First, Disc-SR-2B2Q is defined as the logical union of SR-2B2Q-VZ and SR-2B2Q-Vh, and then the inclusive discovery signal region Disc-SR-Incl is defined as the logical union of SR-4Q-VV and Disc-SR-2B2Q. When evaluating Disc-SR-Incl, both CR0L-4Q and CR0L-2B2Q are included in the simultaneous fit and each has its own floating normalization factor for the reducible backgrounds.

Model-dependent exclusion limits
The results are also interpreted in the context of the specific signal models discussed in Section 2. An exclusion fit is performed for each point in the model space, and a CL s value is assigned based on the hypothesis test. The expected and observed 95% CL exclusion regions correspond to values of CL s < 0.05. Given the large number of models tested, an asymptotic approximation [166] is employed in the CL s calculation instead of the full calculation using pseudo-experiments. The validity is checked and the CL s values of the two methods typically agree within 5%, and maximally within 10%.
The SRs participating in the simultaneous fit vary with the signal model being tested, as summarized in Table 1. CR0L-4Q (CR0L-2B2Q) is included in the fit when at least one SR-4Q (SR-2B2Q) bin is used in deriving the limit, while both of the CR0L bins are included when the SR bins from both SR-4Q and SR-2B2Q participate in the fit. The ( , ) limits are also interpreted for the /ℎ-funnel dark matter model described in Section 2.1. The exclusion limits are shown in Figure 13, overlaid with theoretical predictions of higgsino masses ( 0 2 ) reproducing the observed dark matter relic density (Ωℎ 2 = 0.12) [55]. While the relic density depends drastically on tan , the exclusion limits obtained by the search are assumed to be constant along tan since they do not change when varying B ( 0 2 → 0 1 ), and are interpreted from the exclusion limits for the B ( 0 2 → 0 1 ) = 50% hypothesis in Figure 12(d). For the ℎ-funnel case, where ( 0 1 ) = ℎ /2, the excluded regions are tan > 8.5 for > 0 and 5.5 < tan < 7 for < 0.
The exclusion limits set on the ( , ) and ( , ) models are evaluated in a three-dimensional model space defined by ( 2 , , tan ). For each model point, the mass spectra and the branching ratios are determined using the prescription described in Section 2.1. Figure 14(a) shows the expected limit and observed limit as a function of ( 2 , ) with a fixed tan = 10. Figure 14(b) shows that varying tan or the sign of (sign( )) has very little effect on the sensitivity. The limits are also interpreted as a function of the physical electroweakino masses so that they can be directly compared with the other models. For a given set of (tan , sign( )), a pair of ( 2 , | |) can be projected one-to-one to ( ( ± 2 ), ( 0 1 )) when all the other MSSM parameters are fixed. Figure 14(c) and Figure 14(d) show the limits for the ( , ) model ( 2 > | |) and ( , ) model ( 2 < | |), respectively, assuming tan = 10 and > 0.
The mass exclusion limits are shown to be highly stable with respect to the internal variations within each model. The limits for ( , ) and ( , ) (or for ( , ) and ( , )) are also very similar, despite the branching ratios of heavy being substantially different between the models. This model dependency is small mostly due to the statistical combination of SR-4Q and SR-2B2Q, and the inclusive SR bins (SR-4Q-VV, SR-2B2Q-VZ, and SR-2B2Q-Vh), which are designed to be agnostic with regard to the difference between and bosons.   To summarize, a wino mass between 400 GeV and 1060 GeV is excluded for the wino production models for ( 0 1 ) < 400 GeV; and a higgsino mass between 450 GeV and 900 GeV is excluded for the higgsino production models for ( 0 1 ) < 240 GeV.

Exclusion limits on ( , ) simplified model: ( , )-SIM
The exclusion limits for the ( , ) simplified models (C1C1-WW, C1N2-WZ and C1N2-Wh) are also derived, in order to directly compare the search sensitivity with the previous ATLAS analyses. Figure 15 shows the obtained exclusion as a function of the produced wino mass ( ± 1 / 0 2 ) and the bino LSP mass ( 0 1 ). For C1C1-WW, ( ± 1 ) between 630 GeV and 760 GeV is excluded for ( 0 1 ) < 80 GeV. For C1N2-WZ (C1N2-Wh), ( ± 1 / 0 2 ) between 440 GeV and 960 GeV (400 GeV and 1060 GeV) is excluded for ( 0 1 ) < 300 GeV (420 GeV). 13 The sensitivity to high-mass winos is significantly improved relative to ATLAS searches using the other final states and the same dataset. For example, the expected limits on C1N2-WZ (C1N2-Wh) are typically extended by about 300 (140) GeV in ( ± 1 / 0 2 ), corresponding to exclusion of signals with a 7.5 (2.4) times smaller production cross-section than in the search using final states with three leptons [24] (one lepton and two -jets [29]). This result also sets the most stringent limit on the model to date from the LHC experiments.

Conclusion
A search for electroweakino pair production using final states consisting of miss T and two boosted hadronically decaying heavy SM bosons ( , , or ℎ) is reported. Signatures with large mass-splitting between the produced electroweakino and the lightest SUSY particle (LSP) are targeted.
The use of fully hadronic final states takes advantage of the large SM boson branching ratios, while the backgrounds are efficiently suppressed by reconstructing the , and ℎ bosons using boosted-bosontagging techniques. Using 139 fb −1 of proton-proton collision data at √ = 13 TeV recorded by the ATLAS detector at the LHC, this strategy provides unprecedented sensitivity to the production of heavy electroweakinos.
No excess over the SM background prediction is observed, and 95% CL exclusion limits are set for signal models in various -parity-conserving scenarios. For wino pair production with direct decays into a bino or higgsino LSP, wino masses between 400 GeV and 1060 GeV are excluded when the LSP mass is below 400 GeV and the mass splitting is larger than 400 GeV. For higgsino pair production with direct decays into a bino or wino LSP, higgsino masses between 450 GeV and 900 GeV are excluded when the LSP mass is below 240 GeV and the mass splitting is larger than 450 GeV. The limits are also examined for various wino (higgsino) branching ratio assumptions, by directly scanning over relevant branching ratios (for the bino LSP models) or over the MSSM parameters that dictate them ( 2 , , and tan for the wino/higgsino LSP models). The results are shown to be highly consistent for the variations.
The results are also interpreted in the context of simplified models of wino production with decays into a bino LSP, which are more conventionally explored in electroweakino searches at the LHC. For chargino pair production, with each decaying into a boson and an LSP, a chargino mass between 630 GeV and 760 GeV is excluded for an LSP mass below 80 GeV. For chargino-neutralino pair production involving decays into ( ℎ) and two LSPs, a wino mass between 440 (400) GeV and 960 (1060) GeV is excluded for an LSP mass below 300 (420) GeV. These extend significantly beyond the exclusion limits set by the previous searches at the LHC using different final states.
Finally, exclusion limits are set on higgsino production with decays into a massless gravitino LSP or a massless/massive axino LSP, motivated by the general gauge mediation model (GGM) or new physics models involving the axion, respectively. A higgsino mass between 450 (500) GeV and 940 (850)                    The ATLAS Collaboration