Inclusive search for highly boosted Higgs bosons decaying to bottom quark-antiquark pairs in proton-proton collisions at $\sqrt{s} =$ 13 TeV

A search for standard model Higgs bosons (H) produced with transverse momentum ($p_\mathrm{T}$) greater than 450 GeV and decaying to bottom quark-antiquark pairs ($\mathrm{b\bar{b}}$) is performed using proton-proton collision data collected by the CMS experiment at the LHC at $\sqrt{s}=$ 13 TeV. The data sample corresponds to an integrated luminosity of 137 fb$^{-1}$. The search is inclusive in the Higgs boson production mode. Highly Lorentz-boosted Higgs bosons decaying to $\mathrm{b\bar{b}}$ are reconstructed as single large-radius jets, and are identified using jet substructure and a dedicated b tagging technique based on a deep neural network. The method is validated with Z$\to\mathrm{b\bar{b}}$ decays. For a Higgs boson mass of 125 GeV, an excess of events above the background assuming no Higgs boson production is observed with a local significance of 2.5 standard deviations ($\sigma$), while the expectation is 0.7. The corresponding signal strength and local significance with respect to the standard model expectation are $\mu_\mathrm{H} = $ 3.7 $\pm$ 1.2 (stat)$^{+0.6}_{-0.7}$ (syst)$^{+0.8}_{-0.5}$ (theo) and 1.9 $\sigma$. Additionally, an unfolded differential cross section as a function of Higgs boson $p_\mathrm{T}$ for the gluon fusion production mode is presented, assuming the other production modes occur at the expected rates.


Introduction
The observation of a new boson consistent with the standard model (SM) Higgs boson (H) and the subsequent measurements of its properties [1][2][3] have advanced the understanding of electroweak (EW) symmetry breaking and the origin of the mass of elementary particles [4][5][6][7][8][9][10][11]. The H boson has been observed at the CERN LHC in all of its main expected production modes and several decay modes, including decays to bottom quark-antiquark pairs (bb) when produced in association with a W or Z boson [12,13]. Recently, there has been considerable interest in the measurement of Higgs bosons produced with high transverse momentum, p T , where measurements in the H(bb ) decay channel have better sensitivity than traditional channels because of its large branching fraction, B(H → bb ) = 58.1% [14]. Advances in the identification of largeradius jets [15][16][17][18][19] due to massive color singlet particles with large transverse momentum and decaying to bb pairs have improved the sensitivity of this channel, as demonstrated by the CMS [20,21] and ATLAS [22] Collaborations. The first search for high-p T H(bb ) events by the CMS Collaboration [23] demonstrated the experimental sensitivity of this channel, with an expected significance of 0.7 standard deviations (σ) based on a different theoretical expectation than the latest one used in this paper. Measurements of high-p T H(bb ) events have shown promise in resolving the loop-induced and tree-level contributions to the gluon fusion (ggH) process, provide an alternative approach to study the top quark Yukawa coupling, complementary to associated H production with a top quark-antiquark pair (ttH), and may be sensitive to effects from physics beyond the SM [24][25][26][27][28][29][30][31].
This paper reports the results of an inclusive search for high-p T Higgs bosons decaying to bb pairs in proton-proton (pp) collisions at √ s = 13 TeV. The data set, collected with the CMS detector at the LHC in 2016-2018, corresponds to an integrated luminosity of 137 fb −1 . The search is inclusive in the Higgs boson production mode. The highly Lorentz-boosted H(bb ) candidates are reconstructed as single large-radius jets with the jet mass consistent with that of the observed Higgs boson [19]. The candidate jet is required to have p T > 450 GeV to satisfy restrictive trigger requirements that suppress the large background from jets produced via the strong interaction, referred to as quantum chromodynamics (QCD) multijet events. To further distinguish the H candidates from the background, the jet is required to have a two-prong substructure, as well as displaced tracks and decay vertices consistent with the H(bb ) signal. The events are divided into six adjacent p T categories. The background from QCD multijet production is difficult to model parametrically, and it is therefore estimated in data by inverting the b tagging requirement, which is designed to have reduced correlation with jet mass and p T . The presence of the W and Z boson resonances in the jet mass distribution is used to constrain various systematic uncertainties and to validate the analysis. A separate control region is used to improve the modeling of the tt background. A simultaneous fit to the distributions of the jet mass in all p T categories is performed to determine the normalizations and shapes of the jet mass distributions for the backgrounds and to extract the inclusive H(bb ) signal strength with respect to the SM expectation. The differential cross section for the ggH Higgs boson p T is also extracted under the assumption that H production through other modes occurs at the SM rate.
In contrast with the previous CMS result, the Higgs boson p T spectrum from ggH production is modeled with the HJ-MINLO generator [32][33][34], which includes effects of the finite top quark mass. The predicted cross section is compatible with the latest theoretical calculations [35,36], and is smaller than that used previously [23]. Another major improvement is the development of a b tagging algorithm based on a deep neural network with better H(bb ) signal efficiency. This paper is organized as follows. A brief description of the CMS detector is given in Section 2. Section 3 provides a summary of the various simulated samples used in the analysis. Section 4 describes the event reconstruction and selection criteria used to define the signal and control regions. The background estimation methods are detailed in Section 5. Section 6 lists the sources of systematic uncertainty and their statistical treatment. Section 7 describes the statistical procedure used to derive the results, and reports the results in terms of signal strength modifiers and differential cross sections. Finally, the results are summarized in Section 8.

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 inside its volume. Within the solenoid volume are a silicon pixel and strip tracker, a lead tungstate crystal electromagnetic calorimeter, and a brass and scintillator hadron calorimeter, each composed of a barrel and two endcap sections. Forward calorimeters extend the pseudorapidity (η) coverage provided by the barrel and endcap detectors. Muons are detected in gas-ionization chambers embedded in the steel flux-return yoke outside the solenoid.
Events of interest are selected using a two-tiered trigger system [37]. The first level, composed of custom hardware processors, uses information from the calorimeters and muon detectors to select events at a rate of around 100 kHz within a time interval of less than 4 µs. The second level, known as the high-level trigger, consists of a farm of processors running a version of the full event reconstruction software optimized for fast processing, and reduces the event rate to around 1 kHz before data storage.
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. [38].

