Search for a heavy pseudoscalar boson decaying to a Z and a Higgs boson at s=13Te\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{s}=13\,\text {Te}\text {V} $$\end{document}

A search is presented for a heavy pseudoscalar boson A\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {A}$$\end{document} decaying to a Z boson and a Higgs boson with mass of 125GeV\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\,\text {GeV}$$\end{document}. In the final state considered, the Higgs boson decays to a bottom quark and antiquark, and the Z boson decays either into a pair of electrons, muons, or neutrinos. The analysis is performed using a data sample corresponding to an integrated luminosity of 35.9fb-1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\,\text {fb}^{-1}$$\end{document} collected in 2016 by the CMS experiment at the LHC from proton–proton collisions at a center-of-mass energy of 13Te\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\,\text {Te}\text {V}$$\end{document}. The data are found to be consistent with the background expectations. Exclusion limits are set in the context of two-Higgs-doublet models in the A\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {A}$$\end{document} boson mass range between 225 and 1000GeV\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\,\text {GeV}$$\end{document}.


Introduction
The discovery of a Higgs boson at the CERN LHC [1,2] and the measurement of its mass, spin, parity, and couplings [3,4] raises the question of whether the Higgs boson sector consists of only one scalar doublet, which results in a single physical Higgs boson as expected in the standard model (SM), or whether additional bosons are involved in electroweak (EW) symmetry breaking.
The two-Higgs-doublet model (2HDM) [5] provides an extension of the SM Higgs boson sector introducing a second scalar doublet. The 2HDM is incorporated in supersymmetric models [6], axion models [7], and may introduce additional sources of explicit or spontaneous CP violation that explain the baryon asymmetry of the universe [8]. Various formulations of the 2HDM predict different couplings of the two doublets to right-handed quarks and charged leptons: 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 fermions e-mail: cms-publication-committee-chair@cern.ch and leptons couple to one of the Higgs doublets, while the down-type quarks couple to the other.
The two Higgs doublets entail the presence of five physical states: two neutral and CP-even bosons (h and H, the latter being more massive), a neutral and CP-odd boson (A), and two charged scalar bosons (H ± ). The model has two free parameters, α and tan β, which are the mixing angle and the ratio of the vacuum expectation values of the two Higgs doublets, respectively. If tan β 5, the dominant A boson production process is via gluon-gluon fusion, otherwise associated production with a b quark-antiquark pair becomes significant. The diagrams of the two production modes are shown in Fig. 1. At small tan β values the heavy pseudoscalar boson A may decay with a large branching fraction to a Z and an h boson, if kinematically allowed [5]. These models can be probed either with indirect searches, by measuring the cross section and couplings of the SM Higgs boson [9], or by performing a direct search for an A boson. This paper describes a search for a heavy pseudoscalar A boson that decays to a Z and an h boson, both on-shell, with the Z boson decaying to + − ( being an electron or a muon) or to a pair of neutrinos, and the h boson to bb. The h boson is assumed to be the 125 GeV boson discovered at the LHC. In this search, the candidate A boson is reconstructed from the invariant mass of the visible decay products in events when the Z boson decays to charged leptons, or is inferred through a partial reconstruction of the mass using quantities measured in the transverse plane when the Z boson decays to neutrinos. The signal would emerge as a peak above the SM continuum of the four-body invariant mass (m Zh ) spectrum for the former decay mode and the transverse mass (m T Zh ) for the latter. The signal sensitivity is maximized by exploiting the known value of the h boson mass to rescale the jet momenta and significantly improve the m Zh resolution. In addition, selections based on multivariate discriminators, exploiting event variables such as angular distributions, are used to optimize the signal efficiency and background rejection. This search is particularly sensitive to a pseudoscalar A boson with a mass smaller than twice the top quark mass and for small tan β values. In this region of the 2HDM parameter space, the A boson cross section is larger than 1 pb, and the A boson decays predominantly to Zh [5].
With respect to the CMS search performed at √ s = 8 TeV [10], this analysis benefits from the increased centerof-mass energy and integrated luminosity, includes final states with invisible decays of the Z boson, increases the sensitivity to b quark associated production, and extends the A boson mass (m A ) range from 600 to 1000 GeV. At larger m A , the angular separation between the b quarks becomes small, and the Higgs boson is reconstructed as a single large-cone jet; the corresponding CMS analysis presents limits on the 2HDM from 800 GeV to 2 TeV [11]. The ATLAS Collaboration has published a search probing Zh resonances with similar event selections based on a comparable data set, observing a mild excess near 440 GeV in categories with additional b quarks [12].

