Dark Matter interactions in an $S_4 \times Z_5$ flavor symmetry framework

The interactions of dark matter (DM) with the visible sector are often phenomenologically described in the framework of simplified models where the couplings of quarks to the new particles are generally assumed to be universal or have a simple structure motivated by observational benchmarks. They should, however, a priori be treated as free parameters. In this work we discuss one particular realization of the structure of DM couplings based on an $S_4 \times Z_5$ flavor symmetry, which has been shown to account reasonably well for fermion masses and mixing, and compare their effect on observational signals to universal as well as Yukawa-like couplings, which are motivated by minimal flavor violation. We will also comment on how these structures could be constrained in UV complete theories of DM and how DM observables, such as, e.g., relic density and direct detection, can potentially be used as a smoking gun for the underlying flavor symmetries.

The interactions of dark matter (DM) with the visible sector are often phenomenologically described in the framework of simplified models where the couplings of quarks to the new particles are generally assumed to be universal or have a simple structure motivated by observational benchmarks. They should, however, a priori be treated as free parameters. In this work we discuss one particular realization of the structure of DM couplings based on an S 4 × Z 5 flavor symmetry, which has been shown to account reasonably well for fermion masses and mixing, and compare their effect on observational signals to universal as well as Yukawa-like couplings, which are motivated by minimal flavor violation. We will also comment on how these structures could be constrained in UV complete theories of DM and how DM observables, such as, e.g., relic density and direct detection, can potentially be used as a smoking gun for the underlying flavor symmetries.

