Charming Dark Matter

We have considered a model of Dark Minimal Flavour Violation (DMFV), in which a triplet of dark matter particles couple to right-handed up-type quarks via a heavy colour-charged scalar mediator. By studying a large spectrum of possible constraints, and assessing the entire parameter space using a Markov Chain Monte Carlo (MCMC), we can place strong restrictions on the allowed parameter space for dark matter models of this type.


Introduction
The existence of dark matter (DM) has, since the its early days [1], been established through a wide range of detection techniques, such as galactic velocity curves [2][3][4][5], gravitational lensing [6], and its effects on Big Bang Nucleosynthesis (BBN) and the Cosmic Microwave Background (CMB) [7]. However the interactions of DM outside of its gravitational influence remain elusive, despite concerted efforts to measure its scattering in terrestrial targets (direct detection), its annihilation or decay products in the galaxy or beyond (indirect detection), or through its direct production in colliders [8].
One property of DM that is known to high precision is its abundance in the Universe today. The evolution of the structure of the Universe is well modelled [9] and so the starting point for building a model of a particle DM is to consider how its interactions influence its relic abundance. This leads to the concept of a thermal WIMP (weakly interacting massive particle), in which the DM achieves its relic abundance by decoupling from thermal equilibrium due to its annihilations or decay into standard model (SM) particles.
Under the assumption of a WIMP particle interpretation of DM, we have no concrete indications of its mass, spin or interactions, which leaves tremendous freedom when building models. Although many concrete models, e.g. supersymmetric theories, predict the existence of a DM candidate, so far these theories remain unverified and the phenomenology is often complicated by the large parameter spaces. This represents a top-down approach in which DM arises naturally from a UV complete model.
An alternative approach to DM model building is from the bottom up, where a class of simple low energy models or interactions are considered simultaneously. With no theoretical guiding principle, except gauge symmetry, on which to build such models, one must consider all possible models within a framework of a few assumptions. This is most easily done using a set of EFT (effective field theory) operators. Although an EFT may be perfectly valid for low energy experiments such as direct or indirect detection, they face problems with collider searches where the EFT approximation breaks down when heavy (TeV) states become energetically accessible.
To ensure the model is valid up to high energies and above the reach of colliders, a commonly used tool is simplified models, where often the mediator between the dark sector and the SM is included as a propagating mode. Simplified models arose first in the context of collider searches for missing energy [10][11][12][13][14][15][16], but have recently been applied more widely to indirect and direct detection [10,17,18], they allow for a much more broad study since the models themselves are sufficiently simple to contain only a few parameters which dominate the phenomenology of the DM. This approach is not without criticism, and can at times be too simple, for example neglecting gauge symmetries and perturbative unitarity [19][20][21].
Given the remarkable agreement between the SM and experimentally measured flavour observables it is natural for new physics (NP) models to enforce the minimal flavour violation (MFV) assumption to suppress large NP effects [22,23]. This assumption limits any quark flavour breaking terms to be at most proportional to the Yukawa couplings, which are responsible for the small violation of the flavour symmetry in the SM. This suppresses Flavour Changing Neutral Currents (FCNCs) and avoids strong constraints from rare decays and neutral meson mixing. Nonetheless, some such observables are not reproduced by SM calculations and hence allow room for violations of MFV, for example D 0 mixing which we discuss in section 3.1.
Some recent studies of simplified models have begun to go beyond the MFV assumptions. This has been done in the context of down-type couplings [24], leptonic couplings [25], and more recently top-like [26], or top and charm-like couplings [27]. Such models allow a continuous change from the MFV assumption to strong MFV breaking and can quantify the degree of MFV breaking permitted by the flavour constraints. Similar scenarios have been studied in [28], taking an overview of both lepton and quark flavoured DM and as well as a more focused study on top DM [29], both in the MFV limit.
Our aim in this paper is to extend the work of [26], taking a more general approach to these kinds of beyond MFV models -by placing fewer restrictions on the parameters of the model we include models with dominant up and charm type couplings, which give non-trivially different exclusion regions for different flavours of DM. We note that a similar scenario, except with scalar dark matter and a fermionic mediator has been studied in [30]. We aim to present statistically robust bounds from the entire parameter space based on a Markov Chain Monte Carlo (MCMC) approach.
We consider the following constraints in detail: • Relic Density (section 2): We calculate the relic density of all three DM particles, including their widths and important coannihilation effects.
• Flavour Bounds (section 3): We provide bounds on the model from neutral charm meson mixing, ensuring that the new physics does not exceed 1 σ of the experimental measurement of the mass difference between the heavy and light state of the D 0 . We assess the possibility for constraints on rare decays like D + → π + but find that the NP is relatively unconstrained compared to mixing.
• Direct Detection (section 4): We calculate the event rate for the most excluding DD experiments (LUX and CDMSlite) over a large range of DM masses, including all relevant contributions up to one loop order (including gluon, photon, Z and Higgs exchange) and matching to a full set of non-relativistic form factors.
• Indirect Detection (section 5): We include a large collection of constraints from the literature on the thermally averaged annihilation cross section σv for annihilation into various search targets such as photons, electrons, protons. We also include a study of gamma ray line searches, generated at the one-loop level in our model.
• Collider Searches (section 6): We perform a robust simulation of the dominant signals for a series of monojet, dijet and stop searches for ATLAS and CMS, including the widths of the particles.
We also compute constraints coming from electroweak precision observables, and perturbative unitarity. We calculate the Peskin-Takeuchi parameters [31,32], as these characterise the NP effects in much of the parameter space of our model, and replicate the literature result for a charged singlet scalar [33]. We find that the S, T, U parameters provide no additional constraints beyond those previously described, and similarly perturbative unitarity calculations prove to be unconstraining and so we make no further mention of them.
Including the various constraints named above we can carry out an MCMC scan in order to identify the parameter space left open to the model -our results are collected in section 7. We find that current data can be used to restrict the parameter space where DM of this kind can exist, and go beyond the results of [26] by showing how renormalisation group mixing and running can dramatically improve the direct detection constraints, disfavouring attempts to avoid these limits by predominantly coupling to top quarks.

