Dark Matter and Gauged Flavor Symmetries

We investigate the phenomenology of flavored dark matter (DM). DM stability is guaranteed by an accidental ${\mathcal Z}_3$ symmetry, a subgroup of the standard model (SM) flavor group that is not broken by the SM Yukawa interactions. We consider an explicit realization where the quark part of the SM flavor group is fully gauged. If the dominant interactions between DM and visible sector are through flavor gauge bosons, as we show for Dirac fermion flavored DM, then the DM mass is bounded between roughly $0.5$ TeV and $5$ TeV if the DM multiplet mass is split only radiatively. In general, however, no such relation exists. We demonstrate this using scalar flavored DM where the main interaction with the SM is through the Higgs portal. For both cases we derive constraints from flavor, cosmology, direct and indirect DM detection, and collider searches.


I. INTRODUCTION
The stability of dark matter (DM) is commonly assumed to be due to an exact discrete symmetry, Z n . This can either be imposed by hand or have a dynamical origin. Examples include R-parity in the MSSM [1], and flavor symmetries in the leptonic sector [2][3][4][5]. Here, we explore the intriguing possibility raised in Refs. [6,7] that the stability of DM is due to the Yukawas. The lightest neutral state that is odd under Z 3 is therefore stable and is a DM candidate. This is the idea behind MFV dark matter [6][7][8].
Requiring MFV is sufficient, but not necessary. In this paper we formulate a general condition for flavored DM using flavor triality (see Eq. (3) below). For example, any spurion in the bifundamental of G SM Using this renormalizable model we show below that a thermal relic DM can be in a nontrivial representation of G SM F . There are two conflicting constraints on this setup. On the one hand, Flavor Changing Neutral Current (FCNC) constraints impose lower bounds on the masses of states in nontrivial flavor representations. On the other hand, a DM relic density consistent with observations requires large enough DM annihilation cross section so that some of these same particles need to be sufficiently light. Both of these requirements are satisfied for O(TeV) DM mass. This is low enough that it may be tested by direct and indirect DM detection experiments and searched for at high-energy particle colliders.
While the phenomenology of flavored DM models can be found in Refs. [7,[10][11][12][13][14][15][16][17][18][19][20], the construction of an explicit renormalizable model with inclusion of flavor-gauge interactions is new. Within our framework, the constraints on DM are more severe compared to a generic Effective Field Theory (EFT) analysis [6,8]. In particular, the flavor constraints from new fermionic states, and the fact that the vacuum expectation values (vevs) of the flavon fields need to reproduce the quark masses, makes the structure of the theory much more rigid and predictive.
The paper is structured as follows. In Sec. III we derive general conditions for DM to be stabilized by the exact accidental flavor symmetry of the SM (flavor triality). An explicit realization of this possibility is introduced in Sec. III in the form of a model with fully gauged G SM F . The resulting DM, flavor, and collider phenomenology is analyzed in detail in Sec. IV. We summarize our conclusions in Sec. VI, while more technical details of some of our computations are relegated to the Appendices.

II. STABILITY OF FLAVORED DARK MATTER
We start by formulating the general conditions required for flavored DM to be stable due The G SM F global symmetry is broken by the SM Yukawa terms is identical to a subgroup of U (1) B . This is no longer true in the presence of NP. In MFV for instance, Z U DQ 3 remains exact, while U (1) B can be broken, e.g., by dimension-9 operators [22] (see also [7]).
The Z U DQ 3 may be the underlying reason for the stability of DM. To make this explicit it is useful to introduce Z χ 3 , a diagonal subgroup of Z U DQ 3 × Z c 3 . Here, Z c 3 is the center group of color SU (3) c , under which {U R , D R , Q L } → e −i2π/3 {U R , D R , Q L } . All the SM fields are thus Z χ 3 singlets. In MFV NP Z χ 3 is exact, so that the lightest Z χ 3 odd particle is stable and can be a DM candidate [7].
We generalize this observation beyond MFV. To this end, we introduce the notion of flavor triality [7]. Consider a field X in the G SM F representation X ∼ (n X Q , m X Q ) × (n X U , m X U ) × (n X D , m X D ), where n X i , m X i are the Dynkin coefficients of the corresponding SU (3) i group. We call flavor triality the quantity (n X − m X ) mod 3, where n X = n X Q + n X U + n X D and m X = m X Q + m X U + m X D . The basic requirements for flavored DM to be stable due to Z χ 3 are the following. First of all, G SM F needs to be a good symmetry in the UV. Secondly, G SM F needs to be broken only by spurions Φ with zero flavor triality, (n Φ − m Φ ) mod 3 = 0. This ensures that Z χ 3 is unbroken. (The spurions Φ need to be color singlets in order not to break color.) The lightest Z χ 3 odd state is then stable. If it is a color singlet it is a potential DM candidate. This also means that DM is in a nontrivial flavor representation with nonzero flavor triality, (n χ − m χ ) mod 3 = 0.
The above shows that models with flavored DM can deviate significantly from MFV.
In particular, Z χ 3 is not broken by a vev of any field that is in an adjoint or in a bifundamental of G SM F . Specifically, any function f (y u , y d ) leaves Z χ 3 unbroken. More generally, additional flavor-breaking sources that transform as (8, 1, 1), (1, 3,3), . . . may be present without spoiling DM stability. While the flavor structure of such NP models is not of MFV type in general, the stability of DM is still a consequence of an unbroken flavor subgroup.
DM is in a nontrivial representation of the flavor group, leading to distinct phenomenology depending on the nature of the flavor breaking and on which flavor multiplet χ belongs to.
An important starting point in the above discussion was the assumption that G SM F is a good symmetry in the UV. This is most easily achieved, if G SM F is gauged. We explore this possibility in the remainder of the paper.

III. GAUGED FLAVOR INTERACTIONS AND DARK MATTER
We gauge the full SM quark-flavor symmetry G SM F . The fermionic sector is extended to cancel the anomalies of the new gauge sector. We use the model of Ref. [9] that allows for O(TeV) flavor gauge bosons (FGBs). The SM Yukawas, y u,d , arise from the vevs of new scalar fields transforming as under G SM F . The minimal set of chiral fermions that ensures anomaly cancellation of the new gauged sector is where the index L and R represents their chirality. Together with the SM fermions they, therefore, form vector-like representations of G SM F . The SM gauge quantum numbers of Ψ uR , Ψ dR , Ψ uL , Ψ dL are the same as for U R , D R , U R , D R , respectively, i.e., they are SU (2) L singlets but charged under U (1) Y . Because the new fermions are vector-like under the SM, e.g, Ψ uR transforms like Ψ uL under the SM, all anomalies in the SM sector cancel.
Remarkably, with the above fermionic content also all mixed gauge anomalies cancel. In fact, one could also gauge two additional flavor diagonal U (1)'s, U (1) U and U (1) D , as well as U (1) B−L , a possibility that we do not pursue further, but is discussed in Ref. [9].
The Yukawa and relevant mass terms in the Lagrangian are [9] where λ ( ) u,d are flavor-universal coupling constants and M u,d flavor-universal mass parameters. The mass terms in Eq. (6) mix the states Ψ uL,uR and U L,R forming mass eigenstates u i and u i , where i = 1, 2, 3 is the generation index (and similarly for down-quark states). After electroweak symmetry breaking the masses for the two mass-eigenstate sets are, in the limit [9,23] The mass matrix for the FGBs, A a A , A = Q, U, D and a = 1, . . . , 8, is governed by the vevs of the Y u,d scalar fields and the gauge couplings, [9,23] with λ a the Gell-Mann SU (3) The SM Yukawas, y u,d , are non-analytic functions of the spurions Y u,d , which signals that the theory is strictly speaking not MFV. Analogously, the NP states, u i , d i and A m , have masses that are non-analytic in terms of the SM Yukawas. However, the low-energy observables, with only the SM fields on the external legs can be MFV-like. If the M u,d / Y u,d suppressed terms are kept in Eq. (9), the y u,d become more involved functions of Y u,d −1 .
These are analytic in Y u,d −1 since the effects of NP states decouple in the Y u,d → ∞ limit. Similarly, the NP contributions to the low-energy observables C i take the form with F,F analytic functions. One can thus expand where we assumed for illustration that the transition is due to the left-handed quark current. As long as there are no large flavorconserving ratios, i.e., as long as Taylor expansion can be truncated after a few terms (see Ref. [24] for a more detailed discussion). In this limit, the low-energy effects are of the MFV type, suppressing the FCNCs to acceptable levels already for NP states at the electroweak scale. In a numerical analysis that we perform in Appendix A, we find that the expansion of the effective weak Hamiltonian in terms of SM Yukawa couplings can indeed still be performed for FGB contributions.
Since Y u,d are in the bifundamental representation of G SM F , the Z χ 3 remains unbroken. As argued above the Z χ 3 can be used to make flavored DM candidates stable. We consider two examples: i) a fermionic DM in a vector-like representation of G SM with the visible sector through FGBs, and ii) a scalar flavored DM, that interacts with the visible sector by exchanging FGBs as well as the Higgs.

