Dark Matter Signals at the LHC from a 3HDM

We analyse new signals of Dark Matter (DM) at the Large Hadron Collider (LHC) in a 3-Higgs Doublet Model (3HDM) where only one doublet acquires a Vacuum Expectation Value (VEV), preserving a parity $Z_2$. The other two doublets are \textit{inert} and do not develop a VEV, leading to a {\it dark scalar sector} controlled by $Z_2$, with the lightest CP-even dark scalar $H_1$ being the DM candidate. This leads to the loop induced decay of the next-to-lightest scalar, $H_2 \to H_1 f \bar f$ ($f=u,d,c,s,b,e,\mu,\tau$), mediated by both dark CP-odd and charged scalars. This is a smoking-gun signal of the 3HDM since it is not allowed in the 2HDM with one inert doublet and is expected to be important when $H_2$ and $H_1$ are close in mass. In practice, this signature can be observed in the cascade decay of the SM-like Higgs boson, $h\to H_1 H_2\to H_1 H_1 f \bar f$ into two DM particles and di-leptons/di-jets, where $h$ is produced from either gluon-gluon Fusion (ggF) or Vector Boson Fusion (VBF). However, this signal competes with the tree-level channel $q\bar q\to H_1H_1Z^*\to H_1H_1f \bar f$. We devise some benchmarks, compliant with collider, DM and cosmological data, for which the interplay between these modes is discussed. In particular, we show that the resulting detector signature, $\Et\; f \bar f$, with invariant mass of $ f \bar f$ much smaller than $m_Z$, can potentially be extracted already during Run 2 and 3. For example, the $H_2\to H_1\gamma^*$ and $\gamma^*\to e^+e^-$ case will give a spectacular QED mono-shower signal.


Introduction
The Higgs mechanism of Electro-Weak Symmetry Breaking (EWSB) seems to be the one chosen by Nature to assign mass to fermions and weak gauge bosons. In its minimal realisation, through a single Higgs doublet, it implies the existence of a single Higgs boson, as discovered in 2012 at the Large Hadron Collider (LHC). Indeed such a minimal Standard Model (SM) is compatible with a myriad of experimental results. However, the well known unanswered questions such as the origin of flavour, with three families of quarks and leptons, as well as Dark Matter (DM), suggest that some extension beyond the SM (BSM) is necessary.
Given the existence of three families of quarks and leptons, it is not so far fetched to imagine that there might also be three families of Higgs doublets, where, as for the fermions, the replication is not prescribed by the SM gauge group. Indeed, it is possible that the three families of quarks and leptons could be described by the same symmetries that describe the three Higgs doublets. In such scenarios, this generation/family symmetry could be spontaneously broken along with the EW symmetry, although some remnant subgroup could survive, thereby stabilising a possible scalar DM candidate. For certain symmetries, it is finally possible to find a Vacuum Expectation Value (VEV) alignment that respects the original symmetry of the potential which will then be responsible for the stabilisation of the DM candidate. In such 3-Higgs-Doublet Models (3HDMs), amongst the various symmetries which can govern them [1]- [4], a simple possibility is a single Z 2 , referred to here as Higgs parity, which can prevent Flavour Changing Neutral Currents (FCNCs) and possible charge breaking vacua.
In the present paper, we shall focus on the phenomenology of one of these 3HDMs, namely, the one in which the third scalar doublet is even and the first and second inert 1 doublets are odd under the Z 2 parity. We assume a vacuum alignment in the 3HDM space of (0, 0, v) that preserves the Z 2 symmetry (i.e., the Higgs parity). Thus we are led to consider a model with two inert doublets plus one Higgs doublet (I(2+1)HDM). This model may be regarded as an extension of the model with one inert doublet plus one Higgs doublet (I(1+1)HDM) 2 proposed in 1976 [5] and studied extensively for the last few years (see, e.g., [6]- [8]), by the addition of an extra inert scalar doublet. The lightest neutral scalar or pseudoscalar field amongst the two inert doublets, which are odd under the Z 2 parity, provides a viable DM candidate which is stabilised by the conserved Z 2 symmetry, displaying phenomenological characteristics notably different from the candidate emerging from the I(1+1)HDM case [9], both in the CP-Conserving (CPC) and CP-Violating (CPV) cases, as noted in Refs. [10]- [12]. Within this framework, we study some new SM-like Higgs decay channels offered by the extra inert fields, with the intent of isolating those which would enable one to distinguish between the I(2+1)HDM and I(1+1)HDM, assuming CP conservation throughout. The analysis of the CPV I(2+1)HDM is postponed to a future publication.
In particular, we shall focus on the loop induced decay of the next-to-lightest scalar, H 2 → H 1 ff (f = u, d, c, s, b, e, µ, τ ), mediated by loops involving both dark CP-odd and charged scalars. This decay chain occurs in the I(2+1)HDM but not in the I(1+1)HDM, so it enables the two models to be distinguished. In practice, the loop decay can be observed in the cascade decay of the SM-like Higgs boson into two DM particles and a fermion-antifermion pair, h → H 1 H 2 → H 1 H 1 ff , wherein the h state is produced from gluon-gluon Fusion (ggF) (i.e., gg → h) or Vector Boson Fusion (VBF) (i.e., qq ( ) → qq ( ) h). Notice, however, that this mode competes with the tree-level channel qq → H 1 H 1 Z * → H 1 H 1 ff present also in the I(1+1)HDM. The resulting detector signature, E T ff , with the ff invariant mass well below the Z mass, would indicate the presence of such a loop decay onset by a small difference between H 2 and H 1 which would in turn identify a region of I(2+1)HDM parameter space largely precluded to the tree-level process. Indeed, we will show that such a distinctive signature can possibly be extracted at the LHC during Run 2 and/or Run 3. In fact, amongst the possible ff cases, a particularly spectacular one would be the one in which an electron-positron pair is produced, eventually yielding an isolated mono-shower signal of QED nature, owing to the fact that the dominant component (over the box topologies) of the loop signal is the H 2 → H 1 γ * one, where the photon is (necessarily, because of spin conservation) off-shell, yet eventually producing the e + e − pair in configurations where the fermions are soft and/or collinear. In assessing the scope of the LHC in accessing this phenomenology, we shall consider all available theoretical [13,14] and experimental constraints affecting the I(2+1)HDM parameter space, so as to eventually define some benchmark scenarios which can be tested at the CERN machine.
The layout of the paper is as follows. In the next section we describe the CPC I(2+1)HDM. In Sect. 3, we introduce and discuss the aforementioned loop cascade decays. In Sect. 4 we perform all necessary calculations, both at tree and loop level, including analytic formulae for the H 2 → H 1 ff case. In Sect. 5, we present our results. We then conclude in Sect. 6. Finally, two appendices will collect some key formulae. 2 The CP conserving I(2+1)HDM 2.1 The potential with a Z 2 symmetry It is known [1] that, in a model with several Higgs doublets, the scalar potential which is symmetric under a group G of phase rotations can be written as the sum of V 0 , the phase invariant part, and V G , a collection of extra terms ensuring the symmetry group G.
Here, we study a 3HDM symmetric under a Z 2 symmetry with generator where the doublets, φ 1 , φ 2 and φ 3 , have odd, odd and even Z 2 quantum numbers, respectively. Note that this Z 2 generator forbids Flavour Changing Neutral Currents (FCNCs) and is respected by the vacuum alignment (0, 0, v), since the fermions which only couple to the active scalar doublet, φ 3 , are assigned an even Z 2 charge. The potential symmetric under the Z 2 symmetry in (1) can be written as This potential has only a Z 2 symmetry and no larger accidental symmetry 3 . We shall not consider CP violation in this paper, therefore we require all parameters of the potential to be real.
The full Lagrangian of the model is as follows: where L SM gf is the boson-fermion interaction as in the SM, L scalar describes the scalar sector of the model and L Y (ψ f , φ 3 ) describes the Yukawa interaction with φ 3 the only active doublet to play the role of the SM-Higgs doublet. The kinetic term in L scalar has the standard form of , does not change the phenomenology of the model. The coefficients of these terms, therefore, have been set to zero for simplicity.

