Search for heavy neutral Higgs bosons produced in association with $b$-quarks and decaying to $b$-quarks at $\sqrt{s}=13$ TeV with the ATLAS detector

A search for heavy neutral Higgs bosons produced in association with one or two $b$-quarks and decaying to $b$-quark pairs is presented using 27.8 fb$^{-1}$ of $\sqrt{s}=13$ TeV proton-proton collision data recorded by the ATLAS detector at the Large Hadron Collider during 2015 and 2016. No evidence of a signal is found. Upper limits on the heavy neutral Higgs boson production cross-section times its branching ratio to $b\bar{b}$ are set, ranging from 4.0 to 0.6 pb at 95% confidence level over a Higgs boson mass range of 450 to 1400 GeV. Results are interpreted within the two-Higgs-doublet model and the Minimal Supersymmetric Standard Model.


Introduction
The measured properties of the Higgs boson discovered at the Large Hadron Collider (LHC) by the ATLAS and CMS collaborations [1, 2] with a mass of 125 GeV are consistent with those of the scalar particle that emerges from the mechanism of electroweak symmetry breaking in the Standard Model (SM) with its one doublet of complex scalar fields [3][4][5][6]. Alternative electroweak symmetry breaking models which contain a scalar particle with properties similar to the SM Higgs boson remain viable, however. A simple, well-studied and well-motivated extension of the mechanism of electroweak symmetry breaking in the SM is the two-Higgs-doublet model (2HDM), which contains two doublets of complex scalar fields [7,8]. In the 2HDM there are, assuming negligible CP-violating effects, two CP-even scalar bosons, h and H which satisfy the mass relation m h < m H , one CP-odd pseudoscalar boson, A, and two electrically charged scalar bosons, H ± . The most general renormalizable, electroweak gauge invariant 2HDM contains tree-level Higgs-boson-mediated flavor-changing neutral currents [8] that are in conflict with experimental limits. When symmetries are imposed to naturally suppress flavor changing neutral currents, four model types emerge, distinguished from one another by their Yukawa couplings, as summarized in Table 1 for h, H, and A.
The agreement of SM predictions with measurements of the 125 GeV Higgs boson, assumed in this paper to be the scalar boson h in the 2HDM, is reducing the 2HDM parameter space towards the alignment limit of cos(β − α) ≈ 0, where tan β is the ratio of the vacuum expectation values of the two scalar doublets and α is the mixing angle of the two CP-even scalar bosons [9]. In the alignment limit, decays of the H and A bosons into gauge boson pairs W + W − and Z Z are heavily suppressed, and the fermion coupling pattern simplifies to that of Table 1. The suppression of H/A couplings to W + W − and Z Z, along with ATLAS and CMS limits on new particle production, implies that searches for the heavy neutral Higgs bosons of the 2HDM mainly rely on their couplings to third-generation fermions. Table 1: Tree-level fermion couplings of the 2HDM h, H, and A bosons for model types I, II, X (or lepton-specific), and Y (or flipped). Here U, D, and E refer to up-type quarks, down-type quarks, and charged leptons, respectively, t β ≡ tan β is the ratio of the vacuum expectation values of the two scalar doublets, and = cos (β − α) where α is the mixing angle of the two CP-even scalar bosons [9]. The couplings are normalized to the SM Higgs boson couplings h SMŪ U, h SMD D, and h SMĒ E and are given in the alignment limit cos (β − α) ≈ 0 where the couplings of the light scalar boson h are close to SM expectations. hŪU hDD hĒ E HŪU HDD HĒ E i AŪγ 5 U i ADγ 5 D i AĒγ 5 E The Higgs sector of the Minimal Supersymmetric Standard Model (MSSM) is a Type II 2HDM, which has motivated searches for heavy neutral Higgs bosons at LEP [10] and the LHC [11,12]. These searches use decays of heavy neutral Higgs bosons into τ + τ − , and are sensitive to Type II and lepton-specific 2HDMs. They are not sensitive to flipped 2HDMs at large tan β, however, and they do not cover the entire MSSM Type II 2HDM parameter space since radiative corrections can significantly increase the ratio of the bb and τ + τ − partial widths beyond the tree-level value of 3m 2 b /m 2 τ [13].  d)). In practice, the optimal balance between signal efficiency and background rejection is achieved by requiring that signal events contain at least three b-quark-initiated jets. The search is performed for neutral Higgs bosons in the mass range 450-1400 GeV. A similar search was performed by the CMS Collaboration for the mass range 300-1300 GeV [15].
The kinematic distributions for the production and decay of H and A bosons are nearly identical, and therefore this search is insensitive to the CP properties of the two heavy neutral Higgs bosons of the 2HDM. The φ boson is be used in this paper to represent the CP-even H boson, the CP-odd A boson, or a Higgs boson mass eigenstate with an arbitrary mixture of CP-even and CP-odd eigenstates.

