Search for heavy $ZZ$ resonances in the $\ell^+\ell^-\ell^+\ell^-$ and $\ell^+\ell^-\nu\bar\nu$ final states using proton proton collisions at $\sqrt{s}= 13$ TeV with the ATLAS detector

A search for heavy resonances decaying into a pair of $Z$ bosons leading to $\ell^+\ell^-\ell^+\ell^-$ and $\ell^+\ell^-\nu\bar\nu$ final states, where $\ell$ stands for either an electron or a muon, is presented. The search uses proton proton collision data at a centre-of-mass energy of 13 TeV corresponding to an integrated luminosity of 36.1 fb$^{-1}$ collected with the ATLAS detector during 2015 and 2016 at the Large Hadron Collider. Different mass ranges for the hypothetical resonances are considered, depending on the final state and model. The different ranges span between 200 GeV and 2000 GeV. The results are interpreted as upper limits on the production cross section of a spin 0 or spin 2 resonance. The upper limits for the spin 0 resonance are translated to exclusion contours in the context of Type I and Type II two-Higgs-doublet models, while those for the spin 2 resonance are used to constrain the Randall Sundrum model with an extra dimension giving rise to spin 2 graviton excitations.


Introduction
In 2012, the ATLAS and CMS Collaborations at the LHC discovered a new particle [1,2], an important milestone in the understanding of the mechanism of electroweak (EW) symmetry breaking [3][4][5]. Both experiments have confirmed that the spin, parity and couplings of the new particle are consistent with those predicted for the Standard Model (SM) Higgs boson [6-8] (denoted by h throughout this paper). They measured its mass to be m h = 125.09 ± 0.21(stat) ± 0.11(syst) GeV [9] and reported recently on a combination of measurements of its couplings to other SM particles [10].
One important question is whether the newly discovered particle is part of an extended scalar sector as postulated by various extensions to the Standard Model such as the two-Higgs-doublet model (2HDM) [11]. These extensions predict additional Higgs bosons, motivating searches in an extended range of mass. This paper reports on two searches for a heavy resonance decaying into two SM Z bosons, encompassing the final states Z Z → ℓ + ℓ − ℓ + ℓ − and Z Z → ℓ + ℓ − νν where ℓ stands for either an electron or a muon and ν stands for all three neutrino flavours. These final states are referred to as ℓ + ℓ − ℓ + ℓ − and ℓ + ℓ − νν respectively.
It is assumed that an additional Higgs boson (denoted as H throughout this paper) would be produced predominantly via gluon-gluon fusion (ggF) and vector-boson fusion (VBF) processes, but that the ratio of the two production mechanisms is unknown in the absence of a specific model. For this reason, the results are interpreted separately for the ggF and VBF production modes, with events being classified into ggF-and VBF-enriched categories in both final states, as discussed in Sections 5 and 6. With good mass resolution and a high signal-to-background ratio, the ℓ + ℓ − ℓ + ℓ − final state is well suited to a search for a narrow resonance with mass m H between 200 GeV and 1200 GeV. The ℓ + ℓ − νν search covers the 300 GeV < m H < 1400 GeV range and dominates at high masses due to its larger branching ratio.
These searches look for an excess in distributions of the four-lepton invariant mass, m 4ℓ , for the ℓ + ℓ − ℓ + ℓ − final state, and the transverse invariant mass, m T , for the ℓ + ℓ − νν final state, as the escaping neutrinos do not allow the full reconstruction of the final state. The transverse invariant mass is defined as: where m Z is the mass of the Z boson, p ℓℓ T is the transverse momentum of the lepton pair and ì E miss T is the missing transverse momentum, with magnitude E miss T . In the absence of such an excess, limits on the production rate of different signal hypotheses are obtained from a simultaneous likelihood fit to the two mass distributions. The first hypothesis is the ggF and VBF production of a heavy Higgs boson (spin-0 resonance) under the narrow-width approximation (NWA). The upper limits on the production rate of a heavy Higgs boson are then translated into exclusion contours in the context of the two-Higgs-doublet model. As several theoretical models favour non-negligible natural widths, large-width assumption (LWA) models, assuming widths of 1%, 5% and 10% of the resonance mass, are also studied. The interference between the heavy scalar and the SM Higgs boson as well as between the heavy scalar and the gg → Z Z continuum background are taken into account in this study. Limits are also set on the Randall-Sundrum (RS) model [12,13] with a warped extra dimension giving rise to a spin-2 Kaluza-Klein (KK) excitation of the graviton G KK .
Other searches for diboson resonances decaying into WW or Z Z or W Z have been performed by AT-LAS [14-16] and CMS [17][18][19].
With a significant increase in integrated luminosity and an improved discovery potential from the higher parton luminosities [20] at a centre-of-mass energy of √ s = 13 TeV as compared to √ s = 8 TeV, the results of this paper improve upon previous results published by the ATLAS Collaboration from a search for an additional heavy Higgs boson [21]. Results of a similar search from the data collected at the LHC with √ s = 8 TeV have also been reported by the CMS Collaboration [22].
the MSTW2008LO [30] parton distribution functions (PDF) set. The simulated events are weighted to reproduce the observed distribution of the mean number of interactions per bunch crossing in data (pileup reweighting). The properties of the bottom and charm hadron decays were simulated by the E G v1.2.0 program [31].
Heavy spin-0 resonance production was simulated using the P -B v2 [32] MC event generator. Gluon-gluon fusion and vector-boson fusion production modes were calculated separately with matrix elements up to next-to-leading order (NLO) in QCD. P -B was interfaced to P 8.212 [33] for parton showering and hadronisation, and for decaying the Higgs boson into the H → Z Z → ℓ + ℓ − ℓ + ℓ − or H → Z Z → ℓ + ℓ − νν final states. The CT10 PDF set [34] was used for the hard process. Events from ggF and VBF production were generated in the 300 GeV < m H < 1600 GeV mass range under the NWA, using a step of 100 (200) GeV up to (above) 1000 GeV in mass. For the ℓ + ℓ − ℓ + ℓ − final state, due to the sensitivity of the analysis at lower masses, events were also generated for m H = 200 GeV.
In addition, events from ggF production with a boson width of 5%, 10% and 15% of the scalar mass m H were generated with M G 5_aMC@NLO v2.3.2 [35] interfaced to P 8.210 for parton showering and hadronisation for both final states. For the ℓ + ℓ − ℓ + ℓ − final state, the m 4ℓ distribution is parameterised analytically as described in Section 5.3, and the samples with a width of 15% of m H are used to validate the parameterisation. For the ℓ + ℓ − νν final state, a reweighing procedure as described in Section 6.3 is used on fully simulated events to obtain the reconstructed m T distribution at any value of mass and width tested. To have a better description of the jet multiplicity, M G 5_aMC@NLO was also used to generate events for the process pp → H + ≥ 2 jets at NLO QCD accuracy with the FxFx merging scheme [36].
The fraction of the ggF events that enter into the VBF-enriched category is estimated from the M -G 5_aMC@NLO simulation.
Spin-2 Kaluza-Klein gravitons from the Bulk Randall-Sundrum model [37] were generated with M -G 5_aMC@NLO at leading order (LO) in QCD. The dimensionless coupling k/M Pl , whereM Pl = M Pl / √ 8π is the reduced Planck scale and k is the curvature scale of the extra dimension, is set to 1. In this configuration, the width of the resonance is expected to be ∼ 6% of its mass.
Mass points between 600 GeV and 2 TeV with 200 GeV spacing were generated for the ℓ + ℓ − νν final state. These samples were processed through a fast detector simulation [26] that uses a parameterisation of the response of electromagnetic and hadronic calorimeters [38], while the response of the ID and MS detectors is fully simulated.
The qq → Z Z background for the ℓ + ℓ − νν final state was simulated by the P -B v2 event generator [32] and interfaced to P 8.186 [28] for parton showering and hadronisation. The CT10 PDF set [34] was used for hard-scattering processes. Next-to-next-to-leading-order (NNLO) QCD and NLO EW corrections are included [39-41] as a function of the invariant mass m Z Z of the Z Z system. For the ℓ + ℓ − ℓ + ℓ − final state, this background was simulated with the S v2. EW ) was generated using S , where the process Z Z Z → 4ℓqq is also taken into account.
The gg → Z Z production was modelled by S v2.1.1 at LO in QCD for the ℓ + ℓ − ℓ + ℓ − final state and by 2VV [49] for the ℓ + ℓ − νν final state, both including the off-shell h boson contribution and the interference between the h and Z Z backgrounds. The K-factor accounting for higher-order QCD effects for the gg → Z Z continuum production was calculated for massless quark loops [50][51][52] in the heavytop-quark approximation [53], including the gg → H * → Z Z process [54]. Based on these studies, a constant K-factor of 1.7 is used, and a relative uncertainty of 60% is assigned to the normalisation in both searches.
The WW and W Z diboson events were simulated by P -B , using the CT10 PDF set and P 8.186 for parton showering and hadronisation. The production cross section of these samples is predicted at NLO in QCD.
Events containing a single Z boson with associated jets were simulated using the S v2.  [55].
The triboson backgrounds Z Z Z, W Z Z, and WW Z with fully leptonic decays and at least four prompt charged leptons were modelled using S v2.1.1. For the fully leptonic tt + Z background, with four prompt leptons originating from the decays of the top quarks and Z boson, M G 5_aMC@NLO was used. The tt background, as well as the single-top and Wt production, were modelled using P -B v2 interfaced to P 6.428 [56] with the Perugia 2012 [57] set of tuned parameters for parton showering and hadronisation, to PHOTOS [58] for QED radiative corrections and to T [59,60] for the simulation of τ-lepton decays.
In order to study the interference treatment for the LWA case, samples containing the gg → Z Z continuum background (B) as well as its interference (I) with a hypothetical heavy scalar (S) were used and are referred to as SBI samples hereafter. In the ℓ + ℓ − ℓ + ℓ − final state the MCFM NLO event generator [61], interfaced to P 8.212, was used to produce SBI samples where the width of the heavy scalar is set to 15% of its mass, for masses of 200, 300, 400, 500, 600, 800, 1000, 1200 and 1400 GeV. Background-only samples were also generated with the MCFM event generator, and are used to extract the signal-plus-interference term (SI) by subtracting them from the aforementioned SBI samples. For the ℓ + ℓ − νν final state, the SBI samples were generated with the 2VV event generator. The samples include signal events with a scalar mass of 400, 700, 900, 1200 and 1500 GeV.
The selection used combines the likelihood with the number of track hits and defines two working points (WP) which are used in the analyses presented here. The ℓ + ℓ − ℓ + ℓ − analysis uses a "loose" WP, with an efficiency ranging from 90% for transverse momentum p T = 20 GeV to 96% for p T > 60 GeV. A "medium" WP was chosen for the ℓ + ℓ − νν analysis with an efficiency increasing from 82% at p T = 20 GeV to 93% for p T > 60 GeV. The electron's transverse momentum is computed from the cluster energy and the track direction at the interaction point.
Muons are formed from tracks reconstructed in the ID and MS, and their identification is primarily based on the presence of the track or track segment in the MS [63]. If a complete track is present in both the ID and the MS, a combined muon track is formed by a global fit using the hit information from both the ID and MS detectors (combined muon), otherwise the momentum is measured using the ID, and the MS track segment serves as identification (segment-tagged muon). The segment-tagged muon is limited to the centre of the barrel region (|η| < 0.1) which has reduced MS geometrical coverage. Furthermore, in this central region an ID track with p T > 15 GeV is identified as a muon if its calorimetric energy deposition is consistent with a minimum-ionising particle (calorimeter-tagged muon). In the forward region (2.5 < |η| < 2.7) with limited or no ID coverage, the MS track is either used alone (stand-alone muon) or combined with silicon hits, if found in the forward ID (combined muon). The ID tracks associated with the muons are required to have a minimum number of associated hits in each of the ID subdetectors to ensure good track reconstruction. The stand-alone muon candidates are required to have hits in each of the three MS stations they traverse. A "loose" muon identification WP, which uses all muon types and has an efficiency of 98.5%, is adopted by the ℓ + ℓ − ℓ + ℓ − analysis. For the ℓ + ℓ − νν analysis a "medium" WP is used, which only includes combined muons and has an efficiency of 97%.
Jets are reconstructed using the anti-k t algorithm [64] with a radius parameter R = 0.4 implemented in the F J package [65], and positive-energy clusters of calorimeter cells as input. The algorithm suppresses noise and pile-up by keeping only cells with a significant energy deposit and their neighbouring cells. Jets are calibrated using a dedicated scheme designed to adjust, on average, the energy measured in the calorimeter to that of the true jet energy [66]. The jets used in this analysis are required to satisfy p T > 20 GeV and |η| < 4.5. To reduce the number of jet candidates originating from pile-up vertices, an additional requirement that uses the track and vertex information inside a jet is imposed on jets with p T < 60 GeV and |η| < 2.4 [67].
Jets containing b-hadrons, referred to as b-jets, are identified by the long lifetime, high mass and decay multiplicity of b-hadrons, as well as the hard b-quark fragmentation function. The ℓ + ℓ − νν analysis identifies b-jets of p T > 20 GeV and |η| < 2.5 using an algorithm that achieves an identification efficiency of about 85% in simulated tt events, with a rejection factor for light-flavour jets of about 33 [68, 69].
Selected events are required to have at least one vertex with two associated tracks with p T > 400 MeV, and the primary vertex is chosen to be the vertex reconstructed with the largest p 2 T . As lepton and jet candidates can be reconstructed from the same detector information, a procedure to resolve overlap ambiguities is applied. If an electron and a muon share the same ID track, the muon is selected unless it is calorimeter-tagged and does not have a MS track, or is a segment-tagged muon, in which case the electron is selected. Reconstructed jets which overlap with electrons (muons) in a cone of size ∆R ≡ (∆η) 2 + (∆φ) 2 = 0.2 (0.1) are removed.
The missing transverse momentum ì E miss T , which accounts for the imbalance of visible momenta in the plane transverse to the beam axis, is computed as the negative vector sum of the transverse momenta of all identified electrons, muons and jets, as well as a "soft term", accounting for unclassified soft tracks and energy clusters in the calorimeters [70]. This analysis uses a track-based soft term, which is built by combining the information provided by the ID and the calorimeter, in order to minimise the effect of pile-up which degrades the E miss T resolution. The soft term is computed using the momenta of the tracks associated with the primary vertex, while the jet and electron momenta are computed at the calorimeter level to allow the inclusion of neutral particles. Jet-muon overlap is accounted for in the E miss T calculation. This corrects for fake jets due to pile-up close to muons and double-counted jets from muon energy losses.

