Search for supersymmetry with Higgs boson to diphoton decays using the razor variables at $\sqrt{s} = $ 13 TeV

An inclusive search for anomalous Higgs boson production in the diphoton decay channel and in association with at least one jet is presented, using LHC proton-proton collision data collected by the CMS experiment at a center-of-mass energy of 13 TeV and corresponding to an integrated luminosity of 35.9 fb$^{-1}$. The razor variables $M_\mathrm{R}$ and $\mathrm{R}^2$, as well as the momentum and mass resolution of the diphoton system, are used to categorize events into different search regions. The search result is interpreted in the context of strong and electroweak production of supersymmetric particles. We exclude bottom squark pair-production with masses below 450 GeV for bottom squarks decaying to a bottom quark, a Higgs boson, and the lightest supersymmetric particle (LSP) for LSP masses below 250 GeV. For wino-like chargino-neutralino production, we exclude charginos with mass below 170 GeV for LSP masses below 25 GeV. In the GMSB scenario, we exclude charginos with mass below 205 GeV for neutralinos decaying to a Higgs boson and a goldstino LSP with 100% branching fraction.


Introduction
The discovery of the Higgs boson [1][2][3], the first fundamental scalar particle ever observed, has opened a new window for exploring physics not described by the standard model (SM) of particle physics. Many models of physics beyond the standard model (BSM) postulate the existence of cascade decays of heavy states involving Higgs bosons [4,5]. In the minimal supersymmetric standard model (MSSM) [6], Higgs bosons may be produced in a variety of ways. The bottom squark, the supersymmetric partner of the bottom quark, produced via the strong interaction, may decay to a Higgs boson, quarks, and the lightest supersymmetric particle (LSP). Similarly charginos or neutralinos produced through the electroweak interaction may decay to a Higgs boson and the LSP. Of particular interest are scenarios with gauge-mediated supersymmetry breaking (GMSB), where the lightest neutralino may decay to a Higgs boson and the goldstino LSP ( G) [7,8]. The decay signature depends on whether the chargino and neutralino mixed states are dominated by the wino or higgsino components, the respective supersymmetric partners of the W and Higgs bosons. Diagrams of simplified models [9] for the scenarios considered are shown in Fig. 1. In this paper, we denote the Higgs boson as H to indicate that it is the particle observed by the ATLAS and CMS experiments. In the MSSM, this particle is typically assumed to correspond to the lighter of the two CP-even Higgs particles and is often denoted as h. For the GMSB scenario, we consider simplified models where Higgsinolike charginos and neutralinos are nearly mass-degenerate and both chargino-neutralino and neutralino-pair production result in very similar final state signatures, and are hereafter collectively referred to as chargino-neutralino production in this paper. These examples of BSM production of Higgs bosons motivate an inclusive search for anomalous Higgs boson production that is broadly sensitive to a wide range of supersymmetric (SUSY) scenarios. Similar searches for supersymmetric particles decaying to Higgs bosons have been performed by the ATLAS and CMS collaborations in the past using 8 TeV collision data and can be found in references [10][11][12].
In this Letter, we present an updated search for supersymmetry events with at least one Higgs boson candidate decaying to two photons produced in association with at least one jet produced in 13 TeV proton-proton collisions. The data were collected by the CMS experiment and correspond to an integrated luminosity of 35.9 fb −1 [13]. The diphoton decay mode of the Higgs boson provides a good compromise between branching fraction and background rejection. The transverse momentum of the Higgs boson candidate, the expected mass resolution, and the razor variables M R and R 2 [14,15], explained in detail in Section 4, are used to define event categories which generically enhance BSM signals relative to SM backgrounds. The potential signal is extracted from the dominant nonresonant multijet background through a fit to the diphoton mass distribution. The results of the search are interpreted in terms of simplified models of bottom squark pair production and chargino-neutralino production.

