Search for a massive scalar resonance decaying to a light scalar and a Higgs boson in the four b quarks final state with boosted topology

We search for new massive scalar particles X and Y through the resonant process X $\to$ YH $\to$ $\mathrm{b\bar{b}b\bar{b}}$, where H is the standard model Higgs boson. Data from CERN LHC proton-proton collisions are used, collected at a centre-of-mass energy of 13 TeV in 2016-2018 and corresponding to an integrated luminosity of 138 fb$^{-1}$. The search is performed in mass ranges of 0.9-4 TeV for X and 60-600 GeV for Y, where both Y and H are reconstructed as Lorentz-boosted single large-area jets. The results are interpreted in the context of the next-to-minimal supersymmetric standard model and also in an extension of the standard model with two additional singlet scalar fields. The 95% confidence level upper limits for the production cross section vary between 0.1 and 150 fb depending on the X and Y masses, and represent a significant improvement over results from previous searches.


Introduction
The discovery of a Higgs boson (H) of mass 125 GeV [1][2][3] at the CERN LHC validated the Brout-Englert-Higgs mechanism [4][5][6][7][8][9] of the standard model (SM), yet raised questions of its viability at higher energy scales [10][11][12][13]. Besides, empirical observations such as the measurements of the neutrino masses and the baryon asymmetry in the universe are inconsistent with SM expectations. Beyond the standard model (BSM) theories, including those invoking supersymmetry [14] or extra dimensions [15], seek to address many of the shortcomings of the SM. No BSM phenomena have been observed at the LHC. However, there are unexplored BSM parameter spaces, among which are areas of the scalar sector, the topic of this search.
The minimal supersymmetric extension of the SM (MSSM) [16,17] postulates two complex scalar field doublets with SU(2) gauge symmetry. The next-to-minimal model (NMSSM) [18,19] was proposed to solve the MSSM's "unnaturalness problem" [20], where the Higgs boson mass parameter is many orders of magnitude smaller than the Planck scale. In the NMSSM, an extra complex scalar field gives a total of seven Higgs bosons: three neutral (one would be associated with H) and two charged scalar particles, as well as two neutral pseudoscalars. Searches for a heavier scalar X decaying to SM particles [21-23] have set a lower limit on its mass at M X = 1.5 (1.0) TeV for tan β = 21 (8) [21], where tan β is the ratio of the vacuum expectation values (VEVs) of the two Higgs doublets. The NMSSM favours low tan β, where the current M X bounds are the weakest.
In the NMSSM, the neutral scalar production cross sections may be suppressed because of their small couplings to SM fermions [18]. Enhanced "Higgs-to-Higgs" decays are then possible, such as X → YH [24,25], Y being the lighter scalar. Within the NMSSM, the largest branching fractions for both H and Y (for Y mass M Y less than twice that of the top quark t) are to a b quark-antiquark pair, giving the final state X → YH → bbbb. For higher M Y values the Y → bb branching fraction is ∼10% [25,26]. The second dominant process is X → YH → ττbb, which has been excluded [27] for 0.4 < M X < 0.6 TeV and 50 < M Y < 200 GeV, for specific values of the parameters of the model.
Another interesting model of new physics that motivates this search is the two-real-scalarsinglet extension of the SM (TRSM) [28], which introduces two additional scalar fields. This simplified model, onto which more complicated theories can be mapped, has nine degrees of freedom: the masses and VEVs of the three scalar fields, and three mixing angles. In the scenario where all three VEVs are non-zero, the three fields give rise to three massive scalars, one of which can be associated with the H boson. Depending on their masses and mixing angles, the heaviest scalar can decay to the two lighter scalars. These in turn can decay to SM particles with mass-dependent branching fractions. This Letter describes the search for two new scalar particles, X and Y, the former being more massive and decaying through X → YH. The search uses LHC proton-proton (pp) collision data collected by the CMS experiment in 2016-2018 corresponding to a total integrated luminosity of 138 fb −1 [29][30][31]. The masses of the scalar particles satisfy M X > M Y + M H ; Y may be lighter or heavier than H and both Y and H decay to bb. The search is generic, and X and Y can be associated with the particles predicted in the NMSSM or the TRSM, which are both mentioned above.
This search focuses on the kinematic region where M X is sufficiently larger than both M Y and M H such that Y and H carry considerable momenta and therefore their decay products, i.e. the bb pairs, are highly collimated. We explore the mass ranges 0.9 < M X < 4 TeV and 60 < M Y < 600 GeV, complementing the X → YH → ττbb search [27]. In the highmomentum kinematic regime, special techniques are used to reconstruct the final states containing the collimated bb pairs, in order to increase the signal sensitivity well beyond that covered by previous searches [27].
Tabulated results for this analysis are provided in HEPData [32].