Event selection
Four-lepton events are selected and initially classified according to the lepton flavours: 4µ, 2e2µ, 4e, called "channels" hereafter. They are selected with single-lepton, dilepton and trilepton triggers, with the dilepton and trilepton ones including electron(s)-muon(s) triggers. Single-electron triggers apply "medium" or "tight" likelihood identification, whereas multi-electron triggers apply "loose" or "medium" identification. For the bulk of the data, recorded in 2016, the lowest p T threshold for the single-electron (muon) triggers used is set to 26 (26) GeV, for the dielectron (dimuon) triggers to 15 (10) GeV and for the trielectron (trimuon) triggers to 12 (6) GeV. For the data collected in 2015, the instantaneous luminosity was lower so the trigger thresholds were lower; this increases the signal efficiency by less than 1%. Globally, the trigger efficiency for signal events passing the final selection requirements is about 98%.
In each channel, four-lepton candidates are formed by selecting a lepton-quadruplet made out of two same-flavour, opposite-sign lepton pairs, selected as described in Section 4. Each electron (muon) must satisfy p T > 7 (5) GeV and be measured in the pseudorapidity range of |η| < 2.47 (2.7). The highest-p T lepton in the quadruplet must satisfy p T > 20 GeV, and the second (third) lepton in p T order must satisfy p T > 15 GeV (10 GeV). In the case of muons, at most one calorimeter-tagged, segment-tagged or stand-alone (2.5 < |η| < 2.7) muon is allowed per quadruplet.
If there is ambiguity in assigning leptons to a pair, only one quadruplet per channel is selected by keeping the quadruplet with the lepton pairs closest (leading pair) and second closest (subleading pair) to the Z boson mass, with invariant masses referred to as m 12 and m 34 respectively. In the selected quadruplet, m 12 is required to be 50 GeV < m 12 < 106 GeV, while m 34 is required to be less than 115 GeV and greater than a threshold that is 12 GeV for m 4ℓ ≤ 140 GeV, rises linearly from 12 GeV to 50 GeV with m 4ℓ in the interval of [140 GeV, 190 GeV] and is fixed to 50 GeV for m 4ℓ > 190 GeV.
Selected quadruplets are required to have their leptons separated from each other by ∆R > 0.1 if they are of the same flavour and by ∆R > 0.2 otherwise. For 4µ and 4e quadruplets, if an opposite-charge same-flavour lepton pair is found with m ℓℓ below 5 GeV, the quadruplet is removed to suppress the contamination from J/ψ mesons. If multiple quadruplets from different channels are selected at this point, only the quadruplet from the channel with the highest expected signal rate is retained, in the order: 4µ, 2e2µ, 4e.
The Z + jets and tt background contributions are reduced by imposing impact-parameter requirements as well as track-and calorimeter-based isolation requirements on the leptons. The transverse impactparameter significance, defined as the impact parameter calculated with respect to the measured beam line position in the transverse plane divided by its uncertainty, |d 0 |/σ d 0 , for all muons (electrons) is required to be lower than 3 (5). The normalised track-isolation discriminant, defined as the sum of the transverse momenta of tracks, inside a cone of size ∆R = 0.3 (0.2) around the muon (electron) candidate, excluding the lepton track, divided by the lepton p T , is required to be smaller than 0.15. The larger muon cone size corresponds to that used by the muon trigger. Contributions from pile-up are suppressed by requiring tracks in the cone to originate from the primary vertex. To retain efficiency at higher p T , the track-isolation cone size is reduced to 10 GeV/p T for p T above 33 (50) GeV for muons (electrons).
The relative calorimetric isolation is computed as the sum of the cluster transverse energies E T , in the electromagnetic and hadronic calorimeters, with a reconstructed barycentre inside a cone of size ∆R = 0.2 around the candidate lepton, divided by the lepton p T . The clusters used for the isolation are the same as those for reconstructing jets. The relative calorimetric isolation is required to be smaller than 0.3 (0.2) for muons (electrons). The measured calorimeter energy around the muon (inside a cone of size ∆R = 0.1) and the cells within 0.125 × 0.175 in η × φ around the electron barycentre are excluded from the respective sums. The pile-up and underlying-event contributions to the calorimeter isolation are subtracted event by event [71]. For both the track-and calorimeter-based isolation requirements, any contribution arising from other leptons of the quadruplet is subtracted.
An additional requirement based on a vertex-reconstruction algorithm, which fits the four-lepton candidates with the constraint that they originate from a common vertex, is applied in order to further reduce the Z + jets and tt background contributions. A loose cut of χ 2 /ndof < 6 for 4µ and < 9 for the other channels is applied, which retains a signal efficiency larger than 99% in all channels.
The QED process of radiative photon production in Z boson decays is well modelled by simulation. Some of the final-state-radiation (FSR) photons can be identified in the calorimeter and incorporated into the ℓ + ℓ − ℓ + ℓ − analysis. The strategy to include FSR photons into the reconstruction of Z bosons is the same as in Run 1 [21]. It consists of a search for collinear (for muons) and non-collinear FSR photons (for muons and electrons) with only one FSR photon allowed per event. After the FSR correction, the lepton four-momenta of both dilepton pairs are recomputed by means of a Z-mass-constrained kinematic fit. The fit uses a Breit-Wigner Z boson line-shape and a single Gaussian function per lepton to model the momentum response function with the Gaussian width set to the expected resolution for each lepton. The Z-mass constraint is applied to both Z candidates, and improves the m 4ℓ resolution by about 15%.
In order to be sensitive to the VBF production mode, events are classified into four categories: one for the VBF production mode and three for the ggF production mode, one for each of the three channels. If an event has two or more jets with p T greater than 30 GeV, with the two leading jets being well separated in η, |∆η jj | > 3.3, and having an invariant mass m jj > 400 GeV, this event is classified into the VBF-enriched category; otherwise the event is classified into one of the ggF-enriched categories. Such classification is used only in the search for a heavy scalar produced with the NWA.
The signal acceptance, defined as the ratio of the number of reconstructed events passing the analysis requirements to the number of simulated events in each category, is shown in Table 1, for the ggF and VBF production modes as well as different resonance masses. The contribution from final states with τ leptons decaying into electrons or muons is found to be negligible.

