# The Dark Side of 4321

Diego Guadagnoli, Meril Reboud, Peter Stangl

► To cite this version:

Diego Guadagnoli, Meril Reboud, Peter Stangl. The Dark Side of 4321. Journal of High Energy Physics, 2020, 10, pp.084. 10.1007/JHEP10(2020)084. hal-02863106

HAL Id: hal-02863106

https://hal.science/hal-02863106v1

Submitted on 2 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The dark side of 4321

## Diego Guadagnoli, Méril Reboud and Peter Stangl

LAPTh, Université Savoie Mont-Blanc et CNRS,
Annecy, France

E-mail: diego.guadagnoli@lapth.cnrs.fr, meril.reboud@lapth.cnrs.fr,
peter.stangl@lapth.cnrs.fr

**ABSTRACT:** The evidence of Dark Matter (DM) is one of the strongest observational arguments in favor 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'$. Adding to this setup a ‘minimal’ DM fermionic multiplet, consisting of a 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 favored 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 4 multiplet), and the most important gauge mediator in the DM sector is the $Z'$.

**KEYWORDS:** Beyond Standard Model, Cosmology of Theories beyond the SM, Heavy Quark Physics

ARXIV ePRINT: 2005.10117
# Contents

1 Introduction 1

2 Model setup 4

2.1 Vector bosons 4

2.2 Fermions 5

2.3 Parameter ranges 8

3 Dark-Matter relic abundance 10

3.1 Mass differences 11

3.2 Effective degrees of freedom 12

3.3 Thermally averaged cross-section 13

3.4 Present-day DM abundance 14

4 Dark-Matter direct detection 15

5 Results 17

6 Comments on Dark-Matter indirect detection 20

7 Conclusions 22

A Fermions in 4321 models 23

A.1 The SM fermions 23

A.2 The fermions in the DM sector 26

B Mass splitting 27

C Cross-sections of processes entering the estimation of the relic density 31

# 1 Introduction

After the end of LHC Run 2, no sign of new physics (NP) has been observed in direct searches. There are, however, several indirect hints for NP. Flavor physics experiments have reported a large set of deviations from Standard Model (SM) predictions in B-meson decays, which are also known as B-meson anomalies. They amount to discrepancies both in neutral current $b \to s\bar{l}\bar{l}$ decays [1–11] and in charged current $b \to c\bar{l}\nu$ decays [12–21]. It was realized that the presence of a $U_1$ leptoquark (LQ) with SM quantum numbers (3, 1, 2/3) could simultaneously explain both of these sets of discrepancies [22–30]. While other simultaneous solutions are possible (see e.g. [31–37]), an explanation in terms of a
$U_1$ LQ has become even more consistent with recent data [38–40]. Being a massive vector boson, the $U_1$ requires a UV completion from which it arises either as a gauge boson of a spontaneously broken gauge symmetry or as a composite vector boson (see e.g. [41]). In fact, the $U_1$ is a well-known prediction of Pati-Salam models [42], which extend the SM color $SU(3)_c$ gauge group to $SU(4)$, under which quarks and leptons transform in unified multiplets. Since traditional Pati-Salam models cannot accommodate the flavor structure needed for explaining the $B$-meson anomalies, variants of the Pati-Salam models based on the gauge group $SU(4) \times SU(3)' \times SU(2)_L \times U(1)_X$ have been constructed to this end, the so-called 4321 models [43–52]. In these models, the SM arises after the group $SU(4) \times SU(3)' \times U(1)_X$ is spontaneously broken to its $SU(3)_c \times U(1)_Y$ subgroup. The heavy vector bosons resulting from this symmetry breaking include the $U_1$ LQ but in addition also an $SU(3)_c$ octet $G'$ dubbed “coloron” and a SM singlet $Z'$. Among the 4321 models, those that unify the third family of SM quarks and leptons are of special interest since they imply an approximate global $U(2)^5$ flavor symmetry [53–55]. Such a symmetry is particularly useful for explaining the $B$-meson anomalies without violating other flavor bounds while at the same time reproducing the SM fermion masses and CKM elements [48–52].

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.<fnref n="1" />

To be specific, we consider a fermionic DM candidate $\chi_0$ that fulfills 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 $\Psi_{DM}$ that is vector-like (VL) under the 4321 gauge symmetry,

(v) (co-)annihilation proceeds via $2 \to 2$ processes induced at tree level through the new vector bosons $U_1, G', \text{and } Z'$.

These assumptions put restrictions on the possible 4321 quantum numbers of $\Psi_{DM}$. Conditions (ii) and (iii), i.e. zero electric charge and zero hypercharge of $\chi_0$ require that $\Psi_{DM}$

<footnote n="1">¹Other interesting cases that are beyond the scope of the present article include composite bosonic DM that could naturally arise if the 4321 gauge symmetry is broken by a new strong interaction [52].</footnote>
transforms under an SU(2)$_L$ representation with odd dimension. Condition (iii) further implies that the $\chi_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 + \sqrt{\frac{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<fnref n="2" />

$$ \Psi_{\text{DM}} \sim (\mathbf{4}, \mathbf{1}, \mathbf{N}, +1/2), \quad \mathbf{N} \in \{\mathbf{1}, \mathbf{3}, \mathbf{5}, \dots\} \qquad (1.1) $$

under the 4321 gauge group. After the 4321 symmetry breaking, the $\Psi_{\text{DM}}$ multiplet splits into the two components $\chi$ and $\psi$, which transform under the SM gauge group as

$$ \chi \sim (\mathbf{1}, \mathbf{N}, 0), \quad \psi \sim (\mathbf{3}, \mathbf{N}, 2/3), \quad \mathbf{N} \in \{\mathbf{1}, \mathbf{3}, \mathbf{5}, \dots\}. \qquad (1.2) $$

The Dark Matter candidate $\chi_0$ is then identified with the electrically neutral component of the SU(2)$_L$ N-plet $\chi$. For $\mathbf{N} = \mathbf{1}$ and $\mathbf{N} = \mathbf{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 $\mathbf{N} = \mathbf{1}$ case the field $\psi$ has the same quantum numbers as right-handed up-type quarks and mixing between $\psi_R$ and $u_R^i$ make $\chi$ unstable. In the $\mathbf{N} = \mathbf{3}$ case, a coupling of $\chi$ to the Higgs and lepton doublets is allowed such that the DM candidate could decay to a Higgs and a neutrino. The smallest $\mathbf{N}$ for which the DM candidate is stable because of the absence of renormalizable couplings that would allow it to decay is $\mathbf{N} = \mathbf{5}$ [58].<fnref n="3" />

In the rest of this paper, we will analyze the phenomenology of all the $\mathbf{N} = \mathbf{1}, \mathbf{3}, \mathbf{5}$ cases, bearing in mind that $\mathbf{N} = \mathbf{1}, \mathbf{3}$ require the additional assumption that renormalizable couplings that destabilize DM are absent. In section 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 A. Section 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 B, and the cross-sections relevant for the thermal average are collected in appendix C. In section 4 we then move on to describe our approach towards the estimate of direct-detection signals. Section 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 section 7.

<footnote n="2">²If $\Psi_{\text{DM}}$ is a singlet of SU(4), then its 4321 quantum numbers are fixed to $(\mathbf{1}, \mathbf{1}, \mathbf{N}, 0)$, and couplings to $U_1$, $G'$, and $Z'$ are absent. This corresponds to “Minimal Dark Matter”, discussed in [58].</footnote>

<footnote n="3">³Other less minimal scenarios that even for $\mathbf{N} = \mathbf{1}, \mathbf{3}$ do not contain renormalizable couplings that destabilize DM would require $\Psi_{\text{DM}}$ to transform under larger SU(4) representations or under non-trivial representations of both SU(4) and SU(3)'.</footnote>
## 2 Model setup

We consider a ‘4321’ model [44, 45] based on the gauge group

$$SU(4) \times SU(3)' \times SU(2)_L \times U(1)_X. \qquad (2.1)$$

At a scale<fnref n="4" /> $v_{LQ}$, the spontaneous breaking

$$SU(4) \times SU(3)' \times U(1)_X \rightarrow SU(3)_c \times U(1)_Y \qquad (2.2)$$

yields the SM color times hypercharge factors. Given the $SU(3)_4 \times U(1)_4$ subgroup of $SU(4)$, the spontaneous breaking proceeds such that $SU(3)_c$ is the diagonal subgroup of $SU(3)_4 \times SU(3)'$ and $U(1)_Y$ is the diagonal subgroup of $U(1)_4 \times U(1)_X$.

### 2.1 Vector bosons

Following the notation of [51], we denote the gauge fields of $SU(4)$, $SU(3)'$, and $U(1)_X$ by $H_\mu^\alpha$, $C_\mu^a$, and $B'_\mu$, respectively. The spontaneous breaking yields a massive $U_1$ LQ

$$U_{\mu}^{\pm 1,2,3} = \frac{1}{\sqrt{2}} (H_{\mu}^{9,11,13} \mp i H_{\mu}^{10,12,14}), \qquad (2.3)$$

as well as the massive $Z'_\mu$ and the massive gluon-like ‘coloron’ fields $G'_\mu$, given by the linear combinations

$$Z'_{\mu} = H_{\mu}^{15} \cos \theta_{41} - B'_{\mu} \sin \theta_{41}, \qquad G'_{\mu} = H_{\mu}^{a} \cos \theta_{43} - C_{\mu}^{a} \sin \theta_{43}, \qquad (2.4)$$

whereas the linear combinations orthogonal to $Z'_\mu$ and $G'_\mu$ are the massless hypercharge and QCD gauge bosons $B_\mu$ and $G_\mu^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 $\theta_{41}$ and $\theta_{43}$ are defined analogously to the weak-mixing angle by

$$\cos \theta_{41} = \frac{g_4}{\sqrt{g_4^2 + \frac{2}{3} g_1^2}} = \frac{g_Y}{g_1}, \qquad \cos \theta_{43} = \frac{g_4}{\sqrt{g_4^2 + g_3^2}} = \frac{g_s}{g_3}. \qquad (2.5)$$

Since the couplings $g_s$ and $g_Y$ are the known QCD and hypercharge couplings, eq. (2.5) 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}, \quad g_4. \qquad (2.6)$$

The masses of $U_\mu$, $Z'_\mu$, and $G'_\mu$ are given by

$$\begin{align}
M_U^2 &= \frac{1}{4} g_4^2 v_{LQ}^2, \nonumber \\
M_{Z'}^2 &= \frac{1}{4} \left( g_4^2 + \frac{2}{3} g_1^2 \right) v_{LQ}^2, \nonumber \\
M_{G'}^2 &= \frac{1}{4} (g_4^2 + g_3^2) v_{LQ}^2, \tag{2.7}
\end{align}$$

<footnote n="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.</footnote>
<table><caption>**Table 1.** Quantum numbers of SM-like fermions (upper block) and the vector-like DM multiplet $\Psi_{DM}$ (last row). First and second generation fermions transform under $SU(3)' \times SU(2)_L \times U(1)_X$ like the usual SM fermions; the third generation quarks and leptons are unified into $\Psi_L'^3 = (q_L'^3 \ell_L'^3)^\dagger$, $\Psi_R'^{+3} = (u_R'^3 \nu_R'^3)^\dagger$, and $\Psi_R'^{-3} = (d_R'^3 e_R'^3)^\dagger$.</caption><thead><tr><th>Field</th><th>SU(4)</th><th>SU(3)'</th><th>SU(2)<sub>L</sub></th><th>U(1)<sub>X</sub></th></tr></thead><tbody><tr><td>ℓ'<sub>L</sub><sup>1,2</sup></td><td><b>1</b></td><td><b>1</b></td><td><b>2</b></td><td>-1/2</td></tr><tr><td>e'<sub>R</sub><sup>1,2</sup></td><td><b>1</b></td><td><b>1</b></td><td><b>1</b></td><td>-1</td></tr><tr><td>q'<sub>L</sub><sup>1,2</sup></td><td><b>1</b></td><td><b>3</b></td><td><b>2</b></td><td>+1/6</td></tr><tr><td>u'<sub>R</sub><sup>1,2</sup></td><td><b>1</b></td><td><b>3</b></td><td><b>1</b></td><td>+2/3</td></tr><tr><td>d'<sub>R</sub><sup>1,2</sup></td><td><b>1</b></td><td><b>3</b></td><td><b>1</b></td><td>-1/3</td></tr><tr><td>Ψ'<sub>L</sub><sup>3</sup></td><td><b>4</b></td><td><b>1</b></td><td><b>2</b></td><td>0</td></tr><tr><td>Ψ'<sub>R</sub><sup>+3</sup></td><td><b>4</b></td><td><b>1</b></td><td><b>1</b></td><td>+1/2</td></tr><tr><td>Ψ'<sub>R</sub><sup>-3</sup></td><td><b>4</b></td><td><b>1</b></td><td><b>1</b></td><td>-1/2</td></tr><tr><td>Ψ<sub>DM</sub></td><td><b>4</b></td><td><b>1</b></td><td><b>N</b></td><td>+1/2</td></tr></tbody></table>

and they are related through the angles $\theta_{41}$ and $\theta_{43}$ as<fnref n="5" />

$$M_U = M_{Z'} \cos \theta_{41} = M_{G'} \cos \theta_{43}, \qquad (2.8)$$

implying that, at tree level, the $U_1$ is expected to be the lightest new vector boson.

## 2.2 Fermions

Among the different possibilities for implementing the SM fermions in a 4321 model (see appendix A), a well-motivated and phenomenologically successful variant corresponds to a unification of third-family quarks and leptons [48–52]. In this case, the first and second families of SM-like fermions transform under the $SU(3)' \times SU(2)_L \times U(1)_X$ subgroup of the 4321 symmetry like the usual SM fermions, whereas the third-family quarks and leptons are unified into $\Psi_L'^3 = (q_L'^3 \ell_L'^3)^\dagger$, $\Psi_R'^{+3} = (u_R'^3 \nu_R'^3)^\dagger$, and $\Psi_R'^{-3} = (d_R'^3 e_R'^3)^\dagger$, 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 $\Psi_L'^3$ (cf. table 1) and we denote their left-handed components by $\Psi_L'^{1,2}$. 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

<footnote n="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].</footnote>
DM dynamics as long as their masses are larger than $M_\chi + M_\psi$, 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'_L^{1,2}$, $\ell'_L^{1,2}$ and the new heavy fields $\Psi'_L^{1,2}$. To avoid large flavor 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 flavor-diagonal for the fields

$$q^i = \begin{pmatrix} V_{ji}^* u_L^j \\ d_L^i \end{pmatrix}, \quad \ell^j = \begin{pmatrix} \nu_L^j \\ e_L^j \end{pmatrix}, \quad u_R^i, \quad d_R^i, \quad e_R^i, \quad \nu_R^i, \quad (2.9)$$

where $V$ is the CKM matrix and $u^i, d^i, e^i$, and $\nu^i$ are mass eigenstates. A possible misalignment between the quark and lepton components of the fields $\Psi'_L^i$ is parameterized by embedding the quark and lepton components $\Psi_{qL}^{1,2}$ and $\Psi_{\ell L}^{1,2}$ that have a flavor-diagonal mixing with $q'_L^{1,2}$ and $\ell'_L^{1,2}$, respectively, as

$$\Psi_L'^i = \begin{pmatrix} \Psi_{qL}^{i'} \\ W_{ij} \Psi_{\ell L}^{j'} \end{pmatrix}, \quad (2.10)$$

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

$$W = \begin{pmatrix} 1 & 0 & 0 \\ 0 & \cos \theta_{LQ} & \sin \theta_{LQ} \\ 0 & -\sin \theta_{LQ} & \cos \theta_{LQ} \end{pmatrix}. \quad (2.11)$$

In the absence of additional new heavy fermions that mix with right-handed SM fermions, a possible quark-lepton misalignment in $\Psi_R'^{+3}$ and $\Psi_R'^{-3}$ corresponds to only a phase difference, which we parameterize as

$$\Psi_R'^{+3} = \begin{pmatrix} \Psi_u^3 \\ e^{i\phi_\nu} \Psi_\nu^3 \end{pmatrix}, \quad \Psi_R'^{-3} = \begin{pmatrix} \Psi_d^3 \\ e^{i\phi_e} \Psi_e^3 \end{pmatrix}. \quad (2.12)$$

Consequently, the SM fields in the basis where the down-quark mass matrix is diagonal can be expressed as

$$\begin{align}
q_L^{1,2} &= q_L'^{1,2} \cos \theta_{q1,2} + \Psi_q'^{1,2} \sin \theta_{q1,2}, & q_L^3 &= \Psi_q'^3, \nonumber \\
\ell_L^{1,2} &= \ell_L'^{1,2} \cos \theta_{\ell1,2} + \Psi_{\ell L}'^{1,2} \sin \theta_{\ell1,2}, & \ell_L^3 &= \Psi_{\ell L}'^3, \nonumber \\
u_R^{1,2} &= u_R'^{1,2}, & u_R^3 &= \Psi_u'^3, \tag{2.13} \\
d_R^{1,2} &= d_R'^{1,2}, & d_R^3 &= \Psi_d'^3, \nonumber \\
e_R^{1,2} &= e_R'^{1,2}, & e_R^3 &= \Psi_e'^3, \nonumber \\
\nu_R^3 &= \Psi_\nu'^3. & \nonumber
\end{align}$$
<table><caption>**Table 2.** Constants $κ$ and $ξ$ entering the couplings of fermions to $G'$ and $Z'$ (cf. eq. (2.14)).</caption><thead><tr><th rowspan="2">Field</th><th colspan="2">κ</th><th colspan="2">ξ</th></tr><tr><th>i = 1, 2</th><th>i = 3</th><th>i = 1, 2</th><th>i = 3</th></tr></thead><tbody><tr><td>q<sup>i</sup><sub>L</sub></td><td>sin<sup>2</sup> θ<sub>q1,2</sub> - sin<sup>2</sup> θ<sub>43</sub></td><td>cos<sup>2</sup> θ<sub>43</sub></td><td>sin<sup>2</sup> θ<sub>q1,2</sub> - sin<sup>2</sup> θ<sub>41</sub></td><td>cos<sup>2</sup> θ<sub>41</sub></td></tr><tr><td>u<sup>i</sup><sub>R</sub></td><td>- sin<sup>2</sup> θ<sub>43</sub></td><td>cos<sup>2</sup> θ<sub>43</sub></td><td>-4 sin<sup>2</sup> θ<sub>41</sub></td><td>1 - 4 sin<sup>2</sup> θ<sub>41</sub></td></tr><tr><td>d<sup>i</sup><sub>R</sub></td><td>- sin<sup>2</sup> θ<sub>43</sub></td><td>cos<sup>2</sup> θ<sub>43</sub></td><td>2 sin<sup>2</sup> θ<sub>41</sub></td><td>1 + 2 sin<sup>2</sup> θ<sub>41</sub></td></tr><tr><td>ℓ<sup>i</sup><sub>L</sub></td><td></td><td></td><td>sin<sup>2</sup> θ<sub>ℓ1,2</sub> - sin<sup>2</sup> θ<sub>41</sub></td><td>cos<sup>2</sup> θ<sub>41</sub></td></tr><tr><td>e<sup>i</sup><sub>R</sub></td><td></td><td></td><td>-2 sin<sup>2</sup> θ<sub>41</sub></td><td>1 - 2 sin<sup>2</sup> θ<sub>41</sub></td></tr><tr><td>ν<sup>i</sup><sub>R</sub></td><td></td><td></td><td></td><td>1</td></tr><tr><td>ψ</td><td colspan="2">cos<sup>2</sup> θ<sub>43</sub></td><td colspan="2">1 - 4 sin<sup>2</sup> θ<sub>41</sub></td></tr><tr><td>χ</td><td colspan="2"></td><td colspan="2">1</td></tr></tbody></table>

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

$$ \begin{aligned} \mathcal{L}_{Z'} & \supset \frac{g_4}{2\sqrt{6}\cos\theta_{41}} Z'_{\mu} \left( \xi_q^i \bar{q}_L^i \gamma^{\mu} q_L^i + \xi_u^i \bar{u}_R^i \gamma^{\mu} u_R^i + \xi_d^i \bar{d}_R^i \gamma^{\mu} d_R^i + \xi_{\psi} \bar{\psi} \gamma^{\mu} \psi \right. \\ & \qquad \left. - 3 (\xi_{\ell}^i \bar{\ell}_L^i \gamma^{\mu} \ell_L^i + \xi_e^i \bar{e}_R^i \gamma^{\mu} e_R^i + \xi_{\nu} \bar{\nu}_R^i \gamma^{\mu} \nu_R^i + \xi_{\chi} \bar{\chi} \gamma^{\mu} \chi) \right), \\ \mathcal{L}_{G'} & \supset \frac{g_4}{\cos\theta_{43}} G'_{\mu}^a \left( \kappa_q^i \bar{q}^i \gamma^{\mu} T^a q^i + \kappa_u^i \bar{u}^i_R \gamma^{\mu} T^a u_R^i + \kappa_d^i \bar{d}^i_R \gamma^{\mu} T^a d_R^i + \kappa_{\psi} \bar{\psi} \gamma^{\mu} T^a \psi \right), \\ \mathcal{L}_{U_1} & \supset \frac{g_4}{\sqrt{2}} U_{\mu}^{+} \left( \beta_{q\ell}^{ij} \bar{q}^i \gamma^{\mu} \ell^j + \beta_{de} \bar{d}^3_R \gamma^{\mu} e_R^3 + \beta_{u\nu} \bar{u}^3_R \gamma^{\mu} \nu_R^3 + \bar{\psi} \gamma^{\mu} \chi \right) + \text{h.c.}, \end{aligned} \quad (2.14) $$

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

$$ \beta_{q\ell} = \begin{pmatrix} \sin\theta_{q1} & \sin\theta_{\ell1} & 0 & 0 \\ 0 & \sin\theta_{q2} & \cos\theta_{LQ} & \sin\theta_{q2} \\ 0 & -\sin\theta_{\ell2} & \sin\theta_{LQ} & \cos\theta_{LQ} \end{pmatrix}, \quad \beta_{de} = e^{i\phi_e}, \quad \beta_{u\nu} = e^{i\phi_\nu}. \quad (2.15) $$

Note that in the limit of large $g_4$, where $\cos\theta_{41} \approx \cos\theta_{43} \approx 1$ and $\sin\theta_{41} \approx \sin\theta_{43} \approx 0$, the constants $κ$ and $ξ$ are approximately

$$ \begin{aligned} \kappa_q^{1,2} & \approx \xi_q^{1,2} \approx \sin^2 \theta_{q1,2}, & \xi_{\ell}^{1,2} & \approx \sin^2 \theta_{\ell1,2}, & \kappa_u^{1,2} & = \kappa_d^{1,2} \approx \xi_u^{1,2} \approx \xi_d^{1,2} \approx \xi_e^{1,2} \approx 0, \\ \kappa_q^3 & = \kappa_u^3 = \kappa_d^3 \approx \xi_q^3 \approx \xi_u^3 \approx \xi_d^3 \approx \xi_{\ell}^3 \approx \xi_e^3 \approx \xi_{\nu}^3 \approx \kappa_{\psi} \approx \xi_{\chi} \approx \xi_{\psi} \approx 1, \end{aligned} \quad (2.16) $$

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 $κ \approx ξ \approx 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 flavor 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 $\beta_{de} = -1$ [51], which fixes the phase $\phi_e = \pi$.

* Since the phase $\phi_\nu$ is currently not constrained by any measurement one can use $\phi_\nu = 0$ for simplicity.

* An approximate U(2) symmetry in the quark sector, i.e. $\theta_{q_1} \approx \theta_{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 $\Delta C = 2$ observables would be possible.

* The first-generation lepton doublet can be taken to be purely a singlet of SU(4), i.e. $\theta_{\ell_1} = 0$, to be safe from LFV due to $U_1$ couplings involving the electron.

Making all of the above assumptions and defining $\theta_{q_{12}} = \theta_{q_1} = \theta_{q_2}$, the only remaining free parameters in the fermion sector are

$$ \theta_{q_{12}}, \theta_{\ell_2}, \theta_{LQ}, M_{\chi}, N, \qquad (2.17) $$

where $M_\chi$ denotes, here and henceforth, the mass of the DM-candidate$^6$ $\chi$ and $N$ is the dimension of its SU(2)$_L$ representation.

## 2.3 Parameter ranges

As summarized in eqs. (2.6) and (2.17), within reasonable assumptions the ‘effective’ model parameters are the following:

$$ g_4, v_{LQ}, \theta_{q_{12}}, \theta_{\ell_2}, \theta_{LQ}, M_{\chi}, N. \qquad (2.18) $$

In this section we would like to collect the non-negligible information available on these parameters from a fit to flavor data as well as from constraints due to direct searches. In section 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:

* $g_4, v_{LQ}, M_\chi, N$: the DM phenomenology depends crucially on these parameters.

* $\theta_{q_{12}}$: this parameter is important only for DM direct detection.

* $\theta_{\ell_2}, \theta_{LQ}$: the DM phenomenology is essentially independent of these parameters.

<footnote>$^6$The mass $M_\chi$ is related to the tree-level mass of the multiplet $\Psi_{DM}$ as discussed in appendix B.</footnote>
<figure><img src="image_2.png" /><figcaption>Figure 1. First-generation coupling combinations relevant for $Z'$ searches through the vertex $Z'\bar{u}u$ (left), $Z'\bar{d}d$ (center), or for $G'$ searches through the vertices $G'\bar{u}u$ or $G'\bar{d}d$ (right).</figcaption></figure>

Within our model, a combination of the parameters $v_{LQ}$, $\theta_{q_{12}}$, $\theta_{\ell_2}$, and $\theta_{LQ}$ is constrained by B-physics data alone, in particular by the $R(D^{(*)})$ discrepancies. A global fit of $v_{LQ}$ and the constants $\beta$ (cf. eq. (2.15)) was performed in ref. [51]. Expressed in our notation, the fit prefers values for $v_{LQ}$ in the range $v_{LQ}/\cos\theta_{LQ} \in [3.1, 4.6]$ TeV, with $\cos\theta_{LQ} \approx 0.8-0.9$. While the preferred value of $v_{LQ}$ is correlated with the preferred values of $\theta_{\ell_2}$ and $\theta_{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 $\theta_{\ell_2}$ and $\theta_{LQ}$ to comply with the fit in [51], while fulfilling all DM constraints. The fit to B-physics data leaves some freedom for $\theta_{q_{12}}$, which, however, is constrained by direct searches (see below). In short, for definiteness we take

$$v_{LQ} \in [3, 5] \text{ TeV} \qquad (2.19)$$

as our fiducial range for $v_{LQ}$.

The parameters $g_4$ and $\theta_{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 figure 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\theta_{q_{12}}$. This dependence displays transparently the $g_4$ and $\sin\theta_{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\theta_{q_{12}}$ and large $g_4$. Representative ranges are

$$\sin\theta_{q_{12}} \lesssim 0.2 \quad \text{and} \quad g_4 \gtrsim 3. \qquad (2.20)$$

These two requirements suppress respectively the two terms entering the coupling constants $\kappa_q^{1,2}$, $\xi_q^{1,2}$, and $\xi_{\ell}^{1,2}$ (cf. table 2).$^7$

Bounds from direct searches have been extensively studied in the literature, see [29, 44, 45, 47, 51, 60-64]. It is straightforward to verify that the ranges in eqs. (2.19)-(2.20)

<footnote>$^7$The figure also shows that, for $G'$ couplings, a cancellation between these two terms can be engineered for $g_4 \approx 1.5-2$ and $\sin\theta_{q_{12}} \ge 0.5$. However, the resulting light-quark — $G'$ coupling is not nearly as suppressed as in the case of $\sin\theta_{q_{12}} \lesssim 0.2$ and $g_4 \gtrsim 3$.</footnote>
yield $U_1$, $Z'$ and $G'$ masses and couplings in accord with these bounds. In particular, for $g_4 = 3$ and $v_{LQ} \ge 3$ TeV we have $M_U \approx M_{Z'} \ge 4.5$ TeV and $M_{G'} \ge 4.8$ TeV, such that the bounds found in [63] can be comfortably satisfied.

We conclude that agreement with B-decay discrepancies and with the constraints coming from direct searches select the parameter regions in eqs. (2.19)–(2.20). In section 5 we will confront this parameter space with the constraints imposed by Dark Matter.

## 3 Dark-Matter relic abundance

In this section we shall address the question whether our setup, as introduced around eq. (1.2) and discussed in section 2, can accommodate the relic abundance of DM observed today, $\Omega_0 h^2$.

In addition to the DM candidate $\chi_0$, the DM sector of our model includes also several co-annihilation partners: all the other components, charged under weak isospin, of the $\chi$ and $\psi$ SU(2)$_L$ multiplets. As shown in classic work, even in the presence of co-annihilators an estimate of $\Omega_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 $\Omega_0 h^2$ is the determination of $\Omega 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 \equiv M_{\chi_0}/T$, and to define $x_f \equiv x|_{T=T_f}$. In the case — like ours — where co-annihilators are present, $x_f$ is determined iteratively from the relation [65]

$$x_f = \ln \frac{0.038 g_{\text{eff}} M_{\text{Pl}} M_{\chi_0} \langle \sigma_{\text{eff}} v \rangle}{g_*^{1/2} x_f^{1/2}}, \qquad (3.1)$$

where $M_{\text{Pl}} = 1.22 \times 10^{19}$ GeV, $g_*$ denotes the total number of effectively relativistic d.o.f. at freeze-out, $g_{\text{eff}}$ denotes the number of effective d.o.f. within the DM sector, and the thermally averaged annihilation cross-section $\langle \sigma_{\text{eff}} v \rangle$ 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 = \int_{x_f}^{\infty} \frac{\langle \sigma_{\text{eff}} v \rangle}{x^2} dx. \qquad (3.2)$$

The present-day DM abundance can then be estimated as<fnref n="8" />

$$\Omega_0 h^2 = \sqrt{\frac{45}{\pi} \frac{s_0}{\rho_c} \frac{1}{g_*^{1/2} M_{\text{Pl}} J}} \approx \frac{1.07 \times 10^9 \text{ GeV}^{-1}}{g_*^{1/2} M_{\text{Pl}} J}. \qquad (3.3)$$

The crucial ingredient in the determination of $\Omega_0 h^2$ is the calculation of the thermally averaged annihilation cross-section $\langle \sigma_{\text{eff}} v \rangle$, which enters in $\Omega_0 h^2$ through the post-freeze-out annihilation efficiency $J$. To this end, it is also necessary to determine $g_{\text{eff}}$, which

<footnote n="8">⁸One can derive this relation by using $H(T) = \sqrt{8\pi^3 g_* / 90} T^2 / M_{\text{Pl}}$, $s = 2\pi^2 g_* T^3 / 45$ and $\rho_c = 3H_0^2 / (8\pi G_N)$, with $H(T_0) \equiv H_0 = 100h$ km/ (s Mpc), as customary.</footnote>
enters in both $\langle \sigma_{\text{eff}} v \rangle$ and $x_f$. In turn, $g_{\text{eff}}$ and $\sigma_{\text{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_{\text{eff}}$, and therewith proceed to the estimate of the thermally averaged annihilation cross-section $\langle \sigma_{\text{eff}} v \rangle$.

## 3.1 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 non-trivially 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 B). Applying our generic result, eq. (B.10), to the EW gauge group, we can determine the relative mass difference

$$ \Delta_{\xi\eta} = \frac{M_{\xi} - M_{\eta}}{\hat{M}} \qquad (3.4) $$

between components $\xi$ and $\eta$ of a VL multiplet of hypercharge $Y$ and mass $\hat{M}$. We find

$$ \Delta_{\xi\eta}^{\text{EW}} = \frac{g^2}{16\pi^2} \left\{ ((Q_{\xi}-Y)^2 - (Q_{\eta}-Y)^2) \left[ f\left(\frac{M_W}{\hat{M}}\right) - f\left(\frac{M_Z}{\hat{M}}\right) \right] + s_W^2 (Q_{\xi}^2 - Q_{\eta}^2) f\left(\frac{M_Z}{\hat{M}}\right) \right\}, \qquad (3.5) $$

where $Q_{\xi}$ and $Q_{\eta}$ are the electric charges of $\xi$ and $\eta$, and $f(r)$ is a finite loop function given in eq. (B.6). 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 $\psi$ and $\chi$ $SU(2)_L$ multiplets of our DM sector is between $\mathcal{O}(10^{-3})$ and $\mathcal{O}(10^{-4})$ for $\hat{M} = \mathcal{O}(1 \text{ TeV})$.

Having at hand the generic result, eq. (B.10), it is straightforward to determine the relative mass splitting between our DM candidate $\chi$ and its colored co-annihilation partner $\psi$. This mass splitting is induced by the vector bosons associated with the $43(2)1$ symmetry
<figure><img src="image_3.png" /><figcaption>Figure 2. Mass splitting between $\psi$ and $\chi$ induced by the 4321 gauge bosons.</figcaption></figure>

breaking and is given by

$$ \Delta_{\psi\chi}^{4321} = \frac{g_4^2}{16\pi^2} \left\{ f\left(\frac{M_U}{\hat{M}}\right) + \frac{1}{3}(2\sin^2\theta_{41}+1) f\left(\frac{M_{Z'}}{\hat{M}}\right) + \frac{4}{3}(\sin^2\theta_{43}-1) f\left(\frac{M_{G'}}{\hat{M}}\right) \right\}. \quad (3.6) $$

The value of $\Delta_{\psi\chi}^{4321}$ is around 8–15% for the parameter region of interest (see figure 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 $\Delta^{\text{EW}}(T) \sim (gT/\hat{M})^2$ [58] within $SU(2)_L$ multiplets and $\Delta^{4321}(T) \sim (g_4 T/\hat{M})^2$ within $SU(4)$ multiplets. Recalling that $T_f/\hat{M} \simeq 3\%$, and that such ratio enters quadratically in $\Delta(T)$, we estimate $\Delta^{\text{EW}}(T)$ to lie between $\mathcal{O}(10^{-3})$ and $\mathcal{O}(10^{-4})$ and $\Delta^{4321}(T) \lesssim 1\%$. Since $\Delta^{4321}(T)$ is small compared to the splitting induced by eq. (3.6), we neglect it in the following. The size of $\Delta^{\text{EW}}(T)$ is similar to the splitting induced by eq. (3.5) and we estimate the combination of both contributions to lie between $\mathcal{O}(10^{-3})$ and $\mathcal{O}(10^{-4})$.

## 3.2 Effective degrees of freedom

In our case, the number of effective d.o.f. in the DM sector $g_{\text{eff}}$ is given by [65]

$$ g_{\text{eff}} = \sum_i \left( g_{\chi} (1 + \Delta_{\chi_i})^{3/2} \exp(-x \Delta_{\chi_i}) + g_{\psi} (1 + \Delta_{\psi_i})^{3/2} \exp(-x \Delta_{\psi_i}) \right), \quad (3.7) $$

where the index $i$ runs over the $N$ components of the $SU(2)_L$ multiplets $\chi$ and $\psi$. Besides $g_{\chi} = 4$ and $g_{\psi} = 12$ denote the internal (spin, color, ...) d.o.f. of the components of these multiplets.<fnref n="9" /> The relative mass splittings $\Delta_{\chi_i}$ and $\Delta_{\psi_i}$ are defined as

$$ \Delta_{\chi_i} = (M_{\chi_i} - M_{\chi_0})/M_{\chi_0}, \quad \Delta_{\psi_i} = (M_{\psi_i} - M_{\psi_0})/M_{\psi_0}. \quad (3.8) $$

<footnote n="9">⁹Since $\chi$, $\psi$ belong to complex VL representations of the gauge group, they are Dirac fermions.</footnote>
The relative mass differences within the $\chi$ and $\psi$ multiplets are only at the per mil level (see discussion below eq. (3.5)). We thus neglect them in the following and use a common mass and relative mass splitting for each multiplet,

$$M_{\chi_i} \rightarrow M_{\chi}, \quad M_{\psi_i} \rightarrow M_{\psi}, \quad \Delta_{\chi_i} \rightarrow 0, \quad \Delta_{\psi_i} \rightarrow \Delta_{\psi}, \qquad (3.9)$$

such that $M_{\psi} = (1 + \Delta_{\psi}) M_{\chi}$. Employing this approximation, the number of effective d.o.f. simplifies to

$$g_{\text{eff}} \approx N \left( g_{\chi} + g_{\psi} (1 + \Delta_{\psi})^{3/2} \exp(-x \Delta_{\psi}) \right). \qquad (3.10)$$

We see that, at the freeze-out temperature, $g_{\text{eff}}$ departs appreciably from $N(g_{\chi} + g_{\psi})$ unless $x_f \Delta_{\psi} \ll 1$, i.e. unless the $\psi$-$\chi$ mass splitting is much smaller than the freeze-out temperature. Combining eq. (3.1) with the value of the $\psi$-$\chi$ mass splitting discussed in section 3.1, $\Delta_{\psi} \approx 0.1$, we find

$$(1 + \Delta_{\psi})^{3/2} \exp(-x_f \Delta_{\psi}) \approx 0.06. \qquad (3.11)$$

Eqs. (3.10)-(3.11) provide an accurate determination of $g_{\text{eff}}$ within our setup.

### 3.3 Thermally averaged cross-section

The main dynamical quantity in eq. (3.1) is $\langle \sigma_{\text{eff}} v \rangle$, where [65]

$$\begin{aligned} \sigma_{\text{eff}} = \frac{1}{g_{\text{eff}}^2} \sum_{i,j} & \left( \sigma_{\chi_i \chi_j} g_{\chi}^2 (1 + \Delta_{\chi_i})^{3/2} (1 + \Delta_{\chi_j})^{3/2} e^{-x(\Delta_{\chi_i} + \Delta_{\chi_j})} \right. \\ & \left. + \sigma_{\psi_i \psi_j} g_{\psi}^2 (1 + \Delta_{\psi_i})^{3/2} (1 + \Delta_{\psi_j})^{3/2} e^{-x(\Delta_{\psi_i} + \Delta_{\psi_j})} \right. \\ & \left. + 2 \sigma_{\chi_i \psi_j} g_{\chi} g_{\psi} (1 + \Delta_{\chi_i})^{3/2} (1 + \Delta_{\psi_j})^{3/2} e^{-x(\Delta_{\chi_i} + \Delta_{\psi_j})} \right) \end{aligned} \qquad (3.12)$$

with e.g.

$$\sigma_{\chi_i \psi_j} \equiv \sigma(\chi_i \psi_j \rightarrow XX'), \qquad (3.13)$$

and the other cross-sections defined analogously. Here $X, X'$ denote any particles other than $\chi, \psi$. 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 $\psi$ and $\chi$ multiplets, i.e. employing the replacements in eq. (3.9), the effective cross-section simplifies to

$$\sigma_{\text{eff}} = \frac{1}{g_{\text{eff}}^2} \sum_{i,j} \left( \sigma_{\chi_i \chi_j} g_{\chi}^2 + 2 \sigma_{\chi_i \psi_j} g_{\chi} g_{\psi} (1 + \Delta_{\psi})^{3/2} e^{-x \Delta_{\psi}} + \sigma_{\psi_i \psi_j} g_{\psi}^2 (1 + \Delta_{\psi})^3 e^{-2x \Delta_{\psi}} \right). \qquad (3.14)$$

The cross-sections $\sigma_{\chi_i \chi_j}$, $\sigma_{\chi_i \psi_j}$, and $\sigma_{\psi_i \psi_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, $\sigma_{\chi_i\psi_j}$ and $\sigma_{\psi_i\psi_j}$ are multiplied by one and two powers, respectively, of a factor that is suppressed by the $\psi$-$\chi$ mass splitting, cf. eq. (3.11). Consequently, the $Z'$-mediated cross-section $\sigma_{\chi_i\chi_j}$ is larger than any other contribution to $\sigma_{\text{eff}}$ by one to two orders of magnitude. In view of the overall 10% uncertainty in our analytical estimate, we can therefore approximate

$$ \sigma_{\text{eff}} \approx \frac{1}{N} \sigma(\chi_0\chi_0 \rightarrow Z' \rightarrow XX), \qquad (3.15) $$

where we have used eqs. (3.10)–(3.11) and $\sigma_{\chi_i\chi_j} \approx \delta_{ij} \sigma(\chi_0\chi_0 \rightarrow Z' \rightarrow XX)$.

The thermal average $\langle\sigma_{\text{eff}}v\rangle$ can be determined as an expansion in powers of $1/x$, and this expansion can be related order by order to the expansion of $\sigma_{\text{eff}}$ around $s = 4M_{\chi}^2$ [68]. To this end, we follow the notation of [69]. Neglecting the masses of the annihilation products, we define

$$ \sigma_0(y) = 2\sqrt{y^2 - y} \sigma_{\text{eff}}(y) \qquad (3.16) $$

where we substituted $s$ by the dimensionless variable $y = s/(4M_{\chi}^2)$. The thermal average $\langle\sigma_{\text{eff}}v\rangle$ can then be expressed as [68, 69]

$$ \langle\sigma_{\text{eff}} v\rangle = \sigma_0(1) \left(1 + \sum_{k=1}^{\infty} \frac{c_k}{x^k}\right), \qquad (3.17) $$

where the first coefficients are

$$ c_1 = -3 + \frac{3}{2}\lambda_1, \quad c_2 = 6 - 3\lambda_1 + \frac{15}{8}\lambda_2, \quad c_3 = -\frac{5}{16}(30 - 15\lambda_1 + 3\lambda_2 - 7\lambda_3), \quad (3.18) $$

and we have defined

$$ \lambda_n = \frac{1}{\sigma_0(y)} \left. \frac{d^n \sigma_0(y)}{dy^n} \right|_{y=1} . \qquad (3.19) $$

The above relations allow thus to verify that, by including higher powers in the small-velocity expansion of $\sigma_0(y)$, higher powers in the small-temperature expansion of $\langle\sigma_{\text{eff}}v\rangle$ are smaller and smaller.

## 3.4 Present-day DM abundance

In order to obtain the present-day DM abundance, one convolutes the calculated $\langle\sigma_{\text{eff}}v\rangle$ in the post-freeze-out annihilation efficiency $J$ in eq. (3.2). Using the expansion in eq. (3.17) to order $1/x_f^2$, the present-day DM abundance in eq. (3.3) yields

$$ \Omega_0 h^2 \approx \frac{1.07 \times 10^9 \text{ GeV}^{-1}}{g_*^{1/2} M_{\text{Pl}}} \cdot \frac{x_f}{\sigma_0(1)} \cdot \frac{1}{1 + c_1/(2x_f) + c_2/(3x_f^2)}. \qquad (3.20) $$

We will perform a full numerical study in section 5, following the discussion of the constraints imposed by direct detection in section 4. Here we would like to make a few qualitative considerations around eq. (3.20). Using $\sigma_{\text{eff}}$ as in eq. (3.15), eq. (3.16) yields

$$ \sigma_0(y) = \frac{1}{128 \pi} \left( \frac{g_4}{\cos \theta_{41}} \right)^4 \frac{M_{\chi}^2 (2y^2 + y)}{(4y M_{\chi}^2 - M_{Z'}^2)^2} \frac{1}{N} f(\{\xi^i\}), \qquad (3.21) $$
where $N = 1, 3, 5, ...$ denotes the SU(2)$_L$ size of the $\chi$ and $\psi$ multiplets and for brevity we introduced the flavor function

$$f(\{\xi^i\}) = \sum_{i=1}^{3} (2|\xi_q^i|^2 + |\xi_u^i|^2 + |\xi_d^i|^2 + 3(2|\xi_l^i|^2 + |\xi_e^i|^2 + |\xi_\nu^i|^2)). \quad (3.22)$$

Plugging eq. (3.21) into eq. (3.20), and taking the representative value $g_4 = 3$, we find

$$\Omega_0 h^2 \approx 0.06 \frac{N}{f(\{\xi^i\})} \left(\frac{v_{LQ}}{5 \text{ TeV}}\right)^2 \left(\frac{v_{LQ}}{M_\chi}\right)^2. \quad (3.23)$$

A few remarks are in order. First, eq. (3.23) assumes $4M_\chi^2 \ll M_{Z'}^2$. For $v_{LQ} = 3$ TeV (5 TeV), this approximation implies an error $\lesssim 15\%$ ($\lesssim 40\%$) in eq. (3.23), keeping in mind that the preferred range for $M_\chi$ are $\lesssim 600$ GeV ($\lesssim 1.5$ TeV), see section 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. (3.20)). These terms induce corrections in the per mil ballpark for $g_4 = 3$ and the just mentioned mass ranges for $M_\chi$ and $v_{LQ}$. Besides, these corrections do not depend on the choice of any other of our model's parameters, e.g. $f(\{\xi^i\})$ and $N$, as such dependences cancel in the $\lambda_n$ ratios. The function $f(\{\xi^i\})$ is typically of O(10) in the region satisfying all constraints. For example, taking table 2 with $g_4 = 3$, $\sin\theta_{q12} = 0.2$,<fnref n="10" /> one has $f(\{\xi^i\}) \approx 16$.

The procedure outlined in this section, and leading to eq. (3.20), with $\sigma_{\text{eff}}$ including $U_1$, $Z'$, $G'$ and gluon contributions, will be used in section 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. (3.23), and the discussion around it, demonstrate that $\Omega_0 h^2$ of the order of the observed value can be obtained without effort, in compliance with all other constraints.

# 4 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–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 section 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

$$\mathcal{L}_{\chi q} = \frac{g_{Z'}^2}{12 M_{Z'}^2} \xi_{\chi} (\bar{\chi} \gamma^{\mu} \chi) (\xi_q^i \bar{q}^i \gamma_{\mu} q^i + \xi_u^i \bar{u}_R^i \gamma_{\mu} u_R^i + \xi_d^i \bar{d}_R^i \gamma_{\mu} d_R^i). \quad (4.1)$$

<footnote n="10">¹⁰This choice of parameters follows from the discussion in section 2.3.</footnote>
Furthermore, if we are able to neglect corrections due to the finite momentum transfer between the DM and the nucleons, we may parameterize the matrix elements between vector or axial-vector quark q currents and the external-state nucleons N as

$$ \begin{aligned} \langle N(p') | \bar{q} \gamma^{\mu} q | N(p) \rangle |_{\vec{p}=\vec{p}'=0} &= F_1^{q/N}(0) \bar{u}_N(p') \gamma^{\mu} u_N(p), \\ \langle N(p') | \bar{q} \gamma^{\mu} \gamma_5 q | N(p) \rangle |_{\vec{p}=\vec{p}'=0} &= F_A^{q/N}(0) \bar{u}_N(p') \gamma^{\mu} \gamma_5 u_N(p). \end{aligned} \quad (4.2) $$

For the two form factors at zero momentum transfer we follow conventions common in the literature. In particular, $F_1^{q/N}(0)$ counts the number of valence quarks q in the nucleon N, e.g. $F_1^{d/n}(0) = 2$. eqs. (4.1)–(4.2) yield the following spin-independent cross-section for elastic scattering between DM and a single nucleon N = p or n

$$ \sigma_{\text{SI}}^{N} = \frac{g_{Z'}^{4} \xi_{\chi}^{2} M_{N}^{2}}{144 \pi M_{Z'}^{4}} |C_{V}^{N}|^{2}, \quad (4.3) $$

where

$$ C_V^p = 2C_V^u + C_V^d, \quad C_V^n = C_V^u + 2C_V^d, \quad (4.4) $$

and

$$ C_V^u = \frac{\xi_q^1 + \xi_u^1}{2}, \quad C_V^d = \frac{\xi_q^1 + \xi_d^1}{2}. \quad (4.5) $$

Starting from eq. (4.3), 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. (4.3) with the replacement

$$ |C_V^N|^2 \rightarrow \frac{|ZC_V^p + (A-Z)C_V^n|^2}{A^2}. \quad (4.6) $$

The above procedure is crude in a number of ways, that have been amply discussed in the literature [74–91]. A first outstanding limitation is the fact that the matching scale for the interactions in eq. (4.1) is well above the effective scale for the hadronic matrix elements in eq. (4.2), hence renormalization-group effects are in general non-negligible, and the relativistic-operator basis of eq. (4.1) 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. (4.2) loses validity.

The impact of the above approximations may be explored using the public codes DirectDM [90–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. (4.3)–(4.6) and the numerical estimate obtained within the DirectDM and DMFormFactor codes. The analytic approach of eqs. (4.3)–(4.6) yields a prediction in agreement to $\lesssim 25\%$ with the numerical estimate, provided $M_\chi \lesssim 400$ GeV, and for any $M_\chi \gtrsim 200$ GeV the numerical prediction is lower with respect to the analytic result by a factor of $\lesssim 2$ for $M_\chi \lesssim 1.5$ TeV. This conclusion holds for any choice of $v_{LQ}$ within our fiducial range (see eq. (2.19)).

This comparison is in agreement with the expectation that, for $M_\chi$ light enough, the DM Compton wavelength — of order $v_0/c \times M_\chi \approx 10^{-3} M_\chi$, 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_\chi$ masses in the range required by the relic-density constraint, 300 GeV–1.5 TeV, use of the analytic prediction of eqs. (4.3)–(4.6) will produce direct-detection bounds that are somewhat stronger — yet quite realistic — than those produced with numerical codes. In section 5 we will compare this analytic prediction with the latest bound obtained by the Xenon1T experiment [72].

# 5 Results

The relic-density and direct-detection constraints discussed in the previous sections represent significant phenomenological input for our model. We summarized our parameter space in eq. (2.18) and discussed how B-decay discrepancies and collider constraints lead to the preferred ranges in eqs. (2.19)–(2.20). 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\theta_{q12}$ ranges in eq. (2.20) are also favored by direct-detection constraints, as illustrated in figure 3. Here we show the DM-nucleon cross-section $\sigma_{\text{SI}}^N$ discussed around eq. (4.3), as a function of $g_4$ for increasing values of $\sin\theta_{q12} \ge 0$. As discussed in section 4, experiments yield severe limits, as strong as $\sigma_{\text{SI}}^N < 10^{-45}\text{cm}^2$. As the figure shows, these limits can be comfortably satisfied with the choice $\sin\theta_{q12} \lesssim 0.2$ and $g_4 \gtrsim 3$.<fnref n="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\theta_{q12}$. 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} \text{ cm}^2$.

We next discuss the behavior of $\sin\theta_{q12}$, $g_4$ and $v_{LQ}$ as a function of the DM mass $M_\chi$, when the relic-density and direct-detection constraints are imposed. In figure 4 we show as colored rays the regions selected by the constraint $\Omega_0 h^2$ in the $M_\chi$ vs. $v_{LQ}$ plane.

<footnote n="11">¹¹As an alternative, one may advocate $\sin\theta_{q12} \ge 0.25$ and $g_4 \lesssim 2$, but the constraint would be satisfied in (or very close to) a fine-tuned region.</footnote>
<figure><img src="image_4.png" /><figcaption>Figure 3. The DM-nucleon cross-section (see eq. (4.3)) as a function of $g_4$ for different choices of the mixing angle $\sin\theta_{q12}$ and $v_{LQ} = 3$ TeV. See text for details.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4. DM constraints in the $M_\chi$ vs. $v_{LQ}$ plane. Colored rays fulfill the $\Omega_0 h^2$ constraint within ±15%. Grey regions, with $\sin\theta_{q12}$ fixed at the displayed value, are excluded by Xenon1T [72]. See text for details.</figcaption></figure>

The constraint is imposed with ±15% accuracy, corresponding to the error we attach to its calculation — see discussion in section 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. (2.19), plus the $\Omega_0 h^2$ constraint, allow to identify the following, indicative $M_\chi$ ranges

$$[260, 720] \quad (N = 1)$$

$$M_{\chi}[\text{GeV}] \in [450, 1190] \quad (N = 3). \qquad (5.1)$$

$$[570, 1460] \quad (N = 5)$$
<figure><img src="image_6.png" /><figcaption>Figure 5. Colored rays fulfill the $\Omega_0 h^2$ constraint within $\pm 15\%$. The gray region corresponds to an on-shell intermediate $U_1$ leptoquark (see text for details).</figcaption></figure>

Also shown in different shades of gray depending on $\sin\theta_{q_{12}} \in \{0.2, 0.25, 0.3\}$ are regions excluded by direct detection. It is clear that, for $\sin\theta_{q_{12}} = 0.2$ and $v_{LQ} \ge 3$ TeV (see eq. (2.19)), 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 $\Omega_0 h^2$ constraint on $g_4$, rather than on $v_{LQ}$. We show such dependence in the two panels of figure 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_\chi \ge M_U/2$, that we excluded for simplicity. In fact, as the $U_1$ becomes on-shell, $2 \to 3$ and $2 \to 4$ decay channels open up, whereas we restricted to $2 \to 2$ processes in the calculation of $\sigma_{\text{eff}}$.

The dependence of the $\Omega_0 h^2$ constraint on $M_\chi$ vs. $v_{LQ}$ or $g_4$ shown in figures 4 and 5 can be captured simultaneously in figure 6, that shows this constraint in the plane $v_{LQ}$ vs. $g_4$ for a few reference values of $M_\chi$ and of $N$, represented as isolines. The figure shows at a glance that within the fiducial $v_{LQ}$ range of eq. (2.19), 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, figure 7 displays the allowed parameter space in the $M_\chi$ vs. $\sin\theta_{q_{12}}$ plane, with $g_4 = 3$. The colored ‘curtains’ show the region excluded by direct detection, depending on the choice of $N$. Similarly as figure 3 and the ensuing discussion, this plot shows that direct detection tends to prefer small $\sin\theta_{q_{12}}$, the upper bound for a given $M_\chi$ 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_\chi$ ranges in eq. (5.1). These ranges are not reported in figure 7 to limit clutter.
<figure><img src="image_7.png" /><figcaption>**Figure 6.** The Ω₀h² in the v<sub>LQ</sub> vs. g₄ plane. Along the isolines, the Ω₀h² constraint is fulfilled for the value of M<sub>χ</sub> specified along the line itself.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 7.** DM constraints in the M<sub>χ</sub> vs. sin θ<sub>q12</sub> plane, with g<sub>4</sub> = 3. The colored ‘curtains’ denote regions excluded by Xenon1T [72].</figcaption></figure>

# 6 Comments on Dark-Matter indirect detection

Beside direct-detection searches, our DM candidate may be constrained by so-called indirect-detection signals, namely by cosmic rays produced in DM annihilations in the Milky Way and elsewhere. Out of the numerous channels probed by experiments, three are relevant and potentially constraining for our model, namely $\chi\chi \to \gamma\gamma$, $\tau^+\tau^-$ and $W^+W^-$ (for a recent review see [93]). As already noted below eq. (3.14), the $W^+W^-$ channel has a smaller cross-section than the $\tau^+\tau^-$ one. As the experimental bound is also weaker for the W flux [94], we focus hereafter on the two other channels.
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Left: velocity-averaged DM-annihilation cross section to $\tau\tau$ at present time. The grey area denotes the region excluded by HESS [94], whereas the blue band represents the prediction in our model. The dashed line is the value required for $\langle\sigma v\rangle_0$ to saturate the $\Omega h^2$ constraint. Right: W-$\chi$ contribution to $\sigma(\chi_0\chi_0 \rightarrow \gamma\gamma)$ for $N=3,5$ and corresponding bounds from HESS [95] and Fermi-LAT [96]. See text for further details.</figcaption></figure>

The present-day $\chi_0\chi_0 \rightarrow \tau\tau$ cross-section, velocity-averaged within the Milky Way, may be estimated using eqs. (3.15), (3.17) and the fact that $x$-dependent terms are completely negligible. In this respect, we will denote such cross-section simply as $\sigma v$, following a similar notation used throughout the literature. We obtain

$$ \sigma(\chi_0\chi_0 \rightarrow \tau\tau)v \approx N\sigma_0(1)|_{f(\{\xi^i\})=6|\xi_e^3|^2+3|\xi_e^3|^2}, \quad (6.1) $$

where the choice of $f(\{\xi^i\})$ specializes $\sigma_0(1)$ to the $\chi_0\chi_0 \rightarrow \tau\tau$ case.

This cross-section is displayed in figure 8 (left). The width of the prediction corresponds to the range of $v_{LQ}$ values compatible with the $\Omega_0 h^2$ constraint for a given $M_\chi$ and taking into account the accuracy of 15% we attach to the $\Omega_0 h^2$ calculation — see discussion in section 3. The dashed line shows the value required to obtain the full relic density, if this were the only annihilation channel. Our prediction is below the current HESS bound [94], but close to it.<fnref n="12" /> An improvement of the limit by a factor of a few would offer a valuable probe of our scenario. Actually, an interesting question would be how much more exposure in this channel would lower the bound at the level of the signal.

We next turn to the $\chi_0\chi_0 \rightarrow \gamma\gamma$ channel. Although $\chi_0$ is electrically neutral, this process occurs at one loop via $W-\chi$ or $U_1-\psi$ exchange. As regards the $W-\chi$ contribution alone, given that $M_\chi \gg M_W$, one may consider the well-known large-$M_\chi$ calculation, which

<footnote n="12">¹²The HESS bound we show is the one obtained with the Einasto profile [97], producing the most constraining bound among the different DM-distribution profiles considered by ref. [94].</footnote>
yields (cf. e.g. [58, 98])

$$ \sigma v \approx (N^2 - 1)^2 \frac{\pi \alpha_{\text{em}}^2 \alpha_2^2}{16 M_W^2} \approx (N^2 - 1)^2 \cdot 2 \times 10^{-29} \text{ cm}^3 \text{ s}^{-1}. \qquad (6.2) $$

The W-χ contribution vanishes in the $N=1$ case, but is expected to yield the dominant constraint in the $N=3, 5$ cases. This one-loop result may undergo enhancements due to non-perturbative effects that become important in the $v \ll c$ and $M_χ \gg M_W$ limit, as elucidated in [98–101]. This enhancement is a multiplicative factor $R$ with respect to the 1-loop perturbative result. Defining $\epsilon \equiv M_W/M_χ$ and $\beta \equiv v/c$ with $c$ the speed of light in the vacuum, $R$ may be studied [102] in the two-parameter space of the ratios $\alpha_2/\epsilon$ vs. $\alpha_2/\beta$, where $\alpha_2 = g_2^2/(4\pi)$ is the involved gauge coupling. With typical values $\alpha_2/\beta \sim 30$, $\alpha_2/\epsilon \sim 0.4$, one can see that the perturbative result in eq. (6.2) can undergo a non-perturbative correction $2 \lesssim R \lesssim 3$.

The cross-section in eq. (6.2) — including the $R$ enhancement factor just discussed — may be compared with the 95% CL experimental bound in refs. [95, 96]. This comparison is shown in figure 8 (right), where the displayed uncertainties on the W-χ predictions come from the enhancement factor $R$. Assuming that interference with $U_1 - \psi$ exchange diagrams does not significantly reduce the cross-section, already in the case $N=3$ we obtain a velocity-averaged cross-section in excess of $10^{-27} \text{cm}^3 \cdot \text{s}^{-1}$ for $M_χ$ in the range of interest to us (cf. eq. (5.1)). The figure also shows the observed exclusion lines from HESS [95] and from Fermi-LAT [96] in the respective $M_χ$ ranges, and assuming a NFW [103] distribution for DM. We see that at face value this constraint favors $N=1$, and strongly disfavors $N \ge 5$.

For $N=1$, a detailed analysis of the $U_1 - \psi$ amplitude would be required since this is the only contribution in this case. Such a calculation is also interesting for $N=3, 5$ in order to determine the interference terms. An accurate estimate of the $\chi_0\chi_0 \rightarrow \gamma\gamma$ cross-section in our model along these lines will be interesting to further test its different scenarios against data. Such comparison will also depend in an important way on the model assumed for the distribution of DM in the Milky Way — keeping in mind that different such models imply different ‘best’ regions of interest within the Celestial dataset at each given $M_χ$. We reserve such study to future work.

# 7 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.<fnref n="13" />

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) \times SU(3)' \times SU(2)_L \times 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

<footnote n="13">¹³Such common description has also been studied elsewhere in the literature, see [104–131].</footnote>
the breaking to the SM group, this multiplet gives rise to a DM candidate and a colored co-annihilation partner with mass $\approx 10\%$ larger.

We find that the parameter space selected by collider and B-physics data is also the one favored 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} \in [3, 5]$ TeV, large SU(4) gauge coupling $g_4 \approx 3$ and small light-quark mixing parameter $\sin\theta_{q12} \lesssim 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.

It is also interesting that the model may be constrained by indirect-detection signals, in particular by DM annihilations into two photons. A first look into this channel seems to show a preference for the smallest possible $SU(2)_L$ DM representation — the singlet. Clearly, a complete calculation, along with an accurate modeling of the astrophysical aspects of the problem, could offer a powerful and novel set of tests.

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 section 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 Fermions in 4321 models

The fermion sector of the model contains the SM fermions, the $\psi$ and $\chi$, 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_\psi + M_\chi$, which we assume in the following.

### A.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 $\Psi_{4,2}$ and $\Psi_{1,2}$ shown in table 3, while the SM $SU(2)_L$ singlets are linear combinations of $\Psi_{4,1}$ and $\Psi_{1,1}$. To avoid large flavor violating effects, we
<table><caption>**Table 3.** Quantum numbers of fermions in the model. Fields with the same quantum numbers as the SM fermions are denoted by $\Psi_{IJ}$, where $I \in \{4, 1\}$ and $J \in \{2, 1\}$ are the dimensions of the SU(4) and SU(2)$_L$ representations, respectively. $\Psi_{DM}$ is the multiplet containing $\chi$ and $\psi$.</caption>
  <thead>
    <tr>
      <th></th>
      <th>SU(4)</th>
      <th>SU(3)'</th>
      <th>SU(2)<sub>L</sub></th>
      <th>U(1)<sub>X</sub></th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <td>&Psi;<sub>4,2</sub></td>
      <td><b>4</b></td>
      <td><b>1</b></td>
      <td><b>2</b></td>
      <td>0</td>
    </tr>
    <tr>
      <td>&Psi;<sup>q</sup><sub>1,2</sub> &oplus; &Psi;<sup>l</sup><sub>1,2</sub></td>
      <td><b>1</b></td>
      <td><b>3</b> &oplus; <b>1</b></td>
      <td><b>2</b></td>
      <td>+<sup>1</sup>&frasl;<sub>6</sub> &oplus; -<sup>1</sup>&frasl;<sub>2</sub></td>
    </tr>
    <tr>
      <td>&Psi;<sup>&uarr;</sup><sub>4,1</sub> &oplus; &Psi;<sup>&darr;</sup><sub>4,1</sub></td>
      <td><b>4</b></td>
      <td><b>1</b></td>
      <td><b>1</b></td>
      <td>+<sup>1</sup>&frasl;<sub>2</sub> &oplus; -<sup>1</sup>&frasl;<sub>2</sub></td>
    </tr>
    <tr>
      <td>(&Psi;<sup>u</sup><sub>1,1</sub> &oplus; &Psi;<sup>d</sup><sub>1,1</sub>) &oplus; (&Psi;<sup>&nu;</sup><sub>1,1</sub> &oplus; &Psi;<sup>e</sup><sub>1,1</sub>)</td>
      <td><b>1</b></td>
      <td><b>3</b> &oplus; <b>1</b></td>
      <td><b>1</b></td>
      <td>(+<sup>2</sup>&frasl;<sub>3</sub> &oplus; -<sup>1</sup>&frasl;<sub>3</sub>) &oplus; (0 &oplus; -1)</td>
    </tr>
    <tr>
      <td>&Psi;<sub>DM</sub></td>
      <td><b>4</b></td>
      <td><b>1</b></td>
      <td><b>N</b></td>
      <td>+<sup>1</sup>&frasl;<sub>2</sub></td>
    </tr>
  </tbody>
</table>

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 flavor-diagonal for the fields

$$ q^i = \begin{pmatrix} V_{ji}^* u_L^j \\ d_L^i \end{pmatrix}, \quad \ell^j = \begin{pmatrix} \nu_L^j \\ e_L^j \end{pmatrix}, \quad u_R^i, \quad d_R^i, \quad e_R^i, \quad \nu_R^i, \quad \text{(A.1)} $$

where $V$ is the CKM matrix. A misalignment between the quark and lepton components of $\Psi_{4,2}$ is implemented by embedding the components $\Psi_{4,2}^q$ and $\Psi_{4,2}^l$ that have a flavor-diagonal mixing with $\Psi_{1,2}^q$ and $\Psi_{1,2}^l$, respectively, as

$$ \Psi_{4,2}^i = \begin{pmatrix} \tilde{\Psi}_{4,2}^q \end{pmatrix} \begin{pmatrix} \tilde{\Psi}_{4,2}^l \end{pmatrix} = \begin{pmatrix} \Psi_{4,2}^q \end{pmatrix} \begin{pmatrix} \Psi_{4,2}^l \end{pmatrix}, \quad (W_{ij} \Psi_{4,2}^l)^j, \quad \text{(A.2)} $$

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 ($\Psi_{4,1}^u$ and $\Psi_{4,1}^\nu$), and ($\Psi_{4,1}^d$ and $\Psi_{4,1}^e$) of $\Psi_{4,1}^\uparrow$, and $\Psi_{4,1}^\downarrow$, respectively, i.e.

$$ \Psi_{4,1}^{\uparrow i} = \begin{pmatrix} \Psi_{4,1}^{u i} \\ e^{i\phi_{\nu_i}} \Psi_{4,1}^{\nu i} \end{pmatrix}, \quad \Psi_{4,1}^{\downarrow i} = \begin{pmatrix} \Psi_{4,1}^{d i} \\ e^{i\phi_{e_i}} \Psi_{4,1}^{e i} \end{pmatrix}. \quad \text{(A.3)} $$

Consequently, the SM fields can be expressed as

$$ \begin{align} q^i &= \cos \theta_{q_i} (\Psi_{1,2}^q)^i + \sin \theta_{q_i} (\Psi_{4,2}^q)^i, \notag \\ \ell^i &= \cos \theta_{\ell_i} (\Psi_{1,2}^\ell)^i + \sin \theta_{\ell_i} (\Psi_{4,2}^\ell)^i, \notag \\ u_R^i &= \cos \theta_{u_i} (\Psi_{1,1}^u)^i + \sin \theta_{u_i} (\Psi_{4,1}^u)^i, \tag{A.4} \\ d_R^i &= \cos \theta_{d_i} (\Psi_{1,1}^d)^i + \sin \theta_{d_i} (\Psi_{4,1}^d)^i, \notag \\ e_R^i &= \cos \theta_{e_i} (\Psi_{1,1}^e)^i + \sin \theta_{e_i} (\Psi_{4,1}^e)^i, \notag \\ \nu_R^i &= \cos \theta_{\nu_i} (\Psi_{1,1}^\nu)^i + \sin \theta_{\nu_i} (\Psi_{4,1}^\nu)^i. \notag \end{align} $$
The couplings of the SM fermions to the new vector bosons are given by

$$ \begin{aligned} \mathcal{L}_{Z'} & \supset \frac{g_4}{2\sqrt{6} \cos \theta_{41}} Z'_{\mu} \left( \xi_q^i \bar{q}_L^i \gamma^{\mu} q_L^i + \xi_u^i \bar{u}_R^i \gamma^{\mu} u_R^i + \xi_d^i \bar{d}_R^i \gamma^{\mu} d_R^i \right. \\ & \qquad \left. - 3 \left( \xi_{\ell}^i \bar{\ell}_L^i \gamma^{\mu} \ell_L^i + \xi_e^i \bar{e}_R^i \gamma^{\mu} e_R^i + \xi_{\nu} \bar{\nu}_R^i \gamma^{\mu} \nu_R^i \right) \right), \end{aligned} \quad (A.5) $$
$$ \mathcal{L}_{G'} \supset \frac{g_4}{\cos \theta_{43}} G'_{\mu}^{a} \left( \kappa_q^i \bar{q}^i \gamma^{\mu} T^a q^i + \kappa_u^i \bar{u}_R^i \gamma^{\mu} T^a u_R^i + \kappa_d^i \bar{d}_R^i \gamma^{\mu} T^a d_R^i \right), $$
$$ \mathcal{L}_{U_1} \supset \frac{g_4}{\sqrt{2}} U_{\mu}^{+} \left( \beta_{q\ell}^{ij} \bar{q}^i \gamma^{\mu} \ell^j + \beta_{de}^i \bar{d}_R^i \gamma^{\mu} e_R^i + \beta_{u\nu}^i \bar{u}_R^i \gamma^{\mu} \nu_R^i \right) + \text{h.c.}, $$

where

$$ \begin{aligned} \kappa_q^i &= \sin^2 \theta_{q_i} - \sin^2 \theta_{43}, & \kappa_d^i &= \sin^2 \theta_{d_i} - \sin^2 \theta_{43}, & \kappa_u^i &= \sin^2 \theta_{u_i} - \sin^2 \theta_{43}, \\ \xi_q^i &= \sin^2 \theta_{q_i} - \sin^2 \theta_{41}, & \xi_d^i &= \sin^2 \theta_{d_i} + 2 \sin^2 \theta_{41}, & \xi_u^i &= \sin^2 \theta_{u_i} - 4 \sin^2 \theta_{41}, \\ \xi_{\ell}^i &= \sin^2 \theta_{\ell_i} - \sin^2 \theta_{41}, & \xi_e^i &= \sin^2 \theta_{e_i} - 2 \sin^2 \theta_{41}, & \xi_{\nu}^i &= \sin^2 \theta_{\nu_i}, \\ \beta_{q\ell}^{ij} &= \sin \theta_{q_i} W_{ij} \sin \theta_{\ell_j}, & \beta_{de}^i &= \sin \theta_{d_i} \sin \theta_{e_i} e^{i\phi_{e_i}}, & \beta_{u\nu}^i &= \sin \theta_{u_i} \sin \theta_{\nu_i} e^{i\phi_{\nu_i}}. \end{aligned} \quad (A.6) $$

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 right-handed SM fermions are purely singlets of SU(4). This corresponds to the choice

$$ \sin \theta_{u_i} = \sin \theta_{d_i} = \sin \theta_{e_i} = \sin \theta_{\nu_i} = 0, \quad (A.7) $$

The misalignment matrix W is usually chosen to be CP-conserving and to mix only the second and third generation, i.e.

$$ W = \begin{pmatrix} 1 & 0 & 0 \\ 0 & \cos \theta_{LQ} & \sin \theta_{LQ} \\ 0 & -\sin \theta_{LQ} & \cos \theta_{LQ} \end{pmatrix}. \quad (A.8) $$

Consequently, the only free parameters in the fermion sector are

$$ \theta_{q_1}, \theta_{q_2}, \theta_{q_3}, \theta_{\ell_1}, \theta_{\ell_2}, \theta_{\ell_3}, \theta_{LQ}. \quad (A.9) $$

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. $\theta_{q_1} = \theta_{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 $\Delta C = 2$ observables would be possible.

- The first-generation lepton doublet can be taken to be purely a singlet of SU(4), i.e. $\theta_{\ell_1} = 0$, to be safe from LFV due to $U_1$ couplings involving the electron.
Making both of the two above assumptions and defining $\theta_{q_{12}} = \theta_{q_1} = \theta_{q_2}$, the only free parameters in the fermion sector are

$$ \theta_{q_{12}}, \theta_{q_3}, \theta_{\ell_2}, \theta_{\ell_3}, \theta_{LQ}. \qquad (A.10) $$

If one further maximizes third generation couplings<fnref n="14" /> by taking $\theta_{q_3} = \theta_{\ell_3} = \frac{\pi}{2}$, the set of parameters further reduces to

$$ \theta_{q_{12}}, \theta_{\ell_2}, \theta_{LQ}. \qquad (A.11) $$

* Flavored 4321 models. In “flavored” 4321 models [48–51], all third generation SM fermions are fully unified into **4** representations of SU(4) and only the left-handed first- and second-generation fermions are each a mixture of a **4** and a **1** of SU(4). This corresponds to the choice

$$ \begin{aligned} \sin \theta_{q_3} &= \sin \theta_{\ell_3} = \sin \theta_{u_3} = \sin \theta_{d_3} = \sin \theta_{e_3} = \sin \theta_{\nu_3} = 1, \\ \sin \theta_{u_1} &= \sin \theta_{u_2} = \sin \theta_{d_1} = \sin \theta_{d_2} = \sin \theta_{e_1} = \sin \theta_{e_2} = \sin \theta_{\nu_1} = \sin \theta_{\nu_2} = 0, \quad (A.12) \\ e^{i\phi_{e_3}} &= -1. \end{aligned} $$

The misalignment matrix W is usually chosen as in eq. (A.8). Consequently, the only free parameters in the fermion sector are

$$ \theta_{q_1}, \theta_{q_2}, \theta_{\ell_1}, \theta_{\ell_2}, \theta_{LQ}. \qquad (A.13) $$

Making the above described assumptions to reduce contributions to $\Delta C = 2$ observables and LFV electron couplings, the set of free parameters in the fermion sector is reduced to

$$ \theta_{q_{12}}, \theta_{\ell_2}, \theta_{LQ}. \qquad (A.14) $$

## A.2 The fermions in the DM sector

We consider a DM candidate $\chi$ that, together with its coannihilation partner $\psi$, is part of a vector-like **4** of SU(4) denoted by $\Psi_{DM}$ (cf. table 3). The couplings of $\chi$ and $\psi$ to the new vector bosons and the gluons are thus given by

$$ \begin{aligned} \mathcal{L}_{Z'} &\supset \frac{g_4}{2\sqrt{6} \cos \theta_{41}} Z'_{\mu} (\xi_{\psi} \bar{\psi} \gamma^{\mu} \psi - 3 \xi_{\chi} \bar{\chi} \gamma^{\mu} \chi), \\ \mathcal{L}_{G'} &\supset \frac{g_4}{\cos \theta_{43}} \kappa_{\psi} G'_{\mu}^{a} \bar{\psi} \gamma^{\mu} T^a \psi, \qquad (A.15) \\ \mathcal{L}_{G} &\supset g_s G_{\mu}^{a} \bar{\psi} \gamma^{\mu} T^a \psi, \\ \mathcal{L}_{U_1} &\supset \frac{g_4}{\sqrt{2}} U_{\mu}^{+} \bar{\psi} \gamma^{\mu} \chi + \text{h.c.}, \end{aligned} $$

<footnote n="14">¹⁴Maximizing only the left-handed third-generation couplings, i.e. unifying the third-generation quark and lepton doublets in a pure **4** of SU(4) while keeping the right-handed third-generation fermions pure singlets of SU(4) might be problematic for generating the large Higgs Yukawa coupling in the third generation. In such a case, a “flavored 4321” as described below might be preferable.</footnote>
where

$$ \xi_{\psi} = 1 - 4 \sin^2 \theta_{41}, \quad \xi_{\chi} = 1, \quad \kappa_{\psi} = \cos^2 \theta_{43}. \qquad (A.16) $$

While the above couplings are independent of the representation $N$ of SU(2)$_L$ under which $\Psi_{DM}$ transforms, the couplings to the $W$ and $Z$ bosons are clearly different for different representations. The coupling of $Z$ to a field $\Psi_N$ transforming as a $N$ of SU(2)$_L$ and having hypercharge $Y$ is given by

$$ \mathcal{L}_Z \supset \frac{g}{\cos \theta_W} \bar{\Psi}_N (T_N^3 - \sin^2 \theta_W Q) \gamma^\mu \Psi_N Z_\mu, \qquad (A.17) $$

where $T_N^3$ is the diagonal generator of SU(2) in the $N$ representation and the electric charge is defined by $Q = T_N^3 + Y$. The coupling of $W^{\pm}$ to a field $\Psi_N$ transforming as a $N$ of SU(2)$_L$ is derived from the covariant derivative

$$ \begin{aligned} \bar{\Psi}_N i D_\mu \gamma^\mu \Psi_N &\supset \bar{\Psi}_N (i \partial_\mu + g W_\mu^a T_N^a) \gamma^\mu \Psi_N \\ &\supset \frac{g}{\sqrt{2}} \bar{\Psi}_N (W_\mu^+ T_N^+ + W_\mu^- T_N^-) \gamma^\mu \Psi_N \qquad (A.18) \\ &= \frac{g}{\sqrt{2}} W_\mu^+ \bar{\Psi}_N \gamma^\mu T_N^+ \Psi_N + \text{h.c.} \end{aligned} $$

where $T_N^{\pm} = T_N^1 \pm i T_N^2$.

## B 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 \times H'$ that is spontaneously broken to its subgroup $H$,

$$ G \times H' \to H, \qquad (B.1) $$

where $G$ has a subgroup $H^G \subseteq G$ that is isomorphic to both $H'$ and the unbroken $H$, i.e. $H' \cong H \cong H^G$, and $H$ is the diagonal subgroup of $H^G \times 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 \times H'_2 \times \cdots \times H'_n$.15 An analogous decomposition into simple groups $H_i$ and $H_i^G$ 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 $\theta_i$ by

$$ c_i \equiv \cos \theta_i = \frac{g_{H_i}}{g_{H'_i}}, \quad s_i \equiv \sin \theta_i = k_i \frac{g_{H_i}}{g_G}, \qquad (B.2) $$

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. For non-abelian $H_i$, we set $k_i = 1$. After the spontaneous symmetry breaking $G \times H' \to H$, there are dim($H$) massless vector bosons

<footnote n="15">¹⁵In case of a semi-simple $G = G_1 \times G_2 \times \cdots \times G_m$, one can treat each simple group factor $G_i$ separately. The results given here for a simple group $G$ can therefore be generalized easily to a semi-simple $G$.</footnote>
that are linear combinations of the $H'$ and $H^G$ gauge bosons. The orthogonal linear combinations constitute $\dim(H)$ massive vector bosons with masses $M_{H_i}$. In addition, all vector bosons associated with the coset $G/H$ become massive and have the common mass $M_{G/H}$.

We find it illustrative to show how the above described spontaneously broken gauge group is a generalization both of the EW group in the SM and of the 421 part in 4321 models.

* For EW symmetry breaking, $G = \text{SU}(2)_L$, $H' = \text{U}(1)_Y$, and $H = \text{U}(1)_{\text{em}}$. The couplings are $g_G = g$, $g_{H'} = g_Y$, and $g_H = e$. In this case, there is only a single mixing angle that can be identified with the weak-mixing angle, $\theta = \theta_W$. There is $\dim(H) = 1$ massless vector boson, which can be identified with the photon, $\dim(H) = 1$ massive vector boson, which can be identified with the $Z$, and $\dim(G/H) = 2$ massive coset vector bosons that can be identified with $W^{\pm}$. If one defines the electric charge $Q$ as $Q = T_3+Y$, where $T_3$ is the diagonal generator of $\text{SU}(2)_L$ normalized as $\text{tr}[T_3T_3] = 1/2$ for the fundamental representation, then $k=1$.

* For the 4321 breaking, $G = \text{SU}(4)$, $H'_1 = \text{U}(1)_X$, $H'_2 = \text{SU}(3)'$, $H_1 = \text{U}(1)_Y$, and $H_2 = \text{SU}(3)_c$. The couplings are $g_G = g_4$, $g_{H'_1} = g_1$, $g_{H'_2} = g_3$, $g_{H_1} = g_Y$, and $g_{H_2} = g_s$. In this case, there are two mixing angles $\theta_1 = \theta_{41}$ and $\theta_2 = \theta_{43}$. There are $\dim(H_1) = 1$ plus $\dim(H_2) = 8$ massless vector bosons, which can be identified with the $B$ and the gluons, $\dim(H_1) = 1$ plus $\dim(H_2) = 8$ massive vector bosons, which can be identified with $Z'$ and $G'$, and $\dim(G/H) = 6$ massive coset vector bosons that can be identified with colored $U^{\pm}$. In order to be able to use the conventional normalization of the electric charge defined by $Q = T_3+Y$, one has to set $k_1 = \sqrt{2/3}$.

We now consider a VL fermion multiplet that transforms under a representation $R_G$ of $G$ and under representations $R_{H'_i}$ of the $H'_i$. Note that since the fermion multiplet is VL under $G \times H'$, all of its components transform under $G$ and $H'$ in the same way. However, the components transform differently under the subgroup $H^G$ and under the unbroken group $H$. A given multiplet component $\xi$ transforms under representations $r_{H_i^G}^\xi$ of the $H_i^G$ and under representations $r_{H'_i}^\xi$ of the $H'_i$. The one-loop pole mass can then be expressed in terms of the quadratic Casimir invariants of the different groups and representations involved. In particular, we denote them as follows:

* $C_2^G(R_G)$: quadratic Casimir of the $R_G$ representation of $G$,

* $C_2^{H'_i}(R_{H'_i})$: quadratic Casimir of the $R_{H'_i}$ representation of $H'_i$,

* $C_2^{H_i^G}(r_{H_i^G}^\xi)$: quadratic Casimir of the $r_{H_i^G}^\xi$ representation of $H_i^G$,

* $C_2^{H_i}(r_{H_i}^\xi)$: quadratic Casimir of the $r_{H_i}^\xi$ representation of $H_i$.

Each multiplet component $\xi$ receives contributions to its pole mass from three different kinds of one-loop diagrams.
1. From the massive vector bosons that transform in the adjoint representation of H and have masses $M_{H_i}$,

$$ \Sigma_H = \frac{g_G^2 \hat{M}}{16 \pi^2} \sum_i \left( C_2^{H_i}(r_{H_i}^\xi) + \frac{s_i^2}{c_i^2 k_i^2} C_2^{H'_i}(R_{H'_i}) - \frac{s_i^2}{k_i^2} C_2^{H_i}(r_{H_i}^\xi) \right) \left[ f\left(\frac{M_{H_i}}{\hat{M}}\right) + K(\hat{M}, \mu) \right] \quad (B.3) $$

2. From the massless vector bosons that transform in the adjoint representation of H,

$$ \Sigma_0 = \frac{g_G^2 \hat{M}}{16 \pi^2} \sum_i \frac{s_i^2}{k_i^2} C_2^{H_i}(r_{H_i}^\xi) K(\hat{M}, \mu) \quad (B.4) $$

3. From the massive vector bosons that correspond to the $G/H$ coset and have mass $M_{G/H}$,

$$ \Sigma_{G/H} = \frac{g_G^2 \hat{M}}{16 \pi^2} \left( C_2^G(R_G) - \sum_i C_2^{H_i}(r_{H_i}^\xi) \right) \left[ f\left(\frac{M_{G/H}}{\hat{M}}\right) + K(\hat{M}, \mu) \right] \quad (B.5) $$

In the above expressions, $K(\hat{M}, \mu)$ is a divergent term that depends only on the VL fermion mass and on the renormalization scale $\mu$, while $f(r)$ is a finite loop function given by (cf. [58])

$$ f(r) = r^4 \ln r - r^2 + \frac{r}{2} \sqrt{r^2 - 4} (r^2 + 2) \ln \left( \frac{r^2}{2} - \frac{r}{2} \sqrt{r^2 - 4} - 1 \right) . \quad (B.6) $$

The $\xi$ pole mass $M_\xi$ can then be written as

$$ M_{\xi} = \hat{M} - \Sigma_{H} - \Sigma_{0} - \Sigma_{G/H}. \quad (B.7) $$

It is interesting to note that the contributions from the massless vector bosons, eq. (B.4), exactly cancel the divergent and scale-dependent terms in eq. (B.3) that are proportional to $C_2^{H_i}(r_{H_i}^\xi)$. Similarly, the divergent and scale-dependent terms proportional to $C_2^{H'_i}(r_{H'_i}^\xi)$ cancel between eq. (B.3) and eq. (B.5). All remaining divergent and scale-dependent terms are either proportional to $C_2^G(R_G)$ or $C_2^{H'_i}(R_{H'_i})$, i.e. these terms are the same for any multiplet component and therefore cancel in the mass differences. Consequently, the $\xi$ pole mass can be written as

$$ M_{\xi} = \hat{M} - I - \frac{g_G^2 \hat{M}}{16 \pi^2} \sum_i \left\{ C_2^{H_i}(r_{H_i}^\xi) \left[ f\left(\frac{M_{H_i}}{\hat{M}}\right) - f\left(\frac{M_{G/H}}{\hat{M}}\right) \right] - \frac{s_i^2}{k_i^2} C_2^{H_i}(r_{H_i}^\xi) f\left(\frac{M_{H_i}}{\hat{M}}\right) \right\}, \quad (B.8) $$

where $I$ collects all terms that are the same for each multiplet component and is given by

$$ I = \frac{g_G^2 \hat{M}}{16 \pi^2} \left\{ C_2^G(R_G) \left[ f\left(\frac{M_{G/H}}{\hat{M}}\right) + K(\hat{M}, \mu) \right] + \sum_i \frac{s_i^2}{c_i^2 k_i^2} C_2^{H'_i}(R_{H'_i}) \left[ f\left(\frac{M_{H'_i}}{\hat{M}}\right) + K(\hat{M}, \mu) \right] \right\}. \quad (B.9) $$
Using the result for the one-loop pole mass, we find a generic expression for the relative mass splitting between the components $\xi$ and $\eta$ of a VL fermion multiplet,

$$ \Delta_{\xi\eta} = \frac{M_{\xi}-M_{\eta}}{\hat{M}} = \frac{g_G^2}{16\pi^2} \sum_i \left\{ \left( C_2^{H_i}(r_{H_i}^{\xi}) - C_2^{H_i}(r_{H_i}^{\eta}) \right) \left[ f\left(\frac{M_{G/H}}{\hat{M}}\right) - f\left(\frac{M_{H_i}}{\hat{M}}\right) \right] + \frac{s_i^2}{k_i^2} \left( C_2^{H_i}(r_{H_i}^{\xi}) - C_2^{H_i}(r_{H_i}^{\eta}) \right) f\left(\frac{M_{H_i}}{\hat{M}}\right) \right\}, \quad (B.10) $$

which is finite and scale-independent at one-loop.

Given the generic result, it is straightforward to consider the special cases of the EW gauge group and of the 4321 models.

* 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

$$ C_2^{H_i} = T_3 T_3 = (Q - Y)^2, \quad C_2^H = Q^2, \quad (B.11) $$

and we find

$$ \Delta_{\xi\eta} = \frac{g^2}{16\pi^2} \left\{ \left( (Q_{\xi} - Y)^2 - (Q_{\eta} - Y)^2 \right) \left[ f\left(\frac{M_W}{\hat{M}}\right) - f\left(\frac{M_Z}{\hat{M}}\right) \right] + s_W^2 (Q_{\xi}^2 - Q_{\eta}^2) f\left(\frac{M_Z}{\hat{M}}\right) \right\}, \quad (B.12) $$

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 non-abelian factor. The quadratic Casimir invariants are

$$ \begin{gathered} C_2^{H_1} = T^{15} T^{15} = \frac{3}{2}(Y - X)^2, \quad C_2^{H_1} = Y^2, \quad (B.13) \\ C_2^{H_2} = C_2^{SU(3)_c}, \quad C_2^{H_2} = C_2^{SU(3)_c}, \quad (B.14) \end{gathered} $$

i.e. the Casimir invariants of the abelian groups can be expressed by the U(1) charges Y and X, while those of the non-abelian factors are SU(3) Casimir invariants.

For our DM candidate $\chi$ and its co-annihilation partner $\psi$, the Casimir invariants are given by

$$ \begin{gathered} C_2^{H_1}(\chi) = \frac{3}{8}, \quad C_2^{H_1}(\psi) = 0, \quad C_2^{H_2}(\chi) = 0, \quad C_2^{H_2}(\psi) = 0, \quad (B.14) \\ C_2^{H_1}(\psi) = \frac{1}{24}, \quad C_2^{H_1}(\chi) = \frac{4}{9}, \quad C_2^{H_2}(\psi) = \frac{4}{3}, \quad C_2^{H_2}(\chi) = \frac{4}{3}, \end{gathered} $$

such that the $\psi - \chi$ mass splitting is

$$ \Delta_{\psi\chi} = \frac{g_4^2}{16\pi^2} \left\{ \left( \frac{1}{24} - \frac{3}{8} \right) \left[ f\left(\frac{M_U}{\hat{M}}\right) - f\left(\frac{M_{Z'}}{\hat{M}}\right) \right] + \frac{3}{2} \frac{s_{41}^2}{9} f\left(\frac{M_{Z'}}{\hat{M}}\right) \right. \\ \left. + \frac{4}{3} \left[ f\left(\frac{M_U}{\hat{M}}\right) - f\left(\frac{M_{G'}}{\hat{M}}\right) \right] + s_{43}^2 \frac{4}{3} f\left(\frac{M_{G'}}{\hat{M}}\right) \right\}, \quad (B.15) $$

which can be simplified to

$$ \Delta_{\psi\chi} = \frac{g_4^2}{16\pi^2} \left\{ f\left(\frac{M_U}{\hat{M}}\right) + \frac{1}{3}(2s_{41}^2+1)f\left(\frac{M_{Z'}}{\hat{M}}\right) + \frac{4}{3}(s_{43}^2-1)f\left(\frac{M_{G'}}{\hat{M}}\right) \right\}. \quad (B.16) $$
C Cross-sections of processes entering the estimation of the relic density

Due to the dependence of the cross-sections on the fourth power of the couplings, contri-
butions to the effective cross-section eq. (3.14) boils down to those involving new heavy
gauge bosons and gluons. The contribution of the latter being additionally suppressed by
the mass of the DM candidate, we found

$$
\begin{align}
\sigma^{\chi\chi} &\approx \frac{g_4^4}{512 \pi c_{41}^4} \frac{s(2 M_{\chi}^2 + s)}{\sqrt{s^2 - 4s M_{\chi}^2 (s - M_{Z'}^2)^2}} \tag{C.1} \\
&\qquad \times \sum_{i=1}^3 (2|\xi_q^i|^2 + |\xi_u^i|^2 + |\xi_d^i|^2 + 3(2|\xi_\ell^i|^2 + |\xi_e^i|^2 + |\xi_\nu^i|^2)), \nonumber
\end{align}
$$

$$
\begin{equation}
\begin{split}
\sigma^{\chi\psi} &\approx \frac{g_4^4}{576 \pi} \frac{((M_\chi - M_\psi)^2 - s) ((M_\chi + M_\psi)^2 + 2s)}{\sqrt{s^2 - 2s(M_\chi^2 + M_\psi^2) + (M_\chi^2 - M_\psi^2)^2} (s - M_U^2)^2} \\
&\quad \times \left( 2 \sum_{i,j=1}^3 |\beta_{ij}^{q\ell}|^2 + \sum_{i=1}^3 |\beta_i^{de}|^2 + \sum_{i=1}^3 |\beta_i^{u\nu}|^2 \right),
\end{split}
\tag{C.2}
\end{equation}
$$

$$
\sigma^{\psi\psi} \approx \frac{g_4^4 |\kappa_{\psi}|^2}{324 \pi c_{43}^4} \frac{s(2 M_{\psi}^2 + s)}{\sqrt{s^2 - 4s M_{\psi}^2 (s - M_{\tilde{G}'}^2)^2}} \sum_{i=1}^{3} (2|\kappa_q^i|^2 + |\kappa_u^i|^2 + |\kappa_d^i|^2). \quad (C.3)
$$

**Open Access.** This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.

References

[1] LHCb collaboration, Measurement of Form-Factor-Independent Observables in the Decay $B^0 \to K^{*0}\mu^+\mu^-$, Phys. Rev. Lett. 111 (2013) 191801 [arXiv:1308.1707] [INSPIRE].

[2] LHCb collaboration, Differential branching fractions and isospin asymmetries of $B \to K^{(*)}\mu^+\mu^-$ decays, JHEP 06 (2014) 133 [arXiv:1403.8044] [INSPIRE].

[3] LHCb collaboration, Angular analysis and differential branching fraction of the decay $B_s^0 \to \phi\mu^+\mu^-$, JHEP 09 (2015) 179 [arXiv:1506.08777] [INSPIRE].

[4] LHCb collaboration, Test of lepton universality using $B^+ \to K^+\ell^+\ell^-$ decays, Phys. Rev. Lett. 113 (2014) 151601 [arXiv:1406.6482] [INSPIRE].

[5] LHCb collaboration, Angular analysis of the $B^0 \to K^{*0}\mu^+\mu^-$ decay using 3 fb$^{-1}$ of integrated luminosity, JHEP 02 (2016) 104 [arXiv:1512.04442] [INSPIRE].

[6] ATLAS collaboration, Angular analysis of $B_d^0 \to K^*\mu^+\mu^-$ decays in pp collisions at $\sqrt{s} = 8$ TeV with the ATLAS detector, JHEP 10 (2018) 047 [arXiv:1805.04000] [INSPIRE].

[7] CMS collaboration, Measurement of the $P_1$ and $P'_5$ angular parameters of the decay $B^0 \to K^{*0}\mu^+\mu^-$ in proton-proton collisions at $\sqrt{s} = 8$ TeV, Tech. Rep. CMS-PAS-BPH-15-008 (2017).

[8] CMS collaboration, Angular analysis of the decay $B^0 \to K^{*0}\mu^+\mu^-$ from pp collisions at $\sqrt{s} = 8$ TeV, Phys. Lett. B 753 (2016) 424 [arXiv:1507.08126] [INSPIRE].
[9] LHCb collaboration, Test of lepton universality with $B^0 \to K^{*0}\ell^+\ell^-$ decays, JHEP 08 (2017) 055 [arXiv:1705.05802] [INSPIRE].

[10] LHCb collaboration, Search for lepton-universality violation in $B^+ \to K^+\ell^+\ell^-$ decays, Phys. Rev. Lett. 122 (2019) 191801 [arXiv:1903.09252] [INSPIRE].

[11] BELLE collaboration, Test of lepton flavor universality in $B \to K^+\ell^+\ell^-$ decays at Belle, arXiv:1904.02440 [INSPIRE].

[12] BABAR collaboration, Evidence for an excess of $\bar{B} \to D^{(*)}\tau^-\bar{\nu}_{\tau}$ decays, Phys. Rev. Lett. 109 (2012) 101802 [arXiv:1205.5442] [INSPIRE].

[13] BABAR collaboration, Measurement of an Excess of $\bar{B} \to D^{(*)}\tau^-\bar{\nu}_{\tau}$ Decays and Implications for Charged Higgs Bosons, Phys. Rev. D 88 (2013) 072012 [arXiv:1303.0571] [INSPIRE].

[14] BELLE collaboration, Measurement of the branching ratio of $\bar{B} \to D^{(*)}\tau^-\bar{\nu}_{\tau}$ relative to $\bar{B} \to D^{(*)}\ell^-\bar{\nu}_{\ell}$ decays with hadronic tagging at Belle, Phys. Rev. D 92 (2015) 072014 [arXiv:1507.0323] [INSPIRE].

[15] LHCb collaboration, Measurement of the ratio of branching fractions
$$B(\bar{B}^0 \to D^{*+}\tau^-\bar{\nu}_{\tau})/B(\bar{B}^0 \to D^{*+}\mu^-\bar{\nu}_{\mu})$$, Phys. Rev. Lett. 115 (2015) 111803 [Erratum ibid. 115 (2015) 159901] [arXiv:1506.08614] [INSPIRE].

[16] BELLE collaboration, Measurement of the branching ratio of $\bar{B}^0 \to D^{*+}\tau^-\bar{\nu}_{\tau}$ relative to $\bar{B}^0 \to D^{*+}\ell^-\bar{\nu}_{\ell}$ decays with a semileptonic tagging method, Phys. Rev. D 94 (2016) 072007 [arXiv:1607.07923] [INSPIRE].

[17] BELLE collaboration, Measurement of the $\tau$ lepton polarization and $R(D^*)$ in the decay $\bar{B} \to D^*\tau^-\bar{\nu}_{\tau}$, Phys. Rev. Lett. 118 (2017) 211801 [arXiv:1612.00529] [INSPIRE].

[18] LHCb collaboration, Measurement of the ratio of branching fractions
$$B(B_c^+ \to J/\psi\tau^+\nu_\tau)/B(B_c^+ \to J/\psi\mu^+\nu_\mu)$$, Phys. Rev. Lett. 120 (2018) 121801 [arXiv:1711.05623] [INSPIRE].

[19] LHCb collaboration, Measurement of the ratio of the $B^0 \to D^{*-}\tau^+\nu_\tau$ and $B^0 \to D^{*-}\mu^+\nu_\mu$ branching fractions using three-prong $\tau$-lepton decays, Phys. Rev. Lett. 120 (2018) 171802 [arXiv:1708.08856] [INSPIRE].

[20] BELLE collaboration, Measurement of $R(D)$ and $R(D^*)$ with a semileptonic tagging method, arXiv:1904.08794 [INSPIRE].

[21] BELLE collaboration, Measurement of $R(D)$ and $R(D^*)$ with a semileptonic tagging method, Phys. Rev. Lett. 124 (2020) 161803 [arXiv:1910.05864] [INSPIRE].

[22] R. Alonso, B. Grinstein and J. Martin Camalich, Lepton universality violation and lepton flavor conservation in B-meson decays, JHEP 10 (2015) 184 [arXiv:1505.05164] [INSPIRE].

[23] L. Calibbi, A. Crivellin and T. Ota, Effective Field Theory Approach to $b \to s\ell\ell^{(f)}$, $B \to K^{(*)}\nu\bar{\nu}$ and $B \to D^{(*)}\tau\nu$ with Third Generation Couplings, Phys. Rev. Lett. 115 (2015) 181801 [arXiv:1506.02661] [INSPIRE].

[24] R. Barbieri, G. Isidori, A. Pattori and F. Senia, Anomalies in B-decays and U(2) flavour symmetry, Eur. Phys. J. C 76 (2016) 67 [arXiv:1512.01560] [INSPIRE].

[25] G. Hiller, D. Loose and K. Schönwald, Leptoquark Flavor Patterns & B Decay Anomalies, JHEP 12 (2016) 027 [arXiv:1609.08895] [INSPIRE].
[26] B. Bhattacharya, A. Datta, J.-P. Guévin, D. London and R. Watanabe, *Simultaneous Explanation of the $R_K$ and $R_{D(\pm)}$ Puzzles: a Model Analysis*, JHEP 01 (2017) 015 [arXiv:1609.09078] [INSPIRE].

[27] D. Buttazzo, A. Greljo, G. Isidori and D. Marzocca, *B-physics anomalies: a guide to combined explanations*, JHEP 11 (2017) 044 [arXiv:1706.07808] [INSPIRE].

[28] L. Calibbi, A. Crivellin and T. Li, *Model of vector leptoquarks in view of the B-physics anomalies*, Phys. Rev. D 98 (2018) 115002 [arXiv:1709.00692] [INSPIRE].

[29] A. Angelescu, D. Bečirević, D.A. Faroughy and O. Sumensari, *Closing the window on single leptoquark solutions to the B-physics anomalies*, JHEP 10 (2018) 183 [arXiv:1808.08179] [INSPIRE].

[30] J. Kumar, D. London and R. Watanabe, *Combined Explanations of the $b \to s\mu^+\mu^-$ and $b \to c\tau^-\bar{\nu}$ Anomalies: a General Model Analysis*, Phys. Rev. D 99 (2019) 015007 [arXiv:1806.07403] [INSPIRE].

[31] D. Das, C. Hati, G. Kumar and N. Mahajan, *Towards a unified explanation of $R_{D(\pm)}$, $R_K$ and $(g-2)_\mu$ anomalies in a left-right model with leptoquarks*, Phys. Rev. D 94 (2016) 055034 [arXiv:1605.06313] [INSPIRE].

[32] A. Crivellin, D. Müller and T. Ota, *Simultaneous explanation of $R(D^{(*)})$ and $b \to s\mu^+\mu^-$: the last scalar leptoquarks standing*, JHEP 09 (2017) 040 [arXiv:1703.09226] [INSPIRE].

[33] D. Marzocca, *Addressing the B-physics anomalies in a fundamental Composite Higgs Model*, JHEP 07 (2018) 121 [arXiv:1803.10972] [INSPIRE].

[34] D. Bečirević, I. Doršner, S. Fajfer, N. Košnik, D.A. Faroughy and O. Sumensari, *Scalar leptoquarks from grand unified theories to accommodate the B-physics anomalies*, Phys. Rev. D 98 (2018) 055003 [arXiv:1806.05689] [INSPIRE].

[35] I. Bigaran, J. Gargalonis and R.R. Volkas, *A near-minimal leptoquark model for reconciling flavour anomalies and generating radiative neutrino masses*, JHEP 10 (2019) 106 [arXiv:1906.01870] [INSPIRE].

[36] A. Datta, J.L. Feng, S. Kamali and J. Kumar, *Resolving the $(g-2)_\mu$ and B Anomalies with Leptoquarks and a Dark Higgs Boson*, Phys. Rev. D 101 (2020) 035010 [arXiv:1908.08625] [INSPIRE].

[37] W. Altmannshofer, P.S.B. Dev, A. Soni and Y. Sui, *Addressing $R_{D(\pm)}$, $R_K(\pm)$, muon $g-2$ and ANITA anomalies in a minimal R-parity violating supersymmetric framework*, Phys. Rev. D 102 (2020) 015031 [arXiv:2002.12910] [INSPIRE].

[38] J. Aebischer, W. Altmannshofer, D. Guadagnoli, M. Reboud, P. Stangl and D.M. Straub, *B-decay discrepancies after Moriond 2019*, Eur. Phys. J. C 80 (2020) 252 [arXiv:1903.10434] [INSPIRE].

[39] M. Algueró et al., *Emerging patterns of New Physics with and without Lepton Flavour Universal contributions*, Eur. Phys. J. C 79 (2019) 714 [Addendum ibid. 80 (2020) 511] [arXiv:1903.09578] [INSPIRE].

[40] A. Crivellin, C. Greub, D. Müller and F. Saturnino, *Importance of Loop Effects in Explaining the Accumulated Evidence for New Physics in B Decays with a Vector Leptoquark*, Phys. Rev. Lett. 122 (2019) 011805 [arXiv:1807.02068] [INSPIRE].

[41] R. Barbieri, C.W. Murphy and F. Senia, *B-decay Anomalies in a Composite Leptoquark Model*, Eur. Phys. J. C 77 (2017) 8 [arXiv:1611.04930] [INSPIRE].
[42] J.C. Pati and A. Salam, Lepton Number as the Fourth Color, Phys. Rev. D 10 (1974) 275
[Erratum ibid. 11 (1975) 703] [INSPIRE].

[43] H. Georgi and Y. Nakai, Diphoton resonance from a new strong force, Phys. Rev. D 94 (2016) 075005 [arXiv:1606.05865] [INSPIRE].

[44] B. Diaz, M. Schmaltz and Y.-M. Zhong, The leptoquark Hunter's guide: Pair production, JHEP 10 (2017) 097 [arXiv:1706.05033] [INSPIRE].

[45] L. Di Luzio, A. Greljo and M. Nardecchia, Gauge leptoquark as the origin of B-physics anomalies, Phys. Rev. D 96 (2017) 115011 [arXiv:1708.08450] [INSPIRE].

[46] M. Blanke and A. Crivellin, B Meson Anomalies in a Pati-Salam Model within the Randall-Sundrum Background, Phys. Rev. Lett. 121 (2018) 011801 [arXiv:1801.07256] [INSPIRE].

[47] L. Di Luzio, J. Fuentes-Martín, A. Greljo, M. Nardecchia and S. Renner, Maximal Flavour Violation: a Cabibbo mechanism for leptoquarks, JHEP 11 (2018) 081 [arXiv:1808.00942] [INSPIRE].

[48] M. Bordone, C. Cornella, J. Fuentes-Martín and G. Isidori, A three-site gauge model for flavor hierarchies and flavor anomalies, Phys. Lett. B 779 (2018) 317 [arXiv:1712.01368] [INSPIRE].

[49] A. Greljo and B.A. Stefanek, Third family quark-lepton unification at the TeV scale, Phys. Lett. B 782 (2018) 131 [arXiv:1802.04274] [INSPIRE].

[50] M. Bordone, C. Cornella, J. Fuentes-Martín and G. Isidori, Low-energy signatures of the $\text{PS}^3$ model: from B-physics anomalies to LFV, JHEP 10 (2018) 148 [arXiv:1805.09328] [INSPIRE].

[51] C. Cornella, J. Fuentes-Martín and G. Isidori, Revisiting the vector leptoquark explanation of the B-physics anomalies, JHEP 07 (2019) 168 [arXiv:1903.11517] [INSPIRE].

[52] J. Fuentes-Martín and P. Stangl, Third-family quark-lepton unification with a fundamental composite Higgs, arXiv:2004.11376 [INSPIRE].

[53] R. Barbieri, G. Isidori, J. Jones-Perez, P. Lodone and D.M. Straub, U(2) and Minimal Flavour Violation in Supersymmetry, Eur. Phys. J. C 71 (2011) 1725 [arXiv:1105.2296] [INSPIRE].

[54] G. Blankenburg, G. Isidori and J. Jones-Perez, Neutrino Masses and LFV from Minimal Breaking of $U(3)^5$ and $U(2)^5$ flavor Symmetries, Eur. Phys. J. C 72 (2012) 2126 [arXiv:1204.0688] [INSPIRE].

[55] R. Barbieri, D. Buttazzo, F. Sala and D.M. Straub, Flavour physics from an approximate $U(2)^3$ symmetry, JHEP 07 (2012) 181 [arXiv:1203.4218] [INSPIRE].

[56] G. Bertone, D. Hooper and J. Silk, Particle dark matter: Evidence, candidates and constraints, Phys. Rept. 405 (2005) 279 [hep-ph/0404175] [INSPIRE].

[57] L. Roszkowski, E.M. Sessolo and S. Trojanowski, WIMP dark matter candidates and searches — current status and future prospects, Rept. Prog. Phys. 81 (2018) 066201 [arXiv:1707.06277] [INSPIRE].

[58] M. Cirelli, N. Fornengo and A. Strumia, Minimal dark matter, Nucl. Phys. B 753 (2006) 178 [hep-ph/0512090] [INSPIRE].
[59] M.J. Baker et al., The Coannihilation Codex, JHEP 12 (2015) 120 [arXiv:1510.03434] [INSPIRE].

[60] D.A. Faroughy, A. Greljo and J.F. Kamenik, Confronting lepton flavor universality violation in B decays with high-$p_T$ tau lepton searches at LHC, Phys. Lett. B 764 (2017) 126 [arXiv:1609.07138] [INSPIRE].

[61] M. Schmaltz and Y.-M. Zhong, The leptoquark Hunter's guide: large coupling, JHEP 01 (2019) 132 [arXiv:1810.10017] [INSPIRE].

[62] A. Greljo, J. Martin Camalich and J.D. Ruiz-Álvarez, Mono-$\tau$ Signatures at the LHC Constrain Explanations of B-decay Anomalies, Phys. Rev. Lett. 122 (2019) 131803 [arXiv:1811.07920] [INSPIRE].

[63] M.J. Baker, J. Fuentes-Martín, G. Isidori and M. König, High-$p_T$ signatures in vector-leptoquark models, Eur. Phys. J. C 79 (2019) 334 [arXiv:1901.10480] [INSPIRE].

[64] J. Fuentes-Martín, G. Isidori, M. König and N. Selimović, Vector Leptoquarks Beyond Tree Level, Phys. Rev. D 101 (2020) 035024 [arXiv:1910.13474] [INSPIRE].

[65] K. Griest and D. Seckel, Three exceptions in the calculation of relic abundances, Phys. Rev. D 43 (1991) 3191 [INSPIRE].

[66] E.W. Kolb and M.S. Turner, The Early Universe, Front. Phys. 69 (1990) 1 [INSPIRE].

[67] G. Jungman, M. Kamionkowski and K. Griest, Supersymmetric dark matter, Phys. Rept. 267 (1996) 195 [hep-ph/9506380] [INSPIRE].

[68] M. Srednicki, R. Watkins and K.A. Olive, Calculations of Relic Densities in the Early Universe, Nucl. Phys. B 310 (1988) 693 [INSPIRE].

[69] M. Cannoni, Relativistic $<\sigma v_{rel}>$ in the calculation of relics abundances: a closer look, Phys. Rev. D 89 (2014) 103533 [arXiv:1311.4494] [INSPIRE].

[70] LUX collaboration, Results from a search for dark matter in the complete LUX exposure, Phys. Rev. Lett. 118 (2017) 021303 [arXiv:1608.07648] [INSPIRE].

[71] PANDAX-II collaboration, Dark Matter Results From 54-Ton-Day Exposure of PandaX-II Experiment, Phys. Rev. Lett. 119 (2017) 181302 [arXiv:1708.06917] [INSPIRE].

[72] XENON collaboration, Dark Matter Search Results from a One Ton-Year Exposure of XENON1T, Phys. Rev. Lett. 121 (2018) 111302 [arXiv:1805.12562] [INSPIRE].

[73] P. Salati, Indirect and direct dark matter detection, PoS CARGESE2007 (2007) 009 [INSPIRE].

[74] N. Anand, A.L. Fitzpatrick and W.C. Haxton, Weakly interacting massive particle-nucleus elastic scattering response, Phys. Rev. C 89 (2014) 065501 [arXiv:1308.6288] [INSPIRE].

[75] J. Fan, M. Reece and L.-T. Wang, Non-relativistic effective theory of dark matter direct detection, JCAP 11 (2010) 042 [arXiv:1008.1591] [INSPIRE].

[76] A. Fitzpatrick, W. Haxton, E. Katz, N. Lubbers and Y. Xu, The Effective Field Theory of Dark Matter Direct Detection, JCAP 02 (2013) 004 [arXiv:1203.3542] [INSPIRE].

[77] A. Fitzpatrick, W. Haxton, E. Katz, N. Lubbers and Y. Xu, Model Independent Direct Detection Analyses, arXiv:1211.2818 [INSPIRE].

[78] V. Cirigliano, M.L. Graesser and G. Ovanesyan, WIMP-nucleus scattering in chiral effective theory, JHEP 10 (2012) 025 [arXiv:1205.2695] [INSPIRE].
[79] M. Cirelli, E. Del Nobile and P. Panci, Tools for model-independent bounds in direct dark matter searches, JCAP 10 (2013) 019 [arXiv:1307.5955] [INSPIRE].

[80] G. Barello, S. Chang and C.A. Newby, A Model Independent Approach to Inelastic Dark Matter Scattering, Phys. Rev. D 90 (2014) 094027 [arXiv:1409.0536] [INSPIRE].

[81] R.J. Hill and M.P. Solon, Standard Model anatomy of WIMP dark matter direct detection II: QCD analysis and hadronic matrix elements, Phys. Rev. D 91 (2015) 043505 [arXiv:1409.8290] [INSPIRE].

[82] M. Hoferichter, P. Klos and A. Schwenk, Chiral power counting of one- and two-body currents in direct detection of dark matter, Phys. Lett. B 746 (2015) 410 [arXiv:1503.04811] [INSPIRE].

[83] R. Catena and P. Gondolo, Global fits of the dark matter-nucleon effective interactions, JCAP 09 (2014) 045 [arXiv:1405.2637] [INSPIRE].

[84] R.J. Hill and M.P. Solon, WIMP-nucleon scattering with heavy WIMP effective theory, Phys. Rev. Lett. 112 (2014) 211602 [arXiv:1309.4092] [INSPIRE].

[85] R.J. Hill and M.P. Solon, Universal behavior in the scattering of heavy, weakly interacting dark matter on nuclear targets, Phys. Lett. B 707 (2012) 539 [arXiv:1111.0016] [INSPIRE].

[86] M. Hoferichter, P. Klos, J. Menéndez and A. Schwenk, Analysis strategies for general spin-independent WIMP-nucleus scattering, Phys. Rev. D 94 (2016) 063505 [arXiv:1605.08043] [INSPIRE].

[87] A. Kurylov and M. Kamionkowski, Generalized analysis of weakly interacting massive particle searches, Phys. Rev. D 69 (2004) 063503 [hep-ph/0307185] [INSPIRE].

[88] M. Pospelov and T. ter Veldhuis, Direct and indirect limits on the electromagnetic form-factors of WIMPs, Phys. Lett. B 480 (2000) 181 [hep-ph/0003010] [INSPIRE].

[89] J. Bagnasco, M. Dine and S.D. Thomas, Detecting technibaryon dark matter, Phys. Lett. B 320 (1994) 99 [hep-ph/9310290] [INSPIRE].

[90] F. Bishara, J. Brod, B. Grinstein and J. Zupan, Chiral Effective Theory of Dark Matter Direct Detection, JCAP 02 (2017) 009 [arXiv:1611.00368] [INSPIRE].

[91] F. Bishara, J. Brod, B. Grinstein and J. Zupan, From quarks to nucleons in dark matter direct detection, JHEP 11 (2017) 059 [arXiv:1707.06998] [INSPIRE].

[92] F. Bishara, J. Brod, B. Grinstein and J. Zupan, DirectDM: a tool for dark matter direct detection, arXiv:1708.02678 [INSPIRE].

[93] R.K. Leane, Indirect Detection of Dark Matter in the Galaxy, in 3rd World Summit on Exploring the Dark Side of the Universe, (2020) [arXiv:2006.00513] [INSPIRE].

[94] H.E.S.S. collaboration, Dark matter search in the inner galactic center halo with H.E.S.S, in 51st Rencontres de Morion on Cosmology, pp. 149–152 (2016) [arXiv:1608.08453] [INSPIRE].

[95] HESS collaboration, Search for γ-Ray Line Signals from Dark Matter Annihilations in the Inner Galactic Halo from 10 Years of Observations with H.E.S.S., Phys. Rev. Lett. 120 (2018) 201101 [arXiv:1805.05741] [INSPIRE].

[96] FERMI-LAT collaboration, Updated search for spectral lines from Galactic dark matter interactions with pass 8 data from the Fermi Large Area Telescope, Phys. Rev. D 91 (2015) 122002 [arXiv:1506.00013] [INSPIRE].
[97] J.F. Navarro et al., *The Diversity and Similarity of Cold Dark Matter Halos*, Mon. Not. Roy. Astron. Soc. **402** (2010) 21 [arXiv:0810.1522] [INSPIRE].

[98] J. Hisano, S. Matsumoto, M.M. Nojiri and O. Saito, *Non-perturbative effect on dark matter annihilation and gamma ray signature from galactic center*, Phys. Rev. D **71** (2005) 063528 [hep-ph/0412403] [INSPIRE].

[99] J. Hisano, S. Matsumoto, M.M. Nojiri and O. Saito, *Direct detection of the Wino and Higgsino-like neutralino dark matters at one-loop level*, Phys. Rev. D **71** (2005) 015007 [hep-ph/0407168] [INSPIRE].

[100] J. Hisano, S. Matsumoto and M.M. Nojiri, *Unitarity and higher order corrections in neutralino dark matter annihilation into two photons*, Phys. Rev. D **67** (2003) 075014 [hep-ph/0212022] [INSPIRE].

[101] J. Hisano, S. Matsumoto and M.M. Nojiri, *Explosive dark matter annihilation*, Phys. Rev. Lett. **92** (2004) 031303 [hep-ph/0307216] [INSPIRE].

[102] M. Cirelli, A. Strumia and M. Tamburini, *Cosmology and Astrophysics of Minimal Dark Matter*, Nucl. Phys. B **787** (2007) 152 [arXiv:0706.4071] [INSPIRE].

[103] J.F. Navarro, C.S. Frenk and S.D.M. White, *The Structure of cold dark matter halos*, Astrophys. J. **462** (1996) 563 [astro-ph/9508025] [INSPIRE].

[104] D. Aristizabal Sierra, F. Staub and A. Vicente, *Shedding light on the $b \to s$ anomalies with a dark sector*, Phys. Rev. D **92** (2015) 015001 [arXiv:1503.06077] [INSPIRE].

[105] G. Bélanger, C. Delaunay and S. Westhoff, *A Dark Matter Relic From Muon Anomalies*, Phys. Rev. D **92** (2015) 055021 [arXiv:1507.06660] [INSPIRE].

[106] B. Allanach, F.S. Queiroz, A. Strumia and S. Sun, *$Z'$ models for the LHCb and $g-2$ muon anomalies*, Phys. Rev. D **93** (2016) 055045 [Erratum ibid. **95** (2017) 119902] [arXiv:1511.07447] [INSPIRE].

[107] M. Bauer and M. Neubert, *Flavor anomalies, the 750 GeV diphoton excess, and a dark matter candidate*, Phys. Rev. D **93** (2016) 115030 [arXiv:1512.06828] [INSPIRE].

[108] A. Celis, W.-Z. Feng and M. Vollmann, *Dirac dark matter and $b \to s\ell^+\ell^-$ with U(1) gauge symmetry*, Phys. Rev. D **95** (2017) 035018 [arXiv:1608.03894] [INSPIRE].

[109] W. Altmannshofer, S. Gori, S. Profumo and F.S. Queiroz, *Explaining dark matter and B decay anomalies with an $L_\mu - L_\tau$ model*, JHEP **12** (2016) 106 [arXiv:1609.04026] [INSPIRE].

[110] P. Ko, T. Nomura and H. Okada, *A flavor dependent gauge symmetry, Predictive radiative seesaw and LHCb anomalies*, Phys. Lett. B **772** (2017) 547 [arXiv:1701.05788] [INSPIRE].

[111] P. Ko, T. Nomura and H. Okada, *Explaining $B \to K(*)\ell^+\ell^-$ anomaly by radiatively induced coupling in U(1)$_{\mu-\tau}$ gauge symmetry*, Phys. Rev. D **95** (2017) 111701 [arXiv:1702.02699] [INSPIRE].

[112] J.M. Cline, J.M. Cornell, D. London and R. Watanabe, *Hidden sector explanation of B-decay and cosmic ray anomalies*, Phys. Rev. D **95** (2017) 095015 [arXiv:1702.00395] [INSPIRE].

[113] F. Sala and D.M. Straub, *A New Light Particle in B Decays?*, Phys. Lett. B **774** (2017) 205 [arXiv:1704.06188] [INSPIRE].
[114] J. Ellis, M. Fairbairn and P. Tunney, Anomaly-Free Models for Flavour Anomalies, Eur. Phys. J. C 78 (2018) 238 [arXiv:1705.03447] [INSPIRE].

[115] J. Kawamura, S. Okawa and Y. Omura, Interplay between the $b \to s\ell\ell$ anomalies and dark matter physics, Phys. Rev. D 96 (2017) 075041 [arXiv:1706.04344] [INSPIRE].

[116] S. Baek, Dark matter contribution to $b \to s\mu^+\mu^-$ anomaly in local $U(1)_{L_\mu-L_\tau}$ model, Phys. Lett. B 781 (2018) 376 [arXiv:1707.04573] [INSPIRE].

[117] J.M. Cline, B decay anomalies and dark matter from vectorlike confinement, Phys. Rev. D 97 (2018) 015013 [arXiv:1710.02140] [INSPIRE].

[118] J.M. Cline and J.M. Cornell, $R(K^{(*)})$ from dark matter exchange, Phys. Lett. B 782 (2018) 232 [arXiv:1711.10770] [INSPIRE].

[119] L. Dhargyal, A simple model to explain observed muon sector anomalies and small neutrino masses, J. Phys. G 46 (2019) 125002 [arXiv:1711.09772] [INSPIRE].

[120] C.-W. Chiang and H. Okada, A simple model for explaining muon-related anomalies and dark matter, Int. J. Mod. Phys. A 34 (2019) 1950106 [arXiv:1711.07365] [INSPIRE].

[121] A. Vicente, Anomalies in $b \to s$ transitions and dark matter, Adv. High Energy Phys. 2018 (2018) 3905848 [arXiv:1803.04703] [INSPIRE].

[122] A. Falkowski, S.F. King, E. Perdomo and M. Pierre, Flavourful $Z'$ portal for vector-like neutrino Dark Matter and $R_{K(*)}$, JHEP 08 (2018) 061 [arXiv:1803.04430] [INSPIRE].

[123] G. Arcadi, T. Hugle and F.S. Queiroz, The Dark $L_\mu - L_\tau$ Rises via Kinetic Mixing, Phys. Lett. B 784 (2018) 151 [arXiv:1803.05723] [INSPIRE].

[124] S. Baek and C. Yu, Dark matter for $b \to s\mu^+\mu^-$ anomaly in a gauged $U(1)_X$ model, JHEP 11 (2018) 054 [arXiv:1806.05967] [INSPIRE].

[125] A. Azatov, D. Barducci, D. Ghosh, D. Marzocca and L. Ubaldi, Combined explanations of B-physics anomalies: the sterile neutrino solution, JHEP 10 (2018) 092 [arXiv:1807.10745] [INSPIRE].

[126] B. Barman, D. Borah, L. Mukherjee and S. Nandi, Correlating the anomalous results in $b \to s$ decays with inert Higgs doublet dark matter and muon ($g-2$), Phys. Rev. D 100 (2019) 115010 [arXiv:1808.06639] [INSPIRE].

[127] D.G. Cerdeño, A. Cheek, P. Martín-Ramiro and J.M. Moreno, B anomalies and dark matter: a complex connection, Eur. Phys. J. C 79 (2019) 517 [arXiv:1902.01789] [INSPIRE].

[128] S. Trifinopoulos, B-physics anomalies: The bridge between R-parity violating supersymmetry and flavored dark matter, Phys. Rev. D 100 (2019) 115022 [arXiv:1904.12940] [INSPIRE].

[129] L. Da Rold and F. Lamagna, A vector leptoquark for the B-physics anomalies from a composite GUT, JHEP 12 (2019) 112 [arXiv:1906.11666] [INSPIRE].

[130] J. Fuentes-Martín, M. Reig and A. Vicente, Strong CP problem with low-energy emergent QCD: The 4321 case, Phys. Rev. D 100 (2019) 115028 [arXiv:1907.02550] [INSPIRE].

[131] Z.-L. Han, R. Ding, S.-J. Lin and B. Zhu, Gauged $U(1)_{L_\mu-L_\tau}$ scotogenic model in light of $R_{K(*)}$ anomaly and AMS-02 positron excess, Eur. Phys. J. C 79 (2019) 1007 [arXiv:1908.07192] [INSPIRE].