Search for decays of the 125 GeV Higgs boson into a Z boson and a $\rho$ or $\phi$ meson

Decays of the 125 GeV Higgs boson into a Z boson and a $\rho^0$(770) or $\phi$(1020) meson are searched for using proton-proton collision data collected by the CMS experiment at the LHC at $\sqrt{s}= $ 13 TeV. The analysed data set corresponds to an integrated luminosity of 137 fb$^{-1}$. Events are selected in which the Z boson decays into a pair of electrons or a pair of muons, and the $\rho$ and $\phi$ mesons decay into pairs of pions and kaons, respectively. No significant excess above the background model is observed. As different polarization states are possible for the decay products of the Z boson and $\rho$ or $\phi$ mesons, affecting the signal acceptance, scenarios in which the decays are longitudinally or transversely polarized are considered. Upper limits at the 95% confidence level on the Higgs boson branching fractions into Z$\rho$ and Z$\phi$ are determined to be 1.04-1.31% and 0.31-0.40%, respectively, where the ranges reflect the considered polarization scenarios; these values are 740-940 and 730-950 times larger than the respective standard model expectations. These results constitute the first experimental limits on the two decay channels.


Introduction
In 2012 a boson with a mass near 125 GeV was discovered by the ATLAS and CMS Collaborations at the CERN LHC [1][2][3]. Soon after it was established that the properties of this particle are, within uncertainties, in agreement with those of the Higgs boson (H) in the standard model (SM) [4][5][6][7][8][9]. Decays of the Higgs boson into γγ, ZZ * , W ± W ∓ * , τ + τ − , and bb , as well as Higgs boson production via gluon-gluon fusion (ggH), via vector boson fusion (VBF), in association with a vector boson, and in association with a top quark-antiquark pair, have all been observed [10][11][12][13][14][15][16][17][18][19][20][21][22][23]. While many of the couplings between the Higgs boson and other particles have already been measured, the required sensitivity for measuring Yukawa couplings to second-and first-generation fermions has not yet been reached. Yukawa couplings to secondgeneration fermions are accessible via searches for the decay of the Higgs boson into µ + µ − or cc , both of which have been performed at the LHC [24-27]. The upper limit at the 95% confidence level (CL) for the decay into µ + µ − (cc ) is approximately 3 (70) times the SM expectation. In addition, Yukawa couplings to lighter fermions are also accessible via rare exclusive decays of the Higgs boson. One class of such processes is the decay of the Higgs boson into a photon and a vector meson [28][29][30]. Thus far, the γJ/ψ, γψ(2S), γΥ(nS), γρ, and γφ decays have been searched for [31][32][33]. The 95% CL upper limits on the branching fractions of the Higgs boson into γJ/ψ, γρ, and γφ are 2 orders of magnitude larger than their expected values in the SM. For the γψ(2S) and γΥ(nS) decays, the corresponding upper limits are, respectively, 3 and 5 orders of magnitude larger than the SM expectation.
A related class of rare decays is that of the Higgs boson into a heavy vector boson and a vector meson (V) [34,35]. Up to now only the decays of the Higgs boson into ZJ/ψ and Zη c have been studied experimentally [36]. As indicated in Fig. 1, several processes contribute to the decay of the Higgs boson into a vector boson and a meson. The formation of a vector boson and a meson via H → ZZ * or H → Zγ * decays (Fig. 1, left and middle) are indirect contributions to this process. We refer to the decay of the Higgs boson into light quarks that radiate a vector boson and form a bound meson state (Fig. 1, right) as the direct process. In the SM the indirect processes contribute the most to the decay of the Higgs boson into a heavy vector boson and a vector meson. The direct process is negligible in the SM as it is suppressed by a factor of up to m 2 q /m 2 H relative to the indirect contributions [30]. In that expression m q and m H denote the masses of the quark and of the Higgs boson, respectively. However, in scenarios beyond the SM where the Yukawa couplings to light fermions are enhanced, this direct process could contribute significantly to the Higgs boson branching fraction into a vector boson and a meson [34]. An example of a model beyond the SM with enhanced Yukawa couplings to light fermions is a version of the Giudice-Lebedev model of quark masses [37] that is modified to have two Higgs doublets. In this scenario Yukawa couplings to light quarks could be enhanced by up to a factor of 7 [38]. Other scenarios in which light-quark Yukawa couplings can be larger than predicted in the SM include a single Higgs doublet model with Froggatt-Nielsen mechanism [39,40] and Randall-Sundrum models of warped extra dimensions [41]. In addition, studies of the indirect processes are also of interest as these probe a different phase space from conventional H → WW * and H → ZZ * measurements, and therefore provide complementary information.
This paper describes a search for decays of the 125 GeV Higgs boson into a Z boson and a ρ (770) 0 meson (H → Zρ) or into a Z boson and a φ(1020) meson (H → Zφ). The branching fractions of these processes in the SM are small: B(H → Zρ) = (1.4 ± 0.1) × 10 −5 and B(H → Zφ) = (4.2 ± 0.3) × 10 −6 [34]. The search uses a sample of proton-proton (pp) collisions collected by the CMS experiment at √ s = 13 TeV from 2016 to 2018. The data set corresponds to an integrated luminosity of 137 fb −1 , or 35.9, 41.5, and 59.7 fb −1 collected in 2016, 2017, and 2018, respectively. In this search we select the dimuon and dielectron final states of the Z boson. For the ρ and φ mesons, we select decays containing exactly two charged hadrons, corresponding to the π + π − final state for the ρ meson and the K + K − final state for the φ meson. In the event reconstruction π ± and K ± are not explicitly distinguished. The main source of background events in this analysis is from Drell-Yan production of a Z boson in association with a genuine or misidentified meson candidate. For brevity we do not distinguish between particles and antiparticles in our notations of decay processes in the remainder of this paper.

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, each composed of a barrel and two endcap sections. Forward calorimeters extend the pseudorapidity (η) coverage provided by the barrel and endcap detectors. Muons are detected in gas-ionization chambers embedded in the steel flux-return yoke outside the solenoid.
The silicon tracker measures charged particles within the pseudorapidity range |η| < 2.5. It consists of 1856 silicon pixel and 15 148 silicon strip detector modules. The silicon pixel detector modules are arranged in four layers. In 2016, data were taken with a different detector configuration; at that time there were 1440 silicon pixel detector modules arranged in three layers. For nonisolated particles with transverse momentum in the range 1 < p T < 10 GeV and |η| < 1.4, the track resolution is typically 1.5% in p T [42].
Muons are measured in the pseudorapidity 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 [43].
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 [44].
Events of interest are selected using a two-tiered trigger system [45]. The first level, composed of custom hardware processors, uses information from the calorimeters and muon detectors to select events at a rate of around 100 kHz within a time interval of less than 4 µs. The second level, known as the high-level trigger, consists of a farm of processors running a version of the full event reconstruction software optimized for fast processing, and reduces the event rate to around 1 kHz before data storage.
A more detailed description of the CMS detector, together with a definition of the coordinate system used and the relevant kinematic variables, can be found in Ref. [46].

