The THDMa revisited

The THDMa is a new physics model that extends the scalar sector of the Standard Model by an additional doublet as well as a pseudoscalar singlet and allows for mixing between all possible scalar states. In the gauge-eigenbasis, the additional pseudoscalar serves as a portal to the dark sector, with a priori any dark matter spins states. The option where dark matter is fermionic is currently one of the standard benchmarks for the experimental collaborations, and several searches at the LHC constrain the corresponding parameter space. However, most current studies constrain regions in parameter space by setting all but 2 of the 12 free parameters to fixed values. In this work, we perform a generic scan on this model, allowing all parameters to float. We apply all current theoretical and experimental constraints, including bounds from current searches, recent results from B-physics, in particular B_s ->X_s gamma, as well as bounds from astroparticle physics. We identify regions in the parameter space which are still allowed after these have been applied and which might be interesting for an investigation at current and future collider machines.


I. INTRODUCTION
After the discovery of a scalar boson what complies with predictions of the Standard Model (SM) Higgs boson by the LHC experiments, the quest of additional particles that stem from possible further extensions of the SM scalar sector is an important task for the experimental collaborations. Furthermore, astroparticle observations give evidence for the existance of dark matter. In a standard cosmological description, the particle content of the SM alone does not suffice to explain these observations. Therefore, many beyond the SM (BSM) extensions in addition provide a dark matter candidate. In this work, we concentrate on the model which has been proposed in [1][2][3][4][5]. It extends the scalar sector of the Standard Model by an additional doublet as well as a pseudoscalar singlet, which in turn couples to a fermionic dark matter candidate. The model therefore extends the scalar sector by 5 additional particles, which we label H, A, a, H ± in the mass eigenbasis, as well as a fermionic dark matter candidate χ. In the form discussed here, the model contains 14 free parameters after electroweak symmetry breaking, out of which 2 are fixed by the measurement of the 125 GeV scalar as well as electroweak precision measurements. Furthermore, the models parameter space is subject to a large number of theoretical and experimental constraints, which we will discuss in detail below. Reference [6] contains a detailed discussion as well as benchmark recommendations for the LHC experimental collaborations; see also refs [7][8][9][10][11][12][13][14][15] for more recent work on this model. Finally, refs [16][17][18][19][20][21][22][23][24][25][26][27][28][29][30][31][32][33][34] report on experimental constraints from current collider searches on the parameter space in specific parameter regions 1 . This manuscript is organized as follows. In section II, we briefly review the setup of the model, including the introduction of conventions used in this work. In section III, we discuss current theoretical and experimental bounds that need to be imposed. Section IV introduces the setup of our scan and identifies regions in parameter space that are still allowed after all constraints are taken into account. We present a first glance at possible processes at e + e − colliders in section V, and conclude in section VI.

