Phenomenological consequences of an interacting multicomponent dark sector

We consider a dark sector model containing stable fermions charged under an unbroken $U(1)$ gauge interaction, with a massless dark photon as force carrier, and interacting with ordinary matter via scalar messengers. We study its early Universe evolution by solving a set of coupled Boltzmann equations that track the number density of the different species, as well as entropy and energy exchanges between the dark and visible sectors. Phenomenologically viable realizations include: i) a heavy (order 1 TeV or more) lepton-like dark fermion playing the role of the dark matter candidate, with various production mechanisms active depending on the strength of the dark-visible sector portal; ii) light (few GeV to few tens of GeV) quark-like dark fermions, stable but with suppressed relic densities; iii) an extra radiation component in Universe due to dark photons, with temperature constrained by cosmic microwave background data, and in turn preventing dark fermions to be lighter than about 1 GeV. Extra constraints on our scenario stem from dark matter direct detection searches: the elastic scattering on nuclei is driven by dipole or charge radius interactions mediated by either Standard Model or dark photons, providing long-range effects which, however, are not always dominant, as usually assumed in this context. Projected sensitivities for next-generation detectors cover a significant portion of the viable parameter space and are competitive with respect to the model-dependent constraints derived from the magnetic dipole moments of leptons and cooling of stellar systems.

We consider a dark sector model containing stable fermions charged under an unbroken U (1) gauge interaction, with a massless dark photon as force carrier, and interacting with ordinary matter via scalar messengers. We study its early Universe evolution by solving a set of coupled Boltzmann equations that track the number density of the different species, as well as entropy and energy exchanges between the dark and visible sectors. Phenomenologically viable realizations include: i) a heavy (order 1 TeV or more) lepton-like dark fermion playing the role of the dark matter candidate, with various production mechanisms active depending on the strength of the dark-visible sector portal; ii) light (few GeV to few tens of GeV) quark-like dark fermions, stable but with suppressed relic densities; iii) an extra radiation component in Universe due to dark photons, with temperature constrained by cosmic microwave background data, and in turn preventing dark fermions to be lighter than about 1 GeV. Extra constraints on our scenario stem from dark matter direct detection searches: the elastic scattering on nuclei is driven by dipole or charge radius interactions mediated by either Standard Model or dark photons, providing long-range effects which, however, are not always dominant, as usually assumed in this context. Projected sensitivities for next-generation detectors cover a significant portion of the viable parameter space and are competitive with respect to the model-dependent constraints derived from the magnetic dipole moments of leptons and cooling of stellar systems.

I. MOTIVATIONS AND SYNOPSIS
The existence of a multicomponent dark sector has been extensively discussed in the literature (see [1,2] for two recent reviews). Such framework generally includes many new states with no direct interactions with the Standard Model (SM) particles, but possibly interacting among themselves by means of new forces. Motivations for this construction have been put forward in a variety of different contexts, ranging, e.g., from beyond SM physics in connection to collider data and flavor anomalies, to explaining the nature of the dark matter component of the Universe, and to addressing possible shortcomings in the SM of cosmology.
In particular, regarding the dark matter problem, any non-relativistic stable dark state can potentially contribute to the Universe's matter budget. Because of the secluded nature of the dark sector which prevents large couplings to ordinary matter, these states automatically satisfy observational properties for dark matter, mostly derived under the assumption that the only relevant interaction between dark and ordinary matter is gravity. On the other hand, given the complexity of the dark sector, the phenomenology of dark matter candidates in this context could be richer than simply looking at gravitational effects. For example, dark matter itself could be multicomponent or in composite forms; dark sector interactions may lead to macroscopic effects and, for instance, impact on the paradigm in the SM of cosmology that dark matter should be described as a collisionless fluid.
In this paper we illustrate the interplay among different effects occurring when the dark sector contains several species. More explicitly, we will discuss the early Universe's thermal history in such a scenario and the generation of dark matter and other stable relics. One peculiarity is the fact that there are two reservoirs of states, ordinary and dark, and their temperatures are not necessarily the same. Therefore, a set of coupled Boltzmann equations, tracking at the same time the number density of the different species and the energy exchanges between the two sectors, needs to be considered.
To investigate explicitly this issue, we must first commit ourselves to a specific model of the dark sector (which we do in section II by considering a rather minimal setup). The choice of model provides an explicit spectrum of states within the dark sector, the interaction strengths among the dark states, and the strength of the portal interaction between the dark and SM states. These must be supplied in order to extract definite predictions. In particular, we shall assume that the dark force is long range, that is mediated by an unbroken U (1) gauge interaction. Regarding the particle content, besides the force carrier, a massless dark photon, we introduce a set of stable dark fermions charged under the U (1). One of these may account for most of dark matter in the Universe since it is rather heavy, at the TeV scale or above, and passes upper limits from self-interaction effects [3][4][5]. The others are much lighter, have suppressed relic abundances, but concur in determining the ratio between dark and visible photon temperatures at late times; such ratio is constrained by cosmic microwave background (CMB) data, given that dark photons contribute as an extra radiation component to the Universe's dynamics. In this respect, the role of portal interactions between dark and visible sectors is also important: we consider scalar messengers mediating Yukawa-like interactions. The latter are also crucial for selecting the mechanism for dark matter generation and final relic densities. Such interplay is discussed in detail in section III.
Direct detection, namely the attempt to measure nuclear recoils induces by dark matter scatterings, is one of the main tools to test a given dark matter scenario. In our framework, the direct-detection cross section is mostly driven, via loop induced magnetic dipole and charge radius interactions, by the massless mediators, SM and dark photons. While long-range interactions are present and boost the recoil spectrum at low recoil energies, the correlated contact terms are also contributing to the cross section and may be dominant (contrary to standard lore that contact interactions can be neglected in the presence of long-range effects). These aspects are illustrated in section IV, bridging also between astrophysical, cosmological, and high-energy observables and relative constraints, demonstrating once more the diversity of the phenomenological implications of introducing such a multicomponent dark sector.