Event reconstruction
The products of pp collisions are reconstructed based on a particle-flow algorithm [47], which combines information from all subdetectors to reconstruct individual particle candidates. These particle candidates are classified as muons, electrons, photons, and charged and neutral hadrons.
The candidate vertex with the largest value of summed physics-object p 2 T is taken to be the primary pp interaction vertex (PV). The physics objects are the jets, clustered using the jet finding algorithm [48,49] with the tracks assigned to candidate vertices as inputs, and the associated missing transverse momentum, taken as the negative vector sum of the p T of those jets. Other collision vertices in the event are considered to have originated from additional inelastic pp collisions in each bunch crossing, referred to as pileup (PU). The average number of PU interactions during the 2016 data-taking period was 23, rising to 32 during the 2017 and 2018 data-taking periods. The muons, electrons, and charged hadron tracks used in the search presented in this paper are all required to originate from the PV.
Muons are reconstructed through a simultaneous track fit to hits in the tracker and in the muon chambers [43]. To suppress particles misidentified as muons, additional requirements are applied on the track fit quality and compatibility of individual track segments with the fitted track. Contamination from muons produced within jets is reduced further by requiring the muon to be isolated from hadronic activity in the detector. A relative isolation variable is defined as where ∑ charged p T refers to the scalar sum of the transverse momenta of all charged particles and ∑ neutral p T is the sum of the p T of neutral hadrons and photons. These two sums are calculated within a cone of radius ∆R = 0.4 around the direction of the muon, where ∆R = (∆η) 2 + (∆φ) 2 and ∆η and ∆φ are differences in pseudorapidity and azimuthal angle, respectively. The p T of the muon is excluded from these sums. To reduce the effects from PU, charged particles are only considered in the isolation sum if they are associated with the PV. The term 0.5 ∑ charged p PU T estimates the contributions from neutral particles in PU by summing the p T of charged particles that are within the isolation cone but are not associated with the PV. The factor 0.5 accounts for the ratio of neutral to charged particle production. Muons selected in the analysis must satisfy I µ rel < 0.15. After these identification and isolation requirements are imposed, prompt muons are identified with an efficiency of over 90%. A looser selection, where the isolation requirement is removed, is also used in the analysis to reject events with additional muons.
Electrons are reconstructed by combining clusters of energy deposits in the ECAL with hits in the tracker [44]. To reduce contamination from particles incorrectly identified as electrons, reconstructed electrons are required to pass a multivariate electron identification discriminant. This discriminant, based on the one described in Ref. [44], combines information about the quality of the tracks, the shower shape, kinematic quantities, and hadronic activity in the vicinity of the reconstructed electron. Isolation sums similar to those in Eq. (1) are also included among the discriminant inputs. Therefore no additional isolation requirements are applied. Using the requirements placed on the discriminant in this analysis, the electron identification efficiency is 80%. The rate at which other particles are misidentified as electrons is ≈1%. Looser requirements are used to reject events with additional electrons. Using this looser selection on the multivariate identification discriminant, the electron identification efficiency is 90% and other particles are misidentified as electrons at a rate of 2-5%.
The ρ and φ meson decay products are reconstructed using charged particle tracks measured in the tracker. The tracks are required to originate from the PV and to pass "high purity" reconstruction requirements. These requirements are based on the number of tracker layers with hits, the track fit quality, and the values of the impact parameters relative to their uncertainties. The algorithm is described in more detail in Ref. [42]. In the event selection, described in Section 5, we exploit the known masses of pions and kaons to calculate and restrict the invariant mass of the ρ and φ candidates.

