Search for a heavy Higgs boson decaying into a $Z$ boson and another heavy Higgs boson in the $\ell\ell bb$ final state in $pp$ collisions at $\sqrt{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 $\sqrt{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% confidence-level upper limits on the $pp\rightarrow A\rightarrow ZH$ production cross-section times the branching ratio $H\rightarrow 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 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 the two-Higgs-doublet model.


Introduction
After the discovery of a Higgs boson at the Large Hadron Collider (LHC) [1,2], one of the most important remaining questions is whether the recently discovered particle is part of an extended scalar sector or not. Additional Higgs bosons appear in all models with an extended scalar sector, such as the two-Higgsdoublet model (2HDM) [3,4]. Such extensions are motivated by, and included in, several new physics scenarios, such as supersymmetry [5], dark matter [6] and axion [7] models, electroweak baryogenesis [8] and neutrino mass models [9].
The addition of a second Higgs doublet leads to five Higgs bosons after electroweak symmetry breaking. The phenomenology of such a model is very rich and depends on many parameters, such as the ratio of the vacuum expectation values of the two Higgs doublets (tan β), and the Yukawa couplings of the scalar sector [4]. When CP conservation is assumed, the model contains two CP-even Higgs bosons, h and H with m H > m h , one CP-odd, A, and two charged scalars, H ± . There have been many searches for the heavy neutral Higgs bosons of the 2HDM at the LHC, including H → WW/Z Z [10][11][12][13], A/H → ττ/bb [14][15][16], A → Z h [17,18] and H → hh [19,20]. For the interpretation of these searches it is usually assumed that the heavy Higgs bosons, H and A, are degenerate in mass, i.e. m A = m H .
This assumption of mass degeneracy is relaxed in this Letter by assuming m A > m H . Such a case is motivated by electroweak baryogenesis scenarios [8] in the context of the 2HDM [21,22]. For 2HDM electroweak baryogenesis to occur [8,23], the requirement m A > m H is necessary for a strong first-order phase transition to take place in the early universe. The A boson mass is also bounded from above to be less than approximately 800 GeV, whereas the lighter CP-even Higgs boson, h, is required to have properties similar to those of a Standard Model (SM) Higgs boson and is assumed to be the Higgs boson with mass of 125 GeV that was discovered at the LHC. Under such conditions and for large parts of the 2HDM parameter space, the CP-odd Higgs boson, A, decays into Z H [24]. The production of the A boson in the relevant 2HDM parameter space proceeds mainly through gluon-gluon fusion and b-associated production at the LHC.
This search for A → Z H decays uses proton-proton collision data at √ s = 13 TeV corresponding to an integrated luminosity of 36.1 fb −1 recorded by the ATLAS detector at the LHC. The search considers only Z → , where = e, µ, to take advantage of the clean leptonic final state, and H → bb, because of its large branching ratio. This final state allows full reconstruction of the A boson's decay kinematics. The reconstruction of the A boson's invariant mass uses the assumed value of the mass of the H boson to improve its resolution. The final state is also categorised by the presence of two or three b-tagged jets to take advantage of the b-associated production mechanism. The CMS Collaboration has published a similar search at √ s = 8 TeV [25]. This Letter reports the result of a search at √ s = 13 TeV, which extends the 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. [26]. 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 [27] 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 [28] 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 data-taking 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 datataking 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 gluon-gluon fusion were generated at leading order with M G 5_aMC@NLO 2.3.3 [29,30] using P 8.210 [31] with a set of tuned parameters called the A14 tune [32] for parton showering. For the generation of A bosons produced in association with b-quarks, M G 5_aMC@NLO 2.1.2 [30,33,34] was used following Ref.
[35] together with P 8.212 and the A14 tune for parton showering. The gluon-gluon fusion production used NNPDF2. 3LO [36] as the parton distribution functions (PDF), while the b-associated production used CT10nlo_nf4 [37]. 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 narrow-width 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 S 2.2.1 [38] using the NNPDF3.0NNLO PDF set [39]. Top-quark-pair production was simulated with P -B v2 [40-42] and the CT10nlo PDF set [37], while the electroweak single-top-quark production was simulated with P -B v1 and the fixed four-flavour PDF set CT10nlo_f4 [37]. The parton shower was performed with P 6.428 [43] using the Perugia 2012 set of tuned parameters [44]. The production of top-quark pairs in association with a vector boson was simulated using M G 5_aMC@NLO 2.2.3 and the NNPDF3.0NLO PDF set, whereas P 8.186 was used for the parton shower with the A14 tune. Production of WW, Z Z and W Z pairs was simulated using S 2.2.1 and the NNPDF3.0NNLO PDF set. Finally, SM Higgs boson production in association with a Z boson was generated with P -B v2 and the NNPDF3.0NLO PDF set, whereas the parton shower was performed with P 8.186 using the AZNLO tune [45].
The modelling of bottom-and charm-hadron decays was performed with the EvtGen v1.2.0 package [46] for all samples apart from those simulated with S . The simulated events were overlaid with minimum-bias events, to account for the effect of multiple interactions occurring in the same and neighbouring bunch crossings ('pile-up'). These events were generated using P 8 with the A2 tune [47] 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 θ. and the MSTW2008LO PDF set [48]. 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 G 4-based [49] detector simulation [50] of the ATLAS detector. The ATLFAST2 simulation [50] 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.

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 [63] 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. Events with additional muons or electrons are rejected.
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: Σp 2 T /m bb > 0.4, 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 soft p T spectrum of the associated b-jets.
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 four-momentum 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%. 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 The dominant backgrounds after these selections are from Z+jets and top-quark production. For topquark-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 opposite-charge 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-to-next-to-leading-order (NNLO) cross-sections [64-67]. Single-top-quark production and top-quark-pair production in association with vector bosons are normalised to next-toleading-order (NLO) cross-sections from Refs. [68-70] and Ref. [30], respectively. The normalisation of the Higgs boson production in association with a vector boson follows the recommendations of Ref.

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 ExpGaussExp (EGE) function [71]: 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 [72]: 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 m A . The core width, σ, is observed to monotonically increase with ∆m ≡ m A − m H and is parameterised with a third-degree polynomial. The rest of the parameters are largely constant and are fixed to their average values from the fits, with the exception of mass points with ∆m = 100 GeV. The distributions at mass points with ∆m = 100 GeV correspond to the smallest mass splitting considered in this search and are close to the kinematic cutoff. Their non-core parameters are fixed to the average fit values obtained from signal samples with this mass splitting only. As an example of the performance of this procedure, Figure 1 shows a comparison for the (m A , m H ) = (500, 250) GeV mass point between the simulated distributions and the parametric functions described previously. The m bb distributions are well parameterised by the chosen functional forms.
The previously described parameterisation is adequate for 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 distribution3 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 Figure 2 for the same signal points used in Figure 1. Finally, the signal efficiencies for the interpolated mass points are obtained through separate twodimensional interpolations on the (m A , m H ) plane using thin plate splines [73].

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 described in Refs. [74][75][76]. For a given mass hypothesis of (m A , m H ), the likelihood is constructed as the product of Poisson statistics in m bb bins: Here 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. [10,60], 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 differences between the data and simulation in the Z+jets control regions. For other background processes, they are obtained by varying the factorisation and renormalisation scales, the amount of initial-and 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 of 10 GeV for both the m A range 230-800 GeV and the m H range 130-700 GeV, such that m A − m H ≥ 100 GeV. The step sizes are chosen to be compatible with the detector resolution for m bb and m bb . In the absence of a statistically significant excess, constraints on the production of A → Z H followed by the H → bb decay are derived. The method of Ref. [78] is used to calculate 95% confidence level (CL) upper bounds on the product of cross-section and decay branching ratios, σ × B(A → Z H) × B(H → bb), using the asymptotic approximation [76]. The upper limits are shown in Figure 5 for a narrow-width A boson produced via gluon-gluon fusion and b-associated production. As for the significance calculations above, these limits are derived separately for each production process. For the gluon-gluon fusion limits, only the n b = 2 category is used. For the b-associated production, both the n b = 2 and n b ≥ 3 categories are used.  [79][80][81][82]. For b-associated production a cross-section in the four-flavour scheme is also calculated as described in Refs. [83,84] and the results are combined with the five-flavour scheme calculation following Ref. [85]. The Higgs boson widths and branching ratios are calculated using 2HDMC [86]. The procedure for the cross-section and branching ratio, as well as for the choice of 2HDM parameters, follows Ref. [35].
Since both gluon-gluon fusion and b-associated production are expected, a new signal model weighted by the predicted cross sections of the two processes is built for every point tested in the 2HDM parameter space. Upper limits on σ × B(A → Z H) × B(H → bb) with σ here including contributions from both processes are recalculated and compared with the 2HDM predictions to derive the limits in the 2HDM parameter space. Figure 6 shows the observed and expected limits for Type I, Type II, 'lepton specific' and 'flipped' 2HDMs in the (m A , m H ) plane for various tan β values. Type I and lepton specific 2HDMs share the same couplings at tan β = 1, which leads to an exclusion up to m H = 350 GeV for this tan β value. For tan β > 1, however, the H → ττ decay is preferred in the lepton specific model, resulting in much lower sensitivity compared to Type I. In Type II and flipped 2HDMs the exclusions can reach m H = 400 GeV at lower tan β (less than 10) and m H = 600 GeV at higher tan β (more than 20), due to the enhancement of the couplings of H boson to down-type quarks in this model.

Conclusion
Data recorded by the ATLAS experiment at the LHC, corresponding to an integrated luminosity of 36.1 fb −1 from proton-proton collisions at a centre-of-mass energy 13 TeV, are used to search for a heavy Higgs boson, A, decaying into Z H, where H denotes a heavy Higgs boson with mass m H > 125 GeV. The A boson is assumed to be produced via either gluon-gluon fusion or b-associated production. No significant deviation from the SM background predictions are observed in the Z H → bb final state that is considered in this search. Considering each production process separately, upper limits are set at the 95% confidence level for σ × B(A → Z H) × B(H → bb) of 14-830 fb for gluon-gluon fusion and 26-570 fb for b-associated production of a narrow A boson for the mass ranges 130-700 GeV of the H boson and 230-800 GeV of the A boson. Taking into account both production processes, this search tightens the constraints on the 2HDM in the case of large mass splittings between its heavier neutral Higgs bosons. [