Search for heavy Higgs bosons decaying to a top quark pair in proton-proton collisions at s\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \sqrt{s} $$\end{document} = 13 TeV

A search is presented for additional scalar (H) or pseudoscalar (A) Higgs bosons decaying to a top quark pair in proton-proton collisions at a center-of-mass energy of 13 TeV. The data set analyzed corresponds to an integrated luminosity of 35.9 fb−1 collected by the CMS experiment at the LHC. Final states with one or two charged leptons are considered. The invariant mass of the reconstructed top quark pair system and variables that are sensitive to the spin of the particles decaying into the top quark pair are used to search for signatures of the H or A bosons. The interference with the standard model top quark pair background is taken into account. A moderate signal-like deviation compatible with an A boson with a mass of 400 GeV is observed with a global significance of 1.9 standard deviations. New stringent constraints are reported on the strength of the coupling of the hypothetical bosons to the top quark, with the mass of the bosons ranging from 400 to 750 GeV and their total relative width from 0.5 to 25%. The results of the search are also interpreted in a minimal supersymmetric standard model scenario. Values of mA from 400 to 700 GeV are probed, and a region with values of tan β below 1.0 to 1.5, depending on mA, is excluded at 95% confidence level.


Introduction
The observation of a Higgs boson with a mass of 125 GeV by the ATLAS and CMS experiments [1][2][3] was a milestone in particle physics, confirming the existence of a crucial ingredient of the standard model (SM) of particle physics. Multiple extensions of the SM predict new spin-0 states. These include two-Higgs-doublet models (2HDMs) [4], of which the minimal supersymmetric standard model (MSSM) [5,6] is a particular realization, models predicting a new electroweak singlet [7], and other models with a combination of singlet and doublet fields [8]. The additional bosons may also provide a portal to dark matter, by acting as a mediator between SM and dark matter particles [9,10].
The new states introduced in these extensions of the SM may include charged Higgs bosons, H ± , scalar (CP-even) neutral H and h bosons (here h denotes the lighter of the two states), and a pseudoscalar (CP-odd) neutral A boson. For convenience and depending on the context, a common symbol Φ is used in this paper to represent the H and A bosons.
Top quarks play a key role in searches for new physics because of their high mass and large coupling to the SM Higgs boson. Provided that additional Higgs bosons couple to fermions via a Yukawa interaction, the top quark's high mass suggests the size of the coupling to these new bosons to be large as well. Hence, assuming the masses of the new Φ bosons are sufficiently high, their possible decay to a top quark pair is particularly interesting. Decays of CP-odd Higgs bosons to weak vector bosons, A → WW and A → ZZ, are forbidden (at tree level) if the CP symmetry is assumed. Such decays are also strongly suppressed for H bosons in the vicinity of the alignment limit of 2HDMs, in which the properties of the h boson approach those of the SM Higgs boson [11]. The aligned scenario may naturally occur as a result of a broken symmetry, and such models imply the existence of additional Higgs bosons at relatively low mass ( 550 GeV) [12]. In this paper, however, we do not rely on the assumption of alignment and its naturalness, and probe masses between 400 and 750 GeV. We consider a Yukawa-like coupling between the new spin-0 bosons and the top quark. The corresponding terms in the Lagrangian for the two CP eigenstates read as where m t is the top quark mass, v is the SM Higgs vacuum expectation value, and the strength of the couplings is controlled by real-valued coupling modifiers g Φtt 0.
A special case of a Type-II 2HDM [4] is the Higgs sector in the hMSSM [13], where the h boson is identified with the Higgs boson with mass of 125 GeV. The hMSSM can be fully described by two tree-level parameters: tan β, the ratio of the vacuum expectation values of the two Higgs fields, and m A , the mass of the pseudoscalar boson. The parameter region at low values of tan β is of particular interest, since the coupling of the additional Higgs bosons to top quarks is enhanced in this regime.
We consider the production of new Φ bosons through the gluon fusion process, with only top quarks in the loop. When the heavy Higgs boson decays into a top quark pair, this mode interferes at the quantum level with the SM production of top quark pairs. Example Feynman diagrams are shown in figure 1. As a consequence, the signal consists of a resonant and an interference component. The resonant component corresponds to the square of the amplitude given by the signal diagram, and results in a Breit-Wigner peak in the distribution of the invariant mass of the tt system, m tt . The interference component may be either destructive or constructive, depending on the phase space region and signal model. The sum of the components may result in a peak-dip structure in the m tt distribution [14][15][16]. It is worth noting that the shape and magnitude of the interference depends on the specific signal model, and can be significantly modified by new particles appearing in the loop of the production diagram [17,18].

JHEP04(2020)171
Decays of the scalar and pseudoscalar Higgs bosons produce tt pairs in the 3 P 0 and 1 S 0 states respectively [16], while the SM gg → tt production results in a mixture of states that changes with the partonic center-of-mass energy. Consequently, the signal and the background exhibit different angular properties, providing an additional handle to distinguish them.
A search for H or A bosons decaying to a top quark pair was performed at a centerof-mass energy √ s = 8 TeV by the ATLAS experiment [19]. The results were interpreted within the context of a Type-II 2HDM. The CMS experiment performed a search for top quark associated production of an H or A boson decaying to a top quark pair at √ s = 13 TeV [20]. The ATLAS and CMS collaborations also searched for spin-1 and spin-2 resonances decaying to a top quark pair [21,22], generally probing very high masses and Lorentz-boosted topologies, and without considering quantum interference with SM top quark pair production. In addition, both collaborations performed searches for H ± decaying to a top and a bottom quark [23,24], which are also sensitive to the region of low tan β in the hMSSM parameter space.
This paper describes a search for H or A bosons decaying to a top quark pair in proton-proton collisions at √ s = 13 TeV using the CMS detector at the CERN LHC. The data set analyzed corresponds to an integrated luminosity of about 35.9 fb −1 , collected in 2016. Events are selected in which the top quark pair decays into a final state with one or two leptons, where a lepton refers to an electron or a muon throughout this article. This analysis exploits both m tt and angular variables sensitive to the spin of the heavy Higgs bosons. Constraints on the coupling modifier g Φtt are derived as a function of the boson mass and width. The results are also interpreted in the hMSSM context, putting constraints in the (m A , tan β) plane.

The CMS detector and event reconstruction
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 (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. 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. [25].
The reconstructed vertex with the largest value of summed physics-object p 2 T is taken to be the primary pp interaction vertex. The physics objects are the jets, clustered using the jet finding algorithm [26,27] with the tracks assigned to the vertex as inputs, and the associated missing transverse momentum, taken as the negative vector sum of the p T of those jets.
The particle-flow (PF) algorithm [28] aims to reconstruct and identify each individual particle in an event, with an optimized combination of information from the various ele--3 -

JHEP04(2020)171
in the endcaps. The p T resolution in the barrel is better than 10% for muons with p T up to 1 TeV [33]. Simulation-to-data scale factors that depend on the lepton p T and η are used to correct for small differences in lepton trigger, identification, and isolation efficiency.
We define tight and loose collections of electron and muon candidates, corresponding to the stringency of the lepton identification criteria. For electrons, an updated version of the criteria from ref.
[32] is utilized, while the muon identification is as described in ref. [33]. Tight and loose electrons are furthermore required to satisfy p T > 20 GeV, while tight (loose) muons have p T > 20 (10) GeV. The relative lepton isolation, I rel , is calculated as the sum of the transverse momenta of charged-hadron, neutral-hadron, and photon PF candidates, inside a cone of ∆R = (∆η) 2 + (∆φ) 2 = 0.4 around the lepton, divided by the lepton p T . An estimated contribution from pileup is subtracted in this calculation. Tight electrons must satisfy I rel 0.06, while loose ones are required to have I rel 0.18 (0.16) in the barrel (endcap) region. Tight (loose) muon candidates must satisfy I rel < 0. 15 (0.25).
The variable p miss T , referred to as the missing transverse momentum [34], is defined as the magnitude of the missing transverse momentum vector p miss T , which is the projection on the plane perpendicular to the beams of the negative vector sum of the momenta of all reconstructed PF candidates in the event. The energy corrections applied to the jets are propagated to the p miss T calculation.

Data and simulated event samples
This analysis is performed on a pp collision data set recorded during 2016, at a centerof-mass energy of 13 TeV. The total integrated luminosity of the collected data sample is 35.9 ± 0.9 fb −1 [35].
The single-electron (single-muon) data sample is selected with triggers [36] based on the presence of an isolated electron (muon). The dielectron, electron-muon, and dimuon data samples are selected with triggers that require the presence of two leptons of the corresponding flavors. In order to increase the selection efficiency for dilepton events in which the subleading lepton has a relatively low p T , all the dilepton samples are further extended with events that pass the single-lepton but not the dilepton triggers.
In order to compare the collected data to theoretical predictions, Monte Carlo (MC) samples are produced with events simulating the Φ → tt signal and SM background processes. The signal is simulated at leading order (LO) accuracy in perturbative QCD using a custom model in the MadGraph5 amc@nlo 2.5.1 event generator [37] that implements the top quark loop of the gluon fusion production via an effective coupling between the new boson and gluons [38]. The generator employs the NNPDF3.0 parton distribution functions (PDFs) [39], and is interfaced with pythia 8.212 [40] for fragmentation and hadronization, with the CUETP8M1 underlying event tune [41,42]. Signal event samples are produced for (pseudo)scalar boson masses of 400, 500, 600, and 750 GeV, with relative total decay widths Γ Φ /m Φ of 2.5, 5, 10, and 25% for each mass scenario. The factorization and renormalization scales, µ F and µ R , are set on an event-by-event basis to m tt /2, following the choice in ref. [43]. The top quarks from the heavy Higgs boson decay are decayed in MadGraph5 amc@nlo, preserving their spin correlations. Samples are generated for -5 -events corresponding to the resonant heavy boson signal, and for events corresponding to interference terms in the matrix element calculation between the signal and SM tt background. Events in the interference samples can receive negative weights, which reflects the sign of the corresponding part of the squared matrix element in the presence of a destructive interference. Since the heavy Higgs boson is produced via gluon fusion with a top quark loop, the coupling between the boson and the top quark appears twice in the matrix element. As a result, events originating from the resonance (interference) matrix element terms correspond to a cross section that is proportional to g 4 Φtt (g 2 Φtt ).
We calculate the next-to-next-to-LO (NNLO) cross sections for the resonant part of a given signal using the SusHi 1.6.1 program [44]. The ratio of the NNLO cross section over the LO cross section computed with MadGraph5 amc@nlo determines the K factor, typically of size ≈2, applied to the resonant part of the signal. The K factors applied to the interference component of the signal are obtained as the geometric mean of the K factors of the resonant signal and the SM tt process [43]. The SM tt K factor is 1.6, calculated as the ratio between the NNLO cross section used for the simulated tt sample, as described below, and the LO cross section obtained in a similar setup. The K factors for the resonant part of the signal and the interference are applied throughout this analysis. In the hMSSM interpretation, we also use the 2hdmc program [45] to calculate, for given m A and tan β, the mass of the H boson and the widths of both heavy Higgs bosons, as well as other MSSM parameters.
The main SM background contribution in this analysis originates from tt production. Other background events originate from single top quark production, single boson production (Drell-Yan Z/γ * + jets and W + jets), diboson processes (WW, WZ, and ZZ), tt production in association with a Z or W boson (commonly referred to as ttV), and QCD multijet processes.
The tt process is simulated to next-to-LO (NLO) using the powheg v2 generator [46][47][48][49], assuming a top quark mass of 172.5 GeV. The factorization and renormalization scales are set to m 2 t + p 2 T,t , where m t and p T,t are the mass and the transverse momentum of the top quarks in the underlying Born-level configuration. The NNPDF3.0 PDF set is used, and the events are passed to pythia with the CUETP8M2T4 event tune [50]. The predicted tt production cross section is 831.8 +19. 8 −29.2 (scale) ± 35.1 (PDF + α S ) pb, as calculated with the Top++2.0 program to NNLO in perturbative QCD, including softgluon resummation to next-to-next-to-leading logarithmic (NNLL) order (as discussed in ref. [51] and references therein), and assuming a top quark mass of 172.5 GeV. The first uncertainty comes from the independent variation of the factorization and renormalization scales, while the second one is associated to variations in the PDF and strong coupling α S , following the PDF4LHC prescription with the MSTW2008 68% confidence level NNLO, CT10 NNLO, and NNPDF2.3 5f FFN PDF sets (as discussed in ref. [52] and references therein, and refs. [53][54][55]). The modeling of SM tt production in powheg is known to predict a harder p T spectrum of the top quarks than observed in the data. An empirical reweighting for top quark pairs based on the p T spectrum of generator-level top quarks is applied to obtain a better agreement with the measured differential tt cross section [56,57].

Data analysis
We search for heavy Higgs bosons decaying into a top quark pair, in final states with either one or two leptons, where the lepton is either an electron or a muon. The analysis strategy for single-lepton final states differs from the one for dilepton final states, due to differences in the event selection, SM background composition, kinematic top quark pair reconstruction, and definition of observables that discriminate between the SM background and the signal.

Single-lepton final state
In the single-lepton channel we aim to select events originating from top quark pair decays to a leptonically decaying W boson and a hadronically decaying W boson. The targeted topology is therefore tt → + ν b qq b, where denotes an electron or a muon and the leptonic and hadronic W boson decays may be swapped. Events in the single-electron (single-muon) channel in data and simulation are required to pass a single-electron (singlemuon) trigger, as explained in section 3. Selected events must have exactly one tight electron (muon) with p T > 30 (26) GeV and |η SC | < 2.5 (|η| < 2.4), where η SC is the pseudorapidity of the ECAL supercluster associated with the electron [32]. To suppress the contribution from Z/γ * + jets and other processes in which multiple prompt leptons are produced, an event is rejected if an additional loose electron or loose muon is found.

JHEP04(2020)171
This also ensures orthogonality to the event selection of the dilepton analysis outlined in section 4.2. An event must contain at least four jets with p T > 20 GeV and |η| < 2.4, at least two of which are required to be b tagged. To further suppress the QCD multijet background, only events with m W T > 50 GeV are selected. The transverse mass variable is defined as m W T = 2p T p miss T [1 − cos ∆φ( p T , p miss T )], with p T being the transverse momentum of the only tight electron or muon in the event.
Each event that passes the selection described above is reconstructed under the assumption that it has been produced in the process tt → ν bbqq . The tt system is reconstructed using an approach similar to the one adopted in ref. [56]. All possible ways to assign four reconstructed jets to the four quarks in the final state are considered. To reduce the number of combinations, each of the two b quarks is required to be associated to a b-tagged jet. For each considered choice of a jet assigned to a b quark stemming from the semileptonically decaying top quark, an attempt is made to reconstruct the transverse momentum of the neutrino, p ν T [68]. This is achieved by imposing a constraint on the mass of the semileptonically decaying top quark and on the mass of the leptonically decaying W boson, and choosing the unique solution of the neutrino momentum that minimizes the distance D ν = | p ν T − p miss T |. Next, a likelihood function is constructed from the probability density function of the minimal value of D ν , and the two-dimensional (2D) probability density function of the reconstructed mass of the top quark and the W boson in the hadronic side of the tt decay. The jet assignment with the largest value of this likelihood is chosen to reconstruct the tt system.
The performance of the tt reconstruction algorithm is studied using SM tt simulation, considering only events with the targeted decays at the generator level, tt → ν bbqq . A correct jet-quark assignment exists for 44% of such events that pass the event selection. For those events where a correct assignment exists and a solution can be found for the neutrino momentum, the probability that all four jets are correctly assigned to the quarks varies from around 60 to 80%, depending on the value of the invariant mass of the generator-level top quark pair, m gen tt . The relative m tt resolution, as computed with all selected events with targeted decays, changes from about 17 to 21%, depending on m gen tt . The tt reconstruction results in a solution in about 85% of observed events, and only the events with a solution are considered for further analysis. The search for the Φ → tt signal is performed using two observables. The first one is the invariant mass m tt , as obtained from the tt reconstruction algorithm, and probes the mass of the heavy (pseudo)scalar boson. The second observable is |cos θ * t |, where θ * t denotes the angle between the momentum of the semileptonically decaying top quark in the tt rest frame and the momentum of the tt system in the laboratory frame, as illustrated in figure 2. In a Φ → tt signal process, the heavy Higgs boson would decay into top quarks isotropically, resulting in a flat distribution of cos θ * t at the generator level. The SM tt production, on the other hand, yields a distribution peaking at ±1. As a result, the distribution of |cos θ * t | is relatively enriched in signal events towards |cos θ * t | = 0. The QCD multijet background is estimated from dedicated control regions in the data, independently in the electron and muon channels. The total event yield for this background is obtained from data using a variant of the ABCD method [69]. Four regions are defined, based on the relative isolation of the lepton (smaller or greater than the isolation threshold imposed on the tight lepton), and the m W T variable (smaller or greater than 50 GeV). The three regions complementary to the signal region are relatively enriched in the multijet background. The overall rate of QCD multijet background in the signal region is estimated with a simultaneous fit to the numbers of events observed in the four regions, exploiting the factorization of the distribution of this background in (m W T , I rel ). The shape of the multijet distribution of the observables m tt and |cos θ * t | is determined from data using events with an inverted lepton isolation selection applied, after subtracting the contributions from other backgrounds.
The data and the expected SM background yields are shown in table 1, for selected events that have a solution of the tt reconstruction algorithm. The background predictions are computed with the help of a maximum likelihood fit to the data using the backgroundonly version of the full statistical model that will be described in section 6. The uncertainties obtained in this fit, referred to as post-fit uncertainties, are reported.
The observed and post-fit predicted distributions of m tt in different |cos θ * t | regions are shown in figure 3. The impact of the signal process (including the interference) for the best-fit signal hypothesis, which will be discussed in section 6.2, is shown in the lower panels. It demonstrates the characteristic peak-dip lineshape discussed in section 1. For this benchmark and also in general, the contributions from both the resonant part and the interference are important. The relative importance of the latter increases with the total width or as the coupling modifier decreases.

Dilepton final state
In the dilepton channel we aim to select signal events where both top quarks decay to a leptonically decaying W boson. Hence, the targeted decay topology is tt → + ν b − ν b. Events in the dielectron (ee), electron-muon (eµ), and dimuon (µµ) channel in data and simulation are required to pass a dielectron, electron-muon, and dimuon trigger, respectively, or a single-lepton trigger, as explained in section 3. The subsequent event selection closely follows ref. [70]. Events are required to contain exactly one pair of oppositely charged -9 -JHEP04(2020)171 tight leptons, in which the leading (subleading) lepton has p T > 25 (20) GeV. Events are rejected if they contain additional tight electrons or additional muons that satisfy the tight identification criteria with the exception of a looser selection on isolation, I rel < 0.25. The selected dilepton pair is further required to have an invariant mass of at least 20 GeV, to suppress events from low-mass dilepton resonances. In the ee and µµ channels, events are rejected if they contain a dilepton pair consistent with a decay of a Z boson, namely, if they have an invariant mass in the range 76-106 GeV. Each event must have at least two jets with p T > 30 GeV and |η| < 2.4. Additional jets with p T > 20 GeV and |η| < 2.4 in the event are also considered for further analysis. At least one of the jets is required to be b tagged. In the ee and µµ channels, the p miss T must exceed 40 GeV to further suppress Z/γ * background events.
The contribution from the Z/γ * + jets background process to the selected event yield is estimated using control regions in data, following the procedure described in ref. [71]. The yield outside the Z boson mass window is estimated based on the observed event yield within the window, using the knowledge from simulation of the ratio of Z/γ * + jets inside and outside of the mass window. In this estimation, the non-Z/γ * + jets background contribution within the Z boson mass window is taken from the eµ channel, where the Z/γ * + jets contribution is negligible, and corrected for lepton reconstruction effects before being subtracted from the observed number of events in the Z boson mass window. Using this method, the theoretically predicted Z/γ * + jets event yield is scaled by a factor 1.22, 1.20, and 1.19 in the ee, eµ, and µµ channels respectively.
Each event that passes the selection described above is reconstructed under the assumption that it has been produced in a process tt → + ν b − ν b. A kinematic reconstruction algorithm [72] is applied to reconstruct the tt system. All jets with p T > 20 GeV are considered in the reconstruction of the tt system. Given an assignment of jets and   Observed and expected distributions of m tt in different |cos θ * t | regions in the e + jets (upper) and µ + jets (lower) channels. The expected distributions have been obtained with a background-only fit to the data, and an approximate post-fit uncertainty is shown with a gray band. The impact of the best-fit signal is included in the lower panels for illustration.
leptons to the underlying expected tt decay products, a system of equations is constructed that imposes constraints on the reconstructed W boson mass and the reconstructed top quark and antiquark masses. The transverse momentum imbalance, represented by p miss T , is assumed to originate solely from the two neutrinos. Detector resolution effects are taken into account by sampling both the measured energy and the direction of the leptons and b jet candidates within their respective experimental resolutions. For each sampling, the solution of the equations that results in the smallest m tt is chosen. Per event, 100 sam- (2020) plings are performed, and each is assigned a weight based on the probability density of the invariant mass of the lepton and b jet from the top quark decay. The kinematic properties of the top quark and antiquark are obtained as a weighted average over all samplings. Finally, the assignment of jets resulting in the maximum sum of weights over all samplings is chosen, and preference is given to a jet assignment that contains two b-tagged jets.
The performance of the tt reconstruction algorithm is studied using simulated SM tt events with targeted decays at the generator level, tt → + ν b − ν b. In 75% of the selected events the two generator-level jets are within the acceptance. For those events for which the algorithm has a solution, the probability to correctly match both b jets as chosen by the algorithm to jets originating from b quarks from the top quark decays is 55 to 85%, depending on the value of m gen tt . The m tt resolution, computed using all selected SM tt events with targeted decays, ranges from 20 to 28%, depending on m gen tt . The events for which the tt reconstruction results in a solution, which is the case in about 95% of observed events, are considered for further analysis. The resulting event yields for data and SM background expectations are shown in table 2.
The search for the Φ → tt signal is performed using two observables. The first one is the invariant mass m tt , obtained from the tt reconstruction algorithm. The second observable is a spin correlation variable constructed from the charged leptons in the event. Charged leptons have the highest spin analyzing power amongst the top quark decay products [73], and their properties can be measured precisely. The chosen variable is the cosine of the angle between the charged lepton momenta in their respective helicity frames, and is denoted by c hel . The four-momenta of the leptons in their helicity frames are obtained by first boosting the leptons into the tt rest frame and then boosting them along their parent top quark or antiquark directions in this frame. The distribution of c hel is sensitive to the spin and CP state of the tt system. At the generator level and with no requirements on acceptance, the distribution is linear in shape; its slope is maximally positive for the A resonance and is mildly negative for the H resonance. On the other hand, the slope for the SM tt production, integrated over the entire m tt range, is mildly positive. This allows discriminating both between the signal and background processes and between the H and

Systematic uncertainties
Various sources of uncertainty affect the distributions of the observables used to search for a heavy Higgs boson signal. Below we describe the experimental and theoretical systematic effects considered in the analysis. In the statistical evaluation discussed in section 6, each source of uncertainty corresponds to a nuisance parameter in a binned maximumlikelihood fit to the distributions of the observables in data. Uncertainties that affect only the normalization are modeled using log-normal constraints, while Gaussian constraints are imposed for nuisance parameters that control all other uncertainties. Unless stated otherwise, all uncertainties are evaluated on signal as well as background processes and treated as fully correlated among the processes and lepton channels. The uncertainties are summarized in table 3.
The uncertainty due to the jet p T scale [29] is evaluated by varying the corresponding corrections within their uncertainties. The events are reanalyzed, by reapplying the event selection and recalculating all kinematic quantities. A total of 19 independent jet momentum correction uncertainties affecting jets in the tracker acceptance are considered. As the jet p T resolution in simulation is smeared to match the resolution observed in data, a corresponding uncertainty is evaluated. An uncertainty in the unclustered component of p miss T is computed by shifting the energies of PF candidates not clustered into jets with p T > 15 GeV according to the energy resolution for each type of PF candidate [34]. Uncer--13 -JHEP04(2020)171 tainties in the b tagging efficiency scale factors applied to simulated events are evaluated by varying them within the respective uncertainties [31]. The scale factors for heavy-flavor (b and c) jets are varied independently of those for light-flavor jets. The uncertainties in the trigger scale factors as well as the electron and muon identification scale factors are considered [32,33], where the lepton identification also includes effects originating from the isolation requirement and the track reconstruction. The uncertainties in the trigger efficiency scale factors for the single-electron and single-muon channels are considered not correlated with each other, but each of them is independently correlated with the uncertainty in the trigger scale factors in the dilepton channel, with a 50% correlation coefficient. Effects due to the uncertainty in the distribution of the number of pileup interactions are evaluated by varying the effective inelastic proton-proton cross section in the simulation by 4.6% from its nominal value. The uncertainty in the integrated luminosity amounts to 2.5% [35] and affects the normalization of all simulated processes.
The prediction of the SM tt production, the main background process in the analysis, is affected by various sources of theoretical uncertainties. The overall normalization of the SM tt background is assigned an uncertainty of 6%, from the NNLO + NNLL QCD cross section calculation [51,52] that is used to normalize the events. The effect of the choice of the renormalization and factorization scales, µ R and µ F , in the matrix element is evaluated by varying these scales independently by a factor of 2 and 1/2. The effect on the acceptance is considered, but the effect on the cross section is ignored as it is already included in the considered uncertainty in the cross section. The renormalization scales used in the parton shower simulation of initial-state radiation (ISR) and final-state radiation (FSR) are also varied independently by a factor of 2 in each direction. The effect of the uncertainty in the amount of ISR, as well as in the underlying event tune used in the simulation, was found not to be statistically significant. The uncertainty in the top quark mass is considered by shifting m t in the simulation by ±3 GeV and rescaling the induced variations by a factor of 1/6 to emulate a more realistic top quark mass uncertainty of 0.5 GeV [74]. The uncertainty in the matching scale between the matrix element and the parton shower is evaluated by varying the powheg parameter, h damp , that controls the suppression of radiation of additional high-p T jets [50]. The nominal value of h damp in the simulation is 1.58 m t , and the varied values are 0.99 m t and 2.24 m t . The uncertainty arising from the choice of the PDF set is evaluated by reweighting the simulated tt events using 100 replicas of the NNPDF3.0 set. A principal component analysis is performed on the variations from the PDF replicas to construct two base variations, such that the deviation from the nominal distribution given by each replica can be described as a linear combination of the base variations. The uncertainty in the α S parameter used in the PDF set induces a third independent PDF variation. The uncertainty accounting for the mismodeling of the p T spectrum of top quarks is evaluated by varying the two parameters used in the top quark p T reweighting function.
The renormalization and factorization scale uncertainties in the heavy Higgs boson signal simulation are treated independently for the resonant and interference components. Compared to the alternative of varying the scales for the two components simultaneously, we found this to be the more conservative option. The effect on the acceptance as well as on the cross section is considered. Since the simulated samples have been generated at -14 -

JHEP04(2020)171
LO accuracy, the total effect on the cross section reaches values in excess of 30%. Other theoretical uncertainties in the signal, such as the uncertainties in m t or PDF, are neglected as they are expected to be small compared to this variation.
The expected yields for most of the non-tt background processes are derived using theoretical predictions for the cross sections at NLO or higher accuracy. The uncertainties assumed in the normalization of these processes are conservative and always exceed those of the corresponding theoretical computations. For the single top quark production in the t (tW) channel we assign an uncertainty of 20 (15)%, which is based on the measurements of the cross sections of these processes [75][76][77]. The uncertainty in the ttV production is taken to be 30%, which covers the uncertainties of the experimental measurements [78,79]. To account for the fact that this search probes a restricted region of the phase space of the corresponding processes, we assign uncertainties of 50% for W + jets and Z/γ * + jets production (only in the single-lepton channel) and 30% for the diboson production. Finally, a 20% uncertainty is used for the s-channel single top quark production, for which no measurement at the LHC exists. The adopted conservative normalization uncertainties have little impact on the sensitivity of this search due to the small contribution of these processes.
In cases where the normalization of a background process is estimated using a datadriven method, the corresponding uncertainty is determined by the same method. For the Drell-Yan background in the dilepton channel we assign a 30% uncertainty, from the variation in the scale factors when derived with and without the requirements on p miss T , the jet b tag decisions, and the tt reconstruction. The normalization of the QCD multijet background, which is only relevant in the single-lepton channel, is assigned an uncertainty of +100/−50%, independently in the single-electron and single-muon channels. It covers the statistical uncertainty in the underlying fit and the difference between the data-driven and MC-based estimations.
The nominal background prediction is affected by the limited size of the simulated MC event samples. This statistical uncertainty is evaluated using the "light" Barlow-Beeston method [80], by introducing one additional nuisance parameter per bin of the 2D distribution of the observables.
Several systematic variations in the background, most notably those constructed from dedicated MC samples, are affected by statistical fluctuations. We suppress these fluctuations by smoothing the relative deviations from the nominal background distribution of m tt and the angular variable. The up and down deviations for each independent uncertainty are assumed to be symmetric in shape, but allowed to differ in the overall size. The symmetrized deviation is smoothed using a version of the LOWESS algorithm (LOcally WEighted Scatterplot Smoothing) [81,82]. In the vicinity of each bin of the 2D distribution, the symmetrized deviation is fitted with a plane using a weighted least squares fit, in which nearby bins receive larger weights. The smoothed deviation obtained in this way is rescaled to account for the overall size of the input up or down deviation, and applied to the nominal background expectation in the given bin. A similar procedure is also applied to all signal distributions. The resulting distributions are then used in the subsequent analysis to evaluate the systematic uncertainty under consideration.   Table 3. The systematic uncertainties considered in the analysis, indicating the number of corresponding nuisance parameters (when more than one) in the statistical model, the type (affecting shape or only normalization), the affected processes, and the correlation among the lepton channels. Uncertainties tagged in the last column with "All" are correlated among the single-lepton and dilepton channels. In case an uncertainty is only applicable to the single-electron, the single-muon, the single-lepton, or the dilepton channel, they are indicated with e, µ, , , respectively.

JHEP04(2020)171
In general, the relative importance of different systematic uncertainties depends greatly on the signal hypothesis, especially the mass of the heavy Higgs boson. Typically, among the uncertainties with the largest impact are the signal theoretical uncertainties and some of the jet momentum correction uncertainties. Close to the tt production threshold, the variations in m t and the h damp parameter become important, while for larger m Φ the PDF, µ R , and µ F variations in the SM tt background can have significant impacts. For certain signal hypotheses, the signal distribution partially resembles one of a minor background; in such cases the variation of the normalization of the respective background becomes relevant. In addition, MC statistical uncertainties, when grouped together, often outweigh every other individual uncertainty.

Results
To evaluate the consistency of the observed data with the presence of a signal, we perform a statistical analysis using the 2D binned distribution of (m tt , |cos θ * t |) in the single-electron and the single-muon channels separately, and the 2D binned distribution of (m tt , c hel ) in the combined dilepton channel. The single-lepton and dilepton channels do not overlap as they correspond to orthogonal lepton selection criteria.
The statistical model is defined by the likelihood function with b i denoting the combined background yield in a given bin i, s Φ R,i and s Φ I,i the signal yields in a given bin for the resonant and interference part, respectively, ν the vector of nuisance parameters (on which the signal and background yields generally depend), n i the observed yield, and g Φtt the coupling strength modifiers given by eq. (1.1). The parameters of the signal model (mass m Φ , width Γ Φ , and g Φtt ) are collectively denoted by vector p. Eq. (6.1) is kept generic by including contributions from both CP states. As there is no interference between them, the corresponding signal distributions are trivially added together. We also introduce an auxiliary overall signal strength modifier µ, which rescales the full beyond the SM (BSM) contribution. This allows testing different signal hypotheses in a computationally efficient way, as will be detailed below. The external constraints on the nuisance parameters are taken into account in this likelihood via a product of corresponding probability density functions, G(ν).
The background-only model is constructed by setting µ = 0 in eq. (6.1). To quantify the level of agreement between it and observed data, we perform a goodness-of-fit test based on the so-called "saturated model" [83]. This yields a p-value of 0.43, indicating a good overall agreement.
We perform scans over the parameters of the signal models, p. A variant of the LHC profile likelihood ratio test statisticq µ from refs. [84,85] is utilized:

JHEP04(2020)171
The test statistic is expressed in terms of the auxiliary parameter µ in eq. (6.1), in which the statistical model is linear, while the parameters p are kept fixed at their values being probed in the scan. The likelihood in the numerator is maximized with respect to the nuisance parameters, andν µ,p denotes the vector of their values at the maximum. A similar notation is used in the denominator, where the likelihood is maximized with respect to both µ and ν, under the additional constraint 0 μ p µ. The requirementμ p 0 excludes from the consideration cases in which the shape of the overall BSM contribution gets flipped, resulting in a qualitatively different effect from what is targeted in this search. The conditionμ p µ prevents the exclusion of a signal hypothesis if the data are more compatible with a model that predicts the BSM contribution of a similar shape but a larger overall size.
For each signal hypothesis p, we perform a test according to the CL s criterion [86,87]. This is done for µ = 1 in eq. (6.2), which reproduces the nominal signal expectation. We profit from the known asymptotic approximation [84] for distributions of the adopted test statistic to construct these distributions in a computationally efficient way. If the CL s value computed for µ = 1 and given p is found to be smaller than 0.05, the point p is said to be excluded at 95% confidence level (CL).

Interpolation and extrapolation of signal masses and widths
To construct expected signal distributions for every point encountered in the scans, we apply an interpolation in mass and width of the heavy Higgs boson, starting from the reference generated points. This is done independently for the m tt distribution in each bin of the angular variable. We consider the resonant part of the signal and the interference separately and further split the interference contribution in two according to the sign of the per-event generator weight. A change in the mass results in a horizontal shift of the m tt distribution. The interpolation in this observable is implemented with a nonlinear morphing algorithm [88]. On the other hand, the effect of a change in Γ Φ is evaluated with an independent interpolation in each bin.
The signal model parameter scan may reach values of Γ Φ /m Φ below 2.5%, the lowest value considered in the simulated signal samples. Since the reconstructed m tt resolution is about 17% or worse, the shape of the m tt distribution for a signal with such low widths does not differ from the one corresponding to Γ Φ /m Φ = 2.5%. Hence, for scan points with Γ Φ /m Φ < 2.5% it is sufficient to use the distributions for Γ Φ /m Φ = 2.5% and only scale the cross sections appropriately.

Model-independent interpretation
Constraints on the coupling strength modifier g Φtt are derived as a function of the mass and width of the heavy Higgs boson, for each CP state independently. The coupling modifier for the other CP state in eq. (6.1) is set to zero to exclude it from the statistical model. The scan is performed for m Φ between 400 and 750 GeV and Γ Φ /m Φ between 0.5 and 25%. The mass and width interpolation described in section 6.1 is performed in scan points other than those corresponding to the generated signal samples. Coupling strength values up to 3 are probed to guarantee that the amplitudes preserve perturbative unitarity for all -18 -JHEP04(2020)171 calculations, in accordance with the lower bound tan β = 1/g Att 0.3 given in ref. [4] in the context of 2HDMs.
The constraints obtained on g Φtt are presented in figures 5 and 6 for the scalar and the pseudoscalar scenarios, respectively. Since the total width Γ Φ is kept fixed during the scans and the partial width of Φ → tt is proportional to g 2 Φtt , in some regions the partial width can exceed Γ Φ . These unphysical regions are marked in the figures with hatched lines. In some cases the observed exclusion for a given mass does not extend continuously all the way to the largest probed g Φtt = 3 (e.g., m H ≈ 700 GeV in the panel for Γ H /m H = 25% in figure 5). This is due to the strong dependence of the shape of the signal distribution on the value of the coupling strength modifier. For some values of g Φtt the shape becomes compatible with systematic variations in the background.
As evident from figure 6, there is a signal-like excess for the pseudoscalar hypotheses with low masses. The largest deviation from the SM background is observed for a pseudoscalar Higgs boson with a mass of 400 GeV and a total relative width of 4%, with a local significance of 3.5 ± 0.3 standard deviations. Figure 7 shows scans of −2 ln[L(g Att )/L SM ] for this hypothesis, as a function of the coupling modifier g Att . The likelihoods L(g Att ) and L SM are given by eq. (6.1). They are computed for µ = 1 and 0 respectively and in both cases maximized with respect to all nuisance parameters. The scans are shown for the observed data, as well as for the expectations under the background-only hypothesis and in the presence of the signal. In the latter case the coupling modifier is set to the value obtained in the combined fit, g Att ≈ 0.9. In this case and in general, the expected sensitivity is comparable between the single-lepton and dilepton channels. The single-lepton channel is slightly more sensitive for the scalar hypotheses and in the case of the pseudoscalar hypotheses with larger masses. However, figure 7 demonstrates that the observed excess is driven by the dilepton channel, which is also supported by the comparisons between the observed and expected distributions in figures 3 and 4.
When accounting for the look-elsewhere effect [89] in the mass, total width, and CP state of the heavy Higgs boson, the significance of the excess is 1.9 standard deviations, which corresponds to a p-value of 0.028. We note that higher-order electroweak corrections to the SM tt production can become important in the vicinity of the pair production threshold [90] and may account for the excess.

Interpretation in the hMSSM
For the hMSSM, we perform a scan over the two model parameters, m A and tan β. Both A and H bosons are included. For each point, the coupling strength modifiers g Att and g Htt , the mass m H , and the widths of the two heavy Higgs bosons are determined with the 2HDMC program. The CP-even state is more massive of the two, but the mass separation ∆m = m H − m A decreases with m A and tan β. Typical values of ∆m/m A vary from ≈20% for m A = 400 GeV, tan β = 0.5 to ≈1% for m A = 700 GeV, tan β = 2. The scan is performed for m A between 400 and 700 GeV in steps of 12.5 GeV, and tan β between 0.4 and 5.0 in steps of 0.2. Similarly to what was done above, the lower boundary tan β > 0.4 is imposed to assure perturbative unitarity [4]. The mass and width interpolation described  Figure 5. Model-independent constraints on the coupling strength modifier as a function of the heavy scalar boson mass, for relative widths of 0. 5, 1, 2.5, 5, 10, and 25%. The observed constraints are indicated by the blue shaded area. The inner (green) band and the outer (yellow) band indicate the regions containing 68 and 95%, respectively, of the distribution of constraints expected under the background-only hypothesis. The unphysical region of phase space in which the partial width Γ H→tt becomes larger than the total width is indicated by the hatched lines.
-20 -  Figure 6. Model-independent constraints on the coupling strength modifier as a function of the heavy pseudoscalar boson mass, for relative widths of 0. 5, 1, 2.5, 5, 10, and 25%. The observed constraints are indicated by the blue shaded area. The inner (green) band and the outer (yellow) band indicate the regions containing 68 and 95%, respectively, of the distribution of constraints expected under the background-only hypothesis. The unphysical region of phase space in which the partial width Γ A→tt becomes larger than the total width is indicated by the hatched lines. in section 6.1 is performed in scan points other than those corresponding to the generated signal samples. The expected and observed exclusions in the (m A , tan β) plane are presented in figure 8. The upper boundary of the observed (expected) exclusion in tan β varies from 1.0 (2.3) at m A = 400 GeV to 1.5 (0.8) at m A = 700 GeV. The tension between the observed exclusion and the expectation at low m A is a manifestation of the excess discussed above. These results can be compared to those of the search for H ± → tb/tb in ref. [23], which were also interpreted in the hMSSM benchmark, setting constraints in the (m H ± , tan β) plane. Translating the results from ref. [23] in terms of m A , the present analysis observes a more stringent exclusion in tan β for m A ≈ 700 GeV, while the exclusion for m A ≈ 400 GeV is substantially weaker than in the reference due to the observed signal-like deviation. The expected exclusion is tighter than in ref.
[23] throughout the considered m A range.

Summary
Results are presented for the search for additional heavy Higgs bosons decaying to a pair of top quarks. A data sample recorded with the CMS detector at √ s = 13 TeV, corresponding to an integrated luminosity of 35.9 fb −1 , is analyzed. The final states with one or two leptons are utilized. The invariant mass of the reconstructed tt system as well as angular variables sensitive to the spin of the new boson are used to search for the signal, while taking into account the interference with the standard model tt production. A moderate signal-like deviation is observed for the hypothesis of a pseudoscalar Higgs boson with the mass m A ≈ 400 GeV. After accounting for the look-elsewhere effect, its significance is 1.9 standard deviations. Further improvements of the theoretical description of the standard model tt process in the vicinity of the production threshold will be needed to clarify the origin of this deviation.
Constraints on the strength of the coupling of the sought-for boson to top quarks are reported, separately for the scalar and pseudoscalar cases, for the mass ranging from 400 to 750 GeV and the total relative width from 0.5 to 25%. These are the most stringent constraints on this coupling to date. The results are also interpreted in the hMSSM scenario in the minimal supersymmetric standard model. This search probes the values of m A from 400 to 700 GeV and excludes, at 95% confidence level, the region with values of tan β below 1.0 to 1.5, depending on m A . This extends the exclusion obtained in previous searches.

JHEP04(2020)171
Acknowledgments We congratulate our colleagues in the CERN accelerator departments for the excellent performance of the LHC and thank the technical and administrative staffs at CERN and at other CMS institutes for their contributions to the success of the CMS effort. In addition, we gratefully acknowledge the computing centers and personnel of the Worldwide LHC Computing Grid for delivering so effectively 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: BMBWF and FWF (  Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.