Extending LHC Coverage to Light Pseudoscalar Mediators and Coy Dark Sectors

Many dark matter models involving weakly interacting massive particles (WIMPs) feature new, relatively light pseudoscalars that mediate dark matter pair annihilation into Standard Model fermions. In particular, simple models of this type can explain the gamma ray excess originating in the Galactic Center as observed by the Fermi Large Area Telescope. In many cases the pseudoscalar's branching ratio into WIMPs is suppressed, making these states challenging to detect at colliders through standard dark matter searches. Here, we study the prospects for observing these light mediator states at the LHC without exploiting missing energy techniques. While existing searches effectively probe pseudoscalars with masses between 5 - 14 GeV and above 90 GeV, the LHC reach can be extended to cover much of the interesting parameter space in the intermediate 20 - 80 GeV mass range in which the mediator can have appreciable Yukawa-like couplings to Standard Model fermions but would have escaped detection by LEP and other experiments. Models explaining the Galactic Center excess via a light pseudoscalar mediator can give rise to a promising signal in this regime through the associated production of the mediator with bottom quarks while satisfying all other existing constraints. We perform an analysis of the backgrounds and trigger efficiencies, detailing the cuts that can be used to extract the signal. A significant portion of the otherwise unconstrained parameter space of these models can be conclusively tested at the 13 TeV LHC with 100 fb$^{-1}$, and we encourage the ATLAS and CMS collaborations to extend their existing searches to this mass range.


Introduction
Light, weakly-interacting massive particles (WIMPs) are a particularly compelling class of particle dark matter (DM) candidates. The case for WIMPs with masses close to the electroweak scale has been strengthened by recent observations of an excess in gamma rays originating from the Galactic Center (GC) by the Fermi Large Area Telescope [1][2][3][4][5][6][7][8][9][10]. This signal has garnered much recent attention, since its morphology closely resembles that expected from dark matter pair annihilation into bottom quarks [9,11], though other final states can also provide a good fit when systematics are properly taken into account [12]. Moreover, the signal suggests a WIMP annihilation rate close to that required in the early universe for a thermal relic to saturate the observed dark matter density [9], and the excess is difficult to explain in terms of astrophysical backgrounds alone [9,13]. This has led many to believe that the Fermi GC signal may represent the first (indirect) observation of dark matter to date. A common and well-motivated class of models that can explain the observed excess features dark matter annihilating through a light pseudoscalar with Yukawa-like couplings to Standard Model fermions [14][15][16][17]. For example, these states appear generically in two Higgs doublet models and their extensions [18], as well as pseudo Nambu Goldstone bosons associated with the spontaneous breaking of a new global symmetry [19][20][21]. Their couplings to Standard Model fermions can arise at tree-or loop-level (see e.g. Ref. [22] for an example with heavy vector-like fermions). Since they couple to the visible sector, such pseudoscalars can constitute a portal to the dark sector, mediating the annihilation of dark matter (DM) into SM final states [14,[23][24][25][26][27][28].
Understanding how dark matter interacts with the visible sector is a crucial part of the current dark matter program. Direct detection experiments [29][30][31][32] and the observation of a Standard Model-like 125 GeV Higgs with a small invisible decay width [33,34] have severely constrained Z-and Higgs boson-mediated scenarios [21]. As a result, much recent work has been devoted to studying various possibilities for new mediator particles coupling weakly to the Standard Model degrees of freedom. Of these possibilities, pseudoscalars stand apart for several reasons. For one, they do not predict sizable spin-independent direct detection signals, in contrast with scalar and vector mediators. Furthermore, current collider constraints on new pseudoscalar particles are generally weaker than those on new scalar and vector degrees of freedom [35,36].
If the GC excess is indeed a signal of dark matter annihilation, and if the annihilation is mediated by a new pseudoscalar particle, it is both important and timely to consider how one might probe such scenarios at colliders. Much progress has already been made on this front. Based on the topology and kinematics of the dominant dark matter annihilation channel, scenarios explaining the GC excess with pseudoscalar mediators can be grouped into roughly three types, each with distinct prospects for collider discovery: 1. Models which rely on dark matter annihilating into on-shell mediators [37][38][39][40][41][42]. In this case, the annihilation rate into SM fermions factorizes and the coupling of the pseudoscalar mediator to SM degrees of freedom can be very small. Prospects for direct collider searches are often dim in this case, but there may be other handles on these models provided by direct detection, as well as fixed target and other precision experiments [37][38][39][40].
2. Scenarios featuring a pseudoscalar mediator with a significant invisible branching fraction [22,25,26,[43][44][45][46]. This results in distinctive missing energy signatures at the LHC which can be effectively probed by bb+MET, mono-jet, and other existing and planned LHC searches, as studied in detail in e.g. Refs. [22,25,26,43]. 3. Scenarios in which the pseudoscalar mediator is expected to have a small branching fraction into dark matter particles [14-17, 27, 28]. This can occur when the coupling between the dark matter and the mediator is small relative to the coupling of the mediator to Standard Model degrees of freedom, or when on-shell decays of the mediator into WIMP pairs is not kinematically allowed. Such scenarios can be more difficult to probe directly at the LHC than case 2, since they lack a distinctive missing energy signature [14]. In concrete models of this type, rare Higgs decays can be constraining, however the resulting limits can be straightforwardly avoided in many instances, as can limits from LEP, the Tevatron, and B-physics experiments (see e.g. Refs. [15,16]). While a signal would arise in indirect detection experiments, it has been shown that the dark matter and mediator in this case might avoid detection elsewhere [14]. This rather grim scenario is appropriately known as "Coy dark matter".
In this study we will focus our attention on case 3 above, as it is a generic yet largely unconstrained possibility, as we discuss below. We will restrict our attention to light mediators, with masses below 90 GeV, as pseudoscalars with larger masses are already probed by existing LHC Higgs searches. Furthermore, light pseudoscalars are very attractive from the standpoint of the Galactic Center excess, since they can provide an efficient resonant annihilation channel for the light dark matter masses suggested by the signal and, in some cases, allow for a p-wave annihilation channel into pairs of mediators to drive down the relic abundance without violating constraints from dwarf spheroidal observations [21]. In this situation, on-shell decays of the pseudoscalar to pairs of dark matter particles are suppressed and WIMP production at the LHC through the mediator will be negligible. Our strategy will be to extend LHC coverage to such scenarios by probing the light pseudoscalar directly through its interactions with the Standard Model degrees of freedom. The discovery of such a new particle would constitute a great step forward in our understanding of the dark sector and open up many possibilities for further study, including more dedicated experiments to probe its coupling to dark matter directly.
As we discuss below, the GC excess can suggest an appreciable mediator coupling to down-type fermions. Consequently, we focus on the associated production of the mediator with a b-jet or bb pair. We will assume that the mediator couples to Standard Model fermions with strength proportional to their mass, as in models with minimal flavor violation (MFV). We find that, for a significant range of mediator masses and couplings consistent with the GC excess, a promising signal is predicted in the 1-2 b+a production modes, with a → τ + τ − . We also explore the possibility of a → µ + µ − decays, which is more promising for low masses and likely features lower systematic uncertainties. Existing searches for pseudoscalars motivated by the Minimal Supersymmetric Standard Model (MSSM) and Next-to-MSSM (NMSSM) currently probe mediator masses down to 90 GeV and in the low-mass region between 5-14 GeV. However, we find that coverage can be extended to pseudoscalars in the intermediate mass range (between 20-80 GeV), which are promising for explaining the GC excess and would have evaded detection by LEP. We encourage both ATLAS and CMS to expand their analysis to include this region. In this study, we detail the cuts and kinematic variables that can be used to reduce the large backgrounds and show the extent to which the parameter space in these models can be conclusively tested at the 13 TeV LHC with 100 fb −1 of integrated luminosity. We demonstrate this using a simplified model and show the application of our results to the otherwise unconstrained parameter space of the NMSSM that is consistent with the excess (the NMSSM can be mapped directly onto our model).
It is important to emphasize that, although we will focus on pseudoscalars mediating dark matter annihilation consistent with the GC signal, our study can be applied much more generally to any model featuring light mediators with significant coupling to isospindown Standard Model fermions. Since we assume that the invisible branching fraction of the pseudoscalar is small, our analysis of the predicted collider signatures does not depend on the pseudoscalar's coupling to dark matter, nor on the nature of the dark matter itself.
This study is organized as follows: in Sec. 2, we discuss the simplified model used for our analysis, its relationship to the GC excess, and the existing constraints on light pseudoscalars. The following section (Sec. 3), details the collider signatures of the new mediator, as well as the backgrounds and trigger efficiencies relevant for our analysis. Our results for the LHC discovery potential of light psuedoscalar mediators are presented and discussed in Sec. 4, with further details of the analysis contained in Appendices A, B, and C. We then apply these results to the NMSSM in Sec. 5, showing that the searches we propose here can cover much of the parameter space consistent with the excess and that is currently unconstrained by other experimental searches. Finally, we summarize and conclude in Sec. 6.

A Simplified Model
For our analysis, we follow Ref. [14] and consider a light pseudoscalar that couples to Dirac fermion dark matter, χ, and to Standard Model fermions, with effective Lagrangian where y i = m i /v are the SM Yukawa couplings, with v = 174 GeV. We have assumed that the pseudoscalar a couples to the SM fermions with strength proportional to their masses. The pseudoscalar couplings to up-and down-type fermions are further assumed to depend on the overall scaling factors, g u,d , which we take to be the same for all down-or up-type fermions 1 . These factors appear e.g. in Two Higgs Doublet models (2HDMs) and their extensions; in a Type II 2HDM, g d = 1/g u = tan β, where β is the ratio between the two SU (2) Higgs vacuum expectation values. With the addition of a singlet that mixes with the SU (2) doublets, the effective couplings become g u = cot β cos θ and g d = tan β cos θ, where θ is the mixing angle between the SU (2) and singlet pseudoscalars. Note that Ref. [14] considered the case in which g d = g u = 1. This situation is very difficult to probe at colliders. Explaining the Fermi GC signal with g u = g d = 1 can require rather large values of g DM , unless the annihilation is quite close to the s-channel resonance. Often in ultra-violet (UV) complete models, a sizable value for g DM occurs together with low mass WIMPs in parametric regions featuring a large invisible branching fraction of the Standard Model-like Higgs [23], which is not observed. On the other hand, for pseudoscalar-WIMP couplings that are not too large, the Galactic Center excess suggests enhanced couplings to down-type fermions, as we show below. This situation is much more promising from the standpoint of LHC searches and, in some cases, is not probed by existing searches.

Explaining the Excess
Given the Lagrangian in Eq. 2.1, the zero-temperature s-channel annihilation rate for dark matter through a pseudoscalar into SM fermion pair f ifi is where m i , N C,i are the mass and color factor of the decay states, g i is either g u,d depending on the fermion, and Γ a is the total width of the mediator. Throughout this study we will assume that the dominant DM annihilation channel is χχ → bb. This mode has received the most attention in explaining the GC excess, and although a recent analysis has pointed out that other channels can also explain the signal [12], annihilation into a bb pair still provides a very good fit to the data. There have been several recent developments in determining which annihilation channels, WIMP masses, and annihilation rates best fit the Fermi data. For the bb channel, most previous analyses had suggested that m χ should fall roughly in the range 35 GeV m χ 50 GeV with annihilation rate σv 2-6 × 10 −26 cm 3 /s [9,14] (the required annihilation rate for self-conjugate dark matter would be reduced by a factor of two relative to these values). However, there are large systematic uncertainties associated with the propagation of gamma rays in the Galactic Center that must be taken into account. The impact of these systematics was first studied in Ref. [13], and subsequently by Ref. [12], which performed a detailed analysis incorporating several different models for the diffuse gamma ray background supplied by the Fermi collaboration. The end result is that the range of WIMP masses and annihilation modes statistically consistent with the excess increased significantly once these systematic uncertainties were taken into account. In particular, the range for WIMP masses annihilating primarily into bb was extended to [12] 35 GeV m χ 165 GeV, χχ → bb.
for specific values of the annihilation rate. Across this mass window, the signal from the Galactic Center suggests a clear range of values for the coupling constants of the mediator to SM states for a given g DM in this setup. Since the annihilation cross-section and pseudoscalar width are dominated by down-type interactions (for BR(a → χχ) 1), the only significant parametric dependence is on g DM , g d , m χ and m a . The down-type scale factor g d required to explain the GC excess for m χ = 45 GeV and m χ = 145 GeV (close to the best fit mass for the Fermi model (d) from Ref. [12]) is shown by the bands on the left and right hand side of Fig. 1, respectively, as a function of m a The shaded regions are compatible with the signal, with the red (upper) regions in each band excluded by the recent dwarf spheroidal constraints from Fermi [48], and the yellow (lower) regions with an annihilation rate compatible with both the excess and the constraints. The upper bound on g d from existing LHC searches for a → τ + τ − is shown in blue.
for various values of g DM . The range of annihilation rates allowed in the low mass case (LHS) is taken from Ref. [13], while the allowed values in the high mass case (RHS) are taken from Ref. [12]. In both cases the local dark matter density is assumed to be ρ = 0.4 GeV/cm 3 . The preferred regions depend on J ≡J/J 0 , the ratio of the angularly-averaged integral over the line-of-sight dark matter density ρ DM (r), given bȳ to the canonical valueJ 0 . For the low mass (m χ = 45 GeV) case we take J = 1, while for m χ = 145 GeV we take J = 0.3, which is within the systematic uncertainties discussed in Ref. [12]. The latter choice allows for an annihilation rate close to the canonical thermal freeze-out value ( σv 4.4 × 10 −26 cm 3 /s for Dirac fermion dark matter [49]) and consistent with the Fermi signal while evading the dwarf spheroidal constraints, discussed below.
For reasonable choices for g DM , the value of g d must be quite large to account for the GC excess, unless the masses are tuned to fall very close to the resonance. In addition, reducing the χ abundance has the effect of increasing the preferred value of g d for a given g DM . The regions of parameter space with large g d , in many cases preferred by the signal, predict a significant mediator production cross section at the LHC in association with bottom quarks. Also, the pseudoscalar's invisible branching fraction is small across the entire parameter space, except for low g d and large g DM . For m a < 2m χ an on-shell pseudoscalar cannot decay to a pair of WIMPs, while for m a > 2m χ we find that BR(a → χχ) > 0.1 only for g d 4 for g DM = 0.1 in the m χ = 45 GeV case, since everywhere else g d (g DM ) is too large (small) for this decay to contribute appreciably to the total width.
It is important to note that the Fermi collaboration recently released updated limits on the dark matter annihilation rate from observations of dwarf spheroidal (dSph) galaxies [48]. The resulting constraints 2 are in mild tension with a dark matter explanation of the excess, however there is still a large amount of parameter space capable of explaining the GC excess that survives this constraint. This is shown in Fig. 1, in which the red bands show the impact of the dwarf spheroidal limits (points in these bands could potentially explain the excess but are excluded at 95% C.L.). Meanwhile the yellow bands show points consistent with both the GC excess and dSph limits. Note that in the high mass case all points consistent with the excess are compatible with the dSph constraints for our particular choice of J .
One concern may be that, since the recent dSph constraints disfavor larger annihilation rates, some points with light WIMP masses consistent with the GC excess and dSph limits will tend to produce too large a relic abundance. The dark matter relic density is set by the annihilation rate at finite temperature, which can differ from that at T = 0. In particular, for s-channel annihilation through a pseudoscalar with m a < 2m χ , the annihilation rate at T = 0 is greater than that at freeze-out (T f.o. ∼ m χ /20). The upper limit on the annihilation rate, set by the Fermi dSph results, is below the required annihilation rate at freeze-out for m χ 100 GeV, naively disfavoring this region. However, there are several well-known and straightforward exceptions to this reasoning [53]. For example, p-wave processes with contributions to the total annihilation rate scaling as v 2 DM (with v DM the relative dark matter velocity) will become important at freeze-out, increasing the annihilation rate at T f.o. but not altering the T = 0 prediction. An example of such a process generically expected along with light mediators is χχ → aa (this is another virtue of the light pseudoscalar scenario). Other scenarios allowing for an enhanced annihilation rate at T f.o. relative to that at late times include those with additional co-annihilation channels or featuring m a > 2m χ so that σv T =0 < σv T =T f.o. . Thus, although in some cases the dSph limits may result in requiring some additional tuning or model-building to achieve the correct DM relic abundance, dark matter explanations of the excess, particularly those involving s-channel annihilation through a relatively light pseudoscalar, are alive and well. Note that this tension largely disappears above m χ ∼ 100 GeV, since the dSph upper bound is above the canonical WIMP crosssection in this region (although one should verify that contributions to the annihilation rate at freeze-out from the other states in the theory do not over-dilute the relic density).
In summary, dark matter annihilating through a relatively light pseudoscalar can explain the Galactic Center excess and be compatible with the recent dwarf spheroidal limits from Fermi. In all most discussed and shown above we expect BR(a → χχ) 1, either because m a < 2m χ , g DM 1, or both. This implies a low likelihood of observing the pseudoscalar through missing energy signals at the LHC. In the following subsection, we describe some of the other existing constraints on the parameter space and highlight the need for direct coverage of these scenarios at the LHC.

Existing Constraints
Our goal will be to ascertain to what extent LHC searches can cover the parameter space shown in Fig. 1 that is not currently probed by LHC searches [54][55][56][57][58] . To our knowledge, there are currently no direct constraints on the parameter space of our simplified model with 15 GeV m a 90 GeV. By this, we mean that there do not exist constraints depending only on the pseudoscalar's coupling to SM fermions in this mass range. There are indeed several other indirect constraints, but these are inherently dependent on other degrees of freedom in the UV complete theory and can be straightforwardly avoided in many cases. We will present explicit examples of points evading all of the searches discussed below but still predicting an observable LHC signal in the NMSSM in Sec. 5.
Pseudoscalar mediators with GeV-scale masses predict highly suppressed direct dark matter detection cross-sections. At tree-level, the pseudoscalar only interacts spin-dependently with nuclei. Using the expressions and results found in Refs. [14,59], we find that the spindependent scattering cross-section for dark matter off of nuclei via the pseudoscalar is far below the reach of current and planned experiments (σ SD 10 −48 cm 2 ) across the parameter space we consider. This is thanks to the 1/m 4 a suppression in σ SD in this regime. Also, while spin-independent scattering can occur via one loop diagrams, this contribution is also much too small to be observed. The difficulty in observing dark matter interacting with the visible sector primarily through a pseudoscalar in direct detection experiments in indeed one of the main reasons such models are understood to be coy.
Light pseudoscalars can also be constrained by flavor observables 3 . Loop diagrams involving the pseudoscalar can generate effective flavor-changing vertices [21,61]. The limits are severe for pseudoscalars lighter than the B and Υ meson scale simply because the mediator can be produced on-shell in decays. For m a 10 GeV, the constraints are very significantly relaxed, with the most stringent arising from LHCb [62] and CMS [63] measurements of BR(B s → µ + µ − ). For m a m B , the limits are approximately considering the pseudoscalar contribution alone [21]. This constraint would naively appear to directly constrain some of the parameter space shown in Fig. 1, however, the new contributions to B s → µ + µ − are strongly model-dependent [64]. For example, in supersymmetric UV completions of our model, such as the NMSSM, there are several new contributions which enter with opposite sign to that from the a-induced vertex. Thus, cancellations can occur over large portions of the parameter space allowing for light pseudoscalars with large couplings to SM fermions (i.e. above the naive upper bound of Eq. 2.5) [16], once again highlighting the need for direct probes of this parameter space. For light mediators with 2m a < m h (h is the 125 GeV SM-like Higgs), exotic Higgs decays to pseudoscalar pairs can affect the Higgs width and signal rates [65,66], which are constrained by both ATLAS [67] and CMS [68]. Evidence for h → aa decays was also searched for at LEP [69] and the Tevatron [70]. Such decay modes can also be very effectively probed at the High Luminosity LHC [66]. Indeed, this has long been recognized as an important potential discovery channel of NMSSM pseudoscalars at colliders [71][72][73][74][75]. However, these constraints depend on the haa coupling which, in some cases, can be made appropriately small in realistic models [16,23], especially those in which the pseudoscalar coupling to Standard Model fermions does not arise through mixing with the SM-like Higgs [22]. Alternatively, simply taking m a > m h /2 avoids these constraints altogether.
Another indirect constraint arises from LEP searches for e + e − → ha production [76]. While these results prohibit MSSM-like pseudoscalars lighter than 90 GeV for all values of tan β, these bounds depend on the Zha coupling, which is model-dependent, and can again be straightforwardly avoided [16]. For example, in a Type II 2HDM with an additional singlet (2HDM+S), the Zh i a coupling scales as where S i1 and S i2 are the corresponding entries in the matrix diagonalizing the 3 × 3 CPeven mass matrix with the Higgs bosons ordered in mass (see e.g. Eq. 2.22 in Ref. [77]). Contrasting to g d ∼ cos θ tan β, we see that the simple limit cos θ 1, tan β 1 can result in an appreciable g d with a significantly suppressed g Zha .
Finally, existing MSSM Higgs boson searches at the Tevatron [78][79][80][81][82] and the LHC [54,55,57] constrain g d for m a > 90 GeV, but in an effort to avoid the large backgrounds encountered for lighter masses, and because LEP had already ruled out MSSM-like pseudoscalars with masses below 90 GeV, the published limits do not extend below the Z mass. There are also searches for light (m a 15 GeV) pseudoscalars at CMS [56], motivated by certain limits of the NMSSM. However, the 15 GeV m a 90 GeV mass range remains currently untested 4 .
Although the collider limits on a light pseudoscalar can be avoided, one might also be concerned about the consistency of this scenario once the model is UV-completed. Our Lagrangian is not invariant under SU (2) L × U (1) Y , and so, given a particular UV completion, one should also check that constraints on the other states can be satisfied while demanding a light pseudoscalar. In 2HD+S models, for example, most constraints on the rest of the Higgs sector can be satisfied by simply taking the charged Higgs mass to be moderately heavy (a few hundred GeV) with an appropriate choice of tan β [54,55]. Such requirements are consistent with light pseudoscalars and sizable g d , as shown in e.g. Ref. [16] and in Sec. 5 for the NMSSM.
Perhaps surprisingly, there is a significant gap in coverage for light pseudoscalars with appreciable couplings to SM fermions, as arise in models explaining the GC excess or otherwise. This situation has room for improvement. In the remaining portion of this paper, we will investigate to what extent searches similar to those already existing for heavy MSSM Higgs bosons and for light NMSSM pseudoscalars can directly probe the parameter space motivated by the Galactic Center excess. This task requires a careful treatment of the backgrounds below the Z mass. As we will show below, the backgrounds can be substantially reduced by using a suitably chosen sequence of kinematic cuts.

Production and Signals
Heavy neutral Higgs bosons in two Higgs doublet models are being searched for via a variety of experimental signatures, including gluon fusion (ggF) production, or production in association with top or bottom quarks [54,55]. These canonical Higgs-type searches become much more difficult below the Z threshold, where the backgrounds increase dramatically. Fortunately, as we have shown in Sec. 2 above, light pseudoscalar mediators consistent with the Galactic Center excess can have enhanced couplings to down-type Standard Model fermions relative to those expected for a Standard Model-like Higgs boson of the same mass. This results in an enhanced production cross-section in modes involving b quarks, and (potentially) in the gluon fusion channel relative to the Standard Model-like case. This situation is depicted on the left hand side of Fig. 2, which shows as an example the enhancement of both the inclusive bba (black) and gluon fusion (red) production cross-sections with g d = g −1 u (i.e. cos θ = 1), relative to those with g d = g u = 1 (σ 0 ), as a function of g d 5 .
The enhancement of the bba cross-section is independent of m a , as it only depends on g d for a given m a , while the differently styled red curves correspond to σ ggF /σ ggF,0 for different values of the pseudoscalar mass. The enhancement is substantially larger in the bba mode across the parameter space, which suggests focusing on production processes involving b quarks rather than the gluon fusion process. We consider the branching ratio of the pseudoscalar into various final states, assuming BR(a → χχ) is negligible, on the right hand side of Fig. 2. The pseudoscalar's branching = 1 as a function of g d . The dotted, dash-dotted, dashed, and solid red lines correspond to the enhancement in ggF production for m a = 20, 40, 60, and 80 GeV respectively. The corresponding enhancement for bb associated production is shown by the solid black curve (the enhancement is independent of m a ). Right: Branching fraction of the pseudoscalar into various final states (assuming BR(a → χχ) is negligible). Note that the branching ratios into fermions are nearly independent of g d (since the total width is set primarily by a → bb, τ + τ − decays), while the a → γγ partial width is substantially suppressed for g d > 1.
fraction into photons is small and is further suppressed for g d > 1 which, when combined with the increased backgrounds for m a < m Z , suggests that diphoton searches will likely be unable to probe the low-mass pseudoscalar mediators we are interested in. On the other hand, while the favored decay is into a bb pair, searches for such resonances would contend with large, pure QCD backgrounds to exploit this mode. Thus, to avoid large backgrounds while maintaining a reasonable signal, and to maximize the enhancement of the production cross-section, we propose a search for the pseudoscalar in second and third generation dilepton (τ + τ − and µ + µ − ) pair production in association with one or two b-jets. Of course this strategy requires that the pseudoscalar couples to leptons, which is typical in extended two Higgs doublet models, but need not be the case [21].
Similar searches have been considered by both ATLAS [54] and CMS [55], but are focused on higher mass resonances motivated by two Higgs doublet models and the MSSM, where the mass region of interest is greater than about 90 GeV [76] due to LEP searches and precision constraints on heavy Higgs bosons. Also, previous theoretical studies in the context of the NMSSM have investigated the potential for the LHC to probe light pseudoscalars with somewhat similar searches [83][84][85][86][87][88]. However, these investigations did not incorporate trigger and detector effects, and did not analyze the effects of cuts on the signal and backgrounds in detail, which is a major component of this work and crucial for obtaining an observable signal. While Ref. [87] arrives at largely negative conclusions regarding bba production (at least in the NMSSM with partial universality), our analysis suggests a much more positive picture once appropriate cuts are implemented.
It is worthwhile to point out that the CMS search in Ref. [56] finds sensitivity down to g d ∼ 3 for masses up to m a ∼ 14 GeV in the gluon fusion mode with a → µ + µ − . One might be inclined to conclude that this search channel could simply be extrapolated to larger masses in the scenarios of interest. However, this is unlikely to be the case. Fig. 2 shows that the gluon fusion production cross-section is actually suppressed for 1 < g d 10 as compared to its value with g d = 1, given our assumptions about the couplings. The suppression increases with m a and is due to the decreased top quark loop contribution that is otherwise dominant for heavier masses. In addition, due to the kinematic beta factor 1 − ( 2m i ma ) 2 , the bb branching ratio is suppressed for smaller values of m a , resulting in an increase in the µ + µ − branching fraction. For example, BR(a → µ + µ − ) is enhanced by almost a factor of 2 at m a = 10 GeV versus m a > 20 GeV. Thus, for the scenarios we consider, production modes involving downtype fermions at tree-level would appear more promising than those relying on gluon fusion production and decays to muons, although different assumptions about the coupling structure could alter this conclusion. For a related analysis of the potential LHC reach in the 0b mode in Z models, see e.g. Ref. [89].
In the remainder of this section, we discuss the challenges and strategies for examining low mass pseudoscalars with enhanced couplings to down-type fermions, g d > 1.
We implemented our simplified model in FeynRules 2.0 [90], and generated both our signals and backgrounds at leading order (LO) using MadGraph5+aMC@NLO [91]. We then used Pythia 6.4 [92] to decay the τ leptons and hadronize the b-jets, and incorporated initial and final state radiation, with an appropriate scale used for the MLM matching of hard element and radiated jets. Detector simulation for trigger and tagging was performed using Delphes 3.0 [93]. Trigger effects were implemented as step-function cuts at the analysis stage, though some minimum kinematic requirements were enforced at the generation phase.
Diagrams for some of the primary production modes for the signal are shown in Figure  3. To avoid the appearance of potentially large logarithms arising from the phase space integration over collinear final state quarks, the semi-inclusive b(b)a events were generated with b quarks included in the parton distribution functions (pdfs) of the proton. This is known as the "five flavor scheme" which effectively re-sums the large logarithms [94][95][96]. Exclusive bba events were generated without the inclusion of the b pdfs since the resulting contributions are doubly pdf-suppressed and subleading when compared to the gluon induced processes. To avoid double counting between the two-body, b(b)a, production and the three-body, bba, production mode where one of the b-jets is collinear with the proton beam, the three-body production mode was generated with a minimum p b T > 5 GeV. There are several technical difficulties associated with accurately calculating the twobody b(b)a production cross-section at hadron colliders, which have received much attention in the literature [97][98][99][100][101][102][103][104]. In particular, the leading order production cross-sections are known to exhibit a substantial dependence on the renormalization and factorization scales, µ r and µ f , respectively [99,103]. For our signal generation, we consider dynamic scales defined by Figure 3. Some of the diagrams contributing to the production of the pseudoscalar, a, at the LHC. The two rightmost diagrams arise in the 5FS.
where f is an overall scaling factor, and i refers to the produced b's and a. This is in keeping with previous analyses in the context of Standard Model-like Higgs production [99,[104][105][106][107].
We considered the impact of the scale dependence by varying the overall scaling factor in the range [1/2, 2], which resulted in a 2-20% change in the production cross section, with larger effects occurring for smaller values of m a . This is consistent with the range typically found in the literature [98,99,101].
To further validate the results of our leading order calculation, we have compared our LO result for the dominant (gb(b) → b(b)a) production mode to the next-to-leading order (NLO) result calculated in the five flavor scheme implemented in the program MCFM [108] for several choices for µ f,r (we neglect the difference between scalar and pseudoscalar production which are small [104]). We find that our LO results exhibit reasonable agreement with the NLO result, falling within a factor of 1-2 across the parameter space we consider. Additionally, there are theoretical uncertainties related to the specific choice of parton distribution functions, which have been shown to be of order ∼ 5% for low masses [103], as well as some residual renormalization scheme dependence (MadGraph uses an on-shell scheme, while e.g. MCFM uses MS). To account for these effects, Appendix C takes a conservative approach and explores the effect of a factor of 2 over-estimation in our signal and, separately, a factor of 2 under-estimation in the backgrounds. Our overall conclusions are not significantly affected by this re-scaling, and so we believe them to be quite robust.
For an experimental search, we consider three possible leptonic tagging channels: SR1 requires one electron and one muon; SR2 requires one lepton (e or µ) and one hadronic τ ; SR3 requires two muons. SR1 is motivated by excellent trigger response, while SR2 is motivated by the larger branching ratios and SR3 is motivated by a resonance search methodology in the di-muon invariant mass spectrum that allows for the use of data-driven backgrounds. In all three signal regions, we also require 1-2 b-jet tags, and no light jets, where light jets are defined as p T > 40 GeV. The signals are therefore inclusive for light jets with p T < 40 GeV, such as those that are commonly generated from ISR effects. These tagging requirements significantly suppress fake backgrounds arising from vector boson production in association with light jets.
We assume the default CMS tagging efficiencies that are implemented in Delphes 3.0, which are as follows. For tagging, electrons are required to have p T > 10 GeV and |η| < 2.5. Within the inner region of the detector, |η| < 1.5, we assume a tagging efficiency of e = 0.95, while for the outer region but with |η| < 2.5 we assume e = 0.85. The rate at which jets fake electrons is taken to be j e = 0.0001 and uniform over the whole detector. For muons, we require that candidates have p T > 10 GeV and |η| < 2.4. Since our analysis involves only low p T muons, we take a fixed tagging efficiency of µ = 0.95, which is appropriate for p µ T < 1000 GeV. For the tagging of hadronic taus, we require |η| < 2.5 and take a fixed tagging efficiency of τ = 0.4 with a fake rate for mistagging a light jet as a hadronic tau of j τ = 0.001.

Trigger Effects
Since the signal typically produces very soft jets and leptons, trigger effects are very important to consider. To account for the effect that trigger has on our results, we have implemented a variety of triggers as a step-function cut based on what we believe are reasonable off-line triggers for CMS 6 . The following primary triggers are potentially relevant to our study: • 1e: single electron with p T > 35 GeV; • 1µ: single muon with p T > 25 GeV; • 2µ: di-muon leading with p T > 17 GeV, subleading p T > 10 GeV; • eτ h : electron + hadronic tau with p τ T > 45 GeV, p e T > 19 GeV; • µτ h : muon + hadronic tau with p τ T > 40 GeV, p e T > 15 GeV; • eµ: leading electron + muon with p e T > 23 GeV, p µ T > 10 GeV; • µe: electron + leading muon with p e T > 12 GeV, p µ T > 23 GeV; We also include other triggers, such as those involving photons, jets, τ h plus MET, and b-jets, but these provide a negligible effect on the signal events (i.e. < 0.3% of signal events pass all the non-primary triggers combined) and so are not included in the above list. The nonprimary triggers pass a significant portion of the backgrounds, however, which necessitates their inclusion, but this indicates that these events have distinctive signatures that can be eliminated from the analysis by kinematic cuts.
Due to the low mass of the pseudoscalar in our search, a significant number of the production events will not pass the trigger. Since we are not privy to the details of the final triggers, we consider the effect of varying the muon p T thresholds for the triggers that include a primary muon. These triggers have the greatest likelihood for discretionary variation in a dedicated experimental search, and are the most important due to having the lowest inclusive cross sections and thus p T thresholds. We analyzed the cross section of signal events that pass each of the primary trigger cuts (σ Ty SRx ) as a fraction of the cross section of generated events (σ gen SRx = σ gen × BR(τ + τ − → SRx) × SRx ) with the same tagging signature in each signal region, independently: where SRx refers to the signal region and T y refers to the specific trigger. This ratio can be considered as a sort of trigger efficiency. Of note, we found that the e + τ h and µ + τ h triggers did not pass any of a preliminary 200k generated events, likely due to the hard cut on the p T of the τ h and the low mass of the pseudoscalar. Since the hadronic tau has a large fake background from mistagged light jets, we do not anticipate that the trigger threshold for hadronic tau p T will be improved enough to make these triggers worthwhile to consider. While the fake rate of jets for electrons is smaller than for hadronic taus, we believe it is unlikely that any significant improvement in the electron trigger thresholds will be implemented as there would still be a larger increase in the inclusive cross section than for similar changes in the muon trigger thresholds.
The summary of the trigger efficiency ratios in Eq. 3.2 for the default implemented triggers is shown in Table 1, while an analysis of the effect of varying the threshold for the muon p T in the 1µ, 2µ and µe triggers for each of the three signal regions is shown in Figure  4. A naïve interpretation of this figure suggests that the single muon trigger includes a larger fraction of the signal than µe or 2µ triggers, but it is important to note that the single muon inclusive cross section at the LHC is significantly larger than the muon+electron or dimuon inclusive cross sections, and thus will typically have a higher p T threshold than the other triggers and a lower trigger efficiency, as shown in Table 1.

Backgrounds and Their Reduction
Since the QCD backgrounds at the LHC are significant, the fake rate of jets as electrons, hadronic τ -jets, and b-jets are important to take into account. Additionally, backgrounds with similar kinematics to the signal we examine produce soft leptons that may not be identified as easily or may fall outside of the central region of the detectors where tagging is possible. Thus, backgrounds producing more than two leptons, where one is not tagged, may contribute to the signal regions. To account for these effects, we include backgrounds that produce between one and three leptons (e, µ, τ ), and 0-2 b-jets, in association with 1-3 light jets (with n b + n j ≤ 3), since our signal is inclusive to low p T light jets. The following background processes are generated:  Table 1. The ratio of cross section that passes the trigger cut to the generated cross section for 200k generated events. Kinematic dependent tagging efficiencies are already incorporated into the cross sections. All leptons (e, µ, τ ) are generated with a minimum p T > 10 GeV, but τ decays to leptons can result in a p e,µ T < 10 GeV. The columns in this table are not necessarily independent, as it is possible for an event to simultaneously pass multiple triggers.
• pp → γ * /Z + (0, 1, 2, 3)j, γ * /Z → + − ; • pp → W ± + bb + (0, 1)j, W ± → ± ν (ν ); • pp → W ± + b(b) + (0, 1, 2)j, W ± → ± ν (ν ); • pp → W ± + (0, 1, 2, 3)j, W ± → ± ν (ν ); where = (e, µ, τ ) and j are light jets (u, d, s, c, g) that can come from associated production. Each entry in the list above is produced with the number of quoted jets, and MLM matching and merging is incorporated to avoid double counting of the light jet production with the  . Trigger ratios for each signal region, normalized to the produced and tagged cross section, based on varying the leading muon p T . The µe a trigger is based on a subleading electron p T = 12 GeV, while µe b is based on a subleading electron p T = 17 GeV. The 2µ trigger for SR3 is based on a subleading muon p T = 15 GeV rather than the p T = 10 GeV discussed in the text, as the trigger response for the lower subleading p T is very similar to the single muon rate due to the minimum p T settings in the event generation stage and tagging thresholds. initial state radiation (MLM matching with XQCU T = 15 and QCU T = 20). The two largest contributions to our backgrounds are the inclusive Z production modes (included in the first three entries) and tt (included in the seventh entry), but these are effectively reduced by kinematic cuts. The kinematic distributions of the signal and backgrounds are included in Appendix A.
Based on the kinematic distributions we examined, we have identified a number of possible cuts that improve the signal significance. These cuts are focused on reducing the tt and Z +nj backgrounds. The tt and other backgrounds with W + W − lepton production can be reduced with cuts that involve the E T measurement, including a direct E T cut, as well as the transverse mass m T = 2p 2nd T E T (1 − cos θ). Backgrounds with a Z resonance can be reduced by a cut on the dilepton invariant mass, m . In addition, a large fraction of the backgrounds producing both leptons and jets have a large total p T . Thus, we consider cuts on the scalar sum In the case of SR3, dilepton invariant mass cuts are implemented in a fixed range. While the branching ratio to dimuons is small (< 0.1%), the a → τ + τ − → µ + µ − + E T branching ratio is similarly small, and the invariant mass peak of the direct decay is reconstructible with low smearing. Thus, it may be possible to observe the pseudoscalar with a resonance search methodology. For SR3, we consider only events within a 2-3 GeV invariant mass bin centred at the mass of the pseudoscalar. In contrast, the analysis for SR1 and SR2 are based on a cut-and-count methodology, since the dilepton peak is significant smeared out due to the loss of information from the neutrinos originating from the τ decays. For these signal regions, we do not employ a narrow invariant mass window and instead employ m cuts to exclude backgrounds only.
The cuts for SR1 and SR2 are considered separately in each of two distinct scenarios: hard cuts are better for high luminosity searches and have a greater overall reach, while soft cuts are better for low luminosity searches. Kinematic threshold values for the considered cuts were chosen by maximizing σ sig * L/ σ sig * L + σ bkg * L + 2 sys σ 2 bkg * L 2 , for a systematic uncertainty of sys = 0.2 and luminosity of L = 100/fb, while maintaining σ cut sig /σ tot sig ∼ 0.5(0.8) for hard (soft) cuts for m a > 40 GeV. The dimuon signal region, SR3, is analyzed assuming only a single cut scenario, as background events with m ∼ m a generally have similar acceptance rates to the signal.
The expected search reach using these cuts is given in the next section. Further details about the acceptance rates for each cut are provided in Appendix B. Alternative approaches for determining the cut regions, such as those incorporating repeated algorithmic refinements of the phase space, would optimize cuts for a single mass value and be unable to account for the full range of parameters we explore. Maximizing the acceptance rate for m a = 40 GeV would result in a poorer reach in g d values for m a = 80 GeV, for example. We feel our approach is more appropriate for a general search strategy.

Results
We can now investigate the extent to which the light pseudoscalar parameter space consistent with the Fermi signal can be probed by the searches we propose.  Due to the low pseudoscalar mass region of interest in this study, as well as the cut-andcount search method for SR1 and SR2, systematic uncertainties are a particularly challenging aspect of performing this search. To estimate the effect of systematic uncertainties, we consider two scenarios in addition to our two cut (hard/soft) scenarios -low systematics, with sys = 10%, and high systematics, with sys = 30%. Our analysis of the discovery potential  is based on a signal significance, given by where N s = σ s * L and N b = σ b * L are the number of signal and background events,  respectively, after cuts for a given integrated luminosity, L. Contours of constant luminosity are plotted in Figures 5, 6 and 7. For small enough values of g d , systematic uncertainties dominate the signal, and we expect that greater luminosity will be insufficient to illuminate any signal. Note that we have also verified that each signal data set considered has at least 5 events after cuts. The soft cut scenarios of SR1 and SR2 are optimized for early searches with low luminosity, but suffer from a larger systematics-dominated region, since the total backgrounds are much larger. Thus, their ability to exclude the parameter space ends at approximately L = 10/fb integrated luminosity. Alternatively, hard cuts scenarios have a better reach with exclusions from L = 100/fb, though larger luminosity will be unlikely to push this boundary any further.
As discussed, the expected sensitivites for each case are affected by three primary components: production, trigger and cuts. Production rates decrease with increasing mass, m a , reducing the overall cross section and number of events at the LHC. In contrast to this, trigger response improves for heavier pseudoscalars, but has a significant effect on the lighter pseudoscalar scenario. However, the pseudoscalar is produced in association with b quarks, which results in a boost to the a that allows a large enough fraction of events to pass trigger and thereby make the search viable. Lastly, eliminating backgrounds resulting from the Z peak results in a choice of cut thresholds that has a larger impact on events from heavier pseudoscalar masses, especially for the hard cut scenarios. These issues combined result in the typical shape observed in Figures 5 and 6, with reduced exclusion reach for both the lowest and highest mass scenarios.
The dimuon search uses a different approach, incorporating a pseudo-resonance search methodology. While we do not fit a line-shape over the background and compare the signal, we employ a narrow invariant mass window with a sliding center that effectively estimates the result from such an approach. In practice, an approach that fits a line to the continuum background will reduce systematic uncertainties that are associated with the cut-and-count methodology, which requires simulations to estimate the backgrounds. As a result, we suspect that the low systematics scenario in Figure 7 is potentially a more realistic case, in contrast with the other signal regions, where low systematics may be overly optimistic.
As a result of the relatively large width of the SM Z, combined with detector smearing effects, a dimuon resonance at 80 GeV will contend with increased backgrounds from the Z peak (which is why we do not consider heavier masses). If we assume similar systematic uncertainties for each signal type, then the most promising reach for the high m a region is in the 1e1µ signal regions, while the 1 1τ signal regions are more promising for the low m a regime. Note that the reach in the dimuon signal region is not as promising as the others for any part of the parameter space under the assumption of similar systematics. As mentioned, however, systematic uncertainties in the dimuon search will likely be smaller than in the other modes, and so all signal regions combine to form a complimentary and robust search strategy.
Comparing Figures 1 and 5-7, we see that the searches we propose will cover a significant portion of the otherwise unconstrained parameter space consistent with the Galactic Center excess in scenarios with light pseudoscalar mediators, even with rather low integrated luminosity. This region is both theoretically and phenomenologically well-motivated, and we encourage both ATLAS and CMS to consider searches along the lines of those presented here.

Application to the NMSSM
To illustrate the usefulness of our results in a UV-complete model, we can consider how our searches impact the Z 3 -symmetric NMSSM parameter space consistent with the excess. To set our conventions, we take the superpotential to be with soft supersymmetry breaking terms given by For sizable tan β and cos θ not too small, g d will be larger than 1. Our conventions follow those found in Refs. [23,77], to which we refer the Reader for further details regarding the spectrum.
There have been two scenarios proposed in the Z 3 -invariant NMSSM to explain the GC excess involving neutralino annihilation into SM particles through a light singlet-like pseudoscalar [16,17] (see also Ref. [46] for an analysis of the general NMSSM, which in some cases may also be probed by the searches we present). The first involves a mixed singlino/Higgsino-like neutralino, which, to achieve a Standard Model-(and not singlet-) like 125 GeV Higgs, requires the lightest pseduoscalar to be a nearly pure singlet [16] (i.e. cos θ 1). Since the pseudoscalar couplings to SM fermions are suppressed, to explain the GC excess this scenario requires m a ≈ 2m χ to within about a GeV precision, as well as additional Z-mediated contributions to the annihilation rate in the early universe to drive down the relic abundance. This would seem quite finely tuned, requiring a fortunate conspiracy of parameters to achieve. Instead, we focus on the second possibility, namely that the neutralino is bino/Higgsino-like. In this case, the singlet component of the 125 GeV Higgs is naturally small, and so the lightest pseudoscalar can feature a more significant amount of mixing between the singlet and SU (2) states. As a result, the requirement that the neutralino annihilation is on resonance is relaxed, allowing one to consider a much larger range of masses not precisely tuned to m a ≈ 2m χ [16].
It is worth mentioning that analyses of the NMSSM subsequent to Ref. [16] have found somewhat different results, favoring the singlino/Higgsino scenario [17,109]. However, taking the systematics into account in fitting the Fermi signal [12,13], we find that the bino/Higgsino scenario is fully compatible with both the GC signal and the Fermi dwarf spheroidal limits. Another reason the bino/Higgsino scenario may have been disfavored in Ref. [109] is that the large pseudoscalar couplings to the SM fermions in the bino/Higgsino scenario are constrained by rare meson decays, in particular B s → µ + µ − . As pointed out in Ref. [16], these constraints can be avoided rather straightforwardly by taking advantage of mild cancellations between the various SUSY contributions to BR(B s → µ + µ − ). Such points can be difficult to sample in a large scan of the parameter space, as employed in Refs. [17,109]. However, we have verified that the bino/Higgsino scenario is still in fact viable when taking these constraints into account, as claimed in Ref. [16].
The bino/Higgsino explanation for the GC excess maps directly onto our simplified model (only that the WIMP is a Majorana, instead of Dirac, fermion). To illustrate the effect of our searches on the viable bino/Higgsino parameter space of the NMSSM, we performed a Markov Chain Monte Carlo scan of the parameter space using NMSSMTools 4.4.0 [110], interfaced with micrOmegas 3.1 [111]. Motivated by the parameter space presented in Ref. [16], we  Figure 8. Application of our results to the Z 3 -symmetric NMSSM. The black (gray) contours correspond to the reach at 100 fb −1 (1 fb −1 ) for the hard (soft) cut scenarios and low systematics in the various search channels. The gray points are the result of a Markov Chain Monte Carlo scan of the parameter space (described in the text) consistent with all existing phenomenological constraints with no requirements on the LSP relic abundance or annihilation rate with parameters as in Eq. 5.5 and m A = 550 GeV. The green, blue, and orange points correspond to points capable of explaining the Fermi signal and consistent with the recent dwarf spheroidal constraints for m A = 500, 550, and 600 GeV, respectively. The red band is an example of the NMSSM parameter space found to be consistent with the excess in Ref. [16]. The sample point of Table 2 below is indicated with a star. Note that it may be possible to choose parameters minimizing the haa coupling to fill in the m a < m h /2, g d > 1 region, which we did not attempt in our scan.
with all other soft masses and triscalar couplings at 1 TeV, while varying tan β, κ, and A κ . We required all points to satisfy all existing constraints discussed earlier and implemented in NMSSMTools. The results of the scans are shown, along with our results for the LHC reach across the parameter space, in Fig. 8. The gray points were generated without requiring the lightest supersymmetric particle (LSP) to explain the Galactic Center excess or satisfy constraints on its relic abundance. The green, blue, and orange points correspond to m A = 500, 550, 600 GeV and feature a bino-like LSP with a relic abundance compatible with WMAP  Table 2. Example parameter space point in the NMSSM capable of explaining the GC excess and consistent with the Fermi dwarf spheroidal limits. All dimensionful parameters are in GeV unless otherwise stated. The remaining parameters are set to the values shown in Eq. 5.5. This point would likely be probed by the searches we propose at the 13 TeV LHC with 100 fb −1 of integrated luminosity.
and Planck measurements (including a 2σ theoretical uncertainty) [23] 0.091 ≤ Ωh 2 ≤ 0.138 (5.6) and compatible with both the Galactic Center excess and the dwarf constraints, for self-conjugae dark matter. Points satisfying these constraints typically have small, but non-negligible, p-wave suppressed contributions at freeze-out, such as those involving the Z (but still consistent with limits on the invisible Z width). This slightly reduces the relic abundance relative to the value suggested by χχ → a → bb annihilation alone and allows these points to circumvent the dSph limits. Note that we did not attempt to minimize the haa coupling, and so no points were found with 2m a < m h and g d > 1. However, it might be possible to reach this parametric regime [23] as suggested in Ref. [16], whose results we show along with ours in Fig. 8 by the red band. These values were taken from Fig. 6 of Ref. [16] for m χ = 35 GeV, while our scan was performed assuming m χ ≈ M 1 = 45 GeV. Table 2 provides the detailed spectrum information for an example parameter space point consistent with the GC excess and which would be probed by a → τ + τ − , µ + µ − at the 13 TeV LHC. This point is marked by the black star in Fig. 8. Note also that our scan did not find points with g d > 18. Larger values of g d are typically excluded by LHC limits on the heavy MSSM-like pseudoscalar for the values of tan β sampled. In theories that do not rely on mixing with the SM-like Higgs, these constraints, as well as those from h → aa decays, are often significantly relaxed or absent. The contours in Fig. 8 show the sensitivity of our proposed searches to the NMSSM parameter space consistent with the Galactic Center excess at both 1 fb −1 and 100 fb −1 . A significant portion of the favored region with sizable g d would be probed by the 13 TeV LHC at these luminosities. Even more reach would be expected at the 14 TeV LHC. Our searches are complementary to h → aa observations as well as existing LHC searches for MSSM Higgs bosons and would access regions of the parameter space not currently probed by other experiments, providing a potential window into a dark sector difficult to access otherwise.

Summary and Conclusions
Many dark matter models feature WIMPs that can be very difficult to observe at colliders. Scenarios of this type can be consistent with the Galactic Center excess observed by the Fermi Large Area Telescope. Exploring these "coy dark sectors" at the LHC suggests a shift away from missing transverse energy signals and towards direct signatures of the particle(s) mediating the interaction of the dark matter with the Standard Model.
Models involving pseudoscalar mediators and consistent with the GC excess can be of the coy variety. A good fit to the Fermi signal can be provided by relatively light WIMPs annihilating through a pseudoscalar into b quarks. In many realistic scenarios this suggests substantial couplings of the mediator to down-type Standard Model fermions. The signal favors WIMP masses in excess of ∼ 35 GeV, while current collider bounds often imply pseudoscalar masses below 90 GeV (provided they satisfy constraints from LEP). An interesting and currently untested explanation of the GC signal thus involves a pseudoscalar with mass below about 90 GeV with sizable couplings to down-type fermions and small branching fraction into WIMPs. The latter is generically small in this scenario since the on-shell decay of the mediator into dark matter is often kinematically forbidden and because the pseudoscalar's coupling to WIMPs is relatively small. Our study has attempted to extend LHC coverage to this scenario by taking advantage of the mediator's enhanced couplings to Standard Model fermions (relative to those of a SM-like Higgs boson of the same mass) and studying the production and decays of the pseudoscalar involving down-type final states.
To this end, we explored signals that include one to two b-jets and with either τ or µ lepton pairs in the final state. We employed a simplified model, in which we assumed that the couplings of the pseudoscalar to Standard Model fermions were proportional to their mass, modulo common scaling factors for down-and up-type fermions. While this need not be the case, this situation is common in UV completions involving Type II 2 Higgs doublet models, as in supersymmetry. Our results can be applied to models with different coupling structures by a straightforward re-scaling of the production cross-section and branching ratios.
Due to the rather low pseudoscalar masses we consider, trigger is an important factor in the search reach. We thus performed an analysis of the trigger response of the signal, and explored cuts that were effective in improving the signal significance. Our search strategy comprises a signal excess analysis for the 1e1µ + 1 − 2b and 1 1τ + 1 − 2b modes, including low luminosity (soft cuts) and high luminosity (hard cuts) scenarios, and a dilepton resonance search in the µ + µ − + 1 − 2b signal. Since signal excess searches suffer from large systematics from comparisons to simulated instead of data driven backgrounds, we also analyzed the impact of systematic uncertainties on the LHC reach in all three signal modes.
In the most optimistic scenarios, we find that the LHC should be able to explore values of the reduced pseudoscalar coupling to down-type fermions as low as g d ∼ 8 for 100/fb of integrated luminosity at √ s =13 TeV. Even in more pessimistic scenarios with higher systematics, we find that the LHC should be able to explore down to g d ∼ 10 for some values of m a . This reach, however, is highly dependent on the trigger settings, and so we strongly recommend that the experimental collaborations attempt to account for this type of signal when finalizing their trigger thresholds for leptons, particularly those triggers for muons. The parameter space in the NMSSM not covered by h → aa searches, with m a ∼ 60 − 80 GeV, should be explorable to some extent, and further optimization of the search strategy could focus on this narrow region of masses. More generally, the searches we propose are highly complementary to those already existing at the LHC or elsewhere, highlighting their importance in the interest of fully covering the parameter space in question.
In summary, light pseudoscalars with significant couplings to Standard Model fermions are well-motivated mediators for dark matter annihilation and arise in many models, including those explaining the Fermi Galactic Center excess. In many cases, these new particles would have evaded previous searches but should be testable at the LHC. Significant regions of the parameter space can be explored even with low luminosity, and so this signal presents a possibility for ongoing examination throughout the full LHC program.

B Cut Flow Matrices
When examining the potential of enhancing the visibility of the signal through cuts, we considered a variety of possible kinematic variable distributions, some of which are shown in Appendix A. Of those examined, we chose to consider only those cuts in which the shape of the backgrounds was distinctly different than the shape of the signal for at least one of the signal regions so that cuts on the background had a larger fractional effect on the backgrounds than on the signal. The variables that most effectively improved the signal significance were E T , p T of the leading lepton (p T ), dilepton mass (m ), total scalar sum of visible momenta (H T ), transverse mass of the subleading lepton (m T 2nd ), and the scalar sum of the lepton p T and E T ( H T ).
Variables with E T components were most effective at eliminating backgrounds containing decays of W bosons, including m T 2nd where backgrounds containing intermediate W 's have a longer tail on the distribution. Of note, we found that the transverse mass distribution based on the leading lepton p T had a longer tail for the signal, and so was not quite as effective.
Since E T , for example, is a component of multiple cuts, we consider the correlation between the events passing each pair of cuts in cut flow matrices in Tables 3 through 7. Diagonal entries are the acceptance rate for the single cut labeled in both the column and row headers, where red text indicates background acceptance rates and black text indicates signal acceptance rates. Each off-diagonal entry in these tables represents the acceptance rate (A) of the cut labeled by the row (r) header on the events remaining after performing the cut in the column (c) header, such that each entry is given by For example, the upper right most entry of Table 3 shows an 87.4% acceptance rate for background events and 84.3% acceptance rate for signal events from applying the E T cut to the pool of events that already passed the H T cut. The lower left most entry shows that 17.2% of background events and 87.7% of signal events pass the H T cut after applying the E T cut.
Since the E T cut removes a similar number of events for the signal and background once the H T cut has been applied, the E T cut is superfluous once the H T cut has been applied, and thus should not be included in the final set of cuts. In fact, the E T distribution for signals and backgrounds have similar shapes once the H T cut has been applied, and thus no E T cut value will be effective.  Table 4. Cut flow matrix for SR1: 1e1µ+1−2b+0j signal with soft cuts. The cuts are: E T < 50 GeV, p 1 T < 40 GeV, 12 < m < 45 GeV, H T < 140 GeV, M 2 T < 40, H T < 120. The red/top entry in each cell is the acceptance rate for all backgrounds combined, while the black/bottom entry shows the acceptance rate for the signal. The total cross sections for the 1e1µ + 1 − 2b + 0j signal after applying the trigger cuts are σ bkg = 2187 fb and σ sig = 60.4 fb (m a = 60 GeV, g d = 25).  Table 5. Cut flow matrix for SR2: 1 1τ +1−2b+0j signal with hard cuts. The cuts are: E T < 30 GeV, p 1 T < 40 GeV, 12 < m < 45 GeV, H T < 130 GeV, M 2 T < 25, H T < 100. The red/top entry in each cell is the acceptance rate for all backgrounds combined, while the black/bottom entry shows the acceptance rate for the signal. The total cross sections for the 1 1τ + 1 − 2b + 0j signal after applying the trigger cuts are σ bkg = 742 fb and σ sig = 84.3 fb (m a = 60 GeV, g d = 25).  Table 6. Cut flow matrix for SR2: 1 1τ +1−2b+0j signal with soft cuts. The cuts are: E T < 55 GeV, p 1 T < 55 GeV, 12 < m < 60 GeV, H T < 190 GeV, M 2 T < 45, H T < 140. The red/top entry in each cell is the acceptance rate for all backgrounds combined, while the black/bottom entry shows the acceptance rate for the signal. The total cross sections for the 1 1τ + 1 − 2b + 0j signal after applying the trigger cuts are σ bkg = 742 fb and σ sig = 84.3 fb (m a = 60 GeV, g d = 25  Table 7. Cut flow matrix for SR3: 2µ + 1 − 2b + 0j signal. The cuts are: E T < 60 GeV, p 1 T < 50 GeV, H T < 120 GeV, M 2 T < 45, H T < 120. The red/top entry in each cell is the acceptance rate for all backgrounds combined, while the black/bottom entry shows the acceptance rate for the signal. The total cross sections for the 2µ + 1 − 2b + 0j signal after applying the trigger cuts are σ bkg = 7249 fb and σ sig = 108 fb (m a = 60 GeV, g d = 25).

C Variation of Exclusions
As discussed in section 3, calculations of signal events using the 5FS is quite strongly dependent on the factorization and renormalization scales used. In MadGraph5, we employed a dynamic scale scheme that we then varied by an overall scaling factor between 0.5 and 2.0 (see Eq. 3.1). This factor had the largest effect for low mass pseudoscalar calculations, with a factor of 0.5 reducing the total cross section by approximately 22% for m a = 20 GeV, while only reducing the total cross section by a factor of 4% at m a = 80 GeV. Alternatively, the authors of [103] use a fixed renormalization and factorization scale scheme based on the sum of the masses of the pseudoscalar and the on-shell b quark masses. Variations of this scale by a factor between 0.5 and 2.0 results in a cross section reduced by as much as 50%.
In addition, our calculations of the backgrounds were performed at leading order. Higher order effects, as well as possible unaccounted-for experimental issues, may result in larger backgrounds than we predict. In order to address concerns regarding these two issues, we explore much more conservative contours determined by performing the same calculations but with a factor of 2.0 larger backgrounds, and separately with a factor of 0.5 smaller signal. Figures 18, 19 and 20 give these results. Of note, many regions of parameter space are still explorable at the LHC with 100/fb of integrated luminosity even in the more pessimistic scenarios.  The black lines represent the boundary of the systematics dominated region, the red lines represent the discovery potential at L=10/fb, while the yellow lines represent the discovery potential for L=1/fb.  The black lines represent the boundary of the systematics dominated region, the red lines represent the discovery potential at L =10/fb, while the yellow lines represent the discovery potential for L =1/fb. and conservative factors applied to the signal (dotted) and backgrounds (dashed) (solid lines show the original bounds without any factor applied to the signal or background). Contours correspond to constant values of log(L × fb) needed to achieve k = 3. The black lines represent the boundary of the systematics dominated region, the red lines represent the discovery potential at L =10/fb, while the yellow lines represent the discovery potential for L =1/fb.