Background estimation
The main background component in the H → Z Z → ℓ + ℓ − ℓ + ℓ − final state, accounting for 97% of the total expected background events, is non-resonant Z Z production. This arises from quark-antiquark Table 1: Signal acceptance for the ℓ + ℓ − ℓ + ℓ − analysis, for both the ggF and VBF production modes and resonance masses of 300 and 600 GeV. The acceptance is defined as the ratio of the number of reconstructed events after all selection requirements to the number of simulated events for each channel/category.

Mass
Production annihilation (86%), gluon-initiated production (10%) and a small contribution from EW vector-boson scattering (1%). The last is more important in the VBF-enriched category, where it accounts for 16% of the total expected background. These backgrounds are all modelled by MC simulation as described in Section 3. Additional background comes from the Z + jets and tt processes, which contribute at the percent level and decrease more rapidly than the non-resonant Z Z production as a function of m 4ℓ . These backgrounds are estimated using data where possible, following slightly different approaches for final states with a dimuon (ℓℓ + µµ) or a dielectron (ℓℓ + ee) subleading pair [72].
The ℓℓ + µµ non-Z Z background comprises mostly tt and Z + jets events, where in the latter case the muons arise mostly from heavy-flavour semileptonic decays and to a lesser extent from π/K in-flight decays. The contribution from single-top production is negligible. The normalisations of the Z + jets and tt backgrounds are determined using fits to the invariant mass of the leading lepton pair in dedicated data control regions. The control regions are formed by relaxing the χ 2 requirement on the vertex fit, and by inverting and relaxing isolation and/or impact-parameter requirements on the subleading muon pair. An additional control region (eµµµ) is used to improve the tt background estimate. Transfer factors to extrapolate from the control regions to the signal region are obtained separately for tt and Z + jets using simulated events. The transfer factors have a negligible impact on the m 4ℓ shape of the ℓℓ + µµ background.
The main background for the ℓℓ + ee process arises from the misidentification of light-flavour jets as electrons, photon conversions and the semileptonic decays of heavy-flavour hadrons. The ℓℓ + ee controlregion selection requires the electrons in the subleading lepton pair to have the same charge, and relaxes the identification and isolation requirements on the electron candidate, denoted X, with the lower transverse momentum. The heavy-flavour background is completely determined from simulation, whereas the lightflavour and photon-conversion background is obtained with the sPlot [73] method, based on a fit to the number of hits in the innermost ID layer in the data control region. Transfer factors for the light-flavour jets and converted photons, obtained from simulated samples, are corrected using a Z + X control region and then used to extrapolate the extracted yields to the signal region. Both the yield extraction and the extrapolation are performed in bins of the transverse momentum of the electron candidate and the jet multiplicity.
The W Z production process is included in the data-driven estimates for the ℓℓ + ee final states, while it is added from simulation for the ℓℓ + µµ final states. The contributions from ttV (where V stands for either a W or a Z boson) and triboson processes are minor and taken from simulated samples.

