Search for a pseudoscalar boson decaying into a Z boson and the 125 GeV Higgs boson in llbb final states

Results are reported on a search for decays of a pseudoscalar A boson into a Z boson and a light scalar h boson, where the Z boson decays into a pair of oppositely-charged electrons or muons, and the h boson decays into b anti-b. The search is based on data from proton-proton collisions at a center-of-mass energy sqrt(s)=8 TeV collected with the CMS detector, corresponding to an integrated luminosity of 19.7 inverse femtobarns. The h boson is assumed to be the standard model-like Higgs boson with a mass of 125 GeV. With no evidence for signal, upper limits are obtained on the product of the production cross section and the branching fraction of the A boson in the Zh channel. Results are also interpreted in the context of two Higgs doublet models.


Introduction
The discovery of a scalar boson at the CERN LHC [1-3] with properties in agreement with those predicted by the standard model (SM) raises the question of whether the Higgs sector consists of only one physical state, as expected in the SM, or whether additional bosons are also involved.
An extension of the SM Higgs sector is provided in two Higgs doublet models (2HDM) [4], which introduce a second scalar doublet in addition to the one from the SM. Different formulations of 2HDM predict different couplings of the two doublets to quarks and to leptons. In Type-I 2HDM, all fermions couple to only one Higgs doublet, while in Type-II, up-and down-type quarks couple to different doublets. One example of a Type-II 2HDM is the minimal supersymmetric standard model, despite that supersymmetry is not explicitly required in 2HDM. The second Higgs doublet entails the presence of five physical states: two neutral, CP-even states, representing a light h and a heavy H boson; a neutral, CP-odd A boson; and two charged scalar H ± bosons. The lightest scalar h is assumed to be the boson observed at the LHC at a mass of 125 GeV [5][6][7]. If the masses of the heavier bosons are at or below the TeV scale, they can be accessible at the LHC. Searches for this extended sector can be performed either by measuring the values of the couplings of the discovered h boson to other SM particles [8][9][10], or via direct searches in final states disfavored by the SM [11]. A way to probe this kind of new physics is therefore to search for bosons that decay into final states that contain an SM-like Higgs boson. This paper describes a search for a heavy pseudoscalar A boson that decays into a Z and an h boson, both on-shell, with the Z boson decaying into a pair of + − leptons ( being e or µ), and the h boson into bb. In most 2HDM formulations [4], the A boson is produced predominantly through gluon-gluon fusion and decays to on-shell Z and h bosons, provided that the mass of the A boson satisfies m A m h + m Z ≈ 220 GeV. This channel is expected to be viable for m A smaller than twice the top quark mass (m t ), where the decay A → Zh is generally dominant, but, depending on model parameters, it can also be sensitive at larger values of m A . The h → bb decay has a large branching fraction for most of the parameter space in 2HDM [9]. A similar analysis has been recently published by ATLAS [12].
The analysis strategy is to reconstruct the Z, h, and A boson candidates from the visible decay products in the event. The signal would manifest itself as a peak in the four-body invariant mass (m bb ) spectrum over an expected SM continuum. Irreducible backgrounds correspond to Z boson production with two accompanying b quark jets (genuine or mistagged), and tt events in the dileptonic final state. These backgrounds are evaluated and normalized directly using appropriate control regions in data. The h boson produced in association with a Z boson provides a contribution to the background, but it differs from signal because the m bb mass does not contain a resonant peak. Signal sensitivity is improved by exploiting the known value of the h boson mass, using it to rescale the jet momenta to match the value expected for the dijet invariant mass. In addition, optimal signal efficiency and background rejection is achieved using a multivariate discriminator. Results are extracted through a two-dimensional (2D) fit to m bb and the discriminator output; upper limits are presented on the product of the total cross section and the A → Zh, Z → , and h → bb branching fractions for a pseudoscalar boson, and interpreted within the 2HDM.

