The Dark Side of 4321

The evidence of Dark Matter (DM) is one of the strongest observational arguments in favour of physics beyond the Standard Model. Despite expectations, a similar evidence has been lacking so far in collider searches, with the possible exception of $B$-physics discrepancies, a coherent set of persistent deviations in a homogeneous dataset consisting of $b \to c$ and $b \to s$ semi-leptonic transitions. We explore the question whether DM and the $B$ discrepancies may have a common origin. We do so in the context of the so-called 4321 gauge model, a UV-complete and calculable setup that yields a $U_1$ leptoquark, the by far most successful single mediator able to explain the $B$ anomalies, along with other new gauge bosons, including a $Z^\prime$. Adding to this setup a 'minimal' DM fermionic multiplet, consisting of a ${\bf 4}$ under the 4321's $SU(4)$, we find the resulting model in natural agreement with the relic-density observation and with the most severe direct-detection bounds, in the sense that the parameter space selected by $B$ physics is also the one favoured by DM phenomenology. The DM candidate is a particle with a mass in the WIMP range, freeze-out dynamics includes a co-annihilator (the 'rest' of the ${\bf 4}$ multiplet), and the most important gauge mediator in the DM sector is the $Z^\prime$.

Apart from the hints of NP provided by the B-meson anomalies, there are other observations that suggest an extension of the SM. One of the most solid indications of physics beyond the SM is provided by the strong evidence for the existence of Dark Matter (DM) [56,57]. An immediate question is thus whether any of the new heavy vector bosons in 4321 models could be related to the generation of a DM thermal relic. More specifically, we would like to address the possibility that these vector bosons serve as mediators between SM fermionic currents and a DM current. The latter current may be either bosonic or fermionic. However, a scalar DM candidate can annihilate to SM particles via a Higgs portal such that the DM phenomenology would not rely on the new vector bosons. Therefore, we restrict the discussion to fermionic DM. 1 To be specific, we consider a fermionic DM candidate χ 0 that fulfils the following assumptions (cf. e.g. [58,59]) (i) it is a thermal relic, (ii) it is colorless and electrically neutral, (iii) it has zero hypercharge to avoid direct-detection bounds, (iv) it is the component of a massive fermion multiplet Ψ DM that is vector-like (VL) under the 4321 gauge symmetry, (v) (co-)annihilation proceeds via 2 → 2 processes induced at tree level through the new vector bosons U 1 , G , and Z .
These assumptions put restrictions on the possible 4321 quantum numbers of Ψ DM . Conditions (ii) and (iii), i.e. zero electric charge and zero hypercharge of χ 0 require that Ψ DM transforms under an SU (2) L representation with odd dimension. Condition (iii) further implies that the χ 0 eigenvalue of the hypercharge generator Y vanishes. The definition of Y in terms of the U (1) X charge X and the diagonal SU (4) generator T 15 , Y = X + 2 3 T 15 then fixes X for a given SU (4) representation.
Guided by minimality, we restrict our discussion to singlets of SU (3) and to the smallest non-trivial representation of SU (4), the fundamental 4 representation, which leads to 2 Ψ DM ∼ (4, 1, N , +1/2), N ∈ {1, 3, 5, ...} under the 4321 gauge group. After the 4321 symmetry breaking, the Ψ DM multiplet splits into the two components χ and ψ, which transform under the SM gauge group as The dark matter candidate χ 0 is then identified with the electrically neutral component of the SU (2) L N -plet χ. For N = 1 and N = 3, renormalizable couplings between SM particles and the DM candidate exist, such that the latter is in general not stable on time scales vastly below the age of the Universe [58]. In such cases one has to advocate extra symmetries in order for it to be a viable relic. Within our setup, in the N = 1 case the field ψ has the same quantum numbers as right-handed up-type quarks and mixing between ψ R and u i R make χ unstable. In the N = 3 case, a coupling of χ to the Higgs and lepton doublets is allowed such that the DM candidate could decay to a Higgs and a neutrino. The smallest N for which the DM candidate is stable because of the absence of renormalizable couplings that would allow it to decay is N = 5 [58]. 3 In the rest of this paper, we will analyze the phenomenology of all the N = 1, 3, 5 cases, bearing in mind that N = 1, 3 require the additional assumption that renormalizable couplings that destabilize DM are absent. In Sec. 2 we will describe our model setup, paying particular attention to the fermionic sector, that includes the DM multiplet. An extended discussion about the different possibilities for implementing the SM fermions in such a setup is included in Appendix C. Sec. 3 discusses our analytic approach towards estimating the DM relic within our model, including in particular the impact of mass splittings between the DM and its co-annihilator, and our procedure towards estimating the thermally averaged cross section. Mass splittings are discussed within a more general approach in Appendix A, and the cross sections relevant for the thermal average are collected in Appendix B. In Sec. 4 we then move on to describe our approach towards the estimate of direct-detection signals. Sec. 5 collects our results, addressing the question to what extent B-physics discrepancies and DM phenomenology are compatible with one another within our setup. We conclude in Sec. 6.