Simulated samples
Simulated samples of signal and background events are produced using various Monte Carlo (MC) event generators, with the CMS detector response modeled by GEANT4 [39].
For 2016 running conditions, the QCD multijet and Z+jets processes are modeled at leading order (LO) accuracy using the the MADGRAPH5 aMC@NLO v2.2.2 generator [40]. The W+jets process is modeled at LO accuracy with MADGRAPH5 aMC@NLO v2.3.3. The vector boson (V) samples include decays of the bosons to all flavors of quarks, V(qq ), and include up to 3 (4) extra hard partons at the matrix element level for W+jets (Z+jets). Jets from the matrix element calculation and the parton shower description are matched using the MLM prescription [41]. The tt and single top quark processes are modeled at next-to-LO (NLO) using POWHEG 2.0 [42][43][44][45][46][47]. Diboson processes are modeled at LO accuracy with PYTHIA 8.205 [48].
For 2017 and 2018 running conditions, the same configurations are used, but with newer generator versions. The QCD multijet and V+jets processes are modeled using MAD-GRAPH5 aMC@NLO v2.4.2, and the diboson processes are modeled with PYTHIA 8.226.
For all three years' data collection conditions, the cross sections for the V+jets samples include higher-order QCD and EW corrections, which improve the modeling of high-p T V boson events [49][50][51][52]. The total cross sections for the diboson samples are corrected to next-to-NLO (NNLO) accuracy with the MCFM 7.0 program [53].
The ggH production is simulated using the POWHEG+HJ-MINLO [32,33,43,54] event generator with mass m H = 125 GeV and including finite top mass effects, following the recommendation in Ref. [33]. Additionally, a sample of ggH events is generated with POWHEG and corrected for finite top mass effect using the same procedure as described in Ref. [23], where the NLO to LO ratio of the p T spectrum is approximated by expanding in powers of the inverse square of the top quark mass. The POWHEG generator is used to model Higgs boson production through vector boson fusion (VBF), VH associated production, and ttH channels [54][55][56]. The p T spectrum of the Higgs boson for the VBF production mode is re-weighted to account for next-to-NNLO corrections to the cross section [57,58]. These corrections have a negligible effect on the yield for this process for events with Higgs boson p T > 450 GeV.
For parton showering and hadronization, the POWHEG and MADGRAPH5 aMC@NLO samples are interfaced with PYTHIA 8.205 (8.230) for 2016 (2017 and 2018) running conditions. The PYTHIA parameters for the underlying event description are set to the CUETP8M1 [59] (CP5 [60]) tune, except for the tt sample for 2016, which uses the CUETP8M2T4 tune [61]. For 2016 samples, the parton distribution function set NNPDF3.0 [62] is used, with the accuracy (LO or NLO) corresponding to that used in the matrix element calculations, while for 2017 and 2018 samples, NNPDF3.1 [63] at NNLO accuracy is used for all processes.