The CMS detector
A 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. [13].
The central feature of the CMS apparatus is a superconducting solenoid of 6 m internal diameter, providing a mag-netic 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 (HCAL), 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 1440 silicon pixel and 15,148 silicon strip detector modules. For nonisolated particles with transverse momenta of 1 < p T < 10 GeV and |η| < 1.4, the track resolutions are typically 1.5% in p T and 25-90 (45-150) µm in the transverse (longitudinal) impact parameter [14]. The ECAL provides coverage up to |η| < 3.0, and the energy resolution for unconverted or lateconverting electrons and photons in the barrel section is about 1% for particles that have energies in the range of tens of GeV. The dielectron mass resolution for Z → e + e − decays when both electrons are in the ECAL barrel is 1.9%, and is 2.9% when both electrons are in the endcaps [15]. The muon detectors covering the range |η| < 2.4 make use of three different technologies: drift tubes, cathode strip chambers, and resistive-plate chambers. Combining muon tracks with matching tracks measured in the silicon tracker results in a p T resolution of 2-10% for muons with 0.1 < p T < 1 TeV [16].
The first level of the CMS trigger system [17], composed of custom hardware processors, uses information from the calorimeters and muon detectors to select the most interesting events in a fixed time interval of less than 4 µs. The high-level trigger (HLT) processor farm decreases the event rate from around 100 kHz to about 1 kHz, before data storage.

Event reconstruction
A global event reconstruction is performed with a particleflow (PF) algorithm [18], which uses an optimized combination of information from the various elements of the detector to identify stable particles reconstructed in the detector as an electron, a muon, a photon, a charged or a neutral hadron. The PF particles have to pass the charged-hadron subtraction (CHS) algorithm [19], which discards charged hadrons not originating from the primary vertex, depending on the longitudinal impact parameter of the track. The primary vertex is selected as the vertex with the largest value of summed p 2 T of the PF particles, including charged leptons, neutral and charged hadrons clustered in jets, and the associated missing transverse momentum p miss T , which is the negative vector sum of the p T of those jets.
Electrons are reconstructed in the fiducial region |η| < 2.5 by matching the energy deposits in the ECAL with charged particle trajectories reconstructed in the tracker [15]. The electron identification is based on the distribution of energy deposited along the electron trajectory, the direction and momentum of the track, and its compatibility with the primary vertex of the event. Electrons are further required to be isolated from other energy deposits in the detector. The electron relative isolation parameter is defined as the sum of transverse momenta of all the PF candidates, excluding the electron itself, divided by the electron p T . The PF candidates are considered if they lie within ΔR = √ (Δη) 2 + (Δφ) 2 < 0.3 around the electron direction, where φ is the azimuthal angle in radians, and after the contributions from pileup and other reconstructed electrons are removed [15].
Muons are reconstructed within the acceptance of the CMS muon systems using tracks reconstructed in both the muon spectrometer and the silicon tracker [16]. Additional requirements are based on the compatibility of the trajectory with the primary vertex, and on the number of hits observed in the tracker and muon systems. Similarly to electrons, muons are required to be isolated. The muon isolation is computed from reconstructed PF candidates within a cone of ΔR < 0.4 around the muon direction, ignoring the candidate muon, and divided by the muon p T [16].
Hadronically decaying τ leptons are used to reject W → τ ν background events, and are reconstructed by combining one or three hadronic charged PF candidates with up to two neutral pions, the latter also reconstructed by the PF algorithm from the photons arising from the π 0 → γ γ decay [20].
Jets are clustered using the anti-k T algorithm [21,22] with a distance parameter of 0.4. The contribution of neutral particles originating from pileup interactions is estimated to be proportional to the jet area derived using the FastJet package [22,23], and subtracted from the jet energy. Jet energy corrections, extracted from both simulation and data in multijet, γ +jets, and Z+jets events, are applied as functions of the p T and η of the jet to correct the jet response and to account for residual differences between data and simulation. The jet energy resolution amounts typically to 15-20% at 30 GeV, 10% at 100 GeV, and 5% at 1 TeV [24].
Jets that originate from b quarks are identified with a combined secondary vertex b-tagging algorithm [25] that uses the tracks and secondary vertices associated with the jets as inputs to a neural network. The algorithm provides a b jet tagging efficiency of 70%, and a misidentification rate in a sample of quark and gluon jets of about 1%. The b tagging efficiency is corrected to take into account a difference at the few percent level in algorithm performance for data and simulation [25].