The DMFV Model
The SM (without Yukawa couplings) has a flavour symmetry amongst the quarks -there are no flavour violating effects such as FCNCs at tree level. Minimal Flavour Violation (MFV) is then the statement that the only flavour symmetry breaking terms in the BSM model are the Yukawa terms [23].
In the model of Dark Minimal Flavour Violation (DMFV) originally proposed in [24], the SM quark flavour symmetry is increased by the inclusion of a U(3) symmetry in the dark sector, and the DMFV hypothesis is that this enlarged flavour symmetry is broken only by terms involving the quark Yukawas and a new coupling matrix λ. In the original work [24] λ coupled the DM to right-handed down-type quarks, whereas in this work we couple the DM to up-type right-handed quarks (the choice of right-handed quarks avoids having to introduce any non-trivial SU(2) structure). In this model, we introduce four new particles -a scalar φ that is colour and electrically charged, and a flavour triplet χ i that is a singlet under the SM gauge groups (which allows it to have a standard Dirac mass term). In table 1 we detail the behaviour under various gauge and other symmetry groups of the new particles and the coupling matrix -the transformation of λ under the U(3) flavour symmetries is to be understood in the sense of a spurion field [23]. The new physics Lagrangian reads giving the vertices shown in fig. 1. Note that a coupling between the mediator and the Higgs as well as a mediator self-coupling are allowed by the symmetries of the model, but we neglect them in this work. It was shown in [24] that coupling matrix can be written in the form Figure 1: Feynman rules for the interaction in eq. (2) with the matrices D λ and U λ parametrised as (defining c ij ≡ cos θ ij , s ij ≡ sin θ ij ) where θ ij ∈ [0, π/4] to avoid double counting the parameter space, and we require D ii < 4π for a perturbative theory. The presence of complex couplings (δ ij = 0) creates a violation of CP symmetry (note this is also permissible in the MFV assumption, so long as the complex phases are flavourblind [34]). Due to the stringent constraints from electric dipole moments (EDM) in the presence of CP violation [23] we will set δ ij = 0 throughout. In total we then have a 10 dimensional parameter space Other than those mentioned above, the only other limit we place on our parameters is m χ , m φ 1 GeV, so that the DM is a conventional WIMP candidate and the mediator is sufficiently heavy to decay to at least the up and charm quarks. Although the masses of the DM fields and mediator field are in principle arbitrary free parameters, one must impose m χ,min < m φ +m q (where m q is the lightest quark to which m χ,min couples) to ensure χ cannot decay. Similarly we must have m φ > m χ,min + m q , which ensures the mediator has at least one decay channel and prevents it obtaining a relic abundance itself.
It can be shown additionally that a residual Z 3 symmetry exists in the model [24,35], which prevents either χ or φ decaying into purely SM particles. This useful symmetry argument ensures the relic DM (the lightest of the three) is completely stable even once non-renormalisable effects are considered. It is possible for the heavier χ fields to decay to the lightest χ (DM) -in fact the rate of such decays are always large enough to totally erase the relic density of the heaviest two DM.
Finally, we briefly mention some interesting behaviour of the widths of our new particles. First, the mediator width Γ φ can be shown to be very narrow, with Γ φ /m φ ≤ 9 128π 1 % even in the limit of non-perturbative couplings. Secondly for small mass splittings (m χ i = m χ j (1 + )) the decay rate χ i → χ j + qq scales as 5 , which is important when we consider the relic abundance of the different DM species.

Relic Density with Coannhilations
As mentioned in the introduction, the relic density (RD) of DM is currently measured to a very high accuracy by the Planck collaboration [7], and this must be reproduced by any self-respecting DM model. We will assume that dark matter is produced thermally via a freeze-out mechanism, but the resulting constraints may be alleviated via non-thermal mechanisms as in asymmetric dark matter [36,37]. We leave this possibility to further studies.
In our model with three possible DM candidates, with potentially almost degenerate masses, we follow the results of [38] -Section III in particular deals with the effects of coannihilations (processes with χ i χ j → SM, i = j). In that work, the authors describe how coannihilations can be very important, and can be included in the "standard" computation [39][40][41] of relic density through the use of an effective annihilation cross-section σv eff , defined in eq. (12) of [38]. We will not reproduce all the detail from that paper here, but summarise the key results.
To compute the relic density, one first finds the freeze-out temperature x f ≡ m/T f by solving the equation with g eff an effective number of degrees of freedom of the near-degenerate DM candidates, M pl the Planck mass, g * the total number of relativistic degrees of freedom at freeze-out. The relic density itself can then be written where a ii and b ii are the s-wave and p-wave terms of σv ii (the cross section for the relic, plus any particles with degenerate mass), and I a,b are temperature integrals.

