Search for heavy resonances decaying to ZZ or ZW and axion-like particles mediating nonresonant ZZ or ZH production at $\sqrt{s} =$ 13 TeV

A search has been performed for heavy resonances decaying to ZZ or ZW and for axion-like particles (ALPs) mediating nonresonant ZZ or ZH production, in final states with two charged leptons ($\ell$ = e, $\mu$) produced by the decay of a Z boson, and two quarks produced by the decay of a Z, W, or Higgs boson H. The analysis is sensitive to resonances with masses in the range 450 to 2000 GeV. Two categories are defined corresponding to the merged or resolved reconstruction of the hadronically decaying boson. The search is based on data collected during 2016-2018 by the CMS experiment at the LHC in proton-proton collisions at a center-of-mass energy of 13 TeV, corresponding to an integrated luminosity of 138 fb$^{-1}$. No significant excess is observed in the data above the standard model background expectation. Upper limits on the production cross section of heavy, narrow spin-2 and spin-1 resonances are derived as functions of the resonance mass, and exclusion limits on the production of bulk graviton particles and W$'$ bosons are calculated in the framework of the warped extra dimensions and heavy vector triplet models, respectively. In addition, upper limits on the ALP-mediated diboson production cross section and ALP couplings to standard model particles are obtained in the framework of linear and chiral effective field theories. These are the first limits on nonresonant ALP-mediated ZZ and ZH production obtained by the LHC experiments.


Introduction
Several extensions of the standard model (SM) predict the existence of new particles with large couplings to a pair of bosons, either Higgs bosons H or electroweak vector bosons V = W, Z. One possibility is the Randall-Sundrum model with warped extra spatial dimensions (WED) [1][2][3], which predicts heavy spin-2 gravitons (G). Another example is the spin-1 heavy vector triplet (HVT) model [4], which postulates the existence of Z and W bosons.
Diboson signatures can also be used to search for axion-like particles (ALPs) [5,6], designated in this paper by "a". These are neutral pseudoscalar bosons with derivative interactions with SM particles, and often appear in extensions of the SM. Colliders allow for searches in a wide range of ALP masses and couplings [7][8][9]. An interesting diboson production mechanism is nonresonant ALP-mediated scattering; the ALP is an off-shell mediator in the s-channel 2 → 2 scattering processes [9]. Figure 1 shows the Feynman diagrams for the processes gg → a * → ZZ (left) and gg → a * → ZH (right), where a * is an off-shell ALP. These channels are sensitive to the product of the ALP coupling to gluons and the relevant couplings to dibosons. These channels have been identified as promising because the cross sections are sufficiently large to enable strong constraints to be placed on the theoretical models using the 2016-2018 LHC data set. The derivative interactions enhance the cross section at large diboson invariant masses, √ŝ . In the ALP off-shell regime √ŝ m a , m Z , where m a is the ALP mass, the partonic cross section for gg → a * → ZZ varies asŝ/ f 4 a , where f a is the energy scale characterizing the new physics. The same type of energy behavior holds for gg → a * → ZH. Such energy dependence is only valid as long as the energies probed in the scattering process are smaller than f a . Factoring in the parton distribution functions (PDFs) of the proton, which tame the energy growth, the differential cross sections for the ALP-mediated processes pp → ZZ, ZH diminish more slowly as a function of √ŝ than for the SM backgrounds. The nonresonant cross sections and kinematical distributions are found to be independent of the ALP mass from arbitrarily light masses up to masses of the order of 100 GeV [9]. The strategy of the search is to look for deviations with respect to the SM expectations in the tails of the experimental diboson mass distributions. This paper reports on the results of a search for spin-2 gravitons decaying to ZZ, spin-1 W bosons decaying to ZW, and nonresonant ALP-mediated ZZ or ZH production, in final states with two charged leptons ( = e, µ) produced by the decay of a Z boson, and two quarks produced by the decay of a Z, W, or Higgs boson. These two quarks can be reconstructed as one single "merged" jet in hadronic decays of a Lorentz-boosted boson, or as two individually resolved jets. The analysis is sensitive to resonances with masses in the range from 450 to 2000 GeV. The width of the heavy resonance is taken to be small in comparison to the experimental resolution [10]. This is the first time a search for the effects of nonresonant ALP-mediated ZZ or ZH production has been performed at the CERN LHC. The search analyzes data collected from 2016 to 2018 by the CMS experiment at √ s = 13 TeV.
Previous searches for heavy resonances decaying to two SM bosons at the LHC [11-21] have found no evidence of a signal, providing stringent upper limits on the corresponding production cross sections. Limits on ALP couplings derived from reinterpretations of ATLAS measurements can be found in Ref. [22].
This paper is organized as follows: in Section 2, a description of the signal models, simulated background, and data samples used in the analysis is provided; Section 3 briefly describes the CMS detector; Section 4 provides a description of the event reconstruction; in Section 5, the event selection is discussed; Section 6 describes the estimation of the SM background; the systematic uncertainties affecting the analysis are presented in Section 7; and the results of the search for heavy spin-2 and spin-1 resonances and ALP-mediated diboson production are presented in Section 8. Finally, the results are summarized in Section 9.
Tabulated results are provided in the HEPData record for this analysis [23].

