Search for heavy resonances decaying to a $W$ or $Z$ boson and a Higgs boson in the $q\bar{q}^{(\prime)}b\bar{b}$ final state in $pp$ collisions at $\sqrt{s} = 13$ TeV with the ATLAS detector

A search for heavy resonances decaying to a $W$ or $Z$ boson and a Higgs boson in the $q\bar{q}^{(\prime)}b\bar{b}$ final state is described. The search uses 36.1 fb$^{-1}$ of proton-proton collision data at $\sqrt{s} =$ 13 TeV collected by the ATLAS detector at the CERN Large Hadron Collider in 2015 and 2016. The data are in agreement with the Standard Model expectations, with the largest excess found at a resonance mass of 3.0 TeV with a local (global) significance of 3.3 (2.1) $\sigma$. The results are presented in terms of constraints on a simplified model with a heavy vector triplet. Upper limits are set on the production cross-section times branching ratio for resonances decaying to a $W$ ($Z$) boson and a Higgs boson, itself decaying to $b\bar{b}$, in the mass range between 1.1 and 3.8 TeV; the limits range between 83 and 1.6 fb (77 and 1.1 fb) at 95% confidence level.

Upper limits are set on the production cross-section times branching ratio for resonances decaying to a W (Z) boson and a Higgs boson, itself decaying to bb, in the mass range between 1.1 and 3.8 TeV; the limits range between 83 and 1.6 fb (77 and 1.1 fb) at 95% confidence level. c 2017 CERN for the benefit of the ATLAS Collaboration. Reproduction of this article or parts of it is allowed as specified in the CC-BY-4.0 license.

Introduction
The discovery of the Higgs boson [1,2] confirms the validity of the Standard Model (SM) in the description of particle interactions at energies up to a few hundred GeV. However, radiative corrections to the Higgs boson mass drive its value to the model's validity limit, indicating either extreme fine-tuning or the presence of new physics at an energy scale not far above the Higgs boson mass. It is natural to expect such new physics to manifest itself through significant coupling to the Higgs boson, for example in decays of new particles to a Higgs boson and other SM particles. This Letter presents a search for resonances produced in 36.1 fb −1 of proton-proton (pp) collision data at √ s = 13 TeV which decay to a W or Z boson and a Higgs boson. Such resonances are predicted in multiple models of physics beyond the SM, e.g. composite Higgs [3,4] or Little Higgs [5] models, or models with extra dimensions [6,7].
This search is conducted in the channel where the W or Z and Higgs bosons decay to quarks. The high mass region, with resonance masses m V H > 1 TeV (V = W, Z), where the V and H bosons are highly Lorentz boosted, is considered. The V and H boson candidates are each reconstructed in a single jet, using jet substructure techniques and b-tagging to suppress the dominant background from multijet events and enhance the sensitivity to the dominant H → bb decay mode. The reconstructed dijet mass distribution is used to search for a signal and, in its absence, to set bounds on the production cross-section times branching ratio for new bosons which decay to a W or Z boson and a Higgs boson.
The results are expressed as limits in a simplified model which incorporates a heavy vector triplet (HVT) [8,9] of bosons and allows the results to be interpreted in a large class of models. The new heavy vector bosons couple to the Higgs boson and SM gauge bosons with coupling strength c H g V and to the SM fermions with coupling strength (g 2 /g V )c F , where g is the SM SU(2) L coupling constant. The parameter g V characterizes the interactions of the new vector bosons, while the dimensionless coefficients c H and c F parameterize departures of this typical strength for interactions with the Higgs and SM gauge bosons and with fermions, respectively, and are expected to be of order unity in most models. Two benchmark models are used: in the first, referred to as Model A, the branching ratios of the new heavy vector boson to known fermions and gauge bosons are comparable, as in some extensions of the SM gauge group [10]. In Model B, fermionic couplings to the new heavy vector boson are suppressed, as for example in a composite Higgs model [11]. The regions of HVT parameter space studied correspond to the production of resonances with an intrinsic width that is narrow relative to the experimental resolution. The latter is roughly 8% of the resonance mass. The sensitivity of the analysis to wider resonances is not tested.
Searches for V H resonances, V , have recently been performed by the ATLAS and CMS collaborations. The ATLAS searches (using leptonic V decays) based on data collected at √ s = 8 TeV set a lower limit at the 95% confidence level (CL) on the W (Z ) mass at 1.47 (1.36) TeV in HVT benchmark Model A with g V = 1 [12]. Using the same decay modes, a search based on 3.2 fb −1 of data collected at √ s = 13 TeV set a 95% CL lower limit on the W (Z ) mass at 1