The CMS detector, trigger, and event reconstruction
The central feature of the CMS detector 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. Extensive forward calorimetry complements the coverage provided by the barrel and endcap detectors. Muons are measured in gas-ionization detectors embedded in the magnet 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 2 The CMS detector, trigger, and event reconstruction Figure 1: Diagrams displaying the simplified models that are being considered. Upper left: bottom squark pair production; upper right: wino-like chargino-neutralino production; bottom: the two relevant decay modes for higgsino-like neutralino pair production in the GMSB scenario. Ref. [16].
Signal event candidates are recorded using a diphoton trigger, requiring the transverse energy of the leading and subleading photons to be larger than 30 GeV and 18 GeV, respectively and their invariant mass to be larger than 90 GeV. Additional requirements on the photon shower shape and isolation are imposed to reduce the background rate and improve the signal purity [17]. The efficiency of the trigger with respect to events passing the offline selection is measured to be above 98%.
Physics object candidates are reconstructed using a global event description based on the CMS particle-flow (PF) algorithm [18], which identifies particles through an optimized combination of information from the various detector subsystems. Photon candidates are selected by imposing "loose" requirements on the shower shape in the electromagnetic calorimeter, the ratio of energy measured in the HCAL to the energy measured in the ECAL, and the isolation in a cone around the direction of the photon momentum [19]. The isolation requirement is corrected for the effect of multiple proton collisions in the same or adjacent bunch crossing (pileup) by subtracting the average energy from pileup deposited in the isolation cone, estimated by averaging the energy density over the event. Furthermore, photon candidates are rejected if they match an electron candidate that is not consistent with one leg of a conversion. The photon selection efficiency was measured to be about 90% [20] using tag and probe methods. The measured energy of photons is corrected for clustering and local geometric effects using an energy regression trained on Monte Carlo (MC) simulation [19]. This regression gives a significant improvement in the energy resolution of the photons (about 30%) and provides an estimate of the uncertainty of the energy measurement that is used to separate events into high-and low-resolution categories.
The reconstructed vertex with the largest value of summed physics-object p 2 T is taken to be the primary pp interaction vertex (PV). The physics objects used in this context are the objects returned by a jet finding algorithm [21,22] applied to all charged tracks associated with the vertex under consideration, plus the corresponding associated missing transverse momentum.
The charged PF candidates associated with the PV and the neutral PF candidates are clustered into jets using the anti-k T algorithm [21] with a distance parameter R = 0.4, as implemented in the FASTJET package [22]. The jet momentum is determined as the vectorial sum of all particle momenta in the jet. Jet energy corrections are derived based on a combination of simulation studies, accounting for the nonlinear detector response and the presence of pileup, together with in-situ measurements of the energy balance in dijet and γ+jet events using the methods described in Ref. [23]. For this analysis, jets with |η| < 3.0 that do not overlap with any identified photon are selected by requiring ∆R = √ (∆η) 2 + (∆φ) 2 > 0.5 between photon and jet candidates. The combined secondary vertex (CSVv2) tagging algorithm [24] is used to identify jets originating from the hadronization of b quarks. A loose working point is chosen that yields an efficiency above 80% and a mistag rate for light-flavor jets that is approximately 10%. The negative vector sum of the reconstructed p T of all PF candidates in an event defines the missing transverse momentum p miss T in the event, and its magnitude is referred to as p miss T . Events with detector-and beam-related noise that can mimic event topologies with high energy and large p miss