The CMS detector and event reconstruction
The central feature of the CMS apparatus is a superconducting solenoid of 6 m internal diameter, providing a magnetic field of 3.8 T. Within the solenoid volume are a silicon pixel and strip tracker, a lead tungstate crystal electromagnetic calorimeter (ECAL), and a brass and scintillator hadron calorimeter (HCAL), each composed of a barrel and two endcap sections. Forward calorimeters extend the pseudorapidity (η) coverage provided by the barrel and endcap detectors. Muons are 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. [33]. Events of interest are selected using a two-tiered trigger system. The first level, composed of custom hardware processors, uses information from the calorimeters and muon detectors to select events at a rate of around 100 kHz within a fixed latency of about 4 µs [34]. The second level, known as the high-level trigger (HLT), 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 [35].
The primary vertex is taken to be the vertex corresponding to the hardest scattering in the event, evaluated using tracking information alone, as described in Section 9.4.1 of Ref. [36].
A particle-flow algorithm (PF) [37] aims to reconstruct and identify each individual particle in an event, with an optimized combination of information from the various elements of the CMS detector. The photon energy is obtained from the ECAL measurements. The energy of electrons is determined from a combination of the electron momentum at the primary interaction vertex as determined by the tracker, the energy of the corresponding ECAL cluster, and the energy sum of all bremsstrahlung photons spatially compatible with originating from the electron track. The energy of muons is obtained from the curvature of the corresponding track. The energy of charged hadrons is determined from a combination of their momentum measured in the tracker and the matching ECAL and HCAL energy deposits, corrected for 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.
Jets are clustered from PF candidates using the anti-k T algorithm [38,39] with a distance parameter of either 0.4 (AK4 jets) or 0.8 (AK8 jets). The jet momentum is defined as the vectorial sum of all particle momenta in a jet, and is found from simulation to be, on average, within 5-10% of the true momentum over the whole transverse momentum (p T ) spectrum and detector acceptance [40]. Additional pp interactions within the same or nearby bunch crossings (pileup), averaging 23-32 in 2016-2018, can contribute additional tracks and calorimetric energy depositions to the jet momentum. The effect of pileup is mitigated using the charged-hadron subtraction (CHS) algorithm [41], whereby charged particles identified to be originating from pileup vertices are discarded and an offset correction is applied to correct for remaining contributions. Jet energy corrections are derived from simulation and data to bring the measured response of jets to that of particle level jets on average [40].
For AK8 jets, masses are computed after applying grooming [42] techniques, which remove wide-angled soft and collinear radiation from the jets, in order to mitigate the effects of contamination from initial state radiation, the underlying event, and multiple hadron scattering. The trimming algorithm [43] uses a subjet size parameter of 0.3 and a radiation fraction parameter z = 0.1, which determines the minimum p T fraction that the reclustered jet constituents need to have in order not to be removed. The mass of the resultant jet is referred to as its "trimmed mass". The "soft-drop mass" of the jet is obtained by applying the soft-drop algorithm [44,45]. Here it is obtained using a value z = 0.1 for the radiation fraction parameter. The angular exponent parameter is set as β = 0, so there is no dependence of the p T fraction threshold on the distance between jet constituents.
In case of the soft-drop algorithm, the pileup per particle identification (PUPPI) [41,46] algorithm is used to mitigate the effect of pileup on AK8 jets. In PUPPI, the treatment of charged particles is similar to that in CHS. A weight between zero and one is assigned to neutral particles, larger values indicating higher probability of originating from the primary interaction vertex. The jet mass is computed from the weighted sum of the constituent four-momenta.
The missing transverse momentum vector ( p miss T ) is computed as the negative vector sum of the transverse momenta of all the PF candidates in an event, and its magnitude is denoted as p miss T [47]. The p miss T is modified to account for corrections to the energy scale of the reconstructed jets in the event.