Data and simulated samples
This analysis uses data collected with the CMS detector during proton-proton (pp) collisions at the LHC at √ s = 13 TeV, corresponding to an integrated luminosity of 138 fb −1 . The events were selected online by criteria that require the presence of at least one electron or muon; these criteria are described in Section 5.
Simulated signal samples are used in the analysis to optimize the search for production of heavy spin-2 and spin-1 resonances, and ALP-mediated ZZ and ZH production. Spin-2 and spin-1 resonant signal samples are generated according to the WED and HVT scenarios, respectively.
In the WED scenario [24,25], the main free parameters are the mass of the first Kaluza-Klein graviton excitation (the bulk graviton mass) and the ratio κ = κ/M Pl , where κ is a curvature parameter of the WED metric and M Pl is the reduced Planck mass. The coupling of the bulk graviton to vector bosons is determined by κ. Bulk graviton production cross sections, widths, and branching fractions are taken from Ref. [26].
The HVT model [4] introduces a spin-1 triplet of Z and W bosons. The HVT model generalizes a large number of explicit cases in terms of three parameters: c H determines the coupling to SM vector and Higgs bosons; c F determines the coupling to fermions; and g V is the overall strength of the new vector boson triplet interactions. Benchmark "Model A", with g V = 1, is representative of a model of weakly coupled vector resonances in an extension of the SM gauge group where the HVT bosons have comparable decay branching fractions into SM fermions and vector bosons. Benchmark "Model B", with g V = 3, is representative of a composite model scenario where the HVT boson couplings to fermions are suppressed. The W boson production cross sections, widths, and branching fractions are taken from Ref. [4].
For both scenarios, the samples are generated at leading order (LO) in quantum chromodynamics (QCD) with the MADGRAPH5 aMC@NLO generator [27] versions 2.2.2 (2016) and 2.6.5 (2017 and 2018). Different resonance mass hypotheses are considered in the range 450 to 2000 GeV. Since the resonance width is small in comparison with the experimental resolution, for simplicity, the width is taken to be 1 MeV in the simulation (the narrow-width approximation). The generated spin-2 bulk graviton is forced to decay into two Z bosons, one decaying leptonically into any pair of charged leptons, and the other decaying hadronically into a pair of quarks. In the case of the spin-1 W boson, the resonance is forced to decay into one Z and one W boson; additionally, the Z boson is then forced to decay to a pair of electrons, muons, or tau leptons, while the W boson is forced to decay into a pair of quarks.
The ALP interactions can be parameterized using the model-independent approach of effective field theories (EFTs). There are two EFT implementations: linear (related to weakly coupled new physics models, with a minimal set of parameters) [28,29] and chiral (related to strongly coupled new physics models, with more parameters) [30][31][32][33]. The ALP couplings to SM particles are proportional to a coefficient c and inversely proportional to the new physics energy scale f a . Here, we adopt the definitions of Refs. [7,9]. Classical searches for ALPs at colliders consider their couplings to photons and gluons. The coefficient of the gluon coupling is c G . More recently, interest has extended to consider ALP couplings to electroweak bosons: ZZ, WW, and Zγ [7]. At LO in the linear model, all these and the coupling to photons are related by gauge symmetry to two basic electroweak couplings: c W and c B . The coupling of the ALP to ZZ is where c w and s w denote the cosine and sine of the Weinberg angle, respectively. In this analysis, small couplings of the ALP to fermions are set to zero. In the chiral model, the ALP still couples to gluons with coupling c G . Additionally, the ALP couples to ZH, and there are several coefficients contributing to the ZH vertex [7]. For simplicity, in this chiral analysis, the ALP to ZH coupling coefficient a 2D is studied; small couplings of the ALP to fermions are neglected and the ZZ channel is considered closed. In both the linear and chiral implementations, interference between the SM and ALP-mediated processes is neglected.
Simulated signals of ALP-mediated production are generated with MADGRAPH5 aMC@NLO 2.6.5 using both the linear and chiral EFT models [7]. In the linear model, gg → a * → ZZ events are generated at LO with c G / f a = c Z / f a = 1 TeV −1 and an ALP mass of 1 MeV. The MADGRAPH5 aMC@NLO cross sections are 77 (92, 94, 95) pb for √ŝ < f a = 1 (2, 3, 7) TeV, respectively. One Z boson is forced to decay leptonically and the other hadronically. In the chiral model, gg → a * → ZH events are generated at LO with c G / f a = a 2D / f a = 1 TeV −1 and an ALP mass of 1 MeV; all other chiral coefficients are set to zero. The MADGRAPH5 aMC@NLO cross sections are 63 (78, 80, 81) pb for √ŝ < f a = 1 (2, 3, 7) TeV, respectively. The Z boson is forced to decay leptonically; the Higgs boson is decayed according to the branching fractions predicted in the SM. In both cases, there is a 7% systematic uncertainty in the size of the cross sections related to the renormalization and factorization scales, and a 1.2% systematic uncertainty related to the PDFs.
Several SM processes yielding final states with charged leptons and jets are sources of background events for the analysis, and corresponding Monte Carlo simulated samples have been generated to study them.
The SM production of a Z boson in association with quarks or gluons in the final state (Z + jets) represents the dominant background process for the analysis, having topological similarities to the signal because of the presence of a pair of charged leptons and jets. However, since the quark-and gluon-induced jets are not associated with the decay of a vector boson, the mergedjet and dijet mass spectra are characterized by smooth distributions, and the distribution of the 2 2q system invariant mass falls exponentially. In contrast, the signals are expected to have peaking distributions in the merged-jet and dijet mass spectra, and a distinctive 2 2q mass distribution. The simulated Z + jets samples are produced with MADGRAPH5 aMC@NLO at next to LO (NLO) for 2016 samples, and LO for 2017 and 2018 samples. Both LO and NLO sample events are reweighted to account for electroweak NLO corrections, with the LO events also reweighted for QCD NLO corrections [34].
The SM diboson production of ZV and ZH is an irreducible source of background for the analysis, since the merged-jet and dijet mass spectra will contain a peak from the hadronic decay of H, W, and Z bosons; however, this process produces a smoothly falling 2 2q invariant mass distribution. The SM production of pairs of bosons (ZZ, ZW, and ZH) is simulated at NLO with MADGRAPH5 aMC@NLO.
The parton showering and hadronization is simulated by interfacing the event generators with PYTHIA 8.226 [35] [39] set for 2017 and 2018 resonant samples and backgrounds, as well as all ALP simulated samples. Additional pp interactions occurring in the same or nearby bunch crossings (pileup) are added to the event simulation, with a frequency distribution adjusted to match that observed in data. All samples are processed through a simulation of the CMS detector using GEANT4 [40] and reconstructed using the same algorithms as those for analysis of the data.