Event simulation
Simulated event samples are used to model the SM Higgs backgrounds in the search regions, and to calculate the selection efficiencies for SUSY signal models. Samples of SM Higgs production via gluon fusion, vector boson fusion, associated production with a W or a Z boson, bbH, and ttH are generated using the next-to-leading order (NLO) MADGRAPH5 aMC@NLO v2.2.2 [28] event generator. The Higgs mass is assumed to be 125.0 GeV for the simulated event samples and is within the uncertainty of the currently best measured value [29]. For the gluon fusion production mode, the sample is generated with up to two extra partons to model initial-state radiation (ISR) calculated at the matrix element level with NLO accuracy with the matching scheme described in reference [30]. The SUSY signal MC samples are generated using MAD-GRAPH5 aMC@NLO at leading order accuracy with up to two extra partons in the matrix element calculations with the matching scheme described in reference [31]. In both cases, PYTHIA v8.2 [32] is used to model the fragmentation and parton showering with the CUETP8M1 tune [33]. The NNPDF3.0LO and NNPDF3.0NLO [34] parton distribution functions (PDFs) are used for the LO and NLO accuracy generators respectively. The SM Higgs background and bottom squark pair-production signal samples are simulated using a GEANT4-based model [35] of the CMS detector, while the chargino-neutralino and neutralino-pair production signal samples are simulated with the CMS fast simulation package [36]. While generally providing an accurate description, the fast simulation does sometimes yield inaccurate predictions of the missing transverse momentum tail. These inaccuracies are accounted for by larger systematic uncertainties for the signal yield prediction in the relevant phase space estimated as the difference between signal yields predicted using the generator level missing transverse momentum and the missing transverse momentum reconstructed based on the fast simulation. All simulated events include the effects of pileup and are processed with the same chain of reconstruction programs used for collision data.
To improve the MADGRAPH modeling of ISR in the SUSY signal MC samples, we apply a correction as a function of the multiplicity of ISR jets for sbottom pair production, and as a function of the transverse momentum (p ISR T ) of the chargino-neutralino system for chargino-neutralino production, derived from studies of tt and Z+jets events, respectively. The correction factors vary between 0.92 and 0.51 for ISR jet multiplicity between one and six, and between 1.18 and 0.78 for p ISR T between 125 and 600 GeV. The correction has a small effect on the signal yields at the level of about 1%. The full size of this correction is taken as a systematic uncertainty. The Higgs production cross sections are obtained from the recommendations of the LHC Yellow Report 4 of the LHC Higgs Cross Section Working Group [37]. The SUSY signal production cross sections are calculated to NLO plus next-to-leading logarithmic (NLL) accuracy [38][39][40][41][42][43], assuming all SUSY particles other than those in the relevant diagram to be too heavy to participate in the interaction. These NLO+NLL cross sections and their associated uncertainties [43] are used to derive the exclusion limits on the masses of SUSY particles.

Event selection and search categories
We select events with two photons that satisfy the selection criteria described above. Both photons must be in the barrel region of the electromagnetic calorimeter, with |η| < 1.44, and have p T > 20 GeV. At least one of the two photons must have p T > 40 GeV. If multiple photon pairs are identified, the pair with the largest scalar sum of the transverse momenta is chosen as the Higgs boson candidate in the event. The Higgs boson candidate mass is required to be between 103 GeV and 160 GeV in order to cover a sufficiently large background dominated sideband region.
The Higgs boson candidate and any additional identified jets with p T > 30 GeV and |η| < 3.0 are clustered into two hemispheres (megajets) according to the Razor megajet algorithm [15], which minimizes the sum of the squared-invariant-mass values of the two megajets. To converge, the algorithm requires at least one such identified jet in the event. Next, the razor variables [14] M R and R 2 are computed as follows: where p is the momentum of a megajet, p z is its longitudinal component, and j 1 and j 2 are used to label the two megajets. In the definition of R 2 , the variable M R T is defined as: The razor variables M R and R 2 provide discrimination between SUSY signal models and SM background processes with SUSY signals typically having large values of M R and R 2 , while the SM background exhibits an exponentially falling spectrum in both variables.
The selected events are separated into four mutually exclusive categories. An event is categorized as "HighPt" if the transverse momentum of the selected Higgs boson candidate is larger than 110 GeV. Otherwise, if the event contains two b-tagged jets whose invariant mass is between 76 and 106 GeV, or between 110 and 140 GeV, it is categorized as "H(γγ)-HZ(bb)". The remaining events are categorized as "HighRes" and "LowRes" if the diphoton mass resolution estimate σ M /M is smaller or larger than 0.85%, respectively, where σ M is computed as: where E γ1,2 is the energy of each photon and σ Eγ1,2 is the estimated energy resolution for each photon. A graphical representation of the categorization procedure is shown in Fig. 2.
The "HighPt" category isolates SUSY events producing high-p T Higgs bosons; the "H(γγ)-HZ(bb)" category isolates SUSY signals that produce two Higgs bosons or a Higgs boson and a Z boson in the final state; and the HighRes and LowRes categories further improve the discrimination between signal and background in the remaining event sample. The "H(γγ)-HZ(bb)" category combines events with two Higgs bosons or a Higgs boson and a Z boson in order to achieve a sufficiently large number of events in the sideband for the background estimation method described in Section 5 to remain unbiased.
Each event category is further divided into bins in the M R and R 2 variables, which define the exclusive search regions. A significant excess of events above the SM expectation in one or several bins would provide evidence of BSM physics. The search regions are chosen based on an optimization procedure that maximizes the expected sensitivity to the simplified bottom squark pair production model discussed further in Section 7, and are summarized in Table 1.
The bins in the M R and R 2 variables are kept identical for the HighRes and LowRes categories to allow for simultaneous signal extraction, since the ratio of the event yields in these two categories does not depend on the details of the signal model.