Mass eigenstates
The minimum of the potential is realised for the following point: with The mass spectrum of the scalar particles are as follows.
• The fields from the active doublet The third doublet, φ 3 plays the role of the SM-Higgs doublet, hence, the fields G 0 , G ± are the would-be Goldsone bosons and h the SM-like Higgs boson with mass-squared which has been set to (125 GeV) 2 in our numerical analysis.
• The CP-even neutral inert fields The pair of inert neutral scalar gauge eigenstates, H 0 1 , H 0 2 , are rotated by into the mass eigenstates, H 1 , H 2 , with squared masses • The charged inert fields The pair of inert charged gauge eigenstates, φ ± 1 , φ ± 2 , are rotated by into the mass eigenstates, H ± 1 , H ± 2 , with squared masses m 2 • The CP-odd neutral inert fields The pair of inert pseudo-scalar gauge eigenstates, A 0 1 , A 0 2 , are rotated by into the mass eigenstates, A 1 , A 2 , with squared masses (The model is CP conserving, therefore there is no mixing between CP-even and CP-odd states in the inert sector.) We can separate the inert particles into two families, or generations, with the second generation being heavier than the respective fields from the first generation. We will refer to the set of (H 1 , A 1 , H ± 1 ) as the fields from the first generation and to (H 2 , A 2 , H ± 2 ) as the fields from the second generation.
Each of the four neutral particles could, in principle, be the DM candidate, provided it is lighter than the other neutral states. In what follows, without loss of generality, we assume the CP-even 4 neutral particle H 1 from the first generation to be lighter than all other inert particles, that is: In the remainder of the paper the notations H 1 and DM particle will be used interchangeably and so will be their properties, e.g., m H 1 and m DM .

Simplified couplings in the I(2+1)HDM
Due to the large number of free parameters in the I(2+1)HDM, which makes it impractical to analyse the model in the general case, we focus on a simplified case where the 4 Other neutral scalars could also play the role of DM candidate, e.g., A 1 would be the lightest particle after transformation λ 2,3 → −λ 2,3 . We could also choose H 2 to be the lightest particle with µ 2 12 → −µ 2 12 , or A 2 if both λ 2,3 → −λ 2,3 and µ 2 12 → −µ 2 12 . Hence, the results of our analysis are also applicable to all neutral scalars following suitable sign changes. parameters related to the first inert doublet are n times the parameters related to the second doublet [10]: without introducing any new symmetry to the potential. The motivation for this simplified scenario is that in the n = 0 limit the model reduces to the well-known I(1+1)HDM. We assume no specific relation among the other parameters of the potential. It is important to note that the remaining quartic parameters, (λ 1,11,22,12 , λ 12 ), do not influence the discussed DM phenomenology of the model and thus their values have been fixed in agreement with the constraints discussed in Sect. 2.4 and compliant with the results on unitarity obtained in [14].
With this simplification, it is possible to obtain analytical formulae for the parameters of the potential in terms of chosen physical parameters. In this study, we choose the set (m H 1 , m H 2 , g H 1 H 1 h , θ a , θ c , n) as the input parameters where g H 1 H 1 h is the Higgs-DM coupling. The meaningful parameters of the model are then defined as follows: The mixing angle in the CP-even sector, θ h , is given by the masses of H 1 and H 2 and the dark hierarchy parameter n: Notice that we restore the n = 1 limit of dark democracy discussed in [10,11,12] with θ h = π/4. For the correct definition of tan 2 θ h , the following two relations need to be satisfied: Without loss of generality, we can limit ourselves to n < 1, which will correspond to tan 2θ > 0 for θ h < π/4. Reaching other values of n is a matter of reparametrisation of the potential.

