Minimal Consistent Dark Matter models for systematic experimental characterisation: Fermion Dark Matter

The search for a Dark Matter particle is the new grail and hard-sought nirvana of the particle physics community. From the theoretical side, it is the main challenge to provide a consistent and model-independent tool for comparing the bounds and reach of the diverse experiments. We propose a first complete classification of minimal consistent Dark Matter models, which provides the missing link between experiments and top-down models. Consistency is achieved by imposing renormalisability and invariance under the full Standard Model symmetries. We apply this paradigm to fermionic Dark multiplets with up to one mediator. We also reconsider the one-loop contributions to direct detection, including the relevant effect of (small) mass splits in the Dark multiplet. Our work highlights the presence of unexplored viable models, and paves the way for the ultimate systematic hunt for the Dark Matter particle.

1 Introduction candidate comes in a multiplet of the weak Isospin SU(2), its charged partners may play an important role in the LHC phenomenology often being more important than the neutral state itself. This is the case for charginos in supersymmetry. In addition, simplified models often violate gauge invariance at high scales [41], which is a crucial principle for building a consistent BSM model that incorporates the SM together with new physics. For example, considering simplified models with a new heavy gauge vector boson mediating DM interactions, one should also introduce a mechanism responsible for generating the mediator mass and ensuring gauge invariance for the model [41]. Eventually, this necessarily requires introducing an additional sector into the model that may affect the DM phenomenology [41,43]. These drawbacks strongly indicate the next step in the evolution of the DM investigation, based on building Minimal Consistent Dark Matter (MCDM) models. MCDM models can be still understood as toy models that, however, take in full account the consistency with the symmetries of the SM. In our approach, MCDM models consists of one DM multiplet and at most one mediator multiplet. Furthermore, a particular MCDM model can be easily incorporated into a bigger, more complete and fundamental, BSM model and be explored via complementary constraints from collider and direct/indirect DM search experiments as well as relic density constraints. The exploration of complementarity of the collider and non-collider constraints within the complete models such as MCDM ones is very appealing especially now as we have a large amount of data from the LHC. Combining searches may shed light on the BSM physics in the form of DM, which can be near the corner of the combined collider and non-collider limits. Another attractive feature of the MCDM approach is the minimal but self-consistent parameter space that can be potentially mapped to the parameter space of known (and completely new) BSM models.
Many implementations of MCDM models have been studied in the literature [19-21, 27, 44-46], however there has been no attempt on their systematic classification yet. This is precisely the aim of the present work. In this study we shall: a) perform a complete classification of MCDM models, with at most one mediator and including only renormalisable interactions (with some notable exceptions); b) present the main features for each class of MCDMs constructed using the main building principles we state below.
We believe that this classification, and the MCDM approach in general, will create a solid framework for the consistent exploration of DM models at collider and non-collider experiments for the complementary probe of Dark sectors. The paper is organised as follows: after articulating the main principles behind the MCDM approach in Section 2, we summarise the main properties of models with only a DM candidate in Section 3. Here we also present a detailed calculation of the one-loop cross section for direct detection, which includes for the first time the mass split between components of the DM electroweak multiplet. In Section 4 we classify and characterise models with a single mediator. Finally, in Section 5 we study in detail a new model that emerges from the classification, featuring a Dirac fermionic DM candidate and a CP-odd scalar mediator. In some regions of the parameter space, the scalar mediator can be accidentally stable and contribute to the relic density. We offer our conclusions and outlook in Section 6.
to Refs [47][48][49][50]. The electroweak quantum numbers will be encoded in the weak Isospin, I, and the hypercharge, Y , of the multiplet. Furthermore, we will denote with a tilde the multiplets that belong to the Dark sector, i.e. they cannot decay into purely SM final states. The multiplets we consider, therefore, read: As some mediator multiplets may carry QCD quantum numbers, we will use a superscript c to label this feature.
To construct consistent minimal models, we follow these main building principles: A) We add one Dark multiplet (including the singlet case) and all its renormalisable interactions to SM fields, excluding the ones that trigger the decays of the multiplet, which is therefore stable by construction. The models will automatically include a Dark symmetry, being Z 2 or U (1) depending on the multiplet. The weak Isospin and hypercharge are constrained by the need for a neutral component, therefore we will have the following two cases: -for integer isospin I = n, n ∈ N, then Y = 0, 1 . . . n; -for semi-integer isospin I = (2n + 1)/2, n ∈ N, then Y = 1/2, 3/2 . . . (2n + 1)/2.
Note that the case of negative hypercharge can be obtained by considering the charge conjugate field, thus the sign of Y is effectively redundant, and we will consider Y ≥ 0.
B) We consider models where only one Dark multiplet is present, and mediators are SM fields. 2 While our principle is to be limited to renormalisable interactions, under the assumption that higher order ones are suppressed by a large enough scale to make them irrelevant for the DM properties, in some cases we will consider dimension-5 operators.
C) In additional to point B), we consider adding just one mediator multiplet, characterised by the respective weak Isospin, I , and hypercharge, Y . The mediator multiplet can be odd or even with respect to the Dark symmetry, and its quantum numbers are limited to cases where renormalisable couplings to the Dark multiplet and to the SM are allowed. This leaves open the possibility of multiplets carrying QCD charges, which we label with a superscript c . The mediators are labeled as following: The odd mediator multiplets can also contain a DM candidate if a neutral component is present.
D) We consider all renormalisable interactions allowed by the symmetries of quantum field theory. Our basic assumption for MCDM models is that higher-order operators are suppressed by a scale high enough that the LHC is unable to resolve the physics generating the operators. The effect on the DM properties is also considered negligible (except for dim-5 operators generating mass splits). E) We ensure cancellation of triangle anomalies, so that the MCDM models entails consistent gauge symmetries.
2 Note that this model building approach has been used in [45] to construct models of so-called Minimal Dark Matter, so some of the results we present here can be found in this reference. However, our approach has some differences: in Ref. [45], the symmetry making the DM candidate stable or long lived emerged at low energy, at the level of renormalisable interactions, while decays could be induced by higher dimensional couplings to the Higgs multiplets. In our case, we assume that a parity or global U (1) symmetry is also respected by higher dimensional operators. Henceforth, we do not take into account any constrains on the isospin of the multiple.  With the notations above, following the precepts A) to E), we can classify all MCDM models with up to one mediator multiplet using a 2-dimensional grid in Spin(DM)-Spin(mediator) space, as presented in Table 1. Each specific DM model is denoted by a one-or two-symbol notation, indicating the DM multiplet first, followed by the mediator multiplet. In general, the interactions of the DM candidate to the SM are mediated by SM particles (e.g. by the EW gauge bosons and the Higgs) and other components of the DM multiplet, besides the components of the mediator multiplet. Hence, highly nontrivial interference effects can arise. Furthermore, some couplings entail flavour structure, which need care as they may incur very strong bounds. Eventually, the case with no mediator multiplet is denoted by just one symbol labelling the DM multiplet. In this case the role of mediators can only be played by SM particles and members of the DM multiplet.
In the remainder of this paper, we will focus on spin-1/2 DM multiplets, leaving the other two cases for a future publication.
3 Case of one DM multiplet:F I Y andM I 0 models Models where the DM belongs to a single EW multiplet, while no other light states are present, have been studied in great detail, starting from the seminal paper in Ref. [45]. In this section we briefly review the main properties of these minimal models, and add a detailed discussion of the following novel aspects: i) We provide an improved formula for the mass split induced by EW loops, which is numerically more stable than the one given in Ref. [45].
ii) We discuss in great detail the effect of couplings to the Higgs boson arising as dimension-5 operators. While going beyond renormalisability principles, they are generated by integrating out a single mediator (thus, they can be considered as a limiting case from some of the models discussed in Section 4). Furthermore, a class of these operators have special phenomenological relevance as they help salvage some of the minimal models with non-zero hypercharge.
iii) We provide a detailed and up-to date discussion of direct detection bounds at one-loop level. We include for the first time the effect of mass splits within the DM multiplet, and show their relevance. iv) We discuss the impact of nuclear uncertainties and of the variation of the gluon contribution due to the mass splits. Both generate comparable uncertainties in the total spin-independent cross sections, which emerge as an uncertainty in the DM mass limits of hundreds of GeVs.
This section also serves to fix the notation we will adopt in the rest of the paper. When writing Lagrangians and interactions we will consistently use Ψ = Ψ L + Ψ R for Dirac DM multiplet, with Ψ R = Ψ C L for the Majorana case (where C indicates the charge conjugate field), ψ i for the components of a Dirac multiplet and χ i for the components of a Majorana multiplet. Furthermore, we only consider Y ≥ 0, as the case of negative hypercharge is straightforwardly analogous to the corresponding positive value case. We will use M DM to denote the mass of the neutral component that serves as DM candidate. In the "stand alone" case, only gauge interactions of the EW gauge bosons, W ± , Z and photon, are allowed at renormalisable level. This simple class of models has well established properties [45], which we list below: -A gauge coupling g Zψ 0 ψ 0 is always present for Dirac multiplets with Y = 0, which are thus excluded by direct detection even for under-abundant points (for M DM < m Z /2 the invisible width of the Z also excludes the model). On the contrary, when Y = 0, the coupling g Zψ 0 ψ 0 always vanishes.
-Due to the absence of couplings to the Higgs field, the mass split between the neutral and charged components of the DM multiplet are generated by EW loop corrections and are always small (below a few hundred MeVs, with the precise values depending on the hypercharge of the multiplets). This leads to long lived particles, especially at high mass. The lightest component is not always guaranteed to be neutral: this only occurs for multiplets with Y = 0 and maximal hypercharge, Y = I.
-For Y ≥ 1 and isospin I = Y (hence, I ≥ 2), the mass range with the neutral component being the lightest is excluded by the Z width. Hence, these multiplets in isolation cannot provide a DM candidate.
-For Y = 1/2 and I ≥ 3/2, the lightest component is neutral for M DM 570 GeV. Above this threshold, the charge −1 state becomes the lightest in absence of Higgs couplings.
-For Y = 1/2, a dim-5 operator with the Higgs boson generates a mass that splits the neutral component in two Majorana mass eigenstates (pseudo-Dirac case). This salvages the models from exclusion via the Z interactions.
-Taking into account loop-induced mass splits, the loop-induced cross sections ensures that current and future direct detection experiments can probe multiplets with I ≥ 1, where I ≥ 2 can be completely ruled out, while the case of a doublet I = 1/2 is always below detection. Uncertainties in the nuclear form factors and mass splits for the gluon contribution generate uncertainties of hundreds of GeV in the DM mass limit.
We should finally note that, for DM multiplets with {I, Y } = {0, 0}, {1/2, 1/2}, {1, 0} and {1, 1}, a linear Yukawa coupling with the SM leptons is allowed by gauge symmetries, while larger isospin multiplets are automatically protected at renormalisable level. However, higher order couplings involving the Higgs can always generate decays of the DM multiplets, and it has been the main motivation of Ref. [45] to find multiplets that are long-lived enough to be Cosmologically stable, thus pointing towards multiplets with I = 2. In this work we will be more pragmatic and allow for any multiplet by forbidding implicitly all operators that could mediate the decays of the DM candidate. The origin of such a symmetry is to be searched in the more complete model containing the DM multiplet. Moreover, as the MCDM models are to be considered effective low energy descriptions of the DM phenomenology, we do not consider the upper limit on the isospin value coming from the absence of Landau poles in the renormalisation group running of the EW gauge couplings below the Planck mass.
After reviewing the properties of Dirac and Majorana multiplets in Sec. 3.1 and 3.2 respectively, in Sec. 3.3 we study in detail the effect of dim-5 couplings to the Higgs field. In Sec. 3.4 we provide novel detailed results on one-loop cross sections for direct detection, including for the first time the mass split in the multiplet, and present current exclusion limits and future projections. We also show that, due to delicate cancellations among various amplitudes, both the mass split and nuclear uncertainties have sizeable impact on the cross sections and on the DM mass limits.