Data and simulation 2 CMS detector
A detailed description of the CMS detector, together with a definition of the coordinate system and kinematic variables, can be found in Ref. [13]. The central feature of the CMS apparatus is a superconducting solenoid of 6 m internal diameter, providing a magnetic field of 3.8 T. The field volume contains a silicon pixel and strip tracker, a homogeneous electromagnetic calorimeter (ECAL), and a sampling hadron calorimeter (HCAL), each composed of a barrel and two endcap sections. Muons are measured in gas-ionization detectors embedded in the steel flux-return yoke outside the solenoid.
The silicon tracker measures charged particles within the pseudorapidity range of |η| < 2.5. For non-isolated particles with transverse momenta of 1 < p T < 10 GeV and |η| < 1.4, the track resolutions are typically 1.5% in p T , and between 25-90 and 45-150 µm, respectively, in transverse and longitudinal impact parameters relative to the production vertex [14]. The ECAL consists of lead tungstate crystals that provide a coverage up to |η| < 3.0. The mass resolution for Z → e + e − decays when both electrons are in the ECAL barrel is 1.6%, and is 2.6% when both electrons are in the endcaps. The HCAL has alternating layers of brass as absorber and plastic scintillators, and covers the range of |η| < 3.0, which is extended to |η| 5.2 through forward calorimetry. Muons are measured in the range of |η| < 2.4, with detection planes made using three technologies: drift tubes, cathode-strip chambers, and resistive-plate chambers. Matching muons to tracks measured in the silicon tracker provides a p T resolution of 1.3-2.0% for muons with 20 < p T < 100 GeV in the barrel and better than 6% in the endcaps.

Data and simulation
Data used for this analysis were collected using double-muon and double-electron triggers, with p T thresholds set to 17 and 8 GeV for the highest and next highest p T lepton, respectively, and an isolation requirement used in the electron trigger to maintain an acceptable rate. The analyzed events correspond to an integrated luminosity of 19.7 ± 0.5 fb −1 of pp collisions at √ s = 8 TeV [15].
The SM background processes Z+jets, W+jets, and tt+jets or a vector boson (ttV) are simulated using the MADGRAPH 5.1 [16] Monte Carlo (MC) generator; multijets and dibosons (VV , with V, V being W or Z) are generated using LO PYTHIA 6.4 [17] and the CTEQ6L [18] parton distribution functions (PDF). Single top quark production, and SM Higgs boson production in association with an electroweak vector boson (Vh) are generated using the next-to-leading-order (NLO) POWHEG 1.0 [19][20][21] MC generator and the MSTW2008NLO PDF [22]. Parton showering and hadronization are performed with PYTHIA using the Z2* tune [23]. The generated MC events, including additional pp interactions (pileup) occurring in the bunch crossing containing the high-p T scatter, are processed through a full detector simulation based on GEANT4 [24] and reconstructed with the same algorithms as used for data.
Signal samples are generated with MADGRAPH 5.1. The pseudoscalar A boson is assumed to be produced via the gluon-gluon fusion process and to have a narrow width. The narrowwidth approximation is generally valid for m A 2m t , otherwise the width of the A boson strongly depends on the model parameters. The branching fraction B(A → Zh) is set to 100%, and similarly only the Z boson decays into electrons or muons, and the h boson decays into a pair of bb quarks, are considered. Other decay modes have negligible efficiency for passing the selections. The mass of the light Higgs boson is set to m h = 125 GeV, while the search for A boson is performed in the mass range from 225 to 600 GeV, where the efficiency to reconstruct two resolved jets from the h decay becomes too small.

Event reconstruction
A global reconstruction of the event is achieved using a particle-flow (PF) technique, which reconstructs and identifies individual particles emerging from each collision using information from all CMS subdetectors [25,26].
Electron candidates are reconstructed for |η| < 2.5 by matching energy depositions in the ECAL with reconstructed tracks [27]. The identification relies on a multivariate technique that combines observables sensitive to the amount of bremsstrahlung along the electron trajectory, the compatibility of the measured position, direction and momentum of the reconstructed track in the inner tracker, and the energy deposition reconstructed in the ECAL cluster. Additional requirements are imposed to remove electrons produced by photon conversions in the detector material.
Muons are reconstructed within |η| < 2.4, combining information from both the silicon tracker and the outer muon spectrometer [28], requiring small energy depositions in the calorimeters [29]. Muon candidates have to fulfill restrictive selection criteria based on the quality and the impact parameter of the track, as well as on the number of hits observed in the tracker and muon systems.
An isolation variable is defined for each lepton through the scalar sum over the p T of all PF candidates, excluding the lepton, within a cone of ∆R = √ (∆η) 2 + (∆ϕ) 2 < 0.4 around the lepton direction, where ϕ is the azimuthal angle in radians, subtracting the pileup contribution [30], and then dividing by the lepton p T . The lepton is rejected if the isolation exceeds 0.15 for electrons and 0.12 for muons.
Jets are formed from all particles, charged and neutral, reconstructed through the PF algorithm, and clustered with the anti-k T algorithm [31, 32] with a distance parameter of 0.5. Jet energy corrections are determined from dijet and Z/γ+jet events, and applied to both data and simulation [33]. The imbalance in transverse momentum is calculated as the negative of the vectorial sum of transverse momenta of all the reconstructed particles, and its magnitude is denoted as E miss T . The combined secondary vertex (CSV) b tagging algorithm [34] is used to identify jets that originate from b quarks. This algorithm combines information from track impact parameters relative to the primary and secondary (displaced) vertices into a likelihood discriminant. Two standard working points are set on the discriminant, corresponding to restrictive (tight) and less-restrictive (loose) thresholds that provide on average 50% and 80% b tagging efficiencies, with respective light-jet misidentification probabilities of about 0.1% and 10%. The CSV distribution is corrected in simulation to take into account a difference at the percent level in algorithm performance for data and simulation.