The Generation of Mass Splitting
Almost degenerate DM masses mean the mass splittings (∆m = m χ i − m χ j ) between the different χ i are important to determining the true value of the DM relic density. We can follow two regimes which distinguish the various possibilities by the dominant effect on the signals they generate: 1. The mass splitting is non-zero, the lightest of the χ i survives as the relic. This holds as long as the splitting is large enough to accommodate any kind of decay.
2. The masses are truly degenerate, equivalent to a degeneracy which is sufficiently small to prevent decay, i.e. ∆m ≤ 4 MeV. In this case, the three DM particles obtain equal relic abundances, with the total affected primarily by their coannihilations.
The difference between the effective cross-section method mentioned above and a full solution of the coupled Boltzmann equations, and the effect of degenerate masses is shown on the left of fig. 2. We see that the effective cross section approach correctly reproduces the relic density of the lightest candidate at late times, and that relic density constraints are not hugely sensitive to the mass splitting if it is non-zero. As the final relic density depends sensitively on whether a mass splitting in the candidates exists or not, we briefly talk about how such a splitting can arise. Splittings can arise from two sources -a tree-level contribution where m χ i and m χ j are split by mass terms of the form O(1) × (λ † λ) ii , or a loop-level contribution from renormalisation where the coefficient is instead of the order N c /(16π 2 ) log(µ 2 /Λ 2 ) multiplied by the tree level couplings (λ † λ) ii with Λ some high scale at which the masses are universal, and µ a low scale at which we wish to use the mass (e.g. for direct detection this could well be the nuclear scale of around 1 GeV). Explicitly, the resulting shift in the DM mass will be given by Note that because of our parameterisation of the coupling matrix, λ † λ is diagonal, with elements D 2 ii Relatively large splittings can be generated this way -with a high scale of 100 TeV, then the coefficient of (λλ † ) ii can be as large as ∼ 0.35. We explore the effect of mass splitting in our work by manually setting the mass splitting (∆m/m χ ) to a large (15 %) and small (2 %) value.

Mixing Observables
Since our model introduces couplings to the up-type quarks, we would expect new physics effects in the charm meson sector -in particular in neutral D 0 mesons. Mixing is observed in D, B , and K meson systems, and relates the theoretical quantities Γ 12 and M 12 to the observed decay width differences ∆Γ and mass differences ∆M between the heavy and light mass states of the meson. For D mesons, the current experimental averages from HFLAV are [42], On the theory side however, things are not so well developed. There are two possible ways to calculate the mixing parameters -inclusive, where we assume quark-hadron duality and sum quark level diagrams, or exclusive, where individual decay channels that contribute to D 0 mixing are calculated. In the exclusive approach (e.g. [43,44]), values of x and y on the order of 1 % are believed to be possible. However, currently exclusive D 0 meson decays cannot be calculated from first principles and the estimates in [43,44] were based on phase space arguments and SU(3) F symmetry.
On the inclusive side, we work within the Heavy Quark Expansion (HQE) formalism, see [45] for a review, assuming that the charm quark mass is large compared to the hadronic scale. For charm mixing the three leading dimension six contributions of the HQE suffer, however, from a huge GIM [46] and CKM suppression, leading a prediction that is orders of magnitudes below the experimental values, see e.g. [47], while the individual dimension six contributions are slightly larger than the experimental value. To decide whether the charm quark is heavy enough to apply the HQE one has to study observables that are not affected by any severe cancellations, a prime example for such an observable are lifetimes. First studies [48,49] have suggested that the HQE could hold with corrections of no more than 40 %. Assuming now the applicability of the HQE for the charm system we have to find a mechanism that is violating the severe GIM cancellation. In the literature three possibilities for such a breaking are studied. In [50] it was shown that a small breakdown (O(20 %)) of quark-hadron duality could enhance the predicted value of y up to its experimental value. An older idea [51] is that the GIM cancellation is much less pronounced for higher orders in the HQE. A first estimate of SU(3) breaking dimension nine contributions in the HQE gives x ≈ 6 × 10 −5 , y ≈ 8 × 10 −6 [52]still missing the experimental results by two or three orders of magnitude. Finally there is the possibility that the GIM suppression is lifted by new physics effects, which we will investigate. Because of these difficulties we have some freedom in the treatment of the SM contributions to ∆M and ∆Γ when constraining the allowed BSM contribution by comparison to experiment. One possibility [53] is to require that taking the 1 σ upper limit reported by HFLAV (eq. (8)). This is the limit that would be derived if the NP and SM contributions have roughly the same phase, so that since we know ∆M ≤ 2|M 12 |. The NP contribution to M 12 is given by where we take the decay constant f D from FLAG [54][55][56], the D mixing bag parameter B D from [57], and the loop function F is given by