Simulated samples
Samples of simulated Higgs boson events, produced via the ggH, VBF, W-associated (WH), and Z-associated (ZH) modes, are generated at next-to-leading order (NLO) in quantum chromodynamics using POWHEG 2.0 [50][51][52][53][54][55]. In some of the figures in this paper, and for the evaluation of corrections that account for differences between data and simulation, samples of simulated Drell-Yan Z → events are used. Here, refers to e or µ. These samples are generated at leading order using MADGRAPH5 aMC@NLO 2.2.2 (2.4.2) [56] for the 2016 (2017 and 2018) data-taking periods. All generated samples are interfaced with PYTHIA 8.212 [57] to model parton showering and hadronization. In the signal samples the decays H → Zρ or H → Zφ are also modelled using PYTHIA. These samples are used to build the signal model, which consists of binned templates. The NNPDF3.0 parton distribution functions (PDFs) [58] are used for the 2016 data-taking period. For the samples of signal events NLO PDFs are used, while for the Drell-Yan events leading order PDFs are used. For the 2017 and 2018 data-taking periods the NNPDF3.1 PDFs [59] at next-to-next-to-leading order are used for all samples. The description of the underlying event is provided by the CUETP8M1 tune [60] for the 2016 datataking period and by the CP5 tune [61] for the 2017 and 2018 data-taking periods. Additional PU interactions, generated with PYTHIA, are added to all simulated events in accordance with the expected PU distribution. All generated events are passed through a GEANT4-based [62] simulation of the CMS detector before being reconstructed with the same version of the CMS event reconstruction software as used for data.