Dirac multiplets (F I Y )
In the case of Dirac multiplets, i.e. when both chiralities are present, the lowest order Lagrangian, to be added to the SM one, reads where the covariant derivative includes the EW gauge bosons. It is invariant under a global U(1) DM symmetry, thus an asymmetric contribution to the relic abundance may be present if the complete model preserves this symmetry. Except for the singlet caseF 0 0 , the multiplet contains charged states: The Dirac mass term in Eq. (3.1) gives equal mass to all components of the multiplet. This degeneracy can only be lifted by radiative corrections due to the EW gauge bosons. This contribution has been first computed in Ref. [45], and can be written as where f F (x) is a loop function and x V = m V /m D . This expression explicitly shows that the mass differences vanish in the limit of equal masses for W , Z and photon. For the loop function, we found an alternative form that is numerically more stable than the one given in Ref. [45] (see Appendix A.1 for more details). The result, which is exact, reads This function has been defined in such a way that f F (x γ ) ≡ f F (0) = 0. It is instructive to study how the mass split looks in the limit of DM mass small and large compared to the W and Z masses. For light DM, M DM ≈ m D m W , the leading contribution reads This mass split tends to zero for vanishing DM mass and is proportional to the difference in squared charges, as an indication that it is dominated by the photon exchange. Furthermore, in this limit the lightest component of the multiplet is always the neural one. In the opposite limit, m D m W , the leading term in the expansion reads For Y = 0, the charged states are always heavier than the neutral one as the surviving term is proportional to the difference of squared charges. On the contrary, for Y = 0 the second term, which depends on the sign of the charges (we chose Y > 0 without loss of generality), does not guarantee that the Q = 0 state is always the lightest one. In particular, the state Q = −1 is always lighter than the Q = 0 one in this limit, for any value of Y = 0 and of the isospin of the multiplet. Thus, there exists an upper limit on m D , above which the lightest state in the multiplet is charged, and this value is determined by the Q = −1 state. The values of the mass upper bounds for various Y are shown in the left panel of Fig. 1: the highest value is achieved for Y = 1/2 which gives m max D ≈ 570 GeV (we recall that for Y = 0 there is no limit), while for Y = 1 we find m max D ≈ 42 GeV, which is already below m Z /2. Hence, multiplets with Y ≥ 1 are excluded by the Z-width measurement in the region where the lightest state is neutral, as long as a Q = −1 state exists in the multiplet. In fact, this upper limit is removed for multiplets with maximal hypercharge, Y = I, for which only states with positive charge are present. In the right panel of Fig. 1 we show the mass splits for various charges and for Y = 1/2 as a function of the DM mass, i.e. the mass of the neutral component. This shows that the Q = −1 state is always the lightest above the neutral one for m D 570 GeV, with a mass split always smaller than 100 MeV.
As already mentioned, all models with Y = 0, i.e. b) and c), are excluded by direct detection via the Z exchange. As we will see, however, a dim-5 couplings to the Higgs can salvage the models with Y = 1/2 (see Section 3.3.2).

Pseudo-Dirac multiplets
For completeness, we recall that Dirac multiplets with Y = 0 can be split in two Majorana multiplets M I 0 . This can be effectively described by the addition of a new mass term to the Lagrangian in Eq. (3.1): Without loss of generality, we consider δm to be real and positive. 3 The Lagrangian above effectively describes two Majorana multiplets (see Section 3.2) with masses We highlighted the mass term δm as it breaks the U(1) DM to a Z 2 , hence it may be a small perturbation depending on how this breaking is implemented in the UV completion of the model. Note also that this term is not generated radiatively as long as it is not generated by the complete model. Hence, it may be natural to have a small mass split between the two Majorana multiplets, which leads to a model with two DM candidates, with the relic density dominated by the lighter one for large mass split. We recall that in all pseudo-Dirac models the lightest component is guaranteed to be neutral.

Majorana multiplets (M I 0 )
In the case of a Majorana multiplet,M I 0 , the Lagrangian to be added to the SM one reads: where Ψ C = Ψ and the multiplet can be written in terms of a Weyl spinor Ψ = χ χ with components so that the Majorana DM candidate χ 0 is accompanied by n = I Dirac charged partners. The phenomenology of this multiplet is in large part the same as for aF I 0 Dirac multiplet, in particular the mass split between the various components is given by the same formula given by Eqs. (3.3) and (3.4). Hence, the lightest component is always the neutral one.

Mass split from dim-5 Higgs couplings
In this section we consider minimal couplings to the Higgs field, which can arise at the level of dim-5 operators. While being suppressed by a UV scale, they are relevant because they can induce a mass split between the components of the DM multiplet, potentially competitive with the EW loops, and change drastically the phenomenology of the multiplet. Hence, while they are not renormalisable couplings, we will consider them here as minimal extensions of the single multiplet models. Furthermore, as we shall see in Section 4, they arise by integrating out a heavier fermion or scalar mediator.

