Explaining excesses in four-leptons at the LHC with a double peak from a CP violating Two Higgs Doublet Model

Extended scalar sectors with additional degrees of freedom appear in many scenarios beyond the Standard Model. Heavy scalar resonances that interact with the neutral current could be discovered via broad resonances in the tails of the four-lepton invariant mass spectrum, where the Standard Model background is small and well understood. In this article we consider a recent ATLAS measurement of four-lepton final states, where the data is in excess over the background for invariant masses above 500 GeV. We discuss the possibility that this excess could be interpreted as a"double peak"from the two extra heavy neutral scalars of a CP violating Two Higgs Doublet Model, both coupling to the $Z$ boson. We apply an iterative fitting procedure to find viable model parameters that can match the excess, resulting in a benchmark point where the observed four-lepton invariant mass spectrum can be explained by two scalar particles $H_2$ and $H_3$, with masses of 540 GeV and 631 GeV, respectively, being admixtures of the CP eigenstates. Our explanation predicts additional production processes for $t\bar t$, $W^+W^-$, $4b$ and $\gamma\gamma$, some of which have cross sections close to the current experimental limits. Our results further imply that the electric dipole moment of the electron should be close to the present bounds.


Introduction
The discovery of the scalar resonance with a mass of about 125 GeV, compatible with the predicted Higgs boson [1], completes the Standard Model (SM). While currently most searches for physics beyond the SM do not point at new physics beyond the SM, excesses in final states with leptons [2,3] and di-photons [4][5][6] indicate the possible existence of additional scalar degrees of freedom. The four-lepton final states are often referred to as the 'golden channel' when searching for heavy scalar resonances, due to small and controllable SM backgrounds. The four-lepton analyses from ATLAS [7] and CMS [8,9] show enhanced event rates in final states with high invariant mass. These analyses were combined in ref. [10] to demonstrate the compatibility of the data with a broad resonance structure around 700 GeV. This is compatible with a very recent interpretation of ATLAS data as a second Higgs excitation at 680 GeV [11].
Searches for heavy scalars above 600 GeV in γγ, Zγ, ZZ → 4 , top quarks, pairs of Higgs and W bosons were reviewed and discussed in ref. [12], supporting the claim made in ref. [10] of a resonance in ZZ around 700 GeV. This possible resonant enhancement is visible both in gluon and vector boson fusion channels as reported by ATLAS in ref. [7] The more recent analysis of four-lepton final states by the ATLAS collaboration also shows an enhancement of event rates with invariant masses above about 500 GeV [13]. This observation was considered in ref. [14] to corroborate the possible resonance around 700 GeV. The theoretical framework was the Georgi-Machacek model and a cross section for the heavy resonance was found to be ∼ 160 fb. It is interesting to note that a broader enhancement of the bbbb final state with invariant mass above 500 GeV is visible, which might be compatible with the observed 4 excess from ref. [15]. On the other hand, a recent analysis searching for heavy diboson resonances in semi-leptonic final states is compatible with the SM prediction [16], however here the backgrounds are at a much higher level and might cover up a possible enhancement.
In this article we consider an explanation of the four-lepton excess at invariant masses above 500 GeV by a "double peak" from the two extra heavy neutral scalars of a CP violating Two Higgs Doublet Model (THDM), where the latter extends the scalar sector of the SM by an additional scalar SU (2) L doublet field [17]. Due to CP violation, both extra heavy neutral scalars can couple to the Z boson, thus showing up in the four-lepton invariant mass spectrum. Recently, we studied a class of THDMs with CP violation in ref. [18], exploring the testability of CP violation at the LHC and evaluating the current constraints on the model parameters, pointing out the importance of the the four-lepton final state as "discovery channel". We will make use of the results of [18] to find a viable set of model parameters that can match the current four-lepton excess. The article is structured as follows: in section 2 we briefly review the model framework, in section 3 we describe our analysis and discuss the results, and in section 4 we conclude.