ATLAS detector
The ATLAS detector [20] is a general-purpose particle detector used to investigate a broad range of physics processes. It includes inner tracking devices surrounded by a 2.3 m diameter superconducting solenoid, electromagnetic and hadronic calorimeters and a muon spectrometer with a toroidal magnetic field. The inner detector consists of a high-granularity silicon pixel detector, including the insertable B-layer [21] installed after Run 1 of the LHC, a silicon strip detector, and a straw-tube tracker. It is immersed in a 2 T axial magnetic field and provides precision tracking of charged particles with pseudorapidity |η| < 2.5. 1 The calorimeter system consists of finely segmented sampling calorimeters using lead/liquid-argon for the detection of electromagnetic (EM) showers up to |η| < 3.2, and copper or tungsten/liquid-argon for electromagnetic and hadronic showers for 1.5 < |η| < 4.9. In the central region (|η| < 1.7), a steel/scintillator hadronic calorimeter is used. Outside the calorimeters, the muon system incorporates multiple layers of trigger and tracking chambers within a magnetic field produced by a system of superconducting toroids, enabling an independent precise measurement of muon track momenta for |η| < 2.7. A dedicated trigger system is used to select events [22]. The first-level trigger is implemented in hardware and uses the calorimeter and muon detectors to reduce the accepted rate to 100 kHz. This is followed by a software-based high-level trigger, which reduces the accepted event rate to 1 kHz on average.