Model setup
We consider a '4321' model [44,45] based on the gauge group At a scale 4 v LQ , the spontaneous breaking yields the SM color times hypercharge factors. Given the SU (3) 4 × U (1) 4 subgroup of SU (4), the spontaneous breaking proceeds such that SU (3) c is the diagonal subgroup of

Vector bosons
Following the notation of [51], we denote the gauge fields of SU (4), SU (3) , and U (1) X by H α µ , C a µ , and B µ , respectively. The spontaneous breaking yields a massive U 1 LQ 3 Other less minimal scenarios that even for N = 1, 3 do not contain renormalizable couplings that destabilize DM would require ΨDM to transform under larger SU (4) representations or under non-trivial representations of both SU (4) and SU (3) . 4 If the breaking of the gauge group is due to vacuum expectation values (VEVs) of scalar fields, different scalars can contribute to the breaking at slightly different scales (cf. e.g. [47,51]). In order to reduce the number of parameters, we consider only a single breaking scale (as predicted by the model in [52]), and we do not specify the exact mechanism that triggers the spontaneous breaking.
as well as the massive Z µ and the massive gluon-like 'coloron' fields G a µ , given by the linear combinations whereas the linear combinations orthogonal to Z µ and G a µ are the massless hypercharge and QCD gauge bosons B µ and G a µ . Denoting the gauge couplings of SU (4), SU (3) , U (1) X , SU (3) c , and U (1) Y by g 4 , g 3 , g 1 , g s , and g Y , the angles θ 41 and θ 43 are defined analogously to the weak-mixing angle by Since the couplings g s and g Y are the known QCD and hypercharge couplings, eq. (7) implies that for a given value of g 4 , the other two couplings g 1 and g 3 are fixed. Consequently, the gauge sector of the model can be parameterized by only the two independent parameters v LQ , g 4 .
The masses of U µ , Z µ , and G µ are given by and they are related through the angles θ 41 and θ 43 as 5 implying that, at tree level, the U 1 is expected to be the lightest new vector boson.

Fermions
Among the different possibilities for implementing the SM fermions in a 4321 model (see Appendix C), a well-motivated and phenomenologically successful variant corresponds to a unification of third-family quarks and leptons [48][49][50][51][52]. In this case, the first and second families of SM-like fermions transform under the SU (3) × SU (2) L × U (1) X subgroup of the 4321 symmetry like the usual SM fermions, whereas the third-family quarks and leptons are unified into Ψ 3 , which transform under the 4321 symmetry as shown in table 1. Due to their quantum numbers, the light SM fermions cannot directly couple to the U 1 . However, small but non-vanishing couplings between the U 1 and light SM fermions are required to explain the B-meson anomalies. To realize this, we introduce two massive fermions that couple to the U 1 and mix with the left-handed first and second generation SM-like fermions. In addition to couplings between light fermions and the U 1 , whose sizes are controlled by the mixing, this construction also generates the 2-3 entries in the CKM matrix. The new heavy fermions transform in the same way as Ψ 3 L (cf. table 1) and we denote their left-handed components by Ψ 1,2 L . While the mixing is important for the couplings of the SM fermions to U 1 , Z , and G , we do not further discuss the new heavy fermion mass eigenstates since they are not relevant for the 5 Note that this relation can be modified if the 4321 symmetry is broken by different scalars at slightly different scales, cf. e.g. [47,51]. .
DM dynamics as long as their masses are larger than M χ + M ψ , which we assume in the following. Due to the mixing, the first and second generation SM SU (2) L doublets are in general linear combinations of the SM-like fields q 1,2 L , 1,2 L and the new heavy fields Ψ 1,2 L . To avoid large flavour violating effects, we align the mixings between SM fermions and new heavy fermions in the basis in which the down-quark mass matrix is diagonal (cf. e.g. [47]) such that the mixings are flavour-diagonal for the fields where V is the CKM matrix and u i , d i , e i , and ν i are mass eigenstates. A possible misalignment between the quark and lepton components of the fields Ψ i L is parameterized by embedding the quark and lepton components Ψ 1,2 q L and Ψ 1,2 L that have a flavour-diagonal mixing with q 1,2 L and 1,2 L , respectively, as where W is a unitary matrix parameterizing the misalignment. This matrix is usually chosen to be CP -conserving and to mix only the second and third generation, i.e. we use In the absence of additional new heavy fermions that mix with right-handed SM fermions, a possible quark-lepton misalignment in Ψ + 3 R and Ψ − 3 R corresponds to only a phase difference,  Table 2: Constants κ and ξ entering the couplings of fermions to G and Z (cf. eq. (16)).
which we parameterize as Consequently, the SM fields in the basis where the down-quark mass matrix is diagonal can be expressed as In this basis, the couplings of the new vector bosons U 1 , Z , and G to the SM fermions and to the DM-sector fields ψ and χ are given by where the constants κ and ξ that appear in the G and Z couplings are collected in table 2 and the constants β that appear in the U 1 couplings are given by Note that in the limit of large g 4 , where cos θ 41 ≈ cos θ 43 ≈ 1 and sin θ 41 ≈ sin θ 43 ≈ 0, the constants κ and ξ are approximately i.e. the couplings of left-handed light fermions are proportional to their mixings with the new heavy fermions, the couplings of right-handed light fermions vanish, and all couplings of third-generation SM fermions and DM sector fields satisfy κ ≈ ξ ≈ 1. The constants β, on the other hand, are independent of the value of g 4 and only depend on fermion mixing angles and phases. While the parameterization described above allows explaining the B-meson anomalies and avoids strong constraints from large flavour violating effects, the number of parameters can be further reduced by the following phenomenologically motivated assumptions: • To maximize the agreement with the B-decay measurements that deviate from the SM, one can take β de = −1 [51], which fixes the phase φ e = π.
• Since the phase φ ν is currently not constrained by any measurement one can use φ ν = 0 for simplicity.
• An approximate U (2) symmetry in the quark sector, i.e. θ q 1 ≈ θ q 2 , can be employed to suppress tree-level FCNCs in the up-quark sector that are mediated by the Z and G [47]. Without such a U (2) protection, excessive contributions to ∆C = 2 observables would be possible.
• The first-generation lepton doublet can be taken to be purely a singlet of SU (4), i.e. θ 1 = 0, to be safe from LFV due to U 1 couplings involving the electron.
Making all of the above assumptions and defining θ q 12 = θ q 1 = θ q 2 , the only remaining free parameters in the fermion sector are where M χ denotes, here and henceforth, the mass of the DM-candidate 6 χ and N is the dimension of its SU (2) L representation.

Parameter ranges
As summarised in eqs. (8) and (19), within reasonable assumptions the 'effective' model parameters are the following: In this section we would like to collect the non-negligible information available on these parameters from a fit to flavour data as well as from constraints due to direct searches. In Sec. 5 we will then address the question to what extent these constraints are compatible with those coming from cosmological and direct-detection information about Dark Matter. The above parameters can be grouped into three classes according to their impact on the DM phenomenology: The DM phenomenology depends crucially on these parameters. 6 The mass Mχ is related to the tree-level mass of the multiplet ΨDM as discussed in Appendix A.
• θ q 12 : This parameter is important only for DM direct detection.
• θ 2 , θ LQ : The DM phenomenology is essentially independent of these parameters.
Within our model, a combination of the parameters v LQ , θ q 12 , θ 2 , and θ LQ is constrained by B-physics data alone, in particular by the R(D ( * ) ) discrepancies. A global fit of v LQ and the constants β (cf. eq. (17)) was performed in Ref. [51]. Expressed in our notation, the fit prefers values for v LQ in the range v LQ / cos θ LQ ∈ [3.1, 4.6] TeV, with cos θ LQ ≈ 0.8 − 0.9. While the preferred value of v LQ is correlated with the preferred values of θ 2 and θ LQ , our DM phenomenology is essentially independent of the latter two parameters. Consequently, for any reasonable value of v LQ , we can set the parameters θ 2 and θ LQ to comply with the fit in [51], while fulfilling all DM constraints. The fit to B-physics data leaves some freedom for θ q 12 , which, however, is constrained by direct searches (see below). In short, for definiteness we take v LQ ∈ [3, 5] TeV (21) as our fiducial range for v LQ . The parameters g 4 and θ q 12 enter the definition of the fermionic-currents' couplings to the U 1 , the Z and the G , which are constrained by direct searches. In fig. 1, we show the g 4 dependence of light-quark couplings to the Z as well as to the G , for different values of sin θ q 12 . This dependence displays transparently the g 4 and sin θ q 12 ranges preferred by direct searches. The figure shows at a glance that an efficient suppression of these coupling combinations is achieved for small sin θ q 12 and large g 4 . Representative ranges are sin θ q 12 0.2 and g 4 3 .
We conclude that agreement with B-decay discrepancies and with the constraints coming from direct searches select the parameter regions in eqs. (21)- (22). In Sec. 5 we will confront this parameter space with the constraints imposed by Dark Matter.

Dark matter relic abundance
In this section we shall address the question whether our setup, as introduced around eq. (2) and discussed in Sec. 2, can accommodate the relic abundance of DM observed today, Ω 0 h 2 .
In addition to the DM candidate χ 0 , the DM sector of our model includes also several co-annihilation partners: all the other components, charged under weak isospin, of the χ and ψ SU (2) L multiplets. As shown in classic work, even in the presence of co-annihilators an estimate of Ω 0 h 2 accurate to about 10% may be obtained analytically [65] (see also [66]). This accuracy is satisfactory in our case, in view of several uncertainties inherent in the problem and likewise discussed in the above works.
The first main step towards the estimate of Ω 0 h 2 is the determination of Ωh 2 at the 'freeze-out' temperature T f . It is convenient to introduce the variable x denoting the inverse temperature in units of the DM mass, i.e. x ≡ M χ 0 /T , and to define x f ≡ x| T =T f . In the case -like ours -where co-annihilators are present, x f is determined iteratively from the relation [65] x where M Pl = 1.22 × 10 19 GeV, g * denotes the total number of effectively relativistic d.o.f. at freeze-out, g eff denotes the number of effective d.o.f. within the DM sector, and the thermally averaged annihilation cross section σ eff v is the main dynamical quantity. After freeze-out, the relic abundance is subject to post-freeze-out annihilation processes. The efficiency of this post-freeze-out annihilation is given by [65,67] J The present-day DM abundance can then be estimated as 8 The crucial ingredient in the determination of Ω 0 h 2 is the calculation of the thermally averaged annihilation cross section σ eff v , which enters in Ω 0 h 2 through the post-freeze-out annihilation efficiency J. To this end, it is also necessary to determine g eff , which enters in both σ eff v and x f . In turn, g eff and σ eff depend in an important way on the mass differences between the DM candidate and its co-annihilation partners [65]. In the following, we discuss in great detail the mass differences, the effective degrees of freedom g eff , and therewith proceed to the estimate of the thermally averaged annihilation cross section σ eff v .

Mass differences
For mass differences between the DM candidate and its co-annihilation partners comparable to the freeze-out temperature (which quantifies the average amount of kinetic energy available in the collisions), the co-annihilation partners become nearly as kinematically accessible as the DM candidate. The mass differences are therefore an important ingredient in the determination of the annihilation cross section. In our model, the DM candidate and its co-annihilation partners are components of VL fermion multiplets that transform nontrivially under the SU (2) L and SU (4) gauge groups. The relevant mass differences in this case correspond to the mass splitting inside these multiplets, which are generated after the spontaneous breaking of the SU (2) L and SU (4) symmetries. While it is possible to generate a mass splitting at tree level, e.g. by coupling the VL multiplets to the scalar operators responsible for the spontaneous symmetry breaking, a mass splitting is generated even in the absence of such tree-level terms. At the one-loop order, the gauge bosons associated with the spontaneously broken symmetries, of which some become massive due to the breaking, induce corrections to the fermion masses. Since the components of the VL multiplets correspond to different irreducible representations of the unbroken gauge group, each of them couples differently to the gauge bosons and thus receives a different contribution to its mass. It is possible to obtain a generic result for the one-loop mass splitting among components of a VL multiplet that is applicable to a large set of spontaneously broken gauge groups (see appendix A). Applying our generic result, eq. (62), to the EW gauge group, we can determine the relative mass difference between components ξ and η of a VL multiplet of hypercharge Y and massM . We find where Q ξ and Q η are the electric charges of ξ and η, and f (r) is a finite loop function given in eq. (58). This reproduces the well-known result for the mass splitting in EW VL multiplets (cf. e.g. [58]). For reference, the relative mass splitting within the ψ and χ SU (2) L multiplets of our DM sector is between O(10 −3 ) and O(10 −4 ) forM = O(1 TeV).
Having at hand the generic result, eq. (62), it is straightforward to determine the relative mass splitting between our DM candidate χ and its colored co-annihilation partner ψ. This mass splitting is induced by the vector bosons associated with the 43(2)1 symmetry breaking and is given by The value of ∆ 4321 ψχ is around 8 − 15% for the parameter region of interest (see fig. 2) and its significance is further discussed in sections 3.2 and 3.3.
Additional mass splittings are induced by the non-zero temperature at which our processes of interest take place. We estimate these mass splittings to be of order ∆ EW (T ) ∼