The model
The THDM was first discussed in ref. [19] to discuss the phenomenon of CP violation in the scalar sector. For a comprehensive review we refer the reader e.g. to ref. [20]. In THDMs, the scalar sector contains two SU (2) L -doublet fields, φ 1 and φ 2 , with identical quantum numbers under the SM gauge symmetry group: The components h i , i = 1, ..., 4, are real neutral fields, η + i , i = 1, 2 are complex charged fields, and v i , i = 1, 2 are the vacuum expectation values (vevs). The most general Lagrangian density for the model can be decomposed as where L SM,kin denotes the kinetic terms for SM gauge fields and fermions, L φ,kin denotes the kinetic terms for the two scalar fields φ i , i = 1, 2, V φ denotes the scalar potential, and Y φ contains the Yukawa terms that give rise to the couplings between the SM fermions and the scalar fields.
Since scalar decays into four-lepton final states (henceforth referred to as 4 ) come about from scalar decays into two Z bosons that in turn decay into leptons, we are interested in the interaction of the scalar fields with the neutral current. There are only two interaction eigenstates that can decay into two Z bosons. However, the most general (CP violating) form of the THDM leads to scalar mixing, and in the presence of CP violation all the three neutral mass eigenstates H i of the THDM can mediate the process pp → H i → ZZ → 4 , with i = 1, 2, 3.
The final state from the process pp → H i → ZZ → 4 , with 4 = + α − α + β − β (and where the considered lepton flavors are ± α,β = e ± , µ ± ) features an invariant mass that reflects the mass of the mediating H i . The mass eigenstate H 1 corresponds to the SM-like Higgs boson with m H1 125 GeV, whose four-lepton signal has been studied (see examples for ATLAS [21] and CMS [22] analyses) and H 2 and H 3 are assumed to be heavier. As discussed above, in a CP violating THDM the extra neutral scalars H 2 and H 3 give rise to two additional broad peaks in the 4 invariant mass spectrum -or to a broad "double peak" if the masses are not too separated. This is a feature that is usually not considered when fitting the THDM to the data.

Analysis and results
We consider the measurements of differential cross-sections in 4 events in the 139 fb −1 data set by the ATLAS collaboration [13]. From the results in the ATLAS publication we use the 4 differential cross sections, and in particular the invariant mass spectrum (M 4 ). We digitise the observed event rates, their errors, and the theory prediction.
We use the nine bins from M 4 between 500 GeV and 900 GeV, six of which show event counts in excess of the theory prediction. We create a sample of excess events by subtracting the theory prediction from the observed event rates (which means that "excess events" per bin can also be negative).

Numerical setup
As an explicit example we consider the same class of THDMs with CP violation and type I Yukawa structure as in ref. [18], to which we refer the reader for model details and notation. We calculate testable properties of the THDM via numerical tools that include the following recent constraints: 1 B physics data using FlavorKit [23], Higgs data using HiggsBounds [24] and HiggsSignals [25], and electric dipole moments with the formulae from refs. [29,30].
The process pp → H i → 4 for the heavy scalars with indices i = 2, 3 is calculated in Mad-Graph [31], including the effective gluon-Higgs vertex via SPheno [32,33] and QCD corrections [34]. We note at this point that the inclusive 4 invariant mass spectrum is simulated to include the interference between the scalars. A fast detector simulation is done with 500k events per sample with Delphes [35] using the standard ATLAS detector card. From the reconstructed events we read out the invariant mass spectra.

First approximation
We searched for points with masses m H2 and m H3 around 500 GeV and 700 GeV, respectively, and with total cross sections for the 4 final state that are of similar magnitude. We found the benchmark point P 1 with m P1 H2 = 535 GeV and m P1 H3 = 703 GeV and total 4 cross sections σ P1 H2→4 = 1.3 fb and σ P1 H3→4 = 0.86 fb. For H 2 and H 3 the signal selection efficiency based on the experimental selection criteria is found to be 4 ∼ 0.3 and not too dependent on the scalar mass. We simulated exclusively pp → H 2 → 4 and pp → H 3 → 4 from which we obtained the invariant mass spectra ρ 2 and ρ 3 , respectively. Notice that at this step, the interference between the scalars is not taken into account.
Next we performed a simple χ 2 analysis. We varied the signal peaks of the two spectra with the two parameters δm j , such that the masses are given by m Hj = m P1 Hj + δm j . We also introduce the signal strength multipliers s j which are multiplied with the invariant mass spectra ρ Hj . We construct the χ 2 : where i is the bin number, b i is the measured event rate δ obs, with b SM,i being the SM prediction, ρ Hj the signal distributions, and the parameters δm j and s j are varied to minimise the χ 2 .