Theoretical and experimental constraints
As discussed in [10,11,12], the I(2+1)HDM is subject to various theoretical and experimental constraints.
In [10], we have studied in detail the theoretical constraints, namely the positivity of the mass eigenstates, boundedness of the potential and positive-definiteness of the Hessian. Our parameter choice is also compliant with the EW Precision Test (EWPT) bounds [10,11]. These limits have been taken into account in the present paper. The second set of experimental constraints comes from the relic abundance of DM as well as dedicated direct and indirect searches for DM particles. The Planck experiment provides a DM relic density limit of [15]: In this work, we do not focus on the details of DM annihilation (for detailed discussions see Refs. [10,11,12]). However, we require that the DM candidate of the I(2+1)HDM is in agreement with the upper limit from Planck (24) for all considered points. If relation (24) is exactly satisfied, then H 1 provides 100% of the DM in the Universe. This is a case in benchmark scenario A50 discussed in later sections [10,11]. We also consider cases where H 1 has a subdominant contribution and the missing relic density is to be provided by an extension of the model. This usually happens where mass splittings between H 1 and other inert particles are small, i.e., in the forthcoming benchmarks I5 and I10. In these two cases, the coannihilation channels of H 1 A i → Z → f f are strong and reduce DM relic density to values below the Planck value, even for very small values of Higgs-DM coupling. Benchmark scenario A50 (for 53 GeV m H 1 73 GeV) is in agreement with the most recent direct [16] and indirect [17] detection limits. However, for completeness, we show a larger mass region (40 GeV m H 1 90 GeV) in our cross section plots, and highlight the surviving regions.
For benchmarks I5 and I10, which -as mentioned -correspond to relic density below the Planck value, detection limits should be rescaled, leading to the (relic density dependent) limit of: We ensure this limit is satisfied for all studied points. The detailed analysis of astrophysical signals in benchmarks I5 and I10 is beyond the scope of this paper. However, for all masses in these benchmarks, relic density is within 10% -90% of the observed relic density. The missing relic density can be easily augmented by late-stage decays of an additional particle. The natural candidate here for the completion of the model would be a heavy right handed neutrino in the same vein as the scotogenic model [5], which would decay into DM after the thermal freeze-out of DM and bring back the under-abundant DM relic into the observed range. Finally, we take into account collider data from LEP and the LHC (including the Higgs total decay width [18], Higgs invisible decays [19], direct searches for additional scalars and the Branching Ratio (BR) for h → γ γ [19]), as discussed in [10,11,12]. In all cases, the mass splittings are large enough not to influence the decay widths of the weak gauge bosons, forbidding the on-shell decays Z → H 1,2 A 1,2 and W ± → H ± 1,2 H 1,2 /A 1,2 .
If the Higgs-DM coupling is small enough, i.e., g hH 1 H1 0.02, then both the Higgs invisible decay BR and Higgs total decay width are in agreement with measured values. For benchmark scenario A50, exclusions obtained from applying the LHC constraints are similar to those from dedicated DM experiments, excluding m H 1 53 GeV for a large Higgs-DM coupling. Benchmarks I5 and I10 are in agreement with these constrains for all studied masses.
Charged scalars in all cases are significantly heavier and short-lived than the neutral particles, therefore bounds from long-lived charged particle searches do not apply here. In all benchmarks, in particular I5 and I10, where all mass splittings are of the order of a few GeV, all heavier inert particles decay inside the detector.

Inert cascade decays
In the model studied here, there is one absolutely stable particle, H 1 , as its decays into SM particles are forbidden by the conservation of the Z 2 symmetry. By construction, all other inert particles, which are also odd under the Z 2 symmetry, are heavier than H 1 and hence unstable. The decays of these heavier inert particles may provide striking experimental signals for the I(2+1)HDM.
Access to the inert sector can be obtained through the SM-like Higgs particle, h, and/or the massive gauge bosons, Z and W ± , with the heavy inert particle subsequently decaying into H 1 and on-or off-shell W ± /Z/γ states. In fact, in this model, h can decay into various pairs of inert particles, leading to different signatures. We will consider here h → H 2 H 1 decays. In such a case, as intimated, we will consider Higgs production at the LHC through ggF and VBF.
The interesting production and decay patterns may occur both at tree-and looplevel. In the former case, the colliding protons produce an off-shell gauge boson Z * , which can in turn give us a H 1 A i pair (i = 1, 2), followed by the decay of A i into H 1 Z ( * ) → H 1 ff . In the latter case, one would produce a h state decaying into H 1 H 2 → H 1 H 1 ff , via the loop decay H 2 → H 1 ff . In both cases, one ends up with a E T ff signature (possibly accompanied by a resolved forward and/or backward jet in case of VBF and an unresolved one in ggF), i.e., a di-lepton/di-jet pair, which would generally be captured by the detectors, alongside missing transverse energy, E T , induced by the DM pair. Here, f = u, d, c, s, b, e, µ, τ . For the cases in which the mass difference , only the electron-positron signature would emerge, thus leading to the discussed Electro-Magnetic (EM) shower.
It is important here to notice that the loop decay chain initiated by h → H 1 H 2 is specific to the I(2+1)HDM case, while the one induced by A 1 → H 1 Z ( * ) may also pertain to the I(1+1)HDM case. (In fact, neither H 2 nor A 2 exists in the I(1+1)HDM, unlike A 1 .) Moreover, when the decays are non-resonant, there is no way of separating the two A i (i = 1, 2) patterns. In contrast, the extraction and observation of the decay h → H 1 H 2 (followed by the loop decay H 2 → H 1 ff ) would represent clear evidence of the I(2+1)HDM.
In the upcoming subsections, we will discuss the aforementioned tree-and loop-level decay modes of inert states into the DM candidate in all generality, then we will dwell on the features of the E T ff signature.

Tree-level decays of heavy inert states
CP-odd and charged scalars can decay at tree-level into a lighter inert particle in association with a real(virtual) gauge boson W ±( * ) or Z ( * ) . Assuming the mass ordering m H 1,2 < m A 1,2 < m H ± 1,2 , the following tree-level decays appear (only diagrams with H 1 in the final state are shown in Fig. 1, diagrams (A) and (B)): The leptonic decays(splittings) of real(virtual) massive gauge bosons will result in ff pairs for Z ( * ) and ff for W ±( * ) . The above processes are governed by the gauge couplings and therefore lead to small decay widths, of order 10 −2 − 10 −4 GeV, of heavy inert particles. However, these decay widths could grow if the mass splitting between H 1 and other particles is large. Note that, even if all particle masses are relatively close (of the order of 1 GeV), they all still decay inside the detector. The heavy CP-even scalar, H 2 , cannot couple to H 1 through Z ( * ) , since CP symmetry is conserved in our model. It can decay into the H 1 particle plus a Higgs boson (diagram (C) in Fig. 1), which will then decay via the established SM patterns. Depending on the mass splitting between H 1 and H 2 , the Higgs particle can be highly off-shell (recall that its SM-like nature requires its width to be around 4 MeV), thus leading to a relatively small decay width of H 2 and its relatively long lifetime. However, in all studied points, this width is not smaller than 10 −11 GeV, ensuring the decay of H 2 inside the detector 5 .
Therefore, the H 1 is the only truly invisible dark particle in the benchmark scenarios we consider in the I(2+1)HDM.