Signal and background processes
Monte Carlo simulations of the signal process X → YH → bbbb, with a width of 1 MeV for all the three scalars, are generated at leading order (LO) using the MADGRAPH5 aMC@NLO2.6.5 [48] event generator. The NMSSM model [49,50] is used to produce the simulated samples. However, the kinematic parameters are model-independent, enabling the results to be interpreted using other BSM scenarios.
The two main backgrounds are tt+jets events, where the top quarks decay hadronically, and events with jets arising purely from SM quantum chromodynamics (QCD) interactions (multijet events). Other sources of background like single top quark production, and Higgs boson production in association with a top quark pair or a W or Z boson are found to have negligible contributions.
The tt+jets events with hadronic top quark decays are modelled using POWHEG2.0 [51][52][53][54], at next-to-leading order (NLO). A sample of semileptonic tt decays, with one of the top quarks decaying via t → Wb → νb, being a lepton (electron or muon), is also simulated using the same configuration. These events are used in dedicated tt enriched control regions to derive additional data-to-simulation correction factors. The simulated tt+jets event yields are scaled using a cross section of 832 +46 −52 pb, calculated at next-to-next-to-leading order (NNLO) in QCD with soft gluon resummation at next-to-next-to-leading logarithmic precision [55]. The QCD multijet event samples, containing two to four jets, are simulated at LO using the MAD-GRAPH5 aMC@NLO event generator and are used to develop the tools for the analysis. However, this background is estimated using data-driven techniques.
The signal and semileptonic tt+jets samples are generated using the NNPDF3.1 [56] NNLO parton distribution functions (PDFs) from the LHAPDF6 PDF library [57]. The hadronic tt+jets samples are generated using NNPDF3. The showering and hadronization of partons are simulated with PYTHIA8.226 [59] for 2016 and PYTHIA8.240 for 2017 and 2018 samples. The jet-to-parton matching for all LO samples, i.e. the signal and the multijets background, use the MLM [60] scheme. The CP5 tune [61] is used for all samples, except for the 2016 tt and multijet samples, which use CUETP8M2T4 [62] and CUETP8M1 [63] tunes, respectively.
All generated events are processed through a simulation of the CMS detector based on GEANT4 [64]. The effects of pileup are modelled assuming a total inelastic pp cross section of 69.2 mb [65]. All simulated event samples are weighted to match the distribution of the expected pileup profile of the data.

Event selection
The events are selected in two mutually-exclusive categories: an "all-jets" event sample containing only jets, and a "jets+lepton" sample, containing a lepton (electron or muon). The latter serves to derive corrections to the simulated tt+jets background, in order to match the expectations in the data.

