Comparing 2HDM $+$ Scalar and Pseudoscalar Simplified Models at LHC

In this work we compare the current experimental LHC limits of the 2HDM $+$ scalar and pseudoscalar for the $t \bar{t}$, mono-$Z$ and mono-$h$ signatures and forecast the reach of future LHC upgrades for the mono-$Z$ channel. Furthermore, we comment on the possibility, in case of a signal detection, to discriminate between the two models. The 2HDM+S and 2HDM+PS are two notable examples of the so-called next generation of Dark Matter Simplified Models. They allow for a renormalizable coupling of fermionic, Standard Model singlet, Dark Matter with a two Higgs doublet sector, through the mixing of the latter with a scalar or pseudoscalar singlet.


Introduction
The search of Dark Matter (DM) has become since some years one of the primary objectives of the LHC collaboration and represents, as well, one of the main motivations for the proposal of new collider facilities. From a theoretical perspective, it is crucial to provide an efficient interface for the interpretation of the outcome of collider searches and to enforce the complementarity with the other DM search strategies, namely Direct and Indirect Detection.
The so-called simplified models (see e.g. [1,2] for some reviews) were born to satisfy this need. Their validity might be, however, questionable. The reason is mostly twofold: first of all their extreme simplicity might render these models overconstrained and, moreover, associate them a very limited set of experimental signatures. The second reason is the fact that these models often lack of a theoretical completion, rendering their predictions potentially not reliable (see e.g. [3][4][5][6][7]).
For this reason the community is actively looking for a new generation of DM models [8] which can account for a rich collider phenomenology and solid theoretical predictions still within a not too large set of free parameters.
The so called 2HDM+S and 2HDM+PS are two notable examples of this last category of models. First of all, they allow for a renormalizable coupling between a fermionic gauge singlet DM candidate and the Standard Model (SM). The DM couples, on first instance, with a (pseudo)scalar singlet. A portal with the SM sector is created by the mass mixing of the latter with the CP-even (CP-odd) states of a two doublet Higgs sector. The 2HDM+S and 2HDM+PS provide, moreover, a broad variety of collider signals, beyond simple missing energy signatures (see e.g. [8][9][10][11]).
The aim of this paper is to study some of these particular signatures, namely tt resonances, mono-Z and mono-h. Although, all of these signatures appear in both models, there can be a sizeable difference for the expected signal rates in the pseudoscalar and scalar model. Therefore, in our work we characterize the similarities and differences of the different signatures in the two models, by looking at the limits derived from current LHC data and analyses. In addition, we also look at the reach of the mono-Z channel for future luminosities and comment on the possibility to discriminate between the 2HDM+S and 2HDM+PS models in case of a future signal detection.
The paper is structured as follows. In the next section we will illustrate, from a theoretical perspective, the two models. Section 3 is, instead, devoted to a brief overview of their present phenomenological constraints. The results of our numerical study will be presented in section 4. In Section 5 we finally draw our conclusions.

Scalar Potential
The most general scalar potential we consider is where V 2hdm is the standard 2HDM potential, V S is the potential of the scalar singlet and V S2hdm is the potential involving interaction terms between the scalar singlet and the doublets. We start by reviewing the 2HDM part of the potential. In 2HDMs, one has two doublets with identical charges, therefore having the freedom to choose a specific base Φ 1 , Φ 2 in terms of which to write the potential. A generic SU (2) change of base, where U is an SU (2) matrix, will result in a potential giving the same physics but with different coefficients for the different terms. While it would be possible to get rid of residual freedom through basis independent methods, see e.g. [16][17][18], we will nevertheless consider, in this paper, two reference bases in which to write the potential. The first one is what we will call the flavour base. Assuming the Yukawa interactions of the 2HDM are either of type I, II, X, Y or Inert, the basis Φ 1 , Φ 2 is defined in terms of which of the doublets interacts with the various fermions, and is therefore well defined. This basis is very useful to analyse the interactions of the scalars with fermions. The second basis we choose is the so-called Higgs basis, where one of the two doublets has no vacuum expectation value (vev). When writing the potential, we will label λ i , M ij the coefficients of the potential in the flavour base, and we will useλ i ,M ij to label the coefficients of the potential in the Higgs basis. The doublets in the flavour basis will be labelled Φ 1 , Φ 2 , while the doublets in the Higgs basis will be labelled Φ h , Φ H . The 2HDM potential reads in the flavour basis where a Z 2 symmetry Φ 1 → Φ 1 , Φ 2 → −Φ 2 has been imposed to suppress flavour changing neutral currents, thus removing the additional terms λ 6,7 containing three Φ 1 doublets and one Φ 2 doublet and vice versa. We still allow for the soft-breaking term M 2 12 that is necessary to have a scalar mass spectrum with the desired physical features. We also assume that the potential does not break CP explicitly. This means (M 2 12 ) 2 and λ 5 need to have the same phase (up to a negative relative sign), which can be reabsorbed into Φ 2 by a field redefinition. Therefore we will assume that the parameters λ i , M ij are all real. In this case, the SU (2) transformation freedom for the doublets reduces to SO(2) rotations: The most general expressions for the other two terms in Eq. (2.1) are The Z 2 symmetry might be also extended to these terms of the potential by defining suitable transformation properties for S. By choosing S → −S the terms µ S , µ 11S , µ 22S , would be forbidden, however they would still be allowed in case of soft-breaking of the Z 2 . The λ 12S term on the other hand is forbidden due to the assumed Z 2 symmetry of the Higgs fields. All these terms would have, however, a negligible impact for the phenomenology discussed in this work, since they would affect only the scalar trilinear and quartic interactions. Consequently, in line with [10], we will not include them in our analysis. In Eq. (2.6), all coefficients need to be real, except for λ 12S and µ 12S . As we have set λ 12S = 0, the only remaining physical phase is the one of µ 12S . To conserve CP, this coefficient must be either purely real or imaginary (after rotating the 2HDM doublets to have all phases contained in the 2HDM potential reabsorbed). In the first case, mixing between the CP-even scalar states contained in the Φ 1 and Φ 2 doubles and the singlet scalar would be obtained. In the latter case instead, the singlet state would mix with the neutral CP-odd states contained in the doublets.
The 2HDM+PS, defined in [12], and discussed in [10], is obtained using the latter option. Relabeling S to P to evidence more explicitly its CP-odd nature, the resulting potential is where V 2HDM is the one given in Eq. (2.3) while and v 1,2 are usually parametrized in terms of (2.11) ρ 1,2 are CP-even scalar particles, and η 1,2,3 are CP-odd scalar particles, but not necessarily mass eigenstates.
It is convenient to rotate 14) The two Higgs doublets are then defined by where G ± , G 0 are the SM Goldstone bosons, and H ± is the 2HDM charged scalar. All these particles are mass eigenstates. The particlesρ 1,2 are CP-even scalars that are linear combinations of the standard model Higgs boson h and an additional heavy scalar H.ρ 3 , η 3 are, instead, CP-odd and their combination will originate mass eigenstates which will be labelled a and A with, in general, M a < M A . Note that, in general, the terms λ 6 , λ 7 , λ hHP arise when changing base. Moreover, the coefficient of the term P iΦ † h Φ H + h.c. does not change when changing base. For the scalar case, it is possible to proceed in the same way, assuming µ 12S to be purely real. However, in the case of the scalar, there is also another way to obtain a mixing between the singlet and the CP-even scalars of the doublets. One can in fact assume that the singlet develops a vev. This is the approach taken by [5,11]. In this approach, one assumes that the scalar potential has a spontaneously broken Z 2 symmetry w.r.t. which only the particle S is odd, that gets rid of all cubic terms, independently of the basis in which the potential is written in. This may arise naturally, for example, in the case where S is part of a complex scalar charged under a dark U (1) gauge group. In [11], the authors write the potential directly in the Higgs basis aŝ this time, on top of the SM Goldstone bosons G and the charged scalar H ± , there is a single CP-odd scalar, A, and three CP-even scalars,ρ 1,2,3 , that mix together, and are linear combinations of h and S 1,2 , which are two new CP-even scalar particles. Note that, in this potential written in the Higgs base, the termsλ 6,7 are absent, contrary to the pseudoscalar case, where changing to the Higgs base switches these terms on, in general.