A. Fermionic flavored dark matter
The DM in the first model is a massive Dirac fermion, χ, in a vector-like representation so that no gauge anomalies are induced. Its mass term is Since χ is charged under Z χ 3 , the lightest member of the χ triplet is stable. Note that we could also gauge a larger global group G SM That we identify SU (3) χ with SU (3) U is a model-building choice.
The DM flavor triplet, χ, is split by radiative corrections due to the exchanges of FGBs, see Fig. 1. In the m 0 χ m A limit, the DM mass splitting at one-loop is given by where ∆m χ is a 3 × 3 matrix and so is Ξ ≡ λ a (log M 2 A is given in Eq. (8), while the a, b indices run only over the eight SU (3) U generators. The unphysical µ dependence cancels in the r.h.s. of Eq. (12). The χ i , i = 1, 2, 3, mass eigenstates are obtained by diagonalizing the mass matrix ∆m χ . And Ξ is a function of  In the numerics we use the exact one-loop expressions for the DM mass splitting, so that the lightest state χ 1 is split away from χ 2 and χ 3 , with the latter approximately degenerate, m χ 2 m χ 3 . This is very different from MFV DM [6][7][8] where the DM mass splitting is assumed to be expandable in the SM Yukawa couplings. In that case one has an approximate U (2) symmetry for the first two generations giving m χ 1 m χ 2 , while the topflavored DM, χ 3 , is split away from the first two generations, and can be either significantly heavier or lighter.
The relation ∆m 32 ∆m 31 signifies that our flavored DM is non-MFV. The flavor gauge group SU (3) U is broken by the FGB vevs Y u . This breaking is larger in the first two generations. Since the quark masses are inversely proportional to Y u , this leads to an approximate global U (2) U symmetry in the quark sector. Such an approximate symmetry is not found for the radiative corrections to DM masses, m χ i . The DM multiplet has a chiral symmetry in the limit m χ → 0 ensuring that the radiative corrections are proportional to m χ and only log-dependent on FGB masses. The splitting does not vanish in the Y u → ∞ limit (or, equivalently, y u → 0 limit), since in this limit the SU (3) U gauge group is still completely broken. Numerically, for m χ ∼ TeV the splittings are ∆m 31 ∼ O(10 GeV) and ∆m 32 ∼ O(1 GeV) for g U ∼ 0.4 and can be less than a pion mass for an order of magnitude smaller g U .
The DM multiplet can be split more significantly if there is flavor breaking beyond Y u , Y d . As an example we consider an additional scalar field in the adjoint of SU (3) U , Φ U ∼ (1, 8, 1). The DM mass Lagrangian now reads and yields DM masses that are split already at tree level, ∆m χ = λ χ Φ U . We assume that Then the two diagonalize in the same basis giving O(1) splitting between all three members of the multiplet. The alignment is not needed in general, but does simplify our analysis. For the same reason, we also take the first state to be the lightest one, m χ 1 < m χ 2,3 .
The χ i interact with the SM through FGBs. This also induces the decay of the heavier two states in the DM multiplet, χ 2,3 , to χ 1 . We parametrize the relevant interactions with where P R,L ≡ 1 2 (1 ± γ 5 ) and k, l = 1, . . . , 6. The couplings of χ i to the gauge bosons are where U diagonalizes the m χ mass matrix, U † m χ U =m χ , and W diagonalizes the gaugeboson mass matrix. The explicit form of FGB couplings to exotic and SM quarks, Ĝ u/d L/R kl,m , can be found in Appendix A.2 of Ref. [23].
The partial decay width for χ i → χ j q lqk is, neglecting hadronization effects, where the sum runs over the FGB mass eigenstates m = 1, . . . , 24. Expression (17) is valid in the |∆m ij | m χ i limit, and neglecting the quark masses. The above approximations are valid for all the values of parameter for which the correct relic abundance is obtained and the FCNC, collider and direct DM detection constraints are satisfied.
If the mass splitting is less than the pion mass the decay χ i → χ j q lqk is kinematically not allowed. The heavier χ i states then decay radiatively to χ i → χ j γγ. For our purposes an order of magnitude estimate of the decay width suffices. The naive dimensional analysis estimate gives where Q u = 2/3 and Q d = −1/3 are the electromagnetic charges of up and down quarks.
The sum over m runs over the FGB mass eigenstates, while the sum over f is over the SM quarks and exotic states, of mass m f (for up, down and strange quarks this needs to be replaced with Λ QCD ).
The main difference with the fermionic flavored DM from the previous subsection is that the scalar DM interacts with the visible sector via two different types of interactions. The first is its couplings to the FGBs, which is similar to the case of the fermionic DM. The second is a direct coupling to the Higgs For a thermal relic the DM annihilations proceed predominantly through the Higgs portal.
The interactions via FGBs are subdominant except if m φ m a A /2 for some A a . Unless this is the case, the fact that the DM carries a flavor quantum number is exhibited only through the multiplicity of the states.
After electroweak symmetry breaking, the DM-Higgs interactions are given by and the DM mass term m 2 0 φ † φ is shifted by the Higgs condensate to give the DM mass of The vevs of the flavons, Y u and Y d , split the DM multiplet at tree level through The spectrum is also split by radiative corrections due to FGBs. These are quadratically divergent and proportional to the square of the FGB mass. In principle, it is possible to fine tune the tree-level and loop contributions to give almost degenerate DM flavor multiplet.
However, given the hierarchical FGB masses, it is more likely that the DM multiplet is split completely, and only the lightest state is relevant for DM phenomenology. Depending on the signs of κ i in Eq. (23) the lightest φ component can be either top-quark or up-quark flavored.
We choose the latter option in the numerics for easier comparison with the fermionic case.

IV. DARK MATTER AND NEW PHYSICS PHENOMENOLOGY
We turn next to the phenomenology of the flavored DM models. We perform a scan over the parameters of the models and show that the lowest DM states, both for the fermionic DM, χ, and the scalar DM, φ, can be thermal relics. To make the scan numerically tractable we rely on several approximations in calculating the relic density, which we explain below.
We also discuss the predictions for direct DM detection, and the constraints from FCNCs and collider searches.

A. Scan results
In the scan we fix λ u = 1 and vary λ d ∈ [1/(4π), 1]. The range is chosen with the expectation that one will be able to accommodate both the SM top and bottom quark Yukawas as well as satisfy electroweak precision constraints and direct t and b searches [9]. In addition To a good approximation, the variation of M u effectively compensates the fact that we do not vary λ u as seen from Eq. (9). We have verified that further extending these parameter ranges does not extend the viable DM-model parameter space. For instance, the upper All the points shown in Fig. 3 give the correct relic DM abundance, Ω DM . Different colors denote consecutively applied constraints. The grey points fail the perturbativity require- In the remainder of this section we discuss how the various constraints on the DM model were obtained.

B. Thermal relic
For the calculation of the DM relic density we follow Refs. [25,26]. To speed-up the numerical scan we work in the non-relativistic limit, using the freeze-out approximation.

Fermionic dark matter
In the fermionic DM model the DM annihilation to quarks is dominated by s-channel exchange of the lightest FGB, A 24 , see Fig. 6 (left panel). The χ iχi →ū j u j annihilation cross section is given by In the above expression we have neglected the quark masses for simplicity, with the full expression given in Eq. (B9). The rates for A 24 → χ iχi ,d j d j are obtained by trivial coupling replacements and by correcting the color multiplicity factors.
The correct relic abundance requires resonant annihilation, m χ m A 24 /2, see Fig. 3 (upper panels). This implies an upper bound on the DM mass through the following argument. The thermally averaged DM annihilation cross section scales in the narrow width approximation as Here, we used the approximate scaling for the FGB masses and decay widths, In the limit where only χ 1 contributes to the DM relic abundance we find, using the scan, an upper bound m χ 1 10 TeV. The approximation is valid if χ 2,3 states decay well before χ 1 freezes out (i.e. τ χ 2,3 10 −11 s for m χ ∼ 1 TeV). For purely radiative DM mass splitting this is never the case (c.f. Fig. 5). Instead, if χ 2,3 decay after decoupling, they will also contribute to the final DM relic abundance and one needs to sum all three contributions. In this case, the constraints on the mass spectrum become much more severe. In particular, in order for all χ components to annihilate efficiently their masses need to be within a few decay widths away from the lightest FGB (LFGB) mass. This in turn implies that the (radiative) DM mass splitting has to be of the order of the LFGB width. Since the splitting increases with g U we expect these effects to decrease the effective thermal DM annihilation cross section much before the theory becomes non-perturbative. Indeed we find, using the scan, an upper bound m χ 1 5 TeV.  In the calculation of the thermal relic abundance we include the following annihilation channels: where c W = 1, c Z = 1/2 and The thermally averaged cross sections and relic abundances are computed following the prescription described in Appendix B. The results of the scan are given in Fig. 3  Note that the role of the Higgs portal may be played by other light scalars. In Ref. [20] the flavon field of the Abelian horizontal symmetry was used to enhance the DM annihilation cross section. If the flavons are light, they can also modify the phenomenology of

C. Cosmology
The heavier flavored DM states, both for the fermionic DM, χ 2,3 , and scalar DM, φ 2,3 , are unstable. They decay through the χ i → χ jq q transition when the mass splitting is larger than the pion mass, and through the χ i → χ j γγ otherwise, cf. Eqs. (17), (18). The SM particles in the final state of these decays can have various observable effects in cosmology and astrophysics.
The two relevant sets of parameters are the lifetimes of the two heavy states, τ χ 2,3 , and the related mass splittings of the DM multiplet (with respect to the lightest state), ∆m 31 , ∆m 21 .
The lifetimes determine at which cosmological epoch the heavy states decay. The mass splittings control the released combined electromagnetic and hadronic energy, E vis ∆m 21,31 .
They also determine the relic abundances of the heavy states. Generically, near the degenerate limit each state contributes roughly a third of the total DM relic abundance, Ω DM .
Close to the resonant condition m χ m A 24 /2 the χ 1,2,3 relic abundances may differ from Ω DM /3, depending on the common DM mass and relative mass splittings.
For the scalar DM the mass splitting is expected to be large. The φ 2,3 therefore decay before primordial nucleosynthesis. The decays yield negligible entropy release due to the small φ 2,3 abundances. Such scenarios are basically unconstrained by current cosmological observations. The same is true for the fermionic DM if additional spurions split the DM multiplet at tree level.
If the fermionic DM multiplet is split solely by radiative corrections, the χ 2,3 and χ 1 are generically much more degenerate, cf. In the remainder of this section we consider in more detail the region τ χ 2,3 ∼ 10 −1 s−10 12 s, where the big bang nucleosynthesis (BBN) provides the most stringent constraints [29,33].
The injection of energetic photons or hadrons from χ 2,3 decays during or after BBN adds an additional non-thermal component to the plasma that can modify the abundances of the light elements [34][35][36][37][38]. The bounds differ depending on whether the decays result in hadronic or electromagnetic showers in the plasma. The most stringent bounds are for a relic that produces mostly hadronic showers. This is because the electromagnetically interacting particles such as photons and electrons thermalize very quickly by interacting with the tail of the CMB distribution until the universe is 10 6 s old. In our case, the decays χ 2,3 → χ 1 qq are always kinematically allowed for τ χ 2,3 < 10 12 s. The χ 2,3 decays thus predominantly produce a small number of hadronic jets with a combined released hadronic energy E had E vis .
There are three distinct ranges of lifetimes [39]. For 0.1 s τ χ 2,3 100 s the dominant effect is the inter-conversion between protons and neutrons, which overproduces the 4 He abundance. For longer lifetimes, 100 s τ χ 2,3 10 7 s, hadro-dissociation is the most efficient process and the bounds come from the non-thermal production of Li and D. At late times, 10 7 s τ χ 2,3 10 12 s, photo-dissociation caused by direct electromagnetic showers or by electromagnetic showers from daughter hadrons can lead to overproduction of 3 He.
We impose the 4 He, D and 3 He constraints 1 using the results in Ref. [39]. The visible energy release in the decays is E vis ∼ ∆m 21,31 . For 100 GeV < ∆m 21,31 < 10 TeV the constraints derived from the three relic mass benchmarks in Ref. [39] are well approximated by a power-law scaling with E −η i vis . The exponents for the three constraints are η4 He ≈ 1/3, η D ≈ 1/2 and η3 He ≈ 1. For inter-conversion and hadro-dissociation the power-law scaling is expected to break down at energies below O(10) GeV due to the presence of hadronic thresholds [39]. We thus do not extrapolate the fit results for 4 He and D for ∆m 21,31 below 10 GeV. We assume that the photo-dissociation effects retain approximate power law behavior for E vis large compared to the binding energies of the light nuclei, which is of the order of few tens of MeV. In our model for τ χ 2,3 < 10 2(12) s, the mass splitting, ∆m 21,31 , is always above 10(0.1) GeV. Our approximations are thus always valid for ranges of lifetimes for which the 3 He constraints are the most stringent. For the deuterium bound, on the other hand, the power-law scaling is expected to fail for part of the parameter space where the bound is the most stringent, since ∆m 21,31 can be as low as a few GeV. We have checked using the power-law derived bound that these regions are excluded by several orders of magnitude. This gives us confidence to conclude that they remain excluded even with a more faithful treatment of hadro-dissociation effects. In The Wilson coefficients f p and f n are the couplings to protons and neutrons, respectively, G u,d V are the vectorial couplings of FGBs to quarks, defined in Eq. (25). The spin-independent DM-nucleon cross section as measured by the LUX experiment [42] is where µ χn is the reduced mass of the (χ, n) system. The Xenon atomic and mass numbers are denoted by Z and A, respectively. We thus have Z = 54, while A varies between 128 and 134. With the above relations we calculate the DM-nucleon cross section and compare it with the current best limits reported by the LUX experiment [42]. The results of the scan are shown in Fig. 8. Most of the points lie well below the present LUX bound. This is a consequence of the fact that the relic abundance is given by the s-channel resonant annihilation, while the direct detection scattering is due to t-channel exchanges of FGBs and thus not resonantly enhanced.
For scalar flavored DM the dominant scattering is through t-channel Higgs-boson exchange. This leads to the spin-independent scattering on nucleon N = n, p, [43,44] σ SI χN =  The LUX bound is the brown shaded region. The color coding for the points is as in Fig. 3.
The Higgs-nucleon coupling is where the sum runs over the light quarks, q = u, d, s. f analysis of the πN -scattering data gives σ πN = 59 (7) MeV [46]. This is in agreement with a BχPT fit to world lattice N f = 2 + 1 QCD data, which gives σ πN = 52(3)(8) MeV [47]. Including both ∆(1232) and finite-spacing parametrization in the fit shifts the central value to σ πN = 44 MeV. To be conservative we use σ πN = (50 ± 15) MeV that leads to by using expressions in Refs. [48,49]. This results in f N,h = (29.7 ± 1.3) · 10 −2 where we averaged over Higgs couplings to proton and neutron (the difference is an order of magnitude smaller than the quoted error). The resulting bound from the LUX experiment, assuming correct relic abundance, is shown in Fig. 7 and constrains m φ 1 150 GeV.
Finally, we discuss the constraints from indirect DM searches. The flavored DM annihilates to quarks so that the most constraining indirect DM searches are due to the photon and antiproton cosmic-ray fluxes. The constraints from the antiproton flux are quite depen- thermal DM annihilating to quarks [52] (there are also slightly less stringent constraints from γ-ray emissions from the Large Magellanic Cloud [53], and from isotropic γ-ray background [54]).

E. Searches at the LHC
The partonic differential cross section is given by where, in the partonic center of mass frame, M is the total energy, β f is the velocity of the final-state quarks, z = cos θ * is the cosine of the polar angle of the outgoing quark w.r.t. the direction of the incoming quark, and the couplingsĜ V ,Ĝ A of FGBs to quarks were defined in Eq. (25). We have only included the s-channel contribution that dominates on the FGB resonance peak. Terms odd in z were dropped in the differential cross section since they vanish after integration. The predicted dijet cross sections at the LHC with √ s = 8 TeV are shown in Fig. 9, where the 95% CL exclusion from Ref. [55] is denoted with a solid orange line. This mostly excludes the points where the lightest FGB has large couplings to the quarks. Such points are in fact already mostly excluded either by the perturbativity requirement or from flavor constraints.
The quark partners, u i , d i , have an inverted mass hierarchy w.r.t. the SM quarks so that in most of our scan points the t is the lightest state. The bound on the t mass depends on the t → bW , tZ, and th branching ratios. The respective partial decay widths are given by where x i = m i /m t and s i , c i are the sines and cosines of the mixing angles between the SM and exotic quarks, while V is a unitary matrix describing the misalignment of the Y u and Y d vevs. Their definitions can be found in Ref. [23], where also the relevant Feynman rules are given. (We present the relevant Higgs Feynman rules in App. C, correcting an obvious typographical error of Ref. [23].) In Eq. (38) we took the limit x b → 0 that is justified since m t m b . We use the above expressions for the t → bW, tZ, th branching ratios to obtain the 95% confidence-level bound on the t mass by interpolating between quoted observed-limits table in [56]. In our analysis we include the mass differences in the neutral K 0 , B 0 s , and B 0 d sectors, ∆m K , ∆m B d , and ∆m Bs , respectively. We also include the indirect CP violation in the kaon sector, ε K , and the mixing-induced CP asymmetries S ψKs and S ψφ . The corresponding mixing amplitude where i, j are the flavor indices. The RG evolution is implemented following Ref. [59] (for further details and the dependence of the numerical relevance with the scale see also Ref. [23]). For the non-perturbative inputs, the decay constants and the bag parameters, we use the current lattice averages [60].
The mass difference in the neutral kaon sector, ∆m K , and the CP-violating parameter ε K are given by with ϕ ε = (43.51 ± 0.05) • and κ ε = 0.923 ± 0.006, which includes long-distance effects in both Im M K 0 12 [61] and in the decay, i.e. Im Γ K 0 12 [62]. Our SM expectation for ε K incorporates the known Next-to-Next-to-Leading-Order (NNLO) QCD corrections due to the charm [63] and charm-top [64] contributions. The mass differences in the neutral B sectors are given by The CP violation in the neutral B sector is probed by the time-dependent asymmetries in the decays B 0 d → ψK S and B 0 s → ψφ that define the observables S ψKs = sin(2β + 2φ B d ) and The rate of B → X s γ is also modified by the presence of exotic up-type quarks. These can only enhance the B → X s γ rate with respect to the SM expectations [9]. The contributions of FGBs are loop-suppressed. Even though they may be enhanced by m b they are negligible in models with a seesaw-like mass generation for quarks, like the model we consider [65]. The SM prediction for the rate in our analysis includes the known NNLO corrections [66][67][68].
In our numerical scan we mark parameter space points to have passed the flavor con- a six-point star (4).
"Benchmark 1" is an example of fermionic flavored DM, where the DM multiplet is light, with mass below 1 TeV. The mass of the lightest state in the DM multiplet is m χ 1 520(540) GeV if it lies just below (above) the LFGB resonance. If the mass splitting between χ 1 , χ 2 and χ 3 is solely due to radiative corrections, χ 2 and χ 3 are almost mass degenerate with masses roughly 10 GeV above m χ 1 , and χ 3 is about 100 MeV heavier than χ 2 . The lightest quark partner is t with a mass m t 1.3 TeV. The lightest FGB has a mass m A 24 1.1 TeV.
All the remaining NP states are above 7 TeV. This benchmark point demonstrates that even parameter regions with low lying FGBs can be consistent with both the resonance searches and the FCNC bounds. The most robust constraints in this parameter region come from cosmology (in case of radiatively-split DM masses) and dijet-resonance searches (see Figs. 5 and 9). Note in particular that for the completely (mass) decoupled fermionic DM scenario, in which cosmology bounds are absent, all experimental constraints can be satisfied even for DM (and LFGB) masses below few 100 GeV.
The bottom left panel in Fig. 10 shows the predicted relic abundance for this benchmark, if only the DM mass is varied, which also modifies the splitting within the DM triplet. Relic abundance consistent with observations is obtained for a mass of DM close to half of the mass of the lightest FGB, in which case the annihilation cross section is resonantly enhanced. To saturate the observed DM relic density, two solutions for m χ 1 are obtained, with m χ 1 either above or below the resonant peak. We see that for radiatively split DM masses, where all χ i components contribute to the DM relic density, m χ i need to lie within O(3%) of the resonant peak for the annihilation to be strong enough. For completely decoupled DM multiplet the resonant condition is relaxed and needs to be satisfied to O(5%).
In Fig. 10 (upper panel) we show the spectrum for the lower mass solution and radiative DM multiplet splitting. We see that the quark partners of the lighter generations are heavier than the partners of the third generation quarks. Similarly, the FGBs that couple more strongly to the first two generations are typically heavier than the ones that couple preferably to the third generation. The couplings of the lightest FGB to the light quarks have the form where the relative corrections to this relation are below the percent level. This means that the lightest FGB couples to the light quarks vectorially, case of only radiative mass splitting). The lightest exotic quark is the top partner with mass m t 1.7 TeV, while the mass of the lightest FGB is m A 24 9.2 TeV.
For such high DM masses it is barely possible to obtain the correct relic abundance (see the lower left panel in Fig. 11). Therefore, the DM mass is finely tuned to be exactly on the resonant peak (see the lower left panel in Fig. 11), so m χ i m A 24 /2. Because of the high masses of the NP states the direct searches (direct DM detection, t searches and dijet resonance searches) as well as the indirect flavor constraints are easily avoided, although K −K mixing does receive non-negligible contributions.
Also in this case, the couplings of the lightest FGB to the light quarks have the form  To achieve this we adopted the model of Ref. [9] in which the inverse proportionality is achieved by introducing a set of quark partners also necessary to cancel gauge anomalies.
The same quark partners also mix with the SM quarks and lead to the mass hierarchy of the SM quark masses. We have also considered the possibility that the DM multiplet is split due to an extra source of flavor breaking. Also in this case, the correct relic abundance requires resonant annihilation. The DM mass is in thus still roughly equal to half of the mass of the lightest FGB. However, the heavier DM states decay well before big bang nucleosynthesis so that a much wider range of DM masses is phenomenologically viable. In our scan this includes DM masses as light as 100 GeV (with very small couplings to FGBs) and up to 10 TeV.
A possible signal of the gauged flavor model with fermionic DM at the LHC are mono-jets, where the lightest FGB is produced associated with initial state radiation and decays to χ 1 pairs. The χ 1 s are expected to be non-relativistic in the lightest FGB's rest frame, as their combined mass needs to be close to the FGB mass to fulfill the resonance condition for relic DM abundance. In the event that such a signal would eventually emerge, the corresponding dijet-resonance signal is generically also expected in the model. A final possibility in the case of radiatively split DM mass spectrum is that some of the lightest FGBs decay to slowmoving χ 2,3 . They in turn decay within the detector, leaving (highly) displaced vertices, isolated hits in the calorimeter or in the muon chambers. Unfortunately, in most of the parameter space χ 2,3 are expected to decay well outside the detectors, see In short, we have shown, using an explicit renormalizable model, that it is possible for flavored DM to be a thermal relic. The considered model is not the only choice. One could consider DM in other representations of G SM F . Our analysis can be extended also in other ways: for instance, by enlarging the global symmetry as in Ref. [13] and subsequently gauging it. For instance, with our field content the global group is G SM but other choices could be made. Yet another possibility is to gauge only part of G SM F , for instance a U (2) 3 ⊂ G SM F . Note that for fermionic DM, Z 3 is part of an accidental U (1) χ acting in the dark sector. The U (1) χ can be broken by the dimension-7 operator LHχχχ, but is exact in our renormalizable model. It can in principle be gauged and lead to additional phenomenology. If DM is a scalar, U (1) χ can be broken already at the renormalizable level, leaving only Z 3 exact.
andQ bd 1 ,Q bd 3 , that follow from Q bd 1 , Q bd 3 with P L ↔ P R exchange (the remaining operators can be found in, e.g., Ref. [69]). If the lightest FGB has predominantly left-handed couplings, then the C bd 1 Wilson coefficient is the largest one. If the lightest FGB couples predominantly to the right-handed quarks, thenC bd 1 dominates. For comparable left-and right-handed couplings all four Wilson coefficients, C bd 1,3 ,C bd 1,3 , are important. The analogous discussion applies to the case of FGB contributions to B s mixing obtained with a trivial d → s exchange.
As discussed in Section III the contributions from the gauged flavor model is expandable in terms of the SM Yukawas. The contributions due to FGB exchanges can thus be written where (y d ) ij = diag(y d , y s , y b ), and we set y t = 1 in the second equality. In Eq. (A3) we kept only the two terms relevant for the discussion below. The same expansion applies for δC bs 1 with the replacement d → s in Eq. (A3).
In Fig. 14 (left) we show the ratio δC bd 1 /δC bs 1 , i.e. the NP contribution to the V −A quark current operator due to tree-level FGB exchanges. Note that the ratio δC bd 1 /δC bs 1 can be complex. If c 1 , i.e. the leading MFV (LMFV) term, dominates then δC bd 1 /δC bs 1 (V * td /V * ts ) 2 . This is denoted by a cross in Fig. 14 (left) . The addition of the operators with extra insertions of y d y † d leads to δC bd 1 /δC bs 1 not being equal to (V * td /V * ts ) 2 . We verified that the curve for δC bd 1 /δC bs 1 shown in Fig. 14 Fig. 14 (right), where we plot the cumulative distribution function The function P MFV (n) can be interpreted as the fraction of points that have the ratio |δC bd 1 /δC bs 1 | effectively dominated by operators with up to y n d insertions. That is, the points dominated by the c 1 term contribute to P MFV (0) (and to P MFV (n) with n ≥ 0), while the points dominated by the c 2 term contribute to P MFV (2) (and to P MFV (n) with n ≥ 2). The points with both c 1 and c 2 start contributing to P MFV (n) for n somewhere between 0 and 2, depending on the relative sizes of c 1 and c 2 . Fig. 14 (right) shows that the c 1 term dominates in a subleading (but nonzero) set of points, that about 10% points have both sizeable c 1 and c 2 terms, and that c 2 dominates in about 80% of the points.
The similar analysis can be performed for V + A operators,Q bd 1 andQ bs 1 . We expand the FGB contributions to their respective Wilson coefficients in terms of the SM Yukawas and similarly for δC bs 1 with the d → s replacement. We show in Fig. 15 (left) the ratiõ C bd 1 /C bs 1 for our scan. If thec 1 term dominates, thenC bd 1 /C bs 1 = (m d /m s ) 2 (V * td /V * ts ) 2 , which is denoted by the cross in Fig. 15 (left). The points for which also thec 2 operator (and other operators denoted by ellipses above) is important then lie away from theC bd 1 /C bs 1 = (m d /m s ) 2 (V * td /V * ts ) 2 region. We also define a cumulative functioñ The values forP MFV (n) are shown in Fig. 15 (right). We see that also in this case the points cluster into two groups, with vanishingc 1 term or with bothc 1 andc 2 relevant.
The above analysis demonstrates that the FGB contributions to the Wilson coefficients in the effective weak Hamiltonian can be expanded in terms of the SM Yukawas. This is a hallmark of (general) MFV. In particular, the expansion in terms of m d,s /m b and off-diagonal CKM elements can still be performed and is not ruined by the large ratios of scales present in the problem such as the very disparate FGB masses.