All-jets event selection
A set of triggers based on requirements on jet properties are used for online event selection in the all-jets category.
One trigger criterion required a single AK8 jet with p T > 450 or 500 GeV in 2016 and in 2017-2018, respectively. A second trigger required that the scalar sum (H T ) of the p T of all AK4 jets with p T > 30 GeV and |η| < 2.5 should be greater than 800 or 900 GeV in 2016, depending on the LHC beam instantaneous luminosity. In 2017-2018, H T > 1050 GeV was required.
The third trigger algorithm used required an AK8 jet with a trimmed mass >30 GeV along with p T > 360 GeV (in 2016). In 2017-2018, the AK8 jet p T threshold in this trigger was raised to 400 or 420 GeV, depending on the LHC beam instantaneous luminosity, keeping the same trimmed mass criterion. The fourth trigger required H T > 650 or 700 GeV (in 2016) and H T > 800 GeV (in 2017-2018), together with an AK8 jet having a trimmed mass >50 GeV.
In addition to the above, three trigger algorithms were used in 2016 only, with the following criteria: (1) two AK8 jets with p T > 280 and >200 GeV with one of them having a trimmed mass >30 GeV; (2) having the same requirements as (1) and with one of the AK8 jets passing a loose b tagging criterion using the "combined secondary vertex" algorithm [66] (with efficiency ≈ 81%); (3) H T ≥ 650 GeV with a pair of AK4 jets having an invariant mass >900 GeV with their pseudorapidity separation |∆η| < 1.5.
The combined logical OR of all the triggers improves the overall trigger efficiency, particularly for signals with low values of M X .
Events in the offline all-jets selection are required to have at least two AK8 jets with p T > 350 (450) GeV and |η| < 2.4 (2.5) for 2016 (2017-2018). The higher p T requirement in 2017-2018 reflects the higher trigger thresholds and ensures a trigger efficiency close to 100%. The AK8 jet pairs in multijet backgrounds tend to have a larger separation in pseudorapidity than the signal, for a given M JJ range, and therefore a selection |∆η| < 1.3 is used to reduce such backgrounds.
The two leading-p T jets are considered for H → bb and Y → bb candidates. An H → bb 4.1 All-jets event selection candidate or an "H jet" is a jet whose soft-drop mass is 110 < M H J < 140 GeV. The second jet is designated as the Y → bb candidate or the "Y jet" if its soft-drop mass satisfies M Y J > 60 GeV. When both AK8 jets satisfy the first mass requirement, the Y jet is chosen at random. Events without either an H or a Y jet are rejected. The mass of the Y jet and the invariant mass of the H and Y jets are used to isolate the signal with approximately 15% and 9% resolution in M Y J and M JJ , respectively.
The all-jets event category trigger efficiency is measured in the data requiring a single AK4 jet with p T > 260 GeV by applying the above offline selection, and counting the number of events passing the trigger selection. It is found to be between 94 and 100%. Simulated events are weighted by this efficiency as a function of the invariant mass of the two leading-p T AK8 jets in the event, M JJ .
A graph convolutional neural network algorithm, ParticleNet [67], is employed to discriminate the decays of a boosted massive particle R, which could be an H boson or a Y resonance, to a pair of b quarks against a background of other jets, using the properties of the jet PF constituents as features. As with all heavy-flavour jet classifiers, displaced tracks and vertices are the most important features. The multiclassifier ParticleNet algorithm outputs several variables, each in the range 0-1, and each of which can be interpreted as the probability of a jet having originated from a certain decay, such as from a massive resonance R → bb (P(R → bb )) or from a light-flavoured quark or a gluon (P(QCD)). In this analysis, the ParticleNet score is defined as P(R → bb )/(P(R → bb ) + P(QCD)), where P(R → bb ) is a unified score for jets originating from H or Y decays.
The ParticleNet algorithm is trained [68] on AK8 jets using as the signal simulated Lorentzboosted spin-0 particles decaying to a pair of b quarks, with a wide range of masses. The QCD multijet samples are used for the background. The wide signal mass range in the training sample ensures that the background rejection rate is decorrelated from the mass of the jet [68]. As a consequence, background enriched regions can be defined using low ParticleNet scores on jets that have the same mass spectra as that of the background in the signal region. An accurate background model can therefore be developed.
The ParticleNet scores used for selecting the H → bb and the Y → bb candidates ("signal jets") are either >0.98 (tight requirement) or >0.94 (loose requirement). Depending on the jet p T , the former has an efficiency of 62-72% and a misidentification rate of 0.45%, while the latter has an efficiency of 80-85% and a misidentification rate of 1%.
The efficiency of the ParticleNet classifier is calibrated in data using a sample of jets originating from fragmentation of a gluon to bb (g → bb), which are similar to H → bb and Y → bb jets. Such jets are selected from the data using a boosted decision tree (BDT) classifier, such that their ensemble ParticleNet score resembles that of Y and H jets. Using simulated multijet events, the BDT is trained to separate g → bb jets from jets of other flavours. A systematic uncertainty is assigned to account for different possible choices of such jets. The measurements give a data-to-simulation correction factor of 0.9-1.4 for the ParticleNet selection efficiencies, depending on the jet p T and data-taking year.
The ParticleNet scores of the H and Y jets are used to classify events into either signal, sideband, or validation categories. A layout of the different regions is shown in Fig. 1.
Two signal regions are defined using the tight and the loose ParticleNet scores (Fig. 1). The "signal region 1" (SR1) and the "signal region 2" (SR2) are statistically exclusive. SR1 has a higher signal-to-background ratio and is thus more sensitive to the presence of signal. However, the SR2 improves the sensitivity for signal mass points with low background by increasing the signal efficiency.
Corresponding to the two signal regions, two "sideband regions" are defined for estimating the multijet background from data. They are labelled as "Sideband 1" (SB1) and "Sideband 2" (SB2) in Fig. 1. The SB2 region includes SB1 in order to provide better sideband region characteristics for estimating the multijets background in their respective signal regions.
In addition, six "validation regions" are used to validate the background estimation method. They are grouped into two sets of three regions: VS1, VS2, VB1, and VS3, VS4, VB2 as shown in Fig. 1. The labels VS and VB stand for "signal-like" and "background-like" validation regions, respectively. All these regions are enriched in QCD multijet events, and have much smaller signal-to-background ratios than the signal regions.
The signal selection efficiencies range from 1.7% to 12.6% in SR1 and 1.3% to 5.6% in SR2. Based on simulation, the background composition is about an equal proportion of tt+jets and QCD multijets in the signal regions and in the corresponding validation regions, VS1-VS4. However, the sideband and validation sideband regions are composed ≈ 90% of multijet events.