Event selection
Events are required to have at least two electrons or two muons within the geometrical acceptance regions, and to satisfy the reconstruction, identification, and isolation requirements. The p T threshold is set to 20 GeV for the lepton with highest p T , and to 10 GeV for the lepton with next-highest p T . The Z boson candidate is formed from the two highest-p T , opposite-charge, same-flavour leptons, and must have an invariant mass larger than 50 GeV. In addition, at least two jets are required with p T > 20 GeV, within |η| < 2.4, and with angular separation relative to each lepton of ∆R > 0.5. The h boson candidate is formed from the two jets that have the highest values of the b tagging discriminant. Additional jets in the event are ignored.
Signal events are expected to have at least one jet b-tagged with the tight and one with the loose CSV working points; after b tagging application, the contribution from pileup jets becomes negligible. The invariant masses of the two leptons and of the two jets are required to be compatible with those of the Z boson (75 < m < 105 GeV) and h boson (90 < m bb < 140 GeV), and the value of the E miss T in the event has to be compatible with zero.
Dedicated control regions are defined to check both the normalizations and distributions of the most important backgrounds by inverting the selections used to enhance signal. Drell-Yan backgrounds (from qq → Z/γ → + − production) are considered separately as a function of the number of b jets, distinguishing Z+jets (no b jets), Z + b (1 b jet) and Z + bb (2 b jets). The corresponding control regions are selected by requiring 80 < m < 100 GeV, E miss T < 40 GeV, vetoing of dijet masses close to the Higgs boson mass of 90 < m jj < 140 GeV, and applying different b-tagging selections. In the Z+jets control region, no cutoffs are applied on b tagging discriminators. The Z + b control region contains events with one b jet fulfilling the tight CSV working point, and no other jets passing the loose threshold. Events that enter the Z + bb control region should have at least two b-tagged jets, with one passing the tight and the other the loose working point. The tt control region is defined by inverting the m and E miss T selections, dropping the requirement of the dijet mass being not close to the Higgs boson mass, and requiring at least one tight and one loose b-tagged jet.
The scale factors that are used to correct the normalization of Drell-Yan and tt backgrounds, reported in Table 1, are obtained from a simultaneous likelihood fit to data and simulation in the four control regions and are applied in the following steps of the analysis. Multijet contamination in control and signal regions is evaluated with data by inverting lepton isolation criteria, and its contribution is found to be negligible. The yield of the Vh, ttV, and single top production through the s channel are calculated at NLO [35,36], while the other single top channels and dibosons are normalized to the measured cross sections [37-39]. A mismodelling in simulation, observed in the control regions relative to data, is corrected by reweighting with a linear function the event centrality distribution (defined as the ratio of the sums in scalar p T and the energy of the two leptons and two jets in the rest frame of the four objects) in all Monte Carlo events, improving the agreement between simulation and data for all the variables.

