Branching fraction for $Z$ decays to four leptons and constraints on new physics

The LHC experiments have measured the branching fraction for $Z$ decays to four leptons (electrons or muons). We have combined these measurements with the result $\mathcal{B}(Z \to 4\ell) = (4.58 \pm 0.26) \times 10^{-6}$, allowing a precise comparison to the standard model prediction of $(4.70 \pm 0.03) \times 10^{-6}$. We use a minimal extension of the standard model to demonstrate that this combined value may be used to set stringent limits on new physics.


Introduction
Measurements of the branching fraction for Z decays to four leptons, B(Z → 4 ) where = e, µ, enable tests of standard model (SM) predictions in the rare-decay regime. Such tests are limited, however, by the precision of these measurements. Rare decays such as Z → 4 may contain contributions from physics beyond the standard model (BSM) which so far remain undetected. The Z → 4 process in particular may be linked to observed BSM phenomena such as the muon g−2 anomaly and violations of lepton flavor universality. A more precise measurement of B(Z → 4 ) thus offers not only an improved test of the SM, but also insight into these and other new physics processes by way of improved constraints.
The Feynman diagram for the SM Z → 4 decay in pp collisions is shown in Fig. 1. This process has been observed at the LHC by the CMS Collaboration, which has measured B(Z → 4 ) at √ s = 7 and 13 TeV [1][2][3], and the ATLAS Collaboration, which has performed a combined measurement using data samples at √ s = 7 and 8 TeV [4]. All four of these measurements are in agreement with SM predictions and may be combined to reduce the experimental uncertainty. The combined value may be used to improve existing upper limits on positive contributions to B(Z → 4 ) from hypothesized gauge bosons, such as the Z in L µ − L τ models [5,6] and certain leptophilic dark matter models [7]. We report a combined B(Z → 4 ) determined from the CMS and ATLAS measurements and demonstrate its efficacy in constraining a minimal new physics (NP) model. We begin with a review of the published data and our combination procedure. Next, we describe our NP model and present the resulting constraints. We discuss our findings in terms of how they may be applied to future NP searches and, finally, conclude with a summary of our results.

Measurement combination
We define the phase space region for our combined B(Z → 4 ) by the invariant mass requirements 80 < m 4 < 100 GeV for the four-lepton sum and m + − > 4 GeV for all opposite-sign, same-flavor lepton pairs. This region is chosen to most closely match those used for the published measurements, and is identical to that of the CMS results published in 2016 [2] and 2018 [3]. The region used for the 2012 CMS measurement reported in Ref. [1], however, is further restricted to m > 4 GeV for all lepton pairs (regardless of flavor and charge). This smaller phase space region has a smaller predicted B(Z → 4 ). In order to more accurately include this measurement in our combination, we scaled the central value using a kinematic acceptance correction (approximately 90%) calculated from Monte Carlo (MC) events generated with MadGraph5 aMC@NLO 2.6.1 [8]. The measurement reported by ATLAS in Ref. [4] corresponds to yet another region, defined by m + − > 5 GeV, in which B(Z → 4 ) is smaller by approximately 25%. However, the authors of Ref. [4] also provide a value for B(Z → 4 ) scaled to m + − > 4 GeV, which we use here 1 . This value is listed in Table 1 with the other measurements.
Measurement procedures differ across the analyses. The ATLAS and 2012 CMS analyses determine B(Z → 4 ) by comparing the number of observed signal Z → 4 events to the number of Z → + − events in their respective fiducial regions. By contrast, the 2016 and 2018 CMS analyses use a likelihood fit to measure the pp → Z → 4 cross section, σ(pp → Z → 4 ), and extract B(Z → 4 ) using a theoretical σ(pp → Z → + − ). Hence, these two measurements include separate uncertainties on theory and luminosity, which we do not attempt to correlate with the earlier measurements.
For the combination, we employ a best linear unbiased estimator (BLUE) method [10]. Each source of uncertainty is given a correlation coefficient r. The choices for these coefficients, listed in Table 2, are motivated by knowledge of the measurement process, the experiments, and the LHC. The only nonzero r are among the CMS measurements. The experimental uncertainty coefficient r (syst.) CMS = 0.3 is derived from the contributions tabulated in Ref. [2], and is defined as the correlated fraction of their sum in quadrature (i.e. the variance). The theory coefficient r (theo.) CMS = 1 is chosen to account for CMS analysis procedures which are standard across the measurements, and the luminosity coefficient r (theo.) CMS = 0.5 is chosen as a conservative estimate for correlations among luminosity uncertainty calculations.
This value may be compared to the value (4.59 ± 0.25) × 10 −6 , obtained by setting all r = 0, or to the extreme case r calculated at leading order (LO) with MadGraph5. Here, Γ Z→4 is the partial width for Z → 4 decays subject to m + − > 4 GeV, determined from a sample of 100 000 generated decays, and Γ Z is the total Z width in the program's default SM implementation, which differs slightly from the nominal value. The theoretical uncertainty on Γ Z→4 is evaluated by increasing and decreasing the renormalization scale by a factor of two, and is comparable to the statistical uncertainty.