Event selection
The final states considered in the selection are the µµππ and eeππ decays of the Zρ system, and the µµKK and eeKK decays of the Zφ system. The selection of the µµ and ee pairs, referred to as the dilepton system in what follows, is independent of the meson candidate under consideration. The trigger selection for the µµ final states is based on the presence of at least one isolated muon with p T > 24 GeV in the 2016 and 2018 data-taking periods, and at least one isolated muon with p T > 27 GeV in the 2017 data-taking period. For the ee final states the trig-ger selection requires the presence of at least one isolated electron with p T > 27 GeV in the 2016 data-taking period. In the 2017 (2018) data-taking period this threshold is p T > 35 (32) GeV.
After imposing the trigger requirements, events in the µµ channel are selected by requiring the presence of two oppositely charged muons passing the identification and isolation criteria described in Section 3. At least one of these muons must pass the trigger selection. Both muons must have p T > 20 GeV and |η| < 2.4, while the p T of the muon that satisfies the trigger requirements must be at least 3 GeV above the p T threshold at the trigger level. The ee channel selects events containing two oppositely charged electrons passing the identification criteria described in Section 3. At least one of the electrons must pass the trigger selection. Both electrons must have p T > 20 GeV and |η| < 2.1. The p T of the electron satisfying the trigger requirement must be at least 3 GeV above the trigger-level threshold. The requirement that the p T of the lepton passing the trigger selection is at least 3 GeV above the threshold in the trigger ensures we avoid the part of the phase space where the trigger efficiency increases rapidly. In both the µµ and ee channels, events that contain additional leptons with p T > 5 GeV that pass the loose identification criteria described in Section 3 are rejected. The invariant mass of the dilepton system is required to be in the range 60 < m < 120 GeV.
The ρ (φ) candidate is reconstructed from its decay into π + π − (K + K − ). As the ρ and φ mesons are both light compared with the energy released in the decay, the two charged particles produced in the decay are emitted with small angular separation ∆R, as illustrated in Fig. 2. The events shown in this figure are required to pass the selection criteria described so far. The small separation between the two tracks is exploited in the selection of the ρ and φ candidates. The meson candidate is selected as a pair of oppositely charged particle tracks, both with p T > 1 GeV and separated by ∆R < 0.1. In what follows a pair of oppositely charged particle tracks is also referred to as a ditrack system. The charged particle tracks are required to be separated from each of the Z boson decay products by ∆R > 0.3. In addition, at least one of the tracks must have p T > 10 GeV. Figure 3 shows the p T distribution for the track that has the larger transverse momentum out of the two tracks selected as the meson candidate. This distribution is shown in the H → Zρ and H → Zφ signal events and in the background from Drell-Yan events, illustrating how this requirement helps to reduce the background. If multiple track pairs pass these requirements, we calculate the four-momentum of each ditrack system and select the pair of tracks with the highest p T . This choice maximizes the proportion of signal events in which the correct meson candidate is selected. In all channels, the meson candidate is correctly identified in 98-99% of the signal events.
Furthermore, we require the ditrack system to be isolated. An isolation sum I trk is calculated as in a cone of radius ∆R = 0.3 around the direction of the ditrack system. Only tracks with p T > 0.5 GeV that are associated with the PV are considered, and the tracks forming the ρ or φ candidate are excluded from the sum. Events are selected if I trk < 0.5 GeV, thus with no track around the direction of the ditrack system, as explained previously. Figure 4 shows the distributions of the isolation sum for the data and for the simulated signal, after applying all selection criteria except for the ditrack isolation requirement. The ditrack invariant mass requirement discussed below is also applied. This figure illustrates the reduction in background events because of the isolation requirements. Only events in which the dilepton and ditrack four-body mass is in the range 120-130 GeV are shown. This range is expected to contain 95% of the simulated signal.
The invariant mass of the ditrack system is also used to reduce the contamination from back-   Figure 3: The transverse momentum distribution for the track that has the larger p T out of the two tracks selected as the ρ or φ candidate. The distribution is shown for events that pass the meson candidate selection described in the text, but not the requirement that one of the tracks must have p T > 10 GeV. This distribution is shown for the H → Zρ decay (dashed red), the H → Zφ decay (dotted blue) and the background from Drell-Yan events (solid black). All contributions are normalized to the same area. ground events. Events with a ρ candidate are selected if the invariant mass of the ditrack system is within 0.6 < m π π < 1 GeV, calculated assuming the mass of each particle equals m π ± = 139.6 MeV [63]. The full width at half the maximum of the m π π distribution is approximately 120 MeV in the simulated signal. Figure 5 (left) shows this invariant mass distribution in simulated H → Zρ events. The φ meson has a smaller natural width than the ρ meson, therefore it is possible to use a narrower mass window. This is also reflected in the mass resolu-  Figure 4: The ditrack isolation sum in the ππ (left) and KK (right) channels, combining the µµ and ee channels for all the data-taking years. The distribution in data, as well as in the simulated H → Zρ and H → Zφ signals is shown. A branching fraction of 10 (5)% for the H → Zρ (H → Zφ) signal is assumed. The isolation sum is shown after applying all selection criteria apart from the ditrack isolation requirement. The ditrack invariant mass requirement is also applied. Only events in which the dilepton plus ditrack invariant mass is in the range 120-130 GeV are considered. The dashed line indicates the boundary of the region used in the analysis, for which the isolation sum is required to be smaller than 0.5 GeV.
tion. The full width at half the maximum of the m KK distribution in simulated signal samples is approximately 5 MeV. To select events with a φ candidate, the mass of each particle is taken as m K ± = 493.7 MeV [63] and we require 1.005 < m KK < 1.035 GeV. Figure 5 (right) shows this invariant mass distribution in simulated H → Zφ events.
After these requirements the contribution from H → Zφ events in the ππ channel is smaller than 1% of the total signal in this channel when the SM branching fractions for H → Zρ and H → Zφ are considered. The same is true for contributions from H → Zρ events in the KK channel. After all selection criteria including the ditrack invariant mass requirements are applied, there is no overlap in the events selected by the ππ and KK channels.
The product of signal selection efficiency and acceptance ( A) corresponds to the fraction of simulated signal events that pass the selection. To calculate these values we use the nominal simulated sample, in which the decays of the H and Z bosons are modelled isotropically. On average over the three data-taking years, A in the µµππ (µµKK) channel is 15 (18)%. For the eeππ (eeKK) channel the average A is 8 (10)%. . These masses are calculated assuming the charged particle mass equals the pion mass in the ππ selection and assuming the charged particle mass equals the kaon mass in the KK selection. The events pass all selection criteria described in the text, apart from the requirements on the ditrack invariant mass window. The dashed lines indicate the region selected in the analysis.