II. A MODEL OF THE DARK SECTOR
Several dark sector models have been studied in the literature and they are usually classified [2] according to the portal through which they interact with ordinary matter. We consider a model consisting of dark fermions that are, by definition, singlets under the SM gauge interactions. These dark fermions interact with the visible sector through a portal provided by scalar messengers which carry both SM and dark-sector charges. These scalars are phenomenologically akin to the sfermions of supersymmetric models.
In general, we can have as many dark fermions as there are in the SM; they can be classified conveniently according to whether they couple (via the corresponding messengers) to quarks (q L , u R , d R ) or leptons (l L , e R ): we denote the former (hadron-like) Q and the latter (lepton-like) χ. The Yukawa-like interaction Lagrangian can be written as [6,7]: The L-type scalars are doublets under SU(2) L , while the R-type scalars are singlets under SU(2) L . The S L,R messengers carry color indices (unmarked in (2.1)), while the messengers φ L,R are color singlets. The Yukawa coupling strengths are parameterized by α L,R ≡ g 2 L,R /(4π); they can be different for different fermions and as many as the SM fermions.
In order to generate chirality-changing processes, we must have the mixing terms where H is the SM Higgs boson,H = iσ 2 H , and S 0 a scalar singlet of the dark sector. After both S 0 and H take a vacuum expectation value (VEV) (µ S and v-the electroweak VEV-respectively), the Lagrangian in Eq. (2.2) gives rise to the mixing between right-and left-handed states. Dark sector states interact by means of an unbroken U (1) D gauge symmetry; the corresponding massless gauge boson is the dark photon γ D whose coupling strength we denote by α D ≡ g 2 D /(4π). We assign different dark U (1) D charges to the various dark sector fermions to ensure, by charge conservation, their stability. There is no kinetic mixing between the ordinary and the dark photon [8,9]. The latter is a distinctive feature of models in which the dark photon is, and remains, massless as opposed to those in which the gauge symmetry is broken and the dark photon is massive. While there is no tree-level coupling between dark fermions and SM photons, and between ordinary matter and dark photons, the mixing in Eq. (2.2) leads, through one-loop diagrams and therefore operators of dimension larger than four, to an effective coupling of ordinary matter to the dark photon as well as of the dark fermions to the ordinary photon.
When the dark sector scalar S 0 and the Higgs boson acquire VEVs, the scalar messengers must be rotated to identify the physical states. Considering first the lepton sector, while before the rotation all φ states have the same mass m φ , after the rotation we find the mass eigenstates (labeled by ±) with masses m φ± = m φ √ 1 ± η s , where we defined the mixing parameter: The picture in the hadronic sector is perfectly specular; in the following we will indicate generically with m S the mass for the eigenstates S U Ld and S D Lu before the rotation, and keep η s as mixing parameter for the physical eigenstates: Looking at (2.5), we can see that for χ to be a stable dark-sector species, its mass must be at most m φ− + m e . Similarly, for a dark-sector species Q, the mass must be no heavier than m S− + m q , where m q is the mass of the SM species corresponding to Q. This sets an upper bound for the mixing η s : where M χ,Q stands for the mass of the heaviest stable dark-sector species and m φ,S for the mass parameter of the corresponding messenger. We assume that M χ,Q is much heavier than any SM species. The upper bound in Eq. (2.7) also guarantees that the scalar messengers are heavier than the dark fermion into which can thus decay. This model can be considered as a template for many models of the dark sector with the scalar messenger as stand-in for more complicated portals. It is a simplified version of the model in [6], which might provide a natural solution to the SM flavor-hierarchy problem. It has been used to predict new decays for the Higgs boson [10][11][12], neutral Kaons [13] and the Z-boson [14] as well as invisible decays for the neutral K-and B-mesons [15].
A. Constraining the model Several limits on the parameter space of the model are known from high-energy physics and tests in astrophysical and cosmological environments. We list below the most severe constraints and the relative implications for mass parameters and coupling constants, as a preliminary outline of the regions in parameter space which will be relevant in the analysis of dark matter candidates within this framework. These constraints will be discussed further in Section IV, when examining current limits and projected sensitivities from dark matter direct detection experiments.
Contrary to the case of a massive dark photon, constraints from flavor and precision physics, as well as radiative emission in astrophysical bodies, come from one-loop order corrections providing the coupling to SM fermions. Under the assumption of CP conservation in the dark sector, the limits quoted below are mostly derived from the effective magnetic moment of SM fermions with respect to the dark photon or the ordinary photon, induced by dark fermionscalar messenger loops. Since a change in the chirality of the fermions is required, the limits are strongly dependent to the mixing η s . Depending on the process under consideration, the experimental limits only constrain particular combinations of couplings and masses in the dark sector. At this level, it is then more useful to quote results for Yukawa couplings and dark-sector masses for specific flavors, rather than taking them to be universal as in Eq. (2.1).
• Precision physics: Magnetic dipole moments of leptons provide a deeper insight on the parameter space. From the experimental measurement of the electron magnetic dipole moment [35], we find: where m e χ stands for the mass of corresponding dark fermion. A comparable limit can be found from the experimental measurement of the muon magnetic dipole moment [36]: Since the measurement of the tau magnetic dipole moment is experimentally challenging, the corresponding limit is much less relevant, at about the GeV level.
Except for tau-like dark sector species, these limits point to lepton-like scalar messengers at a heavy scale, say 10 TeV or above, and lepton-like dark fermions significantly lighter, say at 1 TeV or below -unless the couplings α L or α R gets suppressed, or the mixing parameter η s is small.
• Collider physics: Direct searches for charged scalar particles at the LHC [37] set a limit [15] m i S 940 GeV , (2.10) for the messenger mass related to the dark fermions Q U and Q D , while [38] have set constraints on the mass of sleptons, which give the following lower bound on the mass of lepton-like scalar messengers: GeV. (2.11) The limit increases to 1.5 TeV if more families are included. No limits exist for the masses of the dark fermions from events in which they are produced because they are SM singlets and do not interact directly with the detector.
• Astrophysics probes: Dark sector species can change the energy transport in astrophysical environments. Constraints for models with a massless dark photon from astrophysics have been discussed in [39][40][41]. The most stringent limit comes from stellar cooling in globular clusters by dark-photon Bremsstrahlung emission of electrons scattering on 4 He nuclei; for a standard choice of environmental parameters, and an upper value of 10 erg g −1 s −1 on the extra cooling rate by exotic processes [42], we find: This limit applies specifically to the Yukawa coupling to electrons and the corresponding messenger state, and affects regions in parameter space analogous to the limit in Eq. (2.8). When considering, instead, extra cooling effects in supernovae, the most relevant process is the dark photon emission in nucleon-nucleon Bremsstrahlung. From the neutrino signal of supernova 1987A one can deduce: (2.13) The above limit applies to the Yukawa couplings of u and d quarks and the corresponding messenger states. This sets an impact on the parameter space analogous to the leptonic sector, except that, for quark-like dark fermions, we will also explore the possibility of larger mass splittings with respect to the messenger states, with m Q even at the GeV scale.
• Self-interactions for dark matter particles: As already anticipated, our scenario gets severely constrained for light dark matter candidates because of the long-range self-interactions induced by the U (1) D gauge symmetry. The most severe observational limits come from the impact on the dark matter density distribution in collapsed dark matter structures, rather than effects in the early Universe or the early stages of structure formation [3,4,22]. Bounds have been derived from the dynamics in merging clusters, such as the Bullet Cluster [43], the tidal disruption of dwarf satellites along their orbits in the host halo, and kinetic energy exchanges among dark matter particles in virialized halos. Among these limits, the latter turns out to be the most constraining: energy exchanges through dark matter self-interactions tend to isotropize dark matter velocity distributions, while there are galaxies whose gravitational potentials show a triaxial structure with significant velocity anisotropy. A limit has been derived by estimating an isotropization timescale (via hard scattering and cumulative effects of many interactions, with Debye screening taken into account) and comparing that timescale to the estimated age of the object [22]: a refinement of this limit involves tracking the evolution of the velocity anisotropy due to the energy transfer [5]. The ellipticity profile inferred for the galaxy NGC720, according to Ref [5] (see Fig. 4) sets a limit of about: where m χ here stands for the dark matter mass -anticipating that we will focus on a lepton-like dark fermion as dark matter candidate -and the α D scaling quoted this equation is approximate and comes from the leading m χ over α D scaling in the expression for the isotropization timescale. Note that the limit quoted here is subject to a number of uncertainties and assumptions; it is less stringent than earlier results, such as the original bound quoted from [3], as well about a factor of 3.5 weaker than [22] (see also, e.g., [44,45]). On the other hand, results on galaxies from N-body simulations in self-interacting dark matter cosmologies [46], taking into account predicted ellipticities and dark matter densities in the central regions, seem to go in the direction of milder constraints, at about the same level or slightly weaker than the value quoted in Eq. (2.14). This result is also subject to uncertainties, such as the role played by the central baryonic component of NGC720.
As benchmark avoiding self-interaction constraints we will consider cases with dark matter mass about 1 TeV and α D 10 −2 .

B. Reference framework and parameter space
Taking into account the emerging picture, we will consider a scenario with: i) scalar messengers as the heaviest states in the dark-sector, ii) a lepton-like dark fermion χ playing the role of dark matter, lighter than scalar messengers but at a comparable mass scale, and iii) two dark fermions Q U and Q D coupled to the quarks, which are much lighter than χ and representative of the light dark sector (we shall see that the masses of the light dark species turn out to be indirectly constrained by CMB limits on exotic radiation components). Unless comparing to specific observables, to keep the model numerically tractable -but also without losing any of the main trends -we will adopt a set of simplifying assumptions. We restrict ourselves to the case in which all messenger states have a degenerate mass spectrum defined by a single mass parameter m φ = m S and a single mixing parameter η s . For simplicity, the Yukawa couplings of all the dark fermions are also taken to be equal, and with α L = α R . The extra parameters we need to deal with are the mass of the dark matter candidate m χ , the common mass m Q for the two light quark-like dark fermions and the dark photon coupling α D .
The remainder of the paper is devoted to additional constraints coming from the thermal history of the Universe and dark matter searches.

