Searches for heavy $ZZ$ and $ZW$ resonances in the $\ell\ell qq$ and $\nu\nu qq$ final states in $pp$ collisions at $\sqrt{s}=13$ TeV with the ATLAS detector

This paper reports searches for heavy resonances decaying into $ZZ$ or $ZW$ using data from proton--proton collisions at a centre-of-mass energy of $\sqrt{s}=13$ TeV. The data, corresponding to an integrated luminosity of 36.1 fb$^{-1}$, were recorded with the ATLAS detector in 2015 and 2016 at the Large Hadron Collider. The searches are performed in final states in which one $Z$ boson decays into either a pair of light charged leptons (electrons and muons) or a pair of neutrinos, and the associated $W$ boson or the other $Z$ boson decays hadronically. No evidence of the production of heavy resonances is observed. Upper bounds on the production cross sections of heavy resonances times their decay branching ratios to $ZZ$ or $ZW$ are derived in the mass range 300--5000 GeV within the context of Standard Model extensions with additional Higgs bosons, a heavy vector triplet or warped extra dimensions. Production through gluon--gluon fusion, Drell--Yan or vector-boson fusion are considered, depending on the assumed model.


Introduction
The discovery of a Higgs boson h with a mass of approximately 125 GeV in 2012 [1, 2] represents a major milestone in the understanding of electroweak symmetry breaking. Subsequent studies [3][4][5][6] have shown that the properties of the new particle are consistent with those of the Standard Model (SM) Higgs boson. Nevertheless, the possibility that the particle is part of an extended Higgs sector or other extension of the SM cannot be ruled out. Many of these models, motivated by hierarchy and naturalness arguments [7-9], predict the existence of new heavy resonances decaying into dibosons. In models with an extended Higgs sector, such as the two-Higgs-doublet models (2HDM) [10] and the electroweak-singlet model [11], a heavy spin-0 neutral Higgs boson (H) can decay into a pair of Z bosons. In extended gauge models [12], a heavier version of the SM W boson (W ) is predicted to decay into ZW, and in models with warped extra dimensions [13,14], spin-2 Kaluza-Klein (KK) excitations of the graviton (G KK ) are expected to decay into ZZ.
This paper reports searches for heavy resonances X decaying into pairs of vector bosons, ZV (V = W, Z). Production through gluon-gluon fusion (ggF), Drell-Yan (DY) and vector-boson fusion (VBF) processes are considered, depending on the assumed model. Representative Feynman diagrams of these processes are shown in Figure 1. Two ZV decay modes are explored: one in which there is a Z boson decaying into a pair of light charged leptons (electrons or muons, denoted by ), Z → , 1 and the other in which a Z boson decays into a pair of neutrinos, Z → νν. In both cases, the vector boson V is required to decay into a pair of quarks, V → qq, leading to X → ZV → qq and X → ZV → ννqq decay modes. Two different reconstruction techniques for the V → qq decay are considered: resolved and merged. The resolved reconstruction attempts to identify two separate small-radius jets (small-R jet, or j) of hadrons from the V → qq decay, while the merged reconstruction uses jet substructure techniques to identify the V → qq decay reconstructed as a large-radius jet (large-R jet or J). When the resonance mass is significantly higher than the V boson mass, the qq pair from the V boson decay can be collimated. In this case, hadrons from the two quarks overlap in the detector and are more efficiently reconstructed as a single large-R jet. The X → ZV → qq searches utilise both reconstruction techniques for the V → qq decay, whereas the X → ZV → ννqq search uses only the merged reconstruction.
Heavy resonances would manifest themselves as resonant structures above the SM background in the invariant-mass distributions of the qq final state and as broad enhancements in the transverse-mass distributions of the ννqq final state. Thus for the ZV → qq decay mode, the invariant masses of the J system (m J ) from the merged reconstruction and of the j j system (m j j ) from the resolved reconstruction are used as the final discriminants for signal-background separation. For the ZV → ννqq decay mode, the final discriminant used is the transverse mass (m T ) of the large-R jet and the missing transverse momentum (E miss T ). Both the ATLAS and CMS collaborations have reported searches for heavy resonances in ZV decays in proton-proton (pp) collisions at √ s = 8 TeV and 13 TeV [15][16][17][18][19]. In addition to qq and ννqq, these searches include fully hadronic (qqqq), semileptonic (qq ν), and fully leptonic ( , ν, νν) final states. This paper extends previous ATLAS searches for ZV resonances in the qq and ννqq final states at √ s = 13 TeV reported in Ref.
[18] and uses a dataset more than ten times larger.

ATLAS detector
The ATLAS detector [20] at the Large Hadron Collider (LHC) [21] covers nearly the entire solid angle 2 around the collision point, and consists of an inner tracking detector surrounded by a thin superconducting solenoidal magnet producing a 2 T magnetic field, electromagnetic and hadronic calorimeters, and a muon spectrometer incorporating three large toroid-magnet assemblies. The inner detector (ID) consists of a high-granularity silicon pixel detector, including an insertable B-layer [22], and a silicon microstrip tracker, together providing precision tracking in the pseudorapidity range |η| < 2.5, complemented by a transition radiation tracker, providing tracking and electron identification information for |η| < 2.0. A lead/liquid-argon (LAr) electromagnetic calorimeter covers the region |η| < 3.2, and hadronic calorimetry is provided by steel/scintillator-tile calorimeters for |η| < 1.7 and by copper/LAr hadronic calorimeters for 1.7 < |η| < 3.2. The forward region is covered by additional LAr calorimeters with copper and tungsten absorbers. The muon spectrometer (MS) consists of precision tracking chambers covering the region |η| < 2.7, and separate trigger chambers covering |η| < 2.4. A two-level trigger system [23] reduces the event rate to approximately 1 kHz for offline investigations.
3 Data, signal models and simulation