Effective degrees of freedom
In our case, the number of effective d.o.f. in the DM sector g eff is given by [65] where the index i runs over the N components of the SU (2) L multiplets χ and ψ. Besides g χ = 4 and g ψ = 12 denote the internal (spin, color, ...) d.o.f. of the components of these multiplets. 9 The relative mass splittings ∆ χ i and ∆ ψ i are defined as The relative mass differences within the χ and ψ multiplets are only at the per mil level (see discussion below eq. (27)). We thus neglect them in the following and use a common mass and relative mass splitting for each multiplet, such that M ψ = (1 + ∆ ψ ) M χ . Employing this approximation, the number of effective d.o.f. simplifies to We see that, at the freeze-out temperature, g eff departs appreciably from N (g χ + g ψ ) unless x f ∆ ψ 1, i.e. unless the ψ-χ mass splitting is much smaller than the freeze-out temperature. Combining eq. (23) with the value of the ψ-χ mass splitting discussed in Sec. 3 eqs. (32)-(33) provide an accurate determination of g eff within our setup.
9 Since χ, ψ belong to complex VL representations of the gauge group, they are Dirac fermions.

Thermally averaged cross section
The main dynamical quantity in eq. (23) is σ eff v , where [65] with e.g.
and the other cross sections defined analogously. Here X, X denote any particles other than χ, ψ. Since we assume that the DM sector is lighter than any of the U 1 , Z , G mediators or new vector-like fermions, for the X, X we only consider SM particles. Neglecting again the mass differences within the ψ and χ multiplets, i.e. employing the replacements in eq. (31), the effective cross section simplifies to The cross sections σ χ i χ j , σ χ i ψ j , and σ ψ i ψ j are due to the exchange of either SM bosons or the new heavy gauge bosons U 1 , Z , and G . Because of the dependence of the cross sections on the fourth power of the couplings, contributions due to the electroweak sector are negligible compared to those involving the relatively strongly coupled new heavy gauge bosons or gluons. We find that all the cross sections mediated by U 1 , Z , G , and gluons are of comparable size. However, in the effective cross section, σ χ i ψ j and σ ψ i ψ j are multiplied by one and two powers, respectively, of a factor that is suppressed by the ψ-χ mass splitting, cf. eq. (33). Consequently, the Z -mediated cross section σ χ i χ j is larger than any other contribution to σ eff by one to two orders of magnitude. In view of the overall 10% uncertainty in our analytical estimate, we can therefore approximate where we have used eqs. (32)- (33) and The thermal average σ eff v can be determined as an expansion in powers of 1/x, and this expansion can be related order by order to the expansion of σ eff around s = 4M 2 χ [68]. To this end, we follow the notation of [69]. Neglecting the masses of the annihilation products, we define σ 0 (y) = 2 y 2 − y σ eff (y) where we substituted s by the dimensionless variable y = s/(4 M 2 χ ). The thermal average σ eff v can then be expressed as [68,69] where the first coefficients are and we have defined The above relations allow thus to verify that, by including higher powers in the small-velocity expansion of σ 0 (y), higher powers in the small-temperature expansion of σ eff v are smaller and smaller.

Present-day DM abundance
In order to obtain the present-day DM abundance, one convolutes the calculated σ eff v in the post-freeze-out annihilation efficiency J in eq. (24). Using the expansion in eq. (39) to order 1/x 2 f , the present-day DM abundance in eq. (25) yields We will perform a full numerical study in Sec. 5, following the discussion of the constraints imposed by direct detection in Sec. 4. Here we would like to make a few qualitative considerations around eq. (42). Using σ eff as in eq. (37), eq. (38) yields σ 0 (y) = 1 128 π g 4 cos θ 41 where N = 1, 3, 5, ... denotes the SU (2) L size of the χ and ψ multiplets and for brevity we introduced the flavour function Plugging eq. (43) into eq. (42), and taking the representative value g 4 = 3, we find A few remarks are in order. First, eq. (45) assumes 4M 2 χ M 2 Z . For v LQ = 3 TeV (5 TeV), this approximation implies an error 15% ( 40%) in eq. (45), keeping in mind that the preferred range for M χ are 600 GeV ( 1.5 TeV), see Sec. 5. Second, with the above mass ranges at hand, we can discuss the relative size of the corrections due to the c 1,2 terms in the 1/x f expansion (see eq. (42) The procedure outlined in this section, and leading to eq. (42), with σ eff including U 1 , Z , G and gluon contributions, will be used in Sec. 5 to identify the regions of parameter space that are viable in the light of all constraints, including B discrepancies, the relic abundance and also direct-detection constraints, to be discussed in the next section.
The approximate formula in eq. (45), and the discussion around it, demonstrate that Ω 0 h 2 of the order of the observed value can be obtained without effort, in compliance with all other constraints.

Dark-Matter Direct Detection
One of the most straightforward signals one may expect of our model are DM collisions on nuclei. The latter are constrained by a large number of direct-detection experiments, the most stringent bounds for the DM masses of interest to us being Refs. [70][71][72]. Actually, it is precisely in the light of these constraints that we restricted our attention to DM multiplets that transform under SU (2) L representations with odd dimensions, allowing for a Y = 0 multiplet member -the DM candidate -as discussed in Sec. 1.
Even in our case however, the DM-nucleon cross section receives a tree-level contribution mediated by a Z . Since DM is non-relativistic, and since M Z is also much larger than the relevant momentum transfer, the scattering process with the nucleon constituents may be accounted by a local Lagrangian Furthermore, if we are able to neglect corrections due to the finite momentum transfer between the DM and the nucleons, we may parametrize the matrix elements between vector or axial-vector quark q currents and the external-state nucleons N as For the two form factors at zero momentum transfer we follow conventions common in the literature. In particular, F q/N 1 (0) counts the number of valence quarks q in the nucleon N , e.g. F d/n 1 (0) = 2. eqs. (46)-(47) yield the following spin-independent cross section for elastic scattering between DM and a single nucleon N = p or n where and Starting from eq. (48), in order to estimate the matrix element on a nucleus N with mass number A and atomic number Z, one may assume (see e.g. [73,74]) that DM scatters coherently on the A nucleons of the target. In the static limit, the DM -nucleon cross section measured by experiments operating with nuclei N as target material can thus be estimated from eq. (48) with the replacement The above procedure is crude in a number of ways, that have been amply discussed in the literature [74][75][76][77][78][79][80][81][82][83][84][85][86][87][88][89][90][91]. A first outstanding limitation is the fact that the matching scale for the interactions in eq. (46) is well above the effective scale for the hadronic matrix elements in eq. (47), hence renormalization-group effects are in general non-negligible, and the relativistic-operator basis of eq. (46) has to be matched onto the non-relativistic basis relevant for the interaction with the nucleus. A second crucial limitation occurs if the DM momentum is large enough that the pointlike-nucleon approximation inherent in eq. (47) loses validity.
The impact of the above approximations may be explored using the public codes DirectDM [90][91][92] -that accounts for renormalization-group effects from the UV scale to the scale of the (non-relativistic) interaction with the nucleons -and DMFormFactor [74,76,77] -which estimates the non-trivial dynamics due to non-negligible momentum transfers between the DM and the nucleon.
We performed a detailed comparison between the prediction obtained within the analytic approach of eqs. This comparison is in agreement with the expectation that, for M χ light enough, the DM Compton wavelength -of order v 0 /c × M χ ≈ 10 −3 M χ , where v 0 is the typically assumed RMS velocity of the DM halo distribution -is not sufficiently large to resolve the inner nucleon structure, so that the pointlike-nucleon approximation is tenable.
We conclude that, for M χ masses in the range required by the relic-density constraint, 300 GeV−1.5 TeV, use of the analytic prediction of eqs. (48)-(51) will produce directdetection bounds that are somewhat stronger -yet quite realistic -than those produced with numerical codes. In Sec. 5 we will compare this analytic prediction with the latest bound obtained by the Xenon1T experiment [72].

Results
The relic-density and direct-detection constraints discussed in the previous two sections represent significant phenomenological input for our model. We summarised our parameter space in eq. (20) and discussed how B-decay discrepancies and collider constraints lead to the preferred ranges in eqs. (21)- (22). In this section we discuss to what extent such ranges are compatible with those imposed by the DM relic-density and direct-detection constraints.
Quite remarkably, the g 4 and sin θ q 12 ranges in eq. (22) are also favoured by direct-detection constraints, as illustrated in fig. 3. Here we show the DM -nucleon cross section σ N SI discussed around eq. (48), as a function of g 4 for increasing values of sin θ q 12 ≥ 0. As discussed in Sec. 4, experiments yield severe limits, as strong as σ N SI < 10 −45 cm 2 . As the figure shows, these limits can be comfortably satisfied with the choice sin θ q 12 0.2 and g 4 3. 11 Therefore, the constraints from Z and G direct searches on the one side, and from DM direct detection on the other side, identify one and the same region for g 4 and sin θ q 12 . Although a correlation between the suppression required to DY-produced Z and the suppression required to Z -mediated DM -nucleon scattering may be expected just by crossing symmetry, the coupling combinations involved in the two processes are entirely different. Besides, it is non-trivial that couplings in compliance with Z searches would also yield a DM cross section on nucleons as small as 10 −45 cm 2 . . v LQ plane. Coloured rays fulfil the Ω 0 h 2 constraint within ±15%. Grey regions, with sin θ q 12 fixed at the displayed value, are excluded by Xenon1T [72]. See text for details.
We next discuss the behaviour of sin θ q 12 , g 4 and v LQ as a function of the DM mass M χ , when the relic-density and direct-detection constraints are imposed. In fig. 4 we show as colored rays the regions selected by the constraint Ω 0 h 2 in the M χ vs. v LQ plane. The constraint is imposed with ±15% accuracy, corresponding to the error we attach to its calculation -see discussion in Sec. 3. The different rays refer to different choices of N , whereas g 4 = 3 following our above discussion. We see that the v LQ range in eq.
Also shown in different shades of grey depending on sin θ q 12 ∈ {0.2, 0.25, 0.3} are regions excluded by direct detection. It is clear that, for sin θ q 12 = 0.2 and v LQ ≥ 3 TeV (see eq. (21)), only a tiny fraction of the parameter space is excluded by direct detection, and mostly for N = 1. It is interesting to also test the dependence of the region selected by the Ω 0 h 2 constraint on g 4 , rather than on v LQ . We show such dependence in the two panels of fig. 5, corresponding to v LQ set to 3 and 5 TeV, respectively. We see that the dependence on g 4 is actually very weak, especially if this coupling is large. The figure also shows the region where M χ ≥ M U /2, that we excluded for simplicity. In fact, as the U 1 becomes on-shell, 2 → 3 and 2 → 4 decay channels open up, whereas we restricted to 2 → 2 processes in the calculation of σ eff .
The dependence of the Ω 0 h 2 constraint on M χ vs. v LQ or g 4 shown in figs. 4 and 5 can be captured simultaneously in fig. 6, that shows this constraint in the plane v LQ vs. g 4 for a few reference values of M χ and of N , represented as isolines. The figure shows at a glance that within the fiducial v LQ range of eq. (21), the relic-density constraint can be comfortably satisfied whatever the choice of N , and also nearly irrespective of the choice of g 4 -which, as we discussed, is instead constrained by direct detection.
Finally, fig. 7 displays the allowed parameter space in the M χ vs. sin θ q 12 plane, with g 4 = 3. The coloured 'curtains' show the region excluded by direct detection, depending on the choice of N . Similarly as fig. 3 and the ensuing discussion, this plot shows that direct detection tends to prefer small sin θ q 12 , the upper bound for a given M χ becoming stronger as N increases. As discussed earlier, imposing the relic-density constraint plus the v LQ range suggested by B discrepancies yields the indicative M χ ranges in eq. (52). These ranges are not reported in fig. 7 to limit clutter.

Conclusions
We investigate a possible common description of the only hint of physics beyond the SM in collider searches -the B-decay discrepancies -and of one of the strongest phenomenological indications of new physics -the existence of Dark Matter. 12 We adopt the 4321 gauge ansatz, a well-motivated, UV-complete, calculable setup for explaining the B anomalies based on the gauge group SU (4) × SU (3) × SU (2) L × U (1) X . To this ansatz, we add a minimal DM sector, represented by a fermionic multiplet sitting in the fundamental 4 of SU (4) and in a representation of SU (2) L with odd dimension. After the breaking to the SM group, this multiplet gives rise to a DM candidate and a coloured co-annihilation partner with mass ≈ 10% larger.
We find that the parameter space selected by collider and B-physics data is also the one favoured by DM phenomenology, in particular by the constraints imposed by direct detection. These parameter choices include the 43(2)1 symmetry breaking scale v LQ ∈ [3,5] TeV, large SU (4) gauge coupling g 4 ≈ 3 and small fermion mixing parameters of light quarks sin θ q 12 0.2. Within this parameter space, direct-detection signals happen to be below the severe bounds imposed in particular by Xenon1T. For v LQ in the above mentioned range, the requirement of the correct DM relic density is easily fulfilled with a DM mass between about 250 GeV and 1.5 TeV, depending on the SU (2) L representation, and on the v LQ value.
In short, we find our setup neatly compatible with the most accurate DM observations -the relic density and the limits imposed by direct detection. Interestingly, while the new particle in 4321 models that is mainly responsible for B-physics discrepancies is the U 1 LQ, it is the Z that plays the leading role in DM phenomenology.
The study of this class of models thus warrants further scrutiny, especially if the B anomalies will be consolidated by forthcoming measurements.

Acknowledgments
DG would like to thank Marco Cirelli for very useful feedback, and Eugenio Del Nobile for critical comments on Sec. 4. Exchanges with Fady Bishara, Jacopo Ghiglieri and Liam Fitzpatrick are also acknowledged. MR thanks the Université de Montréal for kind hospitality in Spring 2019, and for discussions with David London and Jacky Kumar on related topics. This work is supported by an ANR PRC (contract n. 202650) and by the Labex Enigmass.

A. Mass splitting
In order to determine the one-loop mass splitting, we compute the pole masses of the components of a VL fermion multiplet. It is possible to obtain these pole masses in a way that is applicable to a large set of different spontaneously broken gauge groups. To this end, we consider a gauge group G × H that is spontaneously broken to its subgroup H, where G has a subgroup H G ⊆ G that is isomorphic to both H and the unbroken H, i.e. H ∼ = H ∼ = H G , and H is the diagonal subgroup of H G × H . We take G to be a simple group and H to be semi-simple and given by the direct product of simple groups H i as H = H 1 × H 2 × · · · × H n . 13 An analogous decomposition into simple groups H i and H G i applies to H and H G . We denote the gauge couplings of G, H i , and H i by g G , g H i , and g H i , respectively, and we note that the gauge couplings of all simple group factors of H G are equal to g G . It is convenient to define mixing angles θ i by where k i is a normalization factor relevant in case of an abelian group H i = U (1) and corresponds to the normalization of the U (1) charges. In particular, we denote them as follows: 2. From the massless vector bosons that transform in the adjoint representation of H, • For the EW gauge group, the unbroken group is abelian. Thus, the quadratic Casimir invariants are simply given in terms of squares of U (1) charges. This yields and we find which coincides with the well-known result (cf. e.g. [58]).
• For the 43(2)1 gauge group, the unbroken group contains one abelian and one nonabelian factor. The quadratic Casimir invariants are i.e. the Casimir invariants of the abelian groups can be expressed by the U (1) charges Y and Y , while those of the non-abelian factors are SU (3) Casimir invariants.
For our DM candidate χ and its co-annihilation partner ψ, the Casimir invariants are given by

C. Fermions in 4321 models
The fermion sector of the model contains the SM fermions, the ψ and χ, as well as other heavy vector-like (VL) fermions that mix with the SM fermions. While the mixing of the SM fermions with the heavy VL fermions is important for the couplings of the SM fermions to U 1 , Z , and G , the VL fermions themselves are not relevant for the DM dynamics as long as their masses are larger than M ψ + M χ , which we assume in the following.  (4) and SU (2) L representations, respectively. Ψ DM is the multiplet containing χ and ψ.

C.1. The SM fermions
Due to the mixing with the heavy VL fermions, the SM SU (2) L doublets are in general linear combinations of the fields Ψ 4,2 and Ψ 1,2 shown in table 3, while the SM SU (2) L singlets are linear combinations of Ψ 4,1 and Ψ 1,1 . To avoid large flavour violating effects, we align the mixings between SM fermions and heavy VL fermions in the basis in which the down-quark mass matrix is diagonal (cf. e.g. [47]) such that the mixings are flavour-diagonal for the fields where V is the CKM matrix. A misalignment between the quark and lepton components of Ψ 4,2 is implemented by embedding the components Ψ q 4,2 and Ψ 4,2 that have a flavour-diagonal mixing with Ψ q 1,2 and Ψ 1,2 , respectively, as where W is a unitary matrix parameterizing the misalignment. For simplicity, no misalignment but only a phase difference is introduced for the quark and lepton components (Ψ u 4,1 and Ψ ν 4,1 ), and (Ψ d 4,1 and Ψ e 4,1 ) of Ψ ↑ 4,1 , and Ψ ↓ 4,1 , respectively, i.e.
The above parameterization is general enough to recover the couplings between SM fermions and the heavy vector bosons in several 4321 models in the literature. In particular, the following special cases can be considered.
• Traditional 4321 models: In "traditional" 4321 models [45,47], all three generations of left-handed SM fermions are each a mixture of a 4 and a 1 of SU (4), while all righthanded SM fermions are purely singlets of SU (4). This corresponds to the choice sin θ u i = sin θ d i = sin θ e i = sin θ ν i = 0, The misalignment matrix W is usually chosen to be CP -conserving and to mix only the second and third generation, i.e.
The number of parameters can be further reduced by the following phenomenologically motivated assumptions [47]: -A U (2) symmetry in the quark sector, i.e. θ q 1 = θ q 2 , can be employed to suppress tree-level FCNC in the up-quark sector that are mediated by the Z and G . Without such a U (2) protection, excessive contributions to ∆C = 2 observables would be possible.
-The first-generation lepton doublet can be taken to be purely a singlet of SU (4), i.e. θ 1 = 0, to be safe from LFV due to U 1 couplings involving the electron.
Making the above described assumptions to reduce contributions to ∆C = 2 observables and LFV electron couplings, the set of free parameters in the fermion sector is reduced to θ q 12 , θ 2 , θ LQ .
While the above couplings are independent of the representation N of SU (2) L under which Ψ DM transforms, the couplings to the W and Z bosons are clearly different for different representations. The coupling of Z to a field Ψ N transforming as a N of SU (2) L and having hypercharge Y is given by where T 3 N is the diagonal generator of SU (2) in the N representation and the electric charge is defined by Q = T 3 N + Y . The coupling of W ± to a field Ψ N transforming as a N of SU (2) L is derived from the covariant derivativē where T ± N = T 1 N ± i T 2 N .