Corrections applied to simulated samples
A correction is applied to the simulated events such that the PU distribution in simulation reproduces this distribution in data [64]. Corrections are also applied to the simulation to account for differences in the efficiencies of the trigger selection; of the ditrack isolation requirement; and of the lepton reconstruction, identification, and isolation between simulated events and data. These corrections, deviating from unity by a few percent, are measured using the "tagand-probe" method [65]. The ditrack isolation efficiency correction is determined in Z → µµ events. An uncertainty, described in Section 8, is applied to account for the difference between the phase space where the correction is measured and where it is applied. Energy scale corrections, which are smaller than 1%, are applied to the muons and electrons [43,44]. The event simulations model the decays of the H and Z bosons isotropically, and so do not take into account the impact of particle helicities. However, as there are only a few possibilities for polarizations in the final decay products, we calculate the angular distributions for extreme polarizations and reweight the signal events accordingly following the method described in Ref. [66]. The Z boson and the ρ or φ meson can either both be transversely polarized or both be longitudinally polarized. The two leptons always have opposite helicity in the rest frame of the Z boson. For each possibility the distribution of the polar angle between one of the pions or kaons and the meson, and between one of the leptons and the Z boson, is evaluated analytically. The signal templates are weighted to both of these distributions simultaneously. We ensure that the total normalization of the signal, before event selection, is preserved by the reweighting. However, the reweighting modifies the distribution of the kinematic variables, in particular by changing the lepton p T . Therefore the reweighting reduces (increases) the fraction of signal events that pass the selection criteria in the transversely (longitudinally) polarized case, and so this affects the final results. The change of the signal yield in the two extreme polarizations, relative to the scenario with isotropic decays, is given in Table 1. Table 1: The effect on the signal yields of reweighting to the extreme polarization scenarios, described in more detail in the text, relative to the scenario with isotropic decays. The change in the fraction of signal events that pass the selection criteria affects the final results of the analysis.