Mass Spectrum and Alignment Limit
The M h = 125 GeV Higgs boson is experimentally observed to be very similar to the SM one. To avoid most of the constraints from Higgs physics, we choose to work in the so called alignment limit, i.e. we impose specific relations among the parameters of the scalar potential such that the mixing between theρ 1 andρ 2 is negligible andρ 1 is always identified as the experimentally observed 125 GeV SM-like CP-even Higgs.
In the pseudoscalar model, the alignment limit can be achieved by assuming a specific value for β: Alternatively, one can assume the presence of the CP2 symmetry [19], that gives the following relations between the parameters In the latter approach, one has the advantage that, when switching base, one has λ i =λ i for i = 1, 2, 3, 4, 5 andλ 6,7 = 0.
In case of the scalar model, one also needs the additional assumption (2.28) In our convention, similar to other authors [10,11], we are concentrating on the case M a < M A . The original set of parameters λ 1,2,4,5 , M 2 11 , M 2 22 , M 2 12 , M 2 P P , µ 12P can therefore be expressed in terms of M h , M H , M H ± , M A , M a , θ, tan β, v together with the alignment condition. The parameters λ 3 , λ 11P , λ 22P , λ P remain free if the alignment condition Eq. (2.25) is chosen, while when choosing Eq. (2.26) only λ 11P , λ 22P , λ P remain free, while the value of λ 3 gets fixed.
In the scalar model instead, the mass matrix is made of a diagonal block containing two zero eigenvalues, corresponding to the Goldstone bosons G 0 , G ± , and the masses M 2 h , M 2 H ± , M 2 A , plus a two-dimensional block. Diagonalizing the block, one gets the masses and the mixing angle of the 2 CP-even scalars Also for the scalar case, we will concentrate on the case M S 2 < M S 1 . Again, one can exchange the set of parametersλ To avoid having the couplings λ i varying in an uncontrolled way when varying tan β, throughout the rest of the paper we will adopt the alignment condition Eq. (2.26), together with the mass degenerate assumption M H = M H ± = M A > M a or M A = M H ± = M S 1 > M S 2 , respectively for the pseudoscalar/scalar model. Moreover, to get a reasonable comparison between the pseudoscalar and the scalar model 1 and to avoid variations with tan β, we decide to set λ hhS =λ HHS = 0, (2.30) λ hhP =λ HHP = 0, (2.31) or equivalently for the rest of the paper. This is most relevant for the hAa vertex (hS 1 S 2 vertex for the scalar model), which appears in mono-h processes and depends on Finally, the parameter λ P,S is not relevant for our signatures.