Signal and background modelling
The parameterisation of the reconstructed four-lepton invariant mass m 4ℓ distribution for signal and background is based on the MC simulation and used to fit the data.
In the case of a narrow resonance, the width in m 4ℓ is determined by the detector resolution, which is modelled by the sum of a Crystal Ball (C) function [74, 75] and a Gaussian (G) function: The Crystal Ball and the Gaussian functions share the same peak value of m 4ℓ (µ), but have different resolution parameters, σ C and σ G . The α C and n C parameters control the shape and position of the non-Gaussian tail and the parameter f C ensures the relative normalisation of the two probability density functions. To improve the stability of the parameterisation in the full mass range considered, the parameter n C is set to a fixed value. The bias in the extraction of signal yields introduced by using the analytical function is below 1.5%. The function parameters are determined separately for each final state using signal simulation, and fitted to first-and second-degree polynomials in scalar mass m H to interpolate between the generated mass points. The use of this parameterisation for the function parameters introduces an extra bias in the signal yield and m H extraction of about 1%. An example of this parameterisation is illustrated in Figure 1, where the left plot shows the mass distribution for simulated samples at m H = 300, 600, 900 GeV and the right plot shows the RMS of the m 4ℓ distribution in the range considered for this search.
In the case of the LWA, the particle-level line-shape of m 4ℓ is derived from a theoretical calculation, as described in Ref.
[76], and is then convolved with the detector resolution, using the same procedure as for the modelling of the narrow resonance.
The m 4ℓ distribution for the Z Z continuum background is taken from MC simulation, and parameterised by an empirical function for both the quark-and gluon-initiated processes: The function's first part, f 1 , covers the low-mass part of the spectrum where one of the Z bosons is offshell, while f 2 models the Z Z threshold around 2·m Z and f 3 describes the high-mass tail. The transition between low-and high-mass parts is performed by the Heaviside step function H(x) around m 0 = 240 GeV. The continuity of the function around m 0 is ensured by the normalisation factor C 0 that is applied to the low-mass part. Finally, a i , b i and c i are shape parameters which are obtained by fitting the m 4ℓ distribution in simulation for each category. The uncertainties in the values of these parameters from the fit are found to be negligible. The MC statistical uncertainties in the high-mass tail are taken into account by assigning a 1% uncertainty to c 4 .
The m 4ℓ shapes are extracted from simulation for most background components (ttV, VVV, ℓℓ + µµ and heavy-flavour hadron component of ℓℓ + ee), except for the light-flavour jets and photon conversions in the case of ℓℓ + ee background, which is taken from the control region as described in Section 5.2.

Interference modelling
The gluon-initiated production of a heavy scalar H, the SM h and the gg → Z Z continuum background all share the same initial and final state, and thus lead to interference terms in the total amplitude. Theoretical calculations described in Ref.
[77] have shown that the effect of interference could modify the integrated cross section by up to O(10%), and this effect is enhanced as the width of the heavy scalar increases. Therefore, a search for a heavy scalar Higgs boson in the LWA case must properly account for two interference effects: the interference between the heavy scalar and the SM Higgs boson (denoted by H-h) and between the heavy scalar and the gg → Z Z continuum (denoted by H-B).
Assuming that H and h bosons have similar properties, as postulated by the 2HDM, they have the same production and decay amplitudes and therefore the only difference between the signal and interference terms in the production cross section comes from the propagator. Hence, the acceptance and resolution of the signal and interference terms are expected to be the same. The H-h interference is obtained by reweighting the particle-level line-shape of generated signal events using the following formula: where 1/ s − s H(h) is the propagator for a scalar (H or h). The particle-level line-shape is then convolved with the detector resolution function, and the signal and interference acceptances are assumed to be the same. In order to extract the H-B interference contribution, signal-only and background-only samples are subtracted from the generated SBI samples. The extracted particle-level m 4ℓ distribution for the H-B interference term is then convolved with the detector resolution. Figure 2 shows the overlay of the signal, both interference effects and the total line-shape for different mass and width hypotheses assuming the couplings expected in the SM for a heavy Higgs boson. As can be seen, the two interference effects tend to cancel out, and the total interference yield is for the most part positive, enhancing the signal.