Background prediction
There are two main classes of background events that pass the search selection criteria: SM Higgs production and nonresonant photon production, with either two promptly produced photons or one prompt photon and one jet that is wrongly identified as a photon. The SM Higgs background is estimated from the MC simulation, while the nonresonant background prediction is estimated using a fit to the diphoton mass distribution observed in data.
Within each search bin, we extract a potential signal by performing an unbinned extended maximum likelihood fit to the diphoton mass spectrum. An example of such a fit is shown in Fig. 3. The nonresonant background is modeled with a decaying functional form given in Table 1 for each individual search region bin. All parameters of the function are unconstrained in the fit. The functional form of the model used for each search region bin is selected on the basis of its Akaike information criterion (AIC) score [45], which quantifies the trade-off between goodness-of-fit and model complexity. Each functional form is tested for fit biases 6 5 Background prediction Table 1: A summary of the search region bins in each category is presented. The functional form used to model the nonresonant background is also listed. An exponential function of the form e −am γγ is denoted as "single-exp"; a linear combination of two independent exponential functions of the form e −am γγ and e −bm γγ is denoted as "two-exp"; a modified exponential function of the form e −am b γγ is denoted as "mod-exp"; and a Bernstein polynomial of degree n [44] is denoted by "poly-n". The bin labels 9-13 are used for both the HighRes and LowRes categories because the data in these categories are always fitted simultaneously with potentially different nonresonant background models used. Further details on the simultaneous fit are discussed in Section 5. with respect to a set of alternative models, all of which adequately describe the data in the diphoton mass sideband region (103-121 GeV and 129-160 GeV). The shape of the Higgs boson resonance from SM Higgs production and from decays of SUSY signals is modeled with a double Crystal Ball function [46,47] with two independent tail parameters that is fitted to the diphoton mass distribution obtained from the MC simulation. The parameters of each double Crystal Ball function are held constant in the signal extraction fit procedure, with the exception of the parameter controlling the location of the peak, which is discussed further in Section 6 below. The normalization of the SM Higgs boson background in each bin is predicted from the MC simulation and is constrained to that value in the fit within uncertainties. For the HighRes and LowRes categories, bins in the M R and R 2 variables are fitted simultaneously. For a given search bin, the relative yields in the HighRes and LowRes categories are observed in the simulation to be largely process independent and are therefore constrained according to the simulation prediction. Based on these independent fits in each search bin, we obtain a modelindependent search result, which can be used to evaluate whether the yield in any bin exhibits a statistically significant deviation from the background prediction.
We also perform a combined simultaneous fit using all of the search bins, to test specific SUSY simplified model signal hypotheses. In the combined fit, the yield in each bin for the SM Higgs background and the signal model under test are constrained to the MC simulation predictions within uncertainties. These uncertainties are modeled by use of nuisance parameters that account for various theoretical and instrumental uncertainties that can affect the SM Higgs boson background and SUSY signal normalization, and are profiled in the fit. A more detailed discussion of systematic uncertainties can be found below in Section 6. The MC simulation predictions for the SM Higgs boson background normalization are shown in Table 2 for each bin in the search region.