Data
The data used in the searches were collected with the ATLAS detector in 2015 and 2016 pp collisions at a centre-of-mass energy of 13 TeV, corresponding to a total integrated luminosity of 36.1 fb −1 .
2 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).
Events used in the X → ZV → qq search were recorded with a combination of multiple singleelectron or single-muon triggers with varying transverse energy E T (electron) and transverse momentum p T (muon), quality, and isolation requirements. The lowest E T or p T requirement without trigger prescaling was 26 GeV for both the electrons and muons. Events for the X → ZV → ννqq search were recorded with an E miss T trigger of varying threshold. The lowest threshold without prescaling is 100 GeV. This trigger is fully efficient for events passing the selection described below.
Events are retained for analysis if they were recorded with all detector systems operating normally and pass data-quality requirements [24]. Collision vertices are formed from tracks with p T > 400 MeV. If an event contains more than one vertex candidate, the one with the highest p 2 T of its associated tracks is selected as the primary vertex. All events are required to contain a primary vertex with at least two associated tracks.

Signal models and simulation
Three classes of models of physics beyond the Standard Model are used as benchmarks for the interpretation of the results. Different resonances are predicted by each of these: a neutral heavy spin-0 Higgs boson H using the narrow-width approximation [25], a spin-1 W boson of the heavy vector triplet (HVT) model [26,27], and a spin-2 KK graviton G KK from the bulk Randall-Sundrum model [13,28,29]. The new HVT bosons couple to the SM Higgs boson and 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 characterises the interaction strength of the new vector bosons, while the dimensionless coefficients c H and c F parameterise departures of this typical strength from interactions with the SM Higgs and gauge bosons and with fermions, respectively, and are expected to be of order unity in most models. In Model A with g V = 1, the branching ratios of the new heavy vector boson to known fermions and gauge bosons are comparable, while in Model B with g V = 3, fermionic couplings to the new heavy vector boson are suppressed, giving rise to larger branching ratios for decays into ZW final states than in Model A. In a third model, VBF Model, the couplings to gauge bosons are similar to those in Model A, but the couplings to fermions are set to 0. In the bulk RS graviton model, the G KK couplings to light fermions are suppressed and decays into final states involving heavy fermions, gauge bosons or Higgs bosons are favoured. The strength of the coupling depends on k/M Pl , where k corresponds to the curvature of the warped extra dimension and M Pl is the effective four-dimensional Planck scale M Pl = 2.4 × 10 18 GeV. The cross section and intrinsic width scale as the square of k/M Pl .
Monte Carlo (MC) samples of H → ZZ were generated by Powheg-Box v1 [30][31][32][33] with the CT10 [34] parton distribution functions (PDF) assuming a Higgs boson with width far smaller than the experimental resolution. Both the ggF and VBF production processes are considered. Benchmark samples of HVT W → ZW and RS graviton G KK → ZZ were generated with MadGraph5_aMC@NLO 2.2.2 [35], using the NNPDF23LO [36] PDF set. For the HVT model, two production modes of the charged vector triplet W , DY and VBF, are considered. The W resonance from DY production of Model A has a width approximately 2.6% of its mass, while for VBF production the width is much narrower, since its couplings to fermions are set to 0 in the VBF signal. For HVT Model B, the resonance widths and experimental signatures are similar to those obtained for Model A. Thus results derived from Model A can be directly applied to benchmark Model B by rescaling the relevant branching ratios. The G KK has a mass-dependent width, 3.7% at 500 GeV and 6.4% at 5000 GeV relative to its mass, for k/M Pl = 1. An RS graviton with k/M Pl = 0.5 is also considered and the samples were obtained by reweighting the generated samples for k/M Pl = 1 to account for the resonance width and cross-section differences. The parton showering and