Event selection
The analysis is designed to select Z Z → ℓ + ℓ − νν events (with ℓ = e, µ), where the missing neutrinos are identified by a large E miss T , and to discriminate against the large Z + jets, W Z and top-quark backgrounds. Events are required to pass either a single-electron or a single-muon trigger, where different p T thresholds are used depending on the instantaneous luminosity of the LHC. For the 2015 data the electron and muon triggers had p T thresholds of 24 GeV and 20 GeV respectively, while for 2016 the muon trigger threshold was increased to 24 GeV. For both triggers, the threshold is set to 26 GeV when the instantaneous luminosity exceeds the value of 10 34 cm −2 s −1 . The trigger efficiency for signal events passing the final selection is about 99%.
Events are selected if they contain exactly two opposite-charge leptons of the same flavour and "medium" identification, with the more energetic lepton having p T > 30 GeV and the other one having p T > 20 GeV. The same impact-parameter significance criteria as defined in Section 5.1 are applied to the selected leptons. Track-and calorimeter-based isolation criteria as defined in Section 5.1 are also applied to the leptons, but in this analysis the isolation criteria are optimised by adjusting the isolation threshold so that their selection efficiency is 99%. If an additional lepton with p T > 7 GeV and "loose" identification is found, the event is rejected to reduce the amount of W Z background. In order to select leptons originating from the decay of a Z boson, the invariant mass of the pair is required to be in the range 76 to 106 GeV. Moreover, since a Z boson originating from the decay of a high-mass particle is boosted, the two leptons are required to be produced with an angular separation of ∆R ℓℓ < 1.8.
Events with neutrinos in the final state are selected by requiring E miss T > 120 GeV, and this requirement heavily reduces the amount of Z + jets background. In signal events with no initial-or final-state radiation the visible Z boson's transverse momentum is expected to be opposite the missing transverse momentum, and this characteristic is used to further suppress the Z + jets background. The azimuthal angle between the dilepton system and the missing transverse momentum (∆φ(ℓℓ, ì E miss T )) is thus required to be greater than 2.7 and the fractional p T difference, defined as |p miss,jet T − p ℓℓ T |/p ℓℓ T , to be less than 20%, where p Additional selection criteria are applied to keep only events with E miss T originating from neutrinos rather than detector inefficiencies, poorly reconstructed high-p T muons or mismeasurements in the hadronic calorimeter. If at least one reconstructed jet has a p T greater than 100 GeV, the azimuthal angle between the highest-p T jet and the missing transverse momentum is required to be greater than 0.4. Similarly, if E miss T is found to be less than 40% of the scalar sum of the transverse momenta of leptons and jets in the event (H T ), the event is rejected. Finally, to reduce the tt background, events are rejected whenever a b-tagged jet is found.
The sensitivity of the analysis to the VBF production mode is increased by creating a dedicated category of VBF-enriched events. The selection criteria, determined by optimising the expected signal significance using signal and background MC samples, require the presence of at least two jets with p T > 30 GeV where the two highest-p T jets are widely separated in η, |∆η jj | > 4.4, and have an invariant mass m jj greater than 550 GeV.
The signal acceptance, defined as the ratio of the number of reconstructed events passing the analysis requirements to the number of simulated events in each category, is shown in Table 2, for the ggF and VBF Table 2: Signal acceptance for the ℓ + ℓ − νν analysis, for both the ggF and VBF production modes and resonance masses of 300 and 600 GeV. The acceptance is defined as the ratio of the number of reconstructed events after all selection requirements to the number of simulated events for each channel/category.

Mass
Production production modes as well as for different resonance masses. The acceptance increases with mass due to a kinematic threshold determined by the E miss T selection criteria. Hence the ℓ + ℓ − νν search considers only masses of 300 GeV and above, where its inclusion improves the combined sensitivity.

Background estimation
The dominant and irreducible background for this search is non-resonant Z Z production, which accounts for about 60% of the expected background events. The second largest background comes from W Z production (∼30%) followed by Z + jets production with poorly reconstructed E miss T (∼6%). Other sources of background are the WW, tt, Wt and Z → ττ processes (∼3%). Finally, a small contribution comes from W + jets, tt, single-top-quark and multi-jet processes, with at least one jet misidentified as an electron or muon, as well as from ttV/VVV events. In both the ggF-and in the VBF-enriched signal regions, the Z Z background is modelled using MC simulation and normalised using SM predictions, as explained in Section 3. The remaining backgrounds are mostly estimated using control samples in data.
The W Z background is modelled using simulation but a correction factor for its normalisation is extracted as the ratio of data to simulated events in a dedicated control region, after subtracting from data the non-W Z background contributions. The W Z-enriched control sample, called the 3ℓ control region, is built by selecting Z → ℓℓ candidates with an additional electron or muon. This additional lepton is required to satisfy all selection criteria used for the other two leptons, with the only difference that its transverse momentum is required to be greater than 7 GeV. The contamination from Z + jets and tt events is reduced by vetoing events with at least one b-tagged jet and by requiring the transverse mass of the W boson (m W T ), built using the additional lepton and the E miss T vector, to be greater than 60 GeV. The distribution of the missing transverse momentum for data and simulated events in the 3ℓ control region is shown in Figure 3(a). The correction factor derived in the 3ℓ control region is found to be 1.29 ± 0.09, where the uncertainty includes effects from the number of events in the control region as well as from experimental systematic uncertainties. Since there are few events after applying all the VBF selection requirements to the W Z-enriched control sample, the estimation for the VBF-enriched category is performed by including in the 3ℓ control region only the requirement of at least two jets with p T > 30 GeV. Finally, a transfer factor is derived from MC simulation by calculating the probability of events satisfying all analysis selection criteria and containing two jets with p T > 30 GeV to satisfy the |∆η jj | > 4.4 and m jj > 550 GeV requirements. An additional systematic uncertainty obtained from the comparison of the |∆η jj | distribution between S and P -B generators is included to cover potential mismodellings of the VBF selection. Such systematic uncertainty is included in all background estimations when extrapolating from a control region.
The non-resonant background includes mainly WW, tt and Wt processes, but also Z → ττ events in which the τ leptons produce light leptons and E miss T . It is estimated by using a control sample of events with lepton pairs of different flavour (e ± µ ∓ ), satisfying all analysis selection criteria. Figure 3(b) shows the missing-transverse-momentum distribution for e ± µ ∓ events in data and simulation after applying the dilepton invariant-mass selection but before applying the other selection requirements. The non-resonant background in the e + e − and µ + µ − channels is estimated by applying a scale factor ( f ) to the selected events in the e ± µ ∓ control region, such that: where N bkg ee and N bkg µµ are the numbers of electron-and muon-pair events estimated in the signal region and N data,sub eµ is the number of events in the e ± µ ∓ control sample with Z Z, W Z and other small backgrounds subtracted using simulation. The factor f takes into account the different selection efficiencies of e + e − and µ + µ − pairs at the level of the Z → ℓℓ selection, and is measured from data as f 2 = N data ee /N data µµ , where N data ee and N data µµ are the numbers of events passing the Z boson mass requirement (76 < m ℓℓ < 106 GeV) in the electron and muon channel respectively. As no events survive in the e ± µ ∓ control region after applying the full VBF selection, the background estimation is performed by including only the requirement of at least two jets with p T > 30 GeV. The efficiency of the remaining selection requirements on |∆η jj | and m jj is obtained from simulated events.
The number of Z + jets background events in the signal region is estimated from data, using a socalled ABCD method [78], since events with no genuine E miss T in the final state are difficult to model using simulation. The method combines the selection requirements presented in Section 6.1 (with n b-tags representing the number of b-tagged jets in the event) into two Boolean discriminants, V 1 and V 2 , defined as: V 1 ≡ E miss T > 120 GeV and E miss T /H T > 0.4, V 2 ≡ |p miss,jet T − p ℓℓ T |/p ℓℓ T < 0.2 and ∆φ(ℓℓ, ì E miss T ) > 2.7 and ∆R ℓℓ < 1.8 and n b-tags = 0, with all events required to pass the trigger and dilepton invariant-mass selections. The signal region (A) is thus obtained by requiring both V 1 and V 2 to be true, control regions B and C require only one of the two Boolean discriminants to be false (V 1 and V 2 respectively) and finally control region D is defined by requiring both V 1 and V 2 to be false. With this definition, an estimate of the number of events in region A is given by N est A = N obs C × (N obs B /N obs D ), where N obs X is the number of events observed in region X after subtracting non-Z-boson backgrounds. This relation holds as long as the correlation between V 1 and V 2 is small, and this is achieved by introducing two additional requirements on control regions B and D, namely E miss T > 30 GeV and E miss T / H T > 0.1. The estimation of the Z + jets background was cross-checked with another approach in which a control region is defined by inverting the analysis selection on E miss T /H T and then using Z + jets MC simulation to perform the extrapolation to the signal region, yielding results compatible with the ABCD method. Finally, the estimate for the VBF-enriched category is performed by extrapolating the inclusive result obtained with the ABCD method to the VBF signal region, extracting the efficiency of the two-jet, |∆η jj | and m jj selection criteria from Z + jets simulation.    Figure 3: Missing transverse momentum E miss T distribution (a) for events in the 3ℓ control region as defined in the text and (b) for e ± µ ∓ lepton pairs after applying the dilepton invariant mass requirement, before applying the rest of the control region selection. The backgrounds are determined following the description in Section 6.2 and the last bin includes the overflow. The small excess below 120 GeV in (b) arises from Z + jets background which is here taken from simulation, and lies outside the control region. The error bars on the data points indicate the statistical uncertainty, while the systematic uncertainty in the prediction is shown by the hatched band. The lower panels show the ratio of data to prediction.
The W + jets and multi-jet background contributions are estimated from data using a so-called fake-factor method [79]. A control region enriched in fake leptons or non-prompt leptons from decays of hadrons is designed by requiring one lepton to pass all analysis requirements (baseline selection) and the other one to not pass either the lepton "medium" identification or the isolation criteria (inverted selection). The background in the signal region is then derived using a transfer factor, measured in a data sample enriched in Z + jets events, as the ratio of jets passing the baseline selection to those passing the inverted selection.
Finally, the background from the ttV and VVV processes is estimated using MC simulation.