Background
Scale factor Z + jets 1.069 ± 0.002 An important feature of the signal is that the two b jets originate from the decay of the h boson, whose mass is known with better precision than that provided by the bb invariant mass resolution [40]. The measured jet p T , η, and ϕ values are therefore varied according to their resolution, in a kinematic fit based on Lagrange multipliers, to constrain the dijet invariant mass to m h = 125 GeV. The fit χ 2 is used in subsequent steps of the analysis as a discriminant in place of m bb . The kinematic fit improves the relative four-body invariant mass resolution from 6.3% to 1.2% and 4.0% to 1.9%, respectively, for the smallest and largest values of m A , centering the peaks around their nominal values, as shown in Fig. 1. The effect of the kinematic fit is larger at low m A , where the constraint on the Higgs invariant mass has the largest contribution. Although both the background and signal m bb distributions are modified by the kinematic fit, the signal significance in a mass window close to the investigated A boson mass increases by a factor of two at the lowest mass and by 34% at the highest mass. The resulting jet three-momenta are used to redefine all the kinematic variables in the event.
[GeV]  [41] are trained separately, one for each region. The inputs of each BDT consist of 16 discriminating variables, selected from a list of more than 40 variables: the p T of the Z and h boson candidates, the χ 2 of the kinematic fit, the significance of E miss T [42], the dilepton invariant mass, the ∆R separation and the "twist" angle (defined as tan −1 ∆ϕ/∆η) between the two b jets [43], their CSV discriminator values, the flight directions of the Z boson and of the beam in the rest frame of the A boson (cos θ * ), the decay angle of the Z boson relative to its flight direction in the rest frame of the Z boson (cos θ 1 ), which is sensitive to the transverse polarization of the Z boson along its flight direction, the angle of the pull vector [44] of the highest-p T jet, which exploits the colour connection between the two b quarks originating from the h boson, the scalar sum of E miss T and the p T of jets and leptons in the event (S T ), the number of jets with p T > 20 GeV, the event centrality, and aplanarity [43]. The distribution in BDT outputs for data, for simulated signal (S), and for the expected SM background (B) events are shown in Fig. 2, weighting each entry in the m bb distribution by the expected S/(S+B) ratio of the BDT bin in the signal-sensitive region with BDT > 0.6.

Systematic uncertainties
The sensitivity of the analysis is currently limited by the available data, and not by systematic uncertainties.
The uncertainties in the normalizations of the four main backgrounds (Z+jets, Z + b, Z + bb, and tt) originate both from the fits in the control regions and from the extrapolation to the signal region. The former are reported in Table 1, and the latter is evaluated through a simultaneous likelihood fit to data and to simulated yields in several statistically independent regions, obtained by altering the selections used to define the four control regions. A 13% normalization uncertainty due to this extrapolation is assigned to Z+jets, 12% to Z + b, 2.1% to Z + bb, and 6.2% to tt events. Normalization uncertainties for other SM backgrounds correspond to the ones on their measured or theoretical cross sections.
The uncertainties in lepton reconstruction, identification, isolation, and trigger efficiencies are evaluated through specific studies of events with masses in the region of the Z peak. Uncertainties in background or signal normalization and distribution due to uncertainties in jet energy scale and resolution [33] and b tagging scale factors [34] are estimated by changing the corresponding values by ±1 standard deviation (σ). Additional systematic uncertainties affecting the normalization of backgrounds and signal from the choice of PDF [45,46], contributions from pileup, E miss T fluctuations because of the presence of unclustered energy in the event, and integrated luminosity [15] are also considered in the analysis and are reported in Table 2 (Normalization).
Results are extracted from an analysis based on a binned likelihood fit to the two-dimensional (2D) distribution of m bb versus BDT. Dependence on jet energy scale and resolution, b tagging, as well as factorization and renormalization scales are propagated to the 2D-templates, taking into account the correlations between the two variables. The impact of reweighting the Monte Carlo distributions is considered as an additional source of background in modelling the uncertainty. Finally, the uncertainty from the limited number of simulated events is treated as in Ref. [47]. The sources of systematic uncertainty affecting the forms of the distributions are summarized in the first column of Table 2 (Shape). The second and third columns indicate the respective ranges in the relative impact in percent, obtained by the changes implemented in the background and signal contributions.
The systematic uncertainty with the largest impact on the expected limit is from the reweighting of the background (accounting for a 6% difference on the expected limit, depending on m A ), which is followed by the limitations in the number of simulated events ( 4%), and the factorization and renormalization scales ( 2%). The effect of the other sources is small ( 1%).