The CMS detector
The central feature of the CMS apparatus is a superconducting solenoid of 6 m internal diameter, providing a magnetic field of 3.8 T. Within the solenoid volume are a silicon pixel and strip tracker, a lead tungstate crystal electromagnetic calorimeter (ECAL), and a brass and scintillator hadron calorimeter (HCAL), each composed of a barrel and two endcap sections. Forward calorimeters extend the pseudorapidity η coverage provided by the barrel and endcap detectors. Muons are measured in gaseous detectors embedded in the steel flux-return yoke outside the solenoid.
Events of interest are selected using a two-tiered trigger system. The first level, composed of custom hardware processors, uses information from the calorimeters and muon detectors to select events at a rate of around 100 kHz within a fixed latency of about 4 µs [41]. The second level, known as the high-level trigger, consists of a farm of processors running a version of the full event reconstruction software optimized for fast processing, and reduces the event rate to around 1 kHz before data storage [42].
A more detailed description of the CMS detector, together with a definition of the coordinate system used and the relevant kinematic variables, can be found in Ref. [43].

Event reconstruction
The particle-flow (PF) algorithm [44] aims to reconstruct and identify each individual particle in an event, with an optimized combination of information from the various elements of the CMS detector. The energy of photons is obtained from the ECAL measurement. The energy of electrons is determined from a combination of the electron momentum at the primary interaction vertex as determined by the tracker, the energy of the corresponding ECAL cluster, and the energy sum of all bremsstrahlung photons spatially compatible with originating from the electron track. The energy of muons is obtained from the curvature of the corresponding track. The energy of charged hadrons is determined from a combination of their momentum measured in the tracker and the matching ECAL and HCAL energy deposits, corrected for the response function of the calorimeters to hadronic showers. Finally, the energy of neutral hadrons is obtained from the corresponding corrected ECAL and HCAL energies.
The candidate vertex with the largest value of summed physics-object p 2 T , where p T is the transverse momentum, is taken to be the primary pp interaction vertex. The physics objects are the jets, clustered using the jet finding algorithm [45,46] with the tracks assigned to candidate vertices as inputs, and the associated missing p T , taken as the negative vector sum of the p T of those jets.
The electron momentum is estimated by combining the energy measurement in the ECAL with the momentum measurement in the tracker. The momentum resolution for electrons with p T ≈ 45 GeV from Z → ee decays ranges from 1.7 to 4.5%. It is generally better in the barrel region than in the endcaps, and also depends on the bremsstrahlung energy emitted by the electron as it traverses the material in front of the ECAL [50].
Muons are measured in the range |η| < 2.4, with detection planes made using three technologies: drift tubes, cathode strip chambers, and resistive-plate chambers. The single-muon trigger efficiency exceeds 90% over the full η range, and the efficiency to reconstruct and identify muons is greater than 96%. Matching muons to tracks measured in the silicon tracker results in a relative p T resolution, for muons with p T up to 100 GeV, of 1% in the barrel and 3% in the endcaps, and of better than 7% for muons in the barrel with p T up to 1 TeV [51].
Both electrons and muons are required to be isolated from hadronic activity in the event. An isolation variable is defined as the scalar sum of the p T of charged hadrons originating from the primary vertex, plus the scalar sums of the p T for neutral hadrons and photons, in a cone of ∆R = √ (∆η) 2 + (∆φ) 2 < 0.3 (0.4) around the electron (muon) direction, where φ is the azimuthal angle, corrected to account for the contribution from neutral candidates originating from pileup [50,51].
Hadronic jets are clustered from particles reconstructed by the PF algorithm using the infraredand collinear-safe anti-k T algorithm [45,46] with distance parameters of 0.4 (AK4 jets) and 0.8 (AK8 jets). The jet momentum is determined as the vectorial sum of all particle momenta in the jet, and is found from simulation to be, on average, within 5-10% of the true momentum over the entire p T spectrum and detector acceptance. For AK4 jets, contamination from pileup is suppressed using charged-hadron subtraction, which removes from the list of PF candidates all charged particles originating from vertices other than the primary interaction vertex of the event. The residual contribution from neutral particles originating from pileup vertices is removed by means of an event-by-event jet-area-based correction to the jet four-momentum [52]. For AK8 jets, the pileup per particle identification algorithm (PUPPI) [53] is used to mitigate the effect of pileup at the reconstructed particle level, making use of local shape information, event pileup properties, and tracking information.
Jet energy corrections are derived from simulation studies so that the average measured energy of jets becomes identical to that of particle-level jets. In situ measurements of the momentum balance in dijet, photon + jet, Z + jet, and multijet events are used to determine any residual differences between the jet energy scale in data and in simulation, and appropriate corrections are made [10]. Additional selection criteria are applied to each jet to remove jets potentially dominated by instrumental effects or reconstruction failures [54].
A jet grooming technique is used for AK8 jets in this analysis to help identify and discriminate between merged jets from hadronic decays of boosted V and Higgs bosons, and jets from quarks and gluons. The AK8 jets are groomed by means of the modified mass drop tagger algorithm [55], also known as the soft-drop algorithm, with angular exponent β = 0, soft cutoff threshold z cut < 0.1, and characteristic radius R 0 = 0.8 [56]. Subjets are obtained by reversing the last step of the jet soft-drop clustering. The soft-drop algorithm does not fully reject contributions from the underlying event and pileup. The mass of the AK8 jet (m J ) is therefore defined as the invariant mass associated to the four-momentum of the soft-drop jet, after the PUPPI corrections are applied.
Discrimination between AK8 jets originating from vector boson decays and those originating from gluons and quarks is further improved by the N-subjettiness jet substructure variable [57]. This observable exploits the distribution of the jet constituents found in the proximity of the subjet axes to determine if the jet can be effectively subdivided into a number N of subjets. The generic N-subjettiness variable τ N is defined as the p T -weighted sum of the angular distance of all the k jet constituents from the closest subjet: The normalization factor d 0 is defined as d 0 = ∑ k p T,k R j , with R j the clustering parameter of the original jet. In this analysis, which aims to select V → qq decays, the variable that best discriminates between V boson jets and those from quarks and gluons is the ratio of the 2subjettiness to 1-subjettiness: τ 21 = τ 2 /τ 1 . The τ 21 observable is calculated for the jet before the grooming procedure, and includes the PUPPI algorithm corrections for pileup mitigation.
For the identification of jets originating from the hadronization of bottom quarks, the deep neural network-based combined secondary vertex (DEEPCSV) algorithm [58] is used, either directly on the AK4 jets or on the AK8 soft-drop subjets with PUPPI pileup mitigation applied. A DEEPCSV loose (medium) b tag corresponds to a ≈84 (≈68)% efficiency for b quark identification and a ≈10 (≈1)% light-flavor or gluon jet mistag rate. The b tagging and mistag jet performances are extracted from simulation and then corrected using scale factors measured in data samples enriched in b quark jets and light-flavor or gluon jets [58].
Only those AK4 and AK8 jets passing the tight jet identification requirements [54], and reconstructed in the central region of the detector (|η| < 2.4), are considered in the analysis.