Event reconstruction
Electrons are identified as isolated energy clusters in the electromagnetic calorimeter matched to ID tracks, with requirements that the transverse energy E T > 7 GeV and pseudorapidity |η| < 2.47. A likelihood-based requirement [60] is imposed to reduce the background from misidentified or non-prompt electrons. Electrons are classified as either 'loose', 'medium' or 'tight' according to the likelihood-based identification criteria described in Ref. [60].
Muons are reconstructed by a combined fit to the ID and MS tracks, and are required to have p T > 7 GeV and |η| < 2.5. Muons must pass identification requirements, based on the numbers of hits in the ID and MS subsystems, and the significance of the difference |q/p MS − q/p ID | [61], where q is the charge and p MS (p ID ) is the momentum of the muon measured in the MS (ID). Similar to electrons, muons are classified as either 'loose', 'medium' or 'tight', following the criteria in Ref. [61].
All electrons and muons are required to be isolated using selections on the sum of track p T (excluding the track associated with the lepton) in a p T -dependent cone around their directions. The isolation selection criteria are designed to maintain a constant efficiency of 99% in the p T -η plane for the reconstructed leptons from Z → decays, and to minimise efficiency loss for highly boosted Z bosons in the relevant kinematic range. Furthermore, leptons are required to have associated tracks satisfying |d 0 /σ d 0 | < 5 (3) and |z 0 × sin θ| < 0.5 mm for electrons (muons), where d 0 is the transverse impact parameter with respect to the beam line, σ d 0 is its uncertainty, 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.
Small-R jets are reconstructed using the anti-k t algorithm [62,63] with a radius parameter of R = 0.4. Energy-and η-dependent correction factors derived from MC simulations are applied to correct jets back to the particle level [64]. Jets are required to have p T > 20 GeV for |η| < 2.5 and p T > 30 GeV for 2.5 < |η| < 4.5. Jets with p T > 30 GeV and |η| < 2.5 are called 'signal' jets. To suppress jets from pile-up interactions, a jet vertex tagger [65] is applied to jets with p T < 60 GeV and |η| < 2.4, based on information about tracks associated with the primary vertex and pile-up vertices.
The 'signal' jets containing b-hadrons are identified using a multivariate algorithm (b-tagging) [66] which is based on information such as track impact-parameter significances and positions of explicitly reconstructed secondary decay vertices. The b-tagging is used for identifying Z → bb decays. The chosen b-tagging algorithm has an efficiency of 70% for b-quark jets in simulated tt events, with a light-flavour jet rejection factor of about 380 and a c-jet rejection of about 12 [67]. Correction factors are applied to the simulated event samples to compensate for differences between data and simulation in the b-tagging efficiency for b-, c-and light-jets. The correction for b-jets is derived from tt events with final states containing two leptons, and the corrections are consistent with unity with uncertainties at the level of a few percent over most of the jet p T range.
Large-R jets are reconstructed with the anti-k t algorithm, but with the radius parameter increased to R = 1.0. To mitigate the effects of pile-up and soft radiation, the large-R jets are trimmed [68]. Trimming takes the original constituents of the jet and reclusters them using the k t algorithm [69] with a smaller radius parameter, R subjet , to produce a collection of subjets. These subjets are discarded if they carry less than a specific fraction ( f cut ) of the original jet p T . The trimming parameters optimised for this search are R subjet = 0.2 and f cut = 5%. The large-R jet four-momenta are recomputed from the selected subjets, and the jet energies are calibrated to particle-level using correction factors derived from MC simulations [70]. The mass of a large-R jet (m J ) is computed using a combination of calorimeter and tracking information [71]. The large-R jets are required to have p T > 200 GeV and |η| < 2.0.
A default ATLAS-wide overlap-removal procedure is applied to the selected leptons and jets to prevent double-counting. For nearby electrons and small-R jets, the jet is removed if the separation between the electron and jet is within ∆R < 0.2; the electron is removed if the separation is within 0.2 < ∆R < 0.4. For nearby muons and small-R jets, the jet is removed if the separation between the muon and jet is within ∆R < 0.2 and if the jet has less than three tracks or the energy and momentum differences between the muon and the jet are small; otherwise the muon is removed if the separation satisfies ∆R < 0.4. To prevent double-counting of energy from an electron inside the large-R jet, the large-R jet is removed if the separation between the electron and the large-R jet is within ∆R < 1.0.
Boson tagging [72, 73] is applied to the large-R jets to select those from the V → qq decays. A p Tdependent requirement is applied to the jet substructure variable D 2 , which is defined as a ratio of twoand three-point energy correlation functions [74,75] that are based on the energies and pairwise angular distances of particles within a jet. This variable is optimised to distinguish between jets originating from a single parton and those coming from the two-body decay of a heavy particle. A detailed description of the optimisation can be found in Refs. [72,73]. The V boson jets are then selected by requiring the large-R jet mass m J to be in a p T -dependent window centred around the expected value of the boson mass from simulations. The boson tagging is designed to provide a constant efficiency (working point) independent of the large-R jet p T for the signals studied. Two such working points, 50% efficiency and 80% efficiency, are used, with corresponding misidentification rates for jets from multijet production of ∼ 2% and ∼ 10%, respectively. For the 50% working point, the width of the W (Z) mass window varies from 23 (28) GeV at p T = 500 GeV to 33 (37)  The missing transverse momentum ( E miss T ) is calculated as the negative vectorial sum of the transverse momenta of calibrated electrons, muons, and small-R jets. Large-R jets are not included in the E miss T calculation to avoid double-counting of energy between the small-R jets and large-R jets. Energy depositions due to the underlying event and other soft radiation are taken into account by constructing a 'soft term' from ID tracks associated with the primary vertex but not with any reconstructed object [76-78]. The track-based missing transverse momentum, p miss T , is the negative vectorial sum of the transverse momenta of all good-quality inner detector tracks that are associated with the primary vertex.

X → ZV → qq search
The X → ZV → qq search explores the VBF and ggF production of a heavy Higgs boson H, the VBF and DY production of an HVT W boson, and the ggF production of a bulk RS graviton G KK . It also utilises both the merged and resolved reconstruction for the V → qq decay. The search begins with the identification of the Z → decay, followed by classifying events into the VBF or ggF/DY categories and finally the selection of either the ZV → J or ZV → j j final states. Multiple signal regions (SRs) are defined to enhance the sensitivity of the search. Figure 2 shows a schematic view of the selection. The event selection, the expected signal performance and background estimations are described below.

Selection of Z →
The Z → candidates are identified by requiring two isolated same-flavour leptons (electrons or muons) satisfying the 'loose' criteria. The leading electron (muon) must satisfy E T (p T ) > 28 GeV. Opposite charges are required for the muon pairs but not for the electron pairs. Electrons are more susceptible to charge misidentification due to the conversions of photons from bremsstrahlung, especially at high E T . The dilepton invariant mass m is required to be consistent with the Z boson mass. For electrons, a fixed m window [83, 99] GeV is applied. To account for the effect of dimuon mass-resolution degradation at high transverse momentum (p T ), a p T -dependent mass window [85.6 GeV − 0.0117p T , 94.0 GeV + 0.0185p T ] is used for muons. For both electrons and muons, the mass windows are chosen to ensure that the Z → selection efficiency is approximately 95% and independent of the heavy resonance mass. Events with additional leptons passing the 'loose' criteria are vetoed.

VBF and ggF/DY categories
Signal events from VBF production possess unique kinematic signatures. In addition to the presence of a pair of vector bosons from the resonance decay, VBF events have two additional jets (referred to as tag-jets). These jets typically have a large separation in pseudorapidity and a large dijet invariant mass.
They offer a powerful means for background rejection. Hence for the H → ZZ and W → ZW searches, for which the VBF production is considered, two categories are introduced. Events are subject to the selection designed for VBF production first (VBF category), and if they fail, to the selection designed for ggF or DY production (ggF/DY category). For the G KK → ZZ search, the VBF selection is bypassed as VBF production is not considered.
For the VBF category, events are required to have two tag-jets identified from the small-R jets that fail the b-tagging described in Section 4. The two jets must be in opposite pseudorapidity hemispheres, i.e., η 1 · η 2 < 0, with a pseudorapidity separation |∆η tag j j | > 4.7, and have an invariant mass m tag j j > 770 GeV. Those values are chosen to optimise the search sensitivity to VBF signals for all masses considered. If there is more than one pair of jets satisfying these requirements, the one with the highest m tag j j value is chosen. These jets are not considered in the ZV → qq selection, and large-R jets lying ∆R < 1.5 of either of these two small-R jets are not considered either.
For the ggF/DY category, no tag-jets requirement is applied. Events in the VBF and ggF/DY categories are subject to an identical ZV → qq selection. All events not selected for the VBF category are passed to the ggF/DY selection. This includes events that contain jets satisfying the VBF requirements but which fail the ZV → qq selection, cf. Figure 2. In this case, jets excluded from consideration for the ZV → qq candidate in the VBF category are considered for the ggF/DY category.

Selection of ZV → qq
Identification of ZV → qq decays proceeds by applying the merged ZV → J selection followed by the resolved ZV → j j selection. The order is motivated by a smaller background expectation in the ZV → J final states. These selection criteria are summarised in Table 1 and are explained below.

Merged ZV → J selection
The ZV → J candidates are selected from the Z → events containing at least one large-R jet, among which the one with the highest p T is assumed to be from the V → qq decay. Events are further required to have min(p T , p J T )/m J > 0.3 for the H → ZZ search and > 0.35 for the W → ZW and G KK → ZZ searches. The criteria are optimised to reduce backgrounds while retaining high efficiencies for signals. The looser requirement for the H → ZZ search is motivated by the expected softer p T and p J T spectra resulting from a spin-0 resonance. This requirement suppresses background significantly at large m J while maintaining high efficiencies for signal events.
Figure 3(a) shows the distribution of the large-R jet mass m J of the H → ZZ search, comparing data with the expected backgrounds. Next the boson tagging discussed in Section 4 is applied to select the V → qq decays. Two signal regions are defined, one for events passing the 50% working point of the boson tagging requirement and the other for events failing the 50% but passing the 80% working point requirement. The former is called the high-purity (HP) SR, and the latter the low-purity (LP) SR. Background contributions are mostly from the Z+jets, top quark and diboson production (Section 5.5).

Resolved ZV → j j selection
The ZV → j j candidates are selected from the Z → events that have failed the ZV → J selection, by requiring at least two small-R 'signal' jets. The leading jet p T is required to be greater than 60 GeV.
Similarly to the merged ZV → J selection, the kinematic quantity p T 2 + p j j T 2 m j j is required to be greater than 0.4 for H → ZZ and 0.5 for W → ZW and G KK → ZZ. Here p j j T is the transverse momentum of the dijet system. Signal events are expected to have a dijet invariant mass, m j j , consistent with a V → qq decay. Figure 3(b) shows the m j j distributions of the H → ZZ search. The dijet invariant mass must be in the window [70, 105] GeV for Z → qq and in the window [62,97] GeV for W → qq. The choice of asymmetric windows around m V is motivated by the asymmetry of the m j j distribution expected from V → qq.
About 21% of the Z → qq decays have two b-quark jets, whereas the dominant background, Z+jets, has a smaller heavy-quark content. To exploit this, the ZV → j j candidates are divided into two signal regions: events with two b-tagged jets (b-tagged SR) and events with fewer than two b-tagged jets (untagged SR). Events with more than two b-tagged jets are rejected. No enhancement of b-tagged jets is expected from the W → qq decay, thus the two signal regions are combined for the W → ZW search.
Due to the small number of events, the two regions are also combined for the VBF category.

Signal regions, selection efficiencies and mass resolutions
Signal regions used in the search depend on the signal model under study. There are seven, six and four signal regions, respectively, for the searches for H → ZZ, W → ZW and G KK → ZZ.
Signal selection efficiencies are dependent on the signal model, the production process and the mass of the heavy resonance. As an example, Figure 4 shows the acceptance times efficiency (A × ) of the H → ZZ → qq search as a function of the Higgs boson mass for both ggF and VBF production. The ZV → J selection is more efficient for masses over approximately 500 GeV, while the ZV → j j selection contributes more at lower masses. This reflects the expected large Lorentz boost of the Z bosons from heavy Higgs boson decays and the higher priority given to the ZV → J selection. The A × values of the W → ZW and G KK → ZZ searches are similar with the exception of noticeable decreases, by about 40% from 2 TeV to 4 TeV, due to the merging of the electrons from the Z → ee decays for resonance masses above approximately 3 TeV.
Distributions of the invariant masses m J and m j j are used to search for potential signals. To mitigate the impact of muon momentum-resolution degradation at high p T , a scale factor of m Z /m µµ is applied to the four-momentum of the dimuon system in Z → µµ events, effectively fixing its mass to m Z = 91.  The reconstructed m J and m j j distributions have widths of 2-3% of the resonance mass for the Higgs boson and HVT W boson of VBF production for the entire mass range studied. Given the negligible intrinsic widths of these resonances, the widths largely reflect the detector resolution. For DY-produced HVT W signal, the widths of the m J and m j j distributions are slightly larger at 3-4% since the W boson in this production mode has an intrinsic width of approximately 2.6% of its mass. The widths of the m J and m j j distributions are 4% at 500 GeV for the bulk RS graviton signal with k/M Pl = 1 and rise to 8% at 5000 GeV. The increase can be attributed to the increase in the intrinsic width of the signal. For a given resonance mass, the m J distributions are narrower than those of m j j .

Data control regions and background estimation
The dominant backgrounds to the X → ZV → qq search are the Z+jets, top quark and diboson processes. Their contributions are estimated from a combination of MC and data-driven techniques. In all cases, the shapes of kinematic variables, including those of the final discriminants m J and m j j , are taken from MC simulations. The multijet background is estimated to be negligible.
The Z+jets events are expected to have smooth distributions of m J and m j j , while the signal events should exhibit resonance structures at the mass of the vector-boson V. Thus, a Z+jets control region (CR) is defined for every signal region by reversing the m J or m j j requirement. Events in the control regions are selected in exactly the same way as those in their corresponding signal regions except for the requirement on m J or m j j . For the ZV → J selection, the leading large-R jet mass is required to be outside the large-R jet mass window of the 80% working point of the boson tagging. For the ZV → j j selection, a requirement of 50 < m j j < 62 GeV or 105 < m j j < 150 GeV is applied. These CRs are dominated by the Z+jets contribution, with a purity higher than 96% in all regions, except for the b-tagged CR where the top quark and Z+jets contributions are comparable. They are therefore used to constrain its contribution in signal regions through simultaneous fits as discussed in Section 8.
Top quark production is a significant background source in the b-tagged signal region of the resolved ZV → j j selection. Its contribution is constrained using a top-quark-enhanced control region. Events in this control region must have two different-flavour leptons, eµ, with their invariant mass within [76, 106] GeV, and two b-tagged jets with their invariant mass in the range [50,150] GeV. The leading b-tagged jet is required to have p T > 60 GeV. This selection yields a sample of top quark events with a purity higher than 99%. This top quark control region is used to constrain top quark contributions for all signal regions.
Diboson production, mainly from the SM ZZ and ZW processes, is also a significant background source. The contribution from those processes is estimated completely from MC simulation.
As an example, Table 2 shows the numbers of events observed in the data and estimated from background processes in the seven signal regions of the H → ZZ → qq search. The numbers of background events are extracted through a background-only fit of the signal and control regions discussed in Section 8. Table 2: Numbers of events observed in the data and predicted for background processes from background-only fits to the signal and control regions (Section 8) in the seven signal regions of the H → ZZ → qq search. The numbers of signal events expected from the ggF and VBF production of a heavy Higgs boson with mass of 1 TeV are also shown. The signal yields are calculated using σ×B(H → ZZ) = 20 fb for both processes. The uncertainties combine statistical and systematic contributions. The fit constrains the background estimate towards the observed data, which reduces the total background uncertainty by correlating those from the individual backgrounds.
Events containing a X → ZV → ννqq signal are characterised by a hadronically decaying V boson recoiling against a large missing transverse momentum. Only the merged reconstruction of the V → qq decay is considered. The selection closely follows that of the X → ZV → qq search as illustrated in Figure 2 with additional requirements to remove multijet backgrounds but without the resolved selection.

Selection of ZV → ννqq
An initial selection is made by requiring E miss T > 250 GeV, and vetoing events with electrons or muons passing the 'loose' quality requirements. The multijet background, originating primarily from the presence of mismeasured jets, and non-collision backgrounds are suppressed by using a track-based missing transverse momentum, p miss T , as defined in Section 4. The requirements are p miss T > 50 GeV and the azimuthal separation between E miss T and p miss T directions satisfies ∆φ( E miss T , p miss T ) < 1. An additional requirement is imposed on the azimuthal separation between the directions of E miss T and the nearest small-R jet with min[∆φ( E miss T , small−R jet)] > 0.4. The multijet background is found to be negligible after these selections.
As in the qq search, both the VBF and ggF categories are considered in the signal event selection for the H → ZZ and W → ZW searches; the VBF selection is not motivated by the G KK → ZZ search. The VBF tag-jets are selected in the same manner as in the qq search, but they are required to satisfy m tag j j > 630 GeV and |∆η tag j j | > 4.7. Those values are chosen to optimise the search sensitivity for signals from VBF production, considering all signal masses. Events failing the VBF tag-jets selection are kept in the ggF category. For both the VBF and ggF categories, events are then subject to the selections described below.
Unlike the qq search, this search utilises only the merged selection of the V → qq decay, i.e., considering only the ZV → ννJ final state, since the resolved selection has little acceptance for the given high E miss T trigger threshold. The selection of the large-R jets and the boson tagging are identical to those of the qq search and are described in Section 4. Similarly, candidate events are split into HP and LP signal regions. The H → ZZ search and the W → ZW search use all the four signal regions, while the G KK → ZZ search bypasses the VBF selection, so has only two signal regions. The distributions of the mass and D 2 of the leading large-R jet before the boson tagging are shown in Figure 5. There is a small overprediction of the background at high mass region, but those events are not selected in any signal regions and sidebands. The selection criteria for the ννqq search are summarised in Table 3.
It is not possible to fully reconstruct the invariant mass of the ννJ system due to the presence of neutrinos in the final state, so the transverse mass is used as the final discriminant: where E T,J = m 2 J + p 2 T,J . The resolution of the transverse mass is about 10% of the signal mass and the ratio of the resolution to the signal mass has little dependence on the signal model or the resonance mass.  The signal selection efficiency strongly depends on the signal model and the mass of the resonance. As an example, Figure 6 shows the selection acceptance times efficiency as a function of m(H) for H → ZV → ννqq. Similar results are obtained for the W → ZW and G KK → ZZ searches.

Background estimation
The main backgrounds for this search arise from Z+jets, W+jets and tt production. Data control regions are defined to check the modelling of each contribution. As in the qq search, the shapes of kinematic variables, including the final discriminant m T , are taken from MC simulations.
The normalisation of the Z+jets background is determined using the Z+jets control region from the qq search (Section 5.5). Control regions specific to the ννqq search are also defined to constrain the W+jets and top quark backgrounds. Events in these control regions are selected in the same way as those in the signal regions except that they must have exactly one muon passing the 'tight' quality requirement. These events are then split into the top quark and W+jets control regions according to the number of b-tagged jets that do not overlap with the large-R jet. Events in the top quark control region must have at least one b-tagged jet and events in the W+jets control region must not have any b-tagged jets. For the W+jets control region, the jet mass requirement used in the signal region selection is inverted, to ensure that there is no resonant diboson signal contamination via X → WV → νqq decay, while for the top quark control region, the leading large-R jet is required to be consistent with the mass of the W boson.
For the W+jets and top quark CRs, the muon is treated as a missing particle to make the boson p T acceptance more similar to the signal region, and because the trigger does not include muons in the E miss T calculation. Specifically, E miss T,no µ (p miss T,no µ ) is computed by removing the muon p T contribution from the E miss T (p miss T ) calculation. The requirements on E miss T and p miss T in the signal regions are used for E miss T,no µ and p miss T,no µ respectively for the control regions. Table 4 summarises the numbers of events observed in the data and estimated from background processes in the signal regions from the background-only fit of both the signal and control regions discussed in Section 8. The expected numbers of events from both ggF and VBF production of a 1.6 TeV heavy Higgs boson, assuming σ × B(H → ZZ) = 6 fb, are also included for comparison. Table 4: Numbers of events predicted from background processes and observed in the data in the signal regions of the H → ZZ → ννqq search from the background-only fit of both the signal and control regions (Section 8). The numbers of signal events expected from a Higgs boson with mass of 1.6 TeV are also shown. The signal yields are calculated assuming σ × B(H → ZZ) = 6 fb for both ggF and VBF at 1.6 TeV. The uncertainties combine statistical and systematic contributions.

Systematic uncertainties
Systematic uncertainty sources impacting the search can be divided into three categories: experimental uncertainties related to the detector or to the reconstruction algorithms, uncertainties in the estimations of background contributions, and uncertainties in modelling the signal. Unless explicitly stated, the uncertainties quoted below are the uncertainties in the quantities themselves, not the impact on the search sensitivity.
The uncertainty in the integrated luminosity of the dataset is determined to be 3.2%. It is derived, following a methodology similar to the one detailed in Ref.
[80], from a preliminary calibration of the luminosity scale using a pair of x-y beam-separation scans performed in August 2015 and May 2016. This uncertainty is applied to the normalisation of the signal and also to background contributions whose normalisations are derived from MC simulations. A variation in the pile-up reweighting of MC events is included to cover the uncertainty in the ratio of the predicted and measured inelastic cross section in Ref. [81].
The efficiencies of the lepton triggers for events with selected leptons are high, nearly 100% in the electron channel and approximately 96% in the muon channel. The corresponding uncertainties are negligible. For the selection used in the ννqq search, the efficiency of the E miss T trigger is also close to 100% with negligible associated uncertainty. The modelling of the electron and muon reconstruction, identification and isolation efficiencies is studied with a tag-and-probe method using Z → events in data and simulation at √ s = 13 TeV [82, 83]. Small corrections are applied to the simulation to better model the performance seen in data. These corrections have associated uncertainties of the order of 1%. For the ννqq search, the uncertainty is instead in the efficiency of the lepton veto, and is at the sub-percent level. Uncertainties in the lepton energy (or momentum) scale and resolution are also taken into account.
Uncertainties in the jet energy scale and resolution for small-radius jets are estimated using MC simulation and in situ techniques [64]. For central jets (|η| < 2.0), the total uncertainty in the jet energy scale ranges from about 6% for jets with p T = 25 GeV to about 2% for p T = 1 TeV. There is also an uncertainty in the jet energy resolution [84], which ranges from 10% to 20% for jets with a p T of 20 GeV to less than 5% for jets with p T > 200 GeV. Uncertainties in the lepton and jet energy scales and resolutions are propagated into the uncertainty in E miss T . Uncertainties in the energy scale and resolution of the track soft term are also propagated into the uncertainty in E miss T [78], which is about 2% for the ννqq search. Uncertainties in the efficiency for tagging b-jets and in the rejection factor for light jets are determined from tt samples [85,86].
The uncertainties in the scale of the large-R jet p T , mass and D 2 are of the order of 2-5%. They are estimated using comparisons of data and simulation in Ref. [72]. An absolute uncertainty of 2% is assigned to the large-R jet energy resolution, and relative uncertainties of 20% and 15% are assigned to the resolution of the large-R jet mass and D 2 , respectively.
For the qq search the uncertainty in the spectrum of the final discriminant, m J or m j j , for Z+jets, is assessed by comparing the shape difference between data and MC simulations in the control regions. The uncerainty is found to be approximately 20% for m j j and varies from 3% at low mass to 20% at high mass for m J . For the ννqq search the uncertainty in the modelling of the m T distribution for Z+jets and W+jets production is assessed by comparing the nominal Sherpa MC sample with alternative samples generated by doubling or halving the renormalisation, resummation and factorisation scales independently and varying the matrix-element matching. The uncertainties are of the order of 5-30%. The approach is different from that used in the qq search, since it is difficult to obtain a pure control region for Z+jets or W+jets in the ννqq search.
An uncertainty in the shape of the m J or m T distribution for the tt background is derived by comparing the Powheg-Box sample with the distribution obtained using MadGraph5_aMC@NLO 2.2.2. This uncertainty is found to be approximately 15%. Additional systematic uncertainties are estimated by comparing the nominal sample showered with Pythia 6.428 using the P2012 tune to one showered with Herwig++ 2.7.1 [87] and using the UEEE5 underlying-event tune. The uncertainty varies between 10-40%. Samples of tt events with the factorisation and renormalisation scales doubled or halved are compared to the nominal, and the differences observed are taken as an additional uncertainty. These uncertainties are small, typically less than 3%.
The uncertainties in the diboson cross sections are estimated to be 10% [40,88], and the shape uncertainty is obtained by comparing MC samples generated by Sherpa and Powheg-Box.
The signals acceptance uncertainty due to PDF uncertainties is estimated by taking the uncertainty from the PDF error sets, CT10 for the heavy Higgs boson signal and NNPDF23LO for the HVT and bulk RS signals, and adding it in quadrature to the acceptance difference obtained with the use of alternative PDF sets: MMHT2014LO [89] and NNPDF30NLO [90] for the Higgs boson, and CT10 and MMHT2014LO for the HVT and bulk RS graviton. For the HVT signal, the uncertainty ranges from 1% (2%) to 6% (12%) for the qq (ννqq) search depending on the mass being tested. Similar results are obtained for the heavy Higgs boson signal, while for the bulk RS signal, the uncertainty is generally less than 1%. The signal acceptance uncertainty due to initial-and final-state radiation (ISR/FSR) is estimated by varying relevant parameters in the A14-NNPDF tune [38] for the HVT and bulk RS signals. This uncertainty ranges from less than 1% to 4% (6%) for the qq (ννqq) search depending on the mass being tested. The effect of the QCD scale uncertainty on the heavy Higgs boson signal acceptance is estimated by varying the factorisation and renormalisation scales. It ranges from 1% (2%) to 4% (8%) for the qq (ννqq) search depending on the mass being tested.

Statistical and fit procedures
The statistical analysis is based on the framework described in Refs. [91][92][93]. A binned likelihood function L(µ, θ) is constructed as a product of Poisson probability terms over all bins of the fit templates considered in the search. This function depends on the signal-strength parameter µ, a multiplicative factor applied to the theoretical signal production cross section, and θ, a set of nuisance parameters that encode the effects of systematic uncertainties in the signal and expected backgrounds, described in Section 7. The binning is chosen based on the expected mass resolution and numbers of events. The nuisance parameters are either free to float, or constrained using Gaussian terms defined by external studies. The likelihood function for the combination of the qq and ννqq channels is the product of the Poisson likelihoods of these individual channels. However, only one constraint term per common nuisance parameter is included in the product.
For each signal process considered, a simultaneous maximum-likelihood fit is performed to the observed distributions of the final discriminants in the signal regions, m J or m j j for the X → ZV → qq search and m T for the X → ZV → ννqq search, to extract the signal rate information. The shapes of the signal distributions are interpolated using the moment-morphing method [94] for intermediate mass points where simulated signal samples are not available. The product of acceptance and efficiency for each interpolated signal mass point is obtained from an interpolation with cubic splines between the simulated mass points. The Z+jets, W+jets and top quark control regions are included in the fit's likelihood calculation with one bin per region, i.e., using only their event-count information. Background contributions, including their shapes in the signal regions, are taken from MC simulations. However, they are allowed to vary independently within their uncertainties in each bin. Moreover, normalisation scale factors (SFs) are applied to the MC estimates of the Z+jets, W+jets and top quark contributions. These scale factors are free parameters in the fit and are therefore constrained by the data in both the signal and control regions. The diboson contribution is constrained to the theoretical estimate within the corresponding uncertainties.
In general, one SF is introduced per control region for its intended background component with the following exceptions to take into account different MC modellings in the different phase spaces of the same background component. One common Z+jets SF is used for HP and LP regions of the merged selection of both the ZV → qq and ZV → ννqq searches. However, independent SFs are used for the VBF and ggF categories. For the resolved ZV → qq selection, two independent Z+jets SFs are used for the b-tagged and untagged regions: one for the Z+heavy-flavour component and the other for the Z+light-jet component (Section 3). The ZV → qq search defines a top quark control region and consequently applies one top quark SF to all regions. The ZV → ννqq search shares the Z+jets control regions with ZV → qq and hence applies the same SFs to its Z+jets background contributions. The ZV → ννqq search has its own W+jets and top quark control regions. Similar to Z+jets, one common W+jets SF is used for HP and LP regions, but independent SFs are used for the VBF and ggF categories. All the top quark control regions in the ZV → ννqq search share one common SF and it is distinguished from the top quark SF in the ZV → qq search.
The test statistic q µ is defined as the profile likelihood ratio [95], q µ = −2ln(L(µ,θ µ )/L(μ,θ)), whereμ and θ are the values of the parameters that maximise the likelihood function (with the constraint 0≤μ ≤ µ), andθ µ are the values of the nuisance parameters that maximise the likelihood function for a given value of µ. The test statistic q µ is used to measure the compatibility of the observed data with the background-only hypothesis. In the absence of a signal, constraints on the production of a heavy resonance decaying into ZV pairs are derived. The exclusion limits are calculated with a modified frequentist method [96], also known as CL s , in the asymptotic approximation [97] for resonance masses below 2 TeV. Above 2 TeV, the small number of events makes the asymptotic approximation unreliable and the limits are calculated using pseudo-experiments. All limits are set at the 95% confidence level (CL).

Limits on the production of heavy resonances
The observed distributions of m J or m j j from the H → ZZ → qq search are compared with the background estimates in Figure 7 for the three signal regions of the VBF category and in Figure 8 for the four signal regions of the ggF category. Similarly, Figures 9 and 10 show the m T distributions from the H → ZZ → ννqq search in the signal regions of the VBF and ggF categories, respectively. There are no events in the data beyond the ranges of the distributions shown. Background contributions shown are obtained from background-only fits described previously. The numbers of events observed and estimated in control regions are summarised in Figure 11.
The data distributions are reasonably well reproduced by the estimated background contributions in all these distributions. Similar levels of agreement are observed for the W → ZW and G KK → ZZ searches. The largest difference between the observed data and the SM background prediction is in the H → ZZ → qq search, high-purity signal region, as shown in Figure 8(a). A deficit in data is seen around m J of 800 GeV, with a local significance of 3.0σ, and by taking into account the look-elsewhere effect [98] in the full m J spectrum, a global significance of 1.9σ. Upper limits are set on the product of the production cross section of new resonances and their decay branching ratio to ZV, σ × B(X → ZV).       Limits are presented for both the ggF and VBF productions of H → ZZ in Figure 12 in the resonance mass range between 300 GeV and 3 TeV. The observed limit on σ × B(H → ZZ) varies from 1.7 (0.42) pb at 300 GeV to 1.4 (1.1) fb at 3 TeV for ggF (VBF) H → ZZ production.    Figure 13 shows the limits on σ × B(W → ZW) for DY and VBF production of a W boson in the HVT model. The observed limit ranges from 5.7 pb at 300 GeV to 1.3 fb at 5 TeV for DY production and from 0.98 pb at 300 GeV to 2.8 fb at 4 TeV for VBF production. The theoretical predictions of the HVT Model A, Model B and VBF Model are overlaid for comparison. The observed limits exclude an HVT W boson produced in the DY process lighter than 2.9 TeV for Model A and 3.2 TeV for Model B, while none of the HVT model space can be excluded for the VBF process with the current sensitivity of the analysis. HVT models generally predict a near-by neutral vector boson (Z ) that decays into a Z boson and a SM Higgs boson h. The Z production could contaminate both the signal and control regions of the HVT W search. Assuming the Z → Zh production having the same cross section as the observed upper limit on σ × B(W → ZW), the impacts of the potential contamination to the W → ZW search are found to be negligible.

Effects of systematic uncertainties
The effects of systematic uncertainties are studied for hypothesised signals using the signal-strength parameter µ. The relative uncertainties in the best-fit µ value from the leading sources of systematic uncer-tainty are shown in Table 5 for ggF H → ZZ production with m(H) = 600 GeV and 1.2 TeV. Apart from the statistical uncertainties in the data, the uncertainties with the largest impact on the sensitivity of the searches are from the size of the MC samples, measurements of small-R and large-R jets, E miss T measurement, background modelling and luminosity. For signals with higher mass, the data statistical uncertainty becomes dominant. The effects of systematic uncertainties for the other searches are similar to those shown for the ggF H → ZZ search.

Conclusion
Searches for heavy resonances decaying into ZZ or ZW are performed using proton-proton collision data produced at √ s = 13 TeV and recorded by the ATLAS detector at the LHC in 2015 and 2016. The data correspond to a combined integrated luminosity of 36.1 fb −1 . The searches explore the final states with one Z boson decaying either into a pair of charged leptons or into a pair of neutrinos, and the other Z boson or the W boson decaying into a pair of quarks, identified either as two separate jets or as one largeradius jet. Moreover, the searches are performed in two event categories, targeting VBF and non-VBF production of the resonances.
The data are found to be consistent with the SM background predictions and no evidence of heavy resonance production is observed. Upper limits on the production cross section times branching ratio,  [20] ATLAS Collaboration, The ATLAS Experiment at the CERN Large Hadron Collider, JINST 3 (2008) S08003.