Search for a heavy Higgs boson decaying into a Z boson and another heavy Higgs boson in the (cid:2)(cid:2) bb ﬁnal state in pp collisions at √ s = 13 TeV with the ATLAS detector

A search for a heavy neutral Higgs boson, A , decaying into a Z boson and another heavy Higgs boson, H , is performed using a data sample corresponding to an integrated luminosity of 36.1 fb − 1 from proton– proton collisions at √ s = 13 TeV recorded in 2015 and 2016 by the ATLAS detector at the Large Hadron Collider. The search considers the Z boson decaying to electrons or muons and the H boson into a pair of b -quarks. No evidence for the production of an A boson is found. Considering each production process separately, the 95% conﬁdence-level upper limits on the pp → A → ZH production cross-section times the branching ratio H → bb are in the range of 14–830 fb for the gluon–gluon fusion process and 26–570 fb for the b -associated process for the mass ranges 130–700 GeV of the H boson and 230–800 GeV of the A boson. The results are interpreted in the context of two-Higgs-doublet models.

previous search by considering explicitly the gluon-gluon fusion and b-associated production processes as well as both narrow and wide widths of the A boson.

ATLAS detector
The ATLAS detector is a general-purpose particle detector, described in detail in Ref. [27]. It includes an inner detector surrounded by a 2 T superconducting solenoid, electromagnetic and hadronic calorimeters and a muon spectrometer with a toroidal magnetic field. The inner detector consists of a high-granularity silicon pixel detector, including the insertable B-layer [28] installed in 2014, a silicon microstrip detector, and a straw-tube tracker. It provides precision tracking of charged particles with pseudorapidity |η| < 2.5. 1 The calorimeter system covers the pseudorapidity range |η| < 4.9. It is composed of sampling calorimeters with either liquid argon or scintillator tiles as the active medium. The muon spectrometer provides muon identification and measurement for |η| < 2.7. A two-level trigger system [29] is employed to select events for offline analysis, which reduced the average recorded collision rate to about 1 kHz.

Data and simulation
The data used in this search were collected during 2015 and 2016 from √ s = 13 TeV proton-proton collisions and correspond to an integrated luminosity of 36.1 fb −1 , which includes only datataking periods where all relevant detector subsystems were operational. The data sample was collected using a set of single-muon and single-electron triggers. The lowest-p T trigger thresholds depend on the data-taking period and are in the range of 20-26 GeV for the single-muon triggers and 24-26 GeV for the single-electron triggers. Simulated signal events with A bosons produced by gluongluon fusion were generated at leading order with MadGraph5_aMC@NLO 2.3.3 [30,31] using Pythia 8.210 [32] with a set of tuned parameters called the A14 tune [33] for parton showering. For the generation of A bosons produced in association with b-quarks, MadGraph5_aMC@NLO 2.1.2 [31,34,35] was used following Ref. [36] together with Pythia 8.212 and the A14 tune for parton showering. The gluon-gluon fusion production used NNPDF2.3LO [37] as the parton distribution functions (PDF), while the b-associated production used CT10nlo_nf4 [38]. The signal samples were generated for A bosons with masses in the range of 230-800 GeV and widths up to 20% of the mass and for narrowwidth H bosons with masses in the range of 130-700 GeV.
Background events from the production of W and Z bosons in association with jets were simulated with Sherpa 2.2.1 [39] using the NNPDF3.0NNLO PDF set [40]. Top-quark-pair production was simulated with Powheg-Box v2 [41][42][43] and the CT10nlo PDF set [38], while the electroweak single-top-quark production was simulated with Powheg-Box v1 and the fixed four-flavour PDF set CT10nlo_f4 [38]. The parton shower was performed with Pythia 6.428 [44] using the Perugia 2012 set of tuned parameters [45]. The production of top-quark pairs in association with a vector boson was simulated using MadGraph5_aMC@NLO 2.2.3 and the NNPDF3.0NLO PDF set, whereas Pythia 8.186 was used 1 ATLAS uses a right-handed coordinate system with its origin at the nominal interaction point (IP) in the centre of the detector and the z-axis along the beam pipe. The x-axis points from the IP to the centre of the LHC ring, and the y-axis points upward. Cylindrical coordinates (r, φ) are used in the transverse plane, φ being the azimuthal angle around the beam pipe. The pseudorapidity is defined in terms of the polar angle, θ , as η = − ln tan(θ/2). Transverse momenta are computed from the three-momenta, p, as p T = | p| sin θ .
for the parton shower with the A14 tune. Production of W W , Z Z and W Z pairs was simulated using Sherpa 2.2.1 and the NNPDF3.0NNLO PDF set. Finally, SM Higgs boson production in association with a Z boson was generated with Powheg-Box v2 and the NNPDF3.0NLO PDF set, whereas the parton shower was performed with Pythia 8.186 using the AZNLO tune [46].
The modelling of bottom-and charm-hadron decays was performed with the EvtGen v1.2.0 package [47] for all samples apart from those simulated with Sherpa. The simulated events were overlaid with inelastic proton-proton collisions to account for the effect of multiple interactions occurring in the same and neighbouring bunch crossings ('pile-up'). These events were generated using Pythia 8 with the A2 tune [48] and the MSTW2008LO PDF set [49]. The events were reweighted so that the distribution of the average number of interactions per bunch crossing agreed with the data.
All generated background samples were passed through the Geant4-based [50] detector simulation [51] of the ATLAS detector. The ATLFAST2 simulation [51] was used for the signal samples to allow for the generation of many different A and H boson masses. The simulated events were reconstructed in the same way as the data. Table 1 Summary of the event selection for signal and control regions.
Single-electron or single-muon trigger Exactly 2 leptons (e or μ) (p T > 7 GeV) with the leading one having p T > 27 GeV Opposite electric charge for μμ or eμ pairs; 80 GeV < m , m eμ < 100 GeV, = e, μ At least 2 b-jets (p T > 20 GeV) with one of them having p T > 45 GeV