Event selection
In the online selection, events must satisfy either a single-electron or a single-muon trigger. The electron trigger requires either an electron with p T > 27 (32) GeV and |η| < 2.4 in the 2016 (2017 and 2018) data sample, passing tight identification and isolation requirements, or an electron with p T > 115 GeV and |η| < 2.4 with no isolation requirements. The muon trigger requires a muon with p T > 24 (27) GeV and |η| < 2.4 in 2016 and 2018 (2017), passing loose identification and isolation requirements. The muon trigger efficiency is larger than 99% for the smallest resonance masses and decreases to about 96% for the highest resonance masses.
In the offline selection, two well-identified leptons are first required. The leading and subleading leptons are required to have p T > 40 GeV in the range |η| < 2.4. Electrons (muons) must satisfy tight identification and tight (loose) isolation requirements. In events with more than two leptons selected, the two with largest transverse momentum are considered. Trigger, reconstruction, identification, and isolation efficiencies for leptons and their uncertainties are evaluated from dedicated data samples of leptonic Z boson decays, where one lepton from the decay serves as a tag and the efficiency for the other lepton is measured [59]. A leptonic Z candidate is defined as a pair of leptons with the same flavor, opposite charges, and invariant mass 76 < m < 106 GeV.
Depending on the boson p T in V → qq and H → qq decays, the decay products may either be reconstructed as a single AK8 merged jet (boosted category) or as a pair of AK4 jets (resolved category). In the following, j refers to AK4 jets and J refers to AK8 jets.
In the boosted selection, events must have a leptonic Z candidate with p T ( ) > 200 GeV and one AK8 jet with p T (J) > 200 GeV and soft-drop mass m J > 30 GeV. If several AK8 jets in the event satisfy the requirements, the one with the highest p T is considered. All AK8 jets have two subjets after the soft-drop mass requirement is satisfied. Boosted hadronic V and H candidates are first defined by τ 21 < 0.40 (0.45) in 2016 (2017 and 2018), and then by the soft-drop mass signal regions 65 < m J < 105 GeV (SR1) and 95 < m J < 135 GeV (SR2), respectively. Figure 2 shows the τ 21 and merged p T (J) distributions for boosted hadronic V and H candidates. In addition, boosted hadronic sideband (SB) candidates are defined by the τ 21 requirement and the soft-drop mass regions 30 < m J < 65 GeV and 135 < m J < 300 GeV.
In the resolved selection, boosted events are vetoed; events must have a leptonic Z candidate with p T ( ) > 150 GeV and an AK4 dijet system with p T (j) > 30 GeV, p T (jj) > 150 GeV, m jj > 30 GeV, and ∆R(jj) < 1.5. In events with more than two AK4 jets, the best dijet combination is chosen based on the following criteria: dijets in the b-tagged category are given priority; if more than one dijet combination remains, the pair closest in mass to the target signal boson is chosen. Resolved hadronic V and H candidates are defined by the dijet mass signal regions 65 < m jj < 110 GeV (SR1) and 95 < m jj < 135 GeV (SR2), respectively. In addition, resolved hadronic SB candidates are defined by the dijet mass regions 30 < m jj < 65 GeV and 135 < m jj < 180 GeV.
Boosted (resolved) events belong to the tagged category if the hadronic V or H candidate contains one subjet (jet) with a medium b tag and the other subjet (jet) with a loose b tag. Otherwise, they are placed in the untagged category. Figure 3 shows the distributions of the untagged, exclusive loose, and medium DEEPCSV tags for subjets of the boosted hadronic H candidates in SR2.
In total, eight categories are defined by the combination of electrons/muons, boosted/resolved, and tagged/untagged. Table 1 summarizes the selection and categorization requirements of the analysis. Figure 4 shows the merged jet m J and dijet m jj distributions for the untagged and tagged categories in the signal and sideband regions.