Rare Decays
We consider the semileptonic decay D + → π + µ + µ − , whose short distance contribution comes from the quark level decay c → uµ + µ − . This decay is loop and GIM suppressed in the SM, and so should have good sensitivity to new physics. In our model contributions are no longer GIM suppressed, coming from electroweak penguin diagrams with our new particles in the loop.
Ref. [58] examines rare charm decays to provide limits on the Wilson coefficients of an effective theory -they look at D → µ + µ − as well as D + → π + µ + µ − and find the latter to place the strongest bounds for the coefficients relevant in our model. Matching onto their EFT, and neglecting the Z penguin since the momentum transfer is small, we find only the C 7 , C 9 coefficients are non-zero, corresponding to the operators (our full expressions for the Wilson coefficients can be found in appendix A).
Since the SM branching ratios for the D 0 decay suffer from a strong GIM cancellation, we would expect strong constraints on the flavour breaking terms of the DMFV model. As with the mixing observables, the rare decay process is primarily sensitive to (λλ † ) 12 in the limit of degenerate DM mass. On the right of fig. 3 we show the bounds coming from limits on the Wilson coefficients for (λλ † ) 12 = 1, 2, 4. The bounds on the individual Wilson coefficients are |C i | ∼ 1 (see Table II of [58]). Mediators up to m φ ∼ 50 GeV can be ruled out for couplings D ii ∼ (λλ † ) 12 ∼ O(1). These constraints are therefore substantially weaker than from meson mixing observables.
The rare flavour-changing decays t → u/cγ have been measured by ATLAS [59], but we find that the current limits are again not constraining on our model.

Direct Detection Constraints
Direct detection experiments are one of the most powerful ways of searching for DM, and operate by searching for DM scattering from atomic nuclei. The calculation of the scattering rate is done via an effective theory, where all heavy degrees of freedom (save the DM) have been integrated out, and then amplitudes are matched onto four fermion operators.
We choose to examine data from LUX [60,61] and CDMSlite [62], which together provide the best constraints over the range of DM masses we are looking at. LUX uses liquid xenon as a target, which detects DM with masses above 5 GeV while scattering from DM masses below this is kinematically impossible; CDMSlite is a germanium detector, and best constrains particles with masses between 1.6 GeV and 5.5 GeV. Details of our exact method can be found in appendix B -for now we merely state that we use a Poisson probability distribution for both, comparing the number of observed events in each bin to our predicted signal plus background.
At tree level, the only EFT operator which arises from our model is given by a diagram with t-channel φ exchange. We only consider the scattering amplitudes in which the incoming and outgoing DM (and quark) are the same flavour, as this avoids the computation of (possibly unknown) hadronic matrix elements of quark currentsq i Γq j for i = j. The operator in question is where the Mandelstam variable t has been replaced by its low velocity expansion and we have performed a Fierz transform (see e.g. [63]). Vector and axial-vector currents probe the valence quark content and spin distribution respectively of the scattered nucleon, and so would naively be small for non-valence quarks (i.e. c and t). However, there are 1-loop diagrams (see fig. 4) that mix operators with heavy quarks into those with up and down quarks, and in the case of heavy mediators RG running down to the direct detection scale (µ ∼ 1 GeV) also alters the relative coupling to nuclei. This calculation has been done in [64,65], and we find (see fig. 5) that DM that couples to heavy quarks at the mediator scale will mix into up quark coupling at the low scale with up to 10 % of its high scale coupling strength; tree level scattering is therefore substantial (as can be seen in fig. 6), even in the case of only coupling to heavy quarks. The spin-averaged cross section is parametrised by a series of nuclear form factors F (N,N ) ij [66], which are functions of the local galactic DM velocity squared where we sum over the form factors and the nucleons N, N = p, n. The nucleon coefficients above are related to our Wilson coefficients by where R ju (R jd ) gives the magnitude of the running of operator q j R γ µ q j R onto u γ µ u (d γ µ d ), and we have quoted the i = j = 1 relation since the corresponding form factor has the dominant scaling behavior. i and j run over the DM and quark flavours respectively. The dependence of the R jq parameters on the high scale (which we take to be the mediator mass) is shown in fig. 5.
At loop-level, there are various new operators that arise -in general these are highly suppressed, but we include them both because they can become dominant in particular regions of parameter space (see fig. 7) and for completeness. The operators we consider are photon operators [67,68] which in the non-relativistic limit correspond to the chargeradius, magnetic dipole moment, and anapole moment, Z penguins [67], and those for DM-gluon [69][70][71]. We reproduced the quoted literature results as a check.
The very latest null results from XENON1T [72] and PandaX-II [73] push the constraining potential of direct detection even further -nearly an order of magnitude stronger in cross-section, which translates into a factor of ∼ 2 in mediator mass.

Basics of Indirect Detection
Indirect detection experiments looks for signs of annihilating / decaying DM coming from astrophysical sources, typically the centre of galaxies where DM density is largest. The constraints are based around limits on the annihilation cross-section of DM to SM particles -in our model the main limits come from annihilation to quark pairs There is a bounty of possible search avenues for this annihilation signal; the energetic quarks will hadronize and decay into stable particles (photons, electrons, protons, and their anti-particles, which make up some part of the measured cosmic ray flux), which can be measured directly as they arrive at the earth (in the case of photons especially, which suffer very little energy loss to galactic or inter-galactic material), or indirectly through their influence on cosmic rays (for example photons produced by electrons/protons diffusing through the galaxy). We also have great freedom in where to look; generally anywhere where there is a cosmic overdensity of dark matter, close to home in the galactic centre or further afield in dwarf spheroidal (dSph) galaxies, galaxy clusters or the CMB. Underlying all these is eq. (17) and so ID constraints are frequently quoted as confidence limits on the thermally averaged annihilation cross section σv f f into fermions of the same flavour, covering a mass range m χ ∼ 1 GeV − 100 TeV. The ID signals from heavy quarks (q = c, b, t) are very similar (see Fig. 3 and 4 in [74]), and it is uncommon to find constraints on c, t final states (more common is the b). The primary spectra of electrons, positrons, anti-protons, deuteron and neutrinos are extremely similar between c, b, t quarks, and thus any constraints which look for these particles from DM annihilations will be approximately heavy-flavour independent. The situation is depicted in fig. 8.
It should be noted that the relative strength of these constraints is not robust, different authors use different halo profiles, different astrophysical parameters and are subject to varying degrees of uncertainty, some significantly larger than others, it is beyond the scope of this work to accommodate all these effects and compare constraints on a likefor-like basis and so what we present should be taken as representative but not precise. We will use the b b final state as representative for constraints based on dSph [75] and anti-proton measurements of AMS-02 [76] which dominate other constraints such as those based on other particle targets, such as the positron fraction [77] or neutrinos [78] and also those based on the galactic centre [79], or galaxy clusters [80].