A. General picture
The aim is to compute the cosmological relic density for the stable species in the dark sector. The technical calculation, via a set of coupled Boltzmann equations, is discussed in the next section. However, it is useful to illustrate first a few features characterizing our setup.
The lightest fermions of given dark charge, lepton-like or hadron-like, are stable, and their number density in the early Universe heat bath changes through processes involving pair productions and pair annihilations; initially in equilibrium (chemical equilibrium; see the discussion below for a clarification on this point), they decouple in the non-relativistic regime. Thus, they have a relic density which can be approximated by the celebrated "WIMP miracle" formula: where vσ χχ,QQ is the thermal average of the pair annihilation cross section for either χ or Q, including all kinematically allowed final states. However, there are two elements which make the computation in the case at hand more involved than in other WIMP setups. First, while one usually deals only with SM final states, the pair annihilation may involve both particles belonging to the dark sector and to the SM sector; the leading processes are into two dark photons and a pair of SM fermion-antifermion of the corresponding type, as shown in Fig. 1 for the QQ initial state. Assuming that s-wave processes dominate, the thermal average of the pair annihilation cross sections, in the limit of small temperature corrections and massless final states, are approximately given by: Substituting these approximate expressions into Eq. (3.1), one can find the preferred mass ranges for which Ω χ is at the level of the cosmological dark matter abundance, while Ω Q is instead negligible (fulfilling the scheme emerging from the set of constraints discussed in the previous section). Taking α D and α L to be O(10 −2 ), and messenger scalars lying around 10 TeV, we find that Ω χ h 2 ∼ 0.1 if m χ is in the 1-10 TeV range; χs predominantly annihilate into dark photons (SM fermions) if (α L /α D ) 2 (m χ /m φ ) 4 is much less than (greater than) unity. Requiring that Ω Q is at most 1% of the Universe's matter density, we find as a conservative upper bound on the masses of the hadron-like species m Q 100 GeV; Qs predominantly annihilate into dark photons.
The second point we need to pay attention to is the fact that "thermal bath" effects, neglected so far, can actually have a significant impact on the overall picture. Analogously to the photon in the SM sector, the dark photon is crucial in keeping dark sector particles at a common temperature via, e.g., the large energy exchanges in Compton-like dark fermion -dark photon elastic scatterings. These elastic scattering processes maintain kinetic equilibrium within the dark sector. Moreover, being a stable massless particle, the dark photon can potentially give a sizable contribution to the budget for the energy density in radiation in the Universe, even at epochs, such as recombination, at which extra radiation components are tightly constrained. The general picture is given schematically in Fig. 2. Assuming that the U (1) D coupling α D is perturbative but still sufficiently large, dark photon interactions (or, eventually, a chain of processes involving additional interactions with other mediators/forces in the dark sector) enforce that all dark sector particles in the thermal bath have a common temperature T d . Analogously, Compton scattering between SM photons and SM particles maintains kinetic equilibrium within the visible sector. However, the temperature T of the visible sector may be different from T d . In the regime at which messenger scalars are non-relativistic, and with their number densities suppressed, the communication between visible and dark sectors (both at the level of particle number-changing processes and elastic scatterings) is mostly regulated by the Yukawa-like interactions in Eq. (2.1). Let us first turn off the portal interactions, i.e. α L = α R = 0. In this case the thermal bath in the visible and dark sectors evolve independently, and one can track T and T d by imposing entropy conservation separately in each of the two sectors, see, e.g. [21]. The cooling process goes as the inverse of the scale factor plus a correction due to the change in effective number of relativistic degrees of freedom, when particles becoming non-relativistic transfer their entropy to lighter, relativistic states of the corresponding sector. We define the temperature ratio between dark and visible sectors at a given time t to be and consider some initial time t 0 , with the temperature in the visible sector denoted by T 0 , at which two sectors are already decoupled. Assuming that entropy densities in the dark and visible sectors, which are respectively given by: are separately conserved in a comoving volume, one finds that the temperature ratio at the CMB epoch is given by: where g * Sv (T ) counts the number of internal degrees of freedom (fermionic species are weighted by 7/8) for all SM particles that are relativistic at temperature T , and g * S d (T d ) is the analogous quantity in the dark sector. Evaluating this ratio is relevant since this is the epoch at which extra radiation components are most severely constrained by cosmological observables. The limit is usually given in terms of N ef f , the effective number of neutrino-like species, i.e. fully relativistic fermions with two internal degrees of freedom, and with a temperature which is a factor of (4/11) 1/3 cooler than photons. N ef f is related to the radiation energy density by: The Planck satellite has measured N ef f at the CMB epoch to be [47]: N ef f = 3.27 ± 0.15, 68% CL. Subtracting out the contribution from the three standard model neutrinos [48] N SM ef f = 3.046, and assuming that the dark photon is the only dark sector relativistic state at the CMB epoch, giving rise to the extra radiation component ρ r,d (T d (t)) = ρ γ (T (t)) ξ 4 (t), we can translate the upper limit on N ef f from Planck into a limit on the temperature ratio at the CMB epoch; one finds: The 2σ and 3σ upper limits are, respectively, about 0.59 and 0.63. Our reference dark sector framework consists of: the dark photon, one lepton-like Dirac fermion χ, and N Q light hadron-like Dirac fermions Q being relativistic at the initial time t 0 . From Eq. (3.5) we obtain g * S d (ξ 0 T 0 )/g * S d (ξ CMB T CMB ) = (7N Q + 11)/4. Even for a single family of dark hadrons (N Q = 2) we find ξ CMB ≈ 0.61 ξ 0 , in tension with the limit quoted in (3.7) if ξ 0 = 1 (namely T = T d at t = t 0 ). As we increase the number of light species in the dark-sector, this problem gets more severe. A possible way out is to relax the initial condition. In principle the picture with decoupled sectors can be extrapolated to T 0 as high as, say, the reheating temperature. One can then assume an initial temperature mismatch between the two sectors, with a cooler dark sector (i.e. ξ 0 < 1), and thus the dark photon contribution to the radiation component of the Universe can be made small relative to the visible sector contribution. Similar conclusions (for various implementations of the dark-sector portal) were reached in, e.g., [3,21,40,[49][50][51][52].
On the other hand, when the messenger portal is turned back on, allowing for non-vanishing Yukawa couplings α L and α R , energy (and entropy) can be exchanged between visible and dark sectors. Regardless of what is assumed for ξ 0 , even if the system is not initially in kinetic equilibrium, for couplings sufficiently large, we expect it to relax to a maximum entropy configuration with the two temperature in the two sectors that will tend to become equal. This brings back the problem of satisfying the bound on extra radiation component associated to the dark photon at the CMB epoch, and will effectively translate on an upper bound on the Yukawa couplings. Since α L and α R both enter in the discussion for kinetic and chemical equilibrium, these two aspects have to be considered at the same time, as we will do with the set of coupled Boltzmann equations that we introduce in the next subsection and solve numerically.

B. Boltzmann equations
Having highlighted above that SM and dark sector states may have, in general, different temperatures, T and T d respectively, it is useful to keep track of them separately. Hence, in what follows, we adopt the following notation: i d will generically indicate a species in the dark sector, while species in the visible sector will be denoted by i v ; i will, in general, stand for any species in either sector. To track the distribution function of a state i d , we follow [53,54] and consider the generic Boltzmann equation where f i d is the occupation number for the particle i d , L is the Liouville operator tracking the evolution in the Friedmann-Robertson-Walker (FRW) background, and C is the collision operator. The Liouville operator takes the form where p is the physical momentum of i d and H is the Hubble rate. In the early Universe, the Hubble rate is dominated by radiation components coming from the visible and dark sectors. The first Friedmann equation tells us that (3.10) In the dilute limit, the collision operator acting on f i d is driven by 2 → 2 processes, such as i d + j ↔ k + l. It is then obtained by summing terms of the form: are the usual phase-space integration factors. When tracking chemical equilibrium, i.e. the evolution of the number density of i d , only inelastic processes are relevant. Given the structure of our model, the relevant number changing processes for χ and Q states (for T d not too large) are all in the form of particle-antiparticle pair annihilation or creation (see Fig. 1), namely (3.12) The expression for

and equilibrium distributions with occupation numbers approximated as
and (iii) kinetic equilibrium among dark sector states as enforced by elastic scatterings on the dark photon. Following from (ii), one can safely assume that standard model states follow equilibrium distributions and, using conservation of energy, formally rewrite their occupation numbers in terms of thermal distributions for the dark sector states in the form Note that we have T rather than T d in the last expression. As for (iii), this implies that, for any dark sector state, one may assume that there is an overall scaling -only dependent on time -of the occupation numbers of dark sector species with respect to equilibrium distributions: with n i d and n i d , respectively. To find the evolution equations for the number densities n i d of the relevant dark-sector fermions, we take the zeroth-order moment of the Boltzmann equation to obtaiṅ It is understood that the sum over j d includes the dark photon. The thermally averaged cross section σv in (3.16) is defined, in terms of the corresponding Møller cross section, as .

(3.17)
Looking at (3.16), there are three independent variables: t, T , and T d . In the standard approach, one closes the system by assuming entropy conservation; this leads to a time-temperature relation. In our current set-up, however, the two sectors are allowed to exchange energy and entropy, and thus the entropy of either sector is neither conserved. Nevertheless, the time evolution of the entropies of both sectors will allow us to obtain a well-posed ODE system. In tracking the entropy of both sectors, we first need to introduce the definition of entropy of species i in terms of the occupation number f i . This is given by Its evolution can be obtained by differentiating s i with respect to time, and then using Boltzmann equation. We havė We then take the sum of (3.19) over dark-sector species. Using kinetic equilibrium and dilute limit assumptions, one obtains: (3.20) where A i d has been defined in Eq. (3.15) above, and In the sum over dark-sector/visible sector species, we only take those processes that involve the transfer of entropy from one sector to the other. To proceed further, it is appropriate to digress into the discussion of the elastic part of the collision operator. It encodes the processes of type i + B ↔ i + B, where i is some species scattering from bath particles B, which also contribute to the entropy transfer between the two sectors. As demonstrated in [55], it can be written as where C i+B↔i+B [f i ] is a Fokker-Planck type operator, given by In obtaining this expression, it is assumed that the momentum transfer between i and B is much smaller than the typical momentum of either species. The momentum transfer rate can be shown to be given by: where we identify the thermal average of the momentum transfer rate At this point we would like to emphasize the following: if the species i were non-relativistic, T i m i , one could ignore the dependence of γ iB on energy, and the thermal average may be safely replaced as For instance, this applies for the case of scatterings of non-relativistic DM particles from a bath of relativistic SM species (this is the limit applied, e.g., in [56]). In our scenario, however, we would also like to account for entropy transfers from the dark-sector to the visible sector; this situation corresponds to the case where the dark-sector species act as bath particles for scatterings of visible sector species. When the scattering species are relativistic, one needs to take into account the energy dependence of γ iB and perform the thermal average at each step in the numerical solution of the system of coupled differential equations; further details about the technical implementation of this term are given Appendix B.
We are now in the position to write down the evolution equations for the entropies in the visible and dark sectors; we have (see also the analogous set of equations in [57]) where we have introduced yet another thermal average .
(3.28) From Eq. (3.27) it is transparent that if T = T d at early times the entropy exchange processes balance out, as expected from the condition of thermal equilibrium. Also, once the dark sector particles decouple, the entropies of the two sectors are separately conserved. The approach to kinetic equilibrium between the two sectors will then be relevant if we start with an initial temperature asymmetry and when the heavy dark sector species are still relativistic. We choose to solve the system of coupled differential equations using the scale factor a as the independent variable.
Using Eq. (3.4), we rewrite the evolution equations for the entropies as evolution equations for the temperatures: where we have written explicitly that the Hubble rate H depends both on T and T d , see Eq. (3.10), we have defined and have normalized all number densities to the entropy density in the visible sector, defining Y i ≡ n i /s v , with i being any species -in the visible sector or in the dark sector. For such variables and again using the scale factor a as independent variable, the Boltzmann equation (3.16) takes the form: Equations (3.29) and (3.31) constitute the closed system of differential equations to be solved.

C. Numerical results
As mentioned at the end of Section II, we will consider a dark sector framework with the following fermionic content: (i) one lepton-like dark fermion χ, with mass m χ , which acts as our dark matter candidate, and (ii) two hadron-like states, with masses m Q U and m Q D , that are lighter than χ. The evolution of the number density of each dark sector fermionic species is governed by Eq. (3.31). Regarding scalar messengers, we assume them to be degenerate in mass such that they are specified by a single mass parameter m S , and a universal mixing η S . Meanwhile, the other parameters relevant for the discussion are: (i) the U (1) D gauge coupling α D , and (ii) the Yukawa-like couplings α L and α R , which are taken to be equal for simplicity. Despite the model residing in a seven-dimensional parameter space, main trends can be illustrated on benchmark cases. In particular, unless explicitly stated, we will start illustrating the framework by focusing on the following choice of parameters: We will then vary the Yukawa-like coupling α L and properly adjust m χ , so that the relic density of χ approximately matches the dark matter density in the Universe as measured from cosmological observations. In Fig. 4 we present results for the numerical solution of the Boltzmann code for four different sets of pairs (α L , m χ ). In each panel a solid line shows, as a function of the inverse of the temperature in the visible sector T , the evolution of the number density for χ, Q U and Q D , normalized to the entropy density in the visible sector; such evolution is followed from an initial time t 0 , with initial temperature T 0 = 10 8 GeV, to some low temperature at which all comoving number densities are frozen to their relic values. Y i d for each dark fermion species i d is compared to the corresponding Y i d ,eq (T ), namely the ratio between the equilibrium number density n i d ,eq (T d ) -assuming T d = T -and again s v (T ), which is shown as a dashed line. This comparison is relevant since the case of Y i d tracking Y i d ,eq corresponds to the species i d being in chemical equilibrium as well as kinetic equilibrium between visible and dark sectors. In each panel we also show, with a dash-dotted line, the temperature ratio between dark and visible sectors; values of ξ(t) = T d /T can be read on the vertical scale on the right-hand side -notice that, to show more clearly its variation over time, the range of values displayed is adjusted in each panel (while the displayed range for Y i d , on the left-hand side of each panel, is kept fixed). Following the general discussion in Section III A, for all benchmark models considered in the plot, it is assumed that at t 0 the dark sector is significantly colder than the visible sector, starting the numerical solution with ξ(t 0 ) = 0.1.
In the four panels of Fig. 4, going from left to right and top to bottom, α L is progressively increased from a relatively small value for which the entropy exchanges between dark and visible sector are inefficient at any time, up to a regime at which kinetic equilibrium between the two sectors is reached at the very beginning of the numerical solution and maintained at temperatures lower than the chemical decoupling temperature of the lightest dark fermion. The values of the couplings, the dark fermion mass spectrum, as well as the results of the relic densities of the three dark fermions, and the value ξ CMB of the temperature ratio at the CMB epoch, are given in Table I. To explain trends in Fig. 4, considering the same benchmark cases and focusing on χ, in Fig 5 the effective interaction rates for relevant processes in Eqs. (3.29) and (3.31) are compared to the Universe's expansion rate H (as usual, as a rule of thumb, a given process is efficient only when the ratio is larger than one). The pair annihilation rates into dark photons and/or SM leptons, which are shown separately, drive chemical decoupling; the role of χ in restoring and maintaining kinetic equilibrium can be sketched from the effective energy transfer rate from dark fermion annihilations and χ elastic scattering on SM leptons, i.e. the combinations one obtains when factorizing out   Fig. (4), as well as their corresponding results for the relic densities and temperature ratio at CMB. We have chosen the couplings and masses such that χ would give a relic density that is close to the measured value of the matter density: Ωh 2 = 0.1186. In all cases, we have taken αD = 10 −2 and ms = 10 TeV.
the second and third term on the r.h.s. of Eq. (3.29). In the same plot we also show that, for all benchmark models, the scattering rate of χ on dark photons is much larger than H at any temperature, justifying the assumption of kinetic equilibrium among dark sector states. The four panels in Figs. 4 and 5 correspond to four different regimes in the parameter space. These are: • Region I (top-left plots) This is the regime in which the portal between dark and visible sectors is virtually absent, and the pair annihilation into dark photons enforces chemical equilibrium of fermions in the dark sector at large temperatures. In this case the relic density of χ can be estimated as the thermal freeze-out of a nonrelativistic species from the dark sector, which is analogous to the freeze-out of a standard WIMP from the visible sector: Y χ at freeze-out can be shown to be with the dark-sector freeze-out temperature being about The relic density of χ is then In this regime, the evolution of ξ is obtained by assuming that the entropies of the dark and visible sectors are separately conserved. Note that due to the Universe's expansion, both T d and T decrease; the temperature ratio ξ = T d /T increases whenever T d decreases slower than T . This occurs when a dark species becomes non-relativistic and heats up the dark photon plasma. The ratio reaches a peak at around T = m Q U , and then decreases since SM photons are heated up by SM degrees of freedom becoming non-relativistic, and, especially, at the QCD phase transition when quarks and gluons are transformed into bound-state hadrons.
• Region II (top-right plots): The moderate increase in α L is still insufficient to reach kinetic equilibrium between the two sectors. The effective energy transfer rate and the elastic scattering rates are still smaller than H at all temperatures, see Fig. 5. Nevertheless, the entropy leakage between the two sectors cannot be ignored, as one can see in the partial readjustment of ξ in the top-right panel of Fig. 4. Meanwhile in this regime, the relic density for χ, while still mostly determined by χ pair annihilations into dark photons, is also dictated by pair annihilations of visible sector particles populating the dark-sector with more dark fermions. This scheme is reminiscent of the freeze-in production mechanism for feebly interacting massive particles (FIMPs) [58]. As an approximate expression, Eq. The plot refers to the same four benchmark models displayed in Fig. 6 and specified in Table I, as representative of regions I, II, III, and IV in the parameter space (from left to right and top to bottom).
• Region III (bottom-left plots): This is the regime in which α L is large enough to enforce kinetic equilibrium between the two sectors from the very first steps of the numerical solution, up to the freeze out temperature of the dark matter component (but -for the specific parameter choice displayed -not up to the temperature at which the light fermions become non-relativistic). It is however still too small for the χ pair annihilation into SM leptons to play a role in setting the dark matter relic density; the annihilation into dark photons is still the dominant channel and the standard WIMP formula, Eq. (3.1) applies. Notice that the peak in the temperature ratio exceeds unity, since the light fermions become non-relativistic after kinetic decoupling. It follows that this is the benchmark case with largest ξ CMB , slightly above the 1 σ bound from Planck.
• Region IV (bottom-right plots): This scenario is similar to region III, except that, concerning the relic density of χ, α L is sufficiently large for SM lepton-anti-lepton pairs to be the dominant final state in the annihilation rate driving the WIMP rule-of-thumb formula Eq. (3.1). In the case at hand, α L is also large enough to ensure kinetic equilibrium between dark and visible sectors at all temperatures at which dark fermions are relativistic, hence ξ becomes 1 immediately after t 0 and is not increasing further.
The four regions are also shown in the left panel of Fig. 6, where the relic density of χ is plotted as a function of α L . We have kept m χ , α D , and η S fixed for each curve. As expected from the previous discussion, Ω χ h 2 is not necessarily a monotonic function of α L and so there are multiple values of α L giving the same relic density. In regions I and III, Ω χ h 2 is independent of α L , since in both regimes it is the annihilation to dark photons that determines the relic density of χ. Note that region I is the regime where the dark sector out of kinetic equilibrium with respect to the visible sector at all times, while region III is the regime where kinetic equilibrium holds until, at least, the chemical freeze-out of χ. Region II is the transition region between I and III: since the energy/entropy transfers and freeze-in effects become more efficient as α L increases, Ω χ h 2 increases as well. Region IV is the regime in which annihilations into SM leptons become dominant: following from Eqs. (3.1) and (3.2), we have Ω χ h 2 ∝ α −2 L . The change in m χ produces a vertical shift of regions I, II, and III. This follows from the fact that, for these regimes, the relic density of χ is determined by the annihilation to dark photons, and thus Ω χ h 2 ∝ m 2 χ . The trend changes for region IV; we have Ω χ h 2 ∝ m −2 χ . We also include the case where the left-right mixing between scalar messengers is maximal, i.e. η S = 1 − (m χ /m φ ) 2 ; this makes one of the scalar messengers lighter. A lighter scalar messenger increases the rate of processes enforcing kinetic equilibrium, which slightly changes the transition in region II; it also increases the annihilation rate to SM fermions, leading to the transition from region III to IV at a smaller α L , as well as it leads to a decrease in the relic density in region IV. The value α L * of the transition between regions III and IV can be roughly estimated by imposing that the annihilation cross section to SM species is about the same as the annihilation to dark photons; this leads to For instance, if η s = 0, α D = 10 −2 , and m S /m χ = 2 (as in the blue curve in the left panel of Fig. 6), we have α L * 4 × 10 −2 . In the right panel of Fig. 6 we explore how the relic density of the lighter dark fermions change with their masses. Given the constraints on light particles with long-range interactions in DM halos, such relic densities must be much suppressed compared to Ω χ h 2 , at the level of about 1% or lower. In general, the lighter the dark fermion, the more efficient the pair production/annihilation is into dark photons; since chemical decoupling is regulated by this final state, the relic density decreases accordingly. The right panel of Fig. 6 indeed shows the expected scaling In the left panel of Fig. 7, we show the temperature ratio at the CMB as a function of α L , for fixed α D , m χ , and m φ , while N Q , the number of light dark quarks, and m Q , the common mass of the dark quarks, are allowed to change individually. On the right panel of Fig. 7, we present a contour plot of ξ CMB in the m Q − α L plane, for Each colored regions correspond to 2−σ (green), 3−σ (orange), and > 3-σ (red) bands. The remaining regions correspond to ξCMB that are not excluded at 1σ by the current CMB limit on N ef f . The vertical lines correspond to half the masses of the neutral mesons KL and B 0 , which could decay into a particle-anti-particle pair of dark quarks. (e.g. see [15]).
fixed α D = 10 −2 , m χ = 1 TeV, and N Q = 2. The contour plot has been generated by performing a scan of m Q from 10 MeV to 300 GeV, and α L values from 10 −10 to 10 −1 . All results in Fig. 7 are obtained in numerical solutions of the Boltzmann code assuming as initial temperature ratio ξ 0 = 0.1. As previously mentioned in Sec. III A, bounds on N ef f constrain extra contributions to the amount of radiation energy density. This constraint translates to an upper bound on the temperature ratio at CMB, given by Eq. (3.7).
There are a few features emerging from Fig. 7. As expected, at any given α L , the ratio ξ CMB increases as the number of light species N Q increases. In particular, in the limit of vanishing Yukawa coupling α L , i.e. when the two sectors do not communicate with each other, ξ CMB depends on N Q only. For our reference model, the scaling is ξ CMB ∝ (7N Q + 11) 1/3 . This follows from the fact that entropy is injected into the dark sector bath when dark species become non-relativistic. Since the CMB epoch occurs at relatively late times, ξ CMB does not depend on m Q . Recall also that in this limit, ξ CMB ∝ ξ 0 and we are assuming a rather small ξ 0 .
Starting from a vanishingly small α L , entropy exchanges between visible and dark sectors, that tend to equilibrate the mismatch ξ 0 in the initial temperatures, become more efficient as we increase α L . This leads to increasing ξ CMB . In the left panel of Fig. 7 this is the rising branch at α L 10 −6 . The largest increase is obtained at some intermediate α L for which kinetic equilibrium is reached at early times, but is not maintained at the epoch at which χ or the light dark fermions become non-relativistic. When these particles become non-relativistic, they transfer their entropies mainly to dark photons, which makes ξ(t) become larger than 1 at some intermediate temperatures.
If instead α L is large enough to maintain kinetic equilibrium when dark fermions become non-relativistic, entropy injections are shared by the SM degrees of freedom and the result is a decrease in ξ CMB . At the same time the reverse effect occurs: SM states becoming non-relativistic and injecting entropy into the dark sector, rather than just heating SM photons, with then an increase in ξ CMB . The efficiency in these two-direction exchanges clearly depends on all parameters regulating kinetic equilibrium between the two sectors, including m χ , m Q , and the messenger masses m φ and m S , as well as on the parameters setting the temperatures at which the dark fermions become non-relativistic (regulated also by m χ and m Q ). In the left panel of the figure, we show in particular the α L dependence of ξ CMB for different values m Q and N Q = 1, while the case N Q = 2 is illustrated for a sample value in the left panel and in the full range m Q ∈ (10 MeV, 300 GeV) in the right panel. As the entropy transfer is particularly large at the QCD phase transition, at a temperature of about 150 MeV [59], it is crucial whether, at this epoch, Q are relativistic and/or visible and dark sectors are in kinetic equilibrium.
As seen from Fig. 7, the CMB limits on N ef f turn out to be a very severe constraint on the content of light fermions in the dark sector. Assuming an α L of at least 10 −2 , a favorable situation in order to satisfy the CMB limits at 1-σ level would be to keep N Q ≤ 2 and take m Q to be at least 5 GeV. Future tighter constraints on N ef f will impact on the parameter space even more severely.

IV. DIRECT DETECTION SEARCHES
Direct searches test the interactions of dark matter particles with ordinary matter. As a preliminary step to project direct detection limits into our framework, we need to write down the effective coupling between dark leptons and quarks. Scattering processes are mostly driven by massless mediators: SM and dark photons. Since there is no kinetic mixing between the SM and the dark photon, the leading contributions appear at one-loop order, as shown in Fig. 8.
Computing the diagrams in Fig. (8) yields the following dimension 5 (magnetic dipole) and dimension 6 (chargeradius) effective operators 1 : where F µν and X µν are, respectively, the field strength associated with the SM photon and the dark photon. The dipole and charge-radius couplings, denoted by d M /Λ D and c CR /[Λ CR ] 2 , respectively, carry additional labels. These additional labels specify: the fermion they are associated with, and the massless gauge boson such fermion is coupled to. In the discussion below we will both show results referring to a generic framework in which dipole and chargeradius couplings are treated independently of each other, as well as focus on our specific framework; in the latter case, they are given in terms of our model parameters and strong correlations appear. In particular, assuming universal couplings and g L = g R , we have The exact expressions for the functions F D and F CR , which are either of order 1 or logarithmically-enhanced, are 1 We assume for simplicity that CP invariance is respected in the dark-sector and there are no electric dipole moments.
given in Appendix C. Here we just quote useful approximate expressions assuming that m l m φ− and m Q m S− : If we plug-in typical values of dark-visible couplings and particle masses in the dark sector (where we take m φ+ = m S+ = 10 3 TeV), we have

(4.5)
A further correlation is with the additional contribution to the magnetic dipole moment of SM leptons predicted in our framework; this involves a single class of loop diagrams in which the virtual messenger scalars couple with the SM photon, which is analogous with the bottom-right diagram in Fig. (8) with the particles χ and e exchanged. Such contribution is given by where again F D is order one and can be approximated as (4.7) In the following we will first analyze separately the cases in which DM-nucleus scattering is mediated by: (a) the SM photon, and (b) the dark photon. Case (a) has already been discussed in the literature, and some results are reproduced here (see, e.g., [60][61][62][63][64]; see [65] for another possible long-range interaction for dark matter). Case (b), in which nuclei carry a (dark) magnetic moment, is explored here for the first time. We discuss the differential recoil rates, exclusion curves and projected sensitivities that one obtains considering each of the two massless mediators. Since both cases (a) and (b) involve dipole-vector interactions between DM and nucleons, one expects a term in the scattering amplitude which scales as the inverse of the momentum transfer, giving it an enhancement in the recoil rate at small recoil energies. On the other hand, we demonstrate here that one cannot naively conclude that the latter is the dominant effect and neglect other terms. While the enhancement is indeed present, it may get dominant over other terms only at extremely small recoil energies. It follows that, for what regards the phenomenology of the model, dimension 5 operators are not always playing the main role.

A. Direct detection analysis: an overview
The direct detection differential recoil rate, namely the number of scattering events per unit time, detector mass and recoil energy, can be generally written as In this equation the product of | v |, the modulus of the velocity of the DM particle in the detector frame, times the local DM particle number density, expressed in terms of the ratio between the local DM density ρ 0 and the DM mass m χ , gives the flux of DM particles in the detector at given | v |. Such flux is weighted over the velocity distribution for DM particles in the detector frame f ( v) and convolved with the DM-nucleus differential cross section dσ T /dE R . The sum in the equation is over target nuclear isotopes T , with mass m T and relative abundance c T . The integral includes any | v | large enough to give a recoil energy E R , i.e., larger than v min = | q |/(2µ χT ), where µ χT is the target nucleus-DM reduced mass, µ χT = m χ m T /(m χ + m T ), and the momentum transfer | q | is related to the value of the recoil energy via E R = | q | 2 /(2m T ). In what follows, for the astrophysical dependent quantities ρ 0 and f ( v), we just refer to the standard assumptions in the direct detection community: a local DM halo density of 0.3 GeV/cm 3 and a Maxwellian velocity distribution in the Galactic frame, with standard values of the velocity dispersion, and of the circular and escape velocities at the position of the Sun. While results are mildly dependent on these assumptions, they do not affect in any way the general discussion. The DM-nucleus differential cross section dσ T /dE R is derived in steps. Given the coupling of DM with quarks, one retrieves the effective coupling of DM with nucleons. The general formalism developed to describe non-relativistic EFT interactions goes as follows: the non-relativistic reduction of the Lagrangian density for the elastic scattering of a heavy DM particle on a proton or a neutron at rest can be written in terms of a set of 15 hermitian, leading-order operators (see. e.g., [66,67]), i.e.: is built out of a different contraction of four three-vectors: the momentum transfer q; the transverse component of the DM particle velocity v ⊥ ( v ⊥ · q = 0); the spin of the DM particle and of the nucleon, respectively S χ and S N . The second step is mapping the single-nucleon interactions into nuclear interactions; the general structure for the differential cross section takes the form: where the proton-neutron basis has been replaced by the isospin basis, τ and τ are isospin indices, S α are the darkmatter response functions containing contractions of O (N ) i terms and depend on the coefficients appearing in (4.9), v ⊥ ≡ v + q/(2µ), and W α are the nuclear response functions which are essentially form factors accounting for the composite structure of the nucleus.
Once we have the differential recoil rate, the expected number of direct detection events can be computed using [68]: where M is the mass of the detector, T E is the exposure time, and φ(E R ) is the efficiency curve specific to a particular experiment. We can then use the data on the observed number of scattering events in a direct detection experiment, to constrain DM-nucleon interactions. To obtain the usual exclusion curves with some specified confidence level 1 − α, one must, in principle, obtain the confidence interval [0, N p * ] from the posterior probability distribution of N p , given the observed number of events N o . A fixed value of N p * corresponds to a contour in the space of parameters that we are trying to constrain. Alternatively, to obtain exclusion plots, we use here the likelihood ratio test. First compute the Poisson likelihood functions: 12) where b is the number of background events, and then obtain the test statistic The test statistic λ follows a half-chi-squared distribution with one degree of freedom. The exclusion region will then correspond to those values of N p , which give probabilities above the confidence level: for 90% CL, we reject those values of N p which give λ −1. 64. In what follows we shall use DDCalc [68,69], a package written specifically for dark-matter direct detection calculations, including the calculation of differential recoil rates and likelihoods needed for obtaining parameter constraints at some specified confidence level. We will apply the procedure above to compare against the latest results from the XENON collaboration, which has produced the strongest upper limits in the DM particle mass range of interest for our framework [70], and to infer projected sensitivities of one of the proposed next-generation direct detection experiments, the DARWIN experiment [71], as representative of nearly final target for the direct detection field. In both cases we have checked that our results match closely published results when the DM nucleus interaction is assumed to be mediated by the standard spin independent operator.

