A study of collider signatures for two Higgs doublet models with a Pseudoscalar mediator to Dark Matter

Two Higgs doublet models with an additional pseudoscalar particle coupling to the Standard Model and to a new stable, neutral particle, provide an attractive and fairly minimal route to solving the problem of Dark Matter. They have been the subject of several searches at the LHC. We study the impact of existing LHC measurements on such models, first in the benchmark regions addressed by searches and then after relaxing some of their assumptions and broadening the parameter ranges considered. In each case we study how the new parameters change the potentially visible signatures at the LHC, and identify which of these signatures should already have had a significant impact on existing measurements. This allows us to set some first constraints on a number of so far unstudied scenarios.


Introduction
A Dark Matter (DM) model involving two Higgs doublets and an additional pseudoscalar mediator [1,2] (2HDM+a) has been the subject of several searches at the LHC [3][4][5]. This model provides the simplest theoretically consistent extension of DM simplified models with pseudoscalar mediators [6]. It contains a number of more-or-less free parameters, essentially masses and mixing angles of the bosons and a coupling to the DM candidate. This leads to a particularly rich phenomenology, dominated by the production of the lightest pseudoscalar or the heavier Higgs boson partners via loop-induced gluon fusion, associated production with heavy-flavour quarks or associated production with a Standard Model (SM) Higgs or gauge boson. Depending on the values of the parameters, very diverse characteristic signatures can be produced, including the traditional DM signature of missing transverse momentum (E miss T ), but also a range of SM-only final states. The relative importance of these final states varies strongly as the parameters of the model change. This leads to numerous potentially observable final states at colliders, many of which remain largely unexplored. The purpose of this paper is to explore some of them using existing LHC measurements.
For a given scenario, one expects a well-designed search to give optimal sensitivity to a targeted signature. However, because of the rich particle content of the model there may be contributions from other signatures, even in the benchmark scenarios typically considered in searches. Some of these signatures may contribute to 'SM-like' cross sections which have already been measured 1 . In this study we use Contur [7] to examine the sensitivity of ATLAS, CMS and LHCb particle-level (i.e. unfolded) measurements available in Rivet 3.1 [8] to a two Higgs-doublet model with a pseudoscalar mediator and DM particle.
The paper is structured as follows. First, we discuss the Herwig calculations used, the advantages and limitations of the Contur method, and the default parameters of the model. Then we revisit the benchmark scenarios proposed by the LPCC DM Working Group [2] and report the sensitivity of the measurements to these 2 . We then extend the explored parameter space away from the benchmark region; first, we relax the assumption that the exotic Higgs bosons are all degenerate in mass. Next, we vary the mass of the DM particle and the mediator, with a view to exploring regions where the correct DM relic density can be obtained from this model alone. After that, we vary both the mixing angle between the two neutral CP-odd weak eigenstates (sin θ) and the ratio of the vacuum expectation values of the Higgs doublets (tan β). Finally, we relax the strict imposition of the alignment limit.
In each case we study how the new parameters change the potentially visible signatures at the LHC, and identify which of these signatures should already have had a significant impact on existing measurements. This allows us to set some first constraints on a number of so far unstudied scenarios, and to suggest some future directions for study. These are summarised in the final section.

Methodology
We use Contur to scan over various parameter planes of the model, generating all 2 → 2 processes in which a Beyond the Standard Model (BSM) particle is either an outgoing leg, or an s-channel resonance, with Herwig [10]. This procedure omits some potentially significant 2 → 3 processes, such as pp → ttX and pp → thX, where X is a BSM particle. These processes were studied using dedicated runs with Madgraph5 aMC@NLO [14], and found not to contribute significant additional sensitivity 3 . The calculations are leading order and Herwig includes also leading order loop processes that dominate the production for the DM mediator and the exotics Higgs bosons. It factorises the production of the particles from their decay using the narrow width approximation, so does not include interference terms with SM processes.
No matching or merging between the BSM matrix elements and the parton shower is being used in Herwig 4 . Therefore, because the parton shower can add jets above some minimum transverse momentum k min ⊥ , there is an element of double counting between, for example, gg → H → XY diagrams (where one or both of X, Y are BSM states) with an additional hard gluon radiation from the parton shower, and gg → Hg diagrams where the H decays to XY . The default value of k min ⊥ in Herwig is 20 GeV. A scan of this value over 10 GeV to 160 GeV indicated a relatively rapid fall in sensitivity between 10 GeV and 20 GeV, and a slower fall above this, presumably driven initially by the reduction in double-counted events and then by the loss of valid events as k min ⊥ increases to higher values which the parton shower will not populate. Since for this BSM model we work at leading order, we set k min ⊥ = 50 GeV to be on the conservative side.
All the cross sections quoted have been calculated with Herwig 7 [12,13] and with Mad-graph5 aMC@NLO [14] and were found in good agreement.
We note that interference effects between SM tt production [15][16][17][18][19] have a significant effect on the shape of the mass distribution of the A → tt decay, for example, and as already mentioned these effects are not taken into account by Herwig. However, none of the measurements used are specifically hunting for resonances; in general they are inclusive W +jet or top measurements, and so are unlikely to be very sensitive to differences in shape.
Contur identifies parameter points for which an observably significant number of events would have entered the fiducial phase space of the measurements, and evaluates the discrepancy this would have caused under the assumption that the measured values, which have all been shown to be consistent with the SM, are identical to it. This is used to derive an exclusion for each parameter point, taking into account correlations between experimental uncertainties where available. The speed of the Contur method, the inclusive approach of Herwig in generating all processes leading to BSM particle production, and the variety of measurements available, allow us to broaden the range of parameters considered into some interesting regions, away from the usual benchmarks. Nevertheless, the starting point is the default parameters settings derived from [2], and detailed in Table 1. In each section that follows, the parameters are as given in this Table unless explicitly stated otherwise. , the mass of the DM candidate is M χ , the coupling of the pseudoscalar mediator a to χ is g χ d , θ is the mixing between the two neutral CP-odd weak eigenstates, and sin(β − α) is the sine of the difference of the mixing angles in the scalar potential containing only the Higgs doublets. Finally, λ i are the quartic couplings of the Higgs potential.

Benchmarks with degenerate exotic scalar masses
We first focus on two parameter scans, those of Fig. 19 in Reference [4]: the (M a , M A ) scan and the (M a , tan β) scan. In these scans, which follow the recommendations of the LPCC DMWG [2], the masses of all the exotic Higgs bosons (A, H, H ± ) are degenerate, the mass of the DM candidate M χ = 10 GeV, the coupling of the pseudoscalar mediator a to χ is unity, sin θ = 0.35 where θ is the mixing between the two neutral CP-odd weak eigenstates, and we set sin(β − α), the sine of the difference of the mixing angles in the scalar potential containing only the Higgs doublets, to unity, meaning we are in the aligned limit so that the lightest mass eigenstate has SM Higgs couplings, and the quartic couplings are all set to λ i = 3 (see Table 1). The results are shown in Fig. 1.
In the (M a , M A ) scan (Fig. 1a, 1c) the overall sensitivity is the combined result of relatively marginal (1-2 σ) contributions to a wide range of cross-section measurements. One of the few measurements involving missing energy, which was in fact targetted at ZZ → νν + − [20] and uses only 7 TeV data, is sensitive at low M A and M a < M H − M Z , as can be seen in Fig. 1c. This region overlaps, but is smaller than, the dilepton + E miss T searches discussed in [4], which uses more luminosity and higher-energy collision data. The sensitivity is principally due to the H → aZ decay, which has a branching fraction of ≈ 30% in this region.
Over much of the rest of the 95% excluded region, ATLAS jet substructure measurements [21] are most sensitive, especially in the hadronic W selection; however, several other measurements, notably those aimed at W or Z plus jets, make comparable contributions. These final states generally arise from the production of all the exotic Higgs bosons, with subsequent decays either to top quarks or directly to W bosons 5 .
In the (M a , tan β) scan (Fig. 1b, 1d), most of the plane is excluded for tan β < 1. This is again largely due to analyses involving W boson production, especially now the CMS measurement of semi-leptonic tt decays, which give rise to lepton-plus-E miss T -plus-jet final states [22]. These events again come from decays of the exotic Higgs bosons. For example for tan β = 1, M A = M H = M H ± = 435 GeV and M a = 250 GeV, the production cross-section for the CP-odd Higgs boson A is about 6 pb and its dominant decay channels are in tt (85%) and χχ (15%). In addition, the production cross section for the CP-even Higgs boson H is 3 pb, and it decays into a top-pair with a branching ratio of 88%, while the second dominant decay mode is H → aZ (10%).
The (M a , tan β) scan presented in Fig. 1 has been extended to consider a wider range for the tan β parameter with respect to Ref. [4]. It is worth mentioning that for tan β > 50, the  width of the CP even Higgs boson predicted by this model exceeds 20% of its mass for M a < 350 GeV. As discussed in the previous section, the use of the narrow width approximation means the calculation of the BSM signal will be increasingly unreliable in these regions. Overall, the SM measurements are highly complementary to searches in constraining the parameter space, and extend the coverage for tan β > 10(40) for M a = 100(500) GeV. The sensitivity in this region of the parameter space is driven mainly by E miss T +jet [23,34] final states. This sensitivity comes from bb and gb initiated production, which is enhanced for higher tan β values. There are also contributions from processes such as tH ± → ttb [24] which contribute to top final states, for example boosted top pairs [25].

Non-degenerate masses
The imposition of M H = M A = M H ± is a somewhat arbitrary choice; while one might expect none of the Higgs masses to be far from the electroweak symmetry-breaking scale, some variation may occur and this can have a significant impact on collider signatures. Relaxing this hypothesis was previously studied in the context of exotic Higgs searches in ttZ and tW b final states [26]. There are other constraints on the masses however, as discussed in Ref. [1,26]. Flavour and perturbativity constraints apply at low values of tan β, with a constraint that tan β ≥ 0.8 for M H ± = 750 GeV from B-meson observables, and weaker constraints from LHC searches. In the following, tan β = 1 will be imposed (unless it is varied in a scan), so these constraints are always satisfied. Electroweak precision measurements also imply that H and A can only differ significantly in mass if either M H = M H ± or M A = M H ± ; both of these cases will be studied separately. In the case of M A = M H ± , there are additional constraints on the value of sin θ from the same source, generally favouring lower values. Finally, as discussed in [26], when the A and H masses differ by several 100 GeV, the width of the heavier of the two will get large due to decays to H ± W ∓ .
In this work, we will study the impact of the degeneracy assumptions using three benchmarks: All three benchmarks follow the LPCC DMWG convention for the other parameters (see Table 1).

Case 1
In [1], studies were performed for parameter points in which the degeneracy of the exotic Higgs boson masses was broken. These studies interpreted a number of LHC searches, including mono-jet, mono-Higgs and invisible Higgs searches. Their "Scenario 4" corresponds closely to our first benchmark above, the only difference being the fact that the quartic couplings were set to zero, which has negligible impact on the results. The comparable Contur scan is shown in Figure 2, along with the original figure. Again, the absence of many E miss T measurements, especially mono-Higgs, reduces the sensitivity, visible in this case at low M a and moderate tan β. However, E miss T +jet, top, and W measurements do disfavour the region tan β 1. There is also some sensitivity at high tan β (where our scan is extended to higher values) coming from Z, diphoton and E miss T measurements.

Case 2
Continuing with the charged Higgs mass M H ± equal to the heavy CP-even Higgs mass M H , we now scan over the range 200 GeV to 2200 GeV, and scan M A independently over the same range. This two-dimensional scan was repeated for a few values of M a in the range 100 GeV to 500 GeV, and the results are shown in Fig. 3.
At low M a (Fig. 3a) the model is disfavoured in the M H ± = M H < 500 GeV region for all values of M A , principally due to the + E miss T measurement [20]. This comes from exotic Higgs decays to aZ, with the a then decaying to DM. As M a is increased, this sensitivity is reduced, as these decays are suppressed, and H → bb, a signature with larger SM background, dominates.
At low M A there is reasonable sensitivity for all the M a values considered, coming from a variety of signatures. At low M H ± = M H , dilepton+(b−)jet measurements [28,29] contribute; in this region the A decays to HZ about a third of the time, with the H decaying mostly to bb. There are also contributions from missing energy and a lepton -essentially, top and W +jet cross sections; the charged Higgs bosons decay dominantly to tb at low M A for all the M a values considered. Finally, for M a not too high, low M A and high M H , diphoton analyses have strong sensitivity; in this region the H decays dominantly to ah, with the h → γγ [30] (and also four-lepton final states from h → 4 [31]) becoming visible.
Over the rest of the plane, where the overall sensitivity is low, ATLAS jet substructure measurements [21] and E miss T +jet measurement [23] seem to offer the best chance of an eventual observation.

Case 3
The results of the scan in which the charged Higgs is degenerate with the heavy pseudoscalar A are shown in Fig. 4. The general features of the exclusion are similar to those in Section 4.2, in that the sensitivity is greatest when at least one of the new bosons has a mass of around 500 GeV or below. There are differences in the detailed reasons for this, however.
The dilepton+E miss T signature still plays a role at low M a and M H , although at high M H ± = M H the W or Z+jet signatures become more sensitive. The Higgs-to-diphoton measurements also play a role again, not only at low M A as before, but also at low M H , since A → ah makes a significant contribution.
As in Section 4.2, over the rest of the plane, where the overall sensitivity is low, ATLAS jet substructure measurements [21] and E miss T +jet measurements [23] receive the largest contributions to their cross sections.    Table 1.   Table 1.

Varying the DM mass
In cosmological models where the relic density of DM is determined at thermal freeze-out, and in the absence of further additional scattering mechanisms and particles, the parameters of the model determine this density. Astrophysical observations favour a value of Ωh 2 = 0.12 [32,33]. The most relevant parameters in calculating Ωh 2 from the model are the masses of the DM, M χ , and of the mediator, M a , although other parameters and the relations between them do have an influence. As discussed in [2], for the benchmark choice of M χ = 10 GeV this constraint disfavours the benchmark. However, for other values of M χ the correct relic density can be obtained without significantly influencing the collider phenomenology.
In Figure 5a, we show the benchmark scan in M χ -M a , where in [2] it was shown that for most values of M a , the correct relic density is obtained for M χ values of around 200 GeV. While very little of the plane is currently excluded at 95% c.l., there is sensitivity at the 68% level in vector-boson-plus-jet measurements (including the hadronic decays studied in the ATLAS jet substructure measurement), indicating that future precision measurements will have a significant impact. Figure 5b shows a scan over tan β for M χ values between 1 GeV and 1 TeV. In this plane (for the chosen values of the other parameters, particularly M a = 250 GeV), M χ needs to be between 100 GeV to 150 GeV to achieve Ωh 2 = 0.12. We see that the sensitivity of the measurements shows little dependence on M χ , and values in the range 0.5 < tan β < 30 are still allowed by the data. It is interesting to compare to the M a scan, Fig. 1d, where though the exclusion in tan β is similar, and at low tan β is due to the same final states, at M χ > 100 GeV or so the high tan β exclusion is due to boson+jet rather than E miss T +jet final states.