Iterative analysis
Notice that the signal rate as described in (4) distorts the invariant mass spectrum and thus disconnects it from the underlying benchmark point. However, the distorted spectrum can be used to locate the masses and event rates that are preferred by the fit to the data. Consequently we use the best-fit values for the masses and the total 4 cross sections (converted from the fiducial cross sections using the signal selection efficiency and the integrated luminosity) as selection criteria to find a new benchmark point. From a fine grained scan in the model parameters we select a point P n that has masses m Hj and cross sections σ Hj →4 that are as close to the best-fit results for the masses and event rates of the previous point P n−1 as possible.
For each new point P n we create an inclusive 4 invariant mass spectrum ρ Pn incl with two peaks around m H2 and m H3 , including the interference between H 2 and H 3 . We remark that the interference with H 1 is negligible for the here relevant mass scales of H 2 and H 3 , which we verified by computation. We separate the spectrum into ρ Pn H2 and ρ Pn H3 at the minimum between the two peaks and fit the parameters δm j and s j via the two partial spectra as in eq. (3) to the data.
Once we have a benchmark point with a spectrum that provides a good fit, we carry out a Bayesian fit of the parameters δm j , j = 2, 3 and s j , j = 2, 3 to establish the Bayesian confidence limits on the parameters.
The black dots with error bars denote the difference between observed and predicted data from the four-lepton invariant mass spectrum in ref. [13].
parameters give rise to m P7 H2 = 540 GeV, m P7 H3 = 631 GeV, and σ P7 H2→4 = 0.45 fb, σ P7 H3→4 = 0.42 fb. We call this point the "best-fit benchmark point" as is provides a very good fit to the spectrum with a χ 2 = 5.76, which is to be compared to the SM value for all bins above 500 GeV. For the SM we get a χ 2 SM = 21.0 (16.9) corresponding to an upward fluctuation with a p-value of 0.007 (0.03) considering statistical errors only (all errors). The contribution of the two neutral scalar particles to the four-lepton invariant mass spectrum for the "best fit benchmark point" P 7 is shown in fig. 1. The striking feature of the spectrum is the wide range of M 4 that receives contributions from the two heavy scalars.
In general, scalars that mix with the Higgs doublet can also decay into other SM particles than Z bosons. Thus, after finding a good benchmark point to match the 4 invariant mass spectrum as reported by ATLAS, we explore the possibility of making quantitative predictions for the H i decays into tt, W + W − and γγ. Therefore, we show in the four panels of fig. 2 the projections of the total cross sections for the ∼ 2000 parameter space points from the very fine grained scan over THDM model parameters. The figure shows the 4 final state of H 2 and H 3 into ZZ and the color code of the points denotes the total production cross sections for W + W − (upper panels), tt (middle pannels) and γγ (lower panels).

