Long-lived particle phenomenology in the 2HDM+$a$ model

Higgs decays displaced from the primary interaction vertex represent a striking experimental signature that is actively searched for by the ATLAS, CMS and LHCb collaborations. We point out that signals of this type appear in the context of the 2HDM+$a$ model if the mixing angle $\theta$ of the two CP-odd weak spin-0 eigenstates is tiny and the dark matter (DM) sector is either decoupled or kinematically inaccessible. Utilising two suitable benchmark scenarios, we determine the constraints on the parameter space of the 2HDM+$a$ model that are set by the existing LHC searches for long-lived particles (LLPs) in Higgs decays. We find that depending on the precise mass spectrum of the spin-0 states, mixing angles $\theta$ in the ballpark of a few $10^{-8}$ to $10^{-5}$ can be excluded based on LHC Run II data. This finding emphasises the unique role that searches for displaced signatures can play in constraining the parameter space of the 2HDM+$a$ model. The ensuing DM phenomenology is also discussed. In particular, we show that parameter choices leading to an interesting LLP phenomenology can simultaneously explain the DM abundance observed in today's Universe.


Introduction
The two-Higgs-doublet plus pseudoscalar model (2HDM+a) [1][2][3][4] has by now established itself as a pillar of the LHC dark matter (DM) search programme . It includes a DM candidate in the form of a Dirac fermion which is a singlet under the Standard Model (SM) gauge group, four 2HDM spin-0 states and an additional CP-odd mediator that is meant to provide the dominant portal between the dark and the visible sector. Since in models with pseudoscalar mediators the DM direct detection constraints are weaker compared to models with scalar mediators, the former models are more attractive from an astrophysical point of view since they often allow to reproduce the observed DM relic abundance in a wider parameter space and with less tuning. These features admit a host of missing transverse momentum (E miss T ) and non-E miss T signatures in the 2HDM+a model at colliders which can and have been consistently compared and combined. See for instance [20,22,27] for such combinations.
Beyond the SM (BSM) scenarios in which the hidden and the visible sectors are connected through a Higgs portal are also being actively probed for at colliders. One rather generic feature in such BSM models is the appearance of new electrically neutral longlived particles (LLPs) that give rise to displaced vertex signatures in the LHC detectors -see for example [28][29][30] for detailed reviews of theoretical and experimental aspects of LLPs at the LHC. The main goal of this article is to point out that besides interesting prompt E T,miss and non-E T,miss signatures, the 2HDM+a model can also have an attractive LLP phenomenology. In fact, in the 2HDM+a model the role of the LLP is played by the additional pseudoscalar a, which depending on its mass can be pair produced efficiently in the decays of both the 125 GeV Higgs and the non-SM CP-even Higgs, i.e. h → aa and H → aa. To illustrate the different facets of the LLP phenomenology in the 2HDM+a model we identify two suitable parameter benchmarks. For these benchmark scenarios we determine the bounds on the mixing angle θ of the two CP-odd weak spin-0 eigenstates as a function of the LLP mass that are set by the existing LHC searches for displaced Higgs decays [31][32][33][34][35][36][37][38][39][40][41][42][43][44][45][46][47]. It turns out that depending on the precise mass spectrum of the spin-0 states, mixing angles θ from around a few 10 −8 to about 10 −5 can be excluded with LHC Run II data. To the best of our knowledge, mixing angles θ in this range cannot be tested by any other means, which highlights the special role that LLPs searches play in constraining the parameter space of the 2HDM+a model. In fact, as we will further demonstrate, parameter choices that lead to an interesting LLP phenomenology can in general also correctly predict the measured DM relic density. The regions of 2HDM+a parameter space singled out in our article therefore deserve, in our humble opinion, dedicated experimental explorations in future LHC runs.
This work is structured as follows: in Section 2 we detail the theoretical ingredients that are relevant in the context of this article. Our general findings concerning the LLP phenomenology in the 2HDM+a model will be illustrated in Section 3 by considering two suitable parameter benchmark scenarios as examples. For these two benchmark choices we derive in Section 4 the constraints that the existing LHC searches for LLPs in displaced Higgs decays place on the 2HDM+a parameter space. In Section 5 we discuss the resulting DM phenomenology. Section 6 concludes our work.

