Search for a new resonance decaying into two spin-0 bosons in a final state with two photons and two bottom quarks in proton-proton collisions at √ s = 13 TeV

A search for a new boson X is presented using CERN LHC proton-proton collision data collected by the CMS experiment at √ s = 13 TeV in 2016–2018, and corresponding to an integrated luminosity of 138 fb − 1 . The resonance X decays into either a pair of Higgs bosons HH of mass 125 GeV or an H and a new spin-0 boson Y. One H subsequently decays to a pair of photons, and the second H or Y, to a pair of bottom quarks. The explored mass ranges of X are 260–1000 GeV and 300–1000 GeV, for decays to HH and to HY, respectively, with the Y mass range being 90–800 GeV. For a spin-0 X hypothesis, the 95% confidence level upper limit on the product of its production cross section and decay branching fraction is observed to be within 0.90–0.04 fb, depending on the masses of X and Y. The largest deviation from the background-only hypothesis with a local (global) significance of 3.8 (below 2.8) standard deviations is observed for X and Y masses of 650 and 90 GeV, respectively. The limits are interpreted using several models of new physics. Published


Introduction
With the discovery of a Higgs boson H in 2012 at the CERN LHC by the ATLAS and CMS Collaborations [1][2][3], the standard model (SM) has proven to describe most observed phenomena up to the TeV scale [4][5][6][7][8][9].However, some of the important aspects of nature are left outside of this theory, such as dark matter and the nature of gravity.Beyond-SM (BSM) theories have been proposed to address these shortcomings, and many of these predict the existence of new massive spin-0 and -2 particles.
Among the BSM theories, models with "warped extra dimensions" (WEDs) [10,11] postulate the existence of an additional spatial dimension in which the carriers of a quantized gravitational force (gravitons) propagate.In addition, according to the "Randall-Sundrum bulk model", matter particles are also allowed to propagate in the extra dimension.This model contains massive resonances, the spin-0 radion [12][13][14] and spin-2 bulk gravitons (excited modes of the gravitational field each separated by a characteristic mass scale) [15][16][17], which have sizable branching fractions for decaying into HH [18].For this model, the fundamental scale is the reduced Planck scale M pl with ultra-violet cutoff Λ R [12].In the specific range of model parameters considered for this analysis, a narrow-width approximation is permitted for the resonances [18].Supersymmetry (SUSY) is an extension of the SM which postulates a fermion-boson pairing symmetry (superpartners) [19,20].The minimal supersymmetric SM (MSSM) extends the SM Lagrangian by introducing an extra complex scalar doublet [21,22].The next-to-minimal supersymmetric SM (NMSSM) [23,24] introduces one additional complex singlet field, and is the most straightforward supersymmetric extension to the SM where the electroweak scale originates from the SUSY-breaking scale.The NMSSM could provide a natural solution to the "unnaturalness" problem of MSSM [25], wherein the difference between the m H and the Planck scale is of many orders of magnitude.In the NMSSM Higgs-sector, several spin-0 bosons exist, one of which can be associated with the H discovered so far.More massive spin-0 bosons may decay to the lighter ones in this model.Similar to the NMSSM, a model with two real scalar singlet fields (TRSM), in addition to the SM Higgs field, has been proposed [26].The mixing between the Higgs doublet and these additional scalar fields gives three massive scalar bosons, one of which can be identified as the SM-like H.This model, too, allows the heavier scalar to decay into lighter scalars.Motivated by the above BSM scenarios, this paper presents the search for a massive particle X decaying either to HH or to H and another spin-0 boson Y, using LHC proton-proton (pp) collision data collected by the CMS experiment in 2016-2018 at √ s = 13 TeV and corresponding to a total integrated luminosity of 138 fb −1 .The results are interpreted within the framework of the BSMs presented above.
The search is performed in the final state with a pair of photons γ and bottom b quarks (γγbb).Figure 1 shows a tree-level Feynman diagram of the signal process.In using this channel, the X → HH search benefits from the combination of the H → γγ decay, with its relatively large signal purity, and the H → bb decay, with its large branching fraction.For the H, the SM decay branching fractions [27] are assumed.The HY search is independent of any assumption on the Y → bb branching fraction.
The analysis uses a strategy similar to the search for nonresonant production of HH in the γγbb final state [28,29].One H candidate is reconstructed from a pair of photons γγ, with their invariant mass m γγ compatible with m H .The candidate for the second boson, either H or Y, is reconstructed from a pair of b jets (particle jets originating from the fragmentation This paper is organized as follows: Section 2 provides a brief description of the CMS detector followed by the details of the data and simulated event samples in Section 3. The analysis strategy is discussed in Section 4, including background rejection methods along with the signal and background modeling studies.Section 5 details the systematic uncertainties.The results are presented in Section 6, and the analysis is summarized in Section 7. Tabulated results are provided in the HEPData record for this analysis [30].

The CMS detector
The central feature of the CMS apparatus is a superconducting solenoid of 6 m internal diameter, providing a magnetic field of 3.8 T. Within the solenoid volume, there is a silicon pixel and strip tracker, a lead tungstate crystal electromagnetic calorimeter (ECAL), and a brass and scintillator hadron calorimeter (HCAL), each composed of a barrel and two endcap sections.Forward calorimeters extend the pseudorapidity η coverage provided by the barrel and endcap detectors.Muons are measured in gas-ionization detectors embedded in the steel flux-return yoke outside the solenoid.A more detailed description of the CMS detector, together with a definition of the coordinate system used and the relevant kinematic variables, can be found in Ref. [31].
Events of interest are selected using a two-tiered trigger system [32].The first level (L1), composed of custom hardware processors, uses information from the calorimeters and muon detectors to select events at a rate of around 100 kHz within a fixed latency of about 4 µs [33].The second level, known as the high-level trigger, consists of a farm of processors running a version of the full event reconstruction software optimized for fast processing and reduces the event rate to around 1 kHz before data storage [34].
The global event reconstruction, also called particle-flow (PF) event reconstruction [35], aims to reconstruct and identify each individual particle in an event (PF candidates), with an optimized combination of all subdetector information.In this process, the identification of the particle type (photon, electron, muon, charged or neutral hadron) plays an important role in the determination of the particle direction and energy.
Photons are identified as ECAL energy clusters not linked to the extrapolation of any charged particle trajectory to the ECAL.Electrons are identified as a primary charged particle track and potentially many ECAL energy clusters corresponding to this track extrapolation to the ECAL and to possible bremsstrahlung photons emitted along the way through the tracker material.Muons are identified as tracks in the central tracker consistent with either a track or several hits in the muon system, and associated with calorimeter deposits compatible with the muon hypothesis.Charged hadrons are identified as charged particle tracks neither identified as electrons, nor as muons.Finally, neutral hadrons are identified as HCAL energy clusters not linked to any charged hadron trajectory, or as a combined ECAL and HCAL energy excess with respect to the expected charged hadron energy deposit.
The energy of photons is obtained from the ECAL measurement.The raw ECAL energy is corrected in data using a multivariate analysis (MVA) based regression which is trained on Monte Carlo (MC) simulation of Z → ee events.In the simulated events, the energy of the photons is smeared to match the resolution in data [36,37].For photon identification, an MVA identification method (photon ID) based on photon shower shape, isolation, and kinematic variables is used [36].
The energy of electrons is determined from a combination of the track momentum at the main interaction vertex, the corresponding ECAL cluster energy, and the energy sum of all bremsstrahlung photons attached to the track [38].The energy of muons is obtained from the corresponding track momentum.The energy of charged hadrons is determined from a combination of the track momentum and the corresponding ECAL and HCAL energies, corrected for the response function of the calorimeters to hadronic showers.Finally, the energy of neutral hadrons is obtained from the corresponding corrected ECAL and HCAL energies.In the barrel section of the ECAL, an energy resolution of about 1% is achieved for unconverted or late-converting photons in the tens of GeV energy range.The energy resolution of the remaining barrel photons is about 1.3% up to |η| = 1, changing to about 2.5% at |η| = 1.4.In the endcaps, the energy resolution is about 2.5% for unconverted or late-converting photons, and between 3 and 4% for the other ones [39].The m γγ resolution, as measured in H → γγ decays, is typically in the 1-2% range, depending on the measurement of the photon energies in the ECAL and the topology of the photons in the event [40].
For each event, jets are clustered from these reconstructed particles using the infrared and collinear safe anti-k T algorithm [41,42] with a distance parameter of 0.4.Jet momentum is determined as the vectorial sum of all particle momenta in the jet, and is found from MC simulations to be, on average, within 5-10% of the true value over the whole transverse momentum p T spectrum and detector acceptance.Additional pp interactions within the same or nearby bunch crossings (pileup) can contribute additional tracks and calorimetric energy depositions to the jet momentum.To mitigate this effect, charged particles identified to be originating from pileup vertices are discarded and an offset correction is applied for the remaining neutral particle contributions [43].Corrections, derived from MC simulations, are applied in order to bring the measured jet energies to that of the originating particles, on average.Furthermore, in situ measurements of the momentum balance in dijet, γ+jet, Z+jet, and multijet events are used to account for any residual differences in the jet energy scale between data and MC simulations [44].An additional energy correction algorithm based on machine learning techniques is applied in this analysis to further improve the jet energy scale of b jets [45].
The jet energy resolution amounts typically to 15-20% at 30 GeV, 10% at 100 GeV, and 5% at 1 TeV [44].Additional selection criteria are applied to each jet to remove jets potentially dominated by anomalous contributions from various subdetector components or reconstruction failures.
The missing transverse momentum vector ⃗ p miss T is computed as the negative vector p T sum of all the PF candidates in an event, and its magnitude is denoted as p miss T [46].The ⃗ p miss T is modified to account for the corrections to the energy scale of the reconstructed jets in the event.

Data and simulated samples
The data used in this analysis were collected by the CMS detector from LHC pp collisions at √ s = 13 TeV in 2016-2018.The total analyzed data correspond to an integrated luminosity of 138 fb −1 [47][48][49].Events are collected using the trigger condition for the year 2016 (2017-2018) that requires two photons with thresholds p γ 1 T > 30 GeV and p γ 2 T > 18 (22) GeV and with m γγ > 90 GeV, where γ 1 and γ 2 represent leading-and subleading-p T photons, respectively.The photons in the trigger algorithm have also to satisfy requirements on their calorimeter shower shapes, on the ratio of their energies deposited in the HCAL to that in the ECAL (identification variable), and on the p T carried by charged hadrons in the vicinity of the photon candidate (isolation requirement) [36].
MC simulations of resonant production of X via gluon-gluon fusion (ggF), with subsequent decays, assuming a narrow-width approximation, to either HH or HY were generated at nextto-leading order (NLO) using the MADGRAPH5 aMC@NLO package [50,51].For the WEDmotivated searches for spin-0 and -2 resonances, the MADGRAPH5 aMC@NLO package versions were 2.2.2 and 2.4.2, for the 2016 and 2017-2018 simulated samples, respectively.The signals were generated for m X in range 260-1000 GeV.For the NMSSM-motivated searches, the signal samples were generated using MADGRAPH5 aMC@NLO 2.6.5 for 300 < m X < 1000 GeV and for 90 < m Y < 800 GeV.The samples were generated in steps of 100 GeV in m X , and the kinematic mass distributions for intermediate m X were modeled by interpolating the shape parameters of the above simulated samples.Within signal simulations, a value of 125 GeV is assumed for m H .
The upper bound for m X search at 1 TeV is defined by the transition toward the Lorentz-boosted regime, where particles resulting from the hadronization of both b quarks merge into a single jet.The upper range of m Y is set by the energy-momentum conservation rule for the on-shell particles m Y ≤ m X − m H . Below m Y = 90 GeV, the search capabilities are limited by the position of the kinematic turn-on in the m jj distribution [28,29].
The dominant backgrounds include the SM nonresonant multijet processes with up to two prompt photons, either as the irreducible γγ+jets events or the reducible γ+jets events.The γγ+jets background is modeled with SHERPA 2.2.1 [52] at leading order (LO) and includes up to three additional partons at the matrix element level.The γ+jets background is modeled with PYTHIA 8.212 [53] at LO.The nonresonant background simulations are used only in the optimization and validation of the background estimation method using an MVA-based discriminator.
Single H productions, where the H decays to γγ, are considered as resonant backgrounds.These background events give a signal-like resonant distribution for m γγ , which is one of the observables in the analysis.The ggF H is simulated at NLO using POWHEG 2.0 [54][55][56][57].The vector-boson fusion production (VBF H) and production in association with top quark pairs (ttH) are simulated at NLO using MADGRAPH5 aMC@NLO v2.2.2 (2016) or v2.4.2 (2017-2018).The H production in association with a vector boson (VH) and in association with b quarks (bbH) [58] are simulated using MADGRAPH5 aMC@NLO v2.6.1.The cross sections and decay branching fractions are taken from Ref. [27].The SM nonresonant HH production has a negligible contribution within the analysis phase-space and is not explicitly considered as a separate background source.

Analysis strategy 4.1 Object reconstruction and event preselections
Photons are considered if they are contained within the ECAL and tracker fiducial region of |η| < 2.5, excluding the ECAL barrel-endcap transition region (1.44 < |η| < 1.57), and they also need to pass the selections p T /m γγ > 1/4, and 100 < m γγ < 180 GeV.The primary pp vertex associated with the diphoton system is identified using an MVA technique [68].The efficiency of the correct vertex assignment is 99.9%, aided by the requirement of at least two jets in the γγbb final state.
The jets are required to be within the tracker fiducial volume, so that track information can be used to identify heavy flavor (b and c) jets.These jets must have p T > 25 GeV and |η| < 2.4 (< 2.5) for 2016 (2017-2018), while also passing jet identification criteria, in order to reject spurious jets from calorimeter noise [69].The |η| selection of jets is looser for the 2017-2018 data-taking years because of the installation of the new pixel detector at the beginning of 2017 [70].Then, the jets and photons are required to be separated by ∆R ≡ (∆η) 2 + (∆ϕ) 2 > 0.4, where ∆η and ∆ϕ are the pseudorapidity and the azimuthal angle (measured in radians) differences between a photon and the closest jet, respectively.This selection is used in order to avoid counting of photons as jets.Finally, jet pairs with 70 < m jj < 190 (1200) GeV for HH (HY) analysis are retained for further selection, respectively.
For the identification of b jets, a deep neural network (DNN) score-based tagger DEEPJET [71][72][73] is used.It utilizes information related to tracks and secondary vertex of b jets.The pair with the highest sum of DEEPJET discriminant scores in the event is selected as a H (Y) candidate.If no such pair is present, the event is rejected.The physics object reconstruction and event preselection criteria are summarized in Table 1.The signal efficiency varies from 30-60%.For HH searches, it increases with m X .In HY searches, for a fixed m Y the signal efficiency increases with m X , but, for a fixed m X it decreases with m Y .
The X candidate is reconstructed from the selected pairs of photons and jets.Figure 2 shows the m γγ , m jj , and 4-body invariant mass m γγjj distributions in the data and MC simulations after the full set of selection criteria.The best estimate of the mass of the X candidate is given by the variable M X , which has a better resolution than m γγjj [74] because of the resolution subtraction of m γγ and m jj , defined as: The improvement in resolution by using M X instead of m γγjj ranges from 30% at high m X up to 90% at low m X .

Nonresonant background separation
The dominant nonresonant backgrounds stem from the SM multijet processes with up to two prompt photons.A boosted decision tree (BDT) discriminant is trained using MC simulations to separate signal and nonresonant backgrounds γγ+jets and γ+jets.For this purpose, a gradient boosting algorithm, XGBOOST [75] is used to develop a multiclass BDT classifier.
Various resonant m X and m Y hypotheses with different kinematic distributions are studied within the analysis.The phase space is subdivided into three regions for m X (< 500, 500-700, and > 700 GeV) and for m Y (< 300, 300-500, and > 500 GeV).These ranges are combined into nine m X -m Y domains, three of which are forbidden by the energy conservation constraint m Y < m X − m H .The mass ranges are decided according to the kinematic properties and Lorentzboost factor, estimated by the quantity m X /(m H + m Y ) [76], which defines similar angular properties.One BDT is trained per m X -m Y region.For each mass range in m Y , a dedicated m jj interval is defined: 70-400, 150-560, and 300-1000 GeV, respectively.Each interval is optimized to have the complete m jj shape of the signals lying within one mass range to use as input for signal extraction.
For the purpose of BDT training, all simulated X → HY signal events that fit into a domain are combined into one sample with equal weights to build the signal sample.The simulated nonresonant background events are used as background samples, weighted by the products of their cross sections and the integrated luminosity.Before the training, the signal and background samples are normalized to unity.For HH searches, the training with m jj interval 70-400 GeV is used since it already covers the m jj interval of  GeV.
The signal and background events are intermixed and divided into train and test sets independently for every training.The input hyperparameters of the training are optimized using the 5-fold cross-validation [77] to mitigate over-training.An early-stopping feature of XGBOOST is also used to control the over-training.The BDT is trained for the three data-taking years together.The training strategy is validated by comparing the per year performance from the merged training with the performance obtained by the separate training for each year.The BDT training input variables, based on the two photons and two jets associated with the X decay, are motivated by Ref. [29].They are grouped as follows: 1. Discriminating kinematic variables: • Helicity angles (| cos θ CS HY |, |cos θ jj |, |cos θ γγ |), here CS refers to the Collins-Soper frame [78].
• p T (γ)/m γγ and p T (j)/m jj for leading and subleading photons and jets, respectively.

Object identification variables:
• Leading-and subleading-p T photon ID discriminant to reject misidentified photon contribution • Leading-and subleading-p T jet b tagging discriminant to reject light jets 3. Resolution estimators: • Leading-and subleading-p T photons: Energy resolution (σ E /E) and mass resolution (σ m γ γ /m γγ ); • Leading-and subleading-p T jets: p T resolution (σ p T (j) /p T (j)) and mass resolution (σ m jj /m jj ).
In addition to this, the transverse momentum deposits per unit area of the η-ϕ plane coming from pileup of the event, also known as the pileup transverse momentum density, is added in the inputs of the training to include pileup information which is different in the data-taking periods 2016-2018.The final BDT training output discriminates the signal from the background events by populating the signal events towards high BDT score and background events towards low BDT score.To preserve the signal statistics when defining optimal signal-enriched regions, the BDT outputs are transformed, ensuring that the combined signal BDT output distribution in each m X -m Y domain becomes uniform.It should be noted that, for a given value of m X and m Y , the BDT distribution may exhibit some departure from flatness, especially at very high BDT scores.We have verified for all signal hypotheses that this transformation, made for convenience, does not affect the analysis performance.The distribution of one of the BDT outputs for data and simulated events is shown in Fig. 3.

Event classification
We classify events into three categories based on BDT output along with a selection on M X .The optimization of the BDT categories and M X selection is performed independently aiming for the best signal sensitivity using only the simulations.
The categorization remains the same for the signals within an m X -m Y domain.A constraint is imposed on the background statistics, requiring a minimum of 10 events within each category during the optimization of BDT boundaries.This minimum number of events in category optimization is chosen so that the systematic uncertainties from the background modeling do not exceed the benefits of a better signal over background ratio.Notably, this constraint excludes events falling within the m γγ = [115, 135] GeV region.This approach ensures a sufficient number of background events in the m γγ side-bands for each category, to guarantee a robust background modeling.We use signal and nonresonant background simulations to define categories based on the Punzi figure of merit [79] defined as ϵ s /(1 + √ B), where ϵ s is signal efficiency and B refers to background yield.The category optimization remains independent of assuming any signal cross section since the figure of merit only uses signal efficiency, which prevents Table 2: The BDT based event classification according to defined m X (and m Y ) ranges for HH (and HY) searches.For HH searches, the column with m Y < 300 GeV is used.The number represents the BDT scores showing a decreasing signal purity region from category 0 to 2. categorization bias towards any specific signal hypothesis.The categories for different m X -m Y domains are shown in Table 2.The categories are labeled with 0 to 2 from the highest to lowest purity.For HH searches, categorization in m Y < 300 GeV region is used.The signal efficiency gets reduced by up to a factor of 2 after BDT categorization.
The M X selection differs depending upon m X .For the HH analysis, a selection on M X retaining at least 60% of the selected signal events is defined separately for each m X hypothesis as shown in Fig. 4.This selection remains the same for HH and HY signals because, for low m Y , the M X -resolution of HY signals is similar to HH signals.In contrast, for high m Y , M X holds better resolution, giving higher signal efficiency with the same selection.

Resonant background rejection
The resonant backgrounds originate from the production of single H. Their contribution is strongly reduced by a tight selection applied on the M X for each m X hypothesis as stated in Section 4.3.Nevertheless, for m X < 550 GeV, the presence of ttH, containing genuine b jets from top quark decays, remains significant.A DNN-based discriminant was developed to further reduce it for the nonresonant HH → γγbb analysis [29].Using the same working points for the DNN as in Ref. [29] improved the limits by up to 10%.For m X > 550 GeV, the single H contribution in the total background is negligible, therefore, no selection on the discriminant is applied.Also, it is checked using the MC simulations that this selection does not distort the background mass distributions.

Signal and background modeling
For signal extraction, simultaneous fits are conducted on the m γγ and m jj distributions.Initially, separate parametric fits are executed for each category on these distributions.Subsequently, the final signal and background models are constructed by taking the product of the independent m γγ and m jj models.This independence was verified for signal, background simulations, and data by studying the possible correlations using the 2D distributions of the m γγ and m jj observables which were found to be negligible [29].For signal and resonant backgrounds, MC simulations are used for m γγ and m jj modeling, while the m γγ and m jj distributions in data are used to choose optimal analytical background model.The m γγ fit range remains same for all the signal and background models.The m jj fit range differs depending upon the signal hypothesis.The following m jj fit ranges are used:  GeV in HH searches and for m Y ≤ 150 GeV in HY searches; 70-400 GeV for m Y ∈ (150, 300) GeV; 150-560 GeV for m Y ∈ [300, 500] GeV; 300-1000 GeV for m Y > 500 GeV.During the final fit procedure, the shape parameters are allowed to float within the statistical constraints provided by the MC events used to build them.
For signal modeling, categorized events are fitted with a product of two parametric signal models: A sum of Gaussian distributions for m γγ and a double-sided Crystal Ball (CB) function [80] or the sum of a CB and a Gaussian function for m jj .The m γγ distribution is parametrized using the sum of up to five Gaussian functions without any constraint to have a common mean.
To model the resonant single H backgrounds, the m γγ distribution is constructed from the MC simulations following the same methodology used for the signal m γγ model.The m jj modeling depends on the production mechanism of the single H process, and a parameterization is obtained from the simulated distributions: for the ggF H and VBF H processes, the m jj distribution is modeled with a Bernstein polynomial; for VH production, a CB function is used to model the distribution of the hadronic decays of vector bosons; for ttH, where the two b jets are produced from a top quark decay, a Gaussian function with a mean value of 120 GeV is used.These backgrounds are negligible for m X > 550 GeV, therefore, they are absorbed within the nonresonant background model coming from data to simplify the signal extraction procedure.
The nonresonant background model is extracted from data using the discrete profiling method described in Refs.[68,81].This method makes use of polynomial and exponential functions to decide the analytic functions to fit the background m γγ and m jj distributions.It estimates the systematic uncertainty associated with these functions and treats the choice of background function as a discrete nuisance parameter in the likelihood fit to the data.For background modeling, the fit functions are optimized on data where the events from the signal region with 115 < m γγ < 135 GeV are not taken into account.Considering the negligible correlation between m γγ and m jj , a set of MC pseudo-experiments were generated with and without injecting fake signal and fitted using the factorized 2D model.Negligible bias, as defined in Ref. [2], has been found which validates the robustness of the background model [82].

Systematic uncertainties
The impact of several sources of systematic uncertainties is considered in the analysis.The nonresonant background model is constructed using a data-driven method, and the uncertainties associated with the choice of background fit function are taken into account by the discrete profiling method described in Section 4.5.The systematic uncertainties can affect both the shapes and the normalizations of the m γγ and m jj distributions of the signal and the resonant single H backgrounds.Systematic uncertainties modifying the shape of the m γγ and m jj distributions are built into the models as parametric nuisance parameters.
The theoretical and experimental systematic uncertainties described below impact the normalization of the signal and single H backgrounds.The uncertainties in the cross sections arising from the quantum chromodynamics (QCD) renormalization and factorization scale variations, and PDF uncertainties are considered only for the single H processes.The uncertainties in the branching fractions of the H decay are considered for both the single H backgrounds and signals.Among these, the uncertainty from the QCD scale variations is the dominant one.The uncertainty is computed by varying both scales independently up and down by a factor of 2, while disregarding the two extreme combinations of variation.The theoretical uncertainties are treated as fully correlated across the three data-taking years.
The experimental uncertainties modifying the normalization are studied separately for each category, and they are mostly treated as uncorrelated across the different data-taking years.The dominant sources among them are: • Jet energy scale and resolution corrections: These are measured using the method of balancing of jet p T in dijet, γ+jet, Z+jet, and multijet events [44,83].The uncertainties in these corrections vary from 1-2% depending upon the p T and η of jets.Their impact on the normalization is evaluated by varying the corrections within their uncertainties and propagating them to the final result.
• Jet b tagging: The DEEPJET discriminator distributions in MC simulations are reshaped by data-to-simulation scale factors [71].The impact from different systematic uncertainty sources on the reshaping is studied.The highest impact comes from the misidentification of light flavor jets as b jets.The average magnitude of the uncertainty in b tagging efficiency varies from 4-6%.
• Integrated luminosity: The integrated luminosities for the 2016, 2017, and 2018 datataking years have individual uncertainties of 1.2-2.5% [47][48][49], while the overall uncertainty in the measurement for the 2016-2018 period is 1.6%.The possible correlations from different common sources of luminosity measurements are taken into account.
• Trigger efficiency: This efficiency is calculated with the "tag-and-probe" method using Z → ee events [84].The average size of uncertainty in the efficiency is 1-2% and the impact on the limits is less than 1%.An additional uncertainty is also added to include the trigger inefficiency arising from a gradual shift in the timing of the inputs of the ECAL L1 trigger for |η| > 2.0 in 2016-2017 data-taking years [33].However, it does not impact the results.
• Photon preselection uncertainty: The efficiency after photon preselections is measured in simulations and data on Z → ee events using the "tag-and-probe" method [84].
The ratio of the efficiencies from data and MC simulations is used to derive the correction factor and corresponding uncertainties.Overall, it has an impact of less than 1% on the limits.
• Photon energy scale and resolution: These uncertainties associated with the corrections are estimated using Z → ee events and applied to data for the scale corrections and to MC simulations for resolution corrections [39].The average size of uncertainty is 1-2% with less than 1% impact on the limits.
The impact of other sources of systematic uncertainty is small.The impact of the total systematic uncertainty on the limits for the m X < 500 GeV and m Y < 250 GeV is 1-2%.For other mass ranges, this impact is less than 1%.The dominant experimental systematic uncertainty is associated with the b tagging scale factors.

Results
A likelihood function is defined using signal and background analytic models of m γγ and m jj distributions, optimized as described in Section 4.5, with nuisance parameters related to the uncertainties explained in Section 5. A simultaneous unbinned 2D maximum likelihood fit to the m γγ and m jj distributions is performed on three event categories to extract the signal.
Figure 5 shows the data and signal-plus-background fit for m γγ and m jj observables in CAT 0 defined in Table 2 and with M X selection for two signal hypotheses from HH and HY searches.
For HH searches, no statistically significant deviations in data are observed with respect to the background-only hypothesis.For HY searches, m γγ and m jj distributions in data are shown in the lower panel of Fig. 5 for the signal hypothesis where the largest excess is observed for m X = 650 GeV and m Y = 90 GeV with a local significance of 3.8 σ.
In order to account for the "Look elsewhere effect" [85], a global significance is calculated in the partial phase space of m X within the range 300-1000 GeV and m Y within the range 90-150 GeV, which amounts to 2.8 σ, with a corresponding p-value of 0.0027.This partial phase space is sufficient to provide an upper limit of the global significance, the true value of which would be lower if the full phase space is used.In other words, the chosen phase space is sufficient to demonstrate that the excess is not globally significant.
The observed excess being not statistically significant, we set upper limits at 95% confidence level (CL) on the product of resonant production cross section and branching fraction to γγbb channel using the CL s criterion [86], taking the LHC profile likelihood ratio as a test statistic [87].The asymptotic approximation is used in the limit setting procedure [86,88,89].
Figure 6 shows the upper limits on resonant production cross section as a function of resonance mass m X for HH searches.The expected limits decrease from low to high m X region and observed limits are consistent with them within ±2 standard deviations (±2 σ).Depending on m X , the observed upper limits at 95% CL vary between 0.82-0.07and 0.78-0.06fb for spin-0 and -2 resonant HH searches, respectively.The corresponding expected ranges are 0.74-0.08 and 0.65-0.06fb, respectively.Following the theoretical cross sections from Ref. [18], these results exclude masses up to 600 GeV for spin-0 bulk radion signal at Λ R = 6 TeV and up to 850 GeV for spin-2 bulk KK graviton signal with coupling factor of κ/M pl = 0.5.
Figure 7 describes the upper limits at 95% CL for the pp → X → HY → γγbb signal as a function of m Y with different m X hypotheses.For every value of m X , the expected limits first improve and then worsen with increasing m Y .Indeed, when m Y increases, the signal boost factor decreases, which leads to a drop in the signal selection efficiency.The observed limits remain consistent with the expected limits within ±2 σ of deviation, except for m X = 650 GeV and m Y ≤ 100 GeV.The observed and expected limits vary between 0.90-0.04 and 0.79-0.05fb depending upon the mass ranges in m X and m Y , respectively.Expected and observed 95% CL upper limit on the product of resonant production cross section and branching fraction for spin-0 (upper) and spin-2 (lower) pp → X → HH → γγbb signal hypotheses.The dashed and solid black lines represent expected and observed limits, respectively.The green and yellow bands represent the ±1 and ±2 standard deviations for the expected limit.The red lines show the theoretical predictions with different energy scales and couplings. -1

fb
Figure 7: Expected and observed 95% CL exclusion limit on production cross section for pp → X → HY → γγbb signal.The dashed and solid black lines represent expected and observed limits, respectively.The green and yellow bands represent the ±1 and ±2 standard deviations for the expected limit.The middle plot in the 3rd row shows the largest excess observed for m X = 650 GeV and m Y = 90 GeV.
This search loses its sensitivity above the specified m X range since the bbbb final state becomes more sensitive for HH searches, taking advantage of its larger branching fraction [90].With the improved object identification methods and analysis strategy, the search shows up to 25% improvement on the HH limits when performed with 2016 data only in comparison to the ones previously published by the CMS Collaboration [28].This search reports best-to-date limits for the HH searches in the region with m X < 800 GeV.Furthermore, the observed (expected) limits are better than those reported by the ATLAS Collaboration for HH searches in the γγbb channel using the 2016-2018 data by up to 40% (30%) [91].For HY searches, this is the first result in the γγbb final state, where the analysis reports the best-to-date sensitivity within the m X range of 300-1000 GeV and the m Y range of 200-800 GeV.
The best-fit value of the product of the cross section with the branching fraction for the decay into γγbb at the largest excess is σ(pp → X → HY → γγbb ) = (0.35 ± 0.17 0.13 ) fb.Given that no additional treatment is applied to the single H backgrounds for m X > 550 GeV and they are assumed to be included within the nonresonant background model extracted from data due to their minimal contribution (see Section 4.4 and 4.5), a study was conducted to examine the impact of this simplification of the background modeling on the observed local excess.The study involves quantifying resonant and nonresonant background models separately as a part of the final background model.The results indicate a reduction in the local significance of 8% in units of standard deviation, and 6% in the best-fit value of the cross section.This effect is well within the bias range accepted in this analysis following the Ref. [2].
Figure 8 compares the expected and observed upper limits to the product of the theoretical signal production cross section and the branching fraction pp → X → HY → γγbb.The theoretical predictions are defined as the maximally allowed cross sections for a given m X and m Y hypothesis within the NMSSM framework [92].It is provided by the LHC Higgs Working Group and derived from the codes NMSSMTOOLS 5.5.0 [93] and NMSSMCALC [94].The predictions integrate constraints from measurements of H properties, SUSY searches, B-meson physics, and dark matter searches.For this model, both expected and observed limits are below the theoretical cross sections for some of the signal hypotheses, leading to the exclusion of masses between 400-650 GeV in m X and 90-300 GeV in m Y .

Summary
A search for a new resonance X, decaying either to a pair of Higgs bosons HH or to an H and a new spin-0 boson Y, is presented.The search uses data from proton-proton collisions collected by the CMS experiment at the Large Hadron Collider in 2016-2018 at a center-of-mass energy of 13 TeV, corresponding to 138 fb −1 of integrated luminosity.The search targets beyond standard model particles as predicted by several models of new physics.For X decaying to HH, an m X range of 260-1000 GeV is covered, while for X decaying to HY, the search range is 300-1000 GeV in m X and 90-800 GeV in m Y .Results are presented as the upper limits at 95% confidence level on the product of the production cross section of X and its branching fraction to the γγbb final state, through either HH or HY decays.Depending upon the mass range, the observed limits for a spin-0 resonance X decaying to HH range from 0.82-0.07fb, while the expected limits are 0.74-0.08fb.For X decaying to HY, the observed limits are 0.90-0.04fb, while the expected limits lie in the range 0.79-0.05fb, depending on the masses m X and m Y .The data are found to be compatible with the standard model predictions over the searched domains.The largest deviation from the background-only hypothesis, with a local significance of 3.8 standard deviations, is observed for m X = 650 GeV and m Y =90 GeV.However, this excess is not considered significant as it is reduced to 2.8 standard deviations when a partial look-elsewhere effect is considered.The HY search is performed for the first time in the γγbb channel.The limits from the HH search are the most stringent to date for m X less than 800 GeV.

Figure 1 :
Figure 1: Feynman diagram showing a tree-level gluon-gluon fusion production of a BSM resonance X decaying to a pair of spin-0 bosons (HH or HY), which then decay to the γγbb final state.and hadronization of b quarks), with their invariant mass m jj compatible with m H or with an assumed value of the Y boson mass m Y .Events are selected using machine learning algorithms to suppress background contributions.The signal is extracted from a two-dimensional (2D) maximum likelihood fit spanning the (m γγ ,m jj ) plane.

Figure 2 :
Figure 2: The m γγ (upper left), m jj (upper right) and m γγjj (lower) distributions in data and MC simulations.The signal distributions, shown for different values of m X and m Y with an assumption of 1 fb cross section, have been scaled up by a factor of 10 3 .

Figure 3 :
Figure3: The distribution of the BDT output in data and MC simulations of the single H background processes and a signal, in m X = 500-700 GeV and m Y < 300 GeV range.The signal distribution, with an assumption of 1 fb cross section, has been scaled up by a factor of 10 3 .

Figure 4 :
Figure 4: M X selection for each m X in HH and HY signals.The red and green lines represent the upper and lower boundary of this selection.

Figure 5 :
Figure 5: Invariant mass distributions m γγ (left) and m jj(right)  with the data events (black markers), with M X selection corresponding to an HH signal with m X = 400 GeV (upper panel), and to an HY signal with m X = 650 GeV and m Y = 90 GeV (lower panel).The distributions are shown for the high signal purity category (CAT 0).The red dashed line shows the sum of the fitted signal and background events.The solid black line represents the total background component, including nonresonant and resonant background events.The green and yellow bands represent the ±1and ±2-standard deviations which include the uncertainties in fit to the background-only hypothesis.The lower panel in each plot shows the residual signal yield after the background subtraction.

Figure 6 :
Figure 6:  Expected and observed 95% CL upper limit on the product of resonant production cross section and branching fraction for spin-0 (upper) and spin-2 (lower) pp → X → HH → γγbb signal hypotheses.The dashed and solid black lines represent expected and observed limits, respectively.The green and yellow bands represent the ±1 and ±2 standard deviations for the expected limit.The red lines show the theoretical predictions with different energy scales and couplings.

Figure 8 :
Figure 8: Comparison of the expected (left) and observed (right) limits at 95% CL with the maximally allowed cross sections from the NMSSM model where the area within the red contours indicate the excluded mass regions.The limits are displayed as two-dimensional binned distributions.