Signal and background modelling
The modelling of the transverse mass m T distribution for signal and background is based on templates derived from fully-simulated events and afterwards used to fit the data. In the case of a narrow resonance, simulated MC events generated for fixed mass hypotheses as described in Section 3 are used as the inputs in the moment-morphing technique [80] to obtain the m T distribution for any other mass hypothesis.
The extraction of the interference terms for the LWA case is performed in the same way as in the ℓ + ℓ − ℓ + ℓ − final state, as described in Section 5.3. In the case of the ℓ + ℓ − νν final state a correction factor, extracted as a function of m Z Z , is used to reweight the interference distributions obtained at particle level to account for reconstruction effects. The final expected LWA m T distribution is obtained from the combination of the interference distributions with simulated m T distributions, which are interpolated between the simulated mass points with a weighting technique using the Higgs propagator, a method similar to that used for the interference.

Systematic uncertainties
The systematic uncertainties can be classified into experimental and theoretical uncertainties. The first category relates to the reconstruction and identification of leptons and jets, their energy scale and resolution, and the integrated luminosity. Systematic uncertainties in the data-driven background estimates are also included in this category. The second category includes uncertainties in the theoretical description of the signal and background processes.
In both cases the uncertainties are implemented as additional nuisance parameters (NP) that are constrained by a Gaussian distribution in the profile likelihood ratio, as discussed in Section 8.1. The uncertainties affect the signal acceptance, its selection efficiency and the discriminant distributions as well as the background estimates for both final states. Each source of uncertainty is either fully correlated or anticorrelated among the different channels and categories.

Experimental uncertainties
The uncertainty in the combined 2015 and 2016 integrated luminosity is 3.2%. This is derived from a preliminary calibration of the luminosity scale using x-y beam-separation scans performed in August 2015 and May 2016, following a methodology similar to that detailed in Ref. [81].
The lepton identification and reconstruction efficiency and energy/momentum scale and resolution are derived from data using large samples of J/ψ → ℓℓ and Z → ℓℓ decays. The uncertainties in the reconstruction performance are computed following the method described in Ref.
[62] for electrons. Typical uncertainties in the identification and reconstruction efficiency are in the range 0.5%-3.0% for muons and 1.0%-1.7% for electrons. The uncertainties in the electron energy scale, the muon momentum scale and their resolutions are small, and are fully correlated between the two searches (ℓ + ℓ − ℓ + ℓ − and ℓ + ℓ − νν final states).
The uncertainties in the jet energy scale and resolution have several sources, including uncertainties in the absolute and relative in situ calibration, the correction for pile-up, the flavour composition and response [66]. These uncertainties are separated into independent components, which are fully correlated between the two searches. They vary from 4.5% for jets with transverse momentum p T = 20 GeV, decreasing to 1% for jets with p T = 100-1500 GeV and increasing again to 3% for jets with higher p T , for the average pile-up conditions of the 2015 and 2016 data-taking period.
Uncertainties in the lepton and jet energy scales are propagated to the uncertainty in the E miss T . Additionally, the uncertainties from the momentum scale and resolution of the tracks that are not associated with any identified lepton or jet contribute 8% and 3% respectively, to the uncertainty in the E miss T value.
The efficiency of the lepton triggers in events with reconstructed leptons is nearly 100%, and hence the related uncertainties are negligible.

Theoretical uncertainties
For simulated signal and backgrounds, theoretical modelling uncertainties associated with the PDFs, missing QCD higher-order corrections (via variations of factorisation and renormalisation scales), and parton showering are considered.
For all signal hypotheses under consideration, the largest theoretical modelling uncertainties are due to missing QCD higher-order corrections and parton showering. The missing QCD higher-order corrections for ggF production events that fall into the VBF-enriched category are accounted for by varying the scales in M G 5_aMC@NLO and affect the signal acceptance by 10%. Parton showering uncertainties are of order 10% and are estimated by comparing P 8.212 to H ++ [82].
For the qq → Z Z background, the effect of the PDF uncertainties in the full mass range varies between 2% and 5% in all categories, and that of missing QCD higher-order corrections is about 10% in the ggFenriched categories and 30% in the VBF-enriched category. The parton-shower uncertainties result in less than 1% impact in the ggF-enriched categories and about 10% impact in the VBF-enriched category.
For the gg → Z Z background, as described in Section 3, a 60% relative uncertainty in the inclusive cross section is considered, while a 100% uncertainty is assigned in the VBF-enriched category.

Statistical procedure
The statistical treatment of the data follows the procedure for the Higgs-boson search combination [83,84], and is implemented with RooFit [85] and RooStats [86]. The test statistic employed for hypothesis testing and limit setting is the profiled likelihood ratio Λ(α, θ), which depends on one or more parameters of interest α, and additional nuisance parameters θ. The parameter of interest is the cross section times branching ratio for heavy-resonance production, assumed to be correlated between the two searches. The nuisance parameters represent the estimates of the systematic uncertainties and are each constrained by a Gaussian distribution. For each category of each search, a likelihood fit to the kinematic distribution of a discriminating variable is used to further separate signal from background. The ℓ + ℓ − ℓ + ℓ − final state uses m 4ℓ as the discriminant in each category, while the ℓ + ℓ − νν final state uses m T in each category except for the VBF-enriched one where only the overall event counts are used.
As discussed in Section 7, the signal acceptance uncertainties, and many of the background theoretical and experimental uncertainties, are treated as fully correlated between the searches. A given correlated uncertainty is modelled in the fit by using a nuisance parameter common to all of the searches. The impact of a systematic uncertainty on the result depends on the production mode and the mass hypothesis. For ggF production, at lower masses the luminosity uncertainty, the modelling uncertainty of the Z + jets background and the statistical uncertainty in the eµ control region of the ℓ + ℓ − νν final state dominate, and at higher masses the uncertainties in the electron-isolation efficiency become important, as also seen in VBF production. For VBF production, the dominant uncertainties come from the theoretical predictions of the Z Z events in the VBF category. Additionally at lower masses, the pile-up reweighting and the jet-energy-resolution uncertainties are also important. Table 3 shows the impact of the leading systematic uncertainties on the predicted signal event yield when the cross section times branching ratio is set to the expected upper limit (shown in Figure 6), for ggF and VBF production modes. The impact of the uncertainty in the integrated luminosity, 3.2%, enters both in the normalisation of the fitted number of signal events as well as in the background predicted by simulation. This leads to a luminosity uncertainty which varies from 4% to 7% across the mass distribution, depending on the signal-to-background ratio.

General results
The numbers of observed candidate events with mass above 130 GeV together with the expected background yields are presented in Table 4 for each of the four categories of the ℓ + ℓ − ℓ + ℓ − analysis. The m 4ℓ spectrum for the ggF-enriched and VBF-enriched categories is shown in Figure 4. Table 5 contains the number of observed candidate events along with the background yields for the ℓ + ℓ − νν analysis, while Figure 5 shows the m T distribution for the electron and muon channels with the ggF-enriched and VBF-enriched categories combined.
In the ℓ + ℓ − ℓ + ℓ − search, two excesses are observed in the data for m 4ℓ around 240 and 700 GeV, each with a local significance of 3.6σ estimated in the asymptotic approximation, assuming the signal comes only from ggF production. The global significance is 2.2σ and is calculated, for each excess individually, using the NWA, in the range of 200 GeV< m H < 1200 GeV using pseudo-experiments.
The excess at 240 GeV is observed mostly in the 4e channel, while the one at 700 GeV is observed in all channels and categories. No significant deviation from the expected background is observed in the ℓ + ℓ − νν final state. The excess observed in the ℓ + ℓ − ℓ + ℓ − search at a mass around 700 GeV is excluded at 95% confidence level (CL) by the ℓ + ℓ − νν search, which is more sensitive in this mass range. The excess at 240 GeV is not covered by the ℓ + ℓ − νν search, the sensitivity of which starts from Z Z 297 ± 1 ± 40 480 ± 1 ± 60 193 ± 1 ± 25 15 ± 0.1 ± 6.0 Z Z (EW) 1.92 ± 0.11 ± 0.19 3.36 ± 0.14 ± 0.33 1.88 ± 0.12 ± 0.20 3.0 ± 0.1 ± 2.2 Z + jets/tt/W Z 3.7 ± 0.1 ± 0.8 7.8 ± 0.1 ± 1.1 4.4 ± 0.1 ± 0.8 0.37 ± 0.01 ± 0.05 Other backgrounds 5.1 ± 0.1 ± 0.6 8.7 ± 0.1 ± 1.0 4.0 ± 0.1 ± 0.5 0.80 ± 0.02 ± 0.30 Total background 308 ± 1 ± 40 500 ± 1 ± 60 203 ± 1 ± 25 19.5 ± 0.2 ± 8.0    Figure 4: Distribution of the four-lepton invariant mass m 4ℓ in the ℓ + ℓ − ℓ + ℓ − search for (a) the ggF-enriched category and (b) the VBF-enriched category. The backgrounds are determined following the description in Section 5.2 and the last bin includes the overflow. The simulated m H = 600 GeV signal is normalized to a cross section corresponding to five times the observed limit given in Section 8.3.1. The error bars on the data points indicate the statistical uncertainty, while the systematic uncertainty in the prediction is shown by the hatched band. The lower panels show the ratio of data to prediction.
93 ± 2 ± 4 99.5 ± 2.3 ± 3.2 1.29 ± 0.04 ± 0.27 WW/tt/Wt/Z → ττ 9.2 ± 2.2 ± 1.4 10.7 ± 2.5 ± 0.9 0.39 ± 0.24 ± 0.26 Z + jets 17 ± 1 ± 11 19 ± 1 ± 17 0.8 ± 0.1 ± 0.5 Other backgrounds 1.12 ± 0.04 ± 0.08 1.03 ± 0.04 ± 0.08 0.03 ± 0.01 ± 0.01 Total background 297 ± 4 ± 24 311 ± 5 ± 27 4.6 ± 0.4 ± 0.9 Observed 320 352 9 Events / 50 GeV Wt , WW Other backgrounds Uncertainty  Figure 5: Transverse mass m T distribution in the ℓ + ℓ − νν search for (a) the electron channel and (b) the muon channel, including events from both the ggF-enriched and the VBF-enriched categories. The backgrounds are determined following the description in Section 6.2 and the last bin includes the overflow. The simulated m H = 600 GeV signal is normalized to a cross section corresponding to five times the observed limit given in Section 8.3.1. The error bars on the data points indicate the statistical uncertainty and markers are drawn at the bin centre. The systematic uncertainty in the prediction is shown by the hatched band. The lower panels show the ratio of data to prediction. 300 GeV. When combining the results from the two final states, the largest deviation with respect to the background expectation is observed around 700 GeV with a global significance of less than 1σ and a local significance of about 2σ. The combined yield of the two final states is 1870 events observed in data compared to 1643 ± 164 (combined statistical and systematic uncertainty) for the expected background. This corresponds to a 1.3σ global excess in data. Since no significant excess is found, the results are interpreted as upper limits on the production cross section of a spin-0 or spin-2 resonance.

Spin-0 resonance interpretation
Limits from the combination of the two searches in the context of a spin-0 resonance are described below.