The ATLAS detector
The ATLAS experiment [16] at the LHC is a multipurpose particle detector with a forward-backward symmetric cylindrical geometry and a near 4π coverage in solid angle.1 It consists of an inner tracking detector surrounded by a thin superconducting solenoid providing a 2 T axial magnetic field, electromagnetic and hadronic calorimeters, and a muon spectrometer. The inner tracking detector covers the pseudorapidity range |η| < 2.5. It consists of silicon pixel, silicon microstrip, and transition radiation tracking detectors. The innermost pixel layer [17,18] was added before the start of collisions in 2015. Lead/liquid-argon (LAr) sampling calorimeters provide electromagnetic (EM) energy measurements with high granularity. A hadronic steel/scintillator-tile calorimeter covers the central pseudorapidity range |η| < 1.7. The endcap and forward regions are instrumented with LAr calorimeters for both the EM and hadronic energy measurements up to |η| = 4.9. The muon spectrometer surrounds the calorimeters and features three large air-core toroid superconducting magnets with eight coils each. The field integral of the toroids ranges from 2.0 to 6.0 T·m across most of the detector. It includes a system of precision tracking chambers and fast 1 ATLAS uses a right-handed coordinate system with its origin at the nominal interaction point (IP) in the center of the detector and the z-axis coinciding with the axis of the beam pipe. The x-axis points from the IP towards the center 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 z-axis. The pseudorapidity is defined in terms of the polar angle θ as η = − ln tan(θ/2). Angular distance is defined as ∆R ≡ (∆η) 2 + (∆φ) 2 . Transverse momentum and energy are defined as p T = p sin θ and E T = E sin θ, respectively. detectors for triggering. A two-level trigger system [19] consisting of the level-1 (L1) trigger, implemented in hardware, and the software-based high-level trigger (HLT), selects interesting events. The L1 trigger uses a subset of detector information to reduce the event rate to a design value of at most 100 kHz. The HLT, which can run offline reconstruction algorithms and is used in this analysis for the triggering of b-quark-initiated jets, reduces the event rate to about 1 kHz.

Data
Proton-proton (pp) collision data recorded by the ATLAS detector at the LHC during 2015 and 2016 at a center-of-mass energy of √ s = 13 TeV were used for the analysis described in this paper. For this data, it is required that the LHC operate with stable beam conditions and that all relevant detector systems be fully functional. The data, corresponding to integrated luminosities of 3.2 ± 0.1 fb −1 and 24.5 ± 0.5 fb −1 for 2015 and 2016, respectively, were collected using a combination of HLT triggers, which employ algorithms [19] to identify jets containing b-hadrons (resulting in 'b-tagged jets'). Maximum-likelihood algorithms were utilized in 2015, while the offline multivariate classifier MV2c20 [20,21] was used in 2016. Events were recorded if they passed the L1 single-jet trigger with a transverse energy (E T ) threshold of E T = 100 GeV, and if the HLT identifies either one b-tagged jet with E T > 225 GeV or two b-tagged jets with different thresholds of E T = 150 GeV and E T = 50 GeV. For the single (double) b-tagged jet trigger, the operating points correspond to a b-quark identification efficiency of 79% (72%) in 2015 and 60% (60%) in 2016, as measured with a reference tt sample. An inefficiency in the online vertex reconstruction affected a fraction of the data collected during 2016; events from these running periods were not included in the analysis. The efficiency of the combination of the two HLT b-jet triggers is shown in Figure 2 for events passing the final selection described in Section 4 and ranges from 80% for m φ = 450 GeV to 95% for m φ > 700 GeV.