Appendix B: Thermal relic computation
In this appendix we describe the calculation of relic density that was used in the scans in the main part of the paper. Several approximations to the coupled Boltzmann equations were necessary in order to reduce the evaluation time per benchmark and thus allow adequate coverage of the parameter space. We find the approximate solutions to be in agreement with the full solutions at the O(30%) level. The full numerical solution of the Boltzmann equations was obtained with MadDM [27] using a UFO model file [70], which was generated with the FeynRules package [28].
Denoting the DM multiplet by ϕ, where ϕ is either a Dirac fermion or a complex scalar, the most general set of coupled Boltzmann equations reads [25] dn ϕ i dt + 3Hn ϕ i = − j σ(ϕ i ϕ j ↔ XX) v lab n ϕ i n ϕ j − n eq ϕ i n eq ϕ j − j =i σ(ϕ i X ↔ ϕ j X)v lab n ϕ i − n eq ϕ i n eq ϕ j n ϕ j n eq X − j =i σ(ϕ i ϕ j ↔ ϕ k ϕ )v lab n ϕ i n ϕ j − n eq ϕ i n eq ϕ j n eq ϕ k ϕ eq n ϕ k n ϕ ± j =i Γ(ϕ j,i → ϕ i,j X) n ϕ j,i + σ(ϕ j,i X → ϕ i,j ) n ϕ j,i n eq X , where X denotes a generic SM state. For large mass splittings between the ϕ components it is sufficient to consider the lightest ϕ i state in Eq. (B1). The contributions to the DM relic density from the heavy ϕ components are exponentially suppressed by corresponding Boltzmann factors and can be neglected within our precision. In contrast, when the mass splittings are small the full set of coupled Boltzmann equations in Eq. (B1) needs to be considered. Nevertheless, even in this case several approximations are possible for our model, as we explain below.
First of all, the coannihilation of different ϕ i components into SM particles, ϕ i ϕ j → XX (i = j), can be safely neglected in our model. In benchmarks that survive the experimental constraints the off-diagonal couplings of the lightest FGB to ϕ are much smaller than the diagonal ones, see Fig. 16. Secondly, in the calculation of DM relic density we also neglect the flavor-changing DM scattering off the thermal background, ϕ i X → ϕ j X. The ϕ i X → ϕ j X scattering can be important if σ(ϕ i X → ϕ j X)v lab / σ(ϕ i ϕ j → XX) v lab n ϕ j /n eq X ∼ 10 −9 . In this case the off-diagonal couplings of O(10 −4 ) relative to the diagonal ones are in principle large enough to have O(1) effects on the relic density, and neglecting ϕ i X → ϕ j X may not be justified. Therefore, for the benchmarks with (ĝ 24 χ ) 23 /(ĝ 24 χ ) 33 > 3 × 10 −4 and small mass splittings among ϕ, we explicitly verified using MadDM that neglecting ϕ i X → ϕ j X scattering leads to a change in DM relic density smaller than O(30%).
Finally, we are able to neglect the pure DM scattering process in the third line of Eq. (B1) since σ(ϕ i ϕ j ↔ ϕ k ϕ )v lab σ(ϕ i ϕ j ↔ XX) v lab in our model. The largest contribution to this process is from diagonal couplings between the FGB and DM. This process can couple the evolution of the DM species if σ(ϕ i ϕ j ↔ ϕ k ϕ )v lab ∼ σ(ϕ i ϕ j ↔ XX) v lab . The diagonal FGB couplings to the quarks and to the DM of the same generation are approximately equal. By accounting for the color factors and the multiplicity of channels when annihilating into SM fields one concludes that the pure DM scattering is indeed subleading.
Therefore, for almost mass degenerate ϕ i it is sufficient to consider a set of uncoupled Boltzmann equations The DM relic abundance is in this case the sum of relic abundances for each of the three components obtained from the above set of equations (the heavy ϕ components, in our case ϕ 2 and ϕ 3 , decay after their respective freeze-outs and contribute to the ϕ 1 DM relic abundance). In contrast, for large mass splittings the heavy ϕ components are irrelevant for the calculation of the DM relic abundance. This is then obtained from Eq. (B2) by considering only the lightest DM state, in our case ϕ 1 .
We calculate the DM relic abundance from Eq. (B2) using the freeze-out approximation [26], which gives where M P l = 1.22 × 10 19 GeV, g * is the total number of effectively relativistic degrees of freedom at the time of the freeze-out, and The freeze-out temperature (x f = m ϕ 1 /T f ) is obtained by solving x f = ln 0.038 g eff M P l m ϕ 1 σv lab th √ g * x f , where the thermally-averaged cross section is σv lab th = 2x with v lab = 2 (1 + )/(1 + 2 ) and = s/(2m ϕ 1 ) 2 − 1. The freeze-out approximation is accurate to a few percent with respect to the full numerical solution of the Boltzmann equation [25].
The fermionic flavored DM annihilates through the s-channel exchange of FGBs. In this case, the integration over x can be performed analytically and the double integral in Eq. (B4) reduces to a single one that can be efficiently evaluated numerically. In particular, We evaluate the above integral numerically in the parameter scan.
The rate for A 24 → χ iχi ,d j d j is obtained after trivial replacements for masses and couplings (and dividing by the N c color factor for decays to χ iχi ). The total FGB decay rate is obtained after summing over all kinematically allowed decay channels.