NWA interpretation
Upper limits on the cross section times branching ratio (σ × B(H → Z Z )) for a heavy resonance are obtained as a function of m H with the CL s procedure [87] in the asymptotic approximation from the combination of the two final states. It is assumed that an additional heavy scalar would be produced predominantly via the ggF and VBF processes but that the ratio of the two production mechanisms is unknown in the absence of a specific model. For this reason, fits for the ggF and VBF production processes are done separately, and in each case the other process is allowed to float in the fit as an additional nuisance parameter. Figure 6 presents the observed and expected limits at 95% CL on σ × B(H → Z Z ) of a narrow scalar resonance for the ggF (left) and VBF (right) production modes, as well as the expected limits from the ℓ + ℓ − ℓ + ℓ − and ℓ + ℓ − νν searches. This result is valid for models in which the width is less than 0.5% of m H . When combining the two final states, the 95% CL upper limits range from 0.68 pb at m H = 242 GeV to 11 fb at m H = 1200 GeV for the ggF production mode and from 0.41 pb at m H = 236 GeV to 13 fb at m H = 1200 GeV for the vector-boson fusion production mode. Compared with the results from Run 1 [21], where all four final states of Z Z decays were combined, the exclusion region presented here is significantly extended considering that the ratios of parton luminosities [88] increase by factors of about two to seven for heavy scalar masses from 200 GeV to 1200 GeV.

LWA interpretation
In the case of the LWA, limits on the cross section for the ggF production mode times branching ratio (σ ggF × B(H → Z Z )) are set for different widths of the heavy scalar. The interference between the heavy scalar and the SM Higgs boson, H-h, as well as the heavy scalar and the gg → Z Z continuum, H-B, are modelled by either analytical functions or reweighting the signal-only events as explained in Sections 5.3 and 6.3. Figures 7(a), 7(b), and 7(c) show the limits for a width of 1%, 5% and 10% of m H respectively. The limits are set for masses of m H higher than 400 GeV.

2HDM interpretation
A search in the context of a CP-conserving 2HDM is also presented. This model has five physical Higgs bosons after electroweak symmetry breaking: two CP-even, one CP-odd, and two charged. The model considered here has seven free parameters: the Higgs boson masses, the ratio of the vacuum expectation values of the two doublets (tan β), the mixing angle between the CP-even Higgs bosons (α), and the potential parameter m 2 12 that mixes the two Higgs doublets. The two Higgs doublets Φ 1 and Φ 2 can couple to leptons and up-and down-type quarks in several ways. In the Type-I model, Φ 2 couples to all quarks and leptons, whereas for Type-II, Φ 1 couples to down-type quarks and leptons and Φ 2 couples to up-type quarks. The "lepton-specific" model is similar to Type-I except for the fact that the leptons couple to Φ 1 , instead of Φ 2 ; the "flipped" model is similar to Type-II except that the leptons couple to Φ 2 , instead of Φ 1 . In all these models, the coupling of the heaviest CP-even Higgs boson to vector bosons is proportional to cos(β − α). In the limit cos(β − α) → 0, the light CP-even Higgs boson is indistinguishable from a SM Higgs boson with the same mass. In the context of H → Z Z decays there is no direct coupling of the Higgs boson to leptons, and so only the Type-I and -II interpretations are presented. Figure 8 shows exclusion limits in the tan β versus cos(β − α) plane for Type-I and Type-II 2HDMs, for a heavy Higgs boson with mass m H = 200 GeV. This m H value is chosen so that the assumption of a narrow Higgs boson is valid over most of the parameter space, and the experimental sensitivity is maximal. At this low mass, only the ℓ + ℓ − ℓ + ℓ − final state contributes to this result. The range of cos(β − α) and tan β explored is limited to the region where the assumption of a heavy narrow Higgs boson with negligible interference is valid. When calculating the limits at a given choice of cos(β − α) and tan β, the relative rates of ggF and VBF production in the fit are set to the prediction of the 2HDM for that parameter choice. Figure 9 shows exclusion limits as a function of the heavy Higgs boson mass m H and the parameter tan β for cos(β − α) = −0.1. The white regions in the exclusion plots indicate regions of parameter space which are not excluded by the present analysis. In these regions the cross section predicted by the 2HDM is below the observed cross section limit. Compared with the results from Run 1 [21], the exclusion presented here is almost twice as stringent.

Spin-2 resonance interpretation
The results are also interpreted as a search for a Kaluza-Klein graviton excitation, G KK , in the context of the bulk RS model using the ℓ + ℓ − νν final state because the ℓ + ℓ − ℓ + ℓ − final state was found to have negligible sensitivity for this type of model. The limits on σ × B(G KK → Z Z) at 95% CL as a function of the KK graviton mass, m(G KK ), are shown in Figure 10 together with the predicted G KK cross section. A spin-2 graviton is excluded up to a mass of 1300 GeV. These limits have been extracted using the asymptotic approximation, and they were verified to be correct within about 4% using pseudo-experiments.    Figure 10: The upper limits at 95% CL on cross section times branching ratio σ × B(G KK → Z Z) for a KK graviton produced with k/M Pl = 1. The green and yellow bands give the ±1σ and ±2σ uncertainties in the expected limits. The predicted production cross section times branching ratio as a function of the G KK mass m(G KK ) is shown by the red solid line.

Summary
A search is conducted for heavy resonances decaying into a pair of Z bosons which subsequently decay into ℓ + ℓ − ℓ + ℓ − or ℓ + ℓ − νν final states. The search uses proton-proton collision data collected with the ATLAS detector during 2015 and 2016 at the Large Hadron Collider at a centre-of-mass energy of 13 TeV corresponding to an integrated luminosity of 36.1 fb −1 . The results of the search are interpreted as upper limits on the production cross section of a spin-0 or spin-2 resonance. The mass range of the hypothetical resonances considered is between 200 GeV and 2000 GeV depending on the final state and the model considered. The spin-0 resonance is assumed to be a heavy scalar, whose dominant production modes are gluon-gluon fusion and vector-boson fusion and it is studied in the narrow-width approximation and with the large-width assumption. In the case of the narrow-width approximation, limits on the production rate of a heavy scalar decaying into two Z bosons are set separately for ggF and VBF production modes. Combining the two final states, 95% CL upper limits range from 0.68 pb at m H = 242 GeV to 11 fb at m H = 1200 GeV for the gluon-gluon fusion production mode and from 0.41 pb at m H = 236 GeV to 13 fb at m H = 1200 GeV for the vector-boson fusion production mode. The results are also interpreted in the context of Type-I and Type-II two-Higgs-doublet models, with exclusion contours given in the tan β versus cos(β − α) (for m H = 200 GeV) and tan β versus m H planes. This m H value is chosen so that the assumption of a narrow Higgs boson is valid over most of the parameter space and the experimental sensitivity is maximal. The limits on the production rate of a large-width scalar are obtained for widths of 1%, 5% and 10% of the mass of the resonance, with the interference between the heavy scalar and the SM Higgs boson as well as the heavy scalar and the gg → Z Z continuum taken into account. In the framework of the Randall-Sundrum model with one warped extra dimension a graviton excitation spin-2 resonance with m(G KK ) < 1300 GeV is excluded at 95% CL.