Signal model and background simulations
Signal events for the subprocesses bg → bφ + 0, 1 jet, gg → bbφ, and qq → bbφ with φ → bb were generated at leading order (LO) for fifteen m φ values from 450 to 1400 GeV using the S 2.2.0 [22] Monte Carlo (MC) program in the 5FS with the NNPDF30NNLO [23] set of parton distribution functions (PDF). In order to determine the total width at each value of m φ , a specific MSSM scenario tailored for large values of the branching ratio (B) for φ → bb was used in which tan β = 20, the higgsino mass parameter µ = −800 GeV, the generic soft-SUSY-breaking mass parameter M SUSY = 1000 GeV, the trilinear Higgs-top-squark coupling A t = 2000 GeV, the SU(2) gaugino mass parameter M 2 = 800 GeV, and the SU(3) gaugino mass parameter M 3 = 1600 GeV. These parameters suppress φ boson decays into top quark pairs, top-squark pairs, and electroweak gauginos, while decays into pairs of b-quarks are enhanced through MSSM radiative corrections [13]. The FeynHiggs program [24] was used to calculate the branching ratios and the cross-sections shown in Table 2, where the branching ratio B(φ → bb) > 85% for all m φ values up to 1400 GeV. Given the large values for B(φ → bb) in Table 2, the total widths derived from this set of MSSM parameters also represent the typical total widths in the flipped scenario of the 2HDM in the alignment limit for the same m φ and tan β. The values of the total width in Table 2 are much smaller than the 10-15% experimental bb mass resolution. Although several decay modes are present in this MSSM scenario, only the decay mode φ → bb was simulated in the generated signal samples. The ATLAS detector response to the generated signal events was modeled using the ATLAS full simulation software [25] based on G 4 [26]. The impact of multiple pp collisions in the same or nearby bunch crossings (pileup) was simulated by overlaying minimum-bias events on each generated event. The minimum-bias events were generated with P 8.186 [27], using the A2 set of tuned parameters [28] and the MSTW2008LO PDF sets [29]. Finally, events were processed using the same reconstruction Table 2: Mass, total width, and branching ratios of the φ boson of the MSSM scenario used for signal event generation where tan β = 20, µ = −800 GeV, M SUSY = 1000 GeV, A t = 2000 GeV, M 2 = 800 GeV, and M 3 = 1600 GeV. The pp → bbφ cross-section at √ s = 13 TeV is also shown. Full simulation samples of 300,000 events were produced for each of the mass points. The FeynHiggs program [24] was used to calculate the branching ratios and the cross-sections for pp → bbφ at The background estimate is data-driven, as described in Section 5. Background MC samples (referred to later as 'multi-jet MC samples') served as a guide in developing the background model and consist of a S 2.1.1 simulation of multiple b-jets, a next-to-leading-order (NLO) P [30][31][32] simulation of tt production interfaced to P 6.428 [33], and an LO M [34] simulation of the subprocess j j → Z j j, Z → bb interfaced to P 8.205, where j represents a gluon or a u, d, s, c quark/anti-quark. The full ATLAS detector simulation software was used for tt production. A fast ATLAS detector simulation in which the calorimeter response is parameterized [25,35] was used for the multiple b-jets and j j → Z j j samples.
The dominant background comes from the production of multiple b-jets. In the S 2.1.1 simulation of multiple b-jets all 2 → 2, 3, 4 hard subprocesses with at least one b-quark in the final state were generated at LO. Massive cand b-quarks were used to properly simulate gluon splitting into heavy quarks.