Data and simulation samples
This analysis uses 36.1 fb −1 of LHC pp collisions at √ s = 13 TeV collected in 2015 and 2016. The data were collected during stable beam conditions with all relevant detector systems functional. Events were selected using a trigger that requires a single anti-k t jet [23] with radius parameter R = 1.0 (large-R jet) with a transverse energy (E T ) threshold of 360 (420) GeV in 2015 (2016). The trigger requirement is > 99% efficient for events passing the offline selection of a large-R jet with transverse momentum (p T ) > 450 GeV.
Signal and backgrounds from tt and W/Z + jets production are modelled with Monte Carlo (MC) simulation. While multijet MC events are used as a cross-check, the primary multijet background estimation is performed using data as described in Section 6. The signal is modelled using benchmark Model A with g V = 1. Results derived from this model can be directly applied to benchmark Model B by rescaling the 1 ATLAS uses a right-handed coordinate system with its origin at the nominal interaction point (IP) in the centre of the detector and the z-axis along the beam pipe. The x-axis points from the IP to the centre of the LHC ring, and the y-axis points upward. Cylindrical coordinates (r, φ) are used in the transverse plane, φ being the azimuthal angle around the z-axis. The pseudorapidity is defined in terms of the polar angle θ as η = − ln tan(θ/2). The rapidity is defined relative to the beam axis as y = 1/2 ln((E + p z )/(E − p z )). Track jets clustered using the anti-k t algorithm with R = 0.2 are used to aid the identification of bhadron candidates from the Higgs boson decay [45]. Track jets are built from charged particle tracks with p T > 0.4 GeV and |η| < 2.5 that satisfy a set of hit and impact parameter criteria to minimize the impact of tracks from pile-up interactions, and are required to have track jet p T > 10 GeV, |η| < 2.5, and at least two tracks clustered in the track jet. Track jets are matched with large-R jets using ghost association. The identification of b-hadrons relies on a multivariate tagging algorithm [46] which combines information from several vertexing and impact parameter tagging algorithms applied to a set of tracks in a region of interest around each track jet axis. The b-tagging requirements result in an efficiency of 77% for track jets containing b-hadrons, and a misidentification rate of ∼ 2% (∼ 24%) for light-flavour (charm) jets, as determined in a sample of simulated tt events. For MC samples the tagging efficiencies are corrected to match those measured in data [47].
Muons are reconstructed by combining tracks in the inner detector and the muon system, and are required to satisfy "Tight" muon identification criteria [48]. The four-momentum of the closest muon candidate with p T > 4 GeV and |η| < 2.5 that is within ∆R = (∆η) 2 + (∆φ) 2 = 0.2 of a track jet is added to the calorimeter jet four-momentum to partially account for the energy carried by muons from semileptonic b-hadron decays. This muon correction results in a ∼ 5% resolution improvement for Higgs boson candidate jets (defined in Section 5) [49]. Electrons are reconstructed from inner detector and calorimeter information, and are required to satisfy the "Loose" likelihood selection [50].
Leptons (electrons and muons, ) are also used in a "veto" to ensure the orthogonality of the analysis selection with respect to other heavy V H resonance searches in non-fully hadronic final states. The considered leptons have p T > 7 GeV, |η| < 2.5 (2.47) for muons (electrons), and their associated tracks must have |d 0 |/σ d 0 < 3 (5) and |z 0 sin θ| < 0.5 mm, where d 0 is the transverse impact parameter with respect to the beam line, σ d 0 is the uncertainty on d 0 , and z 0 is the distance between the longitudinal position of the track along the beam line at the point where d 0 is measured and the longitudinal position of the primary vertex. Leptons are also required to satisfy an isolation criterion, whereby the ratio of the p T sum of all tracks with p T > 1 GeV (excluding the lepton's) within a cone around the lepton (with radius dependent on the lepton p T ) to the lepton momentum must be less than a p T -and |η|-dependent threshold I 0 . The value of I 0 is chosen such that a constant efficiency of 99% as a function of p T and |η| is obtained for leptons in events with identified Z → candidates.
The missing transverse momentum ( E miss T ) is calculated as the negative vectorial sum of the transverse momenta of all the muons, electrons, calorimeter jets with R = 0.4, and any inner-detector tracks from the primary vertex not matched to any of these objects [51]. The magnitude of the E miss T is denoted by E miss T .

Event selection
Events selected for this analysis must contain at least two large-R jets with |η| < 2.0 and invariant mass m J > 50 GeV, and cannot have any lepton candidate passing the veto for leptons. The leading and subleading p T large-R jets must have p T greater than 450 GeV and 250 GeV, respectively. The two leading p T large-R jets are assigned to be the Higgs and vector boson candidates, and the invariant mass of the individual jets is used to determine the boson type; the large-R jet with the larger invariant mass is assigned to be the Higgs boson candidate jet (H-jet), while the smaller invariant mass large-R jet is assigned as the vector boson candidate jet (V-jet). In signal MC simulation, this procedure results in 99% correct assignment after the full signal region selections described below. Furthermore, the absolute value of the rapidity difference, |∆y 12 |, between the two leading p T large-R jets must be less than 1.6, exploiting the more central production of the signal compared to the multijet background. To ensure orthogonality with the ZH resonance search in which the Z boson decays to neutrinos, events are rejected if they have E miss T > 150 GeV and ∆φ( E miss T , H-jet) > 120 degrees. The H-jet is further required to satisfy mass and b-tagging criteria consistent with expectations from a Higgs boson decaying to bb [45]. The H-jet mass, m J,H , must satisfy 75 GeV < m J,H < 145 GeV, which is ∼ 90% efficient for Higgs boson jets. The number of ghost associated b-tagged track jets is then used to categorize events. The H-jets with either one or two b-tagged track jets, amongst the two leading p T associated track jets, are used in this analysis. The H-jets with one associated b-tagged track jet are not required to have two associated track jets. The Higgs tagging efficiency, defined with respect to jets that are within ∆R = 1.0 of a truth Higgs boson and its decay b-hadrons, for doubly-(singly-) b-tagged H-jets is ∼ 40% (∼ 75%) for H-jets with p T ≈ 500 GeV and ∼ 25% (∼ 65%) for H-jets with p T ≈ 900 GeV [49]. The rejection factor for jets from multijet production is ∼ 600 (∼ 50) for double (single) tags.
The V-jet must satisfy mass and substructure criteria consistent with a W-or Z-jet using a 50% efficiency working point, similar to the "Medium" working point in Ref. [52]. To be considered a W (Z) candidate, the V-jet must have a mass m J,V within a p T -dependent mass window which varies between m J,V ∈ [67, 95] ([75, 107]) GeV for jets with p T ≈ 250 GeV, and m J,V ∈ [60, 100] ([70, 110]) GeV for jets with p T ≈ 2500 GeV. The jet must also satisfy a p T -dependent D 2 [53, 54] selection (with β = 1) which depends on whether the candidate is a W or a Z boson, as described in Ref. [52]. The variable D 2 exploits two-and three-point energy correlation functions to tag boosted objects with two-prong decay structures. The V-jet tagging efficiency is ∼ 50% and constant in V-jet p T , with a misidentification rate for jets from multijet production of ∼ 2%.
Four signal regions (SRs) are used in this analysis. They differ by the number of b-tagged track jets associated to the H-jet and by whether the V-jet passes a Z-tag or W-tag selection. The "1-tag" and "2tag" SRs require exactly one and two b-tagged track jets associated to the H-jet, respectively. The 2-tag signal regions provide better sensitivity for resonances with masses below ∼ 2.5 TeV. Above 2.5 TeV the 1-tag regions provide higher sensitivity because the Lorentz boost of the Higgs boson is large enough to merge the fragmentation products of both b-quarks into a single track jet. Events in which the V-jet passes a Z-tag constitute the ZH signal regions, while events in which the V-jet passes a W-tag constitute the WH signal regions. While the 1-tag and 2-tag signal regions are orthogonal regardless of the V-jet tag, the WH and ZH selections are not orthogonal within a given b-tag category. The overlap between the WH and ZH selections in the signal regions is approximately 60%.
The final event requirement is that the mass of the candidate resonance built from the sum of the V-jet and H-jet candidate four-momenta, m JJ , must be larger than 1 TeV. This requirement ensures full efficiency for the trigger and jet p T requirements for events passing the full selection. The full event selection can be found in Table 1. The expected selection efficiency for both WH and ZH resonances decaying to qq ( ) (bb + cc) with a mass of 2 (3) TeV in the HVT benchmark Model B is ∼ 30% (∼ 20%).

Background estimation
After the selection of 1-tag and 2-tag events, ∼ 90% of the background in the signal regions originates from multijet events. The remaining ∼ 10% is predominantly tt with a small contribution from V+jets ( 1%). The multijet background is modelled directly from data, while other backgrounds are estimated from MC simulation.
Multijet modelling starts from the same trigger and event selection as described above, but the H-jet is required to have zero associated b-tagged track jets. This 0-tag sample, which consists of multijet events at the ∼ 99% level, is used to model the kinematics of the multijet background in the 1-tag and 2-tag SRs.
To keep the 0-tag region kinematics close to the 1-and 2-tag regions, H-jets in 0-tag events must contain at least one (two) associated track jets when modelling the 1(2)-tag signal region.
The 0-tag sample is normalized to the 1-tag and 2-tag samples and corrected for kinematic differences with respect to the signal regions. These kinematic differences arise from the b-tagging efficiency variations as a function of p T and |η| and because different multijet processes, in terms of quark, gluon, and heavyflavour content, contribute different fractions to the 0-, 1-, and 2-tag samples.
The 0-tag sample is normalized to the 1-and 2-tag samples, separately, using a signal-free high mass sideband of the H-jet defined by 145 GeV < m J,H < 200 GeV. This sideband (SB), illustrated in Figure 1, is orthogonal to the signal region and has similar expected event yield to the signal region. The normalization of the multijet events is set by scaling the number of events in each region of the 0-tag sample by where N Kinematic corrections to the multijet background template are applied by reweighting events from the 0-tag sample. This is performed only for the 2-tag sample, as the modelling of the multijet background in the 1-tag SB and validation regions (described below and depicted in Figure 1) without reweighting is observed to be adequate. The weights are derived in the SB region, from third-order polynomial fits to the ratio of the total background model to data in two distributions that are sensitive to kinematic and b-tagging efficiency differences between the 0-tag and 2-tag samples. The variables are the track jet p T ratio, defined as p lead T /(p lead T + p sublead T ), and p sublead T , both using the p T distributions of the leading two p T track jets associated to the H-jet. The reweighting is performed using one dimensional distributions but is iterated so that correlations between the two variables are taken into account. After each reweighting iteration, the value of µ 1(2)-tag Multijet is recomputed to ensure that the normalization is kept fixed. No explicit uncertainties are associated with this reweighting as these are determined from comparison with validation regions, as described below.
Due to the small number of events in the background template in the high m JJ tail, the backgrounds are modelled by fitting between 1.2 and 4 TeV with power-law and exponential functions. The multijet background in m JJ is modelled using the functional form while the merged tt and V+jets backgrounds are modelled using the functional forms for the 1-tag and 2-tag samples respectively. In these functional forms, x = m JJ / √ s, and p a through p h are parameters determined by the fit. These functional forms are used as they can model changes in the powerlaw behaviour of the respective backgrounds between high and low masses. The exponential function is used for the 2-tag tt and V+jets samples because it was found to model the tail of the distribution well and because a fit to the small statistics of the sample could not constrain a function with more parameters. Fits are performed separately for the 1-tag and 2-tag background estimates, and separately for each background.
The background model is validated in the two regions denoted by VR-SR and VR-SB in Figure 1, each also with two subregions. In all of these, the V-jet is required to have mass 50 GeV < m J,V < 70 GeV but the D 2 selection is only applied in one of the subregions. For the signal-region-like validation regions (VR-SR) the H-jet selection is unchanged, and for the sideband-like validation regions (VR-SB) the H-jet is required to have mass 145 GeV < m J,H < 200 GeV. Both validation regions are kinematically similar to the signal regions but orthogonal to them (and each other). Table 2 compares the observed data yields in the validation regions with the corresponding background estimates. The differences are used as estimators of the background normalization uncertainties, as described in Section 7. The modelling of the m JJ distribution in the signal-region-like validation region is shown in Figure 2 for the 1-tag and 2-tag samples. The data are well described by the background model. Other kinematic variables are generally well described.  Experimental systematic uncertainties affect the signal as well as the tt and V+jets backgrounds estimated from MC simulation. The systematic uncertainties related to the scales of the large-R jet p T , mass and D 2 are of the order of 2%, 5% and 3%, respectively. They are derived following the technique described in Ref. [39]. The impacts of the uncertainties on the resolutions of each of these large-R jet observables are evaluated by smearing the jet observable according to the systematic uncertainties of the resolution measurement [39, 52]. A 2% absolute uncertainty is assigned to the large-R jet p T , and to the mass and D 2 resolutions relative 20% and 15% uncertainties are assigned, respectively. The uncertainty in the btagging efficiency for track jets is based on the uncertainty in the measured tagging efficiency for b-jets in data following the methodology used in Ref. [47]. This is measured as a function of b-jet p T and ranges between 2% and 8% for track jets with p T < 250 GeV. For track jets with p T > 250 GeV the uncertainty in the tagging efficiencies is extrapolated using MC simulation [47] and is approximately 9% for track jets with p T > 400 GeV. A 30% normalization uncertainty is applied to the tt background based on the ATLAS tt differential cross-section measurement [56]. Due to the small contribution of the V+jets background, no corresponding uncertainty is considered.
Systematic uncertainties in the normalization and shape of the data-based multijet background model are assessed from the validation regions. The background normalization predictions in the validation regions agree with the observed data to within ±5% in the 1-tag sample and ±13% in the 2-tag sample. These differences are taken as the uncertainties in the predicted multijet yield. The shape uncertainty is derived by taking the ratio of the predicted background to the observed data after fitting both to a power law. This is done separately for the 1-tag and 2-tag samples. The larger of the observed shape differences in the VR-SR and VR-SB is taken as the shape uncertainty. Separate shape uncertainties are estimated for m JJ above and below 2 TeV in order to allow for independent shape variations in the bulk and tail of the m JJ distribution in the final statistical analysis.
An additional uncertainty in the shape of the multijet background prediction is assigned by fitting a variety of empirical functions designed to model power-law behaviour to the 0-tag m JJ distribution, as described in Ref. [57]. The largest difference between the nominal and alternative fit functions is taken as a systematic uncertainty. Similarly, the fit range of the nominal power-law function is varied, and the largest difference between the nominal and alternative fit ranges is taken as a systematic uncertainty.
The impact of the main systematic uncertainties on event yields is summarized in Table 3.

Results
The results are interpreted using the statistical procedure described in Ref. [1] and references therein. A test statistic based on the profile likelihood ratio [58] is used to test hypothesized values of µ, the global signal strength factor, separately for each model considered. The statistical analysis described below is performed using the m JJ distribution of the data observed in the signal regions. The systematic uncertainties are modelled with Gaussian or log-normal constraint terms (nuisance parameters) in the definition of the likelihood function. The data distributions from the 1-tag and 2-tag signal regions are used in the fit simultaneously, treating systematic uncertainties on the luminosity, jet energy scale, jet energy resolution, jet mass resolution and b-tagging as fully correlated between the two signal regions. Both the multijet normalization and shape uncertainties are treated as independent between the two signal regions. In addition, the multijet shape uncertainties for m JJ above and below 2 TeV are treated as independent. When performing the fit, the nuisance parameters are allowed to vary within their constraints to maximize the likelihood. As a result of the fit, the multijet shape uncertainties are significantly reduced. With the jet mass resolution, jet energy scale and multijet normalization, they have the largest impact on the search sensitivity. Fits in the WH and ZH signal regions are performed separately. The pre-and post-fit m JJ distributions in the signal regions are shown in Figure 3. The number of background events in the 1-tag and 2-tag ZH and WH signal regions after the fit, the number of events observed in the data, and the predicted yield for a potential signal are reported in Table 4. The total data and background yields in each region are constrained to agree by the fit. There is a ∼ 60% overlap of data between the WH and ZH selections for both the 2-tag and 1-tag signal regions, and this fraction is approximately constant as a function of m JJ .

Statistical analysis
To determine if there are any statistically significant local excesses in the data, a test of the backgroundonly hypothesis (µ = 0) is performed at each signal mass point. The significance of an excess is quantified using the local p 0 value, the probability that the background could produce a fluctuation greater than or equal to the excess observed in data. A global p 0 is also calculated for the most significant discrepancy, using background-only pseudo-experiments to derive a correction for the look-elsewhere effect across the mass range tested [59]. The largest deviation from the background-only hypothesis is in the ZH signal region, occurring at m JJ ≈ 3.0 TeV with a local significance of 3.3 σ. The global significance of this excess is 2.1 σ. The data are used to set upper limits on the cross-sections for the different benchmark signal processes. Exclusion limits are computed using the CL s method [60], with a value of µ regarded as excluded at the 95% CL when CL s is less than 5%.  Figure 5 shows the 95% CL limits in the g 2 c F /g V vs. g V c H plane for several resonance masses for both the WH and ZH channels. These limits are derived by rescaling the signal cross-sections to the values predicted for each point in the (g 2 c F /g V , g V c H ) plane and comparing with the observed cross-section upper limit. As the resonance width is not altered in this rescaling, areas for which the resonance width Γ/m > 5% are shown in grey. These may not be well described by the narrow width approximation assumed in the rescaling.     Figure 5: Limits in the g 2 c F /g V vs. g V c H plane for several resonance masses for the (left) ZH and (right) WH channels. Areas outside the curves are excluded. The benchmark model points are also shown. Coupling values for which the resonance width Γ/m > 5% are shown in grey, as these regions may not be well described by the narrow width approximation.

Summary
A search for resonances decaying to a W or Z boson and a Higgs boson has been carried out in the qq ( ) bb channel with 36.1 fb −1 of pp collision data collected by ATLAS during the 2015 and 2016 runs of the LHC at √ s = 13 TeV. Both the vector boson and Higgs boson candidates are reconstructed using large radius jets, and jet mass and substructure observables are used to tag W, Z and Higgs boson candidates and suppress the dominant multijet background. In addition, small radius b-tagged track jets ghost-associated to the large-R jets are exploited to select the Higgs boson candidate jet. The data are in agreement with the Standard Model expectations, with the largest excess observed at m JJ ≈ 3.0 TeV in the ZH channel with a local significance of 3.3 σ. The global significance of this excess is 2.1 σ. Upper limits on the production cross-section times the Higgs boson branching ratio to the bb final state are set for resonance masses in the range between 1.1 and 3.8 TeV with values ranging from 83 fb to 1.6 fb and 77 fb to 1