Basic case for Dirac and Majorana multiplets
The Brout-Englert-Higgs doublet φ H , which has I = 1/2 and Y = 1/2, can only couple to the DM multiplet via higher dimensional operators. The lowest order operators have mass-dimension 5 (dim-5) and read: where T a I are the three SU(2) L generators for the multiplet with isospin I, and Λ is a new scale that we assume being beyond the LHC reach to resolve. For Majorana multiplets, however, the first term is absent as it vanishes identically. The second term generates a common mass contribution for all components, thus it simply shifts the mass of the multiplet and generates a coupling to the Higgs, − κ v Λ hΨΨ, that contributes to direct detection, where v = 246 GeV.
The first one, instead, induces a mass split among the various components, thus it may affect the conclusions about the spectrum of Dirac multiplets we reached in the previous section. We recall that the form of the SU(2) generators for a generic isospin I is with c k = k(2I + 1 − k) , k = 1, . . . 2I, and c 2I+1−k = c k . (3.14) Once the Higgs field develops its VEV, the only non-vanishing component is where µ D = − κv 2 4Λ and we have used the relation T 3 = −Y for the neutral component. In terms of mass splitting, these couplings can be expressed as Together with the EW loops in Eq. (3.3), the master formula for the mass splits reads: (3.20) The condition for maximal hypercharge stems from the fact that the Q = −1 state is absent. The allowed ranges of µ D as a function of m D are shown in Fig. 2 for various values of integer and semi-integer Y , where the upper limit should be removed for multiplets with maximal hypercharge. Hence, µ D allows to salvage multiplets with Y > 1/2. It remains the issue of exclusion by direct detection via the Z coupling: to elude it, one needs to generate a mass split in the neutral state that we discuss in the next subsection.
To connect the feasible values of µ D with the scale at which this interaction is generated, it is useful to compare it with the asymptotic value of the EW loops: This corresponds to the range of µ D allowed asymptotically in the Y = 0 case, and gives a reference for the scale of new physics Λ. Models with Y = 0 are excluded by direct detection via the Z coupling. It is well known that this bound can be avoided if the neutral state is split into two Majorana mass eigenstates via a coupling to the Higgs field. For Y = 1/2, this occurs at dim-5 level via the operator: The operator above is similar in nature to the Weinberg operator in the SM [51] that gives a Majorana mass to the left-handed neutrinos. Note also that it preserves a Z 2 symmetry on the DM candidate, but breaks the U(1) DM . Its most important effect is to split the neutral Dirac state into two Majorana mass states: the Z boson can only couple the two states to each other, without any diagonal couplings. As long as the heavier Majorana state is not Cosmologically stable, the DM candidate is the lightest one and elastic scattering off nuclei mediated by the Z is absent. The price to pay is a new coupling to the Higgs boson, which also contributes to direct detection. As a fist step, we need to determine what is the effect of the new coupling κ M on the mass ordering inside the multiplet. In the operator (3.22), the only non-vanishing component of the Higgs current is which couples toΨT − I Ψ C . The resulting Lagrangian for a generic semi-integer isospin I reads where µ M = κ M v 2 4Λ , and we recall that the neutral state corresponds to k = n = I + 1/2. All states receive a mass correction except the one with the largest electric charge, ψ n+ .  For the neutral state, the mass matrix can be written in a Majorana form as follows: wherem 0 D includes the one-loop EW corrections. The Majorana mass eigenvalues are Note that c I+1/2 is the largest coefficient in the T + I generator and the lightest state always receives a negative contribution to its mass, independently on the sign of κ M . Henceforth, this operator always tends to make one neutral state lighter. For a doublet, I = 1/2, the charged state does not receive a mass correction from κ M , hence the mass split between the charged state and the lightest neutral one can be written as This shows that the presence of a non-zero µ M always enlarges the parameter space where the lightest state is neutral, in particular allowing for larger negative values of µ D compared to the case with µ M = 0.
For larger values of the isospin, I > 1/2, we need to study the correction to the masses of the charged states, whose mass matrix can be written as include the one-loop EW corrections. Using the parametrisation adopted in the previous subsection, the mass eigenstates for charge Q = (n − k) states can be written as The state that receives the potentially largest negative contribution to the mass has charge Q = 1, for which c k → c I−1/2 = c I+3/2 = I + 3 2 I − 1 2 and the mass difference between the lighter charged and neutral states reads The lightest state remains the neutral one as long as −δm This region in the m D -µ D parameter space is represented in Fig. 3 for I = 3/2 (Left) and I = 5/2 (Right), where the curves from solid to dashed correspond to increasing µ M from 0 to 300 MeV. This plot shows that a non-zero µ M always enlarges the allowed band. The same trend occurs for larger isospin values. To have a feeling of the scale involved in the generation of µ M , as a reference the minimal value of µ M above which the neutral state is always the lightest for µ D = 0 and I = 3/2 is: A similar splitting can be obtained also for multiplets with hypercharge larger than 1/2, at the price of higher dimensionality of the operator. For any given semi-integer Y = N + 1/2, the operator contains 2N additional φ H fields, hence having a mass dimension of dim = 5 + 2N = 4 + 2Y . The main issue with this case is that a sizeable µ M would require a relatively low new physics scale: . (3.33) For Y = 3/2, this implies Λ < κ   1/2 model in the (µ M , M DM ) plane: the blue shaded region is excluded by relic density over-abundance; the dark pink region is excluded by PandaX-4T [6] DM direct detection searches; the light pink region presents the region that will be probed by future DM direct detection searches with LZ detector [7]. The narrow green band indicates the allowed region with Ωh 2 = 0.12 that is not accessible to future direct detection experiments.
As an example of how the dim-5 operator in Eq. (3.22) is constrained by relic density and DM direct detection experiments, we show in Fig. 4 the exclusion regions for a doubletF 1/2 1/2 model in the (µ M , M DM ) plane. We set µ D = 0, and recall that the mass splits are given by ∆M 0 = 2|µ M | between the two Majorana mass states, and ∆M + = |µ M |+ the EW loops. The blue shaded region is excluded by relic density over-abundance while the dark pink region is excluded by current direct detection limits from PandaX-4T [6]. The light pink region presents the projected region that future DM direct detection searches with the LUX-ZEPLIN (LZ) detector will be able to probe [7]. We can see that DM masses above 1.1 TeV are excluded by complementary relic density and DM direct detection constraints. For the relic density, increasing |µ M | reduces the co-annihilation via the W and Z gauge bosons, hence requiring a slightly lighter mass, while for µ M 8 GeV the Higgs couplings start dominating, pushing the DM mass to higher values. However, this region is already excluded by direct detection, as PandaX-4T excludes µ M above ∼ 2 GeV. The projected LZ limit will probe µ M down to ∼ 250 MeV, a region where the mass split from EW loops becomes relevant. Direct detection due to the EW loops, however, remains too small to be detected, as we will discuss in the text subsection. The narrow green band indicates the region with Ωh 2 = 0.12 that will not be accessible to direct detection experiments.

Loop-induced Direct Detection
Loop-induced direct detection cross sections in DM models with a single multiplet have been explored in several papers [45,52,53]. In particular, Ref. [53] presents complete results at one-loop (including two-loops for the couplings to gluons via a heavy flavour quark), in the limit where the DM candidate is a Majorana state from a pseudo-Dirac multiplet. Furthermore, the masses of the DM multiplet components are considered to be exactly the same. A cancellation is observed among various amplitudes, leading to a cross section that is significantly smaller than what could be naively expected.
Motivated by this cancellation, in this section we revisit the one-loop calculation and extend the results to cases where the DM candidate is a Dirac state and for Majorana multiplets. We also included the effect of mass splits in the DM multiplet: while the mass splits are numerically small, these effects can alter the delicate cancellation among the various terms, hence changing dramatically the final result. Furthermore, we will discuss the impact of uncertainties in the nucleon form factors and parton density functions, which can be highly enhanced by the cancellations. The one-loop diagrams relevant for direct detection are shown in Figure 5, where q (Q) are external (internal) SM quarks. We do not calculate the two-loop diagram resulting from closing the external quark lines for heavy flavours, instead we employ the results of Hisano et al. [53]. Furthermore, our calculations are done for spin-independent (SI) cross sections in the limit of zero external momenta and assuming that internal quark masses are comparable to the external ones. 4 The amplitudes can be parametrised in terms of the following effective Lagrangian [53]: where D is the DM fermion, which may be pseudo-Dirac, Majorana or Dirac, and the Twist-2 quark current is given by The first term in Eq. (3.34) proportional to f q is the scalar-scalar (SS) operator, the second and third proportional to g (1),(2) q are twist-2 operators, and the last one proportional to f G describes the effectively two-loop coupling to gluons. The coefficients can be explicitly computed (see Appendix A.2) and give: The couplings of the W boson are explicitly given for a multiplet with n = 2I + 1 and hypercharge Y . The vector and axial couplings for the quarks are given by a 0 Fig.5 only contributes to the SS operator and gives rise to the loop function ∆ H , while diagrams A and B generate the loop functions ∆ S , ∆ T 1 , ∆ T 2 . Finally the normalisation factors κ z and κ w depend on the nature of the DM candidate: if D is pseudo-Dirac (κ z , κ w ) = (1, 1), whereas for Majorana (κ z , κ w ) = (0, 1/2) and for Dirac (κ z , κ w ) = (4, 2).
The mass splits within the DM multiplets are represented by the parameters y 0 , y + and y − , which encode the mass split between the two Majorana mass states (y 0 = 0 for a Dirac multiplet, while the whole term vanishes for a Majorana multiplet as Y = 0) and between the charge Q = ±1 states and the DM state, respectively. In the limit of zero mass splits, i.e. y i → 0, our results reproduce the formulas in [53], which we report below for reference: As already mentioned, for the coupling to gluons we use the two-loop computation presented in Ref. [53] for vanishing mass splits of the DM multiplet. The contribution can be expressed in terms of long-distance (LD, dominated by momenta of the order of the light quark masses) and short-distance (SD, dominated by momenta of the order of the W/Z bosons or of the DM states) contributions, as follows: where the LD contribution of the light quarks are taken into account in the SS coefficients f q and NLO corrections in QCD are embedded in the coefficients c Q . Explicit results can be found in Ref. [53].
The SI cross section for DM scattering off target nucleon N is expressed as Here, f T q are the proton form factors for the quarks, while f T G = 1− q=u,d,s f T q applies for gluons (note that charm, bottom and top are considered heavy flavours in this formula, and associated to the gluon form factor), while q(2) andq(2) are second moments (evaluated at µ = m Z ) for quarks and anti-quarks respectively. In the zero mass split case, a strong cancellation has been observed between the contribution of the twist-2 operators (g (1) q and g (2) q ) and the gluon (f G ) contributions, while the SS one (f q ) tends to be smaller [53]. This result is shown in Fig. 6, were we plot the SI cross sections for various cases compared to the current exclusion from PandaX-4T [6] and the projection from the future LZ [7]. The border of the yellow shaded region labelled "Neutrino Floor" corresponds to the sensitivityestimate achievable at each DM mass for a one neutrino event exposure at liquid Xenon detectors [54]. We have digitised data for the PandaX-4T, LZ and neutrino floor limits, and they are now publicly available on the PhenoData platform [55][56][57]. We would like to note that the neutrino floor limit for one neutrino event can be improved (i.e. lowered) by future experiments with lower energy threshold [54] potentially by about one order of magnitude. One can see that only multiplets with I ≥ 2 can be completely probed by LZ up to masses of 10 TeV. The observed cancellation, however, is very sensitive to two important effects: the nuclear uncertainties on the form factors and on the second moments, and the mass splits within the DM multiplet. We will discuss both below, starting from the former.