II. THE MODEL
The model discussed in this work has been introduced in references [1][2][3][4][5], and we refer the reader to these works for a detailed discussion of the model setup. We here just list generic features for brevity. We largely follow the nomenclature of [6].
The field content of the THDMa in the gauge eigenbasis consists of two scalar fields H 1,2 which transform as doublets under the SU(2) × U(1) gauge group, and additional pseudoscalar P transforming as a singlet, as well as a dark matter candidate χ which we choose to be fermionic. The two Higgs Doublet Model (THDM) part of the potential is given by Fields are decomposed according to (see also e.g. [39]) where The scalar potential is Finally, the coupling between the visible and the dark sector is mediated via the interaction L χ = −i y χ Pχγ 5 χ.
Couplings of the scalar sector to the fermionic sector are equivalent to a THDM type II (see e.g. [39]) and follow from where Y i u,d denote the Yukawa matrices, Q and L are left-handed quark and lepton doublets, u R , d R , ℓ R label right-handed uptype, downtype, and leptonic singlets, andH = ǫ H * i , where ǫ = i σ 2 denotes the two-dimensional antisymmetric tensor. We impose an additional Z 2 symmetry on the model, under which the doublets transform as H 1 → H 1 , H 2 → −H 2 , in order to avoid contributions from flavour changing neutral currents. We furthermore set the transformation properties of the additional fields as P → P, χ → −χ. In this work, we concentrate on the case where Y u 1 = Y d 2 = Y ℓ 2 = 0, which corresponds to a type II classification of Yukawa couplings in the THDM notation and can be achieved by requiring u R → −u R under the said symmetry. Note the terms ∼ µ 3 , b P induce a soft Z 2 breaking.
After electroweak symmetry breaking, the model is characterized by in total 14 free parameters. The above potential induces standard mixing for the THDM part of the potential, which we characterize by cos (β − α) , tan β using standard notation. Furthermore, V P introduces a mixing between the pseudoscalar part of the THDM and the new pseudoscalar P , with physical states given by where we use the notation c x ≡ cos x, s x ≡ sin x here and in the following and A 0 denotes the mass-eigenstate of the pseudoscalar from V THDM prior to mixing. We choose as free parameters v, m h , m A , m H , m H ± , m a , m χ , cos (β − α) , tan β, sin θ, y χ , λ 3 , λ P 1 , λ P 2 (6) following the suggestions in [6]. Relations between these and other quantities in the Lagrangian can be found in appendix A.
For the parameters in eqn (6), v and m h are fixed to be ∼ 246 GeV and ∼ 125 GeV, respectively 2 . We are therefore left with 12 parameters which can a priori vary, but are obviously subject to theoretical and experimental constraints. Previous studies of this model usually chose to consider scenarios where all but a few are fixed, in order to facilitate the phenomenological exploration and to provide consistent benchmark scenarios for comparison. In this work, we allow all parameters to float and determine regions in the models parameter space that are highly populated 3 .

III. THEORETICAL AND EXPERIMENTAL CONSTRAINTS
The models parameter space is subject to a list of theoretical and experimental constraints, which we list below.
A. Perturbativity, perturbative unitarity and positivity of the potential In order to guarantee positivity of the potential, several relations between the potential parameters need to be obeyed (see e.g. [6,12,14]). In the notation introduced here, they are given by λ 1,2,P 1,P 2 ≥ 0, λ 3 ≥ −2 λ 1 λ 2 , λ 3 + λ 4 − 2 |λ 5 | ≥ −2 λ 1 λ 2 . Furthermore, we require the self-couplings in the potential to obey an upper limit, which we chose as 4 π. We therefore have Finally, bounds from perturbative unitarity are implemented using the perturbative unitarity bounds derived in [12]; we refer the reader to that work for more details as well as explicit analytic expressions. We furthermore check that eqn (3.11) of [14](taken from [3]) is fulfilled. 2 In principle one can also consider the case that m H ∼ 125 GeV, with m h being a lighter resonance. We discard this scenario in the following. 3 Fine-tuned regions might exist within this model which are missed in a generic scan setup. B. Constraints from electroweak precision observables, flavour constraints, and total widths These bounds are tested via an implementation in SPheno [40,41], using an interface to Sarah [42][43][44][45][46]. Constraints from measurements of electroweak precision observables are taken into account using the oblique parameters S, T, U [47][48][49] as a parametrization. We take the latest fit results from the Gfitter collaboration [50][51][52], with central values given by S = 0.04 ± 0.11, T = 0.09 ± 0.14, U = −0.02 ± 0.11 and the correlation matrix as given in the above references. We demand ∆ χ 2 ≤ 8.02489, which corresponds to a 2 σ region around the central value determined by the SM decoupling.
Two Higgs doublet models and their extensions are also subject to strong bounds from flavour constraints, cf. e.g. [52,53]. Especially the m H ± , tan β plane is strongly constrained by precision flavour observabes as e.g. B → X s γ, B s → µ + µ − , and ∆M s . For these variables, we consider the following allowed regions: The most recent results on the theoretical calculation of this variable, including higher order calculations up to NNLO in QCD, have been presented in [54], superseding previous results [55]. In particular, the charged scalar mass is pushed to even higher masses. We implement the lower bound in this variable using a two-dimensional fit function [56] that reflects the bounds derived in [54].
The ATLAS, CMS and LHCb combination using 2011-2016 LHC data for the branching ratio B s → µ + µ − is given by [57] In principle, theoretical predictions depend on all novel scalar masses (see e.g. [58]).
In order to account for the missing higher order corrections, we apply a multiplicative correction factor to the Spheno result.
Assuming again a similar error in the theory calculation, and taking into account the additional higher order shift to the Spheno result, we then have where we allow for a 3σ deviation.
We use the following derivation for the SM prediction for this value: We start with eqn. (31) of [61] which shows the dependency of this variable on external, experimentally determined parameters. The above expression was updated to [62,63] For input values [64][65][66][67]  where we furthermore assumed an intrinsic error of 0.2 on the central value [63], which mainly stems from uncertainties of the top mass 4 .
Assuming the new physics contributions to be multiplicative, while keeping the error additive 5 as an allowed 2 σ range for this value. 4 Neglecting this leads to a slightly smaller error of ± 1.03 ps −1 . 5 TR wants to thank M. Mikolaj for useful discussions regarding this point.
Finally, although not strictly required, we impose an upper bound on the width over mass ratio on all new scalars. For the 125 GeV candidate, we require that [68] Γ h,125 ≤ 9 MeV. (9) For all other scalars, we impose Γ/M ≤ 0.5.
We consider this a very lenient bound in order to allow a detailed investigation of the parameter space; in practise, ratios larger than 10 − 20% already strongly question the applicability of perturbative treatment for such states.

C. Collider bounds
Agreement with current measurements of the 125 GeV signal strength as well as agreement with null-results from a large number of searches have been tested by interfacing the model with HiggsBounds [69][70][71][72][73][74][75] and HiggsSignals [74, [76][77][78][79]. Furthermore, the model is additionally constrained by dedicated searches by the LHC experiments; we have implemented the corresponding bounds in a pseudo-approximate approach and leave a detailed recast analysis for future work. All production cross sections have been produced using Madgraph5 [80], with the UFO model provided in [4,81].
• Bounds from ℓ + ℓ − + MET searches We consider the experimental bounds presented in [22], which corresponds to a CMS search for this channel mediated via Za production using full Run 2 luminosity. For the parameters specified in figure 8 in that work, we generate parton-level events. In order to marginally mimick the experimental cuts, we have imposed and simulating all parameter points using the cuts specified above. However, this is especially computationally intensive. Therefore, it is useful to consider the main channel leading to the above signature, which is    with subsequent decays a → χχ, Z → ℓ + ℓ − . The value from factorized onshell production times decay branching ratios into the Z χχ final state for the above benchmark points is given in table II.
We therefore could use as a naive condition that A more thorough inspection shows that the above bound has to be modified to σ fac Zχ χ ≤ 0.009085 pb in order to capture all points from our scan which obey eqn (11).
• Bounds from h + MET searches σ res/ unres refer to the resolved/ unresolved region as discussed in the text.
A priori, the dominant production times decay stems from where additional contributions can come from i.e., non-resonant production of A(a) and subsequent decays.
From table IV, we see that we can roughly exclude points for which the conditions are not fulfilled 6 .
For the two points giving the bounds in eqn (15), we have in a factorized approach A more detailed investigation show that for the first parameter point, (13) provides around 62% of the cross-section in the resolved region, while for the second, the processes in (14) mainly contribute. In the following, we require that For our sample points, we have explicitely checked that this value cuts out all parts of parameter space where σ h+ / E ⊥ ≥ 0.002818 pb.
• Bounds from H +t b, H + → tb searches Searches for the THDM sector of the model can be constrained by the above process, which has been presented e.g. by the ATLAS collaboration [25] using an integrated luminosity of 139 fb −1 , which we compare to the upper bound on the cross section as available from Hepdata [82]. Even a comparison with the production cross section for the process prior to the charged scalars decay shows that this search is not yet sensitive in the parameter space investigated here. This is mainly caused by the lower limit on tan β from B → X s γ which preselects point with higher tan β values.
• Bounds from W t + MET searches The ATLAS collaboration has presented results for this model in the above channel using full Run2 data [23]. They give results for bounds derived from single and double top production, divided into several search regions. We here implemented their results into grid-like upper bounds on the total production cross section from single-top like processes, with values taken from figure 28 of the auxiliary material available from [83]. In the mass region considered here, the smallest upper bound is around 32 fb 7 .
• Bounds from tt/bb + MET final states The tt final states have been investigated in [26, 28] using full Run 2 data. We note that the simulation for the signal a different model file was used, and cross sections were rescaled to NLO predictions [84]. In the comparison here, we stick to leading order, as QCD higher order corrections are universal for both models, differing only in the electroweak sector, as long a dominant production is assumed via an s-channel mediator. In table V, we give our predictions for the theory cross section predictions shown in figure 15 of that reference. Note that for generation we used sin θ = 1/ √ 2, as otherwise no significant coupling to both the SM and the dark sector can be achieved. The values in table V have been rescaled accordingly 8 . We use / E ⊥ > 230 GeV. 7 Note that including tt results will improve these bounds. However, here theoretical cross sections which were e.g. used in figure 31 of the above auxiliary material includes detector simulation and cuts; we thank the authors of the above reference for useful discussions regarding this point. Including these would require a recast-type investigation, which is beyond the scope of the current work. 8 The coupling of a to tt is rescaled by sin θ, while the coupling to χχ is rescaled by cos θ. We therefore, in comparison to the calculating in [26] for g = 1 obtain an additional factor 1 4 in the results.
corresponds to a rough estimate of the multiplicative factor from figure 15 of that reference. The last column gives the actual exclusion cross section.
The factor corresponds to a rough estimate of the multiplicative factor from figure 15 of that reference. The last column gives the actual exclusion cross section.
Following the same logic as before, we see that we can naively exclude points where  This leads to bb + / E ⊥ : σ ≤ 0.01420 pb The same channels have been investigated in [16] using an integrated luminosity of 36 fb −1 . In tables VIII and IX, we again show the production cross section for selected points from that reference, as well as resulting upper bounds on these cross sections.
For tt + / E ⊥ , we use / E ⊥ 300 GeV; for bb + / E ⊥ , we take / E ⊥ 180 GeV, p ⊥,b 20 GeV. As in the comparison with results from [26], we used a slightly different way to generate the event samples for consistency reasons.
Following the same logic as before, we now require that σ tt+ / E ⊥ ≤ 0.000703 pb; σ bb+ / E ⊥ ≤ 0.578 pb with the cuts as specified above 9 .

D. Dark matter constraints
The THDMa renders a dark matter candidate and is therefore subject to constraints from astrophysical measurements, such as relic density and direct detection. We have made use of MadDM [88][89][90] for the calculation of relic density. For direct detection, we implemented the analytic expressions presented in [1] 10 . We compare these values to limits from the Planck collaboration [93], and require that Ω h 2 ≤ 0.1224 (16) 10 A more detailed analysis, as e.g. presented in [9,12,91,92], is beyond the scope of the current work. which corresponds to a 2 σ limit. Direct detection bounds are compared to maximal cross section values σ Xenon1T max (m χ ) using XENON1T result [94], which we implemented in terms of an approximation function 11 . Note we rescale bounds from direct detection according to where m χ,i , Ω i refer to the dark matter and relic density of the specific parameter point i tested here.

A. Scan setup
The above constraints are implemented using an interplay of private and publicly available codes in several steps. We discuss these in the following. This also means that, as bounds are implemented successively, subsequent constraints are imposed on the subset of parameter points which have passed previous constraints.
2. Points passing these bounds are then passed on to Spheno, which is used for the calculation of limits from electroweak oblique parameters, B s → µ + µ − , ∆ M s , as well as the total widths for all particles. These values are confronted with the experimental constraints discussed in section III B.
3. The remaining parameter points are now passed through HiggsBounds and HiggsSignals, for comparison with null-results from experimental searches as well as Higgs signal strength measurements which are contained in these tools. 4. We then calculate the relic density as well as direct detection constraints, which are discussed in section III D. 5. Parameter points which have passed all remaining constraints are confronted with the collider bounds presented in section III C.
Our initial scan ranges are determined by a number of prescans to determine regions of parameter space that are highly populated, while still rendering relavitely low masses up to 1 TeV. In particular, we set The values of m h = 125 GeV and v = 246 GeV are set according to measurements of the Higgs boson mass as well as precision measurements. Note the above choices do not imply that values outside the above regions are stricly forbidden (apart from bounds we will discuss explicitely below), but these were found to render a relatively large acceptance rate for the scan discussed above. We furthermore apply the following relations between masses which were superimposed on the original scan ranges 12 . Note that a priori also points outside the above regions were allowed; these were chosen to optimize the selections process imposed by cuts. Note that the fit for B → X s γ implies a lower bound on m H ± of ∼ 800 GeV 13 . 12 We obviously adjusted ranges such that all m i ≥ 0 GeV. 13 We are aware that using these observables in a fit rather than as hard bounds might weaken some of the constraints discussed here, as e.g. discussed in [96] for a normal THDM. This is however beyond the scope of the current work. for m H ± ≤ 850 GeV is set by the bound on B → X s γ as discussed in section III B.

B. Scan results
In the following, we discuss the resulting constraints on the parameter space of the THDMa. Note that not all bounds discussed above lead to a direct limit in a two-dimensional parameter plane. In case the effects can be displayed in such a simple way, we discuss them explicitely.
We generate ∼ 30000 parameter points in the first step of the above scan. Application of bounds in step 2 includes e.g. constraints from flavour-physics in the m H ± , tan β plane. Results from flavour constraints in the m H ± , tan β plane are shown in figure 1, where similar results have been obtained e.g. in [52,53,97] 14 14 Note the newest result from the LHC experimental combination for B s → µ + µ − [57] increases the discrepancy between experimental result and theory prediction, leading to a larger exclusion region. In the mass range considered here, we found the limit can be approximated by requiring tan β ≥ 2.15302 − 0.000930233 m H ± GeV . Furthermore, in [52], the ∆M s bound is not shown in the figure for flavour constraints. However, similar results are obtained as well. TR wants to thank the GFitter authors for clarifying this point.
planes from oblique parameters. We see that regions where both displayed mass differences are large are excluded by the oblique parameters.
After the bounds imposed by theory, the only constraints we find from widths are points which violate (9). This basically rules out all points with m a ≤ 62.5 GeV, as the h → a a decay becomes dominant and can lead to large partial widths 15 . All other mass ratios are Γ/m 28% for all points that fulfill theory constraints, and maximally around 24% for Γ A /m A after width bounds are included. As an example, the point with the largest allowed ratio features a dominant branching ratio for A → ha, mediated via relatively large tan β ∼ 8.1, λ P 2 ∼ 11.  (m a , | sin θ|), where the oblique parameters also set a clear limit, see figure 3. We also investigate how our results compare to the ones obtained in the THDM limit, where a is decoupled. This can be achieved setting sin θ = λ P 1 = λ P 2 = 0. We then compare the limits in the (m H ± − m A , m H − m A ) plane, cf. fig. 4. We see that the admixture of a leads to much weaker constraints on m H ± −m A . For the remaining constraints discussed in section III B, the identification of additional regions in two-dimensional parameter planes that are primarily affected is not straightforward. from oblique parameters; for the latter, sin θ = λ P 1 = λ P 2 = 0. The admixture of a releases the bounds, as expected.
We now turn to the effects of applying collider constraints via HiggsBounds/ HiggsSignals. As expected in THDMs and its extensions, only narrow stripes around cos (β − α) = 0 remain in agreement with current signal strength measurements (see e.g. [98,99]).
Using the tools described above, we find the allowed/ forbidden areas in the (cos (β − α) , tan β) plane as shown in figure 5 16 . Some points are additionally ruled out by H/a → τ τ [103], H → h 125 h 125 [104] and H → a Z [105,106] searches. Values of cos (β − α) > 0.04 and tan β 5 are excluded from h 125 → Z Z [107]. Points which are still allowed after the application of HiggsBounds and HiggsSignals are outside the region sensitive to the full Run 2 A → H Z search [108], assuming on-shell production and decay 17 . Furthermore, the region where cos (β − α) −0.05, tan β 5 is mainly ruled out from both perturbative unitarity and perturbativity.
Constraints from relic density, i.e. the requirement that eqn. (16) is fulfilled, depends on many parameters. In general, however, small values of |m a − 2 m χ | can trigger efficient annihilation processes via the resonant s-channel diagram, and lead to relic density values between 10 −6 and the upper bound. For larger mass differences O (200 GeV) or larger, basically all points are excluded. Dominant annihilation channels are χχ → bb and χχ → tt, where the latter channel opens up for m a 2 m t and tends to lead to smaller values of relic density. We show the values of relic density as a function of the above mass difference in figure 6.
Similarly, for sin θ → 0 annihilation cross sections tend to become small, such that relic density predictions exceed the value measured by the Planck collaboration, cf. eq. (16). We display the dependence on this variable in figure 6, where in addition color coding indicates the mass difference m a − 2 m χ .
Considering only points that pass the relic density constraints, we found that direct detection does not impose any additional constraints. A more detailed analysis, as e.g. presented in [9,12,91,92], could modify this statement. However, from figure 24 in [6], we see that these constraints are effective only for masses m a 60 GeV and additionally depend on the mixing angle. We have verified that the few points that feature this mass range are not excluded when comparing to the above figure. Further comments on the sensitivity of future direct detection experiments for this model can e.g. be found in [12].
Before discussing constraints from direct searches, we will briefly summarize the effects the bounds considered above have imposed on the original parameter space • B-physics constraints set a lower bound on tan β as a function of m H ± ; in general, tan β > 1.
• Furthermore, oblique parameters reduce the allowed mass differences, especially in the THDM scalar sector.
• In the general scan, also the upper limits of m a in dependence on m H , m H ± are reduced by more than 100 GeV. However, this is an artefact of the scan, and one can easily repopulate these regions in a more finetuned setup.
• We also see that m H > 550 GeV, m A > 750 GeV. In both cases, a large number of points are excluded from S, T, U observables, as these variables favour small mass differences between the new scalars. If we force m A ≤ 750 GeV, for example, out of 1000 points only one survives which has very degenerate m H , m H ± masses, differing on the permill level. If we force m H ≤ 550 GeV, more points survive, which feature a small m H ± , m A difference, typically 10% difference, as well as small mass differences |m a − 2 m χ | 120 GeV. In a general, non-finetuned scan, the lower bound on m H ± forces the mass scale of H, A to be 530 − 550 GeV.
All other parameters still populate the original regions set in eqn. (18). Regarding direct collider constraints as discussed in section III C, we find the following results • The strongest constraints stem from h + / E ⊥ searches presented in [32], cutting out around ∼ 28% of the leftover parameter space. As we let all parameters float freely, limits in 2-dimensional planes are not straightforward. Using our approach, only points for which m a − 2 m χ ≥ 0 are affected. Furthermore, points where m χ ≥ 400 GeV are not affected. This however is e.g. due to a strong correlation between m χ and m a , which in turn influences the mass difference m A − m a and the A → h a branching ratio.
• The ℓℓ + / E ⊥ searches cut out roughly 9% of the available parameter space. As in [22], regions with m a 500 GeV are mainly affected. However, as we scan over all parameters, also many points for which m a 500 GeV are still allowed. In addition, the branching ratio for H → a Z swiftly goes to 0 when m H − m a → m Z , and only has reasonably large values when m H − m a 200 GeV, so that only points obeying this mass hierarchy are excluded. Similarly, as the above branching ratio ∼ sin 2 θ, we found that points for which | sin θ| 0.15 are not excluded.
• The bounds from the p p → H +t b, H + →bt searches [25] do not constrain the parameter space here, as production cross sections for the process p p → H +t b are O (10 −2 pb). This can be traced back to the lower limit on tan β from B-physics observables, cf. fig 1. The same holds for searches in the same production mode, but with different final states [30].
• Searches for tt + / E ⊥ cut out a relatively small region of parameter space,∼ 1 %, in the range where m a 300 GeV, m a − 2 m χ ∈ [0; 40] GeV, tan β 2. As before, several parameters contribute, so one also finds allowed points within the regions mentioned above. All of these are also subject to at least one other constraint from the searches mentioned above.
• None of the points is excluded by the W +t + / E ⊥ /W − t + / E ⊥ , bb + / E ⊥ and monojet searches considered here. Out of the remaining parameter space, we now investigate the magnitude of expected rates at possible future e + e − colliders. In the limit where sin θ → 0, we recover the decoupling scenario of a standard THDM. It is therefore interesting whether we can find regions in parameter space where novel signatures, and explicitely final states with missing energy, give the largest rates.
If we focus on the neutral sector and pair-production processes, possible final states are given by ha, hA, HA, Ha. For the first two processes, the sum of masses remains 1 TeV, so colliders with a center-of-mass energy in that range might be relevant. For the latter two, on the other hand, masses can range above 1 TeV, and we can consider cross sections at a 3 TeV collider.

Note that h [H]
A production in a THDM is mediated mainly via Z-boson exchange in the s-channel, where the corresponding coupling is proportional to cos (β − α) [sin (β − α)] respectively. As cos (β − α) is heavily constrained, the corresponding cross-sections are quite small, reaching O (10 −3 fb). Therefore, we will here concentrate on HA (a) production and decay. For this channel, we find that at a 3 TeV collider, production cross-sections can reach up to 1 fb, where largest cross section values are achieved for m A + m H ∼ 1400 GeV. Dominant decay modes for such points are displayed in figure 7. Decays without missing energy lead to tttt and tt tt Z final states, where the first in dominant in large regions of parameter space. The first decay mode that is novel with respect to standard THDMs is the tt + / E final state. In order to identify regions where this state dominates, we show the expected tttt and tt + / E cross sections in figure 8, where in the right plot we additionally include contributions mediated via Ha production. We see that indeed we can identify regions where tt + / E dominates and renders the largest rates. As an example, we here list The total widths are given by Γ H = 2.41 GeV, Γ A = 52.5 GeV, Γ a = 26.5 GeV, Γ H ± = 12.1 GeV, rendering all width/ mass ratios 6 %. The production cross sections at 3 TeVare given by σ HA = 0.512 fb, σ Ha = 0.390 fb. The H decays to tt with a branching ratio of 93%, while χχ final states have branching ratios of 64 %/95 % for A/a, respectively. This leads to an overall estimated production cross section of around 0.65 fb, using factorization and neglecting possible interference effects. The next step in this investigation would be a detailed simulation of signal and background in a CLIC-like environment, as e.g. performed in [109,110] for the Inert Doublet Model. This is in the line of future work.

VI. SUMMARY AND OUTLOOK
In this work, we for the first time presented a scan for the THDMa that lets all 12 free parameters of that model float freely, within ranges that were chosen to optimize scan performance. We have identified regions in parameter space that survive all current theoretical and experimental constraints, and provided a first estimate of possible production cross sections within this model at future e + e − facilities, with a focus on signatures not present in a standard THDM. We have included direct LHC search results using a simplified approach, with maximal cross sections as an upper limit. Such searches typically constrain the low-mass regions of parameter space, where the dark sector masses are 500 GeV, as well as regions with mixings | sin θ| 0.15. In general, however, clear boundaries in twodimensional planes are not easy to identify. An exception are constraints from B-physics observables, that expecially impose a lower limit on the charged mass ∼ 800 GeV that varies with tan β. Furthermore, constraints from electroweak precision observables in the form of oblique parameters pose relatively strong constraints on the mass differences in the THDMa scalar sector for novel scalars, leading to a general lower mass scale ∼ 500 GeV. Requring the relic density to lie below the current experimental measurement furthermore poses strong constraints on |m a − 2 m χ |. We consider this work as a first step towards a more detailed analysis of this model. This includes more detailed recasts of the current LHC experimental results, including possible forecasts for HL-LHC sensitivity, as well as more detailed investigations at future lepton machines, as e.g. CLIC or possible muon-colliders. First estimates for the former show that the most promising novel channel could be tt + / E in a CLIC-like scenario with a 3 TeV center-of-mass energy.
Finally, in view of recent events [111], one might wonder whether the anomalous magnetic momentum of the muon can be explained by the model discussed here. Contributions should be similar to the ones of a general THDM of type II, as e.g. discussed in [112,113], with additional contributions of the second pseudoscalar a, including appropriate rescaling. However, from the above work it becomes obvious that pseudoscalar masses would have to be light, 200 GeV, with additionally relatively large tan β values 15, in order to account for the observed experimental value. Such points have not been included in our analysis here. A quick check using the next-to-leading order approximation given in [114] shows that maximal values of ∆a BSM µ ∼ 2.5 × 10 −12 can be achieved with e.g. the points which pass all other (including experimental LHC) constraints, which is around three orders of magnitude lower than the observed discrepancy. In [52], regions are shown in the (m H ± , tan β) plane; here, for charged masses 800 GeV, tan β 7 seems to be required for 68% confidence level agreement with the result from [115]. This could a priori be an interesting region to investigate, which we will leave for future work.