Gamma rays (and other mono-chromatic lines)
At the one-loop level, the pair production of quarks from annihilating DM can pair produce photons at a fixed energy E γ = m χ /2 via a box diagram. We calculate this cross-section using an EFT where the mediator has been integrated out, in which limit only the axial vector operator (χγ µ γ 5 χ)(q γ µ γ 5 q) contributes to the s-wave annihilation, with cross section where s ≈ 2m 2 χ is the centre of mass energy of the annihilating DM, and C 0 is the scalar integral C 0 (0, 0, s; m 2 f , m 2 f , m 2 f ) in LoopTools notation [81]. As well as γγ final states, there will be γX final states where X = Z, h for example and these also provide constraints. The presence of a massive particle recoiling against the photon shifts the energy to E γ = m χ (1 − m 2 X /4m 2 χ ), but still creates a mono-energetic line signature. We show some results from the indirect searches in fig. 9 -we see that indirect searches can be quite powerful, especially in the case of large coupling to top quarks.

Collider Constraints
Our DMFV model contains a new particle with colour charge, and so we expect there to be significant limits coming from collider experiments. In addition we also have DM which can be searched for in final states with missing energy, and current LHC data can also place limits on the mass of invisible particles. In the past, DM model builders have used effective field theories (EFTs) to analyse NP at colliders, but in recent years it has become clear that the regions of validity of these EFTs at high energy machines such as the LHC are so small as to be almost useless [10,17,18]. We briefly detail in the next section this point for our particular model, before moving on to a more complete analysis.

EFT Limit
In [82] the validity of the EFT approximation for t-channel mediators is quantified by R Λ , which they define as the ratio of the cross section with the constraint t < Λ 2 applied to the total cross section (i.e. the total proportion of the cross section which is valid under the EFT assumption). The lines of R Λ = 0.50 are plotted alongside the EFT limits taken from ATLAS [83] (the R Λ contour assumes |η| < 2 and p T < 2 TeV, the ATLAS results assumed the same range of η, but allow p T 1.2 TeV). It is worth noting that the authors of [82] produce results with the limit g 1, the bounds become significantly weaker by using g 4π which then permit a small region of validity as shown in fig. 10. The EFT breaks down entirely for g 1. Thus the EFT approximation cannot be justified in our analysis and we turn to the simulation of the full cross section.

LHC bounds
To try and cover a large range of constraints, we look at three different LHC processes that could place limits on our model -monojet with missing energy searches, where a single jet recoils off DM pair production; dijet searches with missing energy; and stop searches. The latter are relevant to our model as we have a coloured scalar coupling to top quarks and DM, in analogy with the e.g. stop-top-neutralino vertex in many supersymmetric theories, and provide sensitivity to the φ-t coupling D 33 .
In fig. 11 one example Feynman diagram that generates monojet and dijet signals is shown -in the dijet case the decay of the mediator into quark plus DM is not shown. Other diagrams that contribute can be seen in appendix C.
We produce our collider constraints using MadGraph [84], replicating, except where noted below, the experimental cuts used by the experiments.

Monojet searches
In our analysis, we use the most recent monojet search by ATLAS [85] (which uses the Run 2 data ( √ s = 13 TeV and L = 3.2 fb −1 )), along with a similar analysis performed by CMS [86] with the Run 1 data ( √ s = 8 TeV and L = 19.7 fb −1 ). The total cross section as a function of m φ for a benchmark scenario is shown in fig. 12 with the ATLAS limits overlaid, and the constraints on our model are shown in the top of fig. 13.

Dijet searches
Moving on to dijet searches, we use a Run 1 and Run 2 search by ATLAS [87,88] looking for multiple jets plus missing energy -we restricted our comparison to the 2jet searches which should provide the strongest constraint. In our model, the process pp → φφ →χχjj provides the dominant contribution to this signal. We replicate all the main selection cuts for both analyses, in particular for the Run 1  Table 3: The four relevant signal regions from [89] and the cuts we have implemented.
comparison: E miss t > 160 GeV, p T,(1,2) > 130, 60 GeV, ∆φ > 0.4 (between the jets and missing momentum), and for Run 2 similar cuts are applied (full detail in Table 2 of [88]). The different signals regions (tjl, tjm, tjt) also include a minimum requirement for m eff and E T / √ H T , which are defined as which we implement in MadGraph manually via Fortran code (again, see the respective papers for the cuts in each case). The constraints this places on our model parameters are shown in the bottom left of fig. 13 for the case of no mixing and strong couplings for all DM particles.