Loop-level decays of heavy inert states
Apart from the above tree-level decays there is also the possibility of loop-mediated ones for a heavy neutral inert particle, denoted in Fig. 2 as H 2 , into the lightest inert state, H 1 , and a virtual photon, which then would split into a light ff pair 6 . Figure 2: Radiative decay of the heavy neutral particle The corresponding loops go through triangle and bubble diagrams with H ± i and W ± entering, see Figs 3-4. Note that there are also box diagrams which contribute to the process H 2 → H 1 ff , presented in Fig. 5. Here, the ff pair is produced through the SM gauge-fermion tree-level vertices, without producing an intermediate off-shell photon. The corresponding topologies also see the contribution of inert, both charged and neutral (pseudo)scalars. However, due to the mass suppression, the contribution from the box diagrams is small, of order 10%, and it leaves the results practically unaffected. For reasons of optimisation then, we do not show the results of these box diagrams in the numerical scans and we may refer to this one-loop process as a radiative decay.
Before moving on to study the latter, we would like to stress at this point that one could attempt constructing analogous diagrams to those in Figs. 3-4 with H 2 replaced by A 1 or A 2 , leading to A i → H 1 γ * , i = 1, 2. Notice, however, that this decay would lead to a CPV process, while the model we analyse here is explicitly CPC. Indeed, further notice that spin conservation requires that it is only the scalar polarisation of the virtual photon that contributes to the H 2 → H 1 γ * transition. To check the correctness of the calculations we have explicitly verified this to be the case, as there are cancellations between diagrams that lead to the amplitude being equal to zero otherwise, as discussed in Sect. 4. Also note that the process A i → H 1 Z * does exist at tree-level in both the I(2+1)HDM (for i = 1, 2) and I(1+1)HDM (for i = 1) and contributes to the E T ff signature, as discussed previously. However, in the interesting regions of the parameter space where the invariant mass of the ff pair is small, i.e., << m Z , this process is sub-dominant.
In short, the only (effective) loop-level decay to consider is and this does not exist in the I(1+1)HDM, as CP-conservation prevents the only possibly similar radiative decay in its inert sector (i.e., A 1 → H 1 γ * ). Therefore, as intimated, this signature can be used to distinguish between the I(1+1)HDM and models with extended inert sectors, such as the I(2+1)HDM. Figure 3: Triangle diagrams contributing to the H 2 → H 1 γ * decay, where the lightest inert is absolutely stable and hence invisible, while γ * is a virtual photon that couples to fermion-antifermion pairs. Analogous diagrams cannot be constructed if the initial particle is A 1 or A 2 . Figure 4: Bubble diagrams contributing to the H 2 → H 1 γ * decay, where the lightest inert particle is absolutely stable and hence invisible, while γ * is a virtual photon that couples to fermion-antifermion pairs. Analogous diagrams cannot be constructed if the initial particle is A 1 or A 2 .

The E T ff signature at the LHC
In this subsection, we focus on the possible sources of the aforementioned specific signature that can arise in the I(2+1)HDM, namely, missing transverse energy and a fermionantifermion pair, E T ff . This final state can be produced both at tree-level and through one-loop decays, as previously explained. We dwell further on this here. The first mechanism is related to decays of the SM-like Higgs particle which is produced, e.g., through ggF. The hgg effective vertex is identical to that in the SM, as the gauge and fermionic sectors in the I(2+1)HDM are not modified with respect to the SM. The Higgs particle can then decay into a pair of neutral or charged inert particles, denoted in Fig. 6 by S i,j . Depending on the masses of S i,j , these particles can further decay, providing various final states. Figure 6: The ggF-induced production of the SM-like Higgs particle at the LHC with its decay into inert particles, denoted as S i and S j .
In the CPC I(2+1)HDM, a process contributing to the E T ff signature (and one of our signals) is where the off-shell γ * splits into ff and the H 1 states escape detection 7 .
Notice that there is also a tree-level h decay into two charged scalars with the same signature ( E T ff ), albeit not an identical final state (the two would remain indistinguishable though), following the pattern: where the neutrinos escape detection as (additional) E T . The process in (28) is loop mediated and depends on g H 1 H 2 h , a coupling affecting also DM relic density. Therefore, if this coupling is small, the whole process is suppressed. However, we shall maximise this coupling, while maintaining consistency with DM constraints. We also assume a mass spectrum so that the charged Higgs masses entering the loops are not too heavy, since their large masses would also suppress the loop. In fact, we shall see that there can be parameter configurations for which m H 1 + m H 2 < m h , so that SM-like Higgs production and (loop) decay is resonant, thereby benefiting of an enhancement of O(1/α EM ). The process in (29) is a tree-level one, therefore potentially competitive. However, for the parameter space of interest, maximising the yield of the loop process, this mode becomes negligible, for two reasons: on the one hand, the charged Higgs masses are generally heavy so that there can be no resonant h involved while, on the other hand, the g H ± i H ± i h coupling is generally small. In principle, there is another tree-level signal inducing the E T ff final state in our scenario, see diagrams (A) and (B) in Fig. 7, induced by quark-antiquark annihilation and proceeding via an s-channel off-shell (primary) Z * , wherein the on-shell (secondary) Z eventually decays into an ff pair. However, this is of no concern here. The reason is twofold. On the one hand, as explained, the region of parameter space over which process (28) is interesting for LHC phenomenology is the one where the g H 1 H 2 h strength is maximal and h is possibly resonant: this is when the DM relic density sees a large contribution from H 1 H 2 co-annihilation processes 8 , which in turn means that large g H 1 H 1 h (possibly in presence of a resonant h) and g H 1 H 1 ZZ couplings are forbidden by such data, so that process (30) becomes uninteresting at the LHC. On the other hand, in our construct, process (30) is nothing more than a subleading contribution to the invisible Higgs signature of the SM-like Higgs boson (dominated by ggF and VBF topologies, extensively studied already in Ref. [13]), rather featureless, in fact, as it does not catch any of the heavy scalar states of the model, unlike reaction (28), which is sensitive to all of them, so that one could study the kinematic distributions of the final state attempting to extract their masses by isolating the corresponding thresholds entering the loops 9 . For these reasons, we will not discuss these two topologies any further. Another way of obtaining exactly the H 1 H 1 ff final state is shown in graph (C) of Fig. 7, again produced through s-channel quark-antiquark annihilation into a virtual neutral massive gauge boson, i.e., wherein the DM candidate is produced in association with a pseudoscalar state and the Z may be off-shell. This mode is indeed competitive with the one in (28) over the region of I(2+1)HDM parameter space of interest, so we will extensively dwell with it numerically in the remainder of the paper. Further, diagram (C) in Fig. 7, unlike graphs (A) and (B) herein, because of its heavy pseudoscalar components, may also be isolated in the aforementioned kinematic analysis. Finally, we conclude this subsection by listing, in Fig. 8 (prior to the H 2 → H 1 ff decay), the topologies entering VBF production contributing to the E T ff final state (our second signal) via where q i,j,k,l represents a(n) (anti)quark of any possible flavour (except a top quark).
Here, two aspects are worth noticing. Firstly, there is the additional presence of two forward/backward jets, which may or may not be tagged (we will treat them inclusively). Secondly, not all diagrams proceed via h → H 1 H 2 induced topologies, graph (A), hence unlike the case of ggF, since graphs (B) and (C) are also possible. Clearly, the first diagram dominates when h can resonate while the last two become competitive otherwise.
We shall see how ggF and VBF will compete over the I(2+1)HDM parameter space of interest in being the carrier of its hallmark signature E T ff in a later section.