Jets+lepton event selection
The triggers in the jets+lepton category required events to have either an isolated muon of p T > 24 or 27 GeV; an isolated electron having p T > 27, 32, or 35 GeV, or a photon with p T > 175 or 200 GeV. The thresholds changed between data-taking years. The jets+lepton event trigger efficiencies are measured in a sample of Z → events and are found to be close to 100%.
The loose DeepJet working point, with a mistag rate of 10% and approximately 90% efficiency, is used. The DeepJet score distributions of the AK4 jets are corrected using a weight extracted from measurements in the data [66]. Requirements of p miss T > 60 GeV and H T > 500 GeV are imposed. The lepton, p miss T , and the b-tagged jet provide the signature of the leptonic decay of a top quark. A hadronically decaying top quark candidate is reconstructed from an AK8 jet with p T > 350 (450) GeV and |η| < 2.4 (2.5) for 2016 (2017-2018), a soft-drop mass >60 GeV, and satisfying ∆R(lepton, AK8 jet) > 2. Events in the jets+lepton category, which has a purity of >90%, are split into two regions based on whether the AK8 jet passes the tight or loose ParticleNet scores. Two separate correction factors are derived, one each for SR1 and SR2.

Background estimation
The analysis searches for a narrow signal in the 2-dimensional plane spanned by M JJ and M An initial estimate R init P/F is made using the first set of validation regions, using the data and correcting for the simulated tt+jets component: VS1-to-VB1 and VS2-to-VB1. With the definition R P/F ≡ R init P/F R ratio , only the correction function R ratio needs to be determined directly from the fit to the data. The validation regions provide a good estimate of R P/F , because the pass-to-fail event ratios SR1-to-SB1 and SR2-to-SB2 are close to VS1-to-VB1 and VS2-to-VB1, as borne out in simulations. The values of R ratio are therefore of order unity and lead to stability of the fit of signal and background models to the data.
The values of R init P/F , closely related to the loose and tight misidentification rates of the Parti-cleNet tagger, are determined as functions of M Y J only. The 1-dimensional modelling reduces the statistical uncertainties in the R init P/F . A quadratic function is found to be the best model. Furthermore, the R P/F dependence on M JJ is weaker and is modelled through the R ratio , determined directly from the fit to the data in the signal regions.
The form of the R ratio is chosen to be a product of two polynomials in M The simulated tt+jets event distributions in (M JJ , M Y J ) for the signal regions are corrected for their shape and yield using the jets+lepton event category, which is highly enriched in this background. The AK8 jets from top quark decays fall into three categories, depending on the top quark boost. A high enough boost may result in a fully merged t → Wb → qq b decay, labelled as a bqq jet. At moderate boosts, the W → qq may be merged to form a W jet with the b quark forming its own jet. However, such events are nearly all eliminated in the event selection. Finally, one of the quarks from the W boson decay can merge with the b quark to form a bq jet. Unmerged jets and other combinatorial backgrounds constitute a small fraction outside these three categories.
The masses of the bqq and bq jet components in the jets+lepton event category are fit to the data simultaneously with the all-jets event categories. These two mass distributions are scaled independently, with each being tied to the corresponding jet component from the tt+jets in SR1 and SR2. They are independent for the three years, giving six scales in total ranging from 0.79 to 1.35. Two sets of cross-checks are performed for the background estimation method. The first is to predict the background in the validation regions VS1 and VS2 using the validation region VB1 as a sideband. The R init P/F are estimated from the ratios of events in the regions VS3 to VB2 and VS4 to VB2. The jets+lepton regions are treated as they would be for the true background estimation in the signal regions. Similar to the actual background estimation process, a Fisher's test is used to decide the polynomial order of the R ratio function. Again, the most favoured form for R ratio is found to be the product of linear functions along both M JJ and in M Y J . A goodnessof-fit test confirmed the agreement between the data and the estimated background, with the p-value greater than 0.05.
The second check uses generated toy data sets for SR1 and SR2. A toy QCD multijets background is first obtained for these regions by applying the R P/F of VS1 and VS2, obtained in the first validation exercise, to SB1 and SB2, respectively. The toy multijets background is then combined with the tt+jets sample and different signal strengths to get the toy data sets. The test consists of comparing the estimated and true signal strengths after the full background estimation and signal extraction procedure. The test shows no bias in the estimated signal yields for a wide range of M X and M Y . • ParticleNet scale factor: the uncertainty is 7-37%, depending on the AK8 jet p T , and affects the signal by 15%.