ATLAS 2014 Stop Search
Lastly, a study by ATLAS [89] considers a set of cuts optimized for the detection of stops -the signal consists of a lepton in the final state along with four or more jets. There are four relevant signal regions tN_diag, tN_med, tN_high, tN_boost, each requiring a single lepton with p l T > 25 GeV, and cuts in table 3. 1 We find that the production of the φ pair is dominated by t-channel χ exchange and s-channel gluons; the photon and Z mediated diagrams are neglected. We calculate in MadGraph the cross-section for a single final state ((b b)(d u) + e − ), and then multiply this by four to account for the different top quark decay options (the p T cut means the different masses have a negligible effect). Although the cross section is predominantly controlled by the size of D 33 , the light quark couplings D 11 , D 22 have a mild affect by reducing the branching ratio φ → tχ i and hence suppressing the cross section.
We also examined constraints from a similar ATLAS search for scharms [90] rather than stops, searching for c-tagged jets plus missing energy in the region where the branching ratio φ → cχ i is large. The limits on m φ,χ are similar to the stop search, and thus do not warrant further attention when compared to the dijet searches. 1 We do not include the cuts on the parameters am T 2 and m τ T 2 . From the published cut flows it can be seen that the effect of these cuts is of the order 10 % and 2 % respectively (although the former cut can have a more pronounced effect ∼ 30 % on the tN_med cut choice).

Collider Constraints within DMFV
We have now looked at three classes of analysis: monojet searches, dijet searches, and searches optimised for a stop. Within our model we have couplings to u, c, t (which we denote here by λ u,c,t ) and the relative strengths of these dictate which signals will be dominant. Compared to λ u , the monojet and dijet processes are suppressed by pure λ c (due to the charm parton distribution function (PDF)), but generally are enhanced by mixtures of λ u,c . The coupling λ t reduces the signals since they dominantly come from s-channel φ resonances and thus the branching ratio to u, c jets is ∝ (D 33 ) −2 if λ t λ u,c . The stop search only becomes relevant for large λ t with λ t /λ u,c > 1, and increasing λ u,c suppresses the signal as the branching fraction to top quarks is reduced.
• Mostly up-type: The dominant signal will come from the monojet processes which have the least QCD suppression and which require an up quark in the initial state. Dijet searches are also sensitive but it tends to be the monojet which sets the better constraint.
• Mostly charm-type: The monojet processes are enhanced by the presence of charm couplings, however as the up coupling is reduced the monojet processes become suppressed by the charm PDF by around a factor 10-100. The dijet processes are very similar as for u quarks but the largest contributing diagram is again suppressed by the charm PDF. Both searches provide constraints.
• Mostly top-type: The monojet signal depends primarily on λ u,c , only indirectly on λ t though the widths. λ t can be probed through stop searches with jet multiplicities of ≥ 4.
Colliders provide very powerful exclusions (up to the TeV scale in mediator mass), and cover the full model parameter space in coupling, although these can be significantly weakened by, for example, strong top couplings. The DM is produced on shell, and so the constraints are comparatively weak at high DM mass when compared with searches which depend on the cosmic abundance of DM; on the other hand the fact that the DM is produced in the collider releases any dependence on its abundance in the universe, thus allowing more powerful constraints on DM which has only a fraction of the full relic abundance (or none at all). Similarly, low mass DM is strongly excluded, whereas the most powerful astrophysical probe (direct detection) cannot detect much below the GeV scale due to kinematics.
When compared with the strongest direct detection limits, the collider limits are not as constraining, and this is not likely to change even with more luminosity and higher energy beams.
It is very difficult for a given parameter choice to determine the strongest bound from colliders, except in the extreme cases above, and one should therefore check all available searches as we have done. Due to the interplay between 1 and 2 jet processes, there is no obvious scaling behaviour of the cross section with the coupling parameters, these factors Parameter Range Prior Table 4: Allowed ranges for the parameters used in the MCMC scan, along with the assumed prior likelihood, which is uniform on either a linear or logarithmic scale. make implementing collider searches in an MCMC scan difficult and slow as each cross section must be numerically computed at each point in phase space.