2HDM+a model primer
In order to understand under which circumstances the additional pseudoscalar a in the 2HDM+a model can be an LLP it is useful to recall its partial decay modes -further details on the structure of the 2HDM+a model can be found for instance in [4,9]. In the alignment limit, i.e. cos (β − α) = 0, and choosing for concreteness the Yukawa sector of the 2HDM+a model to be of type-II, one has at tree level (2.1) At the one-loop level the pseudoscalar a can also decay to gauge bosons. The largest partial decay width is the one to gluon pairs. It takes the form Here m a is the mass of the pseudoscalar a, m χ is the mass of the DM particle, y χ is the Yukawa coupling of the pseudoscalar a to a pair of DM particles and sin θ quantifies the mixing of the two CP-odd weak spin-0 eigenstates. Furthermore, N q c = 3, N l c = 1, Examples of tree-level Feynman diagrams representing pp → aa production via gluon-gluon-fusion (ggF) Higgs production (left) and pp → l + l − aa production in associated Zh production (right) in the 2HDM+a model. The possible decay modes of the pseudoscalar a are not shown. Consult the text for further details. η u = cot β, η d = tan β, η l = tan β and y f = √ 2m f /v with m f the mass of the relevant SM fermion, v 246 GeV the Higgs vacuum expectation value (VEV) and α s the strong coupling constant. From the analytic expressions (2.1) and (2.2) it is evident that the pseudoscalar a can only be long-lived if sin θ is sufficiently small, i.e. sin θ → 0, and decays to DM are strongly suppressed/absent which can be achieved either via decoupling, i.e. y χ → 0, or by forbidding the process kinematically, i.e. m χ > m a /2.
Given the strong suppression of the couplings of the pseudoscalar a to SM fermions in the limit sin θ → 0, the only possibility to produce a long-lived a is via the decay of heavier spin-0 state φ into a pair of such pseudoscalars. In the case that the scalar potential is CP conserving the φ has to be a CP-even state which implies that in the 2HDM+a model one can have both decays of the 125 GeV Higgs h and the heavy CP-even Higgs H. The corresponding partial decay widths can be written as with φ = h, H. For sin θ 0 the relevant trilinear couplings are given by [4] g haa − 2v m h λ P 1 cos 2 β + λ P 2 sin 2 β , where m h 125 GeV is the mass of the SM-like Higgs, while λ P 1 and λ P 2 are the quartic couplings that appear in the 2HDM+a scalar potential as follows P 2 λ P 1 H † 1 H 1 +λ P 2 H † 2 H 2 (see for example [4,9] for the complete expression of the scalar potential). Here P denotes the additional pseudoscalar in the weak eigenstate basis which satisfies P a for sin θ 0.
The trilinear couplings entering (2.5) can be constrained phenomenologically. In the case of g haa one can require that the partial decay width Γ (h → aa) does not exceed the total decay width Γ h of the 125 GeV Higgs as measured directly at the LHC. For m a m h this leads to the inequality [22] where in the last step we have employed the latest 95% confidence level (CL) bound of Γ h < 1.1 GeV that follows from the direct LHC measurements of the total SM-like Higgs width [48,49]. Inserting the first expression of (2.5) into (2.6) then leads to the following relation λ P 1 cos 2 β + λ P 2 sin 2 β 0.24 . (2.7) In the case of g Haa one obtains in a similar fashion The above discussion should have shown that in the limit sin θ → 0, an LLP signature can arise in the 2HDM+a model from h or H production followed by the decay of the intermediate Higgs to a pair of pseudoscalars. Representative tree-level graphs showing pp → aa production in ggF Higgs production (left) and pp → l + l − aa production in associated Higgs production (right) that appear in the 2HDM+a model can be found in Figure 1. Notice that in the former case both the 125 GeV Higgs and the heavy CP-even Higgs contribute in the alignment limit. This is not the case for the latter process as the HZZ vertex vanishes identically in the limit cos (β − α) → 0.