Discussion
Our "best-fit benchmark point" P 7 provides an excellent fit to the excesses of events in the 4 spectrum with invariant masses between 500 GeV and 700 GeV observed by the ATLAS collaboration. The fit prefers H 2 and H 3 with similar masses and production cross sections for the 4 final state. Combined with the facts that (i) the THDM adds one CP-even and one CP-odd scalar to the SM field content and (ii) that the CP-odd field does not decay into ZZ this implies that the mass eigenstates must be strongly CP-mixed. The CP properties of the two scalars could be tested at the HL-LHC via correlations in final states with two tau leptons from the processes pp → H i → τ + τ − as discussed in ref. [18].
A comment on our statistical analysis is in order at this point. The central aim of this article is not to establish statistical evidence for the existence of a THDM signal in the M 4 spectrum. We  rather want to point out that if the here discussed signal were to become statistically significant the THDM with CP violation is a suitable model to explain it. We thus decided against combining the data sets from ATLAS with a similar analysis by the CMS collaboration [36], which also observes excesses in 4 final states at high invariant mass. Since the CMS analysis uses a smaller data set, including this result would at best increase the statistical significance of our result very slightly. As mentioned above, our analysis allows us to make quantitative predictions for decays of H j , j = 2, 3 into tt. Our "best-fit benchmark point" has cross sections σ H2→tt = 22 fb and σ H3→tt = 13 fb. This is smaller than the current uncertainty of the recent measurements of the total production cross section σ tt = 830 ± 36(stat) ± 14(syst) pb by ATLAS [26] and 791 ± 25 pb by CMS [27].
Let us now confront our "best-fit benchmark point" with the fact that no enhancement of semi-leptonic final states from W W or ZZ decays has been reported in ref. [16]. For our "best-fit benchmark point", the sum of W W and ZZ cross sections yields: This cross section is comparable with the 2σ upper limits on the production cross section from ref. [16], which are 470 fb for a mass of 500 GeV (for the example of a scalar radion).
The small apparent enhancement of the bbbb = 4b final state for invariant masses above 500 GeV as observed in ref. [15] could be another indication for the process pp → H i → ZZ, i = 2, 3. For the corresponding cross section one would expect that with the b-tagging efficiency b 0.7, and the selection efficiencies 4b ∼ 0.1, and where additional 4b production could come from pp → H i → 2H 1 → 4b. Eq. (7) results in a lower limit of 16.5 additional events in the 4b final state, which matches quite well the observed ∼ 20 events in excess of the background for M 4b ≥ 500 GeV. Next we comment on the H i → γγ channel, for which our "best-fit benchmark point" has cross sections σ Hj →γγ = 1.2 fb and 1.05 fb for j = 2, 3, respectively. The ATLAS search for resonances in diphoton final states places at 2σ upper limits on the production cross section of 1.15 fb and 0.83 fb for resonances with masses corresponding to 540 GeV and 631 GeV, respectively [28]. This implies that our "best-fit benchmark point" is slightly in tension with these limits, but may still be regarded as compatible. We also note that in the current data there indeed exist some upward fluctuations of the observed event counts at 540 GeV and at 680 GeV. In any case, future analyses of the diphoton spectrum with more data should be able to test this prediction of our "best-fit benchmark point".
Last but not least we consider the impact that CP violation in the scalar sector has for lowenergy observables. As mentioned above, large CP-violating phases are an implication of H 2 and H 3 having similar signal strengths in the 4 final state. Large CP phases imply that the THDM fields give rise to the Electric Dipole Moments (EDM) of SM particles, in particular for the electron, cf. e.g. refs. [29,30,37]. In our fine-grained parameter space scan we reject EDM for the electrons that are above the current exclusion limit of |d e | < 1.1 × 10 −29 ecm [38], and we find that the majority ( 90%) of all points has 10 −30 ≤ |de| ecm ≤ 1.09 · 10 −29 .

Conclusions
In this paper we have considered a "double peak" from a CP violating Two Higgs Doublet Model as explanation for the local excess in four-lepton events with invariant masses above 500 GeV as observed by ATLAS and CMS. Within a class of THDMs, we used an iterative fitting procedure to search for model parameters that give rise to heavy Higgs masses and signal strengths towards explaining the excess. The "best-fit benchmark point" (called P7) we found this way provides an excellent explanation for the data, with a χ 2 = 5.76 for 8 bins, predicting events in excess of the SM in the range from 500 GeV to around 700 GeV. It prefers a broad "double peak" in the invariant mass spectrum, with two resonances at 540 GeV and 631 GeV, respectively.
Interpreted in the context of a THDM, this would be an indication for CP violation in the scalar sector. The CP mixing is required to be close to maximal due to the comparable production cross section for the two processes pp → H 2 → 4 and pp → H 3 → 4 . The CP mixing of H 2 and H 3 could in principle be tested at the HL-LHC in the future, e.g. via the di-tau final state (cf. [18]).
Our "best-fit benchmark point" predicts additional tt, V V and γγ production channels with cross sections (summed over H 2 and H 3 ) of about 35 fb for tt, 199 fb for W W and 2.2 fb for γγ. Our results are compatible with present limits, and may be responsible for minor excesses in the 4b and γγ channels. Moreover, the parameter space that leads to an explanation of the observed 4 spectrum gives rise to electron EDMs that are close to the current experimental bounds, providing an example for a complementary way to test the scenario by low energy experiments.