Systematic uncertainties
There are broadly two types of systematic uncertainties. The first and dominant systematic uncertainty is in the shape and normalization of the nonresonant background. This is propagated by profiling the normalization and shape parameters of the nonresonant background functional form in an unconstrained way. The second and subdominant type of systematic uncertainty is in the predictions of the SM Higgs background in the various search bins. These shape uncertainties are propagated through the use of several independent nuisance parameters, where both theoretical and instrumental effects are taken into account. The nuisance parameters are constrained in the fit using log-normal prior functions, whose width reflects the size of the systematic uncertainty. The independent systematic effects considered include missing higher-order corrections, PDFs, trigger and selection efficiencies, jet energy scale uncertainties, b tagging efficiencies, and the uncertainty in the integrated luminosity. The uncertainty due to jet energy resolution uncertainties were also estimated and were found to be negligible. The typical size of these effects on the expected limit is summarized in Table 3. Due to effects of pileup and transparency loss in the ECAL crystals, we observe some simulation mismodeling of the estimated mass resolution, which results in a systematic uncertainty of 10-24% in the prediction of the SM Higgs background and SUSY signal yields in the HighRes and LowRes event categories. The systematic uncertainty in the photon energy scale is implemented as a Gaussian-distributed nuisance parameter that shifts the Higgs boson mass peak position, constrained in the fit to lie within approximately 1% of the nominal Higgs boson mass observed in simulation. The systematic uncertainty for the modeling of the ISR for the signal process is also propagated and is below 1%. For chargino-neutralino and neutralino-pair production signal processes, the fast simulation was used to predict signal yields and an additional systematic uncertainty is propagated for inaccuracies in the modeling of the missing transverse momentum tail. This systematic uncertainty ranges between 1% and 34% depending on the search region bin.

Results and interpretations
The fit results for all search region bins are summarized in Table 4, along with the data yields, fitted background, and signal yields. An example fit result for the search bin with M R > 600 GeV and R 2 > 0.025 in the HighPt category is shown in Fig. 3. The observed signal significance in each bin is summarized in Fig. 4 for all the search region bins, which are statistically independent. None of the 14 bins exhibits a deviation from the background expectation larger than two standard deviations.
We interpret the search results in terms of limits on the production cross section times branching ratio for simplified models of bottom squark pair-production and chargino-neutralino production. Diagrams of these simplified models are shown in Fig. 1. In the case of bottom squark pair production, we consider the scenario where the bottom squark decays to a bottom quark and the next-to-lightest neutralino ( χ 0 2 ), and the χ 0 2 decays to a Higgs boson and the LSP ( χ 0 1 ), and the production cross sections are computed at NLO plus next-to-leading-log (NLL) precision with all the other sparticles assumed to be heavy and decoupled [38][39][40][41][42][43]. We restrict ourselves to the scenario where the mass splitting between the χ 0 2 and the χ 0 1 is 130 GeV, slightly above threshold to produce an on-shell Higgs boson. In the case of chargino-neutralino production, we consider two different scenarios. In the first one, pure wino-like charginos ( χ ± 1 ) [GeV]   Figure 4: The observed significance in units of standard deviations is plotted for each search bin. The significance is computed using the profile likelihood, where the sign reflects whether an excess (positive sign) or deficit (negative sign) is observed. The categories that the bins belong to are labeled at the bottom. The bins in the HighRes and LowRes categories are fitted simultaneously and yield a single combined significance. The yellow and green bands represent the ±1 and ±2 standard deviation regions, respectively. and the next-to-lightest neutralino χ 0 2 are mass-degenerate and are produced together, with the chargino decaying to a W boson and the LSP ( χ 0 1 ) and the χ 0 2 decaying to a Higgs boson and the LSP ( χ 0 1 ). The production cross sections are computed at NLO plus next-to-leading-log (NLL) precision in a limit of mass-degenerate wino χ 0 2 and χ ± 1 , light bino χ 0 1 , and with all the other sparticles assumed to be heavy and decoupled [48][49][50]. In the second scenario, we consider a GMSB [7, 8] simplified model where Higgsino-like charginos and neutralinos are nearly massdegenerate and are produced in pairs through the following combinations: and χ ± 1 χ ∓ 1 . Because of the mass degeneracy, both the χ 0 2 and the χ ± 1 will decay to χ 0 1 and other low-p T (soft) particles, leading to a signature with a χ 0 1 pair. Each χ 0 1 will subsequently decay to a Higgs boson and the goldstino ( G), which is the LSP, or to a Z boson and the goldstino. We consider the case where the branching fraction of the χ 0 1 → H G decay is 100%, and the case where the branching fraction of the χ 0 1 → H G and χ 0 1 → Z G decays are each 50%. The cross sections for higgsino pair production are computed at NLO plus NLL precision in a limit of mass-degenerate higgsino χ 0 2 , χ ± 1 , and χ 0 1 with all the other sparticles assumed to be heavy and decoupled [48][49][50]. Following the convention of real mixing matrices and signed neutralino or chargino masses [51], we set the mass of χ 0 1 ( χ 0 2 ) to positive (negative) values. The product of the third and fourth elements of the corresponding rows of the neutralino mixing matrix N is +0.5 (−0.5). The elements U 12 and V 12 of the chargino mixing matrices are set to 1.
Following the CL s procedure [52][53][54], we use the profile likelihood ratio test statistic and the asymptotic formula [55] to evaluate the 95% confidence level (CL) observed and expected limits on the signal production cross sections. For the bottom squark pair production model, the limits are shown on the left of Fig. 5 as a function of the bottom squark mass and the LSP mass. We exclude bottom squarks with masses below about 450 GeV for all LSP masses below 250 GeV. For the wino-like chargino-neutralino production simplified model, the limits are shown on the right of Fig. 5 as a function of the chargino mass and the LSP mass. We exclude chargino masses below about 170 GeV for all LSP masses below 25 GeV. For the higgsino-like chargino-neutralino production simplified models, the limits are shown in Fig. 6 as a function of the chargino mass for the case where the branching fraction of the χ 0 1 → H G decay is 100% on the left, and for the case where the branching fraction of the χ 0 1 → H G and χ 0 1 → Z G decays are both 50%, on the right. We exclude charginos below 205 GeV and 130 GeV in the former and latter cases, respectively.