Impact of uncertainties on nucleon form factors and parton distribution functions (PDFs)
The proton form factors for light quarks may be calculated [58] in terms of light quark mass ratios, m u /m d = 0.46 ± 0.05 and m s /(m u + m d ) = 13.75 ± 0.15, and quantities associated with nucleonic matrix elements, Σ πN = 46 ± 11 MeV, σ s = 35 ± 16 MeV and z = 1.258 ± 0.081. Explicitly, they are given by In order to combine errors from all sources, we use the Monte Carlo approach whereby we estimate the sampling distribution of the cross section via the generation of points from the sampling distributions of underlying parameters. For the form factors, we sample from a multivariate Gaussian defined by the input parameters given above, assuming that errors are uncorrelated. The distribution of form factor values are computed using Eqs. (3.41).
The uncertainties on the second moments q(2) andq(2) derive from the uncertainties in the parton distribution functions (PDFs). We take into account both the uncertainties in the PDF fitting procedure, and in the scale variation. In practice, we concurrently sample from the CTEQ18NLO [59] PDFs, using the Hessian implementation in LHAPDF [60], before numerically integrating these PDFs to generate the second moment values. We probe the variation from PDF scale by sampling from a log-normal distribution for the PDF scale, µ, with central value µ = m Z such that the 1σ bands fall on µ = m Z /2 and µ = 2m Z (i.e log 2 (µ/m Z ) is normally distributed with mean 0 and standard deviation 1).
The uncertainties propagated on the SI cross section are depicted in Fig. 7, where, for each model, the solid line represents the mean and the band signifies the 95% confidence interval (we show results for pseudo-Dirac case, as the Majorana and Dirac cases can be obtained by a simple numerical scaling). For comparison, the dashed lines show the results of Ref. [53], where the values used are f T u = 0.023, f T d = 0.032, f T s = 0.020 and second moments (evaluated at µ = m Z ) u(2)(ū(2)) = 0.22(0.034), d(2)(d(2)) = 0.11(0.036), s(2)(s(2)) = 0.026(0.026), c(2)(c(2)) = 0.019(0.019), b(2)(b(2)) = 0.012(0.012). The difference is due to the fact that we use a different set of PDFs. To study the effect of scale dependence from the PDFs, we also show in various dotted styles the values of the cross sections for µ = m Z /2, m Z , 2m Z respectively, while keeping all the other parameters fixed. The result shows that the main contribution to the uncertainties derive from the form factors.
For models with I ≥ 1, the uncertainties in the cross section derives in a sizeable uncertainty in the bound on the DM mass, which amounts to several hundred GeV. For the doublet, I = 1/2, instead, we observe that a cancellation may occur for M DM > 450 GeV, hence making a prediction for the direct detection reach impossible in this mass range. Nevertheless, the cross section always remains below the reach of LZ, and will likely escape detection.
. Impact of form factor and PDF uncertainties on the SI cross section for various models with a single DM fermion multiplet. We show results for pseudo-Dirac cases. For Y = 0, Majorana and Dirac cases can be obtained by a simple scaling of the cross sections by a factor of 1/4 and 4 respectively. We recall that only the pseudo-Dirac case is allowed for Y = 0.

Impact of mass splits
As we have shown, the current uncertainties in the nucleon form factors and PDFs produce relevant uncertainties on the SI cross sections, resulting in variations of several hundred GeV in the DM mass bound or producing cancellations in the doublet case. Even if these uncertainties were substantially reduced, the one-loop calculation is sensitive to the mass splits within the DM multiplet, which were not taken into account so far. Here, we extended the loop computation to take into account mass splits, and we re-evaluate the cancellations observed in Ref. [53]. As input parameters, we use the same ones used in [53] and recapped at the end of last section.  Variation in GG Figure 8. Impact of mass splits on the one-loop cross sections for pseudo-Dirac multiplet. On the left, we show the contribution of various operators with and without mass splits, and the impact on the total. The variation is due to a 5% variation in the GG contribution. On the right we show the total contribution with variation for other multiplets with larger isospin, I ≥ 1.
Our calculation takes into account the mass splits in the one-loop results for both SS and twist-2 operators, however this effect is not included yet in the two-loop computation for the gluon couplings. To supply to this, we include a 5% variation to the latter. In the left panel of Fig. 8, we show the impact of the radiative mass split (δm EW ) for the doublet case. We show separately the contribution of SS and twist-2 operators, comparing the result with (dashed) and without (solid) mass splits taken into account. While the effect on each amplitude is small, the total cross section sees a substantial reduction once the mass splits are taken into account. As the two-loop GG contribution does not include the mass splits, we consider an additional 5% variation: as shown by the band, this small effect can cause cancellations to occur for M DM 400 GeV. Numerically, this effect is similar to the impact of the uncertainties, hence showing that the two are comparable.
In the right panel of Fig. 8, we show the total contributions for models with I ≥ 1. An enhanced cancellation is also observed in this case when comparing the solid lines (with the mass splits) to the dashed ones (without). The 5% variation in the GG contribution, instead, generates an uncertainty that derives in an uncertainty of several hundred GeV on the DM mass. The results presented here show that both mass splits and nuclear uncertainties produce similar effects on the SI cross sections, hence motivating further studies in this direction.

Fermionic Dark Matter with one additional multiplet
In this section we present the classification of models that contain one additional multiplet (mediator multiplet), in addition to the DM one. The mediator multiplet can be either odd or even under the symmetry protecting the stability of the DM candidate, and its quantum numbers are limited (and defined) by the requirement of the renormalisability and gauge invariance of its interaction with the DM multiplet (and SM fields). We recall that we use different labels F /F and M /M for Dirac and Majorana fermion multiplets, respectively, since the two choices often lead to rather different models when a mediator is present, as we discuss below.