Yukawa Sector
The Yukawa interactions of the SM fermions with the Higgs doublets can be expressed as As in standard 2HDMs, we shall need to choose Yukawa structures that keep potentially dangerous flavour violating processes under control. We outline the possibilities below, and explore the DM phenomenology of these choices in Sec. 3.1 by determining direct detection constraints. Rewriting Eq. (2.36) in the Higgs basis we have where the matrices Y u,d,l h,ij have to be the SM Yukawa matrices, namely while the Yukawa matrices of the additional doublet are assumed to be proportional to the SM ones where the i are Yukawa scaling factors, with i = u, d, l. This Yukawa structure is the so-called Aligned Yukawa model [20][21][22][23][24], which satisfies Natural Flavour Conservation. In special cases where the i satisfy certain relationships, the Aligned Yukawa structure can correspond to one of the Z 2 symmetric Yukawa structures (type I, II, X or Y), as shown in Tab. 1. As we will only probe values of tan β ≤ 3, our constraints will be valid for all Yukawa structures included in Tab. 1, except the Inert one.
The mass spectrum of the model contains an additional fermion, the DM candidate. Assuming it to be a SM singlet, it will couple, in the interactions basis, only with the scalar, 40) or the pseudoscalar singlet Due to the mass mixing between the SU(2) singlet and doublet scalars, the DM will in any case couple to the CP-even or CP-odd physical states according to whether we will consider the 2HDM+S or 2HDM+PS model.
Type I, II, X and Y To suppress flavour-changing neutral currents in 2HDM, it is possible to assume the presence of a Z 2 symmetry on the Yukawa sector, allowing only one of the two doublets Φ 1,2 to couple to a certain type of quarks and leptons. This hypothesis goes under the name of natural flavour conservation (NFC). The presence of the charged scalar H ± still allows FCNC to appear at loop level. Loop generated FCNC allow to set limits on tan β and the charged scalar mass also in the case of Higgs-alignment, that is weakly constrained by many other Higgs physics observables. Each possible assignment of the Z 2 charges results in a different type of 2HDM. These types are listed in Tab. 1.

Decay Widths and Branching Ratios
In this section we compare the branching ratios (BR) of the four neutral spin-0 states in the 2HDM+PS and 2HDM+S for the later discussion. We give analytic expressions of the dominant decay widths for those and the charged scalars in App. A. For all plots and interpretations we used the parameter values from Eq. (4.4) and fixed M A = 500 GeV and tan β = 1 such that the findings are applicable to all types of 2HDMs (besides the Inert one).
As mentioned, in the alignment limit the couplings of h to the SM fields substantially coincide with the ones expected for the SM Higgs boson. However, its total width can deviate with respect to the SM prediction, because of the eventual presence of additional decay channels. The most relevant, if kinematically allowed, is the one into a pair of a or S 2 states, respectively. Even for mediator masses above half of the SM Higgs mass, the three body decay aχχ, or S 2 χχ respectively, can be sizeable and even dominate for small enough DM masses, see Fig. 1 (for brevity only the two dominating SM decays are shown).  Figure   The width of the additional light (pseudo)scalar shows similar behavior in both models. It is dominated by the χχ channel, if kinematically accessible, even if the decay to tt is allowed, see Fig. 2. The other decay channels arise from mixing with the corresponding doublet state and are therefore suppressed by sin 2 θ.
The decay channels of the heavy scalars H/S 1 , see Fig. 3, and pseudoscalars A, see Fig. 4, are exchanged in the 2HDM+PS and S. All four are dominated by the decay to top quarks. In addition the heavy pseudoscalar A in the 2HDM+S and the heavy scalar H in the 2HDM+PS decay to ZS 2 , or Za respectively, enabling the resonant production of a mono-Z final state discussed below. It can be seen that, in the 2HDM+PS, BR(H → Za) is bigger by roughly a factor of two than BR(A → ZS 2 ) in the 2HDM+S. This is due to the smaller decay width of scalars to quarks compared to pseudoscalars, and therefore a smaller total width, resulting in a larger BR for the 2HDM+S than the 2HDM+PS to the mono-Z final state. In the 2HDM+S the heavy scalar S 1 has a BR of O(10%) to χχ (from mixing with S 2 ) and hS 2 , the latter declining with M S 2 . In the 2HDM+PS the analogue is the heavy pseudoscalar A, where BR(A → ah) and BR(A → χχ) are of

Overview of Model Constraints
In this section we give a brief overview on various constraints on the DM properties of the 2HDM+S/PS mainly from non-collider experiments and on 2HDMs in general as they apply to the studied case and reduce the parameter space.