B. SM photon-mediated processes
We consider first interactions mediated by SM photons (abbreviated as γm in the following). The dipole and charge radius effective coupling between dark leptons and SM quarks can be readily extracted from the effective operators in Eqs. (4.1) and (4.2): where q µ is the transfer four-momentum. We map the quark operators to the nucleon operators by using the form factors in [72]. We havē (4.15) where N = n, p, and the F (q/N ) i coefficients are QCD matrix elements. Applying (4.15) to the quark vector current in (4.14) we get (niσ µα q α n).(4.16) Following the prescription for mapping dark-matter-nucleon operators to their non-relativistic counterparts [66,73], the effective, non-relativistic DM-nucleon interaction is .

(4.17)
Here we have adopted the standard operator numbering where O 1 and O 4 are the operators commonly labelled as, respectively, spin-independent and spin-dependent couplings, and Notice that we have organized the terms in (4.17) in powers of | q |; the first line is of order 1/| q |, while the second line is of order | q | 0 . Looking at (4.17), we see that the γm dipole interaction gives an O 5 contribution, which is   long-range and coherent, a O 1 contribution, which is a contact term and coherent, and other short-range, incoherent contributions. On the other hand, the γm CR interaction gives only a contact, coherent O 1 contribution. We summarize these information in Table II. Coherent terms are likely to provide the largest contributions to the recoil spectrum. Depending on the relative size of the corresponding couplings, the recoil spectrum can either be dominated by dipole or charge-radius interactions. We address this issue by treating first the two couplings, namely d CR,γ ] 2 , as independent free parameters. In the right panel of Fig. 9, assuming that only one of them is non-zero, we show the 90% confidence level exclusion curve from XENON1T data and the projected sensitivity curve for DARWIN as a function of the dark matter mass m χ . Solid lines refer to the case when the γm CR interaction is switched off, with values of γm dipole coupling shown on vertical axis on the left-hand side; on the other hand, dashed lines assume that γm dipole interactions are negligible, with values of the γm CR coupling displayed on the scale on the right-hand side. In the left panel of Fig. 9, we show instead exclusion and sensitivity curves in the plane d . In this plot the solid diagonal line, which runs through the area where exclusion and sensitivity curves bend, approximately marks the separation between the dipole-dominated (region above the line) and the CR-dominated regimes (region below the line). In fact, looking at the expression for the recoil rate contribution from γm dipole interactions, this is mostly driven by the long-range and coherent O 5 operator and can be approximated as: the γm CR contribution is instead of the form: Using Xe as nuclear target, and considering experiments which lose sensitivity below a recoil energy of few keV, we find: which is about the delimiter shown in the plot. There are additional information displayed in Fig. 9. The orange polygonal region in the left panel denotes the pairs of dipole-CR coefficients obtainable in our model assuming α L = 10 −1 , m χ ∈ [200 GeV, 2 TeV], m φ− ∈ [1 TeV, 100 TeV], m φ+ ∈ 11 TeV, 10 3 TeV , and m χ ≤ m φ− ≤ m φ+ . As it can be seen, there are models in our framework that are excluded by XENON1T data, while DARWIN will cut deeper into the parameter space. The full region is within the area delimited by the condition in (4.22). Hence we can infer that within our framework, for what concerns γm interactions, the dipole term contributes more to the direct detection rate than the CR term, although the latter can be relevant as well. Note that this statement depends on the type of the nuclear target and on the range of recoil energies at which the experiment is sensitive.  9: (Left) 90% confidence level exclusion curves from XENON1T data and the projected sensitivity curves for DARWIN for a few values of the dark matter mass mχ in the plane dipole coupling versus CR coupling. The diagonal line gives a visual guidance to separate the regime in which, for a Xenon target and typical detector setups, direct detection rates are dominated by γm dipole interactions or γm CR interactions. The orange region is the area spanned by a large sample of models within our dark sector setup. The horizontal and vertical lines represent a projection of the muon magnetic dipole moment limit into a limit on, respectively, the dipole and CR coefficients, within our framework and for two representative cases: a model with large mixing for scalar messenger (Proj. 1) and one with small mixing (Proj. 2), see the text for details; the intersection points only should be compared with the result for XENON1T and DARWIN. (Right) Exclusion and projected sensitivity curves (90% CL) versus dark matter mass in case of either γm dipole interactions only (solid lines, referring to the vertical scale displayed on the left-hand side) or γm CR interactions only (dashed lines, referring to the vertical scale displayed on the right-hand side); also shown are the limits on the γm dipole coefficient derived within our framework and the same parameter choices as in the left panel.
Finally, in Fig. 9 we try to compare the direct detection limits and projected sensitivities with other constraints. There is no other process in which the operators introduced in Eq. (4.14) are tested at a significant level, and hence a model independent comparison is not possible. On the other hand, as described above, within our framework the loop diagrams giving rise to these interactions are closely related to the loop diagrams contributing to the magnetic dipole moments of leptons, which in turn are providing among the tightest constraints on our model, recall the discussion in Section II A. For reference, we consider the case in which the dark matter particle χ is coupled to muons (stronger constraints would follow in case of coupling to electrons; the limits get essentially irrelevant in case of coupling to tau leptons). The relation between coefficients of the different operators is simply: (4.23) Comparing against the experimental measurement of the muon magnetic dipole moment [36], we find: We project this limit into a limit on the γm dipole and γm CR coefficients (hence comparing at this level against direct detection) choosing two representative set of values for the masses of the corresponding scalar messenger: in the first -to which we refer as projection 1 -we choose a large mixing configuration m φ− , m φ+ = 10 TeV, 10 3 TeV , while in the other -to which we refer as projection 2 -we consider a small mixing case m φ− , m φ+ = (10 TeV, 11 TeV).
In the left panel of Fig. 9 derived limits on γm dipole and γm CR coefficients are shown, respectively, as horizontal and vertical lines; the line-style reflects again the three sample choices for m χ and the position of the crossing point of horizontal and vertical lines for the same model configuration should be compared to the corresponding direct detection curves (crossing points correspond to physical models in our framework, and, as expected, they all lie in the dipole dominated region). We see that in general, within our framework, the muon magnetic dipole moment limit is more constraining than the current direct detection limit. On the other hand, future detectors will be more sensitive to smaller dark matter dipole moments. Note that the effective dipole operator requires a change in the chirality of the external fermion, either through a sizable η s or a mass insertion on the external leg. When η s is sufficiently small, i.e. η s m µ /m χ , the muon dipole is proportional to m µ while the dark matter dipole is proportional to m χ : in this case the muon dipole tends to be much smaller than the dark matter dipole. The projected limits on the γm dipole are also shown in the right panel of Fig. 9; given that physical models in our framework have a direct detection rate mostly driven by γm dipole interactions to a first approximation the displayed limits can be compared to the direct detection curves shown in this plot is the case c (χ) CR,γ ] 2 = 0, reinforcing the picture just described.