Event reconstruction and selection
Event reconstruction is based on a particle-flow algorithm [64], which aims to reconstruct and identify each individual particle with an optimized combination of information from the various elements of the CMS detector. The algorithm identifies each reconstructed particle as an electron, a muon, a photon, or a charged or neutral hadron. The missing transverse momentum vector is defined as the negative vector sum of the transverse momenta of all the particles identified in the event, and its magnitude is referred to as p miss T . The candidate vertex with the largest value of summed physics-object p 2 T is taken to be the primary pp interaction vertex. The physics objects are the jets, clustered using the jet finding algorithm [65] with the tracks assigned to candidate vertices as inputs, and the associated missing transverse momentum, taken as the negative vector sum of the p T of those jets.
Particles are clustered into jets using the anti-k T algorithm with a distance parameter of 0.8 (AK8 jets) or 0.4 (AK4 jets). The larger radius of the AK8 jet better captures the decay products of the high-p T H(bb ) signal. The clustering algorithms are implemented by the FASTJET package [66]. To mitigate the effect from the contributions of simultaneous pp collisions (pileup), the pileup per particle identification algorithm [67,68] assigns a weight to each particle prior to jet clustering based on the likelihood of the particle to originate from the hard scattering vertex. Further corrections are applied to the jet energy as a function of jet η and p T to bring the average measured response of jets to that of jets made directly from the generated particles before simulation of the detector response [69]. These corrections are derived separately for each data collection period. Jet identification criteria are applied to remove spurious jets associated with calorimeter noise as well as those associated with muon and electron candidates that are either misreconstructed or isolated. Specifically, jets are required to have neutral hadron and photon energy fractions less than 90%, nonzero charged hadron energy fractions, muon energy fractions less than 80%, and at least two constituent particles [70]. Additionally, AK8 jets are rejected if a photon with p T > 175 GeV is reconstructed within the jet.
A combination of several event selection criteria is used for the event trigger, all of which impose minimum thresholds on either the AK8 jet p T or the event H T , defined as the scalar p T sum of all jets in the event with |η| < 3.0. For AK8 jets used in the trigger selection, a minimum threshold is also imposed on the trimmed jet mass [71], where remnants of soft radiation are removed before computing the mass, which allows the H T or p T thresholds to be reduced while maintaining manageable trigger rates. The trigger selection efficiency is greater than 95% for events with at least one AK8 jet with |η| < 2.5, mass greater than 47 GeV and p T > 450 (525, 500) GeV for 2016 (2017, 2018) data.
To reduce backgrounds from SM EW processes, events are vetoed if they contain isolated electrons, isolated muons, or hadronically decaying τ leptons with p T > 10, 10, or 18 GeV and |η| < 2.5, 2.4, or 2.3, respectively. For electrons and muons, an isolation variable is calculated as the pileup-corrected p T sum of the charged hadrons and neutral particles surrounding the lepton divided by the lepton p T . For charged particles, only those associated with the primary vertex are considered in the isolation variable. For neutral particles, the pileup correction consists of subtracting the energy deposited in the isolation cone by charged hadrons not associated with the primary vertex, multiplied by a factor of 0.5. This factor corresponds approximately to the ratio of neutral to charged hadron production in pileup interactions [72]. The isolation variable for electrons and muons is required to be less than 15 or 25%, respectively, depending on η [73,74].
For each event, the leading AK8 jet in p T is selected to be the H(bb ) candidate, which is around 60% efficient for the ggH production mode. Alternative H(bb ) candidate jet selection criteria were considered, but were not found to improve the sensitivity. The AK8 jet is required to have |η| < 2.5. To reduce the top quark contamination, events are vetoed if they have p miss T > 140 GeV, or if they contain a b-tagged [20] AK4 jet with p T > 30 GeV located in the opposite hemisphere from the leading AK8 jet (∆φ(AK4, AK8) > π/2). The chosen threshold for the AK4 jet b-tagging algorithm corresponds to a 1% probability to misidentify a jet arising from a light flavor quark or gluon and a 77% probability to correctly identify a jet arising from a b quark in 2017 detector conditions. Approximately 60% of tt events are rejected by this selection.
The soft-drop (SD) algorithm [75] with angular exponent β = 0 and soft radiation fraction z = 0.1 is applied to the Higgs boson jet candidate to remove soft and wide-angle radiation. The parameter β controls the grooming profile as a function of subjet separation; for β = 0, the algorithm is independent of subjet separation, and is equivalent to the modified mass-drop tagger [76]. The resulting SD jet mass, m SD , is strongly reduced for background QCD multijet events, where large jet masses arise from soft gluon radiation. Conversely, the algorithm preserves the mass of jets from heavy boson decays. Corrections to the m SD values from simulation are derived from a comparison of simulated and measured samples in a region enriched with merged W(qq ) decays from tt events [70]. The m SD corrections remove a residual dependence on the jet p T , and match the simulated jet mass scale and resolution to those observed in data.
The resulting m SD distributions are binned from 47 to 201 GeV with a bin width of 7 GeV. The lower bound is sufficiently above the trigger threshold to be insensitive to differences between the online and offline mass calculations, and the bin width corresponds to the m SD resolution near the V resonances. The dimensionless mass scale variable for QCD multijet jets, ρ(m SD , p T ) = 2 ln(m SD /p T ) [76,77], is used to characterize the correlation between the jet b tagging discriminator, jet mass, and jet p T . Its distribution is roughly invariant in different ranges of jet p T . For each p T category, only those m SD bins that satisfy is the upper m SD (p T ) bound and m lo SD (p lo T ) is the lower m SD (p T ) bound. In this restriction, the lower p T bound is weighted more heavily because of the steeply falling QCD multijet p T distribution. The upper bound on ρ is imposed to avoid instabilities at the edges of the distribution due to finite cone limitations from the jet clustering, while the lower bound on ρ avoids the nonperturbative regime of the m SD calculation. This requirement is about 98% efficient for the H(bb ) signal.
The N 1 2 variable [78] is used to determine how consistent a jet is with having a two-prong substructure. It is based on a ratio of 2-point ( 1 e 2 ) and 3-point ( 2 e 3 ) generalized energy correlation functions [79]: where z i represents the energy fraction of the constituent i in the jet, and ∆R ij is the angular separation between constituents i and j. These generalized energy correlation functions v e n are sensitive to correlations of v pairwise angles among n jet constituents [78]. For a two-prong structure, signal jets have a stronger 2-point correlation than a 3-point one. The discriminant variable N 1 2 is defined as The calculation of N 1 2 is based on the jet constituents after application of the SD grooming algorithm to the jet. It provides excellent discrimination between two-prong signal jets and QCD background jets. However, imposing requirements on N 1 2 , or other similar variables, distorts the jet mass distributions differently depending on the jet p T [80]. To minimize this distortion, a transformation is applied to N 1 2 following the designed decorrelated tagger technique [77], reducing its correlation with ρ and p T in multijet events. The transformed variable is defined as is the value corresponding to the 26th percentile of the N 1 2 distribution in simulated QCD events, as a function of ρ and p T . The transformation is derived in bins of ρ and p T . This ensures that the selection N 1,DDT 2 < 0 yields a constant background efficiency for QCD events across the ρ and p T range considered in this search. The chosen efficiency of 26% maximizes the signal sensitivity.
Jets likely to originate from the merging of the fragmentation products of two b quarks are selected using an algorithm based on a deep neural network, composed of multiple layers between input and output, referred to here as the deep double-b tagger (DDBT) [20,21]. The algorithm takes as inputs several high-level observables that characterize the distinct properties of b hadrons and their momentum directions in relation to the two subjet candidate axes, as well as low-level track and vertex observables. Events where the selected AK8 jet is doubleb tagged constitute the "passing," or signal, region, while events failing the DDBT form the "failing" region, which is used to estimate the QCD multijet background in the signal region. Specifically, an AK8 jet is considered double-b tagged if its DDBT discriminator value exceeds a threshold corresponding approximately to a 1% misidentification probability for QCD jets. This threshold corresponds to a 54% efficiency for reconstructed scalar boson resonances with variable masses decaying to bb in the range 40 < m SD < 200 GeV and 450 < p T < 1200 GeV. Compared to the previous double-b tagger (DBT) algorithm [20] used in a prior CMS result [23], the DDBT improves the bb tagging efficiency by a factor of about 1.6 for the same detector conditions and QCD misidentification probability. For SM ggH production specifically, the tagging efficiency is approximately 60%, an improvement over the previous algorithm by a factor of about 1.3. Figure 1 shows the performance curves of misidentification probability for QCD jets versus the identification probability for bb resonance jets for the previous DBT algorithm and the DDBT algorithm in simulation corresponding to the detector conditions in 2017.
After all selections are applied, the Higgs boson candidate jet is categorized into the DDBT passing and failing regions, each with 22 m SD bins evenly dividing the range 47-201 GeV, and  split further into six jet p T categories with bin boundaries of 450, 500, 550, 600, 675, 800, and 1200 GeV. The p T binning is optimized for best signal significance, and the upper m SD bound is due to the requirements imposed on the jet ρ. Specifically, bins that do not satisfy Eq. (1) are removed, resulting in a total of 124 bins each for the passing and failing regions. Namely, the upper m SD bound for the first two p T categories are 166 GeV and 180 GeV, respectively. For the Higgs boson signal processes in the DDBT passing region, the dominant production mode is ggH (56%), followed by VBF (26%), VH (13%), and ttH (5%).