Summary
A search for anomalous Higgs boson production through decays of supersymmetric particles is performed with the proton-proton collision data collected in 2016 by the CMS experiment at the LHC. The sample corresponds to an integrated luminosity of 35.9 fb −1 at the center-ofmass energy √ s = 13 TeV. Higgs boson candidates are reconstructed from pairs of photons in the central part of the detector. The razor variables M R and R 2 are used to suppress Standard Model (SM) Higgs boson production and other SM backgrounds. The non-resonant background is estimated through a fit to the diphoton mass distribution in data, while the SM Higgs background is predicted using simulation. We interpret the results in terms of production cross section limits on simplified models of bottom squark pair production and chargino-neutralino production. We exclude bottom squark masses below 450 GeV for bottom squarks decaying to a bottom quark, a Higgs boson, and the lightest supersymmetric particle (LSP) for LSP masses below 250 GeV and assuming a mass splitting between the χ 0 2 and χ 0 1 of 130 GeV. For wino-like chargino-neutralino production, we improved the search sensitivity by a factor of two with respect to previous results [11] and we exclude charginos with mass below 170 GeV for LSP masses below 25 GeV. In the GMSB scenario, we exclude charginos with mass below 205 GeV for neutralinos decaying to a Higgs boson and a goldstino LSP ( G) with 100% branching fraction. Finally, we exclude charginos with mass below 130 GeV for the case where the branching fractions of the χ 0 1 → H G and χ 0 1 → Z G decays are 50% each.

Acknowledgments
We congratulate our colleagues in the CERN accelerator departments for the excellent performance of the LHC and thank the technical and administrative staffs at CERN and at other CMS institutes for their contributions to the success of the CMS effort. In addition, we gratefully acknowledge the computing centers and personnel of the Worldwide LHC Computing Grid for delivering so effectively the computing infrastructure essential to our analyses. Finally, we acknowledge the enduring support for the construction and operation of the LHC and the CMS detector provided by the following funding agencies: BMWFW and FWF (Aus  [10] ATLAS Collaboration, "Search for direct pair production of a chargino and a neutralino decaying to the 125 GeV Higgs boson in √ s = 8 TeV pp collisions with the ATLAS detector", Eur. Phys. J. C 75 (2015), no. 5, 208, doi:10.1140/epjc/s10052-015-3408-7, arXiv:1501.07110.