Even scalar mediator (F
. The case of a scalar mediator that couples to the SM has been one of the first scenarios considered in simplified models (see e.g. [26,[61][62][63][64][65][66]), however it has been by now established that it is not simple nor minimal to achieve phenomenologically relevant models once the simplified case is included in a fully gauge-invariant model [67,68]. In particular, couplings of a single scalar to SM fermions are hard to obtain without breaking the EW symmetry, while couplings to gauge bosons only arise at dim-5 operator level unless the scalar is allowed to develop a non-zero vacuum expectation value. In the following, we will limit ourselves to the most minimal scenarios and not consider higher dimensional operators.
The models we consider feature a Yukawa coupling connecting two DM multiplets and a (peudo)scalar mediator multiplet. They can be classified as follows (for the sake of minimality we consider multiplets with zero hypercharge, S I 0 , to be real): The properties of all the possible models are summarised in Table 2, where we identified five template scenarios with distinct properties. Note that only integer isospin and hypercharges of the scalar mediator are allowed. In the 4th and 5th columns ("DM sym." and "Etx. sym.", respectively), we list the largest Dark symmetry allowed by the above Yukawa couplings, which could be broken by the couplings of the scalar mediator multiplet Φ to the SM. The last column contains the scalar mediators that can have linear (renormalisable) couplings to the SM.
The most general Lagrangians, for the real and complex scalar multiplets (with integer hypercharges), read: where V is a generic potential for the scalar. Note that, for the real case, only one Higgs portal coupling is allowed due to the fact that Φ has integer isospin. The term V linear contains eventual linear couplings Model D1/D2 D3/D4 DM sym. Ext. sym. Ext. Charges Linear to SM of Φ to a SM operator, which can be made of the Higgs field or leptons. Only 5 such cases occur: where i, j are flavour indices. The linear Higgs portal coupling, allowed only for CP-even S 0 0 and S 1 0 and for the charged iso-triplet S 1 1 , necessarily implies that the scalar mediator acquires a vacuum expectation value Φ = 0 via the Higgs one, thus a universal coupling to SM fermions is generated via the mixing with the physical Higgs boson [69,70]. However, this mixing is strongly suppressed in the triplet cases because of three-level contributions to the ρ parameter, while in the singlet case milder (but still important) bounds derive from the measurement of the 125 GeV Higgs couplings [71]. This shows that the couplings of the scalar mediator to SM fermions and gauge bosons are deemed to be small. If the scalar mediator multiplet is much heavier than the EW scale and the DM mass, the bound on the coupling can be weakened. Integrating out the scalar multiplet generates the dim-5 operators between the DM multiplet and the Higgs boson we introduced in Section 3.3. Depending on the quantum numbers of the scalar multiplet, the following possibilities are realised for Dirac multiplets: The singlet S 0 0 also generates the operator for the Majorana DM multiplet case via the coupling D1: (4.11) The couplings κ, κ and κ M are defined in Eqs (3.11) and (3.22). The other two cases only allow for couplings to leptons. The triplet S 1 1 in Eq. (4.7) corresponds to type-II see-saw models [72][73][74] for neutrino mass generation, and it has been studied in connection to DM in Refs [75][76][77][78]. The doubly-charged scalar S 0 2 in Eq. (4.6) also contributes to neutrino masses, as it breaks lepton number by two units, and it has been studied in Ref. [79] paired with a scalar DM multiplet.
We finally note that a vacuum expectation value for the scalar mediator can be induced in all cases, in particular via the quartic coupling to the Higgs field yielding Φ † Φ = 0, and with all the limitations and bound described above. The phenomenology of such cases follow the analyses done in the simplified models [67,80].
If a vacuum expectation value is not generated, then the presence of the coupling to the scalar mediator multiplet does not affect the mass spectrum of the DM multiplet: this implies that models with Y = 0 are excluded by Z-mediated direct detection, while Y = 0 models are probed at one-loop level, as discussed in Section 3.4, with additional contributions of the scalar mediator if it couples to quarks.
There is, however, a new class of mediators that arise from our classification: scalar mediators that only have bilinear couplings to the SM Higgs field. Such models have new interesting features that we will study in detail in the next section. For now, we content ourselves to classify the relevant models: (a) Accidental stability: the scalar mediator multiplet can be accidentally stable if all linear couplings to the SM are forbidden, and it is lighter than twice the mass of the fermionic DM multiplet. One should note, however, that if the DM multiplet has I = 0, then due to the couplings D1-D4, triangle loops generate couplings of the scalar mediator multiplet to two EW gauge bosons, hence making it unstable. The only model with accidental stability is, therefore,F 0 0 S 0 0 with coupling D2 (i.e. CP-odd mediator). 5 (b) Protected stability by Z 4 : in models with an extended Z 4 symmetry, i.e.F I 0 S I 0 with couplings D3/D4, the stability of the mediator is guaranteed by a discrete charge. Direct detection bounds apply via loop induced couplings for the two DM components, like the ones discussed in Section 3.4. The most minimal surviving model involves gauge singlets,F 0 0 S 0 0 with D4.
(c) Protected stability by U(1): similarly, stability can be guaranteed by a U(1) symmetry inF I Y S I 2Y models. In such cases, however, the fact that Y = 0 requires that a Majorana mass split is generated in the neutral DM fermionic candidate. This can only be achieved in modelsF I 1/2 S 1 1 with the couplings in Eq. (4.5) included: this however explicitly breaks the U(1) symmetry and allows decays of the mediator. By integrating out S 1 1 , or by a small vacuum expectation value, the same mass split induced by the Higgs operator discussed in Section 3.3.2 will arise. The only difference would be the presence of a coupling to a scalar mediator, which can affect the relic density computation.
We note that in all cases, direct detection from the fermionic DM candidate may be avoided if the dominant contribution to the relic density is coming from the stable scalar, however this case will best fit under a scalar DM multiplet study [81]. 5 ForM 0 0 S 0 0 , the scalar mediator is CP-even, thus linear coupling to the Higgs cannot be forbidden.
One should pay a special attention to constraints from EW precision data, in particular from the ρ 0 parameter [82], which is very close to one in the SM and measured with per-mille precision. In the general case of an arbitrary number of SU(2) scalar multiplets, ρ 0 takes the form (see Eq.(10.58) in Ref. [83]): ρ 0 = n [I n (I n + 1) − I 2 3n ]|v n | 2 2 n I 2 3n |v n | 2 , (4.12) where I n , I 3n and v n are the isospin, the third component isospin of the vacuum state and the the vacuum expectation value for the n th scalar multiplet, respectively. Hence, besides the known cases with doublet and singlet, strong bounds from ρ can be avoided in several other non-trivial cases: for example, in the model with a septet scalar (S 3 2 ) [84] that can couple to a DM quintuplet, or the model with custodial combinations like triplets in the Georgi-Machacek model (S 1 0 + S 1 1 ) [85]. The latter would rather be the part of less minimal models, but still possibly quite interesting. From Eq. (4.12) one can see that the caseF I 1/2 S 1 1 described in point (c) can better fit in a Georgi-Machacek scenario, where the triplet VEV is not too constrained. However, for the scenario with a custodial violating triplet only, the coupling may be enough to generate a large enough mass split to avoid constraints from ρ 0 .
To summarise this section, we found a new class of relevant minimal models with a scalar mediator that is (accidentally) stable: this includes a model with two singlets,F 0 0 S 0 0 with D2 or D4, which we study in more details in Section 5.

Odd scalar mediator (F
In this class of models, the DM fermion multiplet Ψ couples to the odd scalarφ and to a SM fermion via a Yukawa coupling: the quantum numbers of the scalar multiplet are, therefore, fixed by the properties of the chosen SM fermion. As the SM fermions are chiral, one can classify two cases, distinguished by their chirality (a SU(2) L doublet f L or a singlet f R ): -for left-handed SM fermions, the respective interactions read: hence,φ f L = {I ± 1/2, Y − Y f } (and an anti-triplet of QCD colour if f L is a quark); -for right-handed SM fermions: (4.14) hence,φ f = {I, Y − Y f } (and an anti-triplet of QCD colour if f R is a quark).
Note that i = 1, 2, 3 is a SM family index, and the two types of couplings cannot co-exist with the same multiplet in minimal models. In other words, the couplings of the mediator necessarily involve one chirality and only one type of SM fermions. The scalar multiplet will also have couplings to the Higgs [46], in a form analogous to that of Eqs (4.1) or (4.2), but in the absence of any linear coupling forbidden by the DM parity. As the cases of quarks and leptons lead to rather different physics, we will discuss them in detail separately.

Quark-type mediators
Firstly, as quark partnersφ q L/R carry QCD charges, they cannot constitute part of the DM relic density and are always required to be heavier than the DM fermion candidate. Besides the effect of EW interactions discussed in the previous section, the scalar mediator will contribute a new tree-level process to direct detection, Dq →φ → Dq, whose rate is determined by the value of the h qL/R coupling and the mass ofφ mediator for any given DM mass. The most minimal, and safe, cases involveF 0 0 andM 0 0 , for which the scalar mediator has the same quantum numbers as the corresponding SM fermion. This case is a template of supersymmetry (φ q L/R being one of the squarks), and has been studied in detail in simplified models withφ q L mediator and Majorana DM [86][87][88].

Lepton-type mediators
In this case, the scalar multiplet may contain a neutral state and therefore also play the role of DM (this case will be covered in a future work [81]). In the case where the DM arises form the fermionic multiplet, direct detection (for the only surviving "safe" cases ofF 0 0 andM 0 0 ) occurs only at one-loop level contrary to the case of the coloured scalar mediators discussed above. DM direct detection rates, in this case, are defined by the respective h lL/R Yukawa coupling and the mass of the scalar multiplet, occurring in the loop. This case also corresponds to the supersymmetry template with sleptons, and has been covered in Refs [90,91].

Even fermion mediator (F
This case does not allow for renormalisable couplings between the mediator and the DM multiplet, however we list it here for completeness and because it leads to interesting new models of leptophilic DM. The only allowed coupling involves one mediator multiplet, Σ, and three DM multiplets Ψ. In turn, the even multiplet Σ needs to couple to the SM via a Yukawa-type coupling to leptons (quarks are excluded to avoid QCD charged DM).
The DM mediator coupling comes from a dim-6 operator: which preserves a Z 3 DM parity [92,93] for a complex Dirac multipletF I Y . 6 Moreover, the hypercharges are related by: (4.16) The last relation imposes a significant constraint on the mediator multiplet, as the hypercharge of the DM one needs to be semi-integer for semi-integer isospin and integer for integer isospin in order to have a neutral component. As a consequence, the only allowed cases (with Yukawa couplings to leptons) are: Class B: ∆L = −ξ RlR φ H Σ + h.c. ;F I=semi-int.
Due to direct detection constraints, and the role played by gauge interactions in the thermal relic abundance (which would make the mediator irrelevant), the only interesting case appears for a singlet DM, F 0 0 F 0 0 , which belongs to class A. Note that Σ is effectively a heavy right-handed neutrino. The relic density will thus be determined by the processes: ΨΨ ↔Ψν , ΨΨ →ΨνH . If the coupling to the SM is very small, being related to neutrino mass generation, then this could be an effective FIMP model. ). Note that in our notation the first multiplet is the one that has the largest component in the DM physical state.
The possible models can be classified based on the form of the Yukawa coupling: where in our convention Ψ indicates the mediator multiplet. Note that the Higgs field may appear as is, or in the form of the complex conjugateφ H . Also, either the mediator or the DM multiplet can be of Majorana nature if either Y = 0 or Y = 0. In general, this class of mediator models have similar features as the simple DM multiplet cases, with an additional coupling to the Higgs boson that could make direct detection more critical.
One point of interest, though, is the fact that in the case of large mediator mass, i.e. M m, by integrating out the mediator multiplet one can generate the dim-5 couplings to the Higgs discussed in Sec. 3.3. In the case of Dirac multiplets, the coefficient of Eq. (3.11) are matched to the Yukawa coupling and mediator mass M as

Even vector mediators (F
Vector mediators are very popular in the simplified model approach to DM phenomenology (see e.g. [94][95][96][97][98][99][100][101][102][103][104][105][106][107][108][109][110][111][112]) mainly because they allow for "gauge invariant" couplings to vector currents of SM fermions. Nevertheless, it is not a simple task to find a consistent, truly gauge invariant, renormalisable model containing vector mediator multiplets. As the vector multiplet couples to a current containing the DM multiplet, the Lagrangian takes the form where P L/R are chirality projectors. As the hypercharge always vanishes (and the isospin is integer), we can always consider real multiplets. For a generic vector field V µ , the most general Lagrangian up to renormalisable couplings reads [113]: where W a µν is the energy-stress tensor of SU(2) L . The second line contains couplings to currents of SM fermions and the Higgs field, compatible with the quantum numbers of the vector multiplet: they are allowed only for the singlet V 0 0 and a triplet V 1 0 . The Lagrangian in Eq. (4.24), which we require to be renormalisable and consistent, needs an additional scalar sector which breaks the gauge symmetry, for which these vector mediators are being gauge bosons (see e.g. [114]). These gauge bosons can come from different theory space, including supersymmetric, extra-dimensional or composite/technicolor origin.
In Ref. [115] it has been shown that the self-interactions of the multiplet can be fixed in order to preserve perturbative unitarity in the scattering amplitude of vector multiplets, however Ref. [116] later showed that violation of perturbative unitarity occurs once the vector multiplet couples to massive gauge bosons (i.e. it is charged under a broken gauge group, like SU(2) L ) and/or to the Higgs: thus new states need to be included in order to restore the consistency of the model. They might affect the low energy properties of the theory by introducing phenomenologically relevant operators. In theories of this kind, the vector mediator may arise as a composite spin-1 meson of a confining strong dynamics, like in models of composite Goldstone Higgs.
One way to avoid these issues is to introduce the vector multiplet as a gauge field: in general, though, a vector carrying isospin needs to come from a model where the gauge symmetry SU(2) L is extended and broken at higher scales. Now, generating the couplings to the SM fermions becomes the challenge, as new fermions are likely to be needed in order to complete multiplets of the extended EW gauge symmetry. Note that here the chiral nature of the SM fermions is the main obstacle, as it may imply the presence of other chiral fermions.
One case that does not suffer from such problem is the singlet, V 0 0 , as it could arise from a broken gauged U(1) symmetry under which some SM fermions are charged. Once again, though, a consistent theory would require an anomaly-free U(1), thus either additional charged heavy states are added, or one has very limited choices, as discussed in various DM Z -portal studies cited above, see e.g. [105] and references therein.

Odd vector mediators (F
In the case of odd vector mediators, the only allowed couplings must involve the DM multiplet and a SM fermion. The classification of mediators, therefore, follows the same as the scalar odd mediators in Sec. 4.2: -for left-handed SM fermions, the coupling reads:

and an anti-triplet of QCD colour if f is a quark);
-for right-handed SM fermions:

and an anti-triplet of QCD colour if f is a quark).
As the mediator typically has non-zero hypercharge, the Lagrangian (4.24) needs to be extended: Similarly to the case of even mediators, the above Lagrangian cannot be complete because of perturbative unitarity violation or the need to extend the gauge symmetries of the SM to generateṼ as a gauge boson. In such a scenarioṼ can play a role of a DM candidate if it is lighter than the fermionic DM candidate. An example of complete model for vector DM involved in the weak interactions has been suggested in [117], where the authors introduce two additional SU(2) triplets -one odd and another even -to make the model consistent.
5 Phenomenology of a new representative model:F 0 0 S 0 0 (CP-odd) In this section we study the modelF 0 0 S 0 0 (CP-odd) with a Dirac fermion singlet (Ψ ≡ ψ) and a pseudoscalar (CP-odd) singlet (Φ ≡ a) -probably the simplest two component DM model introduced in section 4.1. We have reported a preliminary study on this model in [118,119]. During completion of this work, an alternative, partly overlapping, analysis of the same model (without the study of the loop effects) appeared in [120].
The Lagrangian of the dark sector, to be added to the SM one, reads: where φ H is the SM Higgs doublet field. A similar model has been investigated in Ref. [121] where, however, a linear coupling of the pseudo-scalar with the Higgs was also allowed (hence breaking CP), which leads to a very different phenomenology, as this coupling implies that the pseudo-scalar develops a vacuum expectation value. The model contains three new couplings: the Yukawa coupling Y ψ connecting the scalar mediator a to the fermion DM ψ, the self-interaction λ a of the pseudo-scalar a and the quartic coupling to the Higgs λ aH . The latter is the only coupling connecting the new sector to the SM via a Higgs portal. We recall that a linear coupling of a to the Higgs field is forbidden by different CP properties of the Higgs and a. The invariance under CP is preserved as long as a does not develop a vacuum expectation value, i.e. if where m a is the physical mass of the scalar particle, which together with m ψ and three couplings comprises the set of five parameters defining the model: m a , m ψ , Y ψ , λ aH and λ a . (5. 3) The first four parameters only are relevant to the phenomenology we discuss here. We will be working in the region of the parameter space defined by Eq. (5.2), where the phenomenology is very different from the model in Ref. [121] as we have mentioned earlier. As ψ couples exclusively and bi-linearly to a, it is a stable fermionic DM candidate protected by a dark U(1) global symmetry. The pseudo-scalar mediator a can only decay into a pair of DM fermions. Hence, if m a < 2m ψ , a is said to be "accidentally" stable and can contribute to the relic density as a second DM component. In this case the stability of a is protected from its decays to SM particles at all loops because the CP symmetry is conserved in the dark sector. Indeed, a only couples bilinearly to the SM via the Higgs portal and only CP violation can allow for a linear coupling of a to a SM operator. In this sense, it is the CP symmetry itself that prevents a from decaying into SM states.
The interesting dynamics of this model, where a is in touch with the SM via the Higgs portal coupling λ aH while ψ only interacts with a, leads to four distinct regimes of relevance for DM phenomenology, as summarised in Table 3: • In scenario A, both fermion and pseudo-scalar can thermalise with the SM states. If m a ≤ 2m ψ , then a is stable and contributes to the relic abundance. Conversely, if m a > 2m ψ , then it is unstable and merely acts as a mediator for the interactions of the fermionic DM to the SM.
• In scenario B, the relic abundance of ψ is determined by the freeze-in mechanism, driven by the very small value of Y ψ , while a contributes as a thermal DM component for m a < 2m ψ . In the parameter space, where m a > 2m ψ , the smallness of Y ψ can lead to a being long-lived, decaying into a pair of ψ.
• In scenario C, both new particles can freeze-in via their small couplings to the SM sector (including the loop-induced coupling of ψ, as we discuss below), before thermalisation between the two species.
• In scenario D, both particles have very small couplings. While a can freeze-in via its coupling to the Higgs portal, the coupling of the fermion is too small and would lead to a negligible contribution to the total amount of relic density. Depending on its mass, a can be the only significant DM candidate, or decay promptly to the fermion ψ after being produced in the early universe.  Note that any other range of the couplings is excluded by DM over-production or out of control because of perturbativity loss. Furthermore, in scenarios C and D, direct and indirect detection experiments, as well as colliders, would be unable to observe either of these new particles due to the feeble couplings to the SM. In contrast, in scenarios A and B, a may be observable due to the sizeable Higgs portal coupling. In scenario A, the fermion may also be directly observables due to a loop-induced coupling to the Higgs, as we will discuss below. Implementation of this model along with the LANHEP [122] source and libraries required for one-loop calculations have been made publicly available at HEPMDB [123].
Let us start the discussion of the model's phenomenology by presenting some generic features of the new states, ψ and a. If a is stable, its DM fraction can be revealed via direct detection thanks to the following SI elastic cross section on nuclei: where the nucleon effective coupling λ N (N labels the nucleon type) can be written in terms of the nucleon form factors presented in section 3.4 as: The fermion DM, ψ, which is always stable, couples to the SM only via the mediator a. The coupling of ψ to the Higgs boson is, however, generated at one loop level. The complete expression for this coupling is given in Appendix A.3, where δY DD refers to δY (Eq. A.26) evaluated at the direct detection scale, t = 0. In the limit of small m a , the effective Hψψ Yukawa coupling, is given by For larger m a , the loop-induced coupling decreases monotonically, with δY ∝ m −2 a asymptotic for large a masses. This coupling is only relevant when both Y ψ and λ aH are sizeable, and it contributes to direct detection via the following SI cross section of ψ on nucleons: As an illustration, we show in Fig. 9 the SI cross section as a function of the masses, rescaled by the tree-level couplings. Taking into account that the current direct detection limit is in the 10 −10 − 10 −9 pb range, we can infer that this process provides relevant limits only for relatively small ψ masses and couplings of order unity. [pb] Figure 9. Loop-induced direct detection cross section for ψ scattering on nucleons σ SI ψ , scaled by the tree-level couplings (λ aH Y 2 ψ ) 2 , as a function of the masses in GeV.
Another important constraint arises in the region, where the pseudo-scalar and/or the fermion are lighter than half the Higgs mass, i.e. m a , m ψ < m H /2, thanks to the LHC limits on Higgs invisible decays. For the pseudo-scalar, the partial decay width is generated at tree-level: For the fermion, the decay is induced via the one-loop induced coupling in Eq. (5.6). Hence, the loopinduced H → ψψ partial decay width is given by where the effective coupling δY H→ψψ depends on a loop function Υ H→ψψ ≡ Υ(s = m 2 H ) (see Appendix A.3) We recall that a always leads to missing energy, even when it decays promptly. One should also note that the loop-induced coupling δY Relic ≡ δY (s ≈ 4m 2 ψ (1 + 1/(2x))) (where x is mass to temperature ratio) coupling also plays a role for the relic density computation (see Appendix A.3), and is fully taken into account in our numerical results.

Scenario A: 2-component thermal Dark Matter regime
In this scenario (see Table 3), the λ aH and Y ψ couplings are large enough to thermalise both DM components in the early universe.
The relic density in this regime can be evaluated using two coupled Boltzmann equations (see Eq.(5) of Ref. [124]), which are defined by the annihilation and co-annihilation processes, Feynman diagrams of which are shown in Fig. 10. The equations for the two relic densities n a and n ψ read: wheren a andn ψ denote the equilibrium number densities for the two components, and σ v ≡ σv .  We have performed a random scan of the 4-dimensional parameter space of the model and have used MicrOMEGAs [125,126] to evaluate the DM relic density and direct detection rates in the following range of the parameter space: 10 GeV < m ψ < 10 TeV , 10 −1 < Y ψ < 10 , 10 GeV < m a < 1 TeV , 10 −4 < λ aH < 10 . (5.13) The upper limit on the couplings is defined by the loss of perturbativity criteria. We determine the allowed regions surviving after imposing the following constraints: • We use the relic density fit from PLANCK [3] Ω PLANCK h 2 = 0.1186 ± 0.0020 and require Ω 2 h < 0.12 , (5.14) which allows the under-abundant model points.
• We impose the DM direct detection constraints from PandaX-4T [6], which are dominant over the DM indirect detection constraints, as we have explicitly checked.
• We use the invisible Higgs decay constraints at the LHC from ATLAS [127], requiring The results of the scan are presented in Fig. 11, where we show 2D projections of the allowed parameter space for theF 0 0 S 0 0 (CP-odd) model after imposing the constraints listed in the top of each frame. The colour map indicates the relic density normalised to the PLANCK value (Ω PLANCK h 2 = 0.12) for the two DM components a (Ω a /Ω PLANCK , shown in green fading to yellow) and ψ (Ω ψ /Ω PLANCK , shown in magenta fading to cyan), or their sum (Ω tot /Ω PLANCK , shown in black fading to red).
In the top row of Fig. 11 we show the projection in the (m a , λ aH ) plane, where the colour map corresponds to values of Ω a /Ω PLANCK with dark green marking model points that saturate the relic density with a alone. Recall that we keep all points with Ω tot h 2 < 0.12. In Fig. 11(a), no other constraint except the relic density is added: it clearly demonstrates the correlation between Ω a and the value of λ aH , driven by the Feynman diagrams c)/e) and a)/f) of Fig. 10. One can also see the region of the resonant annihilation through the Higgs boson, aa → H, which takes place for m a m H /2. Due to its efficiency, it allows the value of λ aH to go as low as 4 × 10 −4 while being consistent with the Ω PLANCK constraint. Outside of the resonant region, values of λ aH below 10 −2 are excluded by overclosure of the universe. Furthermore, in Fig. 11(b) we present the same 2D projection with points satisfying, in addition, the DM direct detection constraints from PandaX-4T experiment (both on a and on ψ). The plot illustrates how PandaX-4T excludes all points for m a m H , except for a sliver close to the Higgs resonance, which has small couplings or small relic density for the a component, and a few points with very low a relic density (in yellow). One can see that all points with m a m H /2 below the aa → H resonant annihilation region are excluded by PandaX-4T experiment. This happens since in this region the Ωh 2 ≤ 0.12 constraint requires the value of the λ aH coupling to be above 0.1 that, in turn, leads to the SI DM direct detection rates to be above the PandaX-4T limits.
One should also note that, due to the specific set of DM annihilation and co-annihilation diagrams shown in Fig. 10 and their interplay with each other, the relic density constraint requires m a < m ψ in the whole parameter space, except the loop-induced ψψ → H annihilation region (we comment on this region below in more details). This region appears as a vertical strip in Fig. 11(c) for m ψ m H /2. Remarkably, this implies that a is a stable DM component in the whole allowed parameter space, except for the Higgs funnel region for ψ, where a can decay in the fermion DM component.
In Fig. 11(b) we also superimpose the LHC bound on the Higgs invisible decays into a, Br(H → invis) < 0.11, which excludes the Higgs resonant sliver for λ aH 3 × 10 −2 , as shown by the shaded region above the blue line. One can see that this bound is very complementary to the PandaX-4T constraint. Future collider projections are considered as well, showing that the exclusion on λ aH will improve by a factor of about 3 at the High Luminosity LHC run (HL-LHC) (projected bound of Br(H → invis) < 3.8% [128]), as shown by the orange line. The International Linear Collider (ILC) running at √ s = 250 GeV and with an integrated luminosity of 1.15 ab −1 will be able to exclude λ aH 4 × 10 −3 , as indicated by the green line, corresponding to a projected bound Br(H → invis) < 0.4% [129]. One should also note (c) (d) Figure 11. 2D projections of the allowed parameter space forF 0 0 S 0 0 (CP-odd) model (after constraints given at the top of each frame) with the colour map indicating the individual relative relic density of two DM components a (Ω a /Ω PLANCK ), ψ (Ω ψ /Ω PLANCK ) or their sum (Ω tot /Ω PLANCK ). The points with relic density below 10 −3 Ω PLANCK are shown with colour corresponding to the smallest value. that even the ILC will not be able to fully exclude the Higgs resonant region, since λ aH goes below the ILC sensitivity by one order of magnitude.
Besides the Higgs sliver, a second viable region in the parameter space emerges for m a m H , as shown in plot 11(b). It is defined by the interplay of the co-annihilation processes ψψ → aH and aψ → Hψ, involving both new states of the dark sector. This is clearly illustrated by Figs 11(c-d), in the plane defined by the masses and the ψ mass and coupling, respectively. Fig. 11(c), showing a colour map corresponding to the total relic density Ω tot , offers the best view of this region. Besides the Higgs sliver for a, appearing as a horizontal band, the allowed points highlight a vertical strip corresponding to the Higgs resonant region for ψ via the one-loop induced coupling, with An interesting feature is the fact that masses below m H /2 are excluded for both DM candidates: while for a this is due to direct detection and (more marginally) by the Higgs invisible width, for ψ this comes from the fact that for low masses the only efficient annihilation channel is ψψ → aa. This is efficient enough only for m ψ m a , thus, m ψ < m H /2 would result in too much relic density due to the limit on m a . In Fig. 11(d) we show the allowed points projected on the m ψ -Y ψ space, with colour map corresponding to the individual relic density of ψ. We can see a clearly defined triangular shape, which emerges from the ψψ → aa annihilation process and which requires the coupling Y ψ O(1) to be fairly large to avoid overclosure of the universe. On top of this, there is a "leakage" of points for m ψ m H , which emerge from the interplay with the process ψψ → aH, which becomes relevant for m ψ m a ∼ m H . This means that for each value of Y ψ , one can find a value for λ aH that fixes the relic density below the limit. We also observe points with small Y ψ for masses below m H : this is due to an interplay between the two processes aa → H and ψψ → aH above the threshold m ψ   . Direct detection cross-sections (scaled by relevant relic abundance fraction) for the two DM species plotted in the mass plane, with constraints applied from future experiment LZ [7]. Note that small values below the range are shown with colour corresponding to the smallest value.
We remark from Fig. 11(c) that points saturating the measured relic density exist in almost the whole allowed parameter region, thanks to the interplay between the two components a and ψ. In Fig. 12 we show the contribution of each specie to the total relic (left for ψ and right for a) in the (m ψ , m a ) plane. Interestingly, the region with m a ∼ m H /2 contains points with sizeable and dominant relic from ψ, while m ψ ∼ m H /2 is always dominated by relic from ψ. The remaining parameter space contains a region with m a ∼ m ψ where both species can receive competitive relic densities, and regions dominated by a for m a 300 GeV and by ψ for m ψ 1 TeV. Future direct detection experiments will be able to probe most of the remaining points, as demonstrated in Fig. 13, where we impose the projected exclusion by the LZ next generation experiment [7]. The surviving points consist of the Higgs sliver for a, with points dominated by the pseudo-scalar relic, and points with m a ∼ m ψ . The latter ones still have sizeable SI cross-sections, discernible from the neutrino floor at future direct detection experiments.
One should also note that the LZ experiment will be able to almost exclude the whole Higgs resonance region, ψψ → H, which can also be probed, independently, at future colliders via invisible Higgs decays. In this region a is heavier than m H /2, thus contributing a very small fraction to the relic density as shown in Fig. 12. This region of the parameter space can also be efficiently probed by searches for invisible Higgs decays, especially at the ILC that will have the strongest sensitivity. In Fig. 14 we present the comparison of the potential of the current LHC, HL-LHC and the ILC colliders to probe this loop-induced ψψ → H region. In the left panel we fix m ψ = 60 GeV and show the limits as a function of m a . It is remarkable that, for this mass point, the ILC will be able to probe the Higgs invisible decay close to the value corresponding to the ψ relic density saturating the PLANCK limit. The latter corresponds to Br(H → ψψ) 0.24% to be compared to the projected ILC reach of Br(H → ψψ) ≤ 0.4%. In the right panel, instead, we fix m a = 200 GeV and show the limits as a function of m ψ . We can see that the ILC will be able to completely exclude m ψ 59.5 GeV, while a region with the correct relic density will still be allowed for larger masses. Remarkably, the current ATLAS reach excludes m ψ 55 GeV, while the HL-LHC will be able to push the limit to m ψ 56.5 GeV.
To summarise, the viable regions of the parameter space for Scenario A are: • The aa → H annihilation region with m a m H /2 and λ aH 10 −4 , where the right amount of relic density is provided by the diagram in Fig. 10(c). This region can be probed by DM direct detection experiments and collider experiments looking for invisible Higgs decay. The main contribution to DM comes from a. • The ψψ → H annihilation region with m ψ m H /2 and Y 2 ψ λ aH > 1 and coupling generated at one-loop level. The dominant contribution to the relic density comes from ψ. This is the only region where a can be unstable, provided that m a > m ψ /2. This region can be effectively probed and even potentially closed by future ILC searches for invisible Higgs decay channels.