Parameter benchmarks
Besides the phenomenological bounds (2.7) and (2.9) that constrain the sizes of the couplings λ P 1 and λ P 2 , the requirement for the scalar potential to be bounded from below also restricts the quartic couplings as well as other parameters of the 2HDM+a model. Assuming that λ P 1 , λ P 2 > 0 and that sin θ 0, one finds two bounded from below conditions that take the form [9] λ 3 > 2λ , Here the parameter λ 3 denotes the usual quartic coupling from the 2HDM scalar potential and λ = m 2 h /(2v 2 ) 0.13 is the cubic SM Higgs self-coupling. In order to fulfill these relations and to avoid the tight constraints from Higgs and electroweak precision physics, we make the following common parameter choices We furthermore employ a Yukawa sector of type-II throughout this work.
The first 2HDM+a benchmark scenario that we will study as an example to illustrate the possible LLP phenomenology in the 2HDM+a model is: We furthermore treat sin θ and m a as free parameters but require that m a < m h /2 so that the LLP can be pair produced in the decay of the 125 GeV Higgs boson. The precise value of the common heavy Higgs mass is irrelevant in such a situation and we simply set it to m H = 600 GeV in benchmark I. Notice that the quartic couplings λ P 1 and λ P 2 have been chosen such that the constraint (2.7) is easily fulfilled. In fact, in the limit m a → 0 the benchmark I parameter choices imply a value that is very close to the SM prediction of Γ SM h = 4.07 MeV [50]. The corresponding h → aa branching ratio is The proper decay length of the pseudoscalar a for masses in the range m a ∈ [20, 60] GeV can be approximated by a pseudoscalar of m a = 40 GeV has a proper decay length of around 1m. The result (3.6) includes higher-order QCD corrections employing the formulas presented in Appendix A of the paper [51] as implemented in [52]. Also notice that for the choices (3.3) and assuming that sin θ 0, the additional 2HDM Higgses are all narrow, i.e. Γ H /m H 2%, In our second 2HDM+a benchmark scenario that leads to an interesting LLP phenomenology, we consider the following parameter choices λ P 1 , λ P 2 , m χ = 3, 0, 770 GeV , (benchmark II) . (3.8) The parameters sin θ, m H and m a are instead treated as input with the requirements that m a > m h /2 and m a < m H /2 so that the LLP can only be pair produced in the decay H → aa of the heavy CP-even Higgs boson H. Notice that the values λ P 1 and λ P 2 in (3.8) satisfy the constraint (2.9). Taking for example m H = 600 GeV and m a = 150 GeV, the total decay width of the heavy CP-even Higgs is given by which implies that Γ H /m H = 3.7%. The corresponding branching ratios are meaning that the decays of the heavy Higgs to two LLPs does not have the largest branching ratio but that di-top decays are more frequent. Notice that given the structure of the trilinear coupling g Haa in (2.5) this feature will be even more pronounced for heavier CPeven Higgs bosons H. In the range m a ∈ [100, 300] GeV, the proper decay length of the pseudoscalar a is approximately given by where again the results of [51,52] have been used. It follows that for a pseudoscalar of m a = 150 GeV has a proper decay length of about 1m. Notice finally that in the case of (3

LLP constraints
At the LHC, searches for displaced Higgs boson decays into LLPs have been carried out by the ATLAS, CMS and LHCb collaborations in different final states, covering proper decay lengths from around 10 −3 m to 10 3 m [31 -47]. The LLP mean decay length determines the search strategies and reconstruction techniques that are employed -see for instance Section 5 of the review [20] for comprehensive descriptions of the details of the experimental techniques employed in LHC LLP searches.
We first consider the 2HDM+a benchmark I scenario (3.3) with m a < m h /2. In this case the pseudoscalar a can be pair produced in the decay of the 125 GeV Higgs boson. Depending on whether the mass of the a is below or above the bottom-quark threshold, the dominant decay modes of the pseudoscalar are BR (a → cc) 53%, BR (a → τ + τ − ) 38% and BR (a → gg) 10% or BR (a → bb) 85%, BR (a → cc) 4%, BR (a → τ + τ − ) 7% and BR (a → gg) 3%, respectively. LLP searches that target pseudoscalar pair production in ggF Higgs or associated Zh production (cf. Figure 1) leading to multi-jet or four-bottom final states therefore provide the most stringent constraints. Looking for displaced leptons instead leads to significantly weaker restrictions because of the small leptonic branching ratios.
In Figure 2 we show an assortment of LLP constraints in the m asin θ plane that apply in the case of (3.3). All limits result from LHC searches that consider ggF Higgs production. The dotted red exclusion corresponds to the search [34] that considers displaced hadronic jets in the ATLAS calorimeter (CM) and the muon spectrometer (MS) [33], while the dotted blue constraint instead results from the ATLAS search [35] that utilises the inner detector (ID) and the MS. These searches use up to 36 fb −1 and 33 fb −1 of √ s = 13 TeV data, respectively. The dotted green (purple) lines represent an upgrade of the MS (CM) search strategy to 139 fb −1 of luminosity collected in LHC Run II. The corresponding limits are reported in the ATLAS publication [43] and [44], respectively. The dashed yellow contour is finally the exclusion that derives from the CMS search [38] which employs the 95% CL exclusion regions in the m asin θ plane for the 2HDM+a benchmark I scenario (3.3). The dotted red, blue, green and purple lines correspond to the limits following from the ATLAS searches [33,34], [35], [43] and [44], respectively. The dashed yellow curves instead represent the bound that arises from the CMS search [38]. The parameter space between the lines is disfavoured. See main text for further details. muon endcap and 137 fb −1 of √ s = 13 TeV data. From the figure it is evident that in the 2HDM+a benchmark I scenario the existing LHC searches for displaced Higgs decays to hadronic jets allow to exclude values of sin θ between around 10 −7 and 10 −5 with the exact bound depending on the mass of the pseudoscalar a. The excluded parameter space corresponds to proper decay lengths cτ a in the range from around 59 m to 0.08m. Notice that given the smallness of the h → aa branching ratio cf. (3.5) , our benchmark I scenario easily evades the present bounds on the undetected or invisible branching ratios of the 125 GeV Higgs [53] that amount to 19% and 9%, respectively. In fact, even a possible future high-luminosity LHC (HL-LHC) upper limit on the invisible branching ratio of the SM-like Higgs of BR (h → invisible) < 2.5% [54] would not be stringent enough to test (3.3) indirectly. This feature underlines the special role that LLP searches for displaced Higgs decays can play in testing 2HDM+a models with mixing angles θ close to zero.
Let us now turn our attention to the benchmark II scenario (3.8). In this case the parameters are chosen such that an LLP signal may arise from the prompt decay of the heavy CP-even Higgs, i.e H → aa, followed by the displaced decays of the pseudoscalars to a pair of SM fermions a → ff or gluons a → gg. Given our choice of Yukawa sector and tan β, the a dominantly decays to the heaviest SM fermion, which means that depending on the precise value of its mass either a → bb or a → tt provide the largest rate. To illustrate these two possibilities we consider in benchmark II the mass combination m H = 600 GeV   [33,34], [35], [43] and [44], respectively. The dotted red exclusion in the lower panel instead represents the combination of the ATLAS searches [33,34,43,44]. The shaded parameter regions are disfavoured. Further explanations can be found in the main text.
In the upper panel of Figure 3 we display the relevant 95% CL exclusion regions in the m asin θ plane that apply in the case of the 2HDM+a benchmark II scenario with m H = 600 GeV. One observes that taken together the ATLAS searches [33-35, 43, 44] allow to exclude sin θ values between 2 · 10 −8 and 2 · 10 −6 . The corresponding proper decay lengths cτ a range from 53 m to 0.04m. The lower plot in Figure 3 shows the limits on sin θ that a combination of the four ATLAS searches [33,34,43,44] allow to set in the 2HDM+a benchmark II scenario (3.8) assuming m H = 1000 GeV. One observes that the existing LLP searches can exclude mixing angles for the mass points m a = 50 GeV, 150 GeV, 275 GeV and 400 GeV, while the 2HDM+a realisation with m a = 475 GeV remains untested at present. For pseudoscalar masses below the top-quark threshold, sin θ values between around 4·10 −8 and 2 · 10 −6 are excluded, whereas for m a = 400 GeV mixing parameters in the range of about 4 · 10 −9 and 1 · 10 −9 are disfavoured. The excluded parameter space corresponds to cτ a values ranging from around 9.7m to 0.06m. Notice that the order of magnitude improvement of the constraint on sin θ from the point m a = 275 GeV to m a = 400 GeV is readily understood by recalling that in the former case one has Γ a 16 MeV sin 2 θ, while in the latter case Γ a 12 GeV sin 2 θ as a result of the open a → tt channel. We add that improving the limits [43,44] by a factor of four would allow to probe our 2HDM+a benchmark II scenario for m H = 1000 GeV and m a = 475 GeV. A final remark concerns the possibility to search for the heavy CP-even Higgs in the processes pp → H → tt or pp → ttH → 4t. While LHC searches for spin-0 resonances in both di-top [56][57][58] and four-top production [59][60][61][62] have been performed, it turns out that the existing searches do not provide any bound on our 2HDM+a benchmark II model for both m H = 600 GeV and m H = 1000 GeV in the m a ranges considered in Figure 3. This finding again stresses the unique role that LLPs searches can play in the 2HDM+a model in constraining the parameter space.

Relic density
In order to understand the physics of standard thermal relic freeze-out in 2HDM+a realisations with sin θ 0, we first write the cross section for annihilation of DM into a final state X as where v rel is the relative velocity of the DM pair and the coefficient σ 0 X (σ 1 X ) describes the s-wave (p-wave) contribution.
In the alignment limit, the possible DM annihilation channels involving a pseudoscalar a are χχ → a → ff , χχ → a → ZH, χχ → a → ah and χχ → a → aH for what concerns s-channel processes and χχ → aa with DM exchange in the t-channel (cf. also [9]). The annihilation cross sections (5.1) of the former two reactions are, however, proportional to sin 2 θ making them numerically irrelevant in the limit sin θ → 0 unless m a = m χ /2. Such highly tuned solutions to the DM miracle will not be considered in what follows. Similarly, all DM annihilation contributions involving the exchange of a heavy pseudoscalar A are suppressed by at least two powers of the sine of the mixing angle θ, so that only the processes depicted in Figure 4 are relevant for the calculation of the DM abundance in the context of this work.
The annihilation process χχ → a → ah proceeds via s-wave and we find for the corresponding coefficient the following analytic result where the expression for g haa in the limit sin θ → 0 can be found in the first line of (2.5) and Γ a denotes the total decay width of the pseudoscalar a. Since σ 0 ah = 0 we ignore the p-wave coefficient σ 1 ah below by setting it to zero. The result for the s-wave coefficient σ 0 aH describing DM annihilation through χχ → a → aH is simply obtained from (5.2) by the replacements g haa → g Haa and m h → m H . In the case of χχ → aa the annihilation cross section is instead p-wave suppressed (see [63] for the calculation of the t-channel contribution in the simplified pseudoscalar DM model) and the corresponding expansion coefficients take the form σ 0 aa = 0 and Using the velocity expansion (5.1) the DM relic density after freeze-out can be approximated by Here 20,30] with T f the freeze-out temperature and the sum over X in principle includes all possible final states. As we have explained above, for sin θ 0 and away from the exceptional points m a = m χ /2, however, only the channels X = ah, aH, aa Benchmark I Figure 5. Predicted DM relic abundance in the m am χ plane for the 2HDM+a benchmark I parameter choices (3.3). The contour lines indicate the value of Ωh 2 /0.12, meaning that the regions below (above) 1 correspond to a DM underabundance (overabundance) in today's Universe. For additional details we refer the interested reader to the main text.
are numerically important. In the limit of heavy DM, i.e. m χ m a , m h , m H , the velocityaveraged annihilation cross section at the freeze-out temperature can be further simplified: This approximation shows that the s-channel (t-channel) contributions to σv rel f scale as 1/m 4 χ 1/m 2 χ in the limit of infinitely heavy DM. The formulas (5.4) and (5.5) represent useful expressions to estimate Ωh 2 . In the case of the benchmark I scenario (3.3) one has g 2 haa 6 · 10 −5 and g 2 Haa = 0, and it is thus a good approximation to neglect the s-channel contributions to σv rel f . It follows that (5.6) Using x f 25 the relic abundance of Ωh 2 = 0.120 ± 0.001 as determined by Planck [64] is therefore obtained in the case of (3.3) for DM masses m χ 160 GeV while for parameter regions with m χ 160 GeV (m χ 160 GeV) one expects DM underabundance (overabundance). These expectations agree quite well with the results of our exact relic calculation that have been obtained with MadDM [65] and are shown in Figure 5. In fact, the exact computation for (3.3) and m a = 30 GeV leads to Ωh 2 = 0.118, while (5.6) naively predicts a value that is larger by around 15%. The observed difference can be traced back to the fact Benchmark II, m H = 1000 GeV that the MadDM calculation gives x f 21 in the parameter region of interest and correctly takes into account the phase-space suppression present in (5.3) due to the non-zero values of m 2 a /m 2 χ . We add that χχ → a → ah annihilation represents a relative contribution to Ωh 2 of less than about 1% in the part of the m am χ plane that is depicted in the figure.
Neglecting the s-channel contributions in the approximation (5.6) is hence fully justified.
In the case of our 2HDM+a benchmark II parameter scenario (3.8) the coupling g haa is no longer small, in fact g 2 haa 35 and we furthermore have g 2 Haa = 0. On the other hand, y 2 χ /x f 0.04 and thus one can neglect the χχ → aa contribution to the velocity-averaged annihilation cross section at the freeze-out temperature (5.5) to first approximation. Doing this, one obtains the following simple expression  Figure 6. We find that for the parameters (3.8) together with m a = 150 GeV as well as m H = 600 GeV, the exact calculation predicts Ωh 2 = 0.117, a value less than 5% below the naive expectation. Numerically, we furthermore obtain that the relative contribution of χχ → aa to the DM relic density is always below 1% in benchmark II with m H = 600 GeV, showing that it can be safely neglected in the derivation of (5.7).
The result of our MadDM scan in the 2HDM+a benchmark II parameter scenario with m H = 1000 GeV is presented in the lower panel of Figure 6. It is evident from the plot that in this case the above simplistic formula is not able to capture the more intricate behaviour of the contours of constant relic density. In fact, this is not a big surprise because in the derivation of (5.7) we have assumed that m χ m a , m h , m H , however, one has m χ < m H in the entire m am χ plane considered. Still the values m χ ∈ [680, 740] GeV of the DM mass that lead to the correct DM abundance for m a ∈ [50,500] GeV are only less than 10% smaller than what one would expect from (5.7). Like in the case before, it turns out that the relative contribution of s-channel annihilation amounts to more than 99% of the predicted values of Ωh 2 . This shows again that DM annihilation via χχ → aa is phenomenologically irrelevant in 2HDM+a realisations of the type (3.8).
Before concluding, we add that DM direct detection experiments do not set relevant constraints on the 2HDM+a benchmarks (3.3) and (3.8) for the mixing angles θ 0 necessary to have a long-lived a. This is a simple consequence of the fact that the spinindependent DM-nucleon cross section is suppressed by both a loop factor and two powers of sin θ -see the recent articles [8,22,[66][67][68] for explicit formulas and further explanations.

Conclusions
The main lesson that can be learnt from the analytic and numerical results presented in this work is that LHC searches for displaced Higgs decays can provide unique constraints on 2HDM+a model realisations. In fact, LLP signatures appear in the context of the 2HDM+a model if the mixing angle θ of the two CP-odd weak spin-0 eigenstates is very small and the DM sector is either decoupled or kinematically inaccessible. In order to emphasise this generic finding, we have studied two distinct parameter benchmarks and explored the sensitivity of the existing LHC LLP searches by performing parameter scans in the m asin θ plane. The results of these scans can be found in Figures 2 and 3.
In the benchmark I scenario we have chosen the 2HDM+a parameters such that the 125 GeV Higgs boson gives rise to an LLP signature through its prompt decay h → aa followed by the displaced decays of the pseudoscalars to SM fermions such as a → bb or gluons. One important feature that is nicely illustrated in our benchmark I scan is that the LLP searches for displaced hadronic jets that have been performed at LHC Run II can probe regions of parameter space with mixing angles θ in the ballpark of 10 −7 to 10 −5 that are presently not accessible by any other means. In fact, in our benchmark I scenario the predicted h → aa branching ratio turns out to be below the target sensitivity that the HL-LHC is expected to reach on undetected or invisible decays of the 125 GeV Higgs boson. 2HDM+a realisations like (3.3) are therefore unlikely to be testable indirectly at the LHC in searches for both prompt h → aa → 4f production or through signatures involving a significant amount of E T,miss such as h → invisible or mono-jet final states.
The benchmark II parameters were instead chosen such that the LLP signal arises from the prompt decay H → aa of the heavy CP-even Higgs followed by the displaced decays of the pseudoscalars to a pair of SM fermions or gluons. Depending on the precise value of the LLP mass either the a → bb or the a → tt mode turns out to provide the largest rate. In order to illustrate these two distinct possibilities, we have considered in benchmark II the mass combination m H = 600 GeV with m a ∈ [50,275] GeV as well as m H = 1000 GeV with m a ∈ [50,475] GeV. In the former case we found that the existing LHC searches for displaced heavy Higgs decays provide stringent constraints on θ, excluding values of the mixing angle between 2 · 10 −8 and 2 · 10 −6 . The limits on the benchmark II scenario with m H = 1000 GeV turned out to be noticeably weaker than the bounds for m H = 600 GeV due to the order of magnitude smaller inclusive heavy Higgs production cross section. Still, for the three mass values m a = 50 GeV, 150 GeV and 275 GeV, θ values between 4 · 10 −8 and 2 · 10 −6 are excluded, whereas for m a = 400 GeV mixing parameters in the range of 4 · 10 −9 and 1 · 10 −9 are disfavoured by the searches for the displaced heavy Higgs decays performed at LHC Run II. We expect that LLP searches at LHC Run III and beyond will provide sensitivity to 2HDM+a models à la (3.8) with heavy CP-even Higgs masses of the order of 1 TeV and pseudoscalar masses m a above the top-quark threshold. Let us finally note that in future LHC runs it should also be possible to probe 2HDM+a benchmark II models with a heavy CP-even Higgs that satisfies m H > 2m a and m H > 2m t through spin-0 resonance searches in di-top and four-top production. A detailed analysis of this issue is however beyond the scope of this work.
We have furthermore demonstrated that parameter choices that give rise to an interesting LLP phenomenology can simultaneously explain the observed DM abundance without excessive tuning. Our relic density scans have been given in Figures 5 and 6. In our benchmark I scenario, we have seen that due to the smallness of the g haa coupling the value of Ωh 2 is entirely determined by the t-channel process χχ → aa which is p-wave suppressed. For the range m a < m h /2 relevant in benchmark I, we found that the correct relic abundance can be obtained for DM masses of around 170 GeV, rather independently of the precise value of m a . In the case of the benchmark II model instead, only the s-channel annihilation contributions χχ → a → ah and χχ → a → aH that are s-wave turn out to be phenomenologically relevant. While the precise value of the DM mass for which Ωh 2 0.12 is achieved depends on the chosen parameters such as m a and m H , we saw that DM masses in the range m χ ∈ [680, 770] GeV allow to saturate the measured relic abundance in the part of the parameter space most relevant for the LLP signals triggered by H → aa.
The results and parameter benchmarks presented in our article should provide a useful starting point for interpretations of future ATLAS, CMS and LHCb searches for displaced Higgs decays to hadronic jets in the context of the 2HDM+a model. We therefore encourage and look forward to experimental explorations in this direction.