New physics model
To test the constraining power of the combined B(Z → 4 ), we developed a simple, minimal extension of the SM which includes additional contributions to the Z → 4 decay. These contributions arise from the addition of a single new boson resonance, named U , which may appear in the qq → Z → 4 process as shown in Fig. 2. The appearance of U in such a decay is, fundamentally, permitted only in the case that U is uncharged, less massive than the Z boson, and coupled to one or both lepton flavors. These are the only restrictions we consider. We have modeled U as both a scalar and a vector boson. Both NP models were created with FeynRules 2.3.29 [11] and implemented at LO with MadGraph5. The scalar and vector models, respectively, are described by the additions to the SM Lagrangian and where M U is the mass of the U boson and g (= g e , g µ ) the strengths of its lepton couplings. These quantities are free parameters, which we vary over the parameter space defined by 1 ≤ M U ≤ 90 GeV, 0.005 ≤ g ≤ 0.45. The full width of U , Γ U , depends on these parameters and is calculated automatically by MadGraph5. For the scalar (vector) model, 1.9 × 10 −6 ≤ Γ U ≤ 1.5 GeV (5.3 × 10 −6 ≤ Γ U ≤ 3.9 GeV). The largest values of Γ U in these ranges are eventually excluded. For all values considered here, the leptons are prompt.
To explore the models' behavior as a function of the input parameters, we examined the quantity ∆Γ Z→4 , defined as the deviation of the Z → 4 partial width from the SM value. That is, where the quantity denoted "NP" is calculated from the NP model for a given choice of free parameters, and that denoted "SM" is the SM prediction. These calculations were performed with MadGraph5, which gives Γ SM Z→4 = 12.2 eV. Our calculations take interference into account, and ∆Γ Z→4 is almost always positive.
Examples of the behavior of ∆Γ Z→4 as a function of M U and g (≡ g e = g µ ) are shown in Fig. 3. For both models, ∆Γ Z→4 exhibits the expected power-law dependence on each parameter.

Constraints
In order to apply constraints from the measurements to the NP models, we define an acceptance for the fiducial region used for the 2016 and 2018 CMS measurements [2,3], which correspond to our chosen phase space. This fiducial region is defined by-in addition to the phase space T > 20 GeV, p 2 T > 10 GeV, p 3,4 T > 5 GeV; the pseudorapidity requirement |η | < 2.5 for each lepton; and the invariant mass requirement m Z 1 > 40 GeV, where Z 1 is the opposite-sign, sameflavor lepton pair with invariant mass closest to the nominal Z mass. We define our exclusion as the parameter space region such that where σ = σ(pp → Z → 4 ) is the cross section and A the acceptance, and the labels "NP" and "SM" are used as in Eq. 5. The quantity B is our combined branching fraction given in Eq. 1, and B UL is a corresponding upper limit. At the 95% confidence level (CL), for instance, B UL = 5.01 × 10 −6 . Exclusion curves were interpolated from a parameter space sampled in intervals of 10 GeV in M U and 0.05 in g. Approximately 20 000 pp → Z → 4 MC events were generated at each sampled point, and the quantities σ NP and A NP calculated from them at the generator level. For the scalar (vector) model, they fall in the ranges 0.186 ≤ σ NP ≤ 3.66 pb (1.09 × 10 −9 ≤ σ NP ≤ 209 pb) and 9.9 ≤ A NP ≤ 14.3% (3.7 ≤ A NP ≤ 32.6%). The SM quantities σ SM = 194 fb and A SM = 10.5% were obtained by the same procedure. Figure 4 shows the exclusion curves for both models in the case g e = g µ (≡ g), corresponding to B UL at the 95%, 90%, and 68% CLs. Figure 5 compares the g e = g µ case to the case g e = 0 at the 95% CL.

Discussion
As illustrated in Fig. 6, our combined value for B(Z → 4 ) lies well within the range of the measurements listed in Table 1. The total uncertainty on the combined value is comparable to Figure 4: The excluded regions in (M U , g) space when g e = g µ (≡ g), compared at the 95%, 90%, and 68% CLs. The scalar U model is shown in blue and the vector U model in red.   Table 1, with colored error bars corresponding to each source of uncertainty. Also shown are the SM predicted value (thin gray line) and an upper limit on the combined value at the 95% CL (dotted gray line).
the statistical uncertainty on the most recent result published by CMS [3], and is smaller than the statistical uncertainties on all of the other measurements. This reduced uncertainty leads to a tight upper limit and the exclusion of a sizable region in parameter space. The boundary of this excluded region does not depend strongly on the CL of the upper limit, especially for the vector U model, as illustrated in Fig. 4.
Holistic statements about the constraining power of B(Z → 4 ) must consider experimental results. For example, LEP measurements rule out a U coupling to electrons within the energy reach of the LHC [12]. Therefore, bounds addressing a new boson which couples only to muons are of more practical interest in the search for physics to describe observed BSM phenomena such as those mentioned in Sec. 1. Figure 5 demonstrates that the exclusion provided by B(Z → 4 ) is not too diminished in this regime. While B(Z → 4 ) includes contributions from the nonconstructive Z → 4e and Z → 2e2µ channels, it still offers a valuable test of the SM in the Z → 4µ channel which may be applied to future NP searches at the LHC. These searches would complement those already performed at lower-energy experiments, such as BABAR; the collaboration has published limits on light Z bosons with g µ as low as 7 × 10 −4 [13], corresponding to M U < 10 GeV in our parameter space.

Conclusion
We have combined the CMS [1][2][3] and ATLAS [4] measurements of the Z branching fraction to four leptons, B(Z → 4 ) where = e, µ, with the result B(Z → 4 ) = (4.58 ± 0.26) × 10 −6 . This combined value corresponds to the phase space 80 < m 4 < 100 GeV, m + − > 4 GeV and takes into account correlations among the measurements. We have used this combined value to constrain the parameters of a minimal SM extension including a new scalar or vector boson U , which may couple to one or both lepton flavors. These constraints exclude a large region of the parameter space (M U , g ), where M U is the mass of the U boson and g its lepton couplings. These limits demonstrate the potential of the combined B(Z → 4 ) value to constrain theoretically and experimentally informed NP models, as well as searches for BSM physics at the LHC.