Calculation
In this section, we discuss the details of our calculation. In fact, the case of the channel in (31) is easily dealt with, as this is a tree-level process, which we computed numerically using CalcHEP [20]. The bulk of our effort was concentrated upon the loop processes (28) and (32), which we have tackled in factorised form, i.e., by breaking up the two channels into pp → H 1 H 2 X production followed by the H 2 → H 1 ff decay. Here, the ggF and VBF topologies entering at production level are well known in the literature, so we do not discuss them (again, we computed these numerically by exploiting CalcHEP). We therefore address in some detail only the case of the loop decay. This is expressed through a tensor structure appropriate to the I(2+1)HDM particle spectrum and illustrated for the case f = e, so that we can safely take m e = 0 10 . In general, there are two types of one-loop diagrams that contribute to the process namely, those embedding the one-loop effective vertex H 2 H 1 γ * , given by the diagrams in Figs. 3-4 plus the box diagrams shown in Fig. 5. Here, the labels p i and k j identify the external scalar and fermion momenta, respectively. In the following, we use the unitary gauge.
The calculation below is done for the pair of CP-even dark particles H 2 and H 1 , however, all results hold for CP-odd neutral dark particles as well, i.e., A 2 and A 1 , following simple replacements of masses, m H i → m A i , and relevant vertex coefficients, The general expression for the amplitude of the loop calculation is: where is the general structure of the vertex H 1 H 2 γ * obtained in the calculation at one-loop level. However, when we consider the term (p 3 − p 2 ) ν and contract it with γ µ g µν we have: Then the Dirac equation in the limit of m e = 0 gives us: Under these circumstances, we can take which is the same as if the γ were on-shell in the process Therefore, the general structure of the amplitude is: where A(p 3 + p 2 ) µ is related to the contribution of the each diagram in Figs. 3, 4 and 5: where i runs across all diagrams.

Individual contributions to H 2 → H 1 ff
There are six of these, five for the case of the triangle and bubble diagrams of Figs. 3-4 plus two cumulative ones for the box diagrams shown in Fig. 5 11 .
• The first contribution, M µ , comes from a diagram with two charged scalars H ± i (i = 1, 2) and one W ± in the loop, given by diagram (A) in Fig. 3: where 11 Ultraviolet renormalisation is implictly performed for the former.
and m H ± i (i = 1, 2) are the masses of the charged scalars, m H 1 is the mass of the DM candidate and m H 2 is the mass of the next-to-lightest inert particle H 2 . The A ± i s are coefficients related to the vertex structure of the loop diagram whose details are presented in Sect. 4.2. We define m 2 12 = (p 3 −p 2 ) 2 = (k 1 +k 2 ) 2 = 2k 1 ·k 2 , considering the limit m e = 0. Using this tensorial structure we calculate the other diagrams.
• The tensorial amplitude for the diagram with two W ± and one charged scalar H ± i in the loop, given by diagram (B) in the Fig. 3, is: • For the diagram with one H ± i and one W ± particle in the loop, which is (A) in Fig. 4, the tensorial amplitude is • For the diagram with two scalars in the loop, i.e., (B) in Fig. 4, the tensorial amplitude is: However, this last equation is zero because it is an odd function.
• For diagram (C) in Fig. 4, with one H ± i and one W ± in the loop, the tensorial amplitude is given by: • For the box diagrams with W ± in the loop (graphs (B) and (D) in Fig. 5) we obtain: where The structure of M Z−box , i.e., for diagrams with the Z instead of a W ± (graphs (A) and (C) in Fig. 5), is similar, with the replacements (1 − γ 5 ) → (C V − C A γ 5 ) and m W → m Z . When one considers the crossed box diagrams, the ultraviolet divergences cancel. In practice, when performing the loop calculation, one can see that the contribution of the boxes is not important due to the mass suppression and contributes to the aforementioned about 10% of the overall calculation. Hence, as intimated, we shall neglect this from now on.