Background estimation
The dominant background in the signal region is QCD multijet production. The V+jets processes are significant resonant backgrounds. The tt process constitutes a significant nonresonant background across the m SD spectrum. Other EW processes, including diboson, triboson, and ttV, are estimated from simulation and found to be negligible.
The V+jets background is modeled using simulation. Their overall contribution is less than 6% of the total background in the DDBT passing region. The normalizations and shapes of the simulated V+jets background are corrected for NLO QCD and EW effects.
The contribution of tt production to the total background is obtained from simulation, where the normalization and DDBT efficiency are corrected with scale factors derived from a ttenriched control sample. The control sample targets semileptonic tt production, consisting of events with an energetic muon with p T > 55 GeV and |η| < 2.1, a leading AK8 jet with p T > 400 GeV, and an additional b-tagged AK4 jet that is separated from the leading AK8 jet by ∆R > 0.8. The AK8 jet with the highest p T is taken to be the candidate jet. Using the same candidate jet requirements that define the signal selection, DDBT passing and failing regions are constructed in both data and simulation. Due to the relatively low event count in the control sample, the inclusive event counts for 47 < m SD < 201 GeV and p T > 400 GeV are used, totalling 438 (6301) events in the data passing (failing) region. The fraction of tt background relative to the total background expected in this control sample is 72%. Both the absolute normalization and DDBT efficiency of the tt contribution are allowed to vary without constraint from the simulation expectation, but are forced to vary identically in the tt control region and the signal region in the simultaneous fit, thus constraining the background expectation and DDBT mistag probability for this process. The net contribution is about 8% of the total background in the 110 < m SD < 131 GeV range of the DDBT passing region.
The main background in the DDBT passing region, QCD multijet production, has a jet mass shape that depends on p T and is difficult to model parametrically. Therefore, we constrain it using the background-enriched failing region, i.e., events failing the DDBT selection, together with a "pass-fail ratio" function, R p/f , representing the different mass distributions in the two regions. The pass-fail ratio is factorized into two components. First, the DDBT discriminator is designed to be uncorrelated from the jet mass: the training procedure incorporates a penalty term to the loss function for differences in the jet mass distribution between the passing and failing events. Nonetheless, the DDBT exhibits some anticorrelation at high tagger discriminator values and low jet mass, i.e., the mass distributions are different in the passing and failing regions. To account for this, the expected pass-fail ratio is taken from simulated QCD multijet events, by fitting a two-dimensional Bernstein polynomial [81] in ρ and p T , QCD (ρ, p T ), to the distributions in simulation. Second, residual differences arise from the discrepancies in tagger performance between data and simulation, which we parametrize using a Bernstein polynomial in ρ and p T . The complete pass-fail ratio in data is given by the product of these two factors, where n ρ is the degree of the polynomial in ρ, n p T is the degree of the polynomial in p T , a k, is a Bernstein coefficient, and is a Bernstein basis polynomial of degree n.
The pass-fail ratio R p/f is determined from a simultaneous binned fit to the m SD data distributions in the DDBT passing and failing regions across the whole jet mass and p T range, accounting for the contributions from signal and non-QCD backgrounds. In this fit, the coefficients a k, (data correction) are fitted with no external constraints and the QCD coefficients are refitted with constraints from the fit to the QCD simulated data. The p T bin widths, which vary from 50 to 400 GeV, are chosen to provide enough data points to constrain the shape of R p/f . To determine the minimum degree of polynomial necessary to fit the data, a Fisher F-test [82] is performed. As the magnitude of data-to-simulation discrepancies can vary among the data samples and their corresponding simulation samples, an F-test is performed independently for each of the three data taking years. For the 2016 data sample, it is found that a polynomial of order (n ρ , n p T ) = (2, 1) is needed to provide a sufficient goodness of fit with respect to increased orders (p > 0.05), while for 2017 and 2018 data, a residual polynomial of order (n ρ , n p T ) = (1, 1) is found to be sufficient.
The 2017 fitted pass-fail ratio R p/f as a function of m SD and p T under the signal-plus-background hypothesis is shown in Fig. 2. In the absence of correlations between m SD , p T , and the DDBT efficiency, the ratio would be approximately 0.01. The majority of the difference from 0.01 is a result of the expected pass-fail ratio, which ranges from 0.007 to 0.018, while the data residual correction ranges from 0.86 to 1.05. The other data taking periods are similar. As discussed in Section 6, the components of the pass-fail ratio are among the largest sources of uncertainty in the analysis.
In order to validate the background estimation method and associated systematic uncertainties, bias studies are performed using an alternative functional form for pass-fail ratio in the background model. Pseudo-experiment data sets are generated assuming the alternative background model, with and without the injection of signal events, and then fit with the nominal signal-plus-background model. No significant bias in the fitted signal strength is observed; specifically, the means of the differences between the fitted and injected signal strengths divided by the fitted uncertainty are found to be less than 15%. Therefore, no additional systematic uncertainty is assigned for this potential bias from the background modeling.