Results
We have aimed to produce a robust statistical analysis of the eight dimensional parameter space of the DMFV model, using the Bayesian inference tool MultiNest [91][92][93] and its Python interface PyMultiNest [94] with 5000 live points. The motivation for carrying out this analysis is twofold, firstly from a practical standpoint it enables very quick and efficient algorithms for scanning a large dimensional parameter space, allowing us to include all parameters in one analysis. Secondly, a rudimentary "hit-or-miss" analysis leaves a large region of parameter space allowed, which is not surprising given the flexibility of 8 free parameters, with a statistical result we can quantify the regions of parameter space which are allowed but very improbable given the errors of the experimental data. For clarity, we represent the allowed parameters as contours containing credible regions, using the method in [95]; using the posterior probability density function. The 1, 2 σ contours give an indication of the allowed parameter range, with containment probabilities of 68 % and 95 % respectively. Regarding the use of priors: We make one note of caution regarding the results; the credible regions depend sensitively on the choice of priors for the parameters. This is not surprising since our constraints allow large regions of parameter space to be equally well allowed, and so the use of priors which bias the parameters to lower values (i.e. log-uniform compared with linearly uniform) is reflected in the final result. Nonetheless, we are careful to limit the statements made in the text to those which are independent of the choice of priors. In all figures the log-uniform priors have been used for the masses and for D ii , as this represents the more conservative choice. The ranges and priors for the parameters of the scan are summarized in table 4.
Our results are summarized in figs. 14 to 16 as 2 σ contours, and in table 5 as onedimensional 1 σ intervals. We consider three separate samples in which the DM (the lightest χ) is the first, second and third member of the triplet (denoted 'up', 'charm' and 'top' DM). Within each sample we present a low and high mass splitting (2 % and 15 %), which primarily distinguish the effects caused by coannihilation in the calculation of relic density, but affect all other bounds to some extent as we have explicitly included the masses in each.
As we see from see fig. 14, the masses of the DM and mediator are both required to be in the TeV range, with upper limits in the tens of TeV, The DM and mediator masses are strongly correlated with the D ii , as in fig. 14, due to the relic density and mixing bounds which both scale approximately as (D/m) 4 in the high mass limit. Masses in the TeV range favour the D ii to be O(1). The mixing angles are not well constrained in general; θ ij = 0 is favoured, but the full range of angles are usually allowed with 2σ credibility.
The D ii themselves are highly correlated from the mixing constraints (see figs. 15 and 16) which depend on (λλ † ) 12 which is approximately where s ij = sin θ ij and so we see D 11 ∼ D 33 (and less strongly D 11 ∼ D 22 ). Because the correlation between D 22 , D 33 is less pronounced, the RD bound controls the behaviour and produces an anti-correlation, since the annihilation cross section scales like due to coannihilations, as such the trend is most pronounced for small mass splitting. This is seen in the range of D 22 for the small splitting data, fig. 16.   In all cases, increasing the mass splitting reduces the available parameter space of the masses and couplings of the DM since the coannihilations and annihilations of the heavy particles have a reduced effect on the relic density (scaling with a Boltzmann factor exp(−∆m)). This allows less flexibility in the DM parameters whilst potentially opening up the allowed parameters of the heavy particles, since their couplings are out of reach of the astronomical constraints (indirect and direct searches) which are proportional to the relic density of the lightest χ (scaling as Ω 2 and Ω respectively). This effect can be clearly seen in the right panels of fig. 14, where the 2 % splitting allows much smaller DM couplings compared with the 15 % splitting, contrastingly in fig. 15 (middle right panel) the non DM coupling space opens up with a larger splitting. Of course, since we have fixed the mass splitting by hand, the heavy particle parameters are not totally free, and so the parameter space is still reduced by the constraints we consider.
Top quark threshold effects are absent in the MCMC scan, due to the high masses (m χ m t ). Since m χ , m φ m t the three quarks are kinematically equivalent, and so the bounds are not strongly dependent on the flavour of DM. The main differences arise due to the quarks SM interactions which impact the DD and ID limits.
As described in section 6, we have studied collider bounds on our model, but these were not directly incorporated into our MultiNest routine as these bounds are much more computationally intensive than the others. However, as we see from fig. 13, the collider bounds only rule out sub-TeV scale masses, even at large couplings and so we do not expect that a full likelihood function incorporating the LHC constraints would give significantly different results. As a test, we checked a sample of the points inside the 68 % (1 σ) credible regions and found only a small minority (of order 1 %) that would be excluded by collider data. We produce, for each parameter, a marginalized posterior integrated over the remaining 7 parameters. From this distribution we find the 1 σ credible interval. The results are shown in table 5. This contains results for both uniform and log-uniform priors on D ii , m χ and m φ ; when the two cases are discrepant by > 1 σ this is due to a flat posterior, and using the 2 σ band instead the two agree.

Constrained Scenarios
We consider two extensions to the previous results: 1. In section 2.2 we found that the mass splitting which is generated through RG running of the DM self-energy is approximately proportional to D 2 ii , this motivates us to consider a scenario in which the couplings D ii are correlated with the masses (thus introducing a coupling splitting ∆D ii /D χ ∝ ∆m ij /m χ ). The reduced parameter space enforces almost degenerate couplings which leads to two important effects; firstly, it subjects all three χ to the astrophysical constraints of indirect and direct detection, despite the heavier particles having no relic density. By this we mean that, upon fixing the mass splitting, any limits on the coupling strength of the relic particle are translated to restrict the non-relic particles. Secondly, because the D ii are equal the mixing effects are naturally small and as a result the mixing angles are much less constrained as they do not need to be small to counteract flavour effects. This scenario is representative of a model in which MFV is broken only slightly, since the couplings to quark flavours are roughly equal, differing due to the mixing angles and the small differences in the D ii . It is actually only slightly less constrained in both mass and couplings than models in which flavour violation is allowed, which counteracts the naive assumption that without MFV, flavour observables restrict NP very high scales (O(100 TeV)).
2. When compared with the down-type quark sector, flavour bounds are weaker due to D 0 being less well measured and our conservative treatment in which we assume the SM contribution to D 0 mixing is zero and the experimental value comes entirely from the new physics. This is not entirely unreasonable, since short distance calculations of the observable are known to be very discrepant, nor is it completely reasonable, since long distance calculations are able to bring the SM into agreement.
To cover this caveat we consider a future scenario in which the SM calculation reproduces the experimental number (but the precision of the measurement stays at its current value). This is also conservative, since any interference terms between the SM and DMFV amplitude are likely to be large. The constraints on the mixing angles are more pronounced Results for these two further scenarios are shown in fig. 17, and the 1 σ intervals in tables 6 and 7.   are shown as lines and shaded regions respectively, the modal average of the log-uniform prior choice is shown as a dot. The two mass splitting cases are contained in different panels.