Introduction
Motivated by a substantial amount of evidence from cosmological observation, there is a continuing effort by experimental collaborations to obtain evidence for the existence of particle dark matter (DM). This effort is supported by a large number of theoretical analyses investigating the potential characteristics of DM in preparation for a future discovery 1 . There are several options for the choice of framework in which DM can be studied. One can, for example, rely on specific Particle Physics scenarios, like Supersymmetry, and directly confront their predictions with experiments. An orthogonal approach, whose popularity is increased in the recent years, consists of considering simplified models, agnostic to the details of underlying BSM theories, containing only the relevant degrees of freedom, namely the DM and a mediator of its interactions with (typically) SM fermions, and interactions accounting for the DM relic density and possible observable signals at experiments. This kind of approach has been extensively applied to WIMP models [2][3][4][5][6][7][8][9][10][11][12][13][14][15]. As pointed out, e.g. in [8,9,[16][17][18], interpreting the experimental outcome in terms of simplified models requires care in view of their theoretical limitations, as for example the lack of gauge invariance. For this reason the theoretical community is currently working on more solid refinements of simplified models, see e.g. [5,9,13,15,[19][20][21][22][23][24][25] for discussions. By refinements, we mean models in which the coupling of a SM singlet DM candidate is realized in a renormalizable and possibly gauge invariant way but still maintaining a low number of free model parameters.
Along this line of reasoning this manuscript aims to investigate the flavor structure of simplified models. For the latter, one typically relies on the simplifying assumption of flavor-diagonal couplings of the mediator. While this assumption is reasonable and elegant, UV complete scenarios might motivate different choices for the assignations of the couplings of the mediator with SM fermions, see e.g. [26][27][28]. Non trivial flavor structures in the fermion mass matrices can be transferred to the couplings of the SM quarks with the mediator fields, and can be obtained by means of additional symmetries which add to the gauge group of the SM. In recent years, they have been quite popular thanks to their ability to explain with good accuracy the observed pattern of neutrino masses and mixing, generally at the prize of introducing heavy scalar degrees of freedom to mediate the symmetry breaking. Among the various realizations, discrete non-abelian groups G f [29,30] have gained a lot of interest since they have non-trivial multidimensional representations useful to inter-relate different families; in particular, permutation groups as A 4 , S 4 ... were studied in detail as neutrino mixing matrices of tri-bimaximal-mixing [31] and bimaximal-mixing [32] forms naturally emerged from particle assignment to the irreducible tridimensional representations. Quite often such discrete symmetries have been supplemented by U (1) groups a la Froggatt-Nielsen [33] or their discrete counterparts, Z N groups, (N being the order of the group formed from the p-th roots of unity) in order to eliminate unwanted operators preventing the correct description of fermion mass hierarchies. Although not exhaustively discussed in the literature, the compound G f × Z N has been also applied to quarks (see for instance [34][35][36][37][38][39][40]), showing that a leading order (LO) diagonal V CKM can be easily reproduced and higher-order corrections are needed to accommodate the off-diagonal entries. In the present work we adopt a flavor model based on an underlying S 4 × Z 5 flavor symmetry [41] which goes along the lines proposed in [34], making explicit computation of all next-to-leading order (NLO) terms necessary to describe quark masses and mixing.
Our paper is organized as follows: In section 2 we introduce two low energy models, featuring scalar and fermionic DM coupled with a scalar mediator in turn coupled with SM quarks. The flavor structure of the latter coupling will be inspired by the S 4 ×Z 5 symmetry. In section 3 the main phenomenological constraints applied to our study will be illustrated and then transferred to the simplified model in section 4. In section 5, before the conclusions, we will instead consider the phenomenology of a concrete realization of the model, illustrated in detail in appendix A, in which the results are also dependent on a New Physics scale Λ associated to the flavor symmetry. Besides the conventional freeze-out, in this last scenario we will consider also the case of DM production through the so-called freeze-in mechanism. The latter allows to achieve the correct relic density for very high values of Λ, namely 10 11÷12 GeV. In appendix B we finally give some general details about the S 4 group theory.

Low energy simplified model
For what concerns Dark Matter and collider phenomenology, the framework under study can be represented as a portal scenario in which an electrically neutral scalar field φ couples with pairs of SM quarks and pairs of DM particles. For the latter, we will focus, for definiteness, on the cases of scalar S and (dirac 2 ) fermionic χ DM.
The corresponding lagrangians are given by: and for scalar and fermionic DM, respectively. In both lagrangians h 1 = h 2 represent 3 × 3 matrices in the quark flavor space. Since we are considering the mediator φ coupled with all the six SM quarks, there are two copies of h 1 and h 2 associated, respectively, with u-type and d-type quarks. The structure of the matrix elements is specified by the S 4 × Z 5 symmetry [41] (discussed in details in Appendix A) which, at the leading order, are: with ε = 0.05 (and O(1) coefficients simply fixed to unity). The same definition applies also for the matrix acting on the u-type quarks. As a further assumption, we will take, throughout our study, µ 1 = 0, for both fermionic and scalar DM. Two notable differences appear with respect to the simplified models usually considered in the literature. First of all the couplings h 1 , h 2 should be intended as 3 × 3 matrices whose structure is determined by a specific UV model (see next subsection for more detail on this realization) and should comply with specific observational constraints from processes related to flavor physics. A further important feature is represented by the fact that the lagrangians are explicitly CP violating. This is required to properly account for the flavor structure of the SM, and has also relevant phenomenological implications for DM.
The coupling of the mediator with SM quarks originates, at the loop level, an effective coupling with gluons and photons, relevant for collider phenomenology. The latter can be described through the following lagrangian: where [42]:h Here A S q and A P q are loop functions given by: where and withŝ being the gluon-gluon center of mass energy.

Relic Density
Throughout this paper we will mostly consider the standard thermal freeze-out paradigm for the determination of the DM relic density. According to it, the latter is determined by the thermally averaged DM pair annihilation cross-section σv . In the scenario under consideration DM annihilates into SM quark pairs, through s-channel exchange of the mediator φ, and into φφ, through t-channel exchange of the DM (in the case of scalar DM a contact interaction vertex is present as well), if kinematically allowed 3 . For completeness, we have also included in our analysis the annihilation channels into gg and γγ originated by the Lagrangian in 4. Useful analytical expressions can be found e.g. in [15,23,[45][46][47]. Our results will be, however, obtained through precise numerical computations as performed by the package micrOMEGAs [48,49].

Direct Detection
The flux of WIMP DM through a detector volume on Earth can lead to scattering events in a suitable target. This in turn can be detected as recoil energy. This kind of phenomena is at the base of the so-called DM Direct Detection (DD). For the models under consideration, the DM scattering processes are described, at the microscopic level, by interactions of DM pairs with quark (and gluon) pairs through t-channel exchange of the mediator φ. These microscopic interactions lead mostly to Spin-Independent interactions between the DM and the nucleons whose corresponding cross-section is given by: where µ DM,N is the DM-nucleon reduced mass while c DM,N is the effective DM-nucleon coupling. In the case of scalar DM the latter is given by: where f T q are structure functions describing the contribution of the light quarks to the nucleon mass. For these we have adopted the following numerical values [50][51][52] : 3 Notice that the relic density phenomenology can be altered in more realistic setups in which additional BSM states are present. See e.g. [43,44].
In the case of fermionic DM the effective coefficient is a combination of a tree-level and loop-induced contributions [53][54][55][56]: The quantities C box i,q and C box G(G),S(P S) are loop functions which depend on the model parameters, i.e. m χ,φ , λ 1,2 and h 1,2 . Since these expressions are rather complicated, we do not report them here explicitly. The interested reader can find them in [54,56].

Indirect Detection
All the scenarios considered in present work are characterized by DM annihilation cross-sections with sizable s-wave component. Consequently, we can have efficient residual annihilation processes at present times which can be probed through DM Indirect Detection (ID) search strategies. However, while this is a generic feature of the case of scalar DM, for fermionic DM the prospects for ID rely only on the pseudoscalar coupling λ 2 . This is because the scalar coupling χχφ is responsible of a velocity dependent (p-wave) contribution to the annihilation cross-section, whose value at present times is hence suppressed with respect to the one at thermal freeze-out. The pseudoscalar coupling χγ 5 χφ leads, instead, to an s-wave contribution to the DM annihilation cross-section. DM annihilations into SM quarks lead to a γ-ray signal with continous energy spectrum, originating in the quark hadronization process (e.g. decay b → π 0 γ). This kind of signal can be probed, for the range of DM masses considered in our work, by the FERMI-LAT experiment [57,58]. Current constraints exclude cross-sections of the order of the thermally favoured value for DM masses 150 GeV. In addition to the just illustrated continuous γ-ray signal, mono-energetic (lines) γ-rays can be produced as well by DM DM → γγ processes, made possible by the loop induced coupling of the mediator φ with a photon pair. We have adopted the constraints presented in [59]. As shown in figs. 1-2, the latter are nevertheless subdominant.

Collider searches
Given the lagrangians (1,2,4), potential collider signals originate from the resonant production of the mediator φ through gluon fusion or quark-quark fusion. The latter process has, however, negligible impact for the assumed flavor structure and the assignation ε = 0.05. Subsequent decays into visible states lead mostly to a djiet or diphoton signal while, in the case of sizable decay branching fraction into DM, a monojet plus missing energy signal would be originated. The model under consideration cannot, however, be probed through the collider searches mentioned before. As shown, e.g. in [60][61][62][63], current results can probe production of mediators with couplings bigger than the order of the SM Yukawa couplings; the top Yukawa coupling y t is, in particular, crucial to have a sizable production vertex. In the case of S 4 × Z 5 couplings, with ε = 0.05, top and bottom quark loops equally contribute to the production vertex. Given that, nevertheless, ε y t , the production cross-section of the φ state is sensitively more suppressed with respect to the Yukawa simplified model.

Results for the simplified model
We have now all the main ingredients to characterize DM phenomenology within the model specified by the couplings in eq. (3). As a first illustration of our results, we have shown in fig. 1 and fig. 2 the combination of DM constraints in the bidimensional (mass of the mediator, coupling)-planes, for three different assignations of the mediator mass, namely 10, 100 GeV and 1 TeV. In all the plots the parameter spaces corresponding to the correct DM relic density is represented by black isocontours, the regions excluded by DM DD are marked in blue while the regions excluded by indirect searches of DM annihilations in γ-ray continuum (lines) have been marked in red (orange).
In the case of scalar DM we have considered the (m S , g 1 ) bidimensional plane, with the other coupling g 2 set to zero.
This kind of scenario is constrained by both Direct and Indirect detection with the former typically giving the most stringent constraints. As evident, the lightest benchmark with m S = 10 GeV is completely ruled out by the experimental constraints. The latter exclude also most of the parameter space for m S = 100, GeV, ad exception of the pole region, i.e. m S m φ 2 . In the case m S = 1 TeV, the DM relic density is mostly accounted for by the SS → φφ annihilation process so that the corresponding isocontour substantially coincides with an horizontal line. Bounds from Direct Detection are evaded for m φ 500 GeV.
In the case of fermionic DM, shown in fig. 2, we considered the (m φ , λ 1 ), with λ 2 = 0, and (m φ , λ 2 ), with λ 1 = 0. All the considered DM mass assignations are ruled out in the case only the λ 1 coupling is on. This is so because, if only λ 1 = 0, the DM annihilation cross-section is velocity suppressed. Consequently, different from the scalar DM case, larger values are required to match the thermally favored value. This causes, in turn, much more stringent constraints from DM DD. More interesting is the case in which only the coupling λ 2 is different from zero. Indeed it contributes to the DM scattering cross-section only at the loop level and, as evidenced by the right column of fig.2, constrains the parameter space to a negligible extent. The coupling λ 2 is instead sensitive to constraints from  Indirect Detection which, however, exclude only the lightest assignation of the mass of the DM. It can be easily argued that, given the presence of the two couplings λ 1,2 , with different properties, there is a broader available parameter space, with respect to the case of scalar DM.
Instead of assuming some model parameters fixed to constant values, we made our results more systematic by performing a parameter scan over the following ranges: for scalar and fermionic DM, respectively. As already mentioned we have set the tri-scalar coupling µ 1 to zero. The parameter assignations passing all the constraints illustrated in the previous subsections are shown, as blue points, in fig. 3 for scalar DM, and fig. 4 for fermionic DM. For comparison, the same plots show the results of analogous scans conducted for simplified models with Yukawa-like (green point) and flavor-universal (red points). For Yukawa-like simplified model we intend the case in which [64]: for up-type and d-type quarks respectively. y f =u,d,s,c,b,t are the SM yukawa couplings at the EW scale, whose numerical values, have been explicitly reported, for convenience, in eq.18.
By quark-universal model we intend, instead, the assignation: with I being the identity matrix in the flavour space. c φ is common factor which has been varied between 10 −6 and 1 (the reason for this choice of the range of the scan will be clarified below). Focusing, for the moment, on the S 4 × Z 5 model, we see from fig. 3 that the scalar DM scenario is rather constrained. In agreement with the findings of fig. 1, viable model points are found only for m S > m φ or around the m S m φ /2 "pole". Furthermore, only model points with m φ 100 GeV can comply with all the phenomenological constraints. In the case of fermionic DM, on the contrary, we notice viable solutions also for m χ < m φ /2. This is due to the fact that the λ 1 coupling influences mostly DD while it has negligible impact on ID; on the contrary the λ 2 coupling impacts mostly ID. It is then possible to evade experimental constraints by a suitable combination of these two couplings.
Let's now compare the results for the S 4 × Z 5 model with the flavour universal and yukawa cases. In the former case we see that, for both scalar and fermionic DM, there are viable regions of parameters space only for m χ,S > m φ . This is due to the fact that, in the case of universal couplings, the contributions to the DD cross-section from quarks of the first generations are dramatically enhanced, since they are proportional to 1/m q . This does not occur in the other two models since they feature automatically suppressed couplings of the mediator with the first two quark generations. Experimental constraints can be passed only for c φ 0.01. This implies, in turn, that annihilations into SM fermions are too suppressed to ensure the correct DM relic density, which can be achieved only for m χ,S > m φ , when the annihilations into φφ are kinematically accessible 4 . This kind of scenario is often dubbed secluded regime. More interesting is the comparison between the yukawa and the S 4 × Z 5 models. In the case of scalar DM we clearly see that the yukawa model have a larger allowed parameter space. There are, in particular, viable solutions for m S < m φ far from the m φ /2 resonance. This is due to the opening, at high DM masses, of the annihilation channel into tt final state. The same does not occur in the case of the S 4 × Z 5 model since ε y t . Less trivial is, instead, the comparison in the case of fermionic DM. The viable parameter space for the Yukawa simplified model extends at slightly lower mediator masses with respect to the S 4 × Z 5 case. The former model appears to be also slightly more favourable for m χ > m t . This can be explained with the fact that in the S 4 × Z 5 model the mediator Φ has slightly stronger couplings with the first generation quarks as well as the bottom quarks, which then face stronger experimental constraints. At the same time, the stronger coupling with the bottom, with respect to the yukawa case, opens a viable region or parameter space for 35 m χ 120 GeV and m φ 300 GeV.

Towards a more realistic scenario
Up to now we have considered the phenomenology of a simplified model in which the couplings between the mediator φ and the SM quarks are inspired by the structure of a specific flavor group. A more realistic realization can be obtained, in a similar spirit as [63], by extending, through two singlet fields, namely the DM candidate and the mediator φ, an S 4 × Z 5 invariant realization of the SM. The details of the model will be left to appendix A. The most relevant impact on DM phenomenology will be due to the dependence of the couplings h 1,2 to the NP scale Λ associated to the breaking of the flavor symmetry: where v is the vev of the SM Higgs.

Freeze-out regime
Assuming again the thermal freeze-out as generation mechanism for the DM abundance, we can repeat the analysis performed in the previous section. There is, however, an additional parameter influencing DM phenomenology, namely the NP scale Λ at which the flavor symmetry is broken. We have then repeated the parameter scan considering different assignations of the latter parameter, namely Λ = 3, 10, 50 TeV. As shown in fig. 5, our findings in the (m χ,S , m φ )-plane have been compared with the results obtained in the previous section for the simplified model of eq.(3) (which would correspond to the case Λ = 246 GeV). As evident, the viable parameter space, for m S,χ < m φ progressively reduces as the scale Λ increases. This is due to the v Λ suppression factor in the coupling between the mediator and the SM quarks. Because of this suppression, the DM annihilation cross-section into SM quarks cannot match the thermally favored values unless one relies on the resonant enhancement occurring for m S,χ ∼ m φ 2 . For Λ > 50 TeV, viable relic density is obtained only in the secluded regime, for m S,χ > m φ . In such a case, indeed, the relic density is determined essentially by the DM annihilation process into φ pairs, whose cross-section is independent from Λ. Consequently, DM evades most experimental searches ad exception of possible indirect signals from the DM DM → φφ → 4f process, mostly in the case of scalar DM.
In the secluded regime DM observables are not affected by specific value (and flavour structure) of the couplings of the mediator φ with the SM states. The only requirement is that these couplings are not too suppressed so that the DM could exist in thermal equilibrium in the Early Universe and then apply Standard Thermal freeze-out computations. As will be seen in the next subsection, this requirement can be used to set constraints on the scale Λ.

Freeze-in regime
As pointed in the previous subsection, the computation of the DM relic density, based on the thermal freeze-out, is based on the hypothesis that the DM was in thermal equilibrium, at least at temperatures higher than its mass. This might not be the case, however, if the interactions of the DM with the SM states are too suppressed. A rule of thumb to assess whether the DM particle was in thermal equilibrium in the Early Universe consists into comparing the DM annihilation rate Γ ann = n eq DM σv , with n eq for scalar and fermionic DM, respectively. Since the DM can also annihilate into mediator pairs, if the process is kinematically allowed, it can be maintained into thermal equilibrium as long as the mediator is. Having assumed µ 1 = 0, the latter condition can be checked by comparing H with the rate associated to the φ ↔ qq process. From this we can infer the following condition on the scale Λ: Even if the DM never achieved thermal equilibrium in the Early Universe, it can be efficiently produced, through the freeze-in mechanism [65], from qq → χχ(SS) annihilation processes. In such a case the DM relic density can be obtained by solving the following Boltzmann's equation, tracking the time evolution of the DM number density n DM : where X is a SM state in thermal equilibrium (the sum over all the possible annihilation process is implicitly assumed). The right hand side of the equation is formally written as: where σ(s) is the annihilation cross-section as a function of the center of mass energy s while the functionK 1 accounts for the fact that in the freeze-in regime one has to adopt on Fermi-Dirac distribution functions, including a chemical potential, encoded in the parameter η i = η s exp (µ i /T ) with η s = −1(+1) for fermions (boson), for the SM fermions rather than relying on the Maxwell-Boltzmann distribution, as done in the case of WIMP production mechanism. Finally x X = m X /T while g X represents the internal degrees of freedom of the SM state X. By doing the customary change of variables n DM → Y DM = n DM /s, d dt → −HT d dT we can write the solution of the Bolzmann's equation as: where the integral is computed between the present time temperature T 0 and the reheating temperature T R . We remind that the reheating temperature is the temperature at which the radiation dominated epoch in standard cosmology begins. The latter has been assumed to be below the scale of breaking of the flavor symmetry, so that lagrangians 1-2 can be adopted for the computations. g * ,ρ and g * ,s represent, respectively, the number of effective relativistic degrees of freedom contributing to the energy density and the entropy density. Contrary to the case of freeze-out, we see that in the case of freeze-in production the DM relic density is proportional to the strength of the DM interactions with the SM states. This integral, and the related DM relic density have been computed through the package mi-crOMEGAs 5 [66]. We will nevertheless provide below some semi-analytical estimates in order to improve the understanding of the results. Fig 6 shows contours of the correct DM relic density in the (m S , Λ) (left panel) and (m χ , Λ) (right panel) planes, for the three assignations m φ = 10, 100, 1000 GeV. In both scalar and fermionic DM cases, the relic density contours evidence two distinct trends. For m χ,S m φ , the DM relic density increases with the DM mass, so that a comparable increase of the value of the NP scale Λ is needed. In the opposite case, m χ,S m φ , the DM relic density appears to be independent from both the DM and the mediator masses, being just set by the value of Λ. These two different regimes can be explained as follows.
Let us start with the case m S,χ m φ . In such a case, for both fermionic and scalar DM, the cross section σ(s) can be approximated as σ(s) κ s v 2 Λ 2 , where κ is a parameter containing the couplings and numerical factors. We can at this point operate the following change of variables: s → z = √ s/T and T → x = m DM /T . In such a way eq. 25 can be rewritten as: where the argument of the integral does not depend explicitly on the masses of the new particles.
Being Ω DM = m DM Y DM /(s 0 ρ c ), with s 0 being the entropy density at present times while ρ c is the so called critical density, it is easy to see that Ω DM ∝ v 2 Λ 2 , hence without explicit dependence on the DM mass.
In the case in which, instead, the DM is sensitively lighter than the mediator, DM production is in the regime dubbed on-shell in [66], in which the decay of the mediator into DM pairs is the most relevant effect for the DM relic density. In such a case we have that: which simply gives: where y = m φ /T . Using the latter as independent variable in place of the temperature we can write: where we have posed Γ(φ → XX) = κm φ v 2 Λ 2 and Γ(φ→DM DM ) Γtot ∼ 1 since the decay rate of the mediator into DM is not suppressed by the scale Λ. It is then immediate to see that the dependence of the DM relic density on the model parameters is of the form,

Conclusion
In this work we have illustrated the phenomenology of a dark model embedded in a S 4 × Z 5 framework. Focusing at first on a simplified low energy model, in which the flavor symmetry has been just used as ansatz for the structure of the coupling with SM quark of a generic spin-0 mediator, we have shown that, in the cases of both scalar and fermionic DM, it is possible to achieve the correct DM relic density and at the same time comply with constraints from DM searches. Furthermore, we have seen that in the case of fermionic DM, there are specific viable regions of parameter space not present in other simplified models. In a more concrete realization of the scenario under consideration, the DM interactions are suppressed by the scale Λ associated to the breaking of the flavor symmetry. If the latter scale is above 10 (50) TeV for scalar (fermionic) DM, a viable phenomenology is obtained only in the so called secluded regime, in which the relic density is obtained mostly through annihilation in mediator pairs. We have finally considered the possibility in which the NP scale is very large, so that the DM was not existing in thermal equilibrium in the Early Universe. The correct relic density is nevertheless achieved through the freeze-in mechanism for a wide range of values of the DM and mediator masses and for Λ 10 11 − 10 12 GeV.
1 ω 4 ω 4 ω ω ω model assumes the existence of a number of new scalar fields whose transformation properties (together with the SM Higgs) are listed in tab. 1. The flavor symmetry is broken at a generic large energy scale Λ by the vacuum expectation values (vevs) of such fields which, in flavor space, point along the following directions: and In the absence of any dynamical reason, we can assume the same order of magnitude for the vevs v i and introduce the order parameter ε = scalar /Λ which, in the flavor model, governs the ratio of the charged lepton masses as well as the relevant NLO corrections to the neutrino mixing matrix that are needed to shift the (simplistic) LO predictions to a matrix compatible with the current neutrino data [67]. Thus, we expect ε ∼ 0.05. For the quark sector, which is the relevant one for this paper, we report in tab. 2 the transformation properties under the flavor symmetry, where Q is a triplet of SU (2) quark doublets. This assignment is enough to correctly reproduce the ratio among the down-type quarks while a small amount of fine-tuning in the Yukawa couplings is needed to accommodate the large top quark mass as well as the (13) and (23) entries of the CKM. This problem can be prevented by allowing a soft hierarchy among the vevs in eqs. (30) and (31) which, however, is a possibility not contemplated in the present paper.
In order to account for DM in our theory we introduce two states, a scalar (S) of fermionic (χ) DM candidate and a mediator field φ. All these new states are assumed to be singlet both under the SM gauge group and the flavor discrete symmetry. In a similar vein as [63] we can write the following lagrangian: For simplicity we will not consider the possibility of mass mixing between the φ and h states. In addition we will assume isoscalar interactions, i.e., equality between the couplings of up and down quarks in each generation. After flavor and electroweak symmetry breaking the previous Lagrangian generates an effective term that corresponds to the term q c h 1 q of eq. (1) and eq. (2), with where the coefficients a i , b i , c i are linear combinations of the y i parameters and the Higgs VEVs, and ε = ϕ /Λ, assuming a common order of magnitude of the flavon vevs. While the h 1 matrix is flavor diagonal at the leading order, all its entries become not null once the next to leading order effects from the corrections to the vacuum alignment of the flavon fields and from higher order operators of the form: are taken into account: Being of higher order in the ε 1 parameter, the off-diagonal entries of the h 1 matrix have a negligible impact in DM phenomenology. For the latter it is then enough to just take its leading order expression. The coefficients a 1 , b 2 , c 3 are combinations of SM Yukawa couplings and are expected to be of order 1. For simplicity we will assign a 1 = b 2 = c 3 = 1 throughout our study.

B. The Group S 4
The structure of h 1,2 used in our numerical simulations has been obtained adopting the following convention for the generators S and T , according to: In the different representations, they can be written as reported in tab. 3. In the previous basis, the  Clebsch-Gordan coefficients are as follows (α i indicates the elements of the first representation of the product and β i the second one): 1 1 ⊗ η = η ⊗ 1 1 = η with η any representation The multiplication rules with the 2-dimensional representation are the following: The multiplication rules with the 3-dimensional representations are the following: 2 ∼ α 2 β 2 + α 1 β 3 + α 3 β 1 α 3 β 3 + α 1 β 2 + α 2 β 1