Background estimation
Most events in the selected ZV and ZH candidate samples come from Z + jets production. A genuine Z boson decaying into , along with a high-p T merged jet or dijet pair, may easily yield a signal-like combination, even though the jets do not originate from a boson decay. The features of the Z + jets background are studied using simulated samples. The limited number of simulated events remaining in the b-tagged categories results in sizable statistical fluctuations in the predicted distributions of the diboson mass m ZV and m ZH for this background. It has been observed, however, that within simulation uncertainties, the shape of the Z + jets mass distribution is the same for events with and without b-tagged jets. Therefore, the Z + jets shapes in the b-tagged categories are described using the m ZV and m ZH shapes obtained from the simulation without making any b tag requirements.
The modeling of the Z + jets background is obtained from simulation and corrected using data. In the final fit to the data, the Z + jets background normalization in the signal region, SR1 or SR2, is constrained by the observed yield in the SB. The shape of the Z + jets mass distribution is determined from the simulation, with a correction applied to improve the agreement with the data in the SB region. A linear function depending on one single parameter s, corr(m ZX , s) = 1 + s(m ZX − 500 GeV)/(500 GeV), is sufficient to make the invariant mass distributions agree within their uncertainties. This procedure is applied independently to each boosted/resolved and untagged/tagged category. Allowing for different linear corrections in the untagged and tagged categories takes into account possible differences in the shape of the Z + jets background. Figure 5 shows fits to the SB diboson mass distributions including the linear correction with s at its fitted value. Fitted values of s are in the range of −0.05 to +0.10, with s = 0 corresponding to no correction. The relative uncertainties in the SB-fitted slope parameter s are in the range 30-100%, depending on the category. In the final fit of each category, the slope parameter s is allowed to float, constrained by the data in the signal and sideband regions simultaneously. The corrections to the Z + jets shapes derived from the SB and the SR are found to be mutually consistent. Statistical uncertainties associated with the simulated Z + jets distributions are also taken into account in the fit. The fits in the boosted V categories include the peaking region of the background; Fig. 5 shows that the SB data in this particular region are described well by the fit. The systematic uncertainty in the estimation of the Z + jets background is given by the uncertainties in the slope correction parameter s and in the normalization.
Dilepton backgrounds that do not contain a leptonic Z boson decay are estimated from data using eµ events passing the analysis selection. This approach accounts for top quark pair production (tt), WW + jets, Z → ττ + jets, single top quark, and hadrons misidentified as leptons; we collectively refer to this set of backgrounds as t + X. The relative normalization and shape of the diboson mass distributions for ee and µµ events with respect to eµ events have been estimated from data on a top quark-enriched control sample and shown to be consistent with lepton flavor symmetry expectations. The t + X background contamination represents 0. negligible.
The SR1 m ZV and SR2 m ZH distributions for the boosted and resolved categories are shown in Figs. 6 and 7, respectively.