Systematic uncertainties
The systematic uncertainties associated with the jet mass scale, the jet mass resolution, and the N 1,DDT 2 selection efficiency are correlated among the W, Z, and H(bb ) processes. These uncertainties are estimated in data using an independent sample of merged W boson jets in semileptonic tt events, where the hadronically decaying W boson is reconstructed as a single AK8 jet.
For this sample, data events are required to have an energetic muon with p T > 100 GeV and |η| < 2.1 , p miss T > 80 GeV, a high-p T AK8 jet with p T > 200 GeV, and an additional b-tagged AK4 jet separated from the AK8 jet by ∆R > 0.8 with p T > 30 GeV. Using the same N 1,DDT 2 requirement applied in the signal regions, we define two samples, one with events that pass and one with events that fail the N 1,DDT 2 selection, for merged W boson jets in data and simulation. A simultaneous fit to the two samples in m SD is performed in order to extract the selection efficiency of a merged W jet in simulation and in data. The data-to-simulation scale factors for the N 1,DDT 2 selection efficiency are measured separately for the three data taking periods, as listed in Table 1.
The jet mass scale and jet mass resolution data-to-simulation scale factors are extracted from the same fit, and are also shown in Table 1. As the semileptonic tt sample does not contain a large population of very energetic jets, an additional systematic uncertainty is included to account for the extrapolation to very high p T jets. This additional uncertainty is estimated to be 0.5% per 100 GeV, based on a study of fitting the m SD distributions of merged top quark jets in different p T ranges above 350 GeV [83]. In total, the jet mass scale uncertainty increases with jet p T , ranging from 1.2% at 450 GeV to 2.1% at 800 GeV. While the jet mass scale and resolution among the different data collection periods are similar, their data-to-simulation scale factors and uncertainties vary because of the different generator tunes used in the simulations.
The uncertainty on the efficiency of the DDBT is estimated using data and simulation samples enriched in bb pairs from gluon splitting [20]. The gluon splitting samples require that both subjets of an AK8 jet contain a muon, targeting semileptonic decays of the b hadrons. The method is based on yields extracted from fits to the distributions of the jet probability tagger [20,84] discriminant, which uses the signed impact parameter significance of the tracks associated with the jet to obtain a likelihood for the jet to originate from the primary vertex.
Given that the DDBT efficiencies could differ between bb jets from gluon splitting and from color-singlet Z or Higgs boson decays, the efficiencies extracted from the gluon splitting samples are used only to estimate the uncertainty on the DDBT efficiency, and are not used to correct the efficiency. The applied DDBT data-to-simulation scale factor is included in the signal extraction fit as a constrained nuisance parameter, with a nominal value of unity and an uncertainty equal to the difference between the DDBT data-to-simulation scale factor and unity, as shown in Table 1. The scale factor is further constrained via the observed Z boson yield in the passing and failing regions. This strategy differs from that of the previous CMS analysis [23], resulting in an increase in the post-fit systematic uncertainty of the tagger efficiency from 4% to about 14%. The scale factors described above determine the initial distributions of the jet mass for the W(qq ), Z(qq ), and H(bb ) processes. In the fit to data, the jet mass scales and resolutions are treated as constrained nuisance parameters with nominal values and uncertainties as shown in Table 1, and are further constrained by the presence of the V resonances in the jet mass distribution.
The uncertainty associated with the choice of QCD renormalization and factorization scales in the modeling of ggH production is propagated to the overall normalization of the ggH signal via varying each factor by one-half or two around the nominal value and finding the envelope of all combinations of such variations, except those where one scale is multiplied by 0.5 and the other is multiplied by 2 [85,86]. This results in a 30% uncertainty for the POWHEG sample with p T reweighting [23] and as a 20% uncertainty for the HJ-MINLO sample. For the POWHEG sample, the shape of the ggH p T distribution is allowed to vary depending on the Higgs boson p T by up to 30% at 1.2 TeV, without changing the overall normalization. To account for potential p T -dependent deviations due to missing higher-order corrections, uncertainties are applied to the V(qq ) yields that are p T dependent (5-7%) and correlated per p T bin (10%) [49,50,[87][88][89][90][91]]. An additional systematic uncertainty of 2 to 6%, depending on p T , is included to account for potential differences between the higher-order corrections to the W and Z cross sections (EW W/Z decorrelation) [87].
Finally, systematic uncertainties are applied to the W(qq ), Z(qq ), tt, and H(bb ) yields to account for the uncertainties due to the jet energy scale and resolution [92] and the limited simulation sample sizes. Other experimental uncertainties, including those related to the determination of the integrated luminosity [93], variations in the amount of pileup, modeling of the trigger acceptance, and the isolation and identification of leptons are also considered. Table 2 lists the major sources of uncertainty and their observed impact on the Higgs boson signal strength µ H , defined as the ratio of the measured to the SM expected H(bb ) production, in the combined fit. One of the largest sources of statistical uncertainty is the data residual correction to the pass-fail ratio, while the largest source of systematic uncertainty is the expected pass-fail ratio estimated from simulation. Overall, the µ H measurement is limited by statistical sources of uncertainty.