Role of the A ± i s
The coefficients A ± i s, related to the vertex structure of loop diagrams, are the characteristic features of the model. They are sensitive to the CP properties of the decaying particles and they can provide us with the information necessary to cancel the ultraviolet divergences.
For the three neutral scalars we define: where θ h,a,c are the inert mixing angles defined in Sect. 2.2. We use the shorthand A ± i for A ± H ± i ,S where S could be any of the neutral scalars H 2 , A 1 , A 2 . The following relations hold: Despite not being exploited phenomenologically in the remainder of the paper, for completeness, we also describe here the case A 1,2 → H 1 γ * → H 1 e + e − . In the CP conserving I(2+1)HDM, one can distinguish the CP-even inert scalar and CP-odd inert scalar in the diagrams of the Figs. 3 and 4. When considering the amplitude of any diagram plus its crossed companion, one obtains the following results: and as a consequence which is consistent with the observation we made before: CP conservation requires A 1,2 → H 1 γ * → H 1 e + e − to be zero. However, for the box diagrams associated with Fig.  5, A 1,2 → H 1 e + e − decays are possible but their contributions are small. In fact, these decays could also be mediated at one-loop level by an on-or off-shell Z boson, however, the tree-level mode A 1,2 → H 1 Z * → H 1 e + e − (already discussed) is much larger, which is why we concerned ourselves with the latter and not the former. Finally, one can see from (57) that A ± 1 = −A ± 2 , which is crucial for the cancellation of the ultraviolet divergences. In fact, the total contribution of the one-loop calculation is, taking account (61) and (56): Now, taking into account (57), we have One can see then that ultraviolet divergences cancel perfectly.

Partial decay width of H 2 → H 1 ff
When evaluating the tensorial integrals of (41)-(47), these expressions are reduced in terms of Passarino-Veltman scalar functions: where F PV (m H ± i , m W , m 2 12 , m H 1 , m H j ) is given in Appendix A. Then, comparing (38), (39) and (64), we calculate the factor A: One can see that A is a function of the same variables of F PV and the factor A + 1 . Besides, following the notation of [21] for three-body decays, in addition to the variable m 2 12 defined previously, we also introduce m 2 i3 = (k i + p 2 ) 2 = 2k i · p 2 + m 2 H 1 (i = 1, 2). Taking this into account, one can obtain the square amplitude (38) of the loop process (upon the usual final state spin summation): Besides, it is convenient to define Then (dropping henceforth the arguments of F PV ) one has In agreement with Ref. [21], the partial decay width of H 2 → H 1 e − e + is: From (64) and (70), one can observe that the one-loop function F PV contains only the integration variable m 2 12 , so that we can integrate firstly in the variable m 2 23 in the following way: where the integral for m 2 23 is possible to obtain analytically, as With this result we can do the numerical calculation using the LoopTools library [22].

Effective Lagrangian
As it was suggested some years ago [23,24], one can perform a general study of the discussed radiative process in a model independent way using the effective Lagrangian technique, which can parameterise the virtual effects of new physics of a given model. This approach is mandatory in our case, as we will be implementing the effective H 2 H 1 e + e − vertex in CalcHEP, which is otherwise unable to perform the calculation efficiently if using the exact formulae from the previous subsection. The effective Lagrangian for the I(2+1)HDM will be an extension of the SM one [25], following a similar parameterisation to the one used for the case of the 2HDM, when rare decays of neutral CP-odd [23] and charged Higgs [24] bosons were implemented in this way. Following these studies, we use SU (2) × U (1) gauge invariant operators of higher dimension similar to those given in [24]. Adopting this approach, we can define operators that satisfy all symmetries imposed in our model, in particular the discrete symmetry Z 2 . Then the corresponding effective Lagrangian for our model is: where L I(2+1)HDM is the I(2+1)HDM Lagrangian, Λ is the scale of new physics, the O i n s are the higher dimensional operators and the unknown c i n parameters are their dimensionless Wilson coefficients, whose order of magnitude can be estimated since gauge invariance makes it possible to take into account the order of perturbation theory where each operator can be generated in the fundamental theory [26]. This fact allows us to introduce a hierarchy among operators, e.g., when the operators are generated at one-loop level, they must be suppressed by the loop factor (4π) −2 . Using this method, we can study the generic structure of any process.
With the knowledge that the box diagrams and the tree-level diagrams with the offshell Z are sub-dominant, we consider the effective coefficient of the vertex H 2 H 1 e + e − . In practice, we can implement such a vertex in the effective Lagrangian as follows: where Λ ≥ v and c i can be estimated given the order of perturbation theory [25]. In our model, for the full process H 2 → H 1 e + e − , we must consider the following for the coefficient c i : (i) as the process is generated at one-loop level, it must be suppressed by the loop factor (4π) −2 ; (ii) the order in the perturbation theory is proportional to e 2 g 2 , (see (70)). A good approximation is, therefore, c 1 ∝ e 2 g 2 /(4π) 2 . The first and second operators then induce the structure of the loop calculation stemming from the diagrams of Figs. 3-4 while the following operators relate to the structure of the diagrams given in Fig. 5. Given the effective Lagrangian, we can induce the effective vertex H 2 H 1 e + e − as: In this framework, the Wilson coefficient c 1 contains information of the parameters of the Higgs potential of the model, in particular of the mixing angle of the charged sector, which is consistent with the amplitude of loop calculations (see (57), (64)). The Wilson coefficient to this order does not depend on the variables m 2 12 and m 2 23 of (71), and in principle c 1 behaves like a constant in the eyes of these integration variables. In particular, is a function of the charged Higgs masses and their mixing angle [24,27], and the scale of new physics Λ could in general be of order 1 TeV or the energy necessary at the LHC experiments to detect the DM candidate. Now, we can define the effective coefficient K as which we have implemented in CalcHEP as an effective vertex H 2 H 1 e + e − in the following way: In order to relate the K-factor with all numerical results of the previous sections and taking into account the discusion of the Wilson coefficients, we calculate the amplitude of the process H 2 → H 1 e − e + using the effective vertex in (79), which is given by so that the amplitude squared is Thus, the partial decay rate of the H 2 → H 1 e − e + channel, in terms of the K-factor, is where I 3 is given by The K-factor is therefore given by where the width Γ(H 2 → H 1 e + e − ) is calculated using LoopTools. Thus, the K-factor is related directly to the loop calculation through (85). Using this method, we are able to use CalcHEP, since we no longer need to perform any integration externally to the generator itself, as required by the fully fledged computation performed in the previous subsection, thereby by-passing the fact that CalcHEP is actually a tree-level generator.
In order to complete the study of the decay process H 2 → H 1 e + e − , it is necessary to compare the e + e − mode with the others possible final states, H 2 → H 1 ff where f = u, d, c, s, b, µ, τ . Given the effective Lagrangian in (76), one can obtain the contribution of the fermions via the following operators: in agreement with the general structure of the loop calculation, giving the general expression to be M = iev(k 1 ) A( / p 3 + / p 2 ) + (B + C( / p 3 + / p 2 ))P L + (D + E( / p 3 + / p 2 ))P R u(k 2 ), (87) where A, B, C, D and E are form factors associated with the loop structure. This structure helps us to calculate all the aforementioned form factors, taking in account all the contributions of the boxes (factors B, C, D and E) and triangles (factor A given in (67)) in the loops. One can then calculate each form factor separately as they are all individually convergent. Finally, notice that in the channel H 2 → H 1 bb the mass of the top quark appears in the boxes: while this makes the calculation more cumbersome, the mass effects do not contribute significantly to the yield of the total rate. Besides, in the approximation m e = 0, the factors B, C, D are zero and the factor E is small. In appendix B we show the complete expressions of the factors associated with the box diagrams of Fig 5.