Conclusions
In this work, we have analysed a model of dark matter, based on [24] but coupling to up type quarks, that goes beyond MFV in order to allow potentially large new effects in the flavour sector, and have seen how the combination of a wide range of constraints can be used to place limits on models of this type. We approached this task of combining many different constraints using the MCMC tool Multinest, which allowed us to place limits on the high dimensional parameter space of our particular model. As we can see from fig. 14, the MCMC places lower bounds on the new particle masses of at least 1 TeV for Top DM, and a few hundred GeV for Up and Charm DM in certain cases. Our collider bounds ( fig. 13) cannot further exclude Top DM, even in the case of strong couplings, but could remove a small area of allowed parameter space from the bottom end of the mass range in the case of Up/Charm DM.
Ref. [26] considers this model, but examined the region of parameter space with dominant top quark couplings. Our results in general agree with their conclusions if we look at their more focused parameter space. For example, they find strong constraints on θ 12 except in the case of some degeneracy in the D ii , which we replicate. Similarly the strong constraints on DM mass from relic density and direct detection are reproduced. In their work, they explain how loop-level diagrams contributing to direct detection favour the dominant top coupling -however as we explain in section 4, RG effects mean even when DM doesn't couple to up quarks directly, the mixing is substantial enough to weaken this conclusion (as long as the mediator mass is large enough).
Given the current level of data, the model we examine of flavoured DM coupling to up-type quarks has large sections of its parameter space still allowed, so long as one considers large mass new particles. However, even without the complimentary collider results, the lower mass, phenomenologically interesting, regions of parameter space are disfavoured by flavour, relic density, and direct detection considerations.
The MFV assumption is frequently invoked in simplified models in order to evade potentially large flavour-violating effects. The level of robustness of this assumption varies considerably between up-type and down-type quark couplings in the DMFV model; for RH down-type quarks strong flavour bounds do ensure that the assumption is a good one. However for couplings to RH up-type quarks we have seen that in fact the flavour bounds are avoided in a large region of MFV-breaking parameter space.
One particular future development could alter this picture however -if a precise theoretical prediction of D 0 mixing observables could be obtained then either (a) a significant discrepancy requiring new physics is present, or (b) the SM predictions are reproduced with a high precision. The former would motivate the exploration of models which go beyond MFV, and the latter would make the MFV assumption a necessary assumption of the DMFV simplified model if one wants to avoid some fine-tuning. The non-zero Wilson coefficients arise from electroweak penguins (shown in fig. 18), and neglecting Z penguins since the small momentum transfer means they amount to an O(1%) correction, we find where B and C are loop functions using LoopTools [81] notation.

B.1. LUX
For situations where we have both a measured event count, N obs k (binned into energy bins labelled by k) and theoretical background N bck k , we can use the likelihood ratio test, a method based on a hypothesis test between a background only, and background+signal model, with likelihoods L, L bck respectively [96]. The likelihood of observing the data, D, assuming a particular set of parameters {λ}, is denoted L(D|{λ}). The likelihood of each bin is a Poisson distribution Poiss(N obs , N th (λ)) where N th k are the predicted number of signal events (including background), where N th (λ) = N DM (λ) + N bck . The background only model is identical but with N th = N bck . Then the test statistic, follows a χ 2 distribution -the cumulative probability density function of χ 2 (x) represents the probability that we observe the data given the model parameters λ. The value of x such that χ 2 (x) = C (i.e. the C % confidence limit) depends on the number of parameters {λ} -for only one parameter for example one can look up that χ 2 (2.71) = 0.9, which means that the 90 % confidence bounds on λ are given by TS(λ) = 2.71.

B.2. CDMSlite
For CDMSlite, we use a conservative method based on the statement that the 90 % confidence limit is such that there is a probability of 0.9 that if the model were true, then the experiment would have measured more events (n) than have been measured (n obs ).
This equation is numerically solvable for µ giving a required signal µ = 109 +51 −50 , 88 ± 14, 635 ± 37 and 207 ± 20 events for energy bins 1 to 4 respectively. This is conservative since a large portion of the measured events are background, and the resulting limits are slightly weaker than those given by the CDMSlite collaboration.

C.1. Monojet processes
The dominant diagrams contributing to the pure monojet process. Each processes scales as σ ∝ (λλ † )α s and can become extremely large for large λ. The cross section is dominated by the diagrams containing a heavy φ resonance.
u,c u, c χ i χ j φ Figure 19: The above diagram must include initial/final state radiation from external legs or internal bremsstrahlung from the mediator. The contribution is roughly equal amongst these emissions.  Figure 22: The left (right) process has σ ∝ (λλ † ) 2 (α 2 s ) and so the dominance depends on the size of the new couplings -for couplings which are large enough to be excluded it is usually the left diagram which dominates.
The dominant processes contributing to the production of on-shell φ, which decay φ → q i χ j producing a dijet signal. In monojet analyses, this provides a subdominant contribution compared with pure monojet processes figs. 19 and 20 in most of the parameter space.