Data and simulated samples
The data sample analyzed in this search corresponds to an integrated luminosity of 35.9 fb −1 of proton-proton (pp) collisions at a center-of-mass energy of 13 TeV collected with the CMS detector at the LHC. Data are collected using triggers that require either the presence of at least one isolated electron or isolated muon with p T > 27 GeV, or alternatively a p miss T or H miss T larger than 90-110 GeV, the value depending on the instantaneous luminosity. The p miss T is the magnitude of p miss T , and H miss T is defined as the momentum imbalance of the jets in the transverse plane [17].
The pseudoscalar boson signal is simulated at leading order (LO) with the MadGraph5_amc@nlo 2.2.2 matrix element generator [26] in both the gluon-gluon fusion and b quark associated production modes according to the 2HDM [5], assuming a narrow signal width. The h boson mass is set to 125 GeV, and the A boson mass ranges between 225 and 1000 GeV. The A → Zh decay is simulated with MadSpin [27]. The Higgs boson is forced to decay to bb, and the vector boson to a pair of electrons, muons, τ leptons, or neutrinos. In the gluon-gluon fusion production mode, up to one additional jet is included in matrix element calculations, and only the top quark contributes to the loop shown in Fig. 1 (upper). The 2HDM cross sections and branching fractions are computed at next-to-next-to-leading order (NNLO) with 2hdmc 1.7.0 [28] and SusHi 1.6.1 [29], respectively. The parameters used in the models are: m h = 125 GeV, m H = m H ± = m A , the discrete Z 2 symmetry is broken as in the minimal supersymmetric standard model (MSSM), and CP is conserved at tree level in the 2HDM Higgs sector [5]. The branching fractions of the Z boson are taken from the measured values [30].
The SM backgrounds in this search consist of the inclusive production of a vector boson in association with other jets (V+jets, with V = W or Z, and V decaying to final states with charged leptons and neutrinos), and top quark pair production (tt). V+jets events are simulated at LO with MadGraph5_amc@nlo with up to four partons included in the matrix element calculations and using the MLM matching scheme [31]. The event yield is normalized to the NNLO cross section computed with fewz v3.1 [32]. The V boson p T spectra are corrected to account for next-to-leading order (NLO) quantum chromodynamics (QCD) and EW contributions [33]. The tt and single top quark in the t channel and tW production are simulated at NLO with powheg v2 generator [34][35][36]. The number of events for the top quark pair production process is rescaled according to the cross section computed with Top++ v2.0 [37] at NNLO+NNLL, and the transverse momenta of top quarks are corrected to match the distribution observed in data [38]. Other SM processes, such as SM vector boson pair production (VV), SM Higgs boson production in association with a vector boson (Vh), single top quark (t + X) production in the s channel, and top quark production in association with vector bosons, are simulated at NLO in QCD with MadGraph5_amc@nlo using the FxFx merging scheme [39]. The multijet contribution, estimated with the use of samples generated at LO with the same generator, is negligible after analysis selections.
All the simulated processes use the NNPDF 3.0 [40] parton distribution functions (PDFs), and are interfaced with pythia 8.205 [41,42] for the parton showering and hadronization. The CUETP8M1 underlying event tune [43] is used in all samples, except for top quark pair production, which adopts the CUETP8M2T4 tune [44].
Additional minimum bias pp interactions within the same or adjacent bunch crossings (pileup) are added to the simulated processes, and events are weighted to match the observed average number of interactions per bunch crossing. Generated events are processed through a full CMS detector simulation based on Geant4 [45] and reconstructed with the same algorithms used for collision data.

Event selection
Events are classified into three independent categories (0 , 2e, and 2μ), based on the number and flavor of the reconstructed leptons. Events are required to have at least two jets with p T > 30 GeV and |η| < 2.4 to be suitable candidates for the reconstruction of the h → bb decay. If more than two jets fulfill the requirements, the ones with the largest b tagging discriminator value are used to reconstruct the Higgs boson candidate. The efficiency of the correct assignment of the reconstructed jets to initial quarks originating from the Higgs boson decay varies between 80 and 97%, after applying the event selections, depending on the category and final state.
In the 0 category, no isolated electron or muon with p T > 10 GeV is allowed. Events containing isolated hadronic decays of the τ leptons with p T > 18 GeV are vetoed as well. A selection is applied on the reconstructed p miss T , which is required to be larger than 200 GeV, such that the p miss T trigger is at least 95% efficient. In order to select a topology where the Z boson recoils against the Higgs boson, a Lorentz boost requirement of 200 GeV on the p T of the Higgs boson candidate, p bb T , is applied. Multijet production is suppressed by requiring that the minimum azimuthal angular separation between all jets and the missing transverse momentum vector must satisfy Δφ(jet, p miss T ) > 0.4. The multijet simulation is validated in a region obtained by inverting the Δφ selection, finding a good description of data. When the Z boson decays to neutrinos, the resonance mass m A cannot be reconstructed directly. In this case, m A is estimated by computing the transverse mass from the p miss T and the four-momenta of the two jets used to reconstruct the Higgs boson candidate, defined , which has to be larger than 500 GeV. The efficiency of these selections for signal events with m A 500 GeV is small, because the p T of the Z boson is not sufficient to produce a p miss T large enough to pass the selection; thus, the contribution of the 0 category is significant only for large m A .
In the 2e and 2μ categories, events are required to have at least two isolated electrons or muons within the detector geometrical acceptance. The p T threshold on the lepton is referred to as p T , and is set to 30 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-flavor leptons, and must have an invariant mass m between 70 and 110 GeV. The m selection lowers the contamination from tt dileptonic decays, and significantly reduces the contribution from Z → τ τ decays. The reconstructed p miss T also has to be smaller than 100 GeV to reject the tt background. In order to maximize the signal acceptance, no Lorentz boost requirement is applied to the Z and h boson candidates in the dileptonic categories. The A boson candidate is reconstructed from the invariant mass m Zh of the Z and h boson candidates.
If the two jets originate from a Higgs boson, their invariant mass is expected to peak close to 125 GeV. Events with a dijet invariant mass m jj between 100 and 140 GeV enter the signal regions (SRs); otherwise, if m jj < 400 GeV, they fall in dijet mass sidebands, which are used as control regions (CRs) to estimate the contributions of the main backgrounds. Signal regions are further divided by the number of jets passing the b tagging requirement (1, 2, or at least 3 b tags). The 3 b tag category has been defined to select the additional b quarks from b quark associated production. In this region, at least one additional jet, other than the two used to reconstruct the h boson, has to pass the kinematic selections and b tagging requirements. The fraction of signal events passing the m jj selection in the SR is 66-82% and 45-65% in the 1 and 2 b tag categories, respectively. Control regions for the Z+jets background share the same selections as the corresponding SR, except for the m jj mass window.
Dedicated CRs are defined to estimate the tt and W+jets backgrounds, which may enter the 0 SR if the lepton originating from the W decay is outside the detector geometrical acceptance or is not reconstructed. Two W+jets CRs share the same selection as in the 0 categories, but require exactly one electron or one muon passing the same trigger and selections of the leading lepton in the 2 categories. In order to mimic the kinematics of leptonic W decays, where the lepton is outside the geometrical acceptance or is not reconstructed in the detector, the p miss T is recalculated by removing the contribution of the lepton. The min(Δφ) requirement is removed, and the dijet invariant mass selection is not applied, as the signal is absent in 1 final states. Events are required to have three or fewer jets, none of them b tagged, to reduce the tt contribution.
Four different CRs associated with the production of events containing top quarks are defined by inverting specific selections with respect to the SR definition. Dileptonic tt control regions require the same selections as the 2e and 2μ categories with two b tags, but the dilepton invariant mass region around the nominal Z boson mass is vetoed (50 < m < 70 GeV or m > 110 GeV), and the m jj selection is dropped. Two additional top quark CRs are defined specifically for tt events where only one of the two W bosons decays into an electron or a muon, and the lepton is not reconstructed. These events contribute to the tt contamination in the 0 categories. The two single-lepton top quark CRs have the same selections as the two W+jets CRs, but in this case the jet and b tag vetoes are inverted to enrich the tt composition.
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. The measured jet p T values are therefore scaled according to their corresponding uncertainty given by the jet energy scale corrections to constrain the dijet invariant mass to m jj = 125 GeV. The kinematic constraint on the h boson mass improves the relative four-body invariant mass resolution from 5-6 to 2.5-4.5% for the smallest and largest values of m A , respectively. Similarly, in the 2 channels, the electron and muon p T are scaled to a dilepton invariant mass m = m Z . The effect on the m A resolution of the kinematic constraint on the leptons is much smaller than the one of the jets, because of their better momentum resolution.
In the 2e and 2μ categories, the A boson decay chain yields an additional characteristic, which helps distinguish it from SM background. Five helicity-dependent angular observables fully describe the kinematics of the A → Zh → bb decay: the angle between the directions of the Z boson and the beam in the rest frame of the A boson (cos θ * ); the decay angle between the direction of the negatively charged lepton relative to the Z boson momentum vector in the rest frame of the Z boson (cos θ 1 ), which is sensitive to the transverse polarization of the Z boson along its momentum vector; the angle between a jet from the h boson and the h boson momentum vector in the h boson rest frame (cos θ 2 ); the angle between the Z and h boson decay planes in the rest frame of the A boson (Φ); the angle between the h boson decay plane and the plane where the h boson and the beam directions lie in the A boson rest frame (Φ 1 ). The discriminating power and low cross-correlation make these angles suitable as input to a likelihood ratio multivariate discriminator. This angular discriminant is defined as: where the index i runs from 1 to 5 and corresponds to the number N of angular variables x i , and s i and b i are the signal and Z+jets background probability density functions of the i-th variable, respectively. A selection of D > 0.5 is applied in all 2e and 2μ SRs and CRs, except those with three b tags due to the low event count. This working point retains 80% of the signal efficiency and rejects 50% of the Z+jets background.
Considering that top quark pair production may be as large as 50% of the total background in certain regions of the parameter space, a second likelihood ratio discriminator is built specifically to reject the tt events. This discriminator uses only the m and p miss T variables. The background probability density function considers only the top quark background in order to achieve the maximum separation between events with a genuine leptonically decaying Z boson recoiling against a pair of jets and the more complex topologies such as tt decays. Selecting events with a discriminator output larger than 0.5 rejects 75% of the tt events with a signal efficiency of 85%. This selection is applied to the dileptonic SRs and to the Z+jets CRs.
The SRs and CRs selections are summarized in Table 1. The product of the signal acceptance and selection efficiency as a function of m A is presented in Fig. 2 separately for the gluon-gluon fusion and b quark associated production modes.

Systematic uncertainties
The uncertainties in the trigger efficiency and the electron, muon, and τ lepton reconstruction, identification, and isolation efficiencies are evaluated through studies of events with dilepton invariant mass around the Z boson mass, and the variation of the event yields with respect to the expectation from simulation amount to approximately 2-3% for the categories with charged leptons, and 1% in the 0 categories [15,16,20]. The impact of the lepton energy and momentum scale and resolution is small after the kinematic constraint on m . The jet energy scale and resolution [24] affect both the selection efficiencies and the shape of the p miss T and m T Zh distributions, and are negligible in the 2 channels after the kinematic constraint on the dijet mass has been applied. The jet four-momentum is varied by the corresponding uncertainties, and the effect is propagated to the final distributions. The jet energy scale is responsible for a 2-6% variation in the numbers of background and signal events; the jet energy resolution contributes an additional 1-2% uncertainty. The effects of jet energy scale and resolution b-tagged jets 1, 2, or 3 0, 1, 2, or 3 0 ≥ 1 1, 2, or 3 0, 1, 2, or 3 ≥ 2    [25] in the signal yield depends on the jet p T and thus on the mass of the resonance, and the impact on the event yield ranges from 2 to 4% in the 1 b tag category, 4 to 8% in the 2 b tag category, and 8 to 12% in the 3 b tag category.

CMS
The signal and background event yields are affected by the uncertainties on the choice of PDFs [46] and the factorization and renormalization scale uncertainties. The former are derived with SysCalc [47], and the latter are estimated by varying the corresponding scales up and down by a factor of two [48]. The effect of both these uncertainties can be as large as 30% depending on the generated signal mass. The effect of the PDF uncertainties on the signal and background lepton acceptance is estimated to be an average of 3% per lepton. The top quark background is also affected by the uncertainty associated with the simulated p T spectrum of top quarks [38], which results in up to a 14% yield uncertainty. The V+jets backgrounds are affected by the uncertainties on the QCD and EW NLO corrections, as described in Sect. 4.
A systematic uncertainty is assigned to the interpolation between the two mass sidebands to the SR, defined as the difference in the ratio between data and simulated background in the lower and upper sidebands, and ranges between 2 and 10% depending on the channel. The extrapolation to the 3 b tag regions is covered by a large uncertainty (20-46%) assigned to the overall background normalization, and derived by taking the ratio between data and the simulation in the 3 b tag control regions. In the dilepton categories, a dedicated uncertainty is introduced to cover for minor mismodeling effects. The background distribution is reweighted  with a linear function of the event centrality (defined as the ratio between the sums of the p T and the energy of the two leptons and two jets in the rest frame of the four objects) in all simulated events, and the effect is propagated to the m Zh distributions as a systematic uncertainty. Additional systematic uncertainties affect the event yields of backgrounds and signal come from pileup contributions and integrated luminosity [49]. The uncertainty from the limited number of simulated events is treated as in Ref. [50]. A summary of the systematic uncertainties is reported in Table 2.

Results and interpretation
The signal search is carried out by performing a combined signal and background maximum likelihood fit to the number of events in the CRs, and the binned m Zh or m T Zh distributions in the SRs. Systematic uncertainties are treated as nuisance parameters and are profiled in the statistical inter-   pretation [51][52][53]. The asymptotic approximation [54] of the modified frequentist CL s criterion [51,52] is used to determine limits on the signal cross section at 95% confidence level (CL). The background-only hypothesis is tested against the combined signal+background hypothesis in the nine cate- ) + jets  Table 4 Expected and observed event yields after the fit in the signal regions. The dielectron and dimuon categories are summed together. The "-" symbol represents backgrounds with no simulated events passing the selections. The signal yields refer to pre-fit values correspond-ing to a cross section multiplied by B(A → Zh) B(h → bb) of 0.1 pb (gluon-gluon fusion for m A = 300 GeV, and in association with b quarks for m A = 1000 GeV) Signal region 0 , 1 b tag 0 , 2 b tag 0 , ≥3 b tag 2 , 1 b tag 2 , 2 b tag 2 , ≥3 b tag    Table 3, and the overall event yields in the CRs are shown in Fig. 3 before and after the fit. The expected and observed number of events in the SRs are reported in Table 4, and the m Zh and m T Zh distributions are shown in Fig. 4.
The data are well described by the SM processes. Upper limits are derived on the product of the cross section for a heavy pseudoscalar boson A and the branching fractions for the decays A → Zh and h → bb. The limits are obtained by considering the A boson produced via the gluon-gluon fusion and b quark associated production processes separately, in the approximation where the natural width of the A boson Γ A is smaller than the experimental resolution, and are reported in Fig. 5. An upper limit at 95% CL on the number of signal events is set on σ A B(A → Zh) B(h → bb), excluding above 1 pb for m A near the kinematic threshold, ≈0.3 pb for m A ≈ 2m t , and as low as 0.02 pb at the high end (1000 GeV) of the considered mass range. The sensitivity of the analysis is limited by the amount of data, and not by systematic uncertainties. These results extend the search for a 2HDM pseudoscalar boson A for mass up to 1 TeV, which The excluded region is represented by the shaded gray area. The regions of the parameter space where the natural width of the A boson Γ A is comparable to the experimental resolution and thus the narrow width approximation is not valid are represented by the hatched gray areas is a kinematic region previously unexplored by CMS in the 8 TeV data analysis [10]. When m A is larger than 1 TeV, the CMS analysis with merged jets [11] retains a better sensitivity. The sensitivity is comparable to the ATLAS search [12], which observed a mild local (global) excess of 3.6 (2.4) standard deviations corresponding to m A ≈ 440 GeV in final states with 2μ and 3 or more b-tagged jets. A slight deficit is observed by CMS in the corresponding region.
The results are interpreted in terms of Type-I, Type-II, "lepton-specific", and "flipped" 2HDM formulations [5]. In the scenario with cos(β−α) = 0.1 and tan β = 3, an A boson up to 380 and 350 GeV is excluded in 2HDM Type-I and Type-II, respectively, as depicted in Fig. 5. These exclusion limits are used to constrain the two-dimensional plane of the 2HDM parameters [cos(β − α), tan β] as reported in Fig. 6, with fixed m A = 300 GeV in the range 0.1 ≤ tan β ≤ 100 and −1 ≤ cos(β − α) ≤ 1, using the convention 0 < β − α < π. Because of the suppressed A boson cross section and B(A → Zh), the region near cos(β − α)≈0 is not accessible in this search. On the other hand, B(h → bb) vanishes in the diagonal regions corresponding to α close to 0 in Type-II and flipped 2HDM, and α → ±π/2 in Type-I and leptonspecific scenarios. The exclusion as a function of m A , fixing cos(β − α) = 0.1, is also reported in Fig. 7.

Summary
A search is presented in the context of an extended Higgs boson sector for a heavy pseudoscalar boson A that decays into a Z boson and an h boson with mass of 125 GeV, with the Z boson decaying into electrons, muons, or neutrinos, and the h boson into bb. The SM backgrounds are suppressed by using the characteristics of the considered signal, namely the production and decay angles of the A, Z, and h bosons, and by improving the A mass resolution through a kinematic constraint on the reconstructed invariant mass of the h boson candidate. No excess of data over the background prediction is observed. Upper limits are set at 95% confi-dence level on the product of the A boson cross sections and the branching fractions σ A B(A → Zh) B(h → bb), which exclude 1 to 0.01 pb in the 225-1000 GeV mass range, and are comparable to the corresponding ATLAS search. Interpretations are given in the context of Type-I, Type-II, flipped, and lepton-specific two-Higgs-doublet model formulations, thereby reducing the allowed parameter space for extensions of the SM with respect to previous CMS searches.
tively the computing infrastructure essential to our analyses. Finally, we acknowledge the enduring support for the construction and operation of the LHC and the CMS detector provided by the following funding agencies: BMWFW

Data Availability Statement
This manuscript has no associated data or the data will not be deposited. [Authors' comment: Release and preservation of data used by the CMS Collaboration as the basis for publications is guided by the CMS policy as written in its document "CMS data preservation, re-use and open access policy" (https://cms-docdb.cern. ch/cgi-bin/PublicDocDB/RetrieveFile?docid=6032&filename=CMS DataPolicyV1.2.pdf&version=2).] Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecomm ons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Funded by SCOAP 3 .