Signal and background modelling
The dilepton and ditrack four-body mass distribution, corresponding to the reconstructed Higgs boson mass and denoted m hh , where h refers to π or K, is used in the statistical inference. The signal and background are therefore modelled as a function of this observable in the range 118 < m hh < 168 GeV. More than 95% of the expected signal is contained in the range 120 < m hh < 130 GeV; the large tail used at higher masses helps to improve the stability of the background parameterization. As a result of the kinematic selection on the leptons and the meson candidates, the four-body mass distribution for the background changes from rising to falling between 115 < m hh < 118 GeV. For this reason the lower bound of the range is taken as m hh = 118 GeV. The full width at half the maximum of the m hh distribution in samples of simulated signal events amounts to 2-3 GeV, depending on the channel considered. The signal is described through a binned template, built from simulated events. Each bin has a width of 1 GeV in the four-body mass, which matches the binning used for the data.
The background to this search, consisting mainly of Drell-Yan events, is modelled using analytic functions. The values of the parameters of these analytic functions are obtained directly in the final signal extraction fit. Prior to the signal extraction fit we need to determine a set of functional forms that can parameterize the background in the different channels and datataking years. Two sidebands, 118 < m hh < 120 GeV and 130 < m hh < 168 GeV, are used for this. Because the sideband with m hh < 120 GeV is short, we verify that the chosen functional forms also describe the background in a control region where we require 1 < I trk < 2 GeV.
The fitted values of the function parameters in the control region are not required to be the same as those in the analysis phase space. In the control region the full four-body mass range 118 < m hh < 168 GeV is considered.
Chebyshev polynomials are used to describe the backgrounds. The order used depends on the channel and data-taking period, and ranges from 2 to 5. These orders are determined in the sidebands and the control regions described above using an F-test [67]. The results of the fit are shown in Section 9.
Alternative functions can be used to estimate the bias from the choice of a particular background parameterization. As alternatives we choose exponential functions, as well as a function of the form where x = m/ √ s, and p i are parameters of the fit. Here, m represents the four-body mass and √ s = 13 TeV. Such a function has also been used in searches for dijet resonances [68].
The possible bias from the choice of background parameterization is estimated by fitting the alternative function to the four-body mass sidebands. Pseudo-experiments are then drawn from this parameterization, and a signal expectation is added to each pseudo-data set. A maximum likelihood fit of the signal and background models to each pseudo-data set is performed using the nominal background model. This test is performed three times with branching fractions of 0, 2.5, and 5% for H → Zρ or H → Zφ. The difference between the extracted and injected branching fraction is found to be independent of the injected branching fraction. This difference is taken as the uncertainty due to possible bias in the choice of background parameterization. The bias is found to be small and is included in the analysis as a systematic uncertainty.