Scenario B: ψ FIMP regime with thermal a
As we have seen, small values of Y ψ O(10 −1 ) are excluded due to an excessive relic density of the fermionic component ψ. However, for extremely small values, Y ψ O(10 −8 ), ψ will not be in thermal equilibrium at early times and it will freeze-in by means of the scattering of a with the Higgs, aH →ψψ. On the other hand, sizeable values of λ aH would guarantee that a remains thermalised and contributes with a thermal relic component (as a second specie when m a < 2m ψ or by decaying into the fermionic DM).
In Fig. 15 we present the results of this regime for the m a < 2m ψ case, corresponding to twocomponent DM. The first two plots in the top row -Figs 15(a) and (b) -show the Ω a h 2 in the (m a , λ aH ) plane, bearing similarity with Figs 11(a) and (b) and demonstrating that the allowed regions are dominated by the thermal production of a. The only remarkable difference is the absence of "leaking" points, which were due to the co-annihilation processes involving ψ (so the smaller values of λ aH were allowed), which are now suppressed by the small value of Y ψ . The contribution of ψ via freeze-in is shown in the bottom frames of the figure. In Fig. 15(d), in particular, we show the relic density of ψ in the (m ψ , Y ψ ) plane. In this plot we can identify two distinct regions where sizeable values of Ω ψ h 2 can be attained (including saturating the whole DM relic density): one for Y ψ 10 −9 starting from masses m ψ 30 GeV (region BI), and another one for lower couplings, 10 −12 Y ψ 10 −9 , starting at m ψ m H /2 (region BII). These two regions can be better understood by looking at the complementary plane, (m a , m ψ ), shown in Fig. 15(c): the region BI corresponds to points where a is in the Higgs resonant sliver represented by the horizontal band; the region BII corresponds to triangle region at large a mass, where m a > m H . The scenario B can be probed only via the a component of the DM relic and at colliders: BI region is accessible via the Higgs invisible decay searches at colliders, while DM direct detection experiments would be mainly sensitive to the region BII, as one can observe from Fig. 15(b) demonstrating the effect of these searches.
Finally, in Fig. 16 we present numerical results for the region of the parameter space where m a > 2m ψ , region BIII, corresponding to a one-component DM (ψ) that originates from ψ freeze-in as well as a → ψψ decay processes after a freezes out. In general, the correct evaluation of Ωh 2 requires taking into account the fact that a may be long-lived due to the small values of Y ψ . The final relic densities stored in the two species are given by where Ω ψ (t F I ) is the ψ relic density at its freeze-in time and Ω a (t F O ) is the a relic density at its freeze-out time, which is typically much smaller than its life-time τ . In most of the parameter space, τ , is much smaller than the CMB time, τ t CM B 2 × 10 5 years, hence the relic densities in Eq. (5.18) simplify to  while the relic density of a is negligibly small.
There are several important features of the BIII region. One of them is that there is no sensitivity from DM direct detection experiments neither through a, as its relic density is negligibly small, nor through ψ that has a very weak coupling to the SM. This region can be only tested via the invisible Higgs decay search at colliders, as shown in Fig. 16(a). Figure 16(b) presents the interplay between λ aH ans Y ψ couplings: for λ aH 1, the relic density saturating the PLANCK measurement is mainly provided by freezed-in ψ with 10 −12 < Y ψ < 10 −11 . This feature is also clearly visible in Fig. 16(d) via the upper edge in the allowed m ψ values. When Y ψ < 10 −12 , instead, the main contribution to the relic density comes from λ aH 0.1 via the a relic at freeze-out, which then completely decays to ψψ. One can also see one more pattern in Fig. 16(b) represented by "scattered" points with 0.001 < λ aH < 0.1, where aa → H annihilation takes place and provides the right amount of DM via a freeze-out. Figure 16(c) shows the range of m a and m ψ masses viable in this scenario. In particular, it shows that the lower limit on m a is about 50 GeV. This limit comes from the current invisible decay search at the LHC, which extends the m a m H /4 limit which comes from relic density constraints defined by m a and m ψ kinematics, as one can see from Fig. 16(a). One can also see from Fig. 16(a) that the hierarchy between the masses m a and m ψ can be quite large. This means that a small value of the ratio m ψ /m a can provide viable parameter space even if Ω a (t F O ) is too large, as one can see from the second term of Eq. (5.19).