Event selection
The decay A → Z H → bb features a pair of oppositely charged, same flavour leptons and two b-jets. Three resonances can be formed by combining the selected objects: the Z boson ( ), the H boson (bb) and the A boson ( bb). Moreover, additional b-jets may be present if the A boson is produced via the b-associated production mechanism. These features are used to define the event selection as summarised in Table 1.
The events recorded by the single-muon and the single-electron triggers are required to contain exactly two muons or two electrons, respectively. At least one of the leptons must have p T > 27 GeV. Only events that contain a primary vertex with at least two associated tracks with p T > 400 MeV [64] are considered. In the case of muons, they are required to have opposite electric charges. No such requirement is applied to electrons due to their non-negligible charge misidentification rates resulting from conversions of bremsstrahlung photons. The invariant mass of the lepton pair, m , must be in the range of 80-100 GeV to be compatible with the mass of the Z boson.
The H → bb decay is reconstructed by requiring at least two b-jets with the highest-p T one having p T > 45 GeV. When more than two b-jets are present, the two highest-p T b-jets are considered to be from the H decay. Requiring b-jets increases the fraction of top-quark background in the signal region, including top-quark pair and single-top-quark production. This is reduced by requiring E miss T / √ H T < 3.5 GeV 1/2 , where H T is the scalar sum of the p T of all jets and leptons in the event. In addition, a requirement that reduces the Z +jets background is also applied: where m bb is the four-body invariant mass of the two-lepton, two-b-jet system assigned to the A boson and the summation is performed over the p T of these objects.
Subsequently, two categories are defined: the n b = 2 category, which contains events with exactly two b-jets, and the n b ≥ 3 category, which contains events with three or more b-jets. For the gluon-gluon fusion production, 94%-97% of the events passing the above selection fall into the n b = 2 category, depending on the assumed m A and m H . However for the b-associated production, 27%-36% fall into the n b ≥ 3 category. The remaining b-associated produced signal events are categorised as n b = 2 events, even though more than two b-jets are expected, due to the relatively 2 The primary vertex is taken to be the reconstructed vertex with the highest p 2 T of the associated tracks.
soft p T spectrum of the associated b-jets and the geometric acceptance of the tracker. Finally, the invariant mass of the two leading b-jets, m bb , must be compatible with the assumed H boson mass by satisfying the requirement of 0.85 · m H − 20 GeV < m bb < m H + 20 GeV for the n b = 2 category, and 0.85 · m H − 25 GeV < m bb < m H + 50 GeV for the n b ≥ 3 category. The wider window for n b ≥ 3 is motivated by a slightly degraded resolution due to potential b-jet mis-assignments (see later). The overall signal efficiency of the n b = 2 category after this requirement is 5%-11% (3%-7%) for gluon-gluon fusion (b-associated production), depending on the m A and m H values.
Similarly, the efficiency of the n b ≥ 3 category is 2%-4% for the b-associated production. The signal region selection is summarised in Table 1.
The m bb distribution after the m bb requirement is used to discriminate between signal and background. To improve the m bb resolution, the bb system's four-momentum components are scaled to match the assumed H boson mass and the system's fourmomentum components are scaled to match the Z boson mass. This procedure, performed after the event selection, improves the m bb resolution by a factor of two without significantly distorting the background distributions, resulting in an A boson mass resolution of 0.3%-4%.
The dominant backgrounds after these selections are from Z +jets and top-quark production. For top-quark-pair production, a very pure (> 99% of predicted events) control region is used to determine the normalisation of the background, whereas its shape in the signal region is taken from the simulation. This control region is defined by keeping the same selection as discussed previously, apart from an opposite-flavour lepton criterion, i.e., an oppositecharge eμ pair is required instead of an ee or μμ pair (see also Table 1). The shape of the Z +jets background distribution is obtained from simulation and the normalisation is extracted from data together with the signal (see also Section 7). This procedure is possible because of the very different shapes of the m bb distributions from signal and Z +jets events. The normalisation of the Z +jets production is further constrained by a control region defined by inverting the m bb window criterion for each H boson mass hypothesis (see also Table 1). The control regions are distinct for the n b = 2 and the n b ≥ 3 categories, since the accuracy of the background simulation depends on the number of b-jets present in the event. Backgrounds from diboson, single top, and Higgs boson production, as well as top-quark-pair production in association with a vector boson, give a typical contribution of ∼ 5% to the total background. Their shapes are taken from simulation, whereas they are normalised using precise inclusive cross-sections calculated from theory. The diboson samples are normalised using next-tonext-to-leading-order (NNLO) cross-sections [65-68]. Single-topquark production and top-quark-pair production in association with vector bosons are normalised to next-to-leading-order (NLO) cross-sections from Refs. [69-71] and Ref. [31], respectively. The normalisation of the Higgs boson production in association with a vector boson follows the recommendations of Ref. [36] using NNLO QCD and NLO electroweak corrections.