Results
The benchmark scenarios that we study here do not necessarily correspond to regions of the parameter space where our DM candidate accounts for all the observed relic density in agreement with Planck data. In fact, the aim of these benchmark scenarios is to show in which regions of the parameter space the model has a discovery potential at the LHC. Following the discussion in Sect. 2.3, we define three base benchmark scenarios, A50, I5 and I10 in the low DM mass region (m H 1 ≤ 90 GeV) as shown in Tab. 1.
The main distinguishing parameter here is the mass splitting between H 1 and the other CP-even scalar, H 2 . Benchmark A50 (m H 2 − m H 1 = 50 GeV) is taken from the analysis done in [10]. Relatively large mass splittings between H 1 and other neutral scalars leads to a standard DM annihilation in the Universe, providing us with a DM candidate which is in agreement with DM searches for a large part of the parameter space. However, we expect the tree-level decays to dominate over the loop signal through Benchmarks I5 (m H 2 − m H 1 = 5 GeV) and I10 (m H 2 − m H 1 = 10 GeV) have an intermediate mass splitting between H 1 and H 2 of the order of a few GeV. As mentioned in section 2.4, this influences the thermal history of DM, due to the appearance of coannihilation channels.
For the I benchmarks, we expect the tree-level decays to be reduced, since there is a small mass gap between H 1 and A 1,2 . Therefore, the intermediate gauge boson is produced off-shell. Further decreasing of the H 1 -H 2 and H 1 -A 1,2 mass splittings 12 , leads to strengthening of the desired loop signal, with further reduction of all tree-level decays. Note, however, that with increasing the mass splitting, the loop process acquires more phase space and starts seeing the Z * → ll contribution and the partial width grows as a result.
In all cases, differences between m H 1 and masses of both charged scalars are relatively large. This leads to important consequences for the thermal history of DM particles: charged scalars are short-lived and they will not take part in the freeze-out process of H 1 . However, this mass difference is not big enough to suppress the studied loop processes. Increasing this mass difference would lead to a smaller cross-section and, therefore, worse detection prospects. We would also like to stress that the all chosen mass splittings are in agreement with EWPT constraints, which disfavour a significant discrepancy between masses of charged and neutral particles. On the other hand, a significant reduction of this mass splitting would increase the coannihilation effect in the Universe, hence leading to heavily reduced relic density, and thus disfavouring the 3HDM as the model for Dark Matter. 50  75  125  75  125  I5  5  10  15  90  95  I10 10 20 30 90 100 Table 1: Definition of benchmark scenarios with the mass splittings shown in GeV.
Figs. 9-11 show the anatomy of the given scenarios, which include not only the cross sections for leptonic ( E T l + l − ) and hadronic ( E T qq) final states, but also the relevant couplings in each case with the same colour coding. The Higgs-DM coupling is also shown for reference.
For each benchmark scenario, we calculate the cross section for three processes, namely, the ggF process (28), the tree-level process (31) and the VBF process (32) and present the dominant couplings entering in each case. Let us first focus on scenario A50 presented in Fig. 9, which has two special features. First, mass splittings between H 1 and other inert particles are relatively large, as well as the main couplings (in particular the g ZH 1 A 1 ), which leads to large tree-level Z-mediated cross sections (the blue curve). Second, the Higgs-DM coupling, g hH 1 H 1 , is chosen such that the relic density is in exact agreement with Planck measurements. To fulfil that, around the Higgs resonance the coupling needs to be very small, of the order of 10 −4 [10]. As the g hH 1 H 2 coupling is closely related to g hH 1 H 1 , we observe a sudden dip for the orange curve (g hH 1 H 2 ), which then leads to a reduced cross section for the ggF processes, driven by that particular coupling. We also observe that the cross section for the VBF processes, which depend mainly on large mass splittings and relatively constant gauge couplings, are as expected relatively constant for this benchmark. Scenario I5, shown in Fig. 10, differs from the scenario A50 above. Here, the mass splittings are much smaller, but also the Higgs-DM coupling is set to a constant value for all masses, as seen in Fig. 10. This makes the phase space structure more visible. For m H 1 < m h /2 all cross sections are roughly constant, with the ggF processes enhanced through the resonant Higgs production. However, after crossing the Higgs resonance region, with no increase of the Higgs-DM coupling to compensate for that, we observe a rapid decrease of the value of the cross section. For larger masses the cross section are too small to be observed for the current LHC luminosity.  Very similar behaviour is present for scenario I10 depicted in Fig. 11, where, similarly to scenario I5, the Higgs-DM coupling is set to a constant value for all masses. Again we observe the almost constant cross sections, which are rapidly reduced after we cross the Higgs threshold.   In Tabs. 2-4, we show the BR of H 2 → H 1 ff for any ff pair, whose production the m h 2 − m H 1 mass splitting allows for, for an exemplary value of m DM = 54 GeV. For each decay channel, we also show the cross section value (in pb units) for all three discussed processes, the tree-level background as well as the ggF and VBF cross-sections for Higgs h production times the respective branching ratio for h → H 1 H 2 → H 1 H 1 ff . The cross-section for h production and decay into two charged scalars for m DM = 54 GeV is very small as shown in Tab