Object reconstruction and event selection
Primary vertex candidates [36] are reconstructed from tracks in the inner detector, and the vertex with the highest sum of the squared transverse momenta of all associated tracks is selected as the hard-scatter primary vertex. Jets are reconstructed using the anti-k t algorithm [37] with radius parameter R = 0.4 from topological clusters of energy in the calorimeter calibrated at the electromagnetic scale [38]. Jets are then calibrated using correction factors derived from simulation and data [39]. In order to suppress jets arising from pileup, jets with transverse momentum p T < 60 GeV and |η| < 2.4 are removed if they fail to satisfy a requirement imposed by the multivariate jet vertex tagger (JVT) algorithm [40], where the JVT working point provides a 92% selection efficiency for hard-scatter jets. In addition, events with jets consistent with noise in the calorimeter or non-collision backgrounds are vetoed [41].
Jets containing a b-hadron are identified offline using the MV2c20 multivariate classifier [20,21], which combines information from several algorithms. These algorithms are based on impact parameters of tracks, reconstructed secondary vertices, and a multi-vertex fitter which reconstructs the b → c hadron decay chain. A working point with an average b-tagging efficiency of 70%, as determined using simulated tt events, is chosen. The corresponding misidentification rates for c-jets and jets originating from light (u, d, s) quarks or gluons is 8.2% and 0.3%, respectively. Jets tagged as b-jets receive an additional energy correction to account for the presence of muons in the jet [42].
Event preselection begins by requiring that the event pass the trigger selection and that there be at least three jets with p T > 20 GeV and |η| < 2.4. The leading and second-leading jets (ordered in p T ) are then required to have p T1 > 160 GeV and p T2 > 60 GeV, respectively. The two leading jets must also be b-tagged. Events are considered to be in the signal region, and classified as 'bbb', if there exists at least one additional b-tagged jet. Events with only two b-tagged jets are considered to be in the control region, and are classified as 'bbanti'. For bbb events the 'third jet' is defined to be the third-leading b-tagged jet in p T , while for bbanti events the third jet is the third-leading jet in p T . The final preselection requirement is that the minimum ∆R between the third jet and the two leading jets must be greater than 0.8. This requirement reduces the background from gluon splitting into bb in parton showers and subprocesses such as bg → bg * → bbb.
Events are further classified according to the number of jets. The 3-jet, 4-jet and 5-jet regions are defined, where the last one contains events with five or more jets. A larger number of jets often means that significant final-state radiation (FSR) is present in the φ boson decay, making it more difficult to accurately reconstruct m φ from the two highest-p T jets. Consequently, a categorization based on the number of jets improves experimental sensitivity.
Signal sensitivity is enhanced with a transformation of the kinematic variables p T1 , p T2 , and the invariant mass of the two leading b-tagged jets, m bb . Two-dimensional distributions of p T1 versus m bb and of p T2 versus m bb for events with the bbb classification are displayed in Figure 3. As m φ increases the two high-p T jets from the φ boson decay produce additional FSR, but the jet radius parameter remains fixed at R = 0.4. As a consequence, the reconstructed mass distribution based on the two highest-p T jets smears out and it becomes more difficult to distinguish signal from background. However, since FSR occurs stochastically, the φ boson decays with little or no FSR in a subset of the signal events, and these have reconstructed masses close to the true m φ (bottom row of Figure 3). If these events can be isolated from the others, they offer a chance to improve the sensitivity via improved mass resolution and signal-to-background ratio.
To isolate events with small FSR and good m φ resolution, a principal component analysis (PCA) [43] is performed on the three-dimensional distribution of the variables m bb , p T1 , and p T2 using events drawn from the signal MC sample with the bbb classification following preselection. Separate PCAs are performed for each of the fifteen simulated values of m φ and for each of the three n-jet regions. Upon diagonalization of the covariance matrix for m bb , p T1 , and p T2 , the first, second, and third principal components define the variables m bb , p T1 , and p T2 , respectively. The point (m bb , p T1 , p T2 ) = (0, 0, 0) corresponds to the vector of mean values for m bb , p T1 , and p T2 . Two-dimensional distributions of p T1 versus m bb and of p T2 versus m bb are shown in Figure 4.     Figure 4: Two-dimensional distributions of p T1 versus m bb (left column) and of p T2 versus m bb (right column) for events with the bbb classification following preselection, summed over all three n-jet regions, using the PCA for m φ = 1200 GeV. Plots are shown for data (top row), the multi-jet MC sample (middle row), and the m φ = 1200 GeV signal MC sample (bottom row). The multi-jet plots are normalized to an integrated luminosity of 27.8 fb −1 based on the cross-sections provided by the event generators. The signal plots are normalized to σ(pp → bbφ) × B(φ → bb) = 1 pb. The minimum values of p T1 and p T2 in the final event selection are indicated by horizontal lines. Table 3: The squares of the elements c m b b , c p T1 , and c p T2 of the first principal component eigenvectors for the PCAs for m φ = 600, 900, and 1200 GeV and each n-jet region. The eigenvectors are normalized to unity. The principal component is given by where m bb , p T1 , and p T2 are the mean values for m bb , p T1 , and p T2 , respectively. The transformation into the eigenbasis can provide some physical intuition for the relative contributions of m bb , p T1 , and p T2 .  Table 3 for the 3-jet region. The largest component of m bb comes from m bb , regardless of the mass point. However, at larger values of m φ -where there is greater FSR -the p T1 and p T2 components become more important.
The final event selection requirements are p T1 > −10 GeV and p T2 > −50 GeV, independent of n-jet and m φ . These requirements reduce the background while retaining a large fraction of the signal events in regions of high signal density in (m bb , p T1 , p T2 ) space, as shown in Figure 4. Using m bb instead of m bb as the final discriminant leads to increased sensitivity, which becomes more pronounced with increasing values of m φ .

Statistical analysis
The presence of a signal is tested with a binned maximum-likelihood fit to the data using m bb as the final discriminating variable. For each of the considered mass points, a fit is performed simultaneously over six categories corresponding to all combinations of the three jet-multiplicity regions (3-jet, 4-jet, and 5-jet) and of the two b-tag classifications, bbb and bbanti. The shapes and normalizations for the different categories consist of a sum of signal and background contributions. The shapes and normalizations of the signal distributions are obtained from signal MC samples, with the exception of a global normalization factor representing the primary variable of interest, the heavy Higgs bosons production cross-section times branching ratio σ(pp → bbφ) × B(φ → bb). The shapes and normalizations of the background distributions are determined by data. The background shapes are free to take any form satisfying the constraint that the bbb and bbanti shapes for a specific jet-multiplicity region be identical modulo a second-order polynomial correction factor. The six background normalization factors float freely in the fit.
Based on the multi-jet MC sample, the backgrounds for both the bbb and bbanti regions are dominated by the subprocesses gg → gbb and gg → ggbb (such events enter the bbb region via gluon splitting into bb in the parton showering of one of the final-state gluons). However, the subprocesses gb → bbb and gg → bbbb uniquely provide a small but non-negligible contribution to the bbb background, and the polynomial correction factor accounts for this and any other differences between the bbb and bbanti regions. The m bb distributions for both the bbb and bbanti classifications are plotted in Figure 5 for the multi-jet MC sample along with their ratio. The bbb and bbanti shapes for the 3-jet and 4-jet regions in Figure 5 are nearly identical, while the bbb/bbanti ratio for the 5-jet region appears to have an approximately linear dependence on m bb . Application of the χ 2 probability test and the F-test [44] to the simulated multi-jet m bb distributions over all values of m φ demonstrates that a first-order polynomial is sufficient to describe the ratio of the simulated multi-jet bbb and bbanti background shapes for signal masses m φ < 1200 GeV, while a second-order polynomial is needed for m φ ≥ 1200 GeV.
The data bbb and bbanti shapes, as well as their ratio, after applying the selection on p T1 and p T2 are qualitatively similar to the multi-jet MC distributions, as shown in Figure 6. Applying χ 2 probability and F-test criteria to the data m bb distributions over all values of m φ , it is found that a first-order polynomial is sufficient for the 3-jet region with masses m φ < 1200 GeV and for the 4-jet and 5-jet regions with masses m φ < 800 GeV. For all other jet-category/mass combinations, a second-order polynomial is needed. Potential signal contamination does not affect the results of the tests.
The binned maximum-likelihood fit is performed using the RooFit [45] framework and the H F [46] software tool. A product of Poisson probability terms over the bins of the m bb distributions involving the numbers of data events n i, j and the sum of expected signal and background yields ν i, j in each category i and mass bin j forms the binned likelihood function. It accounts for the effects of floating background normalization and systematic uncertainties and is The index t runs over the two flavor categories bbb and bbanti, k runs over the three jet-multiplicity regions, i = (t, k) runs over the six categories, and p runs over the systematic error nuisance parameters. The boldfaced symbols represent vectors of parameters whereas the symbols of the same type in lightface represent individual parameters (usually containing indices). The template histograms, S i, j , are taken directly from signal simulation for the given mass point and are normalized to the event yields expected for a one picobarn signal. Thus, the signal strength parameter µ, which is common to all categories, is σ(pp → bbφ) × B(φ → bb) in picobarns.
The P and L parameters describe the histogram representations of the m 2 bb and m bb functions, respectively. The normalization parameters for these histograms, Q k and A k , correspond to the parabolic and linear parameters, respectively, for the jet-multiplicity region k. The signal strength parameter µ, the six background normalization parameters N i , the background shape parameters B k, j , the linear parameters A k , and the parabolic parameters Q k , are freely floating parameters in the fit, with the exception that, for a fixed jet-multiplicity region k, the sum over bins j of B k, j is constrained to unity.
The fit model contains nuisance parameters accounting for the statistical uncertainty of the MC signal samples and for systematic variations of the shapes and normalizations of the histogram templates used in the fit, as described in Section 6. The variable β i, j represents the systematic variation in the bin content as a function of the nuisance parameters α p . The nuisance parameters α p are constrained within the allowed systematic variations by the f p (a p |α p , σ p ) terms, where a p are auxiliary measurements and σ p denotes the uncertainty in α p . Individual sources of uncertainties are considered uncorrelated. The statistical uncertainty related to the size of MC signal samples is estimated with a variation of the Barlow-Beeston method [47]. In this analysis, each bin in each category is given a single Poissonconstrained nuisance parameter associated with the signal MC prediction for the number of events entering the bin and the total statistical uncertainty in that bin.

Systematic uncertainties
This analysis relies on the prediction of the shapes and normalizations of the discriminating variable m bb for the searched signal. The signal uncertainties are divided into two categories: experimental and those related to the theoretical modeling of the signal. The background model is validated through statistical analysis of the data and tests utilizing the multi-jet MC sample.

Systematic uncertainties of the background model
A cross-check of the background model is performed by applying the full fit procedure to a small multi-jet MC sample with an equivalent integrated luminosity of 6.8 fb −1 for eight m φ hypotheses. The results are summarized in Table 4. The eight separate fits should be uncorrelated given the 10-15% mass resolution for a heavy Higgs boson. With the assumption of no correlation, the χ 2 per degree of freedom is 1.09 for the eight measurements, indicating that no statistically significant spurious signal is found by the analysis. Given this result, as well as the χ 2 probability and the F-test results of Section 5, the systematic uncertainty arising from the choice of function for the ratio of the bbb and bbanti shape distributions is deemed to be much smaller than the statistical uncertainty and is therefore neglected.

Experimental systematic uncertainties of the signal
The dominant experimental systematic uncertainties are related to the calibration of the b-tagging efficiencies in simulation relative to those measured in data for p T < 300 GeV. They are extrapolated to p T > 300 GeV using MC simulation taking into account uncertainties in the jet modeling and detector response affecting b-tagging performance. This calibration is performed separately for b-jets, c-jets, and light-flavor jets and as a function of jet p T and |η| [20]. Uncertainties in the cross-sections for processes used in the b-tagging calibration, modeling of the jet kinematics and flavor composition in the simulated signal samples, detector simulation, and event reconstruction are included [48][49][50]. These uncertainties are decomposed into uncorrelated components resulting in three components for b-jets and c-jets, and five components for light-flavor jets.
Simulation-to-data efficiency differences are also corrected for the trigger, specifically for b-tagged jets.
Since the background estimation is data-driven, this scaling only affects signal. Scale factors obtained by comparing data and simulated dilepton tt events, which are enriched in b-jets, are used to correct simulation-to-data efficiency differences in the b-jet trigger for p T < 240 GeV. For p T > 240 GeV, due to the limited size of the tt data sample, extrapolation based on simulation is used. The systematic uncertainties in these scale factors include mismodeling of the fraction of b-jets in simulation, mismodeling of the b-jet trigger efficiency for non-b-jets, simulation statistical uncertainty, data statistical uncertainty for p T < 240 GeV, uncertainty in the simulation-based extrapolation to p T > 240 GeV, and uncertainties in the dependence of the b-jet trigger efficiency on jet p T and η. The b-jet trigger was calibrated relative to a set of offline b-tagging operating points to correctly take into account correlations between the b-jet trigger and offline b-tagging.
Systematic uncertainties in the jet energy scale and jet energy resolution are based on measurements with data [39,51]. All sources of the jet energy scale uncertainty are decomposed into 21 uncorrelated components that are treated as independent.
The uncertainty in the combined 2015+2016 integrated luminosity is 2.1%. It is derived from the calibration of the luminosity scale using x-y beam-separation scans, following a methodology similar to that detailed in Ref. [52], and using the LUCID-2 detector for the baseline luminosity measurements [53].

Theoretical systematic uncertainties of the signal
The uncertainty related to the choice of generator for the signal hard process and showering model is estimated by comparing the nominal sample with the one obtained by reweighting the nominal sample to the NLO generator M 5_aMC@NLO [54,55] with a 4 flavor scheme (4FS) PDF interfaced to the P 8.205 parton showering model. The particle-level Higgs boson mass, Higgs boson p T , and the p T values of the two leading b-tagged jets are used for sequential reweighting. The uncertainty from the PDF set used in the nominal sample is computed using the standard-deviation method described in Ref. [23].

Results
The search for a single heavy neutral Higgs boson φ produced in association with b-quarks and decaying into a bb pair shows no significant excess above the SM background for any of the analyzed mass points. The post-fit bbb category distributions of the rotated bb invariant mass m bb are shown in Figure 7 together with the m bb distribution extracted from the bbanti region (pre-fit background).

Cross-section limits
Since no significant excess over the background expectation is observed, upper limits on the production of a single heavy neutral Higgs boson φ decaying into a bb pair are set. Figure 8 presents the observed and expected limits for σ(pp → bbφ) × B(φ → bb) at the 95% confidence level (CL). The limits are calculated with the CL s method [56].
The leading sources of uncertainty affecting the best-fit value of σ × B for two of the mass points, 600 and 1200 GeV, are given in Table 5, together with their relative importance. The impact of the given source of uncertainty is obtained by first fixing all the nuisance parameters related to other systematic uncertainties to their best-fit values and then allowing only the nuisance parameters ascribed to the considered source of uncertainty to float in the fit. The uncertainty is dominated by the statistical error, which improves significantly between m φ = 600 and 1200 GeV due to the sharp drop in the background level. The systematic uncertainties with the largest impact on the sensitivity are related to the flavor-tagging calibration of the offline b-tagging algorithm and b-jet trigger and to jet reconstruction.

Model interpretations
The two 2HDM scenarios with enhanced pp → bbφ production and φ → bb decay at large tan β are Type II and Type Y (flipped). The most commonly analyzed scenario is Type II since the Higgs sector of the MSSM is a Type II 2HDM. The results of this search are interpreted in the context of the MSSM for the hMSSM scenario [57] and for the m mod+ h and m mod− h scenarios [58]. The Higgs boson production cross-sections and branching ratios are calculated using the procedures outlined in the LHC Higgs Cross-section Working Group report [59]. The cross-sections for Higgs boson production through bb fusion [14] are determined by matching the 5FS [60,61] and 4FS [62,63] cross-section calculations. For the hMSSM scenario, the Higgs boson masses and branching ratios are calculated using HD [64,65]. For the m mod+ h and m mod− h scenarios, the Higgs boson masses and couplings are calculated with F H [24,[66][67][68][69], and the branching ratios are calculated by combining the most precise results from F H , HD , and PROPHECY4f [70,71].
The 95% CL exclusion limits on tan β as a function of m A are shown in Figure 9 for the hMSSM benchmark. The hMSSM scenario is well defined and broadly representative of the remaining parameter space with SUSY partners too heavy for direct detection at the LHC. As an indication of the sensitivity variation for different MSSM scenarios, Figure 9 also displays the expected sensitivities for the m mod+ [GeV] φ m 600 800 1000 1200 1400

Conclusions
A search for heavy neutral Higgs bosons produced in association with at least one b-quark and decaying into a pair of b-quarks was performed using 27.8 fb −1 of 13 TeV pp collision data recorded by the ATLAS detector at the LHC in 2015 and 2016. The data are compatible with SM expectations, yielding no significant excess of events in the mass range 450-1400 GeV. Upper limits on the cross-section times branching ratio were derived as a function of the mass of the heavy Higgs boson. The 95% CL upper limits are in the range 0.6-4.0 pb. Compared to heavy neutral Higgs boson searches utilizing φ → τ + τ − or A → Z h decays, these limits expand the excluded Type Y (flipped) 2HDM parameter space into regions with | cos(β − α)| ≈ 0 and tan β 20.