Direct Detection
Direct detection (DD) phenomenology and their detection prospects are rather different for the two models. In the case of 2HDM+S model the DM is coupled to scalar mediators. It consequently generates Spin Independent (SI) interactions due to t-channel exchange of the scalar mediators at tree level, described by the following operator [5,11]: where: with the coefficients q as given in Tab. 1 and the coefficients f can be found in [25]. Such interaction is well within the reach of current detectors for , apart from the regions where there are negative interference effects [5]. In the following sections on collider searches we focus on light DM where DD is insensitive, therefore acting as a complementary search.
In the case of the 2HDM+PS model, the DM-nucleon interaction, from tree-level exchange of pseudoscalar mediators is described by the following operator [25]: where the effective coefficientc N is given, for example, in [26,27]. This operator, in the nonrelativistic limit, reduces to: where s χ , s N are the DM and nucleon spins, respectively, while q is the momentum transfer. The corresponding cross-section is hence suppressed by the fourth power of the momentum transfer and experimentally relevant only for very light masses of the mediators [13,26,27]. Spin independent DM nucleon interactions, within the reach of, especially, next generation detectors like XENONnT and DARWIN, emerge however at the one loop level [13,14,[28][29][30] (see also [12,31]).

Indirect Detection and Relic Density
Adopting the WIMP paradigm, the DM relic density, as well as eventual Indirect Detection (ID) signals, rely on DM annihilation processes. The case of the 2HDM+S model has been extensively studied e.g. in [11]. The possible annihilation channels include SM fermion-antifermion pairs as well as two scalars and scalar-gauge bosons final states. The full list of annihilation channels is ff , where f denotes all SM fermions with m f < m χ . All these annihilation channels are characterized by p-wave suppressed annihilation cross-sections. While it is possible to achieve the correct relic density by suitably accommodating the model parameters, we do not expect signals measurable through ID experiments.
The phenomenology related to ID and relic density for the 2HDM+PS model has been extensively reviewed in [32] (see also [12,14,33]). The DM relic density is determined by annihilation processes into ff , XV , with X = h, H, H ± , a, A and V = W ± Z, and the aa, aA, AA final states. The first two types of annihilation channels have s-wave dominated cross-sections, hence also relevant for ID (see below). The annihilation cross-sections into pseudoscalars are, instead, p-wave suppressed. The aa channel can be, nevertheless, especially for very light a, relevant for the relic density.
Concerning ID, possible signals are mostly accounted by the bb, tt and ha channels, the latter giving a four fermion, typically bbbb, signature. For heavy DM, the hA, HZ and H ± W ∓ final states might play a role as well. The aa, aA and AA final states contribute to ID to a negligible extent, because the corresponding cross-sections are p-wave suppressed. As discussed in [32], FERMI-LAT bounds on the bb and tt channels [34] can probe DM masses from 190 GeV up to around 400 GeV, therefore offering a complementary approach to the collider searches. No dedicated studies exist, instead, for annihilation into a gauge and a Higgs boson.

Flavour Constraints
As already mentioned, we are considering specific realizations of the 2HDM in which flavor violation processes are forbidden at tree level through suitable combinations of the couplings of the SM fermions with the different Higgs doublets. FCNC processes can be nevertheless induced at loop level. A comprehensive discussion of the possible bounds has been performed in [35]. Among them, the most relevant ones come from b → s transitions, in particular the B → X s γ process, whose rate is mostly sensitive to tan β and M H ± . b → s transitions mostly constrain 2HDM realizations featuring tan β enhanced interactions of the BSM scalars with d-type quarks, namely the type II and type Y models. For them we have a lower bound M H ± > 570 GeV weakly sensitive to the value of tan β. The bound on M H ± can be translated into bounds for the masses of the other extra scalars in view of the relations imposed by bounds from EWPT and the scalar potential, as illustrated in the subsections below (see also e.g. [15,33]). Additional constraints, from B s → µ + µ − and B → Kµ + µ − processes, affect the type II model for moderate/high values of tan β. Models of type I and type X , having tan β suppressed interactions with the SM quarks, are constrained only in a small region of the parameter space for tan β 2.

Electroweak Precision Constraints
It can be shown that the scalar potential for the 2HDM+S in Eq. (2.1) and (2.18) breaks the custodial symmetry [36][37][38][39][40][41] (in particular the λ 4 , λ 5 and λ hHS terms) and leads to contributions to EW precision observables, unless M A = M H + . The most relevant observable is the ρ parameter, which receives a contribution of One can verify that ∆ρ = 0 for M A = M H + . Similarly, for the 2HDM+PS potential in Eq. (2.7) and (2.12), This, together with the fact that to have resonant enhancement of mono-Z and mono-h one of the (pseudo)scalars needs to be below and the other above the top threshold, is the reason that motivates our choice of concentrating on M S 1 > M S 2 , respectively M A > M a . For more details on EW precision constraints for these models, see [10,11].

Perturbativity and Unitarity Constraints
Perturbativity and unitarity constraints for the 2HDM+S model have been studied in detail in [5] (see also [42]). By requiring that the couplings of the scalar potential are perturbative and that the amplitude of scattering processes of the type φ i φ j → φ k φ l involving scalar states preserves unitarity, we obtain the following constraints: (3.8) These bounds on the quartic couplings can be also re-expressed into upper limits on the mass splittings between the scalar Higgs eigenstates (see, for example, [5]). Together with Eq. (3.8) one should also take into account the following constraints from the requirement that the scalar potential is bounded from below: Unitarity constraints for the 2HDM+PS model have been considered, instead, in [7]. These rely on the amplitudes of the processes aa, aA, AA → W W and can be expressed through the following condition: where Requiring perturbativity on top of unitarity strengthens the limit down to (3.12)

Collider Searches
In this section, we will focus on the collider searches for tt + / E T and mono-jet, whereas tt, mono-Z and mono-h will be discussed in detail in Sec. 4, since they lead to the most stringent limits for the 2HDM+S/PS.

tt + / E T
New spin-0 mediators with large invisible widths can be searched in the tt + / E T (and bb + / E T ) channels. The most recent experimental searches have been reported in [43,44]. Their results have been interpreted in the context of the so-called DM-simplified model (DMF). They can be, however, straightforwardly applied to the 2HDM+PS model, as long as M a M A , by applying the simple scaling relation (a similar relation applies for bb + / E T with the replacement of tan β according to Tab. 1) (3.13) As discussed in [32], the limits of [43,44] can be applied, in analogous manner as illustrated above, also to the 2HDM+S, as long as there is substantial mass splitting between the BSM spin-0 scalars. In case that the new spin-0 scalars have comparable masses, the / E T spectrum features distortions with respect to the DMF case. More refined procedures to map the experimental limits on the models under study should then be applied. The case of the 2HDM+PS model is illustrated in [32].
For the extended 2HDMs we are considering, the tt + / E T exclusions are shown to be subdominant [45]. Therefore, we don't derive explicit bounds from those searches.

Mono-jet
Similarly to what occurs in heavy flavors + missing energy searches, observational constraints are typically interpreted in the context of simplified models. In the limits in which the 2HDM-like neutral bosons are decoupled, in mass, with respect to the singlet-like states, the kinematic distribution of the mono-jet events are substantially the same for the DMF and 2HDM+S/PS models. Experimental limits can be then applied to our setups just by using analogous scaling relations to Eq. (3.13).
As mono-jet events provide the strongest bounds among the initial state radiation signatures [46], we explicitly checked promising points with the help of the CheckMATE [47] implementation of the latest ATLAS search [48] and found no excluded points for tan β = 1. As shown in [10], mono-jet limits only arise for tan β < 1, however, this region is already excluded by other constraints and therefore we do not investigate this further.

Comparison of 2HDM + Scalar/Pseudoscalar and their LHC Signatures
After the general overview of constraints in the last section, we will now turn to the 2HDM+S/PS specific LHC signatures and limits we obtained for the leftover parameter space by dedicated collider simulations.

General Aspects
Before describing the simulations and discussing our parameter space in detail, let us first have a brief look at resonant production processes, which can be described analytically and will become handy for the interpretation of the simulation results. Resonant production The production cross section for a spin-0 scalar (or pseudoscalar) resonance S, with mass M and total decay width Γ tot , subsequently decay into the final state X, can be schematically written, in the narrow width approximation, as [49] σ The sum over i runs over the possible partonic initial states, for example quark or gluon pairs, C i are weight factors that account for the protons PDFs and colour factors, and s = (13 TeV) 2 is the center of mass energy squared. The values of the C i are given by: with q(x) (q(x)) and g(x) being, respectively, the parton distribution functions (PDFs) of quarks (anti-quarks) and gluons inside the proton and x is the conventional Bjorken scaling variable. In general, the two dominant production mechanisms are gluon fusion and bb initial states. We show the production cross sections of the heavy (pseudo)scalar, as a function of the resonance mass for both production modes in Fig. 5. The left panel shows the production cross section for gluon fusion, the right panel for bb initial states. For tan β = 1, the production cross section from gluon fusion results to be about 100 times larger than from bb initial states. This dominance of gluon fusion production depends on the 2HDM Yukawa type in combination with the value of tan β, but holds for all 2HDM Yukawa types in the whole parameter space we take into consideration (0.3 ≤ tan β ≤ 3), so in our simulations we can safely neglect the bb production mode.
From Fig. 5, we can also see that the production of a pseudoscalar particle always results in a larger cross section. This is because the effective coupling of the pseudoscalar mediator is larger, from same assignation of the Yukawa couplings, with respect to the scalar mediator, see e.g. [33,50].
Parameter overview For a better overview we briefly summarize our choice of parameter values, which are motivated by various constraints discussed above and used for the following Monte Carlo simulations and other examples throughout the paper: Signal Generation We simulate events for the signal processes with MadGraph5 aMC@NLO [51][52][53][54][55] at next-to-leading order (NLO) in QCD using the 263000 PDF set (NNPDF3.0) [56] provided through LHAPDF6 [57]. For the parton-showering we use the MadGraph built-in Pythia 8.2 [58]. A fast detector simulation is done with Delphes 3.4.2 [59] using the provided CMS detector card. The final cuts are implemented in MadAnalysis 5 [60,61]. The correct implementation of the program chain and analysis is checked by reproducing the mono-Z and mono-h exclusions for the 2HDM+PS presented in [32,45].
In the following, we will look at the different channels that provide the strongest limits for the 2HDM+S/PS model.

tt Resonances
For masses above the top-threshold, M H/S 1 , A > 2 m t , the additional heavy Higgses dominantly decay into a pair of top quarks, cf. Sec. 2.4. Searches for tt resonances are therefore a powerful tool to test 2HDMs in general and 2HDM+S/PS in particular.
One aspect that complicates the analysis, is that the signal processes interfere non-trivially with the SM background, as pointed out and taken into account by the experimental analysis in [62]. There the ATLAS results are interpreted in a 2HDM of type II, for the case in which the individual contribution of the two mediators to the signal can be singled out as well as for the mass-degenerate case where both mediators contribute. As the later case gives significantly stronger constraints, deriving exclusions for a single mediator serves as a conservative estimate.
The most recent bounds for √ s = 13 TeV and an integrated luminosity of 35.9 fb −1 are provided by CMS in [63]. The limits are presented in terms of simplified models with either a scalar or a pseudoscalar and Yukawa-like couplings to tops and in the hMSSM. Here, a 1.9 σ "signal-like deviation" is observed that would fit to a pseudoscalar with a mass of around 400 GeV, therefore the limits in that mass range are not significantly stronger than the ones from the previous ATLAS analysis with √ s = 8 TeV.
The additional light state a/S 2 can have a non-trivial impact on the limit due to interference effects, cf. Sec. 7.1 in [32]. However, this effect is expected to be small for our choice of small mixing sin θ = 0.3. A detailed analysis of the impact of interference and combining the limits for the two heavy states is beyond the scope of this article. Instead, we recast the recent CMS limits [63] for single mediators to our parameter space, interpolating between the different total width to mass ratios given in the paper. As commented above in the context of [62], this can bee seen as a conservative estimate, since the limits get significantly stronger by taking the contribution of both mass-degenerate states H/S 1 and A into account.
The limits are shown in Fig. 6 in the M A,H/S 1 -tan β-plane for the 2HDM+PS (top) and 2HDM+S (bottom). We choose the setting in Eq. (4.4) and M a/S 2 = 400 GeV. There exists a mild dependency on the mass of the light state due to changes of the total width which is shown in Fig. 13 and 14.
The decay widths of pseudoscalars to quarks and gluons (for M A above the top threshold) are bigger than the ones of scalars, cf. Sec. 2.4, App. A and [64]. Together with the effective coupling approximation for the dominant gluon fusion production, cf. Eq. (4.1) and [65], it can be understood that pseudoscalar resonances are expected to provide a stronger constraint than scalar ones. In the 2HDM+S the couplings of the heavy scalar S 1 are reduced due to the mixing with the singlet S 2 , while the pseudoscalar couplings are untouched. As a consequence the exclusion for S 1 reaches up to tan β = 1, while the one for A clearly exceeds tan β = 1. In the 2HDM+PS instead, the stronger constrained pseudoscalar mixes with the singlet a, weakening the exclusion and leading to overall weaker constraints from tt searches for this model.
In both models, for masses of the heavy Higgs around 500 GeV, tt resonance searches provide a strong lower limit on tan β and essentially exclude values below tan β 1.

Mono-Z
For a strong mono-Z signal, which is a general feature for 2HDMs extended by pseudoscalar mediated DM [66], the heavy neutral spin-0 particle of the doublet which doesn't mix with the singlet has to be produced, meaning the scalar H in the 2HDM+PS and the pseudoscalar A in the 2HDM+S. This resonantly produced state can decay to a Z boson and the light state a or S 2 , respectively, which further decay to χχ with a high branching ratio, cf. Sec. 2.4. The Feynman diagrams for the processes are shown in Fig. 7 for the dominant gluon fusion production.
Searches for the relatively clean final state, where the Z boson decays leptonically (meaning to electrons or muons), in association with / E T are performed by the ATLAS and CMS collaborations [67,68].
/ E T Spectrum The maximum value of the missing transverse energy / E T can be obtained from kinematics and is given by [10] / E where λ(m 1 , m 2 , m 3 ) is given in Eq. (A.17) and the first (second) subscript is used for the 2HDM+S (PS). The missing energy spectrum is given by  Two example / E T spectra after detector simulation for the 2HDM+PS in the e + e − + / E T final state are shown in Fig. 8 together with the predicted SM backgrounds and the observed number of events as provided by ATLAS [67]. Nearly identical / E T spectra exist for the µ + µ − + / E T final state. As described below both final states are combined to determine the exclusion limits. The signal features a peak a bit below / E T = M A /2 on top of the smoothly falling SM background.
Backgrounds For mono-Z searches the main irreducible background is ZZ production with one Z decaying to neutrinos. Another important background is W Z production, where one lepton from  the W -decay escapes detection or a τ decays hadronically, cf. Fig. 8. Minor backgrounds are Z+jets processes with poor / E T reconstruction and non-resonant production. In [67], the backgrounds are estimated from simulations and data-driven methods. For both the electron and muon final state the background uncertainty is dominated by systematic errors which are driven by the uncertainty on the Z+jets background.

Current Constraints
To determine model constraints we use the ATLAS results [67] as they are also used by the LHC-DMWG and easier to reproduce than the (slightly stronger) CMS results [68]. To be explicit, for our exclusion bounds, we use the expected number of background events b and the corresponding uncertainty σ b from [67], together with our simulated event numbers s for various parameter points and the sensitivity formula from [32,69]: . (4.7) The index i refers to the different bins and the values for Z i are added up quadratically to find the (square of the) overall sensitivity Z 2 , where one expects to exclude parameter points with Z > 2 at 95% confidence level. In this way we obtain the expected exclusion limits, however, as they are very similar to the observed limits, cf. [45], we can use the expected mono-Z limits to compare them to the observed tt and mono-h ones in Sec. 4.5.
The constraints from mono-Z searches are similar in shape and reach for both models, see Fig. 9. This can be understood with help of the approximation in Eq. (4.1) and its corresponding discussion: on the one hand, in the 2HDM+S a heavy pseudoscalar A needs to be produced, which happens with a production cross section that is bigger by roughly a factor of two than the one for the production of the scalar H in the 2HDM+PS, cf. Fig. 7 for the Feynman diagrams. On the other hand, the relevant branching ratios for the mono-Z final state is bigger in the 2HDM+PS by a factor of approximately two, compensating for the smaller production cross section. The difference in the branching ratios is related to the different decay width of scalars and pseudoscalars into top quarks, as mentioned in Sec. 2.4. The general features of the constraints for the mono-Z channel depicted in Fig. 9 can be readily understood by considering kinematics and the couplings. In the M a/S 2 -tan β-plot, the weakening of the exclusion limit for larger tan β values is due to the fact that the top-coupling scales like (tan β) −1 and the top-coupling is essential for the production of the intermediate heavy Higgs, cf. Fig. 7. For the M a/S 2 -M A -plot, there are three different features that can be understood: first, the "diagonal" lower bound to the exclusion region is due to the fact that for M A M a/S 2 + M Z , resonant production is not allowed with on-shell a/S 2 . Second, the upper bound of the exclusion limit stems from the heavy Higgs being the harder to produce the heavier it is and third, a heavier M a/S 2 leaves less energy available for the Z, so the Z production gets kinematically suppressed, thereby smoothing out the transition between the first two features.

Projected Sensitivity
With the detailed data of the experimental analysis at hand [67], we estimate the reach of the high luminosity phase of the LHC (HL-LHC). To do so, we assume integrated luminosities of 300 fb −1 and 3000 fb −1 and took a projected reduction of systematic uncertainties by 50% into account called YR18 scenario [70]. The projections for the pseudoscalar model are shown in Fig. 10 and are nearly identical to the ones for the scalar model, as can be seen for the current bounds in Fig. 9. From  Fig. 10, one can see that the increased integrated luminosity will lead to a substantial increase in sensitivity, strengthening the 5 σ limits by roughly a factor of two in terms of the masses for 3000 fb −1 . In contrast to that, the influence of the improved systematic uncertainties will likely be small.

Mono-h
The most recent searches for mono-h with h → bb by ATLAS and CMS can be found in [71,72]. Furthermore, in [73] CMS performs a first search of mono-h with h → bb, γγ, τ + τ − , W + W − , ZZ, where the exclusions are dominated by the bb channel. For our analysis, we use the model-independent upper limits on the h + / E T cross section provided by ATLAS [71] and compare it to the one we find from our simulation, as also done by the LHC-DMWG.
Backgrounds The main backgrounds for mono-h with h → bb are top quark pair + single top as well as leptonically decaying vector bosons + jets processes [71]. In case of top quark pair + single top, a top decaying to a bottom can look like it came from a Higgs decay in combination with another b-tagged jet, while the missing energy is in neutrinos. Similarly, decaying vector bosons also produce / E T via decays to neutrinos or missed charged leptons and the jets can by chance look like a Higgs decay. To generate those backgrounds Monte-Carlo simulations are used in [71].

Current Constraints
The non-deviation of the observed events from the SM background is translated into (model-independent) upper limits on the h(bb) + / E T cross section, which we use in the following by comparing our simulated cross-section to this upper limit, as also done by [32]. In [71], only one / E T bin is used at a time to minimize the model dependency, which implies that the derived limits are conservative estimates. Furthermore, the dependency of the limit on the kinematics within one bin and on the acceptance and efficiency is calculated for several parameter points of a benchmark model and the least stringent limits are given -again, this leads to conservative estimates for the exclusion limits of the 2HDM+S/PS, which are shown in Fig. 12.
Similarly to the mono-Z exclusion plots (cf. Fig. 9), we can see that for larger tan β, M A and M a/S 2 the exclusion limits in each case weaken, as they also do in the case of M a/S 2 ∼ M A . The otherwise most prominent feature in the M a/S 2 -M A -plane (cf. left panel of Fig. 12), is the dip in the exclusion limits around M A ∼ 700 GeV. It originates from the single / E T bin used at a time to derive the upper limit on the cross section in the ATLAS search, or in other words that one cannot be model-independent and use the spectral shape information of the / E T -spectrum at the same time. Hence, splitting the events between different / E T bins leads to weaker limits. For M A > 700 GeV a significantly higher fraction of the signal events reaches / E T > 350 GeV (as can be seen from the simulated data) and therefore ending up in the stronger constrained bin reaching from / E T = 350 to 500 GeV. This (over-)compensates the decrease in production cross section by making the mediator heavier, leading to a comparably strong or even stronger limit. The effect of more events ending up in the bin of / E T = 350 to 500 GeV can also be understood in the light of the / E T spectrum discussion in case of the mono-Z (cf. Sec. 4.3 and Fig. 8). Since the / E T spectrum does not qualitatively change by replacing the Z with an h, for values of M A > 700 GeV the peak of the / E T spectrum starts to shift from the / E T = 200 to 350 GeV bin to the / E T = 350 to 500 GeV bin because the peak is located a bit below / E T = M A /2. In contrast to the mono-Z searches, the exclusions limits from mono-h differ significantly between the 2HDM+S and 2HDM+PS. The upper bound on the mass of the light pseudoscalar a is stronger by approximately a factor two in the 2HDM+PS than the one for the corresponding scalar S 2 in the 2HDM+S. This is due to the fact that, to resonantly produce mono-h events in the 2HDM+S, one needs to produce the heavy scalar S 1 , rather than the the heavy pseudoscalar A, and the production cross section for scalars is smaller than the one for pseudoscalars, cf. Eq. (4.1) and Fig. 5. The opposite happens for the 2HDM+PS, where again the roles of the heavy scalar and pseudoscalar are switched going from mono-Z to mono-h, see also the corresponding Feynmann diagrams in Fig. 7 and 11. Opposite to the mono-Z case, the branching ratios relevant for mono-h events are similar in both models, cf. Sec. 2.4: Therefore the 2HDM+PS is expected to have a higher mono-h cross section in the resonant region and as the kinematics do not change significantly, this results in stronger exclusion bounds.

Combined Constraints
Finally, we summarize our results by comparing all derived limits for the two models in the M a/S 2 -M A -plane and the M a/S 2 -tan β-plane in Fig. 13 and Fig. 14, respectively. In addition to the tt, mono-Z and mono-h constraints discussed above, we also include results from Higgs to invisible  for the observed (expected) upper limit at 95% confidence level. As a and S 2 dominantly decay to DM, the 3-body final state a/S 2 χχ is also invisible. Using the decay widths from Sec. 2.4 we find a lower limit for M a/S 2 of about 100 GeV for m χ = 10 GeV with a residual dependence on the mass of the heavy Higgs M H/S 1 = M A = M H ± . This is also the reason why we start the parameter scans at 100 GeV in terms of the light new Higgs mass M a/S 2 .
The differences between the pseudoscalar and scalar model in the single searches have been discussed in the corresponding last chapters, so here we focus on how the different constraints compare to each other within one model. Starting with the 2HDM+PS model in the M a -M A -plane with tan β = 1, cf. left panel of Fig. 13, we can see that the dominant limit comes from mono-Z searches, excluding light Higgs masses up to M a ∼ 320 GeV and heavy Higgs masses between 200 and 1000 GeV, while for large masses of the heavy Higgses and small masses of the light Higgs the mono-h limit leads to slightly stronger bounds. Light masses below ∼ 100 GeV for the new pseudoscalar state are excluded by Higgs to invisible searches, while tt resonance searches can exclude small parts of the parameter space nearly independent of the light new pseudoscalar mass around M H = M H ± = M A ∼ 400 and 500 GeV.  The same holds true for the 2HDM+S model, cf. right panel of Fig. 13, especially for the mono-Z and Higgs to invisible searches. In contrast to that, the mono-h limit for the scalar model is significantly weaker compared to the mono-Z limits and thus never the strongest limit, as explained in Sec. 4.4. Furthermore, the tt limits exclude a larger band of heavy Higgs masses from 400 to 500 GeV, however this is just a reflection of slightly stronger limits compared to the pseudoscalar case with the limits changing from just below 2 σ to just above 2 σ.
For the M a/S 2 -tan β-plane and at heavy Higgs masses of 500 GeV, cf. Fig. 14, we have again similar limits for both models, except for the mono-h limit being weaker in the scalar model. As in the M a -M A -plane, the dominant limit is given by mono-Z searches, specifically for values of tan β > 1. They begin to weaken for larger values of tan β of around three. While these limits, for 0.3 tan β 3 apply to 2HDMs of all Yukawa sector types, cf. Tab. 1, for tan β 3 the bb production mode starts to become relevant for type II and Y, leading to stronger limits for such Yukawa sectors, while for type I and X bb production never becomes relevant, and limits will continue to get weaker for larger tan β values. The tt resonance searches provide a lower limit on tan β of around 1, being slightly stronger in the 2HDM+S and nearly independent of the light mediator mass, therefore being the strongest limit for M a/S 2 > 270 GeV.
Finally, another interesting aspect which is accessible via the comparison plots is the question of how to distinguish the two models. Here, the weaker mono-h limits for the 2HDM+S model can come in handy. This discrepancy in terms of sensitivity to the mono-h channel could be exploited to distinguish between the 2HDM+PS and 2HDM+S, since the ratio of the signal strength in mono-Z and mono-h is characteristic for the model. So if signals are detected in both channels, their signal strength ratio could be used to discriminate between the the models.

Conclusions
In this paper, we have investigated the collider phenomenology of the 2HDM+S and 2HDM+PS, in particular focusing on the tt, mono-Z and mono-h signatures.
The 2HDM+S and 2HDM+PS are the new standard paradigm used by the CMS and ATLAS collaborations to interpret experimental results in the context of DM searches and Simplified Models. These models feature an extended scalar sector, with a second Higgs doublet, and an additional singlet and differ from the previous generation of Simplified Models and other DM models by tending to generate a rich collider phenomenology. In particular, it is possible to have resonant production of the heavy Higgs bosons that can decay to a SM boson and the additional singlet, generating resonantly-enhanced mono-Z and mono-h signatures.
In our analysis, we made use of a few assumptions that have been established by the experimental collaborations [8] to reduce the number of dimensions of the parameter space, and simplify the comparison between the two models and the different experimental signatures. These assumptions are essentially the mass degeneracy of the heavy scalars, motivated by the bounds on the EW precision observables, together with the Higgs alignment limit, motivated by Higgs measurements. The alignment limit simplifies the Higgs sector by turning one of the doublets into the SM Higgs doublet, so that most constraints from Higgs physics are avoided. Moreover, while we have described the possible Yukawa sectors of the 2HDM, our results are universal due to the range of parameters we considered.
We have reviewed all the principal constraints of the two models. Some of them, like Direct and Indirect Detection, are mostly complementary to collider searches, because they tend to require different mass spectra. Some other, like flavour and EW precision observables, perturbativity and unitarity constraints, instead are very relevant also for our signatures, and are used as a guidance to select the appropriate ranges for the various parameters.
Using the data provided by different LHC analyses, we derived limits for the 2HDM+S and 2HDM+PS for the tt, mono-Z and mono-h signatures and discussed how they compare to each other and between the two models. We found that the mono-Z limits are in general the most constraining ones, while the tt limits are nearly independent on the mass of the additional singlet and especially relevant for constraining tan β. A lower bound of ∼ 100 GeV for the additional singlet mass is given by Higgs to invisible searches, giving us a natural starting point for this parameter in our scans.
We found that in principle the two models could be discriminated at a collider by detecting both mono-Z and mono-h signatures, from the ratio of their strengths. Also, the absence or appearance of mono-jet signals would give further insights into the nature of the dark sector. However, depending on the DM mass, other probes, mostly astrophysical, would be powerful tools to help discriminating between the two models. A Direct Detection signal for DM is several orders of magnitude stronger for the 2HDM+S, while Indirect Detection signals would be several orders of magnitude stronger for the 2HDM+PS. So detecting one of these two astrophysical signatures would also give a clear indication towards the nature of the mediators to the dark sector and discriminate between 2HDM+S and 2HDM+PS. In this case, collider studies could help to further investigate the inner workings of the dark sector.

A Formulae for the Decay Widths
In this appendix we give analytic expressions for the dominant branching ratios of the spin-0 states in the 2HDM+S/PS as partially shown in Sec. 2.4.
We focus on the mass hierarchy M A , M H/S 1 , M H ± > M a/S 2 , M h and tan β = O(1). For the value of f , denoting the ratio between the Yukawa coupling of the fermion f in the different types of 2HDMs and the SM value y f = √ 2m f /v, as given in Tab. 1. In the following we use the abbreviation τ i,j = 4M 2 i /M 2 j .

A.1 Scalar Model
Higgs Boson h As mentioned above, in the decoupling limit the couplings of h with the SM states substantially coincide with the ones for the SM Higgs boson. However its total width can deviate, with respect to the SM prediction, because of the eventual presence of additional decay channels. The most relevant, if kinematically allowed, is the one into a pair of S 2 states. As the Higgs width is small, also three-body decays to S 2 χχ can be relevant and the additional widths are given by Light Scalar S 2 The light scalar S 2 mostly decays into gg, ff and χχ (direct couplings with gauge boson are forbidden in the alignment limit), depending on the mass. We quote below the corresponding decay widths and the loop-induced one into gluons which is useful for the interpretation of the collider studies: Γ(S 2 → gg) = α 2 s 16π 3 M S 2 sin 2 θ q 2 q y 2 q F S (τ q,S 2 ) , (A.6) Heavy Scalar S 1 Besides an additional cos 2 θ due to the mixing with the additional singlet the couplings of the heavy scalar to SM fields remain similar to the ones known from 2HDMs. Additional decay channels are χχ, suppressed by sin 2 θ,S 2 S 2 , which is very small for our parameter choice and