Results and their interpretation
Results are obtained from the combined signal and background fit to the binned two-dimensional distribution of the four-body invariant mass m bb and the BDT output in the signal-sensitive (BDT >0.6) region. With no evidence of significant deviation from background expectations, the asymptotic modified frequentist method is used to determine the limit at the 95% confidence level (CL) on the contribution from signal, treating systematic uncertainties as nuisance parameters that are integrated over in the fit [48][49][50][51].
The observed limit, as well as the expected limit and its relative ±1 and ±2σ bands of uncertainty, are reported as a function of the A boson mass in Fig. 3 for σ A B(A → Zh → bb), i.e. the product of the cross section and the A → Zh, h → bb, and Z → branching fractions, with = e or µ. The limits are obtained by considering the A boson produced via the gluon-gluon fusion process in the narrow-width approximation. Interpolated mass points are obtained as in Ref. [52], and numerical values are reported in Table 3. A signal upper limit at 95% CL is set on σ A B(A → Zh → bb), excluding from 10 to 30 fb for m A near the kinematic threshold, ≈8 fb for m A ≈ 2m t , and up to ≈3 fb at the high end (600 GeV) of the considered mass range. Comparable limits have been recently obtained for the same channel by the ATLAS Collaboration [12]. The most significant excess at m bb = 560 GeV has a local significance of 2.6σ, which For m A > 2m t , the width of the A boson depends strongly on the model parameters. Different limits are provided by taking into account the natural width of the A boson (Γ A ) in the reconstructed m bb , leaving the BDT unchanged. Figure 4 shows the exclusion limit above 2m t for an average width of 30 GeV, and the dependence of the observed limit on Γ A for m A = 500 GeV. The local and global significance at m bb = 560 GeV become, respectively, 2.9σ and 1.5σ, assuming Γ A = 30 GeV.
As an independent cross-check of the 2D fit, the signal is extracted by applying two complementary strategies, based on one-dimensional fits. The first one consists of fitting the m bb distribution, after selecting events in a signal-enriched region by applying a BDT > 0.8 selection. The second relies on fits to the BDT distributions, after selecting events within the resolution of the signal m bb peak. The two methods give upper limits compatible with those from the 2D fit, but 10 to 20% less stringent.
The results are interpreted in terms of Type-I and Type-II 2HDM formulations [4]. The B(A → Zh) and B(h → bb) branching fractions and the signal cross sections are computed at next-tonext-to-leading-order (NNLO) with SUSHI 1.2.0 [54] and 2HDMC 1.6.4 [55], respectively, using the MSTW2008LO, NLO, and NNLO sets of PDF. The B(Z → ) branching fraction, with = e or µ, is taken from the measured value [56]. Both gluon-gluon fusion and associated production with b quarks have been considered. The latter is rescaled to the fusion process, taking account of the difference in acceptance for signal, as well as the efficiency for selecting dijet pairs in the presence of combinatorial contributions from additional b quarks in the event. The parameters used for the models are: m h = 125 GeV, m H = m H ± = m A , m 2 12 = m 2 A [tan β]/[1 + tan 2 β], λ 6,7 = 0 [57], while m A = 225-600 GeV, 0.1 ≤ tan β ≤ 100, and −1 ≤ cos(β − α) ≤ 1, using the convention 0 < β − α < π, where tan β and α are, respectively, the ratio of the vacuum expectation values, and the mixing angle of the two Higgs doublets [4].
The observed limit, together with the expected limit and its relative ±1 and ±2σ uncertainty bands, are shown in Fig. 5, interpreted in Type-I and Type-II 2HDM for m A = 300 GeV. A sizeable fraction of the 2HDM phase space is excluded at a 95% CL with respect to the previous CMS searches [11].    : Observed and expected (together with ±1, 2σ uncertainty bands) exclusion limit for Type-I (left) and Type-II (right) models, as a function of tan β and cos(β − α). Contours are derived from the projection on the 2HDM parameter space for the m A = 300 GeV signal hypothesis; the observed limit is close to 1σ above the expected limit, as shown in Fig. 3.
in the context of Type-I and Type-II 2HDM formulations, thereby reducing the parameter space for extensions of the standard model.   [27] CMS Collaboration, "Performance of electron reconstruction and selection with the CMS detector in proton-proton collisions at √ s = 8 TeV", (2015). arXiv:1502.02701. Submitted to J. Instrum.