Systematic uncertainties
The systematic uncertainties influence both the normalization and shape of the background and signal distributions in the analysis.
The uncertainty in the background shape, which is the dominant systematic effect, is evaluated by comparing data and simulation after the fit in the signal and SB regions, and the residual shape difference is treated as an uncertainty. The background shape correction procedure has been explained in Section 6. The impacts of these uncertainties are quantified by calculating the change on the fitted signal cross section when a given parameter is displaced by ±1 standard deviation from its post-fit value, divided by the total uncertainty in the fitted signal cross section. The impacts of the background shape uncertainty in the boosted untagged, boosted tagged, resolved untagged, and resolved tagged categories are 11, 13, 3, and 3%, respectively, for a bulk graviton with a mass of 1 TeV. For an ALP linear ZZ (ALP chiral ZH) signal, the corresponding impacts are 42 (9), 42 (44), 16 (7), and 16 (23)%, respectively.
Uncertainties associated with the description in simulation of the trigger efficiencies, as well as the uncertainties in the efficiency for lepton reconstruction, identification, and isolation, are extracted from dedicated studies of events with leptonic Z boson decays, and amount to 1.5% for muons and 2% for electrons. The uncertainties in the lepton momentum and energy scales are propagated to the signal shape and normalization; the largest impact on the normalization is 0.2%.
Uncertainties in the jet energy scale and resolution [54] affect both the normalization and the shape of the signal samples. The momenta of the reconstructed jets are varied according to the uncertainties in the jet energy scale and resolution, and the selection efficiencies and signal shapes are reevaluated using these modified samples. This procedure results in a change of 1-2% in the normalization of the resolved categories; the largest change in the normalization of the boosted categories is 0.3%. These sources of uncertainty are responsible for the largest signal shape systematic uncertainty.
The uncertainty in the boosted V boson identification efficiency and the groomed mass scale and resolution [60] are responsible for a 7% signal normalization uncertainty in the merged category of the analysis. These are measured in data and simulation in an almost pure selection of semileptonic tt events where boosted W bosons produced in top quark decays are separated from the combinatorial tt background by means of a simultaneous fit to the soft-drop mergedjet mass. The uncertainties in the soft-drop mass scale and resolution are propagated to the groomed jet mass, and the impact on the selection efficiency of signal and ZV background is taken into account. The simulated τ 21 distributions have been found to be similar for boosted W, Z, and H bosons, allowing the same signal normalization uncertainty to be applied to all boosted bosons. An additional uncertainty affecting the signal normalization is included to account for the extrapolation of the uncertainties extracted from a tt sample at typical jet p T of 200 GeV to higher regimes, yielding an uncertainty of 2.6-6%. Extrapolation uncertainties are estimated by comparing the results from simulated event samples using HERWIG [61] and PYTHIA.
The scale factors (SFs) for the b tagging efficiencies and their corresponding systematic un-Entries / 50 GeV  certainties [58] are measured in data, as functions of p T , using samples enriched in b quark content, and their propagation to the signal region of the analysis produces an uncertainty of 12% in the tagged boosted category and 4% in the tagged resolved category. The uncertainties in the mistag efficiency SFs are also considered; the uncertainties in the b tagging and mistag efficiencies are treated as anticorrelated between the tagged and untagged categories.
For the t + X background, a 4% uncertainty is estimated by comparing the yield of eµ events with ee + µµ events in a top quark-enriched control region.
The different shapes of the m J (m jj ) distributions in the NLO and LO Z + jets simulations induce a 3 (5)% change in the ratio of the SR to SB normalizations in the boosted (resolved) categories. This difference is taken as a systematic uncertainty.
The effect of varying the factorization and renormalization scales by a factor of 2 is propagated to the signal shape and normalization. The impact on the signal acceptance is evaluated to be 0.1-0.3%, depending on the signal mass and analysis category. Varying the factorization and renormalization scales in the Z + jets simulation changes the shape of the distribution in the same way as the linear correction, by an amount that is less than the correction itself. The effect is therefore absorbed by the fit, and an explicit associated uncertainty is not required.
The systematic uncertainty associated with the PDFs used to generate the simulated samples is evaluated by varying the NNPDF 3.0/3.1 PDF set within its uncertainties, following the prescription in Ref. [62], and its effect is propagated to the signal shape and normalization, resulting in a measured uncertainty of 0.3% in the resolved categories and 1.5% in the boosted categories.
The effect of varying the scales and PDFs on the shape of the ALP distribution is also evaluated, but as it is found to be small (at most 3% in the tail of the distribution), it is neglected in the final fit to the data.
Additional systematic uncertainties affecting the normalization of backgrounds and signal from the contributions of pileup events and the integrated luminosity [63][64][65] are also considered.
The complete list of background and signal normalization systematic uncertainties in the analysis is reported in Table 2.