Results
A binned maximum likelihood fit to the observed m SD distributions is performed using the sum of the signal and background contributions. The fit is performed simultaneously in the DDBT passing and failing regions of the six p T categories, as well as in the DDBT passing and failing components of the tt-enriched control region. The fit is performed separately for the three year periods. A combined fit over the three periods is performed for the final result.
The systematics related to theoretical uncertainties are correlated between different years. The chosen test statistic, used to determine how signal or background-like the data are, is based on the profile likelihood ratio [94]. Systematic uncertainties are incorporated into the analysis via nuisance parameters and treated according to the frequentist paradigm. The best-fit value of each signal strength parameter and an approximate 68% confidence level (CL) interval are extracted following the procedure described in Section 3.2 of Ref. [95]. Figure 3 shows the m SD distributions in the combined data set for the DDBT passing and failing regions with the fitted background. The bottom panels of Fig. 3 show the difference between the data and the prediction from the background, divided by the statistical uncertainty in the data. These highlight the contributions from Higgs and V boson production in the failing and passing regions. The W boson contribution in the passing region is due to the misidentification of W (qq) decays by the DDBT. The agreement between the data and the signal-plusbackground model is quantified with a Kolmogorov-Smirnov goodness-of-fit test [96], which yields a p-value of 17%. In Fig. 4, the m SD distributions are shown for each p T category in the passing region. The nuisance parameters related to the jet mass scale uncertainties, whose values ranges up to 2 GeV in the case of the Z boson as discussed in Section 6, do not significantly deviate from their pre-fit expectations.
To validate the substructure and b tagging techniques employed in this search, a maximum likelihood fit is performed using a model where the Z (qq) signal strength (µ Z ) and µ H are left unconstrained. In the DDBT passing region, decays of the Z boson to bb constitute 79% of all Z decays. The measured µ Z value is 1.01 ± 0.05 (stat) +0.20 −0.15 (syst) +0.13 −0.09 (theo). This demonstrates that the Z boson is clearly separable from the background. In this measurement, the dominant source of systematic uncertainty is the DDBT scale factor. For the remainder of results, µ Z is fixed to its expectation, with the corresponding uncertainties, as described in Section 6. Thus, the Z boson resonance is used to further constrain the DDBT scale factor in the Higgs boson measurements.
To extract the Higgs boson signal, three maximum likelihood fits are performed to the data, each with a different degree of reliance on the modeling of the Higgs boson p T spectrum: the nominal inclusive fit using one µ H parameter for all H production modes and all jet p T categories, an alternative fit using an independent µ H parameter for each p T category for all H production modes to assess the compatibility among the p T categories, and a fit which unfolds detector effects to present results for the ggH production mode at the generator level.
In the inclusive fit using the HJ-MINLO sample as the ggH signal model and including the contributions from the other production modes, the measured µ H value is 3.7 ± 1.2 (stat) +0. 6 −0.7 (syst) +0.8 −0.5 (theo). Upper limits at 95% CL using the CL s criterion [97,98] are obtained using asymptotic formulae [99]. The corresponding observed and expected upper limits on µ H at a 95% CL are 3.7 and 2.9, respectively, while the observed and expected significances [100] with respect to the background-only hypothesis are 2.5 σ and 0.7 σ. The measurement exhibits an excess over the SM expectation (µ H = 1), with a significance of 1.9 σ. Table 3 summarizes the measured signal strengths and significances for the Higgs and Z boson pro- The shaded blue band shows the systematic uncertainty in the total background prediction. The bottom panel shows the difference between the data and the total background prediction, divided by the statistical uncertainty in the data. In the failing region, the background model includes a free parameter for each m SD bin, ensuring the nearly perfect agreement between the model and the data. Thus, the statistical uncertainty in the data gives rise to the systematic uncertainty in the background prediction. This is reflected in the fact that the error bar for the data and the uncertainty band for the background are approximately equal in size.  cesses. The primary results using the ggH p T spectrum from HJ-MINLO [32,33] are shown, alongside results using the ggH p T spectrum from Ref. [23] for ease of comparison. The prediction used for the ggH p T spectrum in Ref. [23] is different from that of HJ-MINLO in both shape and total cross section, which is primarily due to different accuracy of finite top mass correction included in the simulation. In particular, the number of ggH signal events predicted by HJ-MINLO in the fiducial region of the analysis is approximately a factor of two smaller than that of Ref. [23], which is reflected in the fitted µ H values and their uncertainties. The fitted signal strength value and its uncertainty are sensitive to the ggH theoretical prediction and associated uncertainty, which are challenging to obtain in the high-p T regime. To assess the compatibility between the observed signal strengths in the different jet p T categories, an alternative fit to the data is performed. In this fit, an independent µ H is assigned to each of the six reconstructed jet p T bins. These signal strengths are unconstrained in the fit and are varied simultaneously. All other parameters are profiled, as in the original fit. Figure 5 (left) illustrates the compatibility in the best fit signal strengths between the different p T categories, showing an excess with respect to the SM expectation for categories with jet p T above 550 GeV. Separately, the same exercise is performed with an independent µ Z in each p T category. The fitted signal strengths, shown in Fig. 5 (right), are consistent with the SM expectation.
To facilitate comparisons with theoretical predictions, we isolate and remove the effects of limited detector acceptance and response to the ggH production cross section using a maximumlikelihood unfolding technique as described in Section 5 of Ref. [24]. In our treatment, the remaining Higgs boson production modes are assumed to occur at SM rates, as predicted by  [101]. As the minimum reconstructed jet p T is 450 GeV, a negligible signal contribution is expected from events with p H T < 300 GeV. The folding matrix M ji , defined as the product of the acceptance and the efficiency for an H(bb ) event in p H T bin j to be found in jet p T bin i, is shown in Fig. 6 for the ggH HJ-MINLO simulation. This matrix is found to be well-conditioned. Therefore, we omit any regularization in the unfolding procedure [102].
The ggH fiducial cross section in each STXS p H T bin is then extracted by scaling the cross section found in simulation, imposing no selection requirements other than those on p H T , by the corresponding signal strength parameter. The uncertainty in this value is taken from the correspondingly scaled signal strength uncertainty. For the theoretical uncertainties, only those that affect the acceptance of signal events into the reconstructed selection are taken into account. Based on the envelope of acceptance values from varying the renormalization and factorization scales by factors of two, this theoretical acceptance uncertainty is estimated to be 2%. We verify that this unfolding procedure is unbiased through signal injection studies.
The result of this unfolding procedure is shown in Fig. 7 and Table 4, along with the predicted cross sections from Ref. [33] and the predictions of the signal event generators described in Section 3. The correlation coefficients among the three p H T bins are shown in Table 5. The measured cross section uncertainty in the first p H T bin is larger because of limited acceptance. The first and second p H T bins have a mild anti-correlation, primarily because of the imperfect jet energy response of the detector, which inflates the corresponding per-bin uncertainties in the unfolded cross section. The observed cross section in the third p H T bin has a smaller relative uncertainty than that in the second bin because of the larger magnitude of the central value in that bin. With respect to the SM, the upward deviation of the cross section in the third p H T bin, when profiling the other two, corresponds to a local significance of 2.6 σ. When considering all three cross section parameters of interest simultaneously, the total deviation from the SM corresponds to a significance of 1.9 σ. Table 4: Measured and predicted ggH differential fiducial cross section as a function of Higgs boson p T . All cross sections are in units of fb. The cumulative cross section predictions from Ref. [33] are converted to differential cross section predictions by subtraction assuming the cumulative cross section uncertainties are fully correlated.