Signal modelling
The good m bb mass resolution together with the fact that theory models often predict A bosons with large widths inflates the number of signal mass and width hypotheses that need to be considered. For this reason, the m bb distributions are taken from simulation of a limited number of (m A , m H ) mass points and an interpolation using analytical functions is employed for the rest.
The m bb distributions for A bosons produced by gluon-gluon fusion and with negligible widths compared with the experimental resolution are found to be adequately described by the ExpGauss- On the other hand, m bb distributions for A bosons from b-associated production, also with negligible widths compared with the experimental resolution, are better described by a double-Gaussian Crystal Ball (DSCB) function [73]: Both functions consist of a Gaussian core with mean a and variance σ 2 , whereas the rest of the parameters (k L , k H , n 1 , n 2 ) describe the tails. The DSCB function describes better than the EGEE function the slightly longer tails of the mass distribution of the b-associated production compared to gluon-gluon fusion. This is due to the few cases in which the b-quark produced in association with the Higgs boson is taken to be one of the b-quarks from the Higgs boson decay. The values of the function parameters are extracted from unbinned maximum-likelihood fits to the simulated m bb distributions. The core mean, a, is parameterised using a linear function of The cores of the m bb distributions are well-parameterised by the chosen functional forms. The small differences seen in the tails of some distributions between the functional forms and the simulations have only negligible effects on the final results, and moreover they are included as a source of systematic uncertainty. The previously described parameterisation applies to signal samples generated with narrow-width A bosons. In some regions of 2HDM parameter space relevant to this analysis, the A boson's width is significant compared with the detector resolution while the H boson's width remains negligible. In order to model the m bb shape of A bosons with large natural widths, a modified Breit-Wigner distribution 3 is convolved with the EGE and DSCB functions. The procedure is validated by comparing the results of the convolution with those of the simulated samples of A bosons with large natural widths. Widths of up to 20% of the A boson mass are considered. An example of signal distributions with large natural widths is shown in Fig. 2 for the same signal points used in Fig. 1.

Fit model and systematic uncertainties
The m bb distribution is expected to exhibit a resonant structure if signal events are present, while background events result in a smooth shape. Therefore m bb is chosen as the final signal and background discriminating variable. The shape differences in the m bb distribution between the signal and background contributions are exploited through binned maximum-likelihood fits of the signal-plus-background hypotheses to extract potential signal contributions. The fits are based on the statistical framework de- Here N i is the number of observed events and S i (m A , m H , θ) and B i ( α, θ) are the expected number of signal and estimated background events in bin i. The vector α represents free background normalisation scale factors (described later) and vector θ denotes all non-explicitly listed parameters of the likelihood function such as nuisance parameters associated with systematic uncertainties.
The function G( θ) represents constraints on θ . The parameter of interest, μ, is a multiplicative factor to the expected signal rate and is called the signal-strength parameter. The m bb bin widths are chosen according to the expected detector resolution and taking into account the statistical uncertainty related to the number of background Monte Carlo events. The bin centres are adjusted such that at least 68% of the test signal is contained in one bin.
Only the n b = 2 category is considered for the gluon-gluon fusion production while both the n b = 2 and n b ≥ 3 categories are included in the likelihood calculation for the b-associated production.
For each bin, S i is calculated from the total integrated luminosity, the theoretical cross-section for the signal and its selection efficiency. The sum of all background contributions in the bin, B i , is estimated from simulation. However, the tt and Z +jets control regions are included in the likelihood calculation as one bin each, to help constrain their respective contributions in the signal regions. This is achieved by introducing two free normalisation  Systematic uncertainties are incorporated in the likelihood as nuisance parameters with either Gaussian or log-normal constraint terms. They include both the experimental and theoretical sources of uncertainty. Experimental uncertainties comprise those in the luminosity measurement, trigger, object identification, energy/momentum scale and resolution as well as underlying event and pile-up modelling. These uncertainties, discussed in Refs. [61,10], impact the simulations of signal and background processes. Theoretical uncertainties include both the signal and background modelling. For the signal modelling, uncertainties due to the factorisation and renormalisation scale choice, the initial-and final-state radiation treatment and the PDF choice are considered. Additional systematic uncertainties are assigned to cover the differences in signal efficiencies and m bb parameter values between the interpolations and the simulations. For the background modelling, the most important sources of systematic uncertainty are from the modelling of the m bb and the p T distributions of Z +jets. They are taken to be the difference between the data and simulation of the selected samples before the event categorisation and the m bb requirement. The samples are dominated by the Z +jets contribution, and any potential signal contamination is expected to be negligible. For other background processes, they are obtained by varying the factorisation and renormalisation scales, the amount of initialand final-state radiation, and the choices of PDF parameterisations.
The effect of these systematic uncertainties on the search is studied using the signal-strength parameter μ for hypothesised signal production. Uncertainties having the largest impact depend on the choice of (m A , m H ) signal point. Table 2 shows the relative uncertainties in the best-fit μ value from the leading sources of systematic uncertainty for two example mass points of both gluon-gluon fusion and b-associated production of a narrow-width A boson. The leading sources of systematic uncertainty are similar for other mass points studied and for larger A boson widths. For all cases, the limited size of the simulated samples has the largest impact on the search sensitivity among all the sources of systematic uncertainty. While systematic uncertainties and the statistical uncertainty of the data have comparable impact at low masses, the search sensitivity is mostly determined at high masses by the limited size of the data sample.

Results
The m bb distributions from different m bb mass windows are scanned for potential excesses beyond the background expectations through signal-plus-background fits. The scan is performed in steps   The last bin is used as a reference for normalisation, and its width is noted in the y-axis label. In this case, the content displayed in a bin is the number of events shown in that bin multiplied by the ratio of widths of the last bin relative to the bin shown.

Conclusion
Data recorded by the ATLAS experiment at the LHC, corresponding to an integrated luminosity of 36.