1 Introduction

The recent discovery of the Higgs boson [1, 2] marked a new era for fundamental physics. For the first time an electroweak-scale scalar resonance has been discovered, supposedly a remnant of the mechanism underlying electroweak symmetry breaking [3–6].

While elementary scalar particles have been observed in nature for the first time, they are often an integral part of Standard Model (SM) extensions, e.g. supersymmetry or general N-Higgs Doublet Models. When these extensions contain complex scalar fields, as a result, CP-odd scalars are introduced in the spectrum of the theory. Hence, since their existence would be evidence for physics beyond the SM, searches for CP-odd scalars are at the core of the current LHC program.

Recently, CP-odd scalars as mediators between Dark Matter (DM) and SM particles have received attention as possible explanations of the diffuse gamma-ray excess from the Galactic Centre [7–10] in the contexts of the so-called Coy Dark Matter models [11–14], and the Next-to-Minimal-Supersymmetric-Standard-Model (NMSSM) [15, 16] (for subsequent discussions in the NMSSM, see [17–21]). Hence, they are included as mediators in simplified models by the ATLAS and CMS collaborations to recast searches for jets and missing transverse energy (‘monojet’) during the upcoming LHC runs [22, 23].

In contrast to the widely accepted paradigm that new physics particles have to be heavy, i.e. with masses of \(\mathcal {O}\)(1) TeV or beyond, the mass of CP-odd scalars is almost unconstrained by direct searches. As interactions between gauge bosons and CP-odd scalars are only induced via higher-dimensional operators, e.g. \(\frac{1}{2} A \epsilon _{\mu \nu \sigma \rho } V_{\sigma \rho } V^{\mu \nu }\), limits from LEP are fairly weak. The main collider sensitivities may primarily arise from bottom quark or top quark-associated productions (for recent explorations, see [24, 25]). Further, due to the predicted velocity suppression in direct detection experiments for CP-odd scalar mediators, even light CP-odd scalars are still in agreement with experimental observations. Some constraints from flavour physics exist but limits are again weak if \(m_A \gtrsim 5\) GeV [26], assuming the CP-odd scalar interacts with fermions in agreement with the hypothesis of minimal flavour violation [27].

Therefore, indirect detection experiments and direct searches at the LHC appear to be the most sensitive ways to search for the existence of electroweak-scale CP-odd scalar particles. Most previous studies of direct searches for CP-odd scalars at the LHC have been focussed on the gluon-fusion production mechanism [28, 29]. Instead, in this paper we explore the direct production of such particles in association with a top quark pair and subsequent decay into a bottom quark pair, \(p p \rightarrow t \bar{t} A \rightarrow t \bar{t} b \bar{b}\). Thus, we derive limits on the mass and coupling strength of the CP-odd scalar in a process with unsuppressed fermion couplings only.Footnote 1

This paper is organised as follows. In Sect. 2 we briefly outline the way we incorporate the CP-odd scalar into the theory, using a simplified model approach. The event generation and details of the final state reconstruction are described in Sects. 3 and 4. In Sect. 5 we derive limits on the mass of the CP-odd scalar and its couplings to top quarks. Such limits can be applied to models where the CP-odd scalar arises as part of a Higgs multiplet. In Sect. 6 we recast these limits in the context of a two-Higgs-doublet model (2HDM) and the Next-to-Minimal Supersymmetric Standard Model (NMSSM). Finally, in Sect. 7 we offer conclusions.

2 Simplified model

While CP-odd scalars are present in many extensions of the SM, for simplicity and generality of our results, we use a simplified model approach [30] to parameterise the contribution of this particle to the process \(p p \rightarrow t \bar{t} A \rightarrow t \bar{t} b \bar{b}\). More precisely, we add couplings of the CP-odd scalar with the bottom and top quarks to the full SM Lagrangian

$$\begin{aligned} \mathcal {L}_{\mathrm {UV}} \supset \mathcal {L}_{\mathrm {SM}} + \mathcal {L}_{\mathrm {CP-odd}}, \end{aligned}$$
(1)

where

$$\begin{aligned} \mathcal {L}_{\mathrm {CP-odd}} = i \frac{g_{t} y_t}{\sqrt{2}}~\bar{t} \gamma _5 t A + i \frac{g_{b} y_b}{\sqrt{2}}~\bar{b} \gamma _5 b A, \end{aligned}$$
(2)

and \(g_i\) (\(i=t,b\)) parameterises the deviation from the SM Yukawa coupling \(y_i = m_i/v\).

Recently, a similar approach was proposed to recast monojet searches at the LHC in terms of scalar mediators between the SM and a secluded sector [31–34]. In a similar way, we will focus on the minimal set of free parameters relevant to the process considered. Throughout this paper we will assume A to be a narrow resonance with \(2 m_b \le m_A < 2m_t\). Hence, in our approach the CP-odd scalar decays exclusively into bottom quarks with \(\mathrm {BR}(A\rightarrow b\bar{b})=1\) and its width \(\Gamma _A\) is completely determined by the value of \(g_b\). For a narrow resonance, the kinematic distributions are expected to remain largely independent of the value of \(\Gamma _A\) and interference effects with other electroweak-scale particles, necessarily present in a UV-complete model, are expected to be small [35, 36].

The purpose of our simplified model approach is to explicitly include the degrees of freedom of a UV-complete model that contribute to the process of interest, i.e. the process \(pp \rightarrow t\bar{t} A \rightarrow t\bar{t}b\bar{b} \). Hence, the CP-odd scalar A and its couplings to top and bottom quarks are to be understood as residual low-energy states and interactions of well-defined UV models.Footnote 2 To minimise potential contributions of electroweak-scale resonances to our signal region we use a “bump hunt” inspired reconstruction approach, as discussed in Sects. 3–4.

3 Event generation and simulation details

3.1 Signal and background modelling

Signal and background samples corresponding to pp collisions at \(\sqrt{s}=14\) TeV are generated using the Madgraph5 2.1.1 [38] leading-order (LO) generator and the CTEQ6L1 [39] set of parton distribution functions (PDF), interfaced to Pythia v6.427 [40] for parton showering and fragmentation and using the Perugia2011C [41] underlying event tune. In all cases, a top quark mass of 172 GeV is assumed and top quarks are decayed inclusively by Pythia.

Table 1 Leading-order cross section for \(t\bar{t}A\) production in pp collisions at \(\sqrt{s}=14\) TeV as a function of the A boson mass \(m_A\). As discussed in the text, this LO cross section is obtained assuming \(g_t=1\) and will be multiplied by a k-factor of 1.3 to approximate the NLO cross section
Full size table
Fig. 1
figure 1

a Comparison of the leading-order cross section for \(t\bar{t}A\) (solid line) and \(t\bar{t}h\) (dashed line) in pp collisions at \(\sqrt{s}=14\) TeV as a function of Higgs boson mass. In both cases a value of \(g_t = 1\) is assumed. b Comparison of the Higgs boson \(p_\mathrm{{T}}\) between \(t\bar{t}A\) (black) and \(t\bar{t}h\) (red) for two different values of the Higgs boson mass, 20 GeV (solid) and 100 GeV (dashed)

Full size image

Samples of \(t\bar{t}A\) signal events are generated for different values of the A boson mass, \(m_A = 20, 30,40, 60, 80\) and 100 GeV, and assuming \(g_t=1\) and \(\mathrm {BR}(A \rightarrow b\bar{b})=1\). A model corresponding to the Lagrangian shown in Eq. (1) is implemented using Feynrules 2.1 [42] and is imported as an UFO model [43] in Madgraph5. The LO signal cross section predicted by Madgraph5 (see Table 1) is scaled by a k-factor of 1.3. This k-factor is obtained as the ratio of the NLO to LO cross sections for \(t\bar{t}h\) production, where h is a CP-even Higgs boson. It has been checked that this k-factor is rather constant as a function of \(m_h\), varied between 20 and 125 GeV. Figure 1a compares the production cross section between \(t\bar{t}h\) and \(t\bar{t}A\) as a function of the Higgs boson mass, in both cases assuming \(g_t=1\). The ratio between both cross sections varies significantly versus mass, with the \(t\bar{t}h\) cross section being about a factor of 20 larger than the \(t\bar{t}A\) cross section at a mass of 20 GeV, and only about a factor of 2 larger at a mass of 120 GeV [44]. This difference results from the presence of the extra \(\gamma _5\) factor in the interaction between a CP-odd Higgs boson and the top quark, compared to the case of a CP-even Higgs boson. Another consequence of the different interaction is that a CP-odd Higgs boson has a substantially harder \(p_\mathrm{{T}}\) spectrum compared to the CP-even case, particularly at low mass, as illustrated in Fig. 1b. This is a key feature exploited in this analysis, as discussed in Sect. 4.

A large sample of \(t\bar{t}\)+jets background events is generated including tree-level diagrams with up to two additional partons in the 5F scheme (i.e. including b- and c-quarks). To avoid double-counting of partonic configurations generated by both the matrix-element calculation and the parton shower, a parton–jet matching scheme (“MLM matching”) [45] is employed. The sample is normalised to a cross section of 990 pb obtained using Top \(++\) v2.0 [46] at next-to-next-to-leading order (NNLO) in QCD, including resummation of next-to-next-to-leading logarithmic (NNLL) soft gluon terms [47–51], and using the MSTW 2008 NNLO [52, 53] PDF set. The \(t\bar{t}\)+jets sample is generated inclusively, but events are categorised depending on the flavour content of additional particle jets in the event (i.e. jets not originating from the decay of the \(t\bar{t}\) system). Particle jets are reconstructed with the anti-\(k_t\) [54–56] algorithm with a radius parameter \(R=0.4\) and are required to have \(p_\mathrm{{T}}>15\) GeV and \(|\eta |<2.5\). Events where at least one such particle jet is matched within \(\Delta R<0.4\) to a b-hadron with \(p_\mathrm{{T}}>5\) GeV not originating from a top quark decay are generically labelled as \(t\bar{t}+{\ge }1b\) events. Similarly, events where at least one such particle jet is matched within \(\Delta R<0.4\) to a c-hadron with \(p_\mathrm{{T}}>5\) GeV not originating from a W boson decay, and that are not labelled already as \(t\bar{t}+{\ge }1b\), are labelled as \(t\bar{t}+{\ge }1c\) events. Events labelled as either \(t\bar{t}+{\ge }1b\) or \(t\bar{t}+\ge 1c\) are generically referred to below as \(t\bar{t}+ \)HF events, where HF stands for “heavy flavour”. We do not apply dedicated corrections to the normalisation of \(t\bar{t}+ \)HF events, since Run 1 searches at the LHC [57] showed that the LO prediction from Madgraph5 using the same settings as us is consistent with data within \({\sim } 20\,\%\), and a larger systematic uncertainty will be assumed in this study. As in Ref. [57], a finer categorisation of \(t\bar{t}+ \)HF events is considered for the purpose of assigning systematic uncertainties associated with the modelling of heavy-flavour production in different topologies. In this way, a distinction is made between events with only one extra heavy-flavour jet satisfying the above cuts (referred to as \(t\bar{t}+b\) or \(t\bar{t}+c\)), events with two extra heavy-flavour jets (referred to as \(t\bar{t}+b\bar{b}\) or \(t\bar{t}+c\bar{c}\)) and events with one extra heavy-flavour jet containing two b- or c-hadrons (referred to as \(t\bar{t}+B\) or \(t\bar{t}+C\)). The remaining events are labelled as \(t\bar{t}+\) light-jet events, including those with no additional jets.

Additional background samples corresponding to \(t\bar{t}W\), \(t\bar{t}Z\) and \(t\bar{t}h_\mathrm{SM}\) production, where \(h_\mathrm{SM}\) is the SM Higgs boson, are also produced. The \(t\bar{t}W\) sample is generated requiring at least one W boson in the event to decay leptonically, and is normalised to the corresponding LO cross section, 0.404 pb, times a k-factor of 1.4 [58]. The \(t\bar{t}Z\) sample is generated requiring \(Z \rightarrow q\bar{q}\) decays and is normalised to the corresponding LO cross section, 0.353 pb, times a k-factor of 1.3 [58]. Finally, the \(t\bar{t}h_\mathrm{SM}\) sample is generated assuming \(m_h=125\) GeV and requiring \(h \rightarrow b\bar{b}\) decays. It is normalised to the NLO cross section [59–61], 0.611 pb, times the \(h_\mathrm{SM} \rightarrow b\bar{b}\) branching ratio of 57.7 % [62–65], collected in Ref. [66]. In these samples \(Z \rightarrow q\bar{q}\) and \(h_\mathrm{SM} \rightarrow b\bar{b}\) decays are performed by Madgraph5 and top quarks and W bosons are decayed by Pythia.

3.2 Event reconstruction

The generated samples at the particle level are processed through a simplified simulation of the detector response and object reconstruction.

Isolated leptons (electrons or muons) are required to have \(p_\mathrm{{T}}>25\) GeV and \(|\eta |<2.5\). Furthermore, they are required to not overlap with jets, as discussed below. A typical per-lepton identification efficiency of 80 % is assumed.

Stable particles from Pythia, except for muons and neutrinos, are processed through a simplified simulation of a calorimeter. The four-momenta of particles falling within the same window in \(\eta -\phi \) space of size \(\Delta \eta \times \Delta \phi = 0.1\times 0.1\) are added together to simulate the finite granularity of calorimeter cells. For each cell, the total three-momentum is rescaled such as to make the cell massless. Cells with energy larger than 0.1 GeV and \(|\eta |<5.0\) become the inputs to the jet algorithm. Several types of jets are considered in this analysis.

The anti-\(k_t\) algorithm is used to reconstruct jets with two different radius parameters, \(R=0.2\) and \(R=0.4\), referred to as AKT2 and AKT4 jets, respectively. The minimum jet \(p_\mathrm{{T}}\) threshold for reconstruction is 5 GeV. During jet reconstruction, no distinction is made between identified electrons and jet energy deposits, and so every electron is also associated with a reconstructed jet. In order to remove this double counting, if any of the jets in the AKT2 and AKT4 collections lie within \(\Delta R=0.2\) of a selected electron, the closest jet from each jet collection is discarded. Since this analysis has a large number of b-quark initiated jets, for which a significant fraction of energy is carried away by muons in semi-muonic b-hadron decays, the four-momenta of all reconstructed muons with \(p_\mathrm{{T}}>4\) GeV that are ghost-associated [67, 68] to a jet are added to the calorimeter jet four-momentum. After this correction, a minimum \(p_\mathrm{{T}}\) requirement of 15 and 25 GeV is made for AKT2 and AKT4 jets, respectively. All jets are required to satisfy \(|\eta |<2.5\). Finally, any electron or muon within \(\Delta R=0.4\) of a selected AKT4 jet is discarded. In this analysis AKT4 jets are used to define the minimum jet multiplicity required in the event selection, while AKT2 jets are used to define the b-tag multiplicity of the event. The latter is particularly important since at low \(m_A\) values the b-quarks from the \(A \rightarrow b\bar{b}\) decay emerge with small angular separation. The flavour of an AKT2 jet is determined by matching it within \(\Delta R=0.15\) with a b-hadron or a c-hadron (not originating from a b-hadron decay), resulting in the jet being labelled as b-jet or c-jet, respectively. The rest of the jets are taken to originate from the fragmentation of a light quark or gluon and are labelled as “light jets”. Heavy-flavour tagging is modelled in a probabilistic fashion by assigning a per-jet efficiency of 70 % to b-jets, 20 % to c-jets, and 0.7 % to light jets.

In addition, jets are reconstructed with the Cambridge–Aachen (C/A) algorithm [69, 70] for the purpose of reconstructing the \(A \rightarrow b\bar{b}\) decay, taking advantage of the boost with which A bosons are produced in the \(t\bar{t}A\) process (see Fig. 1b). Two radius parameters are considered for C/A jets, \(R^\mathrm{C/A}=0.6\) and 0.8, referred to as CA6 and CA8 jets, respectively. The choice of radius for C/A jets is optimised in order to optimally reconstruct the \(t\bar{t}A\) signal depending on the value of \(m_A\). To minimise the impact of soft radiation and pile-up (the latter not modelled in this analysis), the mass-drop (a.k.a. BDRS) filtering algorithm [71, 72] is applied to the reconstructed C/A jets. For BDRS filtering, the following parameters are used: \(\mu _\mathrm{frac}=0.67\) and \(y_\mathrm{cut}=0.09\) [73]. A semi-muonic energy correction is also applied to the C/A jet four-momentum, as in the case of AKT2 and AKT4 jets.

4 Experimental analysis

4.1 Analysis strategy and event selection

This search is focussed on the \(t\bar{t}A \rightarrow W^+b W^- \bar{b}b\bar{b}\) process, with one of the W bosons decaying leptonically and the other W boson decaying hadronically. The resulting final state signature is thus characterised by one electron or muon, and high jet and b-jet multiplicities that can be exploited to suppress the background, dominated by \(t\bar{t}\)+jets production. Therefore, the following preselection requirements are made: exactly one electron or muon, \(\ge \)5 AKT4 jets and \(\ge \)3 AKT2 b-tagged jets, in the following simply referred to as \(\ge \)5 jets and \(\ge \)3 b-tags. In order to optimise the sensitivity of the search, the selected events are categorised into two separate channels depending on the number of b-tags (3 and \(\ge \)4). The channel with \(\ge \)5 jets and \(\ge \)4 b-tags has the largest signal-to-background ratio and therefore drives the sensitivity of the search. It is dominated by \(t\bar{t}\)+HF background. The channel with 3 b-tags has significantly lower signal-to-background ratio and the background is enriched in \(t\bar{t}+ \)light jets. The simultaneous analysis of both channels is useful to calibrate in-situ the \(t\bar{t}+\) jets background prediction (including its heavy-flavour content) and constrain the related systematic uncertainties, as will be discussed in Sect. 4.3. This is a common strategy used in many experimental searches in the ATLAS and CMS collaborations [57, 74, 75], which we mimic here in order to obtain more realistic projected sensitivities.

Fig. 2
figure 2

Distribution of \(\Delta R\) between the two b-quarks from the \(A \rightarrow b\bar{b}\) decay prior to any selection requirements, for different values of \(m_A\): a \(m_A=20, 30\) and 40 GeV, and b \(m_A=60, 80\) and 100 GeV

Full size image

An extra handle is provided by the significant boost of the A boson in a fraction of signal events, which results in the two b-jets from the \(A \rightarrow b\bar{b}\) decay emerging with small angular separation between them. This is particularly relevant for low \(m_A\) values, as shown in Fig. 2. As a result, the A boson decay products can be reconstructed into a single fat jet, whose mass distribution would show a resonant structure peaked at the correct \(m_A\) value. This feature is also very powerful to discriminate against the background. Therefore, a further requirement is made to have at least one C/A BDRS-filtered jet with radius parameter \(R^\mathrm{CA}\) and minimum \(p_\mathrm{{T}}\) depending on the \(m_A\) hypothesis being tested. In order to correctly reconstruct a significant fraction of the signal while rejecting as much background as possible, CA6 jets are used for \(m_A \le 40\) GeV, while CA8 jets are used for higher \(m_A\) values (up to 100 GeV). The minimum \(p_\mathrm{{T}}\) requirements on the C/A jets are 60, 100, 120, 150, 200 and 250 GeV for \(m_A=20, 30, 40, 60, 80\) and 100 GeV, respectively. As shown in Fig. 2, for high values of \(m_A\) only a small fraction of signal events would have the A decay products contained within the CA8 jet. The small signal acceptance comes with the benefit of improved background rejection and the ability to reconstruct the A boson mass, desirable in such simple analysis. However, it is expected that a dedicated multivariate analysis focussed on the sample rejected by this analysis, similar in spirit to the ATLAS and CMS searches for the SM Higgs boson in \(t\bar{t} h\), \(h \rightarrow b\bar{b}\) [57, 74], could also achieve significant signal sensitivity at high \(m_A\). Evaluating this possibility is beyond the scope of this study. The number of b-tags inside the C/A jet is determined by matching the b-tagged AKT2 jets within a cone of radius \(\Delta R =0.75 R^\mathrm{C/A}\). Finally, a requirement is made that the C/A jets have \(\ge \)2 b-tags inside. In the case of more than one selected C/A jet, the leading \(p_\mathrm{{T}}\) one is chosen.

Table 2 presents the expected yields for signal and the SM backgrounds per fb\(^{-1}\) of integrated luminosity as a function of the selection cuts applied in each of the analysis channels under consideration: (\(\ge \)5j, 3b) and (\(\ge \)5j, \(\ge \)4b). In the case of the (\(\ge \)5j, 3b) channel, the dominant background after final selection is \(t\bar{t}\)+light jets, where typically the two b-quarks from the top quark decays, as well as the c-quark from the \(W \rightarrow c\bar{s}\) decay, are b-tagged. In contrast, in the (\(\ge \)5j, \(\ge \)4b) channel half of the background is \(t\bar{t}+{\ge }1b\), with \(t\bar{t}\)+\(b\bar{b}\) being its leading contribution. The rest of the background is approximately equally split between \(t\bar{t}+{\ge }1c\) and \(t\bar{t}+\) light jets. In this table the expected contribution from \(t\bar{t}A\) signal is obtained under the assumptions that \(g_t=2\) and \(\mathrm {BR}(A\rightarrow b\bar{b})=1\). Both analysis channels have approximately the same amount of signal, while the background is about a factor of 4 higher in the (\(\ge \)5j, 3b) channel than in the (\(\ge \)5j, \(\ge \)4b) channel. Together with the different composition of the background, the very different signal-to-background ratio between both channels is the primary motivation for analysing them separately.

The final discriminating variable is the invariant mass of the selected C/A jet, referred to as “BDRS jet mass”. Figures 3 and 4 show the expected distribution of the BDRS jet mass for signal and background in each of the analysis channels, for the different \(m_A\) values considered. The distributions correspond to \(\sqrt{s}=14\) TeV and are normalised to an integrated luminosity of 30 fb\(^{-1}\). For the assumed values of \(g_t=2\) and \(\mathrm {BR}(A\rightarrow b\bar{b})=1\), the signal is clearly visible on top of the background.

Table 2 Expected signal and SM backgrounds at \(\sqrt{s}=14\) TeV per fb\(^{-1}\) of integrated luminosity as a function of the selection cuts applied in each of the analysis channels under consideration (see text for details). The signal prediction is obtained under the assumptions that \(g_t=2\) and \(\mathrm {BR}(A\rightarrow b\bar{b})=1\). Several background categories have been merged for readability. The sum of \(t\bar{t}W\), \(t\bar{t}Z\) and \(t\bar{t}h_\mathrm{SM}\) is denoted as \(t\bar{t}+X\). The yields shown correspond to the optimised selections for two different values of \(m_A\), 30 and 80 GeV. Shown in bold are the signal and backgrounds expectations after full selection in each of the analysis channels considered
Full size table
Fig. 3
figure 3

Distribution of the BDRS jet mass in the two analysis channels considered after final selection: (top) (\(\ge \)5j, 3b) and (bottom) (\(\ge \)5j, \(\ge \)4b), for different values of \(m_A\) (20, 30 and 40 GeV). The prediction corresponds to \(\sqrt{s}=14\) TeV and an integrated luminosity of 30 fb\(^{-1}\). Several background categories have been merged for visibility. The expected contribution from the \(t\bar{t}A\) signal under the assumptions that \(g_t=2\) and \(\mathrm {BR}(A\rightarrow b\bar{b})=1\) is also shown (red histogram), stacked on top of the SM background. The dashed red line shows the \(t\bar{t}A\) signal distribution normalised to the background yield to better compare the shape to that of the background. The bottom panel displays the expected total systematic uncertainty on the total prediction prior to the fit to the pseudo-data

Full size image
Fig. 4
figure 4

Distribution of the BDRS jet mass in the two analysis channels considered after final selection: (top) (\(\ge \)5j, 3b) and (bottom) (\(\ge \)5j, \(\ge \)4b), for different values of \(m_A\) (60, 80 and 100 GeV). The prediction corresponds to \(\sqrt{s}=14\) TeV and an integrated luminosity of 30 fb\(^{-1}\). Several background categories have been merged for visibility. The expected contribution from the \(t\bar{t}A\) signal under the assumptions that \(g_t=2\) and \(\mathrm {BR}(A\rightarrow b\bar{b})=1\) is also shown (red histogram), stacked on top of the SM background. The dashed red line shows the \(t\bar{t}A\) signal distribution normalised to the background yield to better compare the shape to that of the background. The bottom panel displays the expected total systematic uncertainty on the total prediction prior to the fit to the pseudo-data

Full size image

4.2 Systematic uncertainties

Several sources of systematic uncertainty are considered that can affect the normalisation of signal and background and/or the shape of the BDRS jet mass distribution. Individual sources of systematic uncertainty are considered uncorrelated. For each systematic uncertainty, correlations are maintained across processes and analysis channels. The choices of what uncertainties to consider and their magnitude are inspired by recent \(t\bar{t}h_\mathrm{SM}\), \(h_\mathrm{SM}\rightarrow b\bar{b}\) searches at the LHC [57].

A 15 % normalisation uncertainty is assigned to the \(t\bar{t}\)+light-jet background corresponding to the modelling of the jet multiplicity spectrum. A 30 % normalisation uncertainty is assigned to each of the \(t\bar{t}+ \)HF background components (\(t\bar{t}+b\), \(t\bar{t}+b\bar{b}\), \(t\bar{t}+B\), \(t\bar{t}+c\), \(t\bar{t}+c\bar{c}\), \(t\bar{t}+C\)), and taken to be uncorrelated among them. These uncertainties are expected to be conservative given the recent progress in NLO predictions for \(t\bar{t}\) production with up to two jets merged with a parton shower [76], as well as NLO predictions for \(t\bar{t}+{\ge } 1b\) production in the 4F scheme matched to a parton shower [77]. Cross section uncertainties for \(t\bar{t}W\), \(t\bar{t}Z\) and \(t\bar{t}h_\mathrm{SM}\) backgrounds are taken to be 30 % for each process. Uncertainties associated with jet energy and jet mass calibrations are taken to be 5 % per jet, fully correlated between energy and mass and across all jets in the event. Finally, uncertainties on the b-, c- and light-jet tagging efficiencies are taken to be 3, 6 and 15 %, respectively. These uncertainties are taken as uncorrelated between b-jets, c-jets and light jets. As shown in Figs. 3 and 4, the resulting total background normalisation uncertainty is about 20 %, although the different uncertainty components have different shape in the final distribution.

4.3 Statistical method

The BDRS jet mass distribution in the two analysis channels under consideration (see Figs. 3, 4) are tested for the presence of a signal. To obtain the most realistic possible sensitivity projection, a sophisticated statistical analysis is performed, following very closely the strategy adopted in the experimental searches at the LHC.

For each \(m_A\) hypothesis, 95 % CL upper limits on the \(t\bar{t}A\) production cross section times branching ratio, \(\sigma (t\bar{t}A) \times \mathrm {BR}(A \rightarrow b\bar{b})\), are obtained with the CL\(_\mathrm{{s}}\) method [78, 79] using a profile likelihood ratio as test statistic implemented in the RooFit package [80, 81]. The likelihood function \(L(\mu ,\theta )\) depends on the signal-strength parameter \(\mu \), a multiplicative factor to the theoretical signal production cross section, and \(\theta \), a set of nuisance parameters that encode the effect of systematic uncertainties in the analysis. The likelihood function is constructed as a product of Poisson probability terms over all bins of the distributions analysed, and of Gaussian or log-normal probability terms, each corresponding to a nuisance parameter. For a given assumed value of \(\mu \), the profile likelihood ratio \(q_\mu \) is defined as:

$$\begin{aligned} q_\mu = -2\ln (L(\mu ,\hat{\hat{\theta }}_\mu )/L(\hat{\mu },\hat{\theta })), \end{aligned}$$
(3)

where \(\hat{\hat{\theta }}_\mu \) are the values of the nuisance parameters that maximise the likelihood function for a given value of \(\mu \), and \(\hat{\mu }\) and \(\hat{\theta }\) are the values of the parameters that maximise the likelihood function (with the constraint \(0\le \hat{\mu } \le \mu \)). The maximisation of the likelihood function over the nuisance parameters allows variations of the expectations for signal and background in order to improve the agreement with (pseudo-)data, yielding a background prediction with reduced overall uncertainty and thus resulting in an improved sensitivity. For a given \(m_A\) hypothesis, values of the production cross section (parameterised by \(\mu \)) yielding CL\(_\mathrm{{s}}<\)0.05, where CL\(_\mathrm{{s}}\) is computed using the asymptotic approximation [82], are excluded at \(\ge \)95 % CL.

5 Estimated limits on a light CP-odd scalar

Following the analysis steps and the statistical method outlined in Sects. 2–4, we estimate expected 95 % CL upper limits on the production cross section times branching ratio, \(\sigma (t\bar{t}A) \times \mathrm {BR}(A\rightarrow b\bar{b})\), as a function of \(m_A\) (see Fig. 5). Table 3 summarises the 95 % CL upper limits on \(\sigma (t\bar{t}A) \times \mathrm {BR}(A\rightarrow b\bar{b})\) as a function of \(m_A\) for different values of the integrated luminosity. Under the assumption that \(\mathrm {BR}(A\rightarrow b\bar{b})=1\), the upper limits on \(\sigma (t\bar{t}A) \times \mathrm {BR}(A\rightarrow b\bar{b})\) can be translated into upper limits on \(|g_t|\), which are summarised in Table 4.

Using the reconstruction strategy outlined in Sect. 4.1, a CP-odd scalar that couples with \(g_t=1\) can be excluded for \(20 \le m_A \le 90\) GeV with only \(30~\mathrm {fb}^{-1}\) of data (see Fig. 5). With an increased statistics of \(300~\mathrm {fb}^{-1}\) couplings as low as \(g_t \simeq 0.5\) can be constrained over a large mass range, i.e. \(30 \le m_A \le 80\) GeV.

Fig. 5
figure 5

Expected 95 % CL upper limits on \(\sigma (t\bar{t}A) \times \mathrm {BR}(A\rightarrow b\bar{b})\) as a function of \(m_A\) in pp collisions at \(\sqrt{s}=14\) TeV for an integrated luminosity of a 30 fb\(^{-1}\) and b 300 fb\(^{-1}\). The green and yellow bands correspond to 1 and 2 standard deviations, respectively, around the median expected limit, obtained under the background-only hypothesis. Also shown are the theoretical cross sections for \(\sigma (t\bar{t}A)\) for different assumed values of \(g_t\) (0.5, 1.0 and 2.0) and \(\mathrm {BR}(A\rightarrow b\bar{b})=1\)

Full size image

6 Interpretation of limits

A light CP-odd Higgs boson (\(m_A < 125\) GeV), which may or may not be related to global symmetries being present, exists in many extensions of the SM. Its couplings with gauge bosons are generically suppressed, yielding weak bounds from LEP. If \(m_A<m_{h_\mathrm{SM}}/2\), it may be searched via the decay \(h_\mathrm{SM} \rightarrow AA\). Although such decay sometimes has a large branching ratio, being in conflict with current Higgs precision data, there do exist scenarios, in both supersymmetric and non-supersymmetric theories, where the \(\mathrm {BR}(h_\mathrm{SM}\rightarrow AA)\) is suppressed. Therefore, new strategies for collider searches that could cover as large as possible model parameter space with a light CP-odd Higgs boson, are necessary. Next, we will interpret our collider analysis of \(t\bar{t} A\) in several representative beyond-SM scenarios.

6.1 2HDM

In the Minimal Supersymmetric Standard Model (MSSM), a supersymmetric extension of a type-II 2HDM, a scenario with a light CP-odd Higgs boson is hard to achieve, given constraints from precision Higgs data. This is not surprising since there are only two free parameters at tree level in the Higgs sector, due to supersymmetric interrelations. The picture, however, is changed in the 2HDM without supersymmetry. With a softly broken \(Z_2\) symmetry (\(\Phi _1 \rightarrow \Phi _1\), \(\Phi _2 \rightarrow -\Phi _2\)), which is often introduced to suppress scalar-mediated flavour-changing processes, the Higgs potential of the 2HDM is given by

$$\begin{aligned} V(\Phi _1,\Phi _2)= & {} m^2_1 \Phi ^{\dagger }_1\Phi _1+m^2_2 \Phi ^{\dagger }_2\Phi _2 + (m^2_{12} \Phi ^{\dagger }_1\Phi _2+{\mathrm {h.c.}}) \nonumber \\&+\frac{1}{2} \lambda _1 (\Phi ^{\dagger }_1\Phi _1)^2 +\frac{1}{2} \lambda _2 (\Phi ^{\dagger }_2\Phi _2)^2\nonumber \\&+\lambda _3 (\Phi ^{\dagger }_1\Phi _1)(\Phi ^{\dagger }_2\Phi _2) + \lambda _4 (\Phi ^{\dagger }_1\Phi _2)(\Phi ^{\dagger }_2\Phi _1)\nonumber \\&+ \frac{1}{2}\lambda _5[(\Phi ^{\dagger }_1\Phi _2)^2+{\mathrm {h.c.}}], \end{aligned}$$
(4)

where \(\Phi _{1,2}\) are complex \(SU(2)_L\) doublets. Assuming no CP-violation, the model has two CP-even and one CP-odd spin-0 neutral eigenstates, denoted as h, H, and A, respectively. Such a setup contains seven free parameters at tree level (including all Higgs masses), yielding a large parameter space that can accommodate a light CP-odd Higgs boson.

Table 3 Expected 95 % CL upper limits on \(\sigma (t\bar{t}A) \times \mathrm {BR}(A\rightarrow b\bar{b})\) as a function of \(m_A\) in pp collisions at \(\sqrt{s}=14\) TeV for different integrated luminosities
Full size table
Table 4 Expected 95 % CL upper limits on \(|g_t|\) as a function of \(m_A\) in pp collisions at \(\sqrt{s}=14\) TeV for different integrated luminosities, under the assumption that \(\mathrm {BR}(A\rightarrow b\bar{b})=1\)
Full size table

Theoretically, the SM-like Higgs boson \(h_\mathrm{SM}\), with mass of 125 GeV, could be either the light CP-even Higgs boson (h) or the heavy one (H). If \(m_A < m_{h_\mathrm{SM}}/2\), the decay \(h_\mathrm{SM}\rightarrow AA\) is kinematically allowed. Often the partial width for \(h_\mathrm{SM}\rightarrow AA\) becomes comparable or ever dominant over that of \(h_\mathrm{SM}\rightarrow b\bar{b}\), given that the latter is suppressed by the small value of the b-quark mass. Therefore, \(h_\mathrm{SM} \rightarrow AA\) decays become a good probe for these light bosonic particles. However, as discussed recently [83],Footnote 3 in the alignment limit [\(\cos (\beta -\alpha )=0\) if \(h_\mathrm{SM}=h\), and \(\sin (\beta -\alpha )=0\) if \(h_\mathrm{SM}=H\)], which is favoured by current precision Higgs measurements, the Higgs coupling \(g_{h_\mathrm{SM} AA}\) is reduced to:

$$\begin{aligned} \left| g_{h_\mathrm{SM} AA}\right| = \left| -\frac{2 m_A^2 + m_{h_\mathrm{SM}}^2 - 4 m_{12}^2/\sin 2\beta }{v}\right| . \end{aligned}$$
(5)

In the case that \(2 m_A^2 + m_{h_\mathrm{SM}}^2 \sim 4 m_{12}^2/\sin 2\beta \), the decay \(h_\mathrm{SM}\rightarrow AA\) would be greatly suppressed. Therefore, collider strategies are needed to probe these scenarios with \(m_A < m_{h_\mathrm{SM}}/2\), as well as the scenarios with \(m_A > m_{h_\mathrm{SM}}/2\).

We should note that the perturbation requirement for Higgs couplings yields bounds on \(\tan \beta \). Particularly, the coupling \(\lambda _1\) is related to the Higgs boson mass via the relation [83]:

$$\begin{aligned} \lambda _1 = \frac{m_h^2 + m_H^2 \tan ^2\beta - m_{12}^2 (\tan \beta +\tan ^3 \beta )}{v^2}. \end{aligned}$$
(6)

Assuming \(g_{h_\mathrm{SM}\rightarrow AA}=0\), it becomes

$$\begin{aligned} \lambda _1 = \frac{m_h^2 + \tan ^2\beta (m_H^2 - m_{h_\mathrm{SM}}^2/2 - m_A^2) }{v^2}. \end{aligned}$$
(7)

Given that \(m_H^2 - m_{h_\mathrm{SM}}^2/2 - m_A^2 > 0\) for \(m_A < m_{h_\mathrm{SM}}/2\), the perturbativity condition \(\lambda _1 < 4\pi \) immediately sets an upper bound on \(\tan \beta \) in this region:

$$\begin{aligned} \tan \beta< & {} \sqrt{\frac{4 \pi v^2 - m_h^2}{m_H^2 - m_{h_\mathrm{SM}}^2/2 - m_A^2 }} \nonumber \\< & {} \sqrt{\frac{4 \pi v^2 }{m_{h_\mathrm{SM}}^2/2 - m_A^2 }} \sim 10-20. \end{aligned}$$
(8)

These features are illustrated in Fig. 6. Additionally, the perturbation requirement for top Yukawa couplings can bound the \(\tan \beta \) value from below. So we will limit our discussions for \(\tan \beta > 0.1\).

Fig. 6
figure 6

Parameter region with \(\mathrm {BR}(h_\mathrm{SM}\rightarrow AA)<30\, \%\) in the 2HDM (blank region). The blank belts in the region with \(m_A < m_{h_\mathrm{SM}}/2\), which are characterised by different boundary colours, are yielded by different \(m_{12}^2\) values. The black dashed line represents a universal upper limit on \(\tan \beta \) due to the perturbation requirement for \(\lambda _1\)

Full size image
Fig. 7
figure 7

Sensitivity reach (at 95 % CL) of the \(t\bar{t}A\) and \(b\bar{b} A\) channels within a the type-I 2HDM and b the type-II 2HDM. The green bands represent a region where the recently observed gamma-ray excess from the Galactic Centre can be explained, yielding a DM annihilation cross section of \(\langle \sigma v \rangle \simeq 1-2.5 \times 10^{-26}\,\mathrm{cm}^3\,\mathrm{s}^{-1}\). Here the DM particle mass and the coupling between the mediator and the DM particles are assumed to be \(m_{\chi }=50\) GeV and \(y_\chi = 0.3\), respectively

Full size image

The expected sensitivities for probing these scenarios in the 2HDM via \(b\bar{b}A\) and \(t\bar{t}A\) production are presented in Fig. 7. The \(b\bar{b}A\) reach is estimated based on the projections from Ref. [85], neglecting systematic uncertainties. For illustration, we focus on type-I and type-II 2HDMs. Within a type-II 2HDM, the \(t\bar{t}A\) and \(b\bar{b}A\) channels are complementary to each other in searching for light CP-odd Higgs bosons, since the coupling \(g_{bbA}\) is \(\tan \beta \)-enhanced whereas \(g_{ttA}\) is \(\cot \beta \)-enhanced. With integrated luminosities in excess of 300 fb\(^{-1}\), the whole parameter region can be covered except a corner with relatively large \(m_A\) and moderate \(\tan \beta \). This is interesting given that low \(\tan \beta \) is particularly favoured by perturbativity. In contrast, within a type-I 2HDM, the coupling \(g_{bbA}\) would also be \(\cot \beta \)-enhanced, so both search channels are no longer probing complementary \(\tan \beta \) regions. As a matter of fact, in such scenario the \(t\bar{t}A\) channel provides better sensitivity to search for the light CP-odd Higgs boson over the whole mass range of \(20~\mathrm{GeV} < m_A < 100~\mathrm{GeV}\), although the high-\(\tan \beta \) region remains difficult to probe.

Searches for \(t\bar{t}A\) and \(b\bar{b} A\) also provide a probe for DM physics. For example, consider a Dirac fermion \(\chi \) that is a DM candidate, with mass \(m_\chi \), and coupling to the CP-odd scalar A via:

$$\begin{aligned} \mathcal {L} \supset y_\chi A \bar{\chi }i \gamma ^5 \chi . \end{aligned}$$
(9)

Integrating out A yields a dimension-six effective operator:

$$\begin{aligned} \mathcal {L}_\mathrm{eff} \sim \frac{- y_b y_\chi m_b }{\Lambda ^3} \bar{\chi }\gamma ^5 \chi \bar{b} \gamma ^5 b. \end{aligned}$$
(10)

Such an operator implies s-wave DM annihilation \(\chi \chi \rightarrow b\bar{b}\) with

$$\begin{aligned} \langle \sigma v \rangle = \frac{3}{8 \pi } \frac{y^2_\chi g_b^2 y_b^2 m_{\chi }^2}{(m_A^2 - 4 m_{\chi }^2)^2 + m_A^2 \Gamma _A^2} \sqrt{1-\frac{m_b^2}{m_\chi ^2}}, \end{aligned}$$
(11)

allowing an explanation for the recently observed diffuse gamma-ray excess from the Galactic Centre [7, 86]. In Fig. 7, the \(\tan \beta -m_A\) values consistent with an explanation of the gamma-ray excess are indicated, yielding a DM annihilation cross section of \(\langle \sigma v \rangle \simeq 1-2.5 \times 10^{-26} \mathrm{cm}^3 \mathrm{s}^{-1}\), under the assumptions that \(m_{\chi }=50\) GeV [87] and \(y_\chi = 0.3\). This scenario results in a spin-dependent and p-wave-suppressed direct detection signal, resulting in a weak bound from current direct detection searches. Monojet searches at the LHC would also be insensitive since the decay \(A \rightarrow \chi \chi \) would be kinematically forbidden, while the \(t\bar{t}A\), \(A \rightarrow b\bar{b}\) search would provide an effective probe.

6.2 NMSSM

Another class of benchmark scenarios for light CP-odd Higgs bosons arise in the NMSSM, with the superpotential and soft supersymmetry-breaking terms of its Higgs sector given by

$$\begin{aligned} \mathbf {W}&= \lambda \mathbf {S} \mathbf {H_u} \mathbf {H_d} + \frac{1}{3}\kappa \mathbf{S}^3, \nonumber \\ V_{\text {soft}}&= {m^2_{H_d}} |H_d|^2 + {m^2_{H_u}} |H_u|^2 + {m^2_S}|S|^2 \nonumber \\&\quad - (\lambda A_{\lambda } H_u H_d S + \text { h.c.})+\left( \frac{1}{3} \kappa A_{\kappa } S^3 + \text { h.c.} \right) , \end{aligned}$$
(12)

where \(H_d\), \(H_u\) and S denote the neutral Higgs fields of the \(\mathbf{H_d}\), \(\mathbf{H_u}\) and \(\mathbf{S}\) supermultiplets, respectively. For convenience, let us define its CP-even and CP-odd mass eigenstates as \(H_i\), \(i=1,2,3\), and \(A_j\), \(j=1,2\), respectively.

In contrast with the 2HDM case, the light CP-odd Higgs boson in the NMSSM often results from breaking an approximate global symmetry spontaneously, serving as an axion or a pseudo-Goldstone boson. Its appearance is thus less “artificial” than in a 2HDM (including the MSSM). Let us start with the tree-level mass matrix of the CP-odd Higgs bosons in the NMSSM:

$$\begin{aligned} \mathcal {M}^2_{P}= & {} \left( \begin{array}{ll} m_A^2 &{}\quad \lambda v \left( \frac{m_A^2 }{2\mu }\sin 2 \beta -\frac{3 \kappa \mu }{\lambda }\right) \\ &{}\lambda ^2 v^2 s_{2\beta } \left( \frac{m_A^2}{4\mu ^2} \sin 2 \beta +\frac{3 \kappa }{2\lambda }\right) -\frac{3 \kappa A_{\kappa } \mu }{\lambda }\\ \end{array} \right) ,\nonumber \\&m_A^2 = \frac{2 \mu (A_\lambda + \kappa s) }{\sin 2 \beta } \end{aligned}$$
(13)

which yields a determinant

$$\begin{aligned} \det (\mathcal {M}^2_{P}) = 9 \kappa \lambda v^2 \mu A_\lambda - \frac{6 A_\kappa \kappa \mu ^2}{\lambda \sin 2\beta } \left( A_\lambda + \frac{\kappa \mu }{\lambda } \right) . \end{aligned}$$
(14)

Necessarily, the scenarios with a light \(A_1\) (\(A_1\) denotes the lightest CP-odd Higgs boson) or \(m_{A_1} \rightarrow 0\) yield \(\det (\mathcal {M}^2_{P}) \rightarrow 0\) and vice versa, if such a stable vacuum exists. Among various possibilities, two have been studied extensively: R-symmetry (or R-limit) and Peccei–Quinn (PQ)-symmetry (or PQ-limit), both of which yield a vanishing determinant at tree level.

Another difference with a 2HDM is that the light CP-odd Higgs boson in the NMSSM is typically singlet-like. This can be understood since the Goldstone boson of a spontaneously broken global U(1)-symmetry is manifested as

$$\begin{aligned} A_1 \sim \sum _i \frac{q_i v_i}{v_{U(1)}} \Phi _i, \quad \Phi _i = S, H_u, H_d. \end{aligned}$$
(15)

Here \(v_{U(1)} = \sqrt{\sum q_i^2 v_i^2}\) is the U(1) breaking scale and \(q_i\) is the U(1) charge of \(\Phi _i\). An effective parameter \(\mu = \lambda \langle v_S\rangle \) of the electroweak scale with \(\lambda \sim \mathcal O(0.1)\) naturally yields \(v_S \gg v_u, v_d\), and hence a singlet-like pseudo-Goldstone boson. This feature renders such a light boson much more difficult to probe at colliders, compared to the 2HDM case.

Fig. 8
figure 8

Sensitivity reach (at 95 % CL) of the a \(t\bar{t}A_1\) and b \(b\bar{b} A_1\) channels within the R-limit scenario in the NMSSM. The scan is performed over all parameters, in the ranges \(0.1 \le \lambda \le 0.6\), \(0.1 \le \kappa \le 0.6\), \(-7 \le A_\lambda \le 7\) GeV, \(-7 \le A_\kappa \le 0\) GeV, \(0.1 \le \tan \beta \le 5\), and \(100 \le \mu \le 500\) GeV. We have assumed soft squark masses of 2 TeV, slepton masses of 200 GeV, \(A_{u,d,e} = -3.5\) TeV, and bino, wino and gluino masses of 100, 200, and 2000 GeV, respectively. The hue of the scatter points represents the corresponding \(\tan \beta \) values

Full size image
Fig. 9
figure 9

Sensitivity reach (at 95 % CL) of the a \(t\bar{t}A_1\) and b \(b\bar{b} A_1\) channels within the PQ-limit scenario in the NMSSM. The scan is performed over all parameters, in the ranges \(0.06 \le \lambda \le 0.6\), \(5 \le \kappa / \lambda \le 100\), \(|\varepsilon '| = \left| A_\lambda /\mu \tan \beta -1\right| \le 0.25\), \(-100 \le A_\kappa \le 0\) GeV, \(0.1 \le \tan \beta \le 5\), and \(100 \le \mu \le 500\) GeV. We have assumed soft squark masses of 2 TeV, slepton masses of 200 GeV, \(A_{u,d,e} = 3.5\) TeV, and bino, wino and gluino masses of 100, 200, and 2000 GeV, respectively. The hue of the scatter points represents the corresponding \(\tan \beta \) values

Full size image

Next, we evaluate the collider constraints on the two NMSSM scenarios discussed above:

  1. 1.

    R-limit: \(A_\lambda \rightarrow 0\), \(A_\kappa \rightarrow 0\), where the theory is approximately invariant under the transformation

    $$\begin{aligned}&H_u \rightarrow H_u \exp (i \phi _\mathrm{R}), \quad H_d \rightarrow H_d \exp (i \phi _\mathrm{R}),\nonumber \\&S \rightarrow S \exp (i \phi _\mathrm{R}), \end{aligned}$$
    (16)

    and the tree-level couplings of the R-axion \(A_1\) with the top and bottom quarks are given by

    $$\begin{aligned} y_{A_1 tt} = \frac{2 \lambda v \cos ^2 \beta }{\mu },\quad y_{A_1 bb} = \frac{2 \lambda v \sin ^2 \beta }{\mu }. \end{aligned}$$
    (17)

    In this scenario, both \(\lambda \) and \(\kappa \) can be large, yielding a sizeable contribution to the mass of the SM-like Higgs boson at tree level. Hence, a large value for \(\tan \beta \) is unnecessary. A scan in the parameter space in this scenario is performed using NMSSMTools 4.2.1 [88–93] including all built-in constraints, such as from Higgs searches, superpartner searches, muon \(g-2\), flavour physics, invisible Z-boson decay and the constraints from \(\Upsilon \) decays (with the exception of the Landau pole test and DM related-constraints, which are not considered). The resulting values for the \(y_{A_1tt}\) and \(y_{A_1bb}\) couplings are compared to the expected collider bounds in Fig. 8a. Depending on the parameter values, the magnitude of \(y_{A_1tt}\) in this scenario can be up to \(\sim 0.5\). Only for an integrated luminosity of 3000 fb\(^{-1}\) the LHC can probe a coupling \(y_{A_1tt}\) as small as 0.5 via the \(t\bar{t} A_1\), \(A_1 \rightarrow b\bar{b}\) channel. Therefore this scenario is difficult to probe, even at the HL-LHC.

  2. 2.

    PQ-limit: \(\frac{\kappa }{\lambda } \rightarrow 0\), \(A_\kappa \rightarrow 0\), where the theory is approximately invariant under the transformation

    $$\begin{aligned}&H_u \rightarrow H_u \exp (i \phi _\mathrm{PQ}),\quad H_d \rightarrow H_d \exp (i \phi _\mathrm{PQ}), \nonumber \\&\quad S \rightarrow S \exp (-2i \phi _\mathrm{PQ}), \end{aligned}$$
    (18)

    and the tree-level mass of the PQ pseudo-Goldstone boson \(A_1\) are given by

    $$\begin{aligned} m_{A_1} = - \frac{3 \kappa A_\kappa \mu }{\lambda }. \end{aligned}$$
    (19)

    This scenario has been proposed as a supersymmetric benchmark for sub-electroweak scale (singlino-like) DM [84], since its lightest neutralino is generically singlino-like and lighter than the electroweak scale. Particularly, in this scenario \(A_1\) can serve as the mediator for DM annihilation into a bottom quark pair and explain the diffuse gamma-ray excess from the Galactic Centre [15, 16]. In this limit, the tree-level couplings of \(A_1\) with the top and bottom quarks are given by

    $$\begin{aligned} y_{A_1 tt} = \frac{\lambda v \cos ^2 \beta }{\mu }, \quad y_{A_1 bb} = \frac{\lambda v \sin ^2 \beta }{\mu }, \end{aligned}$$
    (20)

    and so they are smaller by a factor of 2 than the corresponding couplings in the R-limit. Furthermore, a smaller \(\lambda \) is favoured in this limit and a relatively large \(\tan \beta \) is needed to generate a mass of 125 GeV for the SM-like Higgs boson. Therefore, the coupling \(y_{A_1 tt}\) tends to be smaller than in the R-limit scenario. The resulting values for the \(y_{A_1tt}\) and \(y_{A_1bb}\) couplings are compared to the expected collider bounds in Fig. 8b. For most of the points, the magnitude of \(y_{A_1 tt}\) is below 0.1, which renders this scenario extremely difficult to probe at the LHC using the \(t\bar{t}A_1\) channel.Footnote 4

Finally, we stress that the \(b\bar{b} A_1\) channel does not help much in probing the R- and PQ-limit scenarios. The sensitivities of both searches are suppressed by the mixture with the singlet. Even worse, the mixing is approximately \(\tan \beta \) enhanced, further suppressing the sensitivity of the \(b\bar{b} A_1\) in probing the large \(\tan \beta \) region in both scenarios (Fig. 9).

7 Conclusions

Searches for CP-odd scalars, as predicted by many extensions of the Standard Model and motivated by some recent astroparticle observations, are part of the core program of upcoming LHC runs at \(\sqrt{s}=13\) and 14 TeV. Searches at LEP and during Run 1 of the LHC at \(\sqrt{s}=7\) and 8 TeV have placed only weak constraints on the coupling strengths of CP-odd scalars with top and bottom quarks, or in their allowed mass range.

Using a simplified model approach for the signal, we have carried out a detailed study to evaluate the prospects at the LHC for probing scenarios with a CP-odd scalar with mass \(20 \le m_A < 100\) GeV, via the process \(pp \rightarrow t\bar{t}A\) with subsequent decay \(A\rightarrow b\bar{b}\). To separate the signal from the large background from \(t\bar{t} +\) jets production, we apply jet-substructure techniques, reconstructing the mass of the CP-odd scalar as the mass of a large-radius jet containing two b-tagged subjets. The chosen method allows for a so-called “bump hunt” over a fairly smooth background, and it may be the most promising strategy for searching for a CP-odd scalar with mass \({\lesssim } 50\) GeV, i.e. about twice the typical minimum \(p_\mathrm{{T}}\) cut for narrow jets used in standard LHC searches. A significant effort has been made to develop a semi-realistic experimental analysis, including a fairly complete description of systematic uncertainties. For more realistic estimates, sophisticated statistical tools are used to constrain in-situ the effect of systematic uncertainties, thus limiting their impact on the search sensitivity. We then derive expected upper limits on the production cross section times branching ratio using the CL\(_\mathrm{{s}}\) method.

In specific models, e.g. 2HDM or NMSSM, the coupling of the A boson with the top quark is related to other couplings in a well-defined way. Hence, the upper limits obtained on this coupling for a given mass \(m_A\) can be used to bound other couplings of these models indirectly or as input for a global coupling fit. We find that in a type-I and type-II 2HDM the LHC can constrain a large fraction of the \((m_A, \tan \beta )\) parameter space, including the region preferred to explain the diffuse gamma-ray excess from the Galactic Centre as dark-matter annihilation via a CP-odd scalar mediator which decays into \(b\bar{b}\). However, in the case of the NMSSM with a light CP-odd scalar, a Goldstone boson of either a spontaneously broken R- or PQ-symmetry, the LHC appears to have very limited sensitivity in probing these models. Hence, depending on the concrete embedding of the scalar sector into a UV-complete theory, the LHC can provide complementary information, not accessible at either indirect detection experiments or electron–positron colliders, on the existence of CP-odd scalars, their mass, and their couplings to third-generation fermions.