Scenarios C and D: 2-component FIMPs
These two scenarios are characterised by a very small coupling of the Dark sector to the SM, i.e. a tiny λ aH . As such, they are very difficult to test while they can provide the right amount of relic density. For this reason, we do not present any numerical scan, instead we will qualitatively discuss the main features of the two scenarios.
In case C, λ aH is very small while Y ψ is sizeable. Hence, both a and ψ can be produced via freeze-in via the couplings to the Higgs (for ψ loop induced). A large Y ψ would simply reshuffle the relic density of the two components at later times. When m a > 2m ψ , then a would promptly decay resulting in ψ saturating the relic density.
In case D, the smallness of both relevant couplings would lead to ψ occupying an insignificant part of the relic as freeze-in for this species would be hampered doubly by the small couplings. This means that in this scenario ψ decouples from a and the model is effectively reduces to the well-explored scalar portal model with freeze-in scenario.

Conclusions and outlook
After the discovery of the Higgs boson, the search for a Dark Matter particle has become the new grail and hard-sought nirvana of the particle physics community. The diversity in the experimental techniques and the remarkable progress achieved in each one call for more sophisticated theoretical studies, especially when trying to combine and compare various experimental bounds. The main difficulty stands in the large array of energies probed by the experiments: from low energy interactions in direct and indirect detection, to high energies at colliders like the LHC and the future FCC-ee/hh, ILC and CEPC. Moreover, additional constraints come from Cosmology via the relic density, precisely determined via the cosmic microwave background measurements, and from precision measurements in the electroweak sector of the Standard Model.
Complete models that contain a Dark Matter candidate, like supersymmetry or composite Higgs models, provide a consistent comparison at the price of specificities that are hard to disentangle from the phenomenology and generic features of the Dark Matter sector itself. Exploration of Dark Matter properties independent of quite a few details of the complete model is, therefore, a challenge. In this work, we propose a systematic classification of minimal consistent Dark Matter models, which are required to respect the complete symmetries of the Standard Model, as summarised in Table 1. They provide the missing link between effective field theory approaches and complete models, and offer a consistent and model-independent comparison between various experimental constraints. Moreover, because of their consistency, MCDM models can serve as a complete theory by themselves or be used as a building blocks within a bigger framework. This approach allows to create a convenient basis for the DM model space which can be used for a systematic DM exploration at various experiments.
In our framework the Dark Matter particle is embedded in an electroweak multiplet, characterised by its weak Isospin and hypercharge. Similarly, a mediator multiplet is included with all renormalisable interactions. The only exception to the latter is given by dimension-5 couplings to the Higgs, which can split masses and, therefore, crucially influence direct detection bounds. We consider in this work fermionic Dark multiplets and discover that many models are still allowed by all constraints, beyond the simplified models currently considered in the literature. We also revisit one-loop contributions to direct detection, including for the first time the mass splits in the dark multiplet. Due to the presence of a fine cancellation among various contributions to the amplitudes, the presence of a small mass split affects significantly the total spin-independent cross section.
Our main results can be summarised as follows: • Dark multiplets with hypercharge equal of above 1 are excluded by the presence of a charged lightest component and Z decay bounds.
• The loop-induced direct detection excludes multiplets with Isospin equal or above 3 (sextet), while other multiplets are probed by current or future experiments. The doublet escapes detection thanks to a cancellation among various contributions to the elastic scattering amplitude.
• Dimension-5 couplings of the Higgs, potentially generated by a heavy scalar or fermion mediators (i.e. described by MCDMs with a mediator multiplet), play a crucial role in splitting Dirac multiplets in Majorana mass eigenstates, hence removing the strong constraints from Z-mediated direct detection. On the other hand, the value of the mass split of the neutral states of the order of few GeV is being tested by DM direct detection experiments at present, while future DM direct detection experiments will be able to test it at sub-GeV level.
• Besides the important role of the mass split effects for DM direct detection, we have also shown the role of the PDF and the QCD scale uncertainties, which can be similarly important to provide cancellations for a loop-induced direct detection amplitudes.
We also study in detail a new model with a Dirac singlet Dark multiplet and a CP-odd singlet scalar mediator. While the mediator is even under the Dark parity protecting the fermion, it can be accidentally stable if it is lighter that twice the fermion mass. Thanks to the interplay between the two components, the parameter space can be probed by the synergy between future direct detection and the measurement of the invisible Higgs decay width at colliders. Furthermore, in the small coupling regime, both fermion and scalar can be produced in the universe by freeze-in. This is one example of interesting models, neglected in the literature, which is highlighted by our complete MCDM classification.
In this paper, we provide a first complete and consistent classification of effective models for Dark Matter that allows for a consistent and systematic comparison between all constraints on the Dark Matter particle candidate. We focus here on fermionic spin-1/2 Dark multiplets, while the paradigm can be applied to any spin option. We leave for a future work to compile a classification for spin-0 and spin-1 Dark multiplets and their minimal one-mediator extensions. Many models are still allowed and viable, beyond the simplified cases analysed in the literature. A systematic study of all the cases can help us establish the feasibility of a Dark Matter candidate around the electroweak scale, which seems to be under siege by the non-discovery of the historical WIMP candidates. Furthermore, as our classification requires full invariance under the Standard Model symmetries, the models we present can be easily embedded into more complete models and they can be further UV-completed in a consistent way. Henceforth, we believe that our classification can provide the required missing link between experimental searches and the model building required to obtain the new Standard Model that includes Dark Matter.
First we discuss the box diagrams (A and B in Fig. 5). For Dirac and pseudo-Dirac cases, D + and D − are distinct particles (D ± is used to refer to their respective antiparticles), which can have different masses. As such, a DM particle(antiparticle) will couple to up(down)-type quarks via diagram A containing a D + (D + ) or diagram B containing a D − (D − ) and down(up)-type quarks via diagram A containing a D − (D − ) or diagram B containing a D + (D + ), or to all quarks by diagrams A and B when two Z bosons are exchanged. In the Majorana case, the charged state in the two loops is the same. Hence, these diagrams can be combined when in the Majorana case, or when the masses are the same (as it is for Y = 0 with mass splits generated by EW loops).