Systematic uncertainties
• Jet energy scale and resolution: the uncertainties are applied to both AK4 and AK8 jets, and are fully correlated between the two sets of jets. The signal is affected by 5%.
• Jet mass scale: this is modelled as a ±5% shift in the AK8 jet soft-drop mass. It is uncorrelated between the bqq, the bq, and the signal jets. It affects the signal by 13%. The jet mass scale uncertainty in the tt+jets background is reduced by including the jets+lepton control region.
• Jet mass resolution: simulated AK8 jet masses are smeared to match their distributions in the data, based on studies using Lorentz-boosted W → qq (W boson jets). The nominal simulated unsmeared jet mass resolution is taken as the downward uncertainty while applying a 20% larger smear [41] is used to estimate the upward uncertainty in the AK8 jet mass resolution. The resultant impact on the signal yield is an uncertainty of 4%.
The following uncertainties affect only the backgrounds.
• tt normalization: the uncertainties in the bqq and bq jet scale factors range from 6% to 16%.
• Top quark p T modelling in Monte Carlo simulations: an uncertainty is assigned to the tt+jets simulation process [73], resulting in a 2% uncertainty in this background.
• Multijets background uncertainty: the uncertainty derives mainly from the uncertainty in the R init P/F (M Y J ), which is driven by the sample sizes in the sideband regions VS1, VS2, and VB1. It corresponds to a 7-11% change in the multijet background yields.
Other systematic uncertainties with minor impact are the following.
• Trigger efficiency: the difference between the jet energy scale at the HLT and in the offline reconstruction [74] is appreciable for M JJ < 1100 GeV, resulting in the trigger efficiency dropping below 100%. An uncertainty equal to half of the difference between unity and the measured trigger efficiency is assigned. It is larger than the statistical uncertainty and is expected to cover jet energy scale uncertainties in the trigger efficiencies. Its maximum value is 3%.
• Trigger timing correction: during the 2016-2017 data taking, a gradual shift in the timing of the inputs of the ECAL hardware level trigger in the region of |η| > 2 caused a specific trigger inefficiency. To take this effect into account, a 2% normalization uncertainty is applied to tt events and signal for these years.
• Pileup: the value of the pp total inelastic cross section that is used in the simulation of pileup events is varied upwards and downwards from its assumed value of 69.2 mb by its uncertainty of 4.6% [65].
• PDF and scale uncertainties: the impact of the PDF and the QCD factorization µ F and renormalization µ R scale uncertainties in the signal acceptance and selection is estimated to be 1%. The former is derived using the PDF4LHC procedure [75] and the NNPDF3.1 PDF sets. The latter is evaluated by separately changing µ R and µ F in simulation by factors of 0.5 and 2.
• Sample size of sideband regions: the effects of the limited sizes of the SB1 and SB2 samples are included as statistical uncertainties in the multijets background predicted in SR1 and SR2, using the Barlow-Beeston Lite prescription [76,77]. These uncertainties are small compared to the uncertainties in R init P/F . • Lepton ID and isolation efficiencies: the data-to-simulation correction factors for the efficiencies have uncertainties that affect the event yields by 1-2% in the jets+lepton selection.
• AK4 jet b-tagging data-to-simulation scale factor uncertainty: this uncertainty amounts to about 2% and affects the semileptonic tt+jets event yields.
All uncertainties affecting the signal and the tt+jets samples are uncorrelated among years, except those associated with the PDF choice, the renormalization and factorization scales, the pileup correction, the integrated luminosity, and the top quark p T modelling.