Varying the mixing parameters
Varying the mixing angle, sin θ, between the 2HDM pseudoscalar A and the mediator a can change the interplay between different signatures, in particular those that couple to top quarks. The DMWG benchmark parameters focus mostly on an intermediate mixing (sin θ = 0.35). In this section we investigate instead the case of maximal mixing (sin θ = 1/ √ 2), as well as scanning down to low values of mixing, where the DM candidate gradually decouples.
First we repeat the scans discussed in Section 3, keeping all parameter choices fixed, with the exception of the mixing angle. The (M a , M A ) scan is shown in Figures 6a,c. The overall sensitivity of SM model measurements for this slice of the parameter space is increased towards higher M A masses for small and intermediate M a masses. The ZZ → 2 + E miss T measurement [20] still plays an important role for M H = M A < 500 GeV, especially due to the fact that the H → aZ branching ratio increases with sin θ. For example, at M H = 435 GeV, M a = 100 GeV it approximately doubles as sin θ changes 0.35 → 1/ √ 2. In addition the production cross section for the H boson increases from 3.4 pb (sin θ = 0.35) to 3.6 pb (sin θ = 1/ √ 2). A second important analysis that dominates the sensitivity for large M A − M a values is the h → γγ analysis [30]. This analysis selects events where a Higgs boson is radiated from the pseudoscalar mediator a, a process whose cross-section is proportional to the mixing and the M A − M a mass difference.  Table 1.
The last but not least important signature that acquires additional exclusion sensitivity in the maximal mixing case is the jets+E miss T final state. Particularly strong sensitivity is obtained by the search for supersymmetric quarks in the 0-lepton channel [34] around M A ∼ 1 TeV. This sensitivity arises from the increased gluon-initiated production cross section for a (up to a factor 3-4 for M a = 100 GeV, assuming M A = 1 TeV), as well as the increase in branching ratio of many decays that contribute to the jets+E miss T final state. In this category we have not only the classic a/A + j → χχ + j but also H → aZ, H ± → aW ± or A → ah(bb).
The (M a , tan β) scan is shown in Figure 6b,d. Also in this case, the jets+E miss T final states play an important role for small a masses for all tan β values but the overall exclusion sensitivity in this region is a combination of multiple final states with similar sensitivities.
In Figure 7, we show a series of scans from minimal to maximal mixing, for various values of M a and tan β, with M a = 400 GeV.
For intermediate tan β, the overall sensitivity shows little dependence on sin θ, presumably because for these parameters it is not driven by the E miss T signature of DM. The E miss T +jet final state, which is more sensitive at higher BSM Higgs masses, does decline in prominence as the mixing, and hence the DM production cross section, declines. For high and low tan β, the sensitivity extends to larger values of M A , and increases with sin θ, again driven by E miss T +jet and jet substructure measurements.

Away from the alignment limit
Finally, the SM Higgs measurements and electroweak fits do still allow some mixing between the SM Higgs and the H. Figure 8 shows a scan where sin(β − α) is allowed to move from the alignment limit; for some limited parameter settings (around 0.5 < tan β < 10 and M H ± ≥ 600 GeV) values as low as 1 √ 2 are still allowed by the studies in [2]. The combined measurement sensitivity again shows little dependence on sin(β − α), although as we move away from the alignment limit, the diphoton [30,35] and four-lepton measurements [31,36] play an increasing role, since these decay channels open up for the H. The impact of the measurements of the SM Higgs branching fractions, and of any implied decrease in the SM Higgs cross sections, is not considered here.

Conclusions
There is interesting sensitivity to the two Higgs doublet plus pseudoscalar DM model across several final states already measured by ATLAS and CMS at the LHC. This is due to the quite complex phenomenology of the model, which can change a lot when the parameters change. The Contur approach allows us to efficiently scan a wider range of parameters than usually considered, relaxing some of the current assumptions imposed in benchmark scenarios.
Future searches for this model should also consider final states involving top and/or W production, even in the absence of a large missing energy signature. Producing cross section measurements together with such searches would increase the impact and generality of the results, as well as allowing predictions for the relevant SM processes to be probed and tested. If future W W measurements can be made less model-dependent, for example by including all relevant kinematic cuts in th fiducial cross section definitions, they may also make a significant     [4] in the (left) (M a , M A ) and (right) (M a , tan β) planes, but considering only 2 → 3 processes generated with Mad-graph5 aMC@NLO. The top row shows the sensitivity map in terms of CLs, with 95% (68%) excluded contours indicated by solid (dashed) white lines. The bottom row shows the analysis which contributes the most to the sensitivity of each signal point.