Results and interpretation
The signal selection efficiencies obtained from simulation are presented in Table 3. The efficiencies are relative to simulated events in which one Z boson decays to an electron or muon pair, the W boson or second Z boson decays to quarks, and the H boson decays according to the branching fractions of the SM.
Results are extracted from a combined maximum likelihood fit of signal and background to the m ZV or m ZH distribution in the 400-3000 GeV interval, simultaneously in all the categories used in the analysis. The systematic uncertainties discussed in Section 7 are included as nuisance parameters in the maximum likelihood fit, and the background-only hypothesis is tested against the combined background and signal hypothesis [66,67].
No statistically significant excess is observed with respect to the SM expectations. The lowest single-sided p-values are of the order of 2.0%, corresponding to a local significance at the level of 2 standard deviations, and occur for bulk gravitons with masses in the vicinity of 575, 775 and 1350 GeV, and for a W boson with a mass in the vicinity of 775 GeV. Upper limits at 95% confidence level (CL) on the signal production cross sections are set using the modified frequentist (CL s ) approach in the asymptotic approximation [66][67][68][69].
The observed and expected upper limits on the resonance cross section, multiplied by the branching fraction for the decay into one Z boson and a W or Z boson, σ G B(G → ZZ) or σ W B(W → ZW), are reported as a function of the resonance mass in Fig. 8, assuming a W or G produced in the narrow-width approximation. A WED bulk graviton is excluded up to masses of 1200 (expected 1150) GeV for κ = 0.5. This mass limit is similar to the one obtained by the CMS Collaboration in Ref. [16], and the present analysis also sets cross section limits on gravitons with a mass down to 450 GeV. The strongest W mass exclusion limits in HVT models A and B are those quoted in Refs. [15,16,20]. Figure 9 shows the observed and expected 95% CL upper limits on the ALP linear |c G c Z | (left) and the ALP chiral |c G a 2D | (right) coupling coefficients, as functions of the mass scale f a . For each black dot in Fig. 9, events with m ZV or m ZH > f a have been excluded from the fit, restrict- ing the analysis to the EFT mass consistency region. Limits on the coupling coefficients for f a > 3 TeV are calculated from the same fit result as f a = 3 TeV. We define the ALP mass range where our results are considered valid as the region where the cross section diverges by less than 10% from its asymptotic value, m a < 100 GeV. The expected and observed upper limits on σ(gg → a * → ZZ/ZH) for f a = 3 TeV are reported in Table 4. In the ALP linear model with f a ≥ 3 TeV, the observed (expected) 95% CL limit on |c G c Z |/ f 2 a is 0.0415 (0.0400) TeV −2 . The limit on |c Z |/ f a and |c G |/ f a in the case where the two are equal is 0.204 (0.200) TeV −1 ; for |c G |/ f a = 0.25 TeV −1 , the limit on |c Z |/ f a is 0.166 (0.160) TeV −1 .
In the ALP chiral model with f a ≥ 3 TeV, the observed (expected) 95% CL limit on |c G a 2D |/ f 2 a is 0.0269 (0.0281) TeV −2 . The limit on | a 2D |/ f a and |c G |/ f a in the case where the two are equal is 0.164 (0.168) TeV −1 ; for |c G |/ f a = 0.25 TeV −1 , the limit on | a 2D |/ f a is 0.108 (0.113) TeV −1 . These limits are valid for a Higgs boson with SM decays. The experimental upper limit on B(H → BSM) is 0.34 at 95% CL [70]; thus, relaxing the SM condition would increase the limits on |c G a 2D |/ f 2 a by a factor of 1.23 at most.

Summary
A search has been presented for heavy resonances decaying to ZZ or ZW, and nonresonant ZZ or ZH production (where H is the Higgs boson) mediated by axion-like particles (ALPs). The analysis is sensitive to resonances with masses in the range from 450 to 2000 GeV. Two categories are defined based on the merged or resolved reconstruction of the hadronically de-  No significant excess is observed in the data above the standard model expectations. Depending on the resonance mass, upper limits of 2-90 and 5-120 fb have been set on the product of the cross section of a spin-2 bulk graviton and the ZZ branching fraction, and on a spin-1 W signal and the ZW branching fraction, respectively. Upper limits on the nonresonant ALPmediated ZZ and ZH production cross sections for a new physics energy scale f a = 3 TeV and ALP masses m a < 100 GeV have been established at 162 and 57 fb, respectively. Depending on the value of the scale f a , upper limits on the product of the ALP coupling to gluons with the relevant coupling to ZZ or ZH of 0.02-0.09 TeV −2 have been set, valid for ALP masses m a < 100 GeV. These are the first limits based on nonresonant ALP-mediated ZZ and ZH production obtained by the LHC experiments.

Acknowledgments
We congratulate our colleagues in the CERN accelerator departments for the excellent performance of the LHC and thank the technical and administrative staffs at CERN and at other CMS institutes for their contributions to the success of the CMS effort. In addition, we gratefully acknowledge the computing centers and personnel of the Worldwide LHC Computing Grid and other centers for delivering so effectively the computing infrastructure essential to our analyses. Finally, we acknowledge the enduring support for the construction and operation of the LHC, the CMS detector, and the supporting computing infrastructure provided by the follow- [13] CMS Collaboration, "A multi-dimensional search for new heavy resonances decaying to boosted WW, WZ, or ZZ boson pairs in the dijet final state at 13 TeV", Eur. Phys. J. C 80 (2020) 237, doi:10.1140/epjc/s10052-020-7773-5, arXiv:1906.05977.