Results
The   The signal hypothesis with M X = 1.6 TeV and M Y = 90 GeV gives the highest observed local significance of 3.1σ, which becomes 0.7σ after accounting for the look-elsewhere effect [78]. However, the excess is not apparent in Fig. 2  The upper limits are computed with a modified frequentist approach, using the CL s criterion [79,80] with the profile likelihood ratio used as the test-statistic and with the asymptotic approximation [81]. As the signal distributions only assume that they originate from spin-0 particle decays, the limits are model-independent. The expected and observed limits at 95% CL as a function of M X and M Y are shown in Fig. 4, and range from 0.1 fb to 150 fb.
The cross section limits are compared with the maximally allowed cross sections in the NMSSM and TRSM. In the NMSSM, no mass range is excluded by the median expected limits. However, the observed limits exclude an area with M X range of 1.00-1.15 TeV and M Y range of 101-145 GeV. For TRSM, an expected exclusion area with the bounds 0.90 < M X < 1.26 TeV and 100 < M Y < 126 GeV is found while the observed exclusion range spans 0.95 < M X < 1.33 TeV and 110 < M Y < 132 GeV.

Summary
A search for massive scalar resonances X and Y, where X decays to Y and the standard model Higgs boson H, has been performed using proton-proton collision data collected at the LHC by the CMS detector between 2016 and 2018, and corresponding to an integrated luminosity of 138 fb −1 . Events are selected using jet substructure and neural network based boosted H/Y → bb identification algorithms. Upper limits at 95% confidence level are set on the cross section of the process pp → X → YH → bbbb for assumed masses of X in the range 0.9-4 TeV and Y between 60-600 GeV. The expected and observed cross section limits for the considered process, set between 0.1 and 150 fb, are the most stringent to date over much of the explored mass range. These limits are interpreted as exclusion of possible M X and M Y within the frameworks of the next-to-minimal supersymmetric model and the two-real-scalar-singlet extension of the standard model.

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 centres and personnel of the Worldwide LHC Computing Grid and other centres 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, the CMS detector, and the supporting computing infrastructure provided by the following funding agencies: