Search for new neutral Higgs bosons through the H $\to$ ZA $\to$ $\ell^{+}\ell^{-} \mathrm{b\bar{b}}$ process in pp collisions at $\sqrt{s} =$ 13 TeV

This paper reports on a search for an extended scalar sector of the standard model, where a new CP-even (odd) boson decays to a Z boson and a lighter CP-odd (even) boson, and the latter further decays to a b quark pair. The Z boson is reconstructed via its decays to electron or muon pairs. The analysed data were recorded in proton-proton collisions at a center-of-mass energy $\sqrt{s} = $ 13 TeV, collected by the CMS experiment at the LHC during 2016, corresponding to an integrated luminosity of 35.9 fb$^{-1}$. Data and predictions from the standard model are in agreement within the uncertainties. Upper limits at 95% confidence level are set on the production cross section times branching fraction, with masses of the new bosons up to 1000 GeV. The results are interpreted in the context of the two-Higgs-doublet model.


Introduction
The CMS and ATLAS experimental programmes are focusing efforts on the measurement of the properties of the Higgs boson discovered in 2012 [1][2][3], which has a mass of about 125 GeV [4][5][6]. All measurements to date are consistent with the expectations for a standard model (SM) Higgs boson within the experimental uncertainties.
Additional Higgs bosons are predicted in several extensions of the SM. Examples of these extensions are the two-Higgs-doublet model (2HDM) [7], whose phenomenology is based on the presence of an additional scalar Higgs doublet, and the minimal supersymmetric extension of the SM (MSSM) [8], which is a particular realisation of the 2HDM. The two Higgs doublets entail the presence of five physical states: two neutral and CP-even bosons (h and H); a neutral and CP-odd boson (A); and two charged scalar bosons (H ± ). Under particular theoretical assumptions [7], the model is often described by the following parameters: the mass of the CPeven boson H, m H ; the mass of the pseudoscalar A, m A ; the ratio of the vacuum expectation values of the two doublets, tan β; the mixing angle α between the two CP-even bosons; and the soft-breaking term, m 2 12 . Different couplings of the two doublets to right-handed quarks and charged leptons are predicted in various formulations of the 2HDM: in the Type-I formulation, all fermions couple to only one Higgs doublet; in the Type-II formulation, the up-type quarks couple to a different doublet than the down-type quarks and leptons; in the "lepton-specific" formulation, the quarks couple to one of the Higgs doublets and the leptons couple to the other; and in the "flipped" formulation, the up-type quarks and leptons couple to one of the Higgs doublets, while the down-type quarks couple to the other.
Different models and assumptions also alter the mass hierarchies, as shown in Fig. 1. There, and in the rest of the paper, h is identified with the observed Higgs boson. Two scenarios are possible. In the conventional scenario, the pseudoscalar is degenerate in mass with the charged scalars and is heavier than the scalar H, thus allowing for the A → ZH process; while in the twisted [9] scenario, the scalar H is degenerate in mass with the charged scalars and is heavier than the pseudoscalar, thus allowing for the H → ZA process. Moreover, in the parameter space region where cos(β − α) approaches 0, the CP-even h has properties indistinguishable from a SM Higgs boson with the same mass. In this region, known as the alignment limit, the branching fraction of the heavy scalar H to a Z boson and a lighter pseudoscalar A is the largest. The branching fractions for several decay channels of the H and A bosons for m H = 300 GeV and m A = 200 GeV are shown in Fig. 2 as functions of cos(β − α) (left) and tan β (right). This paper reports on a search for a new CP-even (odd) neutral Higgs boson decaying into Z and a lighter CP-odd (even) neutral Higgs boson, where the Z decays into an opposite-sign electron or muon pair, and the light Higgs boson into a b quark pair. The analysis is performed under the assumption of the twisted mass hierarchy scenario, and subsequently extended to the conventional scenario by interchanging the masses of the two bosons. While this is not crucial when setting model independent upper limits, it appears particularly interesting for theoretical interpretation of the results in the context of the 2HDM, since the theoretical cross sections differ in the two scenarios. The search is based on LHC proton-proton (pp) collision data at a center-of-mass energy √ s = 13 TeV collected by the CMS experiment during 2016, corresponding to an integrated luminosity of 35.9 fb −1 . The analysis exploits the invariant mass distributions of the bb, with being an electron or muon, and bb systems to search for a resonant-like excess of events compatible with the H and A masses.

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-ionisation chambers embedded in the steel flux-return yoke outside the solenoid. A two-level trigger system [16] is used to reduce the rate of recorded events to a level suitable for data acquisition and 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. [17].

Event simulation and background predictions
Background samples for this search are produced for Z boson production through the Drell-Yan (DY) process, top quark pair production (tt), single top quark, diboson, triboson, ttV (V = W, Z), W+jets, and SM Higgs boson production. They are generated at next-to-leading order (NLO) precision in perturbative quantum chromodynamics (QCD  [26] for parton shower and hadronisation. The underlying event tune is CUEPT8M1 [27], derived from the MONASH tune [28]. The NNPDF 3.0 [29] LO and NLO parton distribution functions (PDF) are used.
Signal samples of 207 different mass hypotheses are produced for the process H → ZA → bb, with m H and m A ranging from 120 to 1000 GeV and from 30 to 1000 GeV, respectively. The choice of the mass hypotheses is strongly motivated by the need of achieving a complete coverage of the parameter space. The spacing between two adjacent mass hypotheses is chosen so as to take into account the worsening of the signal resolution as the mass increases, such that the signal shape can be interpolated with good accuracy over the whole search region. These samples are produced using MADGRAPH5 aMC@NLO version 2.3.2 [18] interfaced with PYTHIA 8.212 [26] for simulation of the parton shower, hadronisation, and underlying event using the CUEPT8M1 tune [27]. The PDF set used is NNPDF 3.0 [29] at leading order (LO) in the four-flavour scheme, and the factorisation and renormalisation scales are estimated dynamically. The total widths assumed for the H and A resonances are computed with 2HDMC [30]. Finally, the signal simulation does not account for production of the scalar H in association with b jets.
For all processes, the detector response is simulated using a detailed description of the CMS apparatus, based on the GEANT4 package [31]. Additional pp interactions in the same and or neighbouring bunch crossings (pileup) are generated with PYTHIA 8.212 [26], and overlapped with the simulated events of interest in order to reproduce the pileup measured in data.
All background processes are normalised to their most accurate theoretical cross sections. The tt, DY, single top quark, W + W − , and W+jets samples are normalised to next-to-next-to-leading order (NNLO) precision in QCD [32][33][34][35], while the remaining diboson, triboson and ttV processes are normalised to NLO precision in QCD [18,36]. The SM Higgs boson production cross section is computed at NNLO QCD precision and NLO electroweak precision [37]. We indicate the SM Higgs boson, the ttV, and the W+jets backgrounds with Other in the figures.

Event reconstruction and selection
Events considered for this search are selected by a trigger based on the dilepton signature. The leading and subleading transverse momentum (p T ) thresholds applied by the triggers are channel dependent, and vary from 17 to 23 GeV (8 to 12 GeV) for the leading (subleading) lepton. Trigger efficiencies are measured with a "tag-and-probe" method [38] as a function of lepton p T and η in a data control region consisting of Z → events. Events with two oppositely charged leptons (e ± e ∓ , µ ± µ ∓ ) are selected using asymmetric p T requirements, chosen to be above the corresponding trigger thresholds, for the leading and subleading leptons. These requirements are 25 and 15 GeV, respectively, for e ± e ∓ events; and 20 and 10 GeV, respectively, for µ ± µ ∓ events. Electrons in the range |η| < 2.5 and muons in the range |η| < 2.4 are considered. Events with different-flavour leptons (e ± µ ∓ ) are also selected. The p T requirement for the leading lepton is 25 and 15 (10) GeV for the subleading electron (muon). These events mostly arise from tt production, and this region is used in the final template fit described in Section 7 to obtain an estimate of the normalisation of the non-resonant background processes (tt, single top quark, diboson, and triboson) and of the shape of the tt process only. For simplicity, we will refer to events with two electrons, muons, and mixed-flavour leptons as ee, µµ, and eµ, respectively, throughout this paper.
A particle-flow (PF) algorithm [39] aims at reconstructing all particles (PF candidates) in an event by combining information from all subdetectors. The PF candidates include photons, electrons, muons, neutral hadrons, and charged hadrons. The candidate vertex with the largest value of summed physics-object p 2 T is taken to be the primary pp interaction vertex. The physics objects used in this determination are the jets, clustered using the jet finding algorithm [40, 41] 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. Electrons, reconstructed by associating tracks with ECAL clusters, are identified by a sequential selection using information on the cluster shape in the ECAL, track quality, and the matching between the track and the ECAL cluster. Additionally, electrons from photon conversions are rejected [42]. Muons are reconstructed from tracks found in the muon system, associated with tracks in the silicon tracking detectors. They are identified based on the quality of the track fit and the number of associated hits in the various tracking detectors [43]. The Rochester correction [44] is applied to the muon momenta to correct for misalignments of the detector or uncertainties in the magnetic field. The lepton isolation, defined as the scalar p T sum of all PF candidates in a cone of radius ∆R = 0.4 around the lepton, excluding the lepton itself and corrected for contributions from particles not coming from the primary vertex, divided by the lepton p T , is required to be <0.06 for electrons and <0.15 for muons. Here, ∆R is defined in terms of the track separation in η and azimuthal angle (φ, in radians) as ∆R = √ (∆φ) 2 + (∆η) 2 . Moreover, the lepton tracks are required to be connected to the primary vertex. Lepton identification and isolation efficiencies in the simulation are corrected for residual differences with respect to data. These corrections are measured in a data sample enriched in Z → events, using the "tag-and-probe" method, and are parameterised as a function of lepton p T and η. The DeepCSV algorithm [48] is used to identify jets originating from b quarks. Jets are considered as b tagged if they have p T > 20 GeV and they pass the medium working point of the algorithm, which provides around 70% efficiency with a mistag rate of less than 1%, while the mistag rate for c jets is around 10%. Correction factors are applied in the simulation to the selected jets to account for the different response of the DeepCSV algorithm between data and simulation [48]. Among all possible dijet combinations fulfilling the previous criteria, we select the two jets with the highest DeepCSV algorithm outputs.
The final object selection consists of two opposite-sign leptons and two b-tagged jets, after which a requirement of 70 < m < 110 GeV is applied to enhance the presence of Z → events. In addition, the events are required to have a p miss T < 80 GeV in order to reduce the background contributions from processes with large p miss T , such as tt production. Both requirements have negligible impacts on the signal efficiency.
The main background processes, in decreasing order of importance, are DY in association with b quarks and tt production where both top quarks decay leptonically (fully leptonic tt). The contribution from QCD multijet events with jets misidentified as leptons constitutes a negligible background after requiring a pair of well-identified leptons, as described in Section 4.

Signal extraction
We search for the process H → ZA → bb by fully reconstructing its final-state objects and applying selection requirements in order to remove as many background events as possible, as explained in Section 4. From the reconstructed objects, we search for resonances in the invariant masses. Specifically, the invariant mass of the A can be reconstructed from the b jet pair; and that of the H from the b jet pair and the lepton pair. Two categories are defined based on the lepton flavours considered: ee and µµ. The Z mass, reconstructed from two oppositesign leptons, is used in the selection criteria described in Section 4 since it is common to all signals studied in this paper. The masses of the other two particles, H and A, vary according to the signal scenarios considered. Therefore, a simple and effective model independent approach to isolate the signal is to search for an excess of events in the reconstructed m jj and m jj distributions centered around the H and A candidate mass for each signal hypothesis. These distributions for µµ + ee events are shown in Fig. 3, where the background shapes and normalisations are obtained from simulation.
Since the m jj and m jj distributions are inherently positively correlated under a particular signal hypothesis, an elliptical signal region is chosen in order to optimize the sensitivity of the search. Figure 4 (left) shows the reconstructed mass distributions for three different signals in the m jj vs. m jj plane along with their defined elliptical signal regions. Because the shape of the signal is driven by the energy resolution of the final-state objects, ellipses take different sizes and tilt angles, depending on the masses being considered. A parametrisation is therefore performed in order to guarantee a good description of the signal shape for each signal hypothesis. For each ellipse, it provides the center, the major and minor semi-axes, and the tilt angle. Since each ellipse must be well-centered around the maximum of the two-dimensional (2D) mass distribution, the reconstructed center is extracted from a one-dimensional Gaussian fit in both m jj and m jj . The diagonalisation of the covariance matrix of the 2D distribution provides the axes of the ellipse and its tilt angle.
Since the shape of the signal is not exactly Gaussian, concentric elliptically shaped regions are defined in the parameter space using a parameter called ρ. Specifically, an ellipse with ρ = i contains roughly the fraction of signal events expected within i standard deviations in a 2D distribution. Selected events in the m jj vs. m jj plane are classified in six regions around the center of the ellipse defined for each signal point. The regions are built in ρ steps of 0.5, from 0 to 3, as illustrated in Fig. 4 (right), and lead to a template containing six bins used to perform the statistical analysis. By construction, the bulk of the signal is located at small values of ρ. The yield in data and the expected yields in simulation are reported in Table 1 for each elliptical bin under the mass hypothesis m H = 500 GeV and m A = 200 GeV. The ee and µµ categories are summed. Table 1: Expected and observed event yields prior to the fit in the signal region with m H = 500 GeV and m A = 200 GeV for each elliptical bin. The signal is normalised to its theoretical cross section for the Type-II 2HDM benchmark tan β = 1.5 and cos(β − α) = 0.01. The ee and µµ categories are summed. The reported uncertainties are statistical only.

Systematic uncertainties
We consider different sources of systematic uncertainties that may affect the statistical interpretation of the results, through their modification of both the normalisation and the shape of the distributions for the signal and background processes.
Theoretical uncertainties in the cross sections of the background processes estimated using simulation are considered as systematic uncertainties in the yield predictions. The uncertainty in the total integrated luminosity is determined to be 2.5% [49].
The signal region contains events that have at least two b-tagged jets. A control region is built by requiring events to pass the selection, as described in Section 4, but with no b tag requirement for the jets. In that region, a discrepancy between data and simulation of up to 10% is observed in the shape of the mass distributions, which hints at a mismodeling of the DY + heavyflavour jets background in some specific regions of the reconstructed mass plane. To account for this mismodeling, given the similar precision of the DY + light-flavour and DY + heavy-flavour jets background samples at simulation level, the observed data-MC discrepancy is fitted with a polynomial function, which is used to reweight each DY + heavy-flavour jets simulated event in the signal region, and a shape uncertainty equal to 100% of the correction is applied. In order to avoid assigning only one shape uncertainty to regions characterised by very different values of the above-mentioned correction, this uncertainty is considered independently in 42 regions of approximately 150 × 150 GeV 2 in the m jj vs. m jj plane. This procedure ensures enough degrees of freedom in the maximum likelihood fit (used to extract the best fit signal cross section, as explained in Section 7) to properly account for the mismodeling of the DY + heavy-flavour jets background shape.
The following sources of systematic uncertainties that affect the normalisation and shape of the templates used in the statistical evaluation are considered: • Trigger efficiency, lepton identification and isolation: uncertainties in the measurement of trigger efficiencies, as well as electron and muon isolation and identification efficiencies, are considered. These are evaluated as a function of lepton p T and η, and their effect on the analysis is estimated by varying the corrections to the efficiencies by ±1 standard deviation. • Renormalisation and factorisation scale uncertainty: this uncertainty is estimated by varying the renormalisation (µ R ) and the factorisation (µ F ) scales used during the generation of the simulated samples independently by factors of 0.5, 1, or 2. Cases where the two scales are at opposite extremes, are not considered. An envelope is built from the 6 possible combinations by keeping maximum and minimum variations for each bin of the distributions, and is used as an estimate of the scale uncertainties for all the background and signal samples.
• PDF uncertainty: the magnitudes of the uncertainties related to the PDFs and the variation of the strong coupling constant for each simulated background and signal process are obtained using variations of the NNPDF 3.0 set [29], following the PDF4LHC prescriptions [34].
• Drell-Yan additional uncertainty: additional shape uncertainties are applied to DY events to correct for mismodeling of this background as explained above. Their values range up to 10%, depending on the region of the reconstructed mass plane.
• Simulated sample size: the finite nature of simulated samples is considered as an additional source of systematic uncertainty. For each bin of the distributions, one additional uncertainty is added, where only the considered bin is altered by ±1 standard deviation, keeping the others at their nominal value.
The variations that these uncertainties induce on the total event yields in the analysis selection are summarised in Table 2 for a specific signal hypothesis, where the ee and µµ categories are combined together, yielding, for some uncertainties, a range of variations.

Methodology and results
A binned maximum likelihood fit is performed in order to extract best fit signal cross sections. The fit is performed using the six binned templates mentioned above in the ee and µµ channels. An additional six-bin template is included in the fit containing the eµ selection to further constrain the tt process, which is the major background component in this region, together with the minor non-resonant background processes. The systematic uncertainties are introduced as nuisance parameters in the fit. For each systematic uncertainty affecting the shape (normalisation) of the templates, a nuisance parameter is constrained with a Gaussian (log-normal) prior. The best fit values for all the nuisance parameters, as well as the corresponding uncertainties, are extracted by performing a binned maximum likelihood fit to the data. Figure 5 shows final distributions of events after a background-only fit in bins of ρ under two different mass hypotheses for the µµ + ee and eµ categories with all the nuisance parameters set to their best fit values. The corresponding signals are also displayed and normalised to 1 pb for illustrative purposes.
No significant deviations from the SM expectations are observed. The highest asymptotic local significance observed corresponds to 3.9 standard deviations for the signal hypothesis with m H = 627 GeV and m A = 162 GeV, which globally becomes 1.3 standard deviations once accounting for the look elsewhere effect [51], evaluated with the method described in Ref. [52]. The local p-value in the m H vs. m A plane is displayed in Fig. 6. Figure 7 shows expected (left) and observed (right) model independent upper limits at 95% confidence level (CL), σ 95% , on the product of the production cross section and branching fraction (σB) for H(A) → ZA(H) → bb, evaluated using the CL s criterion [53,54] in the asymptotic approximation [55] as a function of the H and A and mass hypotheses. Model dependent exclusion regions at 95% CL in the m H vs. m A plane can be obtained by comparing σ 95% to the theoretical cross section predicted by a particular model. Figure 8 shows the expected and observed 95% CL exclusion regions for the Type-II 2HDM benchmark scenario tan β = 1.5 and cos(β − α) = 0.01, while Fig. 9 shows the 95% CL exclusion region in the tan β vs. cos(β − α) plane for m H = 379 GeV and m A = 172 GeV.

Summary
This paper reports on a search for a new CP-even (odd) neutral Higgs boson, decaying into a Z boson and a lighter CP-odd (even) neutral Higgs boson, where the Z decays into an electron or muon pair, and the light Higgs boson into a b quark pair. The search is based on LHC protonproton collision data at a center-of-mass energy √ s = 13 TeV collected by the CMS experiment during 2016, corresponding to an integrated luminosity of 35.9 fb −1 . We consider decays such as H → ZA → bb, where H and A are the additional CP-even and -odd Higgs bosons above-mentioned, respectively, in the context of the two-Higgs-doublet model (2HDM). They are searched for in the mass range from 120 to 1000 GeV for H and 30 to 1000 GeV for A. The search is subsequently extended to the A → ZH → bb process via interchanging the two mass parameters.        [15] CMS Collaboration, "Search for a heavy pseudoscalar boson decaying to a Z and a Higgs boson at √ s = 13 TeV", Eur. Phys. J. C 79 (2019) 564, doi:10.1140/epjc/s10052-019-7058-z, arXiv:1903.00941.