Conclusion and outlook
In this paper, we have assessed the sensitivity of the LHC to Higgs signals in the E T ff channel, f = u, d, c, s, b, e, µ, τ , with invariant mass of the ff pair much smaller than the Z mass. This signature would in fact point towards an underlying 3HDM structure of the Higgs sector, with one active and two inert doublets (so that the scenario can evocatively be nicknamed as I(2+1)HDM), induced by the decay H 2 → H 1 ff , where H 1 represents the lightest CP-even neutral Higgs state from the inert sector (thereby being a DM candidate) and H 2 the next-to-lightest one. The decay proceeds via loop diagrams induced by the propagation of both SM weak gauge bosons (W ± and Z) and inert Higgs states (H ± 1,2 and A 1,2 ) in two-, three-and four-point topologies, wherein the leading contribution comes from the intermediate decay step H 2 → H 1 γ * , involving a very low mass virtual photon scalarly polarised, eventually splitting in a collimated ff pair, which would be a distinctive signature of this Higgs construct. In fact, the corresponding 2HDM version, with one inert doublet only, i.e., the I(1+1)HDM, contains only one CP-even and only one CP-odd neutral Higgs state, so that no such a decay is possible owing to CP conservation.
This signature would emerge from SM-like Higgs boson production, most copiously via ggF and VBF, followed by a primary h → H 2 H 1 decay, so that the complete particle final state is H 1 H 1 ff , wherein the two DM candidates would produce missing transverse energy, accompanied by some hadronic activity in the forward and backward directions, originating by initial state gluon radiation or (anti)quark remnant jets, respectively, for ggF and VBF. In fact, amongst the possible fermionic flavours f , the cleanest signature is afforded by the leptonic ones (f = l), in view of the overwhelming QCD background. While the muon and tauon cases are the cleanest, the latter being larger than the former (assuming only leptonic decays of the τ 's), the electron case is potentially the one giving raise to the most spectacular signal, which, owing to parton distribution imbalances, so that the h state would be boosted, would appear at detector level as a single EM shower with substantial E T surrounding it.
However, there is a substantial tree-level contribution, due to qq → Z * H 1 H 1 topologies (a first one involving single h-strahlung followed by h → H 1 H 1 splitting, a second one via a Z * Z * H 1 H 1 vertex and a third one through A 1,2 H 1 production followed by A 1,2 → H 1 Z * decay), which is potentially much larger than the aforementioned loop diagrams, thereby acting as an intrinsic background. In fact, even though the Z * ought to be significantly off-shell in its transition to ff pairs to mimic the γ * → ff splitting, this can happen with substantial rates, because of the rather large value of the total Z decay width. It is therefore clear that the H 2 → H 1 ff signal can only be established in presence of a rather small mass gap between H 2 and H 1 . To this effect, we have then defined a few benchmarks on the I(2+1)HDM parameter space where the mass difference m H 2 − m H 1 is taken to be increasingly small, varying from 50, to 10 to 5 GeV. Correspondingly, we have seen the relevance of the loop processes growing with respect to the tree-level one, with ggF dominating VBF, to the point that the former become comparable to the latter for cross sections and BRs directly testable at Run 2 and/or Run 3 of the LHC. This is particularly true over the DM mass region observable at the CERN machine, i.e., for small values of the DM candidate mass, typically less than m h /2. In this case, the cumulative signal can be almost within an order of magnitude or so of such an intrinsic background.
Furthermore, other (irreducible) background processes can be present. The first one is the tree-level h decay into two charged scalars with the same signature ( E T ff ), albeit containing two (invisible) additional neutrinos, which has a very small cross section, as shown in Tab. 5 for each of our benchmarks for the usual illustrative value of m DM = 54 GeV. A second one is due to gg → h → V V (via resonant h production) and qq → V V (gauge boson pair production), where V V = W + W − or ZZ. These two subprocesses have inclusively very large cross sections, of O(10 pb) (prior to V decays), compared to our signals, and a significant amount of (differential) kinematical selection ought to be employed to reduce these noises, which is clearly beyond the scope of this paper. However, a few handles can be clearly exploited. For the case V = Z, a veto m ll = m Z can always be adopted. For the case V = W ± , a requirement of the kind m ll << m W , combined with the request of identical lepton flavours, can be used.
We have obtained these results in the presence of up-to-date theoretical and experimental constraints, including amongst the latter those from colliders, DM searches and cosmological relic density. Therefore, we believe that the advocated discovery channel might serve as smoking-gun (collider) signature of the I(2+1)HDM, that may enable one to distinguish it from the I(1+1)HDM case, in a few years to come. In fact, once this signal is established and some knowledge of the H 2 and H 1 masses gained, the latter can be used to extract additional manifestations of the prevalent H 2 → H 1 γ * decay, by considering the selection of additional splittings γ * → ff , where f can be identified with q = u, d, c, s, b, depending upon the relative value of m H 2 − m H 1 and 2m f . Finally, in reaching these conclusions, we emphasise that we have done a complete one-loop calculation of the H 2 → H 1 ff decay process, including all topologies entering through the same perturbative order, i.e., not only those proceeding via H 2 → H 1 γ * → H 1 ff , which was never attempted before, so that we have collected the relevant formulae in this paper for future use.
In conclusion, the 3HDM with two inert doublets, provides a well motivated dark matter model with distinctive LHC signatures in certain regions of parameter space arising from novel Higgs decays, the most spectacular being e + e − + E T mono-shower. A The Passarino-Veltman functions, F PV In this appendix, we present detailed formulae for the functions F PV in a useful shorthand notation. In order to describe the cancellations of the ultraviolet divergences we first define the differences of scalar functions B 0 s and C 0 s in the following way 13 : We also define the following parameters: With these definitions, we can write the function F PV as:

B The functions of the box diagrams
Firstly, we define the following Passarino-Veltman functions which can obtained in the FeynCalc package [29]: (C A + C V ) 2 /2 inside B, C, D and E defined below.