Signal extraction and systematic uncertainties
The results of this analysis are presented as upper limits on B(H → Zρ) and on B(H → Zφ). All limits quoted in what follows are set at the 95% CL. Limits are set using the modified frequentist CL s criterion [69,70], in which the profile likelihood ratio modified for upper limits [71] is used as the test statistic. In the limit setting procedure we make use of the asymptotic approximation [72].
Several systematic uncertainties are incorporated in the likelihood as nuisance parameters. They are described in this section and summarized in Table 2.
Most of the systematic uncertainties affect only the normalization of the simulated signal templates: i. The uncertainties in the integrated luminosity measurements are, respectively, 2.5, 2.3, and 2.5% for the 2016, 2017, and 2018 data-taking periods [73][74][75].
ii. Uncertainties in the muon identification, isolation, and trigger efficiency measurements arise from the method used to measure the efficiency, from the kinematic phase space in which the measurement is performed, and from the limited size of the simulated samples used for the measurement in simulation [43]. These uncertainties affect the normalisation of the simulated processes by ≈1% for all the data-taking periods.
iii. Uncertainties in the electron reconstruction, identification, and trigger efficiency measurements range from 2 to 3%, depending on the data-taking period. These uncertainties mainly arise from the method used for the efficiency measurement [44].
iv. The uncertainty in the tracking efficiency amounts to 4.6-4.8% (corresponding to 2.3-2.4% per track), depending on the data-taking period. This uncertainty is determined by comparing ratios of D * meson decay chains in data and simulation. The dominant components of the uncertainty come from limited sample sizes and the uncertainties in the SM predictions of these ratios.
v. The uncertainty in the ditrack isolation efficiency measurement is 2% for all three datataking periods. This uncertainty mainly arises from the method used to measure the efficiency.
vi. Theoretical uncertainties in the ggH production cross section amount to 3.9%, with uncertainties in the VBF, WH, and ZH production cross sections being, respectively, 0.4, 0.7, and 3.8% [34].
vii. Uncertainties from the choice of PDF and the value of the strong force coupling constant (α S ) depend on the Higgs boson production mode and range from 1.6 to 3.2% [34].
Four systematic uncertainties affect both the shape and normalization of the simulated signal templates: i. Uncertainties in the lepton energy scales are typically less than 0.3% for both muons and electrons [43,44].
ii. An additional uncertainty in the ditrack isolation efficiency measurement is applied. This uncertainty is taken as the difference between the ditrack isolation efficiency in the phase space where the correction is measured, and the efficiency as evaluated in the simulated signal. This uncertainty is in the range 1-6%, depending on the data-taking period.
iii. The uncertainty in the total inelastic cross section, used for correcting the PU profile in simulation to the profile in data, is 4.6% [64]. The overall effect on the normalisation of the simulated signal templates ranges from 0.5 to 1.5%, depending on the data-taking period and the channel considered.
iv. Uncertainties due to the limited number of simulated events are taken into account by allowing each bin of the signal template to vary within its statistical uncertainty, independently from the other bins. The largest possible bias from the choice of the function modelling the background is included in the likelihood as a modification of the number of expected events. The number of expected events in a given bin i is obtained as (B + ∆ bias )s i + b i , where s i is the number of signal events and b i is the number of background events. The parameter B is the branching fraction of the Higgs boson and the parameter on which we set limits. The parameter for the bias from the choice of background function is ∆ bias . It is subject to a Gaussian constraint with a mean of 0 and a width equal to the largest possible bias due to the choice of background function, which ranges from 0.01 to 0.20%. These values are obtained using the method described in Section 7.
Theoretical uncertainties in the production cross sections, and the uncertainties due to the choice of PDF and the value of α S are treated as correlated between the different data-taking periods. The uncertainty in the integrated luminosity measurement is treated as partially correlated between the different data-taking periods. The other experimental uncertainties are treated as uncorrelated between the different data-taking periods.

Results
To present results in terms of B(H → Zρ) and B(H → Zφ), the signal templates are normalized by taking into account the ggH, VBF, WH, and ZH production cross sections. The ggH cross section is calculated at next-to-next-to-next-to-leading order as 48.58 pb [34]. The cross sections for the other production modes are calculated at next-to-next-to-leading order in quantum chromodynamics and next-to-leading order in electroweak accuracy, and amount, respectively, to 3.78, 1.37, and 0.88 pb [34]. In addition, SM branching fractions of 3.37% are assumed for each of the Z → decays [63].
In the limit setting procedure we do not take into account potential contributions of Higgs boson decays into a Z boson and other vector mesons.
The four-body mass distributions in data and the background model are shown in Fig. 6. The expected H → Zρ (H → Zφ) signal, in the isotropic-decay scenario, at a branching fraction of 3.0 (0.7)% is also shown. In this figure the µµ and ee channels, as well as all three data-taking periods, are combined for visualization. In the statistical inference these channels are considered separately in a simultaneous fit. No significant excess above the background expectation is observed in either of the two searches.  Table 3 for B(H → Zρ) and in Table 4 for B(H → Zφ). While these limits are set on the total branching fractions B(H → Zρ) and B(H → Zφ), the results mainly probe the indirect process via the H → ZZ * decay as the direct decay process (Fig. 1, right) is greatly suppressed in the SM.  [15] ATLAS Collaboration, "Study of (W/Z)H production and Higgs boson couplings using H → WW * decays with the ATLAS detector", JHEP 08 (2015) 137, doi:10.1007/JHEP08(2015)137, arXiv:1506.06641.
[17] ATLAS and CMS Collaborations, "Measurements of the Higgs boson production and decay rates and constraints on its couplings from a combined ATLAS and CMS analysis of the LHC pp collision data at √ s = 7 and 8 TeV", JHEP 08 (2016) 045, doi:10.1007/JHEP08(2016)045, arXiv:1606.02266. [20] ATLAS Collaboration, "Observation of Higgs boson production in association with a top quark pair at the LHC with the ATLAS detector", Phys.