C. Dark photon-mediated processes
The same procedure outline above can be applied to compute the recoil rate in case of processes that are dark photon-mediated (in the following: γ D m); we start with the effective SM quark-dark lepton interaction: (4.25) Borrowing the terminology from the previous section, we identify the first and second terms in (4.25) as γ D m dipole and charge-radius (CR) interactions, respectively. We then map the quark vector and tensor currents to nucleonic operators. The non-relativistic reduction of the effective DM-nucleon interaction yields where O 1 , O 4 , and O 6 are defined in Eq. (4.18), and Using the numerical values of QCD matrix elements obtained from lattice calculations [72], the coefficients in Eq. (4.26) are in the form with f T = 0.59 ± 0.023, and angle brackets denoting weighted averages that can be safely removed if γ D m dipole and CR coefficients are about the same for all light quarks. Analogously to the previous case, we organized the terms in Eq. (4.26) in powers of | q |, with the first line of order 1/| q | and the second of order | q | 0 . In the non-relativistic reduction, the γ D m dipole interaction has led to: (i) a long-range, incoherent O 3 term, (ii) a contact, coherent O 1 term, and (iii) other short-range, incoherent terms; the γ D m CR interaction has generated only one leading operator corresponding to a contact, coherent O 1 term. A summary with relativistic operators and the corresponding non-relativistic reductions is given in Table II.
Similarly to what has been done above for the γm case, we consider first the γ D m dipole and CR couplings to a quark as two independent coefficients, without any reference to our scheme. In the left panel of Fig. 10, 90% confidence level exclusion curve from XENON1T data and projected sensitivity curves for DARWIN are shown in the plane CR,γ D ] 2 for a few sample values of the dark matter mass: m χ = 200 GeV, 1 TeV, 2 TeV. In the right panel they are shown instead versus mass, assuming that only one among the two coefficients are different from zero. The solid diagonal line in the left panel marks again the separation between FIG. 10: 90% CL exclusion curves from XENON1T data and projected sensitivity curves for DARWIN in case of γD-mediated DM-nucleus scatterings, in the plane dipole-CR coefficients (left panel) or when assuming that only one of the two coefficients is non-zero (right panel, solid lines and the vertical scale on the left side refer to the γDm dipole operator, while dashed lines and the vertical scale on the right side refer to the γDm CR operator). Model-independent supernova cooling limit on the γDm dipole for quarks are displayed as horizontal solid lines for two representative values of αD: 10 −2 (black) and 5 × 10 −2 (green). Also shown in the left panel are two results specific for our dark sector framework: vertical dashed lines represent the projection of the supernova limit on the γDm CR coefficient -the intersection points with horizontal lines should be compared against direct detection results; the coloured regions correspond to two representative scans in the model parameter space, see the text for details.
the dipole-dominated region and CR-dominated region, as we can check looking at the expressions for the differential recoil rate. As in the previous case, CR interactions contributes with the coherent term in the form where A is the atomic number of the target nucleus. On the other hand and contrary to the previous case, for γ D m dipole interaction we cannot a priori assume that the long-range effects dominate: given that the long-range O 3 term is incoherent, we need to keep also the short-range coherent O 1 term, getting where S N is the spin operator for the valence nucleon (which is usually relevant for odd-even nuclei) and L is the angular momentum associated with the internal motion of the valence nucleon. Among the three contributions on the right-hand-side, although the first has a m N /E R enhancement, this has to compete with the large A 3 and 1/v 2 enhancements in the third term; moreover, the second term is most often sub-leading compared to the third given that ( L · S N ) 2 ≈ l 2 max , where l max is the maximum angular quantum number attained by the valence nucleon, typically much less than A. Comparing first and third contributions, one finds that the long-range 1/E R enhancement takes over only at recoil energies i.e. in a range which is irrelevant for a Xe target (as well as any target presently considered) and current detector technologies. Hence the third term is the leading one, and when taking the ratio between the rate in Eq. (4.31) and that in Eq. (4.30) one finds (4.33) It follows that: (4.34) which is the delimiter shown as a solid diagonal line in the left panel of Fig. 10. Turning now to constraints competing with direct detection results, contrary to the γm case, there is a strong model-independent bound impacting directly on the first operator in Eq. (4.25). In fact, the γ D m dipole for quarks can be responsible for enhancing the cooling rate in supernovae, allowing for nucleon-nucleon Bremsstrahlung emission of dark photons; as discussed in Section II A, there is a tight constraint one can extrapolate from the observed neutrino flux from SN1987A. The detailed derivation of the limit is rather involved and beyond the scope of this paper; we consider instead an extrapolation from analogous scenarios. In particular, Raffelt [42] computed the energy loss rate due to nucleon-nucleon Bremsstrahlung with axion emission, with the axion entering through a derivative coupling with the nucleon axial current. More recently an improved calculation has been implemented in [74]. For the case of γ D emission, Dobrescu [40] assumed that the rate of energy loss is two times larger than in the case of axion emission, given that the dark photon has two propagating degrees of freedom. If one writes the effective nucleon-γ D interaction as following Raffelt and imposing that the extra energy loss rate per unit mass induced by the novel Bremsstrahlung process cannot exceed 10 19 erg g −1 s −1 , we find: Here, f is a fudge factor accounting for the deviation from the Dobrescu assumption on the cooling rate when actually using (4.1) (in the following we will just set it to 1). Mapping the quantity g N γ D to the quark dipole moment in (4.1), and then mapping to the γ D m dipole coefficient constrained by direct detection, we have This limit is shown with horizontal solid lines in the left and right panels of Fig. 10, for α D = 10 −2 (black) and α D = 5 × 10 −2 (green). As it can be seen, at face value, the supernova limit is constraining γ D m dipole of quarks at a comparable level with respect to current direct detection data, while, regarding future sensitivities, direct detection experiments are going to be more competitive. On the other hand, the validity of the supernova limit has been recently questioned [75] since it relies on a mainstream picture for the explosion mechanism of core-collapse supernovae which is still, to a large extent, not well-established; in alternative scenarios the limit in Eq. (4.36) simply does not apply.
In this respect, information on the γ D m dipole of quarks derived from direct detection searches seem more reliable. While the γ D m CR operator does not contribute the dark photon emission via nucleon-nucleon Bremsstrahlung, a constraint can be indirectly derived within our framework implementing In the left panel of Fig. 10, limits on the γ D m CR coefficient, as derived from the supernova limit on the γ D m dipole, are shown with dashed vertical lines; these projections are obtained assuming m Q = 10 GeV, m S− = 10 TeV, and m S+ = 10 3 TeV. Note that vertical and horizontal lines cross in the dipole-dominated regime, hence the relevant comparison with direct detection rates is still in the limit of vanishing CR coefficient. Dipole dominance is typical for the parameter space in our scheme. In the left panel of Fig. 10 we show the regions in the dipole-CR plane corresponding to a scan with α L = 10 −1 , m Q ∈ [10 GeV, 50 GeV], m S− ∈ 10 TeV, 10 3 TeV , m S+ = 1.001 × 10 3 TeV, and either α D = 10 −2 (grey region) or α D = 5 × 10 −2 (green region); in scanning the model space, we ensured that m Q ≤ m S− ≤ m S+ . Most models are in the dipole-dominated area, with only tails extending into the CR-dominated regime in case the γ D m dipole gets severely suppressed when S − and S + are very close in mass and hence the mixing η s is very small. In Fig. 11, we plot recoil spectra in case of γm interactions (left panel) and γ D m interactions (right panel) for sample models in our dark sector framework. For each mediator, we have chosen two representative points such that "Pt. 1" lies in the corresponding dipole-dominated region, while "Pt. 2" in the CR-dominated regime: the corresponding model parameters are specified in Table III. Contributions to the differential rate of the dipole and CR operators are shown separately. Notice the qualitatively different shapes of the dipole contribution in the two cases: the 1/E R scaling due to long-range interactions can be seen in the γm case, while contact interactions dominate in the γ D m case. Notice also that Pt. 2 in the γm case is rather peculiar, since to find a model within the CR-dominated regime we were forced to consider a relatively small m χ , below the range considered for the scan displayed in Fig. 10 and what we expect typically in our framework.