A.3
Loop induced h-ψ-ψ coupling inF 0 0 S 0 0 (CP-odd) representative model Here we present details of the calculation the loop which induces h-ψ-ψ coupling, key to the phenomenology of the model studied in Section 5.  The loop-induced h-ψ-ψ coupling is generated by a loop containing the pseudo-scalar, as shown in Fig. 17. This loop must be evaluated at three different scales: for DM direct detection the external Higgs momentum squared is q 2 = t ≈ 0; for Higgs invisible decays to a pair of DM fermions, q 2 = s = m 2 H ; finally, for the relic abundance computation, the scale varies with the temperature, q 2 = s = 4m 2 ψ (1 + 1 2x ) where x ≈ 20 around freeze-out temperature (in this region the loop factor is relatively insensitive to x). It is useful to define the quantities Υ DD ≡ Υ(s = 0) , The amplitude for this diagram is given by. where p 1 (p 2 ) is the incoming(outgoing) momentum of fermion ψ and q is the incoming momentum of the Higgs boson. The functions C i correspond to the 3-point Passarino-Veltman integrals [133,134]. The derivation of this function using Package-X is provided at [132] and a c library used by a LANHEP implementation of this model are given at [123].   which is a function of a single variable, β ≡ m 2 a /m 2 ψ . The value of this function over the range of β is presented in Fig. 19.