Summary
An inclusive search for the standard model (SM) Higgs boson decaying to a bottom quarkantiquark pair and reconstructed as a single large-radius jet with transverse momentum p T > 450 GeV has been presented. The search uses a data sample of proton-proton collisions at √ s = 13 TeV, corresponding to an integrated luminosity of 137 fb −1 . The associated production of a Z boson and jets is used to validate the method and is measured to be consistent with the SM prediction. The inclusive Higgs boson signal strength is measured to be µ H = 3.7 ± 1.2 (stat) +0. 6 −0.7 (syst) +0.8 −0.5 (theo) = 3.7 +1.6 −1.5 , based on the theoretical prediction from the HJ-MINLO generator for the gluon fusion production mode. The measured µ H corresponds to an observed significance of 2.5 standard deviations (σ) with respect to the background-only hypothesis, while the expected significance of the SM signal is 0.7 σ. The significance of the observed excess with respect to the SM expectation is 1.9 σ. With respect to the previous CMS result, the relative precision of the µ H measurement improves by approximately a factor of two  [33], shown in red, and HJ-MINLO [32], shown in blue. The two predictions are nearly identical. The larger gray band shows the total uncertainty in the measured cross section while the red and blue hatched bands show the uncertainties in the predictions of Ref. [33] and HJ-MINLO, respectively. In the bottom two panels, the dotted line corresponds to a ratio of one. The relative uncertainties in the predictions of Ref. [33] and HJ-MINLO are approximately 10 and 20%, respectively. because of the increased integrated luminosity, an improved b tagging technique based on a deep neural network, and smaller theoretical uncertainties. Finally, the differential cross section for the p T of a Higgs boson produced through gluon fusion, assuming the other production modes occur at the SM rates, in the phase space regions recommended by the LHC simplified template cross section framework has also been presented. An excess is seen for Higgs boson p T > 650 GeV with a local significance of 2.6 σ with respect to the SM expectation including the Higgs boson.

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: