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

TL;DR

This paper reports a search for heavy neutral Higgs bosons produced with b-quarks and decaying into b-quark pairs at 13 TeV, setting upper limits on production cross-sections and interpreting results within theoretical models.

Contribution

First search for heavy neutral Higgs bosons with b-quark associated production and decay into b-quark pairs at 13 TeV with ATLAS, providing new constraints on models.

Findings

No evidence of Higgs signal found.

Upper limits on cross-section range from 0.6 to 4.0 pb.

Results constrain two-Higgs-doublet and MSSM models.

Abstract

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.

Tables

Table 1. Table 1 : Tree-level fermion couplings of the 2HDM h ℎ h , H 𝐻 H , and A 𝐴 A bosons for model types I, II, X (or lepton-specific), and Y (or flipped). Here U 𝑈 U , D 𝐷 D , and E 𝐸 E refer to up-type quarks, down-type quarks, and charged leptons, respectively, t β ≡ tan ⁡ β subscript 𝑡 𝛽 𝛽 t_{\beta}\equiv\tan\beta is the ratio of the vacuum expectation values of the two scalar doublets, and ϵ = cos ⁡ ( β − α ) italic-ϵ 𝛽 𝛼 \epsilon=\cos{(\beta-\alpha)} where α 𝛼 \alpha 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 ¯ ​ U subscript ℎ SM ¯ 𝑈 𝑈 h_{\mathrm{SM}}\bar{U}U , h SM ​ D ¯ ​ D subscript ℎ SM ¯ 𝐷 𝐷 h_{\mathrm{SM}}\bar{D}D , and h SM ​ E ¯ ​ E subscript ℎ SM ¯ 𝐸 𝐸 h_{\mathrm{SM}}\bar{E}E and are given in the alignment limit cos ⁡ ( β − α ) ≈ 0 𝛽 𝛼 0 \cos{(\beta-\alpha)\approx 0} where the couplings of the light scalar boson h ℎ h are close to SM expectations.

	$h\overline{U}U$	$h\overline{D}D$	$h\overline{E}E$	$H\overline{U}U$	$H\overline{D}D$	$H\overline{E}E$	$iA\overline{U}{\gamma }_{{}_5}U$	$iA\overline{D}{\gamma }_{{}_5}D$	$iA\overline{E}{\gamma }_{{}_5}E$
I	$1+\frac{ϵ}{t_{\beta }}$	$1+\frac{ϵ}{t_{\beta }}$	$1+\frac{ϵ}{t_{\beta }}$	$-\left(\frac{1}{t_{\beta }}-ϵ\right)$	$-\left(\frac{1}{t_{\beta }}-ϵ\right)$	$-\left(\frac{1}{t_{\beta }}-ϵ\right)$	$-\frac{1}{t_{\beta }}$	$\frac{1}{t_{\beta }}$	$\frac{1}{t_{\beta }}$
II	$1+\frac{ϵ}{t_{\beta }}$	$1-ϵt_{\beta }$	$1-ϵt_{\beta }$	$-\left(\frac{1}{t_{\beta }}-ϵ\right)$	$t_{\beta }+ϵ$	$t_{\beta }+ϵ$	$-\frac{1}{t_{\beta }}$	$-t_{\beta }$	$-t_{\beta }$
X	$1+\frac{ϵ}{t_{\beta }}$	$1+\frac{ϵ}{t_{\beta }}$	$1-ϵt_{\beta }$	$-\left(\frac{1}{t_{\beta }}-ϵ\right)$	$-\left(\frac{1}{t_{\beta }}-ϵ\right)$	$t_{\beta }+ϵ$	$-\frac{1}{t_{\beta }}$	$\frac{1}{t_{\beta }}$	$-t_{\beta }$
Y	$1+\frac{ϵ}{t_{\beta }}$	$1-ϵt_{\beta }$	$1+\frac{ϵ}{t_{\beta }}$	$-\left(\frac{1}{t_{\beta }}-ϵ\right)$	$t_{\beta }+ϵ$	$-\left(\frac{1}{t_{\beta }}-ϵ\right)$	$-\frac{1}{t_{\beta }}$	$-t_{\beta }$	$\frac{1}{t_{\beta }}$

Table 2. Table 2 : Mass, total width, and branching ratios of the ϕ italic-ϕ \phi boson of the MSSM scenario used for signal event generation where tan ⁡ β = 20 𝛽 20 \tan\beta=20 , μ = − 800 𝜇 800 \mu=-800   GeV , M SUSY = 1000 subscript 𝑀 SUSY 1000 M_{\textrm{SUSY}}=1000   GeV , A t = 2000 subscript 𝐴 𝑡 2000 A_{t}=2000   GeV , M 2 = 800 subscript 𝑀 2 800 M_{2}=800   GeV , and M 3 = 1600 subscript 𝑀 3 1600 M_{3}=1600   GeV . The p ​ p → b ​ b ¯ ​ ϕ → 𝑝 𝑝 𝑏 ¯ 𝑏 italic-ϕ pp\rightarrow b\bar{b}\phi cross section at s = 13 𝑠 13 \sqrt{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 p ​ p → b ​ b ¯ ​ ϕ → 𝑝 𝑝 𝑏 ¯ 𝑏 italic-ϕ pp\rightarrow b\bar{b}\phi at s = 13 𝑠 13 \sqrt{s}=13   TeV .

$m_{\varphi }$	${\Gamma }_{\varphi }\left(\text{total}\right)$	$ℬ\left(\varphi \to b\overline{b}\right)$	$ℬ\left(\varphi \to {\tau }^{+}{\tau }^{-}\right)$	$ℬ\left(\varphi \to t\overline{t}\right)$	$\sigma \left(b\overline{b}\varphi \right)$
$[$GeV$]$	$[$GeV$]$				$[$fb$]$
450	4.7	0.91	0.07	0.02	2792
500	5.1	0.91	0.07	0.02	1694
550	5.5	0.91	0.07	0.02	1066
600	5.9	0.91	0.07	0.02	693
650	6.3	0.91	0.07	0.02	463
700	6.7	0.91	0.07	0.02	317
750	7.1	0.90	0.07	0.02	222
800	7.5	0.90	0.08	0.02	158
850	7.9	0.90	0.08	0.02	115
900	8.3	0.90	0.08	0.02	85
950	8.7	0.90	0.08	0.02	66
1000	9.1	0.90	0.08	0.02	48
1100	9.8	0.90	0.08	0.02	29
1200	10.7	0.89	0.08	0.02	18
1300	11.6	0.87	0.08	0.02	12
1350	12.1	0.87	0.08	0.02	10
1400	12.5	0.86	0.08	0.02	8

Table 3. Table 3 : The squares of the elements c m b ​ b subscript 𝑐 subscript 𝑚 𝑏 𝑏 c_{m_{bb}} , c p T1 subscript 𝑐 subscript 𝑝 T1 c_{p_{\textrm{T1}}} , and c p T2 subscript 𝑐 subscript 𝑝 T2 c_{p_{\textrm{T2}}} of the first principal component eigenvectors for the PCAs for m ϕ subscript 𝑚 italic-ϕ m_{\phi} = 600, 900, and 1200  GeV  and each n 𝑛 n -jet region. The eigenvectors are normalized to unity. The principal component is given by m b ​ b ′ = c m b ​ b ​ ( m b ​ b − ⟨ m b ​ b ⟩ ) + c p T1 ​ ( p T1 − ⟨ p T1 ⟩ ) + c p T2 ​ ( p T2 − ⟨ p T2 ⟩ ) superscript subscript 𝑚 𝑏 𝑏 ′ subscript 𝑐 subscript 𝑚 𝑏 𝑏 subscript 𝑚 𝑏 𝑏 delimited-⟨⟩ subscript 𝑚 𝑏 𝑏 subscript 𝑐 subscript 𝑝 T1 subscript 𝑝 T1 delimited-⟨⟩ subscript 𝑝 T1 subscript 𝑐 subscript 𝑝 T2 subscript 𝑝 T2 delimited-⟨⟩ subscript 𝑝 T2 m_{bb}^{\prime}=c_{m_{bb}}(m_{bb}-\langle m_{bb}\rangle)+c_{p_{\textrm{T1}}}(p_{\textrm{T1}}-\langle p_{\textrm{T1}}\rangle)+c_{p_{\textrm{T2}}}(p_{\textrm{T2}}-\langle p_{\textrm{T2}}\rangle) , where ⟨ m b ​ b ⟩ delimited-⟨⟩ subscript 𝑚 𝑏 𝑏 \langle m_{bb}\rangle , ⟨ p T1 ⟩ delimited-⟨⟩ subscript 𝑝 T1 \langle p_{\textrm{T1}}\rangle , and ⟨ p T2 ⟩ delimited-⟨⟩ subscript 𝑝 T2 \langle p_{\textrm{T2}}\rangle are the mean values for m b ​ b subscript 𝑚 𝑏 𝑏 m_{bb} , p T1 subscript 𝑝 T1 p_{\textrm{T1}} , and p T2 subscript 𝑝 T2 p_{\textrm{T2}} , respectively. The transformation into the eigenbasis can provide some physical intuition for the relative contributions of m b ​ b subscript 𝑚 𝑏 𝑏 m_{bb} , p T1 subscript 𝑝 T1 p_{\textrm{T1}} , and p T2 subscript 𝑝 T2 p_{\textrm{T2}} .

	3-jet			4-jet			5-jet
$m_{\varphi }$ $[$GeV$]$	$c_{m_{bb}}^2$	$c_{p_{\text{T1}}}^2$	$c_{p_{\text{T2}}}^2$	$c_{m_{bb}}^2$	$c_{p_{\text{T1}}}^2$	$c_{p_{\text{T2}}}^2$	$c_{m_{bb}}^2$	$c_{p_{\text{T1}}}^2$	$c_{p_{\text{T2}}}^2$
600	0.74	0.12	0.14	0.76	0.12	0.13	0.82	0.10	0.08
900	0.55	0.20	0.26	0.64	0.16	0.20	0.74	0.14	0.13
1200	0.45	0.25	0.30	0.59	0.20	0.21	0.72	0.15	0.14

Table 4. Table 4 : Best-fit values for σ ​ ( p ​ p → b ​ b ¯ ​ ϕ ) × ℬ ​ ( ϕ → b ​ b ¯ ) 𝜎 → 𝑝 𝑝 𝑏 ¯ 𝑏 italic-ϕ ℬ → italic-ϕ 𝑏 ¯ 𝑏 \sigma(pp\rightarrow b\bar{b}\phi)\times\mathcal{B}(\phi\rightarrow b\bar{b}) in the multijet MC sample with a total uncertainty Δ ​ ( σ × ℬ ) Δ 𝜎 ℬ \Delta(\sigma\times\mathcal{B}) for each mass point. The multijet MC sample has an equivalent integrated luminosity of 6.8 fb -1 .

Mass $[$GeV$]$	$\sigma \times ℬ$ $[$pb$]$	$\Delta \left(\sigma \times ℬ\right)$ $[$pb$]$
450	–0.99	2.42
550	0.77	1.44
650	–0.75	0.80
750	0.42	0.62
850	–0.31	0.50
1000	–0.24	0.49
1200	0.71	0.35
1400	–0.29	0.20

Table 5. Table 5 : Grouped contributions to the systematic uncertainty of the best-fit value of σ × ℬ 𝜎 ℬ \sigma\times\mathcal{B} . The best-fit values for m ϕ = 600 subscript 𝑚 italic-ϕ 600 m_{\phi}=600 and 1200  GeV  are 0.76 and − 0.1 0.1 -0.1  pb, respectively.

	$m_{\varphi }=600$ GeV	$m_{\varphi }=1200$ GeV
Source of uncertainty	$\Delta \left(\sigma \times ℬ\right)$ [pb]	$\Delta \left(\sigma \times ℬ\right)$ [pb]
Total	0.80	0.29
Statistical	0.77	0.26
Systematic	0.20	0.11
Experimental uncertainties
Jet-related	0.05	0.05
Flavor-tagging	0.12	0.05
Trigger	0.04	0.05
Luminosity	0.02	0.01
Theoretical and modeling uncertainties
Generator	0.03	0.03
PDF	0.08	0.04
MC statistical	0.09	0.04

Equations

$$P(\mathbf{n},\mathbf{a}|\mu ,\mathbf{N},\mathbf{\gamma },\mathbf{\alpha })=\prod _{i\in \text{categories}}\prod _{j\in \text{bins}}\mathrm{Pois}(n_{i,j}|{\nu }_{i,j})\prod _{p\in \text{sys. nuis. param.}}{f}_p(a_p|{\alpha }_p,{\sigma }_p),$$

$${\nu }_{i,j}$$

$${\gamma }_{i,j}$$

Full text

1. 1 Introduction

2. 2 The ATLAS detector

3. 3 Data and simulated samples

4. 3.1 Data

5. 3.2 Signal model and background simulations

6. 4 Object reconstruction and event selection

7. 5 Statistical analysis

8. 6 Systematic uncertainties

9. 6.1 Systematic uncertainties of the background model

10. 6.2 Experimental systematic uncertainties of the signal

11. 6.3 Theoretical systematic uncertainties of the signal

12. 7 Results

13. 7.1 Cross section limits

14. 7.2 Model interpretations

15. 8 Conclusions

1 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^{\pm}$. 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 toward the alignment limit of $\cos(\beta-\alpha)\approx 0$, where $\tan\beta$ is the ratio of the vacuum expectation values of the two scalar doublets and $\alpha$ 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 $ZZ$ are heavily suppressed, and the fermion coupling pattern simplifies to that of Table 1. The suppression of $H/A$ couplings to $W^{+}W^{-}$ and $ZZ$, 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.

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 $\tau^{+}\tau^{-}$, and are sensitive to Type II and lepton-specific 2HDMs. They are not sensitive to flipped 2HDMs at large $\tan\beta$, however, and they do not cover the entire MSSM Type II 2HDM parameter space since radiative corrections can significantly increase the ratio of the $b\bar{b}$ and $\tau^{+}\tau^{-}$ partial widths beyond the tree-level value of $3m_{b}^{2}/m_{\tau}^{2}$ [13].

This paper presents a search for heavy neutral Higgs bosons produced in association with one or two $b$-quarks and decaying into $b$-quark pairs using 27.8 fb*-1* of $\sqrt{s}=13$ TeV proton-proton collision data recorded by the ATLAS detector at the LHC during 2015 and 2016. The search is sensitive to the Type II and flipped scenarios of the 2HDM in the regime where $\tan\beta\gg 1$. In the 5-flavor scheme (5FS) [14], processes such as those shown in Figure 1 lead to the production of heavy neutral Higgs bosons in association with one $b$-quark (Figures 1(a) and 1(b)) or two $b$-quarks (Figures 1(c) and 1(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 $\phi$ boson is 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.

2 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\pi$ coverage in solid angle.111 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 toward the center of the LHC ring, and the $y$-axis points upward. Cylindrical coordinates ($r$,$\phi$) are used in the transverse plane, $\phi$ being the azimuthal angle around the $z$-axis. The pseudorapidity is defined in terms of the polar angle $\theta$ as $\eta=-\ln\tan(\theta/2)$. Angular distance is defined as $\Delta R\equiv\sqrt{(\Delta\eta)^{2}+(\Delta\phi)^{2}}$. Transverse momentum and energy are defined as $p_{\text{T}}=p\sin\theta$ and $E_{\text{T}}=E\sin\theta$, respectively. It consists of an inner tracking detector surrounded by a thin superconducting solenoid providing a $2\text{\,}\mathrm{T}$ axial magnetic field, electromagnetic and hadronic calorimeters, and a muon spectrometer. The inner tracking detector covers the pseudorapidity range $|\eta|<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 $|\eta|<1.7$. The endcap and forward regions are instrumented with LAr calorimeters for both the EM and hadronic energy measurements up to $|\eta|=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 across most of the detector. It includes a system of precision tracking chambers and fast 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.

3 Data and simulated samples

3.1 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 $\sqrt{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\pm 0.1$ fb*-1* and $24.5\pm 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_{\textrm{T}}$) threshold of $E_{\textrm{T}}=100$ GeV, and if the HLT identifies either one $b$-tagged jet with $E_{\textrm{T}}>225$ GeV or two $b$-tagged jets with different thresholds of $E_{\textrm{T}}=150$ GeV and $E_{\textrm{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 $t\bar{t}$ 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_{\phi}=450$ GeV to 95% for $m_{\phi}>700$ GeV. The efficiency of the double $b$-tagged jet trigger falls with increasing $m_{\phi}$ since the per-jet $b$-tagging efficiency is dropping with increasing $E_{\textrm{T}}$. The single $b$-tagged jet trigger efficiency increases with increasing $m_{\phi}$ because a larger fraction of the second leading jets satisfy $E_{\textrm{T}}>225$ GeV.

3.2 Signal model and background simulations

Signal events for the subprocesses $bg\to b\phi+0,1$ jet, $gg\to b\bar{b}\phi$, and $q\bar{q}\to b\bar{b}\phi$ with $\phi\to b\bar{b}$ were generated at leading order (LO) for seventeen $m_{\phi}$ values from 450 to 1400 GeV using the Sherpa 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_{\phi}$, a specific MSSM scenario tailored for large values of the branching ratio ($\mathcal{B}$) for $\phi\rightarrow b\bar{b}$ was used in which $\tan\beta=20$, the higgsino mass parameter $\mu=-800$ GeV, the generic soft-SUSY-breaking mass parameter $M_{\textrm{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 $\phi$ 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 $\mathcal{B}(\phi\rightarrow b\bar{b})>85\%$ for all $m_{\phi}$ values up to 1400 GeV. Given the large values for $\mathcal{B}(\phi\rightarrow b\bar{b})$ 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_{\phi}$ and $\tan\beta$. The values of the total width in Table 2 are much smaller than the 10%–15% experimental $b\bar{b}$ mass resolution. Although several decay modes are present in this MSSM scenario, only the decay mode $\phi\to b\bar{b}$ is simulated in the generated signal samples. Since they are ignored in the analysis, additional $\phi$ decay modes that happen to leak into the $\phi\to b\bar{b}$ acceptance will only make any limits on $\phi\to b\bar{b}$ more conservative.

The ATLAS detector response to the generated signal events was modeled using the ATLAS full simulation software [25] based on Geant4 [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 Pythia 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 software as in data.

The background estimate is data-driven, as described in Section 5. Background MC samples (referred to later as ‘multijet MC samples’) served as a guide in developing the background model and consist of a Sherpa 2.1.1 simulation of multiple $b$-jets, a next-to-leading-order (NLO) Powheg [30, 31, 32] simulation of $t\bar{t}$ production interfaced to Pythia 6.428 [33], and an LO Madgraph [34] simulation of the subprocess $jj\to Z\ jj$, $Z\to b\bar{b}$ interfaced to Pythia 8.205, where $j$ represents a gluon or a $u,d,s,c$ quark/antiquark. The full ATLAS detector simulation software was used for $t\bar{t}$ production. A fast ATLAS detector simulation in which the calorimeter response is parametrized [25, 35] was used for the multiple $b$-jets and $jj\to Z\ jj$ samples.

The dominant background comes from the production of multiple $b$-jets. In the Sherpa 2.1.1 simulation of multiple $b$-jets all $2\rightarrow 2,3,4$ hard subprocesses with at least one $b$-quark in the final state were generated at LO. The $c$- and $b$-quarks masses are set to their running Yukawa values to properly simulate gluon splitting into heavy quarks.

4 Object reconstruction and event selection

Primary vertex candidates [36] are reconstructed using 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_{\textrm{T}}<60$ GeV and $|\eta|<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 noncollision 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\to c$ hadron decay chain. A working point with an average $b$-tagging efficiency of 70%, as determined using simulated $t\bar{t}$ 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 reconstructed 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_{\textrm{T}}>20$ GeV and $|\eta|<2.4$. The leading and second-leading jets (ordered in $p_{\textrm{T}}$) are then required to have $p_{\textrm{T1}}>160$ GeV and $p_{\textrm{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_{\textrm{T}}$, while for bbanti events the third jet is the third-leading jet in $p_{\textrm{T}}$. The final preselection requirement is that the minimum $\Delta 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 $b\bar{b}$ in parton showers and subprocesses such as $bg\to bg^{*}\to bb\bar{b}$.

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 $\phi$ boson decay, making it more difficult to accurately reconstruct $m_{\phi}$ from the two highest-$p_{\textrm{T}}$ jets. For example, the root-mean-square values of the reconstructed signal MC mass distributions for the two highest-$p_{\textrm{T}}$ jets for the 3-jet, 4-jet and 5-jet regions are 56 GeV, 65 GeV, and 74 GeV, respectively at $m_{\phi}$=600 GeV, and 163 GeV, 196 GeV, and 210 GeV, respectively at $m_{\phi}$=1200 GeV. Studies performed with and without jet multiplicity categorization demonstrate improvement in the expected limits from this categorization.

Signal sensitivity is enhanced with a transformation of the kinematic variables $p_{\textrm{T1}}$, $p_{\textrm{T2}}$, and the invariant mass of the two leading $b$-tagged jets, $m_{bb}$. Two-dimensional distributions of $p_{\textrm{T1}}$ versus $m_{bb}$ and of $p_{\textrm{T2}}$ versus $m_{bb}$ for events with the bbb classification are displayed in Figure 3. As $m_{\phi}$ increases the two high-$p_{\textrm{T}}$ jets from the $\phi$ 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_{\textrm{T}}$ jets smears out and it becomes more difficult to distinguish signal from background. However, since FSR occurs stochastically, the $\phi$ boson decays with little or no FSR in a subset of the signal events, and these have reconstructed masses close to the true $m_{\phi}$ (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_{\phi}$ resolution, a principal component analysis (PCA) [43] is performed on the three-dimensional distribution of the variables $m_{bb}$, $p_{\textrm{T1}}$, and $p_{\textrm{T2}}$ using events drawn from the signal MC sample with the bbb classification following preselection. Separate PCAs are performed for each of the seventeen simulated values of $m_{\phi}$ and for each of the three $n$-jet regions. Upon diagonalization of the covariance matrix for $m_{bb}$, $p_{\textrm{T1}}$, and $p_{\textrm{T2}}$, the first, second, and third principal components define the variables $m_{bb}^{\prime}$, $p_{\textrm{T1}}^{\prime}$, and $p_{\textrm{T2}}^{\prime}$, respectively. The point $(m_{bb}^{\prime},p_{\textrm{T1}}^{\prime},p_{\textrm{T2}}^{\prime})=(0,0,0)$ corresponds to the vector of mean values for $m_{bb}$, $p_{\textrm{T1}}$, and $p_{\textrm{T2}}$. Two-dimensional distributions of $p_{\textrm{T1}}^{\prime}$ versus $m_{bb}^{\prime}$ and of $p_{\textrm{T2}}^{\prime}$ versus $m_{bb}^{\prime}$ are shown in Figure 4.

The variables $m_{bb}$, $p_{\textrm{T1}}$, and $p_{\textrm{T2}}$ that are used to calculate the PCA variables are reconstructed signal MC quantities, and are subject to the experimental systematic uncertainties, such as jet energy scale uncertainty, discussed in Sec. 6. These uncertainties can lead to jet-multiplicity migration and other effects which alter the event populations in the 3-jet, 4-jet and 5-jet regions and change the PCA variables. Full error propagation is performed to account for these effects, from variation in $m_{bb}$, $p_{\textrm{T1}}$, and $p_{\textrm{T2}}$ to jet-multiplicity migration, PCA variable calculation, and, finally, experimental sensitivity to the Higgs boson cross section times branching ratio.

Examples of the relative contributions of $m_{bb}$, $p_{\textrm{T1}}$, and $p_{\textrm{T2}}$ to the rotated variable $m_{bb}^{\prime}$ are shown in Table 3 for the 3-jet region. The largest component of $m_{bb}^{\prime}$ comes from $m_{bb}$, regardless of the mass point. However, at larger values of $m_{\phi}$ - where there is greater FSR - the $p_{\textrm{T1}}$ and $p_{\textrm{T2}}$ components become more important.

The final event selection requirements are $p_{\textrm{T1}}^{\prime}>-10$ GeV and $p_{\textrm{T2}}^{\prime}>-50$ GeV, independent of $n$-jet and $m_{\phi}$. These requirements reduce the background while retaining a large fraction of the signal events in regions of high signal density in $(m_{bb}^{\prime},p_{\textrm{T1}}^{\prime},p_{\textrm{T2}}^{\prime})$ space, as shown in Figure 4. Using $m_{bb}^{\prime}$ instead of $m_{bb}$ as the final discriminant leads to increased sensitivity, which becomes more pronounced with increasing values of $m_{\phi}$.

5 Statistical analysis

The presence of a signal is tested with a binned maximum-likelihood fit to the data using $m^{\prime}_{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 $\sigma(pp\rightarrow b\bar{b}\phi)\times\mathcal{B}(\phi\rightarrow b\bar{b})$. 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 multijet MC sample, the backgrounds for both the bbb and bbanti regions are dominated by the subprocesses $gg\rightarrow gb\bar{b}$ and $gg\rightarrow ggb\bar{b}$ (such events enter the bbb region via gluon splitting into $b\bar{b}$ in the parton showering of one of the final-state gluons). However, the subprocesses $gb\rightarrow bb\bar{b}$ and $gg\rightarrow b\bar{b}b\bar{b}$ 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}^{\prime}$ distributions for both the $bbb$ and $bbanti$ classifications are plotted in Figure 5 for the multijet 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^{\prime}_{bb}$. Application of the $\chi^{2}$ probability test and the $F$-test [44] to the simulated multijet $m^{\prime}_{bb}$ distributions over all values of $m_{\phi}$ demonstrates that a first-order polynomial is sufficient to describe the ratio of the simulated multijet $bbb$ and $bbanti$ background shapes for signal masses $m_{\phi}<1200$ GeV, while a second-order polynomial is needed for $m_{\phi}\geq 1200$ GeV.

The data $bbb$ and $bbanti$ shapes, as well as their ratio, after applying the selection on $p_{\textrm{T1}}^{\prime}$ and $p_{\textrm{T2}}^{\prime}$ are qualitatively similar to the multijet MC distributions, as shown in Figure 6. Applying $\chi^{2}$ probability and $F$-test criteria to the data $m^{\prime}_{bb}$ distributions over all values of $m_{\phi}$, it is found that a first-order polynomial is sufficient for the 3-jet region with masses $m_{\phi}<1200$ GeV and for the 4-jet and 5-jet regions with masses $m_{\phi}<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 HistFactory [46] software tool. A product of Poisson probability terms over the bins of the $m^{\prime}_{bb}$ distributions involving the numbers of data events $n_{i,j}$ and the sum of expected signal and background yields $\nu_{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

[TABLE]

[TABLE]

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 $\mu$, which is common to all categories, is $\sigma(pp\rightarrow b\bar{b}\phi)\times\mathcal{B}(\phi\rightarrow b\bar{b})$ in picobarns.

Within the HistFactory framework, the second-order polynomial correction is implemented with the histograms $L_{j}$ and $P_{j}$, which are binned in $m^{\prime}_{bb}$ and have bin contents given by the bin center value and bin center value squared, respectively. The normalization parameters for these histograms, $A_{k}$ and $Q_{k}$ , correspond to the linear and parabolic parameters, respectively, for the jet-multiplicity region $k$. The signal strength parameter $\mu$, 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 $\beta_{i,j}$ represents the systematic variation in the bin content as a function of the nuisance parameters $\alpha_{p}$. The nuisance parameters $\alpha_{p}$ are constrained within the allowed systematic variations by the $f_{p}(a_{p}|\alpha_{p},\sigma_{p})$ terms, where $a_{p}$ are auxiliary measurements and $\sigma_{p}$ denotes the uncertainty in $\alpha_{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 Poisson-constrained nuisance parameter associated with the signal MC prediction for the number of events entering the bin and the total statistical uncertainty in that bin.

6 Systematic uncertainties

This analysis relies on the prediction of the shapes and normalizations of the discriminating variable $m^{\prime}_{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 multijet MC sample.

6.1 Systematic uncertainties of the background model

The background model is validated through the $\chi^{2}$ probability and $F$-test analyses described in Sec. 5. As a cross-check, the fit procedure is applied to a small multijet MC sample with an equivalent integrated luminosity of 6.8 fb*-1* for eight $m_{\phi}$ hypotheses. Events in MC samples are weighted to reflect data-based corrections for pileup, flavor tagging and trigger efficiency. In order to use the data fit procedure without modification, a special unweighted version of the MC sample is produced using acceptance-rejection sampling and the total weight of each event. 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 $\chi^{2}$ per degree of freedom is 1.09 for the eight measurements, indicating that no statistically significant spurious signal is found by the analysis.

6.2 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_{\textrm{T}}<300$ GeV. They are extrapolated to $p_{\textrm{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_{\textrm{T}}$ and $|\eta|$ [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[51]. Since the background estimation is data-driven, this scaling only affects signal. Scale factors obtained by comparing data and simulated dilepton $t\bar{t}$ events, which are enriched in $b$-jets, are used to correct simulation-to-data efficiency differences in the $b$-jet trigger for $p_{\textrm{T}}<240$ GeV. For $p_{\textrm{T}}>240$ GeV, due to the limited size of the $t\bar{t}$ data sample, extrapolation based on simulation is used, as described in detail in Ref. [51]. 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_{\textrm{T}}<240$ GeV, uncertainty in the simulation-based extrapolation to $p_{\textrm{T}}>240$ GeV, and uncertainties in the dependence of the $b$-jet trigger efficiency on jet $p_{\textrm{T}}$ and $\eta$. 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, 52]. 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. [53], and using the LUCID-2 detector for the baseline luminosity measurements [54].

6.3 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 Madgraph5_aMC@NLO [55, 56] with a 4 flavor scheme (4FS) PDF interfaced to the Pythia 8.205 parton showering model. The particle-level Higgs boson mass, Higgs boson $p_{\textrm{T}}$, and the $p_{\textrm{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]. Variations in PDF parameters can modify the $\phi$ boson detection efficiency through changes to the relative contribution of the four Higgs production Feynman diagrams of Figure 1, as well as the distributions of kinematic variables, such as $p_{\textrm{T}}$ and $\eta$ for the third $b$-jet. The PDF uncertainty is likely overestimated, but remains small compared to the statistical uncertainty.

7 Results

The search for a single heavy neutral Higgs boson $\phi$ produced in association with $b$-quarks and decaying into a $b\bar{b}$ pair shows no significant excess above the SM background for any of the analyzed mass points. The postfit bbb category distributions of the rotated $bb$ invariant mass $m^{\prime}_{bb}$ are shown in Figure 7 together with the $m^{\prime}_{bb}$ distribution extracted from the bbanti region (prefit background).

7.1 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 $\phi$ decaying into a $b\bar{b}$ pair are set. Figure 8 presents the observed and expected limits for $\sigma(pp\rightarrow b\bar{b}\phi)\times\mathcal{B}(\phi\rightarrow b\bar{b})$ at the $95\%$ confidence level (C.L.). The limits at specific mass points are calculated with the $\text{C.L.}_{\text{s}}$ method [57]. The expected and observed limits between mass points are given by straight lines connecting the mass point limits. Linear interpolations of the PCA parameters for adjacent mass points are used to calculate the observed $95\%$ C.L. limits for masses between the mass points. The straight lines in Figure 8 are consistent with these observed limits to within the accuracy of the method, which grows from 0.00 to 0.56 times the $+1\sigma$ expected limit band as the mass is varied from the mass point itself to a position midway between adjacent mass points.

The leading sources of uncertainty affecting the best-fit value of $\sigma\times\mathcal{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_{\phi}=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.

7.2 Model interpretations

The two 2HDM scenarios with enhanced $pp\rightarrow b\bar{b}\phi$ production and $\phi\rightarrow b\overline{b}$ decay at large ${\tan\beta}$ 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 [58] and for the $m_{h}^{\mathrm{mod}+}$ and $m_{h}^{\mathrm{mod}-}$ scenarios [59]. The Higgs boson production cross-sections and branching ratios are calculated using the procedures outlined in the LHC Higgs Cross-section Working Group report [60]. The cross-sections for Higgs boson production through $b\bar{b}$ fusion [14] are determined by matching the 5FS [61, 62] and 4FS [63, 64] cross-section calculations. For the hMSSM scenario, the Higgs boson masses and branching ratios are calculated using HDecay [65, 66]. For the $m_{h}^{\mathrm{mod}+}$ and $m_{h}^{\mathrm{mod}-}$ scenarios, the Higgs boson masses and couplings are calculated with FeynHiggs [24, 67, 68, 69, 70], and the branching ratios are calculated by combining the most precise results from FeynHiggs, HDecay, and PROPHECY4f [71, 72].

The 95% C.L. exclusion limits on ${\tan\beta}$ 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_{h}^{\mathrm{mod}+}$ and $m_{h}^{\mathrm{mod}-}$ scenarios, in which top squark mixing parameters are chosen to allow a wide range of ${\tan\beta}$ values while maintaining compatibility with $m_{h}=125$ GeV. The hMSSM limits obtained by this search are comparable to the limits obtained in the charged Higgs boson search in the $H^{+}\rightarrow\tau\nu$ decay channel by ATLAS and CMS Collaborations [73, 74] but not as stringent as the hMSSM limits obtained in the ATLAS and CMS searches for heavy neutral Higgs bosons decaying via $\phi\rightarrow\tau^{+}\tau^{-}$ [11, 12].

The 95% C.L. ${\tan\beta}$ exclusion limits for this search assuming the Type Y (flipped) 2HDM scenario are presented, first in Figure 10 for the specific $\phi$ mass of 450 GeV as a function of $\cos(\beta-\alpha)$ and then in Figure 11 as a function of $m_{\phi}$ in the alignment limit $\cos(\beta-\alpha)=0$. For these limits, it is assumed that the flipped 2HDM is *CP-*conserving with $m_{h}=125$ GeV and $m_{H}=m_{H^{\pm}}=m_{A}$. The model grid points are generated using SusHi [62] and 2HDMC [75]. These limits complement the flipped 2HDM limits obtained from the searches for $A\rightarrow Zh$ in ATLAS [76], which exclude regions with $|\cos(\beta-\alpha)|\gtrsim 0.2$ or $\tan\beta\lesssim 4$, and from the search for $A\rightarrow ZH$ in ATLAS [77], which excludes regions with $m_{A}-m_{H}\gtrsim 100$ GeV.

8 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% C.L. upper limits are in the range 0.6–4.0 pb. Compared to heavy neutral Higgs boson searches utilizing $\phi\to\tau^{+}\tau^{-}$ or $A\to Zh$ decays, these limits expand the excluded Type Y (flipped) 2HDM parameter space into regions with $|\cos(\beta-\alpha)|\approx 0$ and $\tan\beta\gtrsim 20$.

——————————————–

Acknowledgments

We thank CERN for the very successful operation of the LHC, as well as the support staff from our institutions without whom ATLAS could not be operated efficiently.

We acknowledge the support of ANPCyT, Argentina; YerPhI, Armenia; ARC, Australia; BMWFW and FWF, Austria; ANAS, Azerbaijan; SSTC, Belarus; CNPq and FAPESP, Brazil; NSERC, NRC and CFI, Canada; CERN; CONICYT, Chile; CAS, MOST and NSFC, China; COLCIENCIAS, Colombia; MSMT CR, MPO CR and VSC CR, Czech Republic; DNRF and DNSRC, Denmark; IN2P3-CNRS and CEA-DRF/IRFU, France; SRNSFG, Georgia; BMBF, HGF and MPG, Germany; GSRT, Greece; RGC and Hong Kong SAR, China; ISF and Benoziyo Center, Israel; INFN, Italy; MEXT and JSPS, Japan; CNRST, Morocco; NWO, Netherlands; RCN, Norway; MNiSW and NCN, Poland; FCT, Portugal; MNE/IFA, Romania; MES of Russia and NRC KI, Russia Federation; JINR; MESTD, Serbia; MSSR, Slovakia; ARRS and MIZŠ, Slovenia; DST/NRF, South Africa; MICINN, Spain; SRC and Wallenberg Foundation, Sweden; SERI, SNSF and Cantons of Bern and Geneva, Switzerland; MOST, Taiwan; TAEK, Turkey; STFC, United Kingdom; DOE and NSF, United States of America. In addition, individual groups and members have received support from BCKDF, CANARIE, Compute Canada and CRC, Canada; ERC, ERDF, Horizon 2020, Marie Skłodowska-Curie Actions and COST, European Union; Investissements d’Avenir Labex, Investissements d’Avenir Idex and ANR, France; DFG and AvH Foundation, Germany; Herakleitos, Thales and Aristeia programmes co-financed by EU-ESF and the Greek NSRF, Greece; BSF-NSF and GIF, Israel; La Caixa Banking Foundation, CERCA Programme Generalitat de Catalunya and PROMETEO and GenT Programmes Generalitat Valenciana, Spain; Göran Gustafssons Stiftelse, Sweden; The Royal Society and Leverhulme Trust, United Kingdom.

The crucial computing support from all WLCG partners is acknowledged gratefully, in particular from CERN, the ATLAS Tier-1 facilities at TRIUMF (Canada), NDGF (Denmark, Norway, Sweden), CC-IN2P3 (France), KIT/GridKA (Germany), INFN-CNAF (Italy), NL-T1 (Netherlands), PIC (Spain), ASGC (Taiwan), RAL (UK) and BNL (USA), the Tier-2 facilities worldwide and large non-WLCG resource providers. Major contributors of computing resources are listed in Ref. [78].