D. Comparison with relic density limits
We are now ready to combine direct detection results with the constraints on our framework obtained by imposing that the relic density of χ matches the observed abundance of DM in the Universe, Ω DM h 2 = 0.1200 ± 0.0012 [47]. We refer to our minimal 6-parameter setup, slicing the parameter space along the m χ − α L plane for reference values of the dark photon coupling α D , of the common scalar messenger mass parameter m φ = m S and mixing η s , and of the mass m Q for light quark-like dark fermions. In Fig. 12, along the curves labelled "relic" the dark matter relic density matches the observed dark matter density. In the "south-east" direction, i.e. towards larger m χ and smaller α L , the χ relic density exceeds the observed dark matter density, assuming that α D is fixed. In the opposite direction, the χ relic density is a fraction of the observed dark matter density. These portions of the parameter space could be, in principle, recovered referring to, e.g., non-thermal production of dark matter or non standard cosmological frameworks, see, e.g., [76,77]).
In the top-left panel a maximal scalar mixing η s = 1 − (m χ /m φ ) 2 has been considered, while in the top-left panel it is tuned to zero; results for two representative values of α D are displayed, namely 10 −2 (solid lines) and 5 × 10 −2 (dashed lines), while the other parameters are fixed to m Q = 10 GeV and m S = 10 TeV. Each isolevel curve for Ω χ h 2 exhibits the features described by the decoupling regimes for χ discussed in Sec. III C. The upper branch corresponds to region IV, where the χχ annihilation to SM leptons controls the final relic density of χ. The vertical branch corresponds to region III where the annihilation to γ D determines the final relic density of χ: note that, in order to have the same relic density, increasing α D requires increasing m χ as well, which is consistent with the expectation from Eq. (3.2). The remaining branch corresponds to region II, where the final relic density is still determined by the γ D channel, but in the relic density regime given by Eq. (3.35), where a larger α L leads to larger ξ f.o. and thus m χ must decrease accordingly (since σv γ D goes as m −2 χ ). The branch for region I is not shown in these plots, but would simply correspond to a vertical line at lower values of α L .
Concerning direct detection limits and projected sensitivities, in the top-left panel of Fig. 12, the XENON1T and DARWIN curves (solid lines corresponding again to α D = 10 −2 , while dashed lines to 5 × 10 −2 ) are driven by the γ D m dipole operator, given that in the large mixing scenario this gives a larger event rate than the γ D m CR operator: we find that a large part of the upper branches with correct value of the relic density is already excluded by current direct detection limits, while a larger portion of region III will be tested with DARWIN. On the other hand, in the top-right panel η s = 0 suppresses the role of the γ D m dipole operator and the γ D m CR operator provides instead the bulk of the direct detection events: while current esperiments do not test this regime, DARWIN will be able to probe the branch with correct relic density in the region IV snd a portion of the one in region III. Note however that these results depend to some extent on the assumption of universality in the scalar messenger sector: the displayed direct detection curves would shift to larger values of α L in case some hierarchy between φ and S is assumed, with a larger m S relaxing the direct detection limits, without any significant impact on the result for the relic density of χ, given that the S scalars only interact with SM quarks.
A tighter connection appears with limits and projected sensitivities when including γm interactions; having artificially switching off γ D m interactions, in the bottom panel of Fig. 12 we show the curves stemming from the γm dipole coupling, together again with the relic density isolevel curves in case of both η s = 1 − (m χ /m φ ) 2 and η s = 0 (in this plot we are always the dipole-dominated region in Fig. 9, hence the γm CR coupling plays a minor role). Such results look less constraining than for γ D m interactions, however they do not depend on m S or η s , and hence can be more solidly compared against the relic density lines: we find that XENON1T data are in fact excluding part of the upper branch at small η s which was not tested via the γ D m operators, as well as that DARWIN will have a sizable impact in probing our scenario.

V. SUMMARY AND CONCLUSIONS
It is plausible that the solution to the dark matter problem may be in a context in which, on top of one or more particles accounting for dark matter, there are several extra states and/or extra forces. In this work, we have considered a toy model realization of a multicomponent dark sector, with an additional unbroken U(1) gauge interaction, mediated by a massless dark photon, and with portal interactions between dark fermions and SM fermions through scalar messengers. The model is characterized by: (i) the dark U(1) coupling α D , (ii) the Yukawa-like portal couplings α L,R , and (iii) the masses of the scalar messengers and the dark fermions. Despite its simplicity, this model has a rich dynamics and several phenomenological consequences. Its stable relics can provide a significant additional radiation component, as well as match the measured dark matter density in the Universe, with dark components having sizable interactions with ordinary matter as well as non-negligible self-interactions.
To characterize these features and have reliable estimates of final particle densities and temperatures, we have introduced a properly extended system of coupled Boltzmann equations, which track simultaneously the number density of several particle species, as well as entropy and energy exchanges between the dark and visible sectors. We have solved it numerically, implementing a few procedures allowing for fast -but very accurate -solutions. The target is to have a lepton-like dark fermion χ as the dark matter candidate, a requirement which selects viable regions in the model parameter space, without however singling out a definite scheme for the dark matter generation in the early Universe. In fact, depending on the strength of the Yukawa portal between visible and dark sector, we have identified four different regimes for the χ production, ranging from the limit of a WIMP-like scenario in a totally decoupled dark sector, to a FIMP-like generation in case of intermediate coupling, and up to a standard WIMP framework when the two sectors come and stay in kinetic equilibrium all the way through the chemical decoupling of all dark sector species. As a consequence, our framework is not very predictive regarding the mass scale of the dark matter candidate, which we can only point to be rather heavy, in the range, say, 500 GeV-10 TeV, with portal couplings all the way from about α L 10 −2 -1, down to around 10 −9 -10 −7 .
The result on the relic density for lepton-like dark fermion χ is weakly dependent on the choice for the masses of the lighter quark-like dark fermions Q. The latter can have a negligible contribution to the Universe matter density, say FIG. 12: The lines labelled "relic" correspond to model parameters αL and mχ for which the relic density of χ matches the abundance of dark matter in the Universe. In the top-left panel the case of maximal mixing for scalar messengers is considered, while in the top-right a case with ηs = 0 is displayed; in these two panels solid lines refer to the choice of αD = 10 −2 and dashed lines to αD = 5 × 10 −2 , while the other parameters the model are fixed to sample values, see the text for details. Also displayed in two top panels are XENON1T limits and DARWIN sensitivity curves due to γDm interactions, mainly due to the dipole operator in case of large mixing and the CR operator in case of zero mixing; solid and dashed lines refer again to the two sample values of αD. In the bottom panel the four relic density isolevel curves from the top panels are reproduced to be shown against XENON1T and DARWIN results in case γm interactions are included while γDm interactions are switched off (by e.g. raising the mass scale for scalar messengers in the quark sector); the dipole term is dominant and results do not depend on ηs or αD.
below 1% with respect to the heavy dark lepton contribution, if there is a sizable mass splitting between χ and Q, say m Q 100 GeV for m χ 1 TeV. On the other hand, the presence of light dark fermions enters critically in setting the temperature ratio ξ between dark photons and SM photons at the kinetic decoupling between the two sectors; this is one of the most critical observables in our model, since the CMB constraint on the amount of extra radiation in the Universe (usually given in terms of the effective number of neutrino-like species N ef f ) limits ξ to be at most 0.6 (at the 3-σ level). For a given portal coupling, constraints on the number and masses of light quark-like dark fermions follow: E.g., in a scenario with two light dark quarks (N Q = 2) and the early-time temperature ratio initialized to ξ 0 = 0.1, the limit on ξ CMB is satisfied at 1-σ level for any dark quark mass m Q if α L 10 −7 ; for α L 10 −3 , m Q lighter than 10 GeV (2 GeV) are excluded at 1-σ (at 2-σ); if α L 10 −1 , Q lighter than about 0.5 GeV are excluded at more than 3-σ.
Regarding other constraints on our scenario, we checked its testability with direct detection searches. The elastic scattering of the dark matter candidate χ on a nucleus can be mainly driven by dipole (dimension 5) or charge radius (dimension 6) interactions mediated by either the SM photon or the dark photon. We have analyzed on general grounds the interplay among the different operators, discussing features in the recoil spectrum and enlightening that long-range effects are not always predominant (as usually assumed in this context). After deriving current limits and projected sensitivities for next-generation detectors in terms of generic dipole and charge radius couplings, we have applied the results to our specific toy model, showing, e.g., that the DARWIN experiment will cover a significant portion of the parameter space in which χ is a viable dark matter candidate, as well as it will be competitive against the tightest (but model-dependent) constraints at present, including extra contributions to the magnetic dipole moments of leptons and extra cooling of stellar systems.
This exploratory work on a particular realization of a multicomponent dark sector model, can be extended further by investigating more general early-time initial conditions as well as a further extension of the particle content or more general particle interactions. Furthermore aspects are also not discussed here, such as its mapping on precision and accelerator physics, or further cosmological and astrophysical implications, including, e.g., the level of dark matter self-interactions. Some of this directions will be investigated in future work.
In writing the collision term in the right-hand side of the Boltzmann equations in section III, we need the amplitude squared of the relevant annihilation and elastic scattering amplitudes. Regarding the annihilation processes, the dark fermions can annihilate to dark photons or SM fermions (see Fig. 1). The squared amplitudes in case of χ (the expressions for Q are specular) are given by where s, t and u are the standard Mandelstam variables. When computing the pair annihilation cross section of dark fermions we need to include the Sommerfeld enhancement induced by the long-range attractive force mediated by dark photons [78] (the importance of this non-perturbative effect in the context of dark matter annihilations was first pointed out by [79]); for such Coulomb term, the enhancement can be computed analytically and added as a multiplicative factor to the cross section σ 0 accounting for contact interactions where v is the velocity of each annihilating species in the center-of-mass frame. As for the elastic scattering processes, the dark fermions can either undergo Compton-like processes with dark photons, or scatter on SM fermions (see Fig. 3). The squared amplitudes in case of χ are In the Boltzmann code developed in Section III there are several quantities involving thermal averages. Starting with pair annihilation cross sections, a method to efficiently compute σv (T ), as defined in Eq. (3.17), was detailed in [53]: Assuming equilibrium distribution functions with occupation numbers approximated by the exponential in (3.13), you can manipulate the numerator by performing a change of integration variables from the two momenta p 1 and p 2 to E + ≡ E 1 + E 2 , E − ≡ E 1 − E 2 and s, with the integral in the first two that can be performed analytically, giving where K 1 (z) and K 2 (z) are the modified Bessel functions of order 1 and 2, respectively. The same method can be applied to σvE (T ), see the definition in Eq. (3.28), obtaining Both these expressions are very convenient when coming to their numerical implementation: for any given particle physics model, you can first tabulate the cross sections σ as a function of s, and then link to such tabulations for a fast computation of thermal averages at any T in the temperature evolution equations. On the other hand, an analogous shortcut cannot be implemented in thermal averages for momentum transfer rates. Referring generically to the scattering process i + B → i + B, the thermally averaged momentum transfer rate γ iB (T i , T B ) is a function of the temperature of both the species i and bath particles B, see Eqs. (3.26) and (3.24), and such dependences cannot be simply factorised, making the implementation in the numerical Boltzmann code CPU-demanding. When investigating the kinetic decoupling of massive dark matter particles i from the heat bath, since this typically occurs in the regime at which the temperature T i is small compared to the particle mass m i , [56] and [80] noticed that the dependence on the particle momentum p i in γ iB can be approximately dropped, thereby allowing to replace γ iB (T i , T B ) with γ iB (E i = m i , T B ). When the scattering species is relativistic, this is not a fair estimate; on the other hand, we can still use it as a guideline for a better approximation: At T i m i the occupation number in the integrand at the numerator of the l.h.s. of Eq. (3.26) is sharply peaked at the stationary point E i = m i and γ iB (E i , T B ) simply picks up the contribution coming from the stationary point. On the other hand, when T i m i , the occupation number has a relatively longer tail at higher energies. The pre-factor (E 2 i − m 2 i ) 3/2 in the integral cannot be neglected, and, to extract the peak contribution, one has to search for the stationary point of the function which is now at Note that going back to the limit T i m i , you correctly retrieve E * = m i +O(T i ). The agreement between γ (T i , T B ) and γ(E i = E * (T i ), T B ) is very good, as shown in a sample case in Fig. (13).

Appendix C: Loop calculations
We report here a few details regarding the computation of the γ and γ D vertex functions represented by the loop diagrams in Fig. 8. For the γ D vertex function, involving a SM quark q on the external legs and a quark-like dark fermion Q (with U (1) D charge Q Q ) and scalar messengers S ± in the loop, we havē q(k ) iΓ µ γ D q(k) = g D Q Qq (k ) λ=± g 2 L + g 2 R 2 I µ a,γ D (m Q , m S λ , k, q) + (λ) g L g R m Q I µ b,γ D (m Q , m S λ , k, q) q(k) , (C1) where q µ is the momentum transfer. For the γ vertex function, with a lepton-like dark fermion χ on the external legs and the corresponding lepton (with U (1) em charge Q l ) and messengers scalars φ ± in the loop, we havē χ(k ) iΓ µ γ χ(k) = eQ lχ (k ) λ=± g 2 L + g 2 R 2 I µ a,γ (m l , m φ λ , k, q) + (λ) g L g R m l I µ b,γ (m l , m φ λ , k, q) χ(k) .
where we introduced the function D(p, m) ≡ p 2 − m 2 , while s γ = 1 and s γ D = −1. The s V sign structure is motivated by the form of the interaction Lagrangian in Eqs. (2.1) and (2.2): a dark fermion and its corresponding messenger scalar must have opposite U (1) D charges, while a SM fermion and its corresponding messenger scalar must have the same U (1) em charge. The additional contribution to the magnetic dipole moment of SM leptons predicted in our dark sector framework is computed from the γ vertex function having SM leptons as external legs and a loop with χ and φ ± . We find l(k ) iΓ µ γ l(k) = eQ ll (k ) λ=± g 2 L + g 2 R 2 J µ a (m χ , m φ λ , k, q) + (λ) g L g R m χ J µ b (m χ , m φ λ , k, q) l(k),