Multi-component Dark Sectors: Symmetries, Asymmetries and Conversions

We study the relic abundance of several stable particles from a generic dark sector, including the possible presence of dark asymmetries. After discussing the different possibilities for stabilising multi-component dark matter, we analyse the final relic abundance of the symmetric and asymmetric dark matter components, paying special attention to the role of the unavoidable conversions between dark matter states. We find an exponential dependence of the asymmetries of the heavier components on annihilations and conversions. We conclude that having similar symmetric and asymmetric components is a natural outcome in many scenarios of multi-component dark matter. This has novel phenomenological implications, which we briefly discuss.


Introduction
Several observations indicate the existence of non-visible matter termed as Dark Matter (DM), that comprises roughly 25% of the total energy budget of the universe, and is therefore almost five times more abundant than the visible matter [1]. Despite this significant contribution, hardly anything is known about its composition, nature or dynamics. Among existing efforts and speculations made regarding the particle nature of DM, Weakly Interacting Massive Particles (WIMPs) are one of the simplest and most popular particle DM candidates, contributing to the energy density once they undergo thermal freeze-out. On the other hand, Asymmetric Dark Matter (ADM) models are motivated by the close relation between DM and baryonic energy densities, with ρ DM ∼ 5ρ B [2]. Given that the baryon density is set by the baryon asymmetry of the universe, η B ∼ 6 × 10 −10 , the DM abundance could stem from an asymmetry in the dark sector in a similar manner (see Refs. [3,4] for a review). In the case of WIMPs, the value of the abundance is set by the thermally averaged annihilation cross section, whereas for ADM it is set by the initial asymmetry and the DM mass. Moreover, as discussed in Ref. [5], there is also the possibility of an intermediate regime between the two extremes discussed above, where the DM is partially asymmetric and there is an exponential dependence on the annihilation cross section.
However, most of these scenarios assume that the entire dark sector is composed of just one DM state. Given the plethora of stable particles in the Standard Model (SM) (electrons, protons, nuclei, neutrinos, photons...), it is reasonable to think that the dark sector could be made up of more than one stable component. Efforts are being made towards the study of such multi-component DM models, where the dark sector is comprised by several stable particles. For example, in Refs. [6][7][8] the influence of DM conversions (annihilation of one DM state into another) in the final relic abundance was emphasised. Similarly, semiannihilations can affect the final abundance [9]. Finally, multi-component DM, with its many degrees of freedom, leads to a richer phenomenology (see for instance Ref. [10][11][12]), and can help to relax the existing bounds from direct detection, indirect detection and collider searches. Another proposal for multi-component dark sectors is Dynamical Dark Matter [13]. In this framework, the dark states come in large numbers and are unstable, with lifetimes proportional to their individual abundance. For a recent review of their possible collective effects see Ref. [14].
The numerous states and interactions with SM and among themselves in a multicomponent DM setup allow us to treat it as a complex system, which exhibits a behaviour that can be characterised by different scalings: power law and exponential. Such scaling laws are common in Physics and in other branches of science. Understanding this behaviour of multi-component DM subject to different interactions (annihilations and conversions) and conditions (symmetric/asymmetric) is one of the main goals of this paper. For that purpose, we will extend the analysis performed for one component asymmetric WIMP in Ref. [5] to two or more components with conversions among them.
The paper is structured as follows. In Section 2, we discuss the symmetries and combinations that are required to stabilise multiple states in the dark sector. In Section 3, we briefly discuss how DM asymmetries in different components can be generated from a B −L asymmetry. In Section 4, we develop the multi-component framework for both symmetric and asymmetric components. First, we study the final relic abundance of multi-component DM taking into account the unavoidable DM conversions. Then, in Section 5, we parameterise the possible particle physics models using effective interactions. We discuss the phenomenological implications of these scenarios in Section 6 and provide some conclusions in Section 7.

Dark Sector Stabilisation Symmetries
In this section, we briefly discuss the simplest possible symmetries to have multi-component DM. In DM models, usually a discrete and/or continuous global symmetry is imposed to stabilise the DM candidate. As pointed in Ref. [15], continuous gauge symmetries alone cannot stabilise the DM and need to be spontaneously broken. However, they can motivate the origin of discrete gauge symmetries as a remnant of spontaneously broken gauge symmetries [16], such as matter parity (and R-parity in supersymmetry) from gauged B −L [17]. In this work, we consider a multi-component DM framework that can be stabilized by adopting different symmetries, with the presence of both symmetric and asymmetric components. For symmetric DM, it is more natural to have a Majorana fermion or real scalar, stabilised by discrete symmetries. On the other hand, for asymmetric DM, the candidate needs to be a Dirac fermion or a complex scalar, stabilised by continuous symmetries.
Let us start with the case of discrete symmetries. Most popular DM models employ a single Z 2 symmetry which guarantees the stability of the lightest Z 2 odd particle, allowing only one DM candidate in the model. In order to have multi-component DM, we would require a larger discrete symmetry or consider products of more than one discrete symmetry. We first investigate the number of stable DM candidates that are allowed for Z N >2 , in particular for N = 3, 4, and then for the product of two or more discrete symmetries.
Under an abelian Z N symmetry, a field transforms as where q = 0, 1, 2, . . . , N − 1 are the Z N charges. Depending on how the fields transform under a given charge, we can divide them into classes. Class 0 corresponding to q = 0 (implying that the field is neutral under the given symmetry) is by default common to all Z N symmetries. For example, consider the case of Z 2 , where q = 0, 1 and the field transforms as φ → φ (even) and φ → −φ (odd) respectively. Thus, there are only two classes for Z 2 : {0, 1}. The lightest Z 2 even particle will be from the SM and the odd, the DM candidate.
In the case of Z 3 , we can assign charges q = 0, 1, 2 to the fields. However, the charge q = 2 corresponding to ω 2 is same as assigning q = −1 since ω 2 = ω * , implying that the particle with charge q = 2 is nothing but the antiparticle of the particle with q = 1. Therefore, in the case of Z 3 , we also have only two classes {0, 1} and there is only one DM candidate, the lightest particle transforming under a non-trivial Z 3 charge [18].
For Z 4 , ω q = 1, i, −1, −i for q = 0, 1, 2, 3 respectively. Here, note that q = 3 is equivalent to q = −1, and they would thus be placed in the same class. So, the classes for Z 4 are {0, 1, 2}. The lightest particle of class 1 (i, −i), let's say φ 1 is always stable and a DM candidate. The lightest particle of class 2 (−1), say φ 2 can also be stable, provided m φ 2 < 2m φ 1 . Thus, Z 4 guarantees the stability of at least one DM candidate, but can also allow for a second DM candidate in case its decays are kinematically forbidden to the particles belonging to class 1 [19]. For DM models with Z 3 and Z 4 symmetries, see Refs. [9,[20][21][22].
In general, if N is a prime number, the lightest particle transforming with a non-trivial Z N charge is always stable, and a potential DM candidate. One can also have heavier particles with Z N charges as DM if their decay modes are kinematically forbidden. If N is a composite number, such that N = p s 1 1 p s 2 2 . . . p s k k , where p i is a prime number (p i = p j for i = j) and s i is natural, Z N can be decomposed as a product of smaller groups [23,24] and there can be at most s 1 + s 2 + . . . + s k stable particles depending on the kinematic conditions. Therefore, in order to realize multi-component DM with just one imposed discrete symmetry Z N , we require N ≥ 4. For studies of multi-component scalar DM in larger Z N groups, see Ref. [25]. Now we turn our discussion to the product of two or more of the same discrete symmetries. The simplest and most commonly used symmetry is Z 2 × Z 2 [7,8,11,[26][27][28]. Such a choice of symmetry allows the existence of at least two stable particles and at most three, subject to kinematics. For the case of an n Z 2 product: Z 2 × Z 2 × · · · × Z (n) 2 , there are at least n stable DM particles and at most 2 n − 1 stable particles can be DM candidates. It is easy to generalise this for the case of Z 3 and Z 4 where c denotes the number of classes corresponding to the discrete symmetry Z N . For N = 3, since there are only two classes, the number of DM candidates for Z 3 × Z 3 is equal to that of Z 2 × Z 2 (see Ref. [29] for a model with a Z 3 × Z 3 symmetry). On the other hand, in the case of Z 4 which has three classes, Z 4 × Z 4 allows for at max 3 2 − 1 = 8 DM candidates. The '−1' term guarantees stability of the lightest particle in the SM. We can also consider products of different discrete symmetries, for example Z 2 × Z 3 (see Ref. [30]), and generalise Eq. (2.3) as 4) where N i corresponds to the discrete symmetry used and n i denotes its multiplicity in the product.
In the following, we will also be interested in the case of multi-component DM with states that can be asymmetric. In this case, the most natural framework is the stabilisation of the DM components due to continuous global symmetries, which may be accidental. The possibilities we will implicitly consider are: i) a single global continuous U (1) symmetry: with the DM χ i having charge q i such that no linear terms can be written in χ i at least at the renormalizable level. There could be linear terms at the non-renormalizable level, possibly suppressed by Planck-scale physics and therefore sufficiently suppressed [31].
(ii) the direct product of several U (1)s, such that each particle only transforms under one of them and is also the lightest particle charged under it, with DM i transforming with the global charge χ i ∼ (0, . . . , i = 1, . . . , 0). For multicomponent DM models stabilised by the product of global symmetries, see Refs. [6,7].

Asymmetries in the Dark Sector
In general, DM may have both symmetric and asymmetric components. One of the simplest options is that a dark asymmetry can be generated from an asymmetry in the visible sector [2]. The idea is that an asymmetry in the SM sector is transferred to the DM sector charges generate a mixed operator, which is active at high energies and transfers any initial B − L asymmetry in the SM to the DM sector. For definiteness, let us consider the case of Dirac DM, χ i . We define n as the number densities of the χ i particle (anti-particle), with n + i > n − i without loss of generality, and similarly for baryons. The baryon number density, [32]. The DM number density is X i ≡ n + i − n − i , which can also be expressed in terms of B − L as X i ≡ x i (B − L). It is useful to express the number densities in terms of the yields Y i ≡ n i /s, where s is the entropy density given by s(T ) = (2π 2 /45) h eff (T )T 3 and h eff (T ) are the effective degrees of freedom related to entropy. We also define η Depending on whether the interaction rate generated by the transfer operator, Eq. (3.1), is active below or above the EW scale, one obtains [32,33]: If there are strong dark interactions (large annihilations), the symmetric component may be totally erased. In that case, in order to reproduce the DM abundance, we require If this is the case, there is a prediction for the combination of DM masses. It only depends on the energy at which the operator is active, and on how the asymmetry is transferred, i.e., on the B − L charge. For example, in the case of 2DM there are different possibilities: 1. Identical DM particles with same charges, e.g. B − L → X 1 , X 2 = X 1 : the same operator couples to both DM particles. For instance, if the operator is active above (below) the EW scale, one obtains In the case of N components, we get Therefore, the larger the number of asymmetric states, the lower their masses. Indeed, if all masses are similar, the presence of N -asymmetric components would translate into an N -suppressed mass-scale compared to the 1-DM scenario, Notice that this can be understood as an upper limit on the mass scale in order for the DM not to be overabundant. Of course, as for 1DM, heavier DM is allowed if the DM asymmetry is unrelated or suppressed compared to the baryonic one.
2. Different DM particles with different charges, e.g., B − L → X 1 , X 2 = X 1 : different operators couple to both DM particles. For example, if operator one (two) with B − L charge Q 1 (Q 2 ) is active above (below) the EW scale, one gets In the case of N = (N 1 +N 2 ) DM components, where N 1 (N 2 ) components are coupled to operators active above (below) the EW scale, the above equation generalizes to Note how Eq. 3.4 (Eq. 3.5) is a particular case of Eq. 3.7 (Eq. 3.8) for identical charges above (below) the EW scale.
3. Different DM particles with different charges, B − L → X 1 → X 2 . Only one of them, χ 1 , is coupled to the SM; the asymmetry is first transferred to one DM state and then there is a further sharing of the asymmetry with another DM state. For example, this can be generated by the operators with φ being a scalar DM particle. This implies x 2 = 2 x 1 and therefore 1 m 1 6.5 GeV + m 2 3.25 GeV = 2 . (3.10) In the next section we consider the general case where both asymmetric and symmetric components may be present. If the asymmetric component dominates, the previous mass relations are very good approximations.

Multi-component Framework
In the following, we set up the framework for multi-component DM with both symmetric and asymmetric components. Let us consider an N -component DM scenario stabilised by a certain choice of continuous U (1) symmetries, as discussed in Section 2, with m 1 > m 2 > . . . > m N > m SM . Within this setup we may have one or more of the following reactions: 1. Annihilations χ iχi → SM SM. We use the notation σ χ iχi →SM SM ≡ σ ann, i .
4. Semi-annihilations χ i χ j → χ k SM, for N > 2. These processes depend on the symmetry imposed. In this case, DM needs to carry some SM charge.
In the following, we consider only the case of annihilations and conversions. We are interested in conversions of the form χ iχi ↔ χ jχj , which are not forbidden by any symmetry that allows co/annihilations and are therefore typically always present. This is the case if both χ i and χ j are charged under Z 2 and Z 4 (class 2) symmetries respectively. Other conversions such as χ i χ i χ i ↔ χ j χ j , for example, are phase-space suppressed and allowed only if χ i is charged under Z 3 , and χ j under Z 2 . Similarly, conversions of the form χ i χ i χ j → χ iχjχj are also phase-space suppressed.
Furthermore, we always have annihilations that are significant and comparable to conversions. This avoids significant reheating in the dark sector if DM states become nonrelativistic and their temperatures increase exponentially [39]. The non-negligible coupling to the SM also ensures that they are thermalised with the SM plasma at freeze-out. Therefore, considering only 2 → 2 processes, assuming kinetic equilibrium with the SM of all species and neglecting the quantum statistics for the species involved, the evolution equation for the number densities of the particles read wheren ± i = n eq e ±ξ denote equilibrium number densities of the species, for which we use Maxwell-Boltzmann distributions, n eq = g(mT /2π) 3/2 e −m/T , with µ = ξT the chemical potential and g the number of internal degrees of freedom. Throughout this work, we denote by σ ann, i v and σ conv, ij v the thermally-averaged annihilation and conversion cross sections, respectively. The Hubble parameter is given in the radiation-dominated epoch by 19 GeV, and g eff (T ) the effective relativistic degrees of freedom. 3 Furthermore, the temperature evolution of the Universe can be expressed by scaling it with the mass of the heaviest component (m 1 in our case) using x = m 1 /T , so that Eq. (4.1) can be written in terms of the yields as For ease of analysis and to highlight the qualitative aspects of the setup, we first consider a two-component DM system and discuss the cases where both components are symmetric or asymmetric. Later on, we extend the analysis to more complex sectors.

Symmetric components
For 2DM with symmetric components, Eq. (4.2) can be written as If conversions are neglected, the Boltzmann Equations in Eqs. (4.3) decouple, and in order to reproduce the observed relic abundance, we require (for m i ≥ O(GeV) and assuming s-wave annihilations), Therefore, the abundance is dominated by the particle with the smallest annihilation cross section, i.e., it is like 1DM, unless σ ann, 1 v 1 σ ann, 2 v 2 . If both particles have similar annihilation cross sections, the yield is dominated by the lightest particle, and the contributions to the relic abundance from both states are of similar size.
On the other hand, if conversions are significant, both components may end up with similar abundances, i.e., DM may be multi-component. For instance, assuming the hierarchy σ ann, 2 v σ conv v σ ann, 1 v , the abundance of χ 1 (χ 2 ) gets reduced (increased) with respect to the case without conversions [6-8, 11, 28]. This behaviour is illustrated in Fig. 1, where we show the evolution of abundances for χ 1 and χ 2 for the assumed hierarchy of cross sections in presence and absence of conversions. With conversions, Eq. (4.4) is roughly modified as It is clear from this that conversions reduce the total relic abundance, and so, in comparison to the single-DM-component case, smaller annihilations into SM states are necessary. This is turn allows to have smaller couplings between the visible and the dark sector, and therefore, to relax existing constraints. We can extend the above discussion to N components. Let us consider the case where there are N DM components that have similar cross-sections for annihilations and conversions, e.g. σ ann,i v = σ conv,ij = σv for all i, j. This simplifies the analysis so that the only parameters of interest are the individual DM masses m i and σv . For example, let us consider the case of N = 10 and take m i − m i+1 = 100 GeV, with m 1 = 1000 GeV.
In Fig. 2, we plot the individual asymptotic yields (left panel) and total abundance (right panel), as well as the fits to power laws of DM mass. Since the cross section is equal for all components, the yield goes as ∼ m −1 for the case without conversions and ∼ m −2 for the case with conversions. The inverse dependence on mass is well known from the analytic expression for the yield in the standard freeze-out case without conversions, as all the Boltzmann Equations are decoupled. The additional mass suppression in the case of conversions can be understood as follows: the heavier the mass of the DM component, the larger the number of channels for conversions. For example, χ 1 , being the heaviest component, can undergo conversions to all the other 9 components whereas for χ 9 , the only allowed conversion is χ 9 χ 9 → χ 10 χ 10 . This implies that the individual abundances are almost independent of mass in the case without conversions, whereas in presence of conversions, they are inversely proportional to mass, as Ω i ∝ m i Y i . For the lightest component, the yields are identical in both cases since there are no channels for its conversion and the excess due to its production from all the other components gets annihilated to SM. In Fig. 2 (right panel), we can see that the total abundance is inversely proportional to the cross section in both the cases. Furthermore, the presence of conversions allows us to reproduce the correct abundance for smaller values of the annihilation cross-section (roughly a factor of 2), and thus in general of the DM couplings to SM. This is a generic feature, which as argued before translates in a relaxation of the experimental bounds.
Finally, let us mention that for cold thermal relics, there exists a bound on how heavy the DM particle can be from unitarity in order not to overclose the universe: m χ 110 TeV [40,41]. In the case of multi-component DM, this translates into a bound on the masses of the DM components for s-wave annihilations, Therefore, as for the case of asymmetric DM, the existence of a complex dark sector points to lower DM masses. Again, if all masses are similar, the presence of N -symmetric components translates into an √ N -suppressed mass-scale compared to the 1-DM scenario, where in the last step we assume s-wave annihilations and conversions. Notice how, now, the suppression with N is milder as compared to the asymmetric case, c.f. Eq. 3.6. The above bounds may be relaxed to 190 TeV (220 TeV) for p-wave (for a combination of sand p-wave) annihilations and conversions [41].

Analytical expressions
Now we study how the picture changes if there are initial DM asymmetries. We follow closely the work of Ref. [5] and extend it to the case of 2DM with initial asymmetries Of course, the asymmetries remain constant and equal to the initial value throughout the evolution. Following their notation, we can define the asymmetric where the upper (lower) limit is for totally (a)symmetric DM. In this way, the co-moving number densities Y ± i can be completely expressed in terms of η i and r i as Given that the chemical potentials of particles and their antiparticles are equal and opposite, we definer i ≡ e −2ξ(x) with ξ = sinh −1 (η i /(2Ȳ i )). Hence, Eq. (4.2) can be written for the asymmetric ratios r i in the suggestive forms Notice that if the initial asymmetries η i are zero, the ratios r i would remain constant and equal to one, like in the standard freeze-out case.
It is interesting to derive some analytical solutions to better understand the behavior of the system. Considering the evolution of r i (i = 1, 2), we have two different regimes: where x f i denotes the freeze-out temperature for the species. We can expand the annihilation and conversion cross sections in the relative velocity and express the thermally averaged cross sections as where one can consider different combinations of s-wave (k i/c = 0) and p-wave (k i/c = 1) cross sections. 6 We can write where λ a,i = (π/45) 1/2 M Pl m 1 σ i , f c ≡ σ c /σ 1 , and In the following, for simplicity, we consider the case: k 1 = k 2 = k c ≡ k, η 1 = η 2 ≡ η, λ a,2 λ a,1 ≡ λ a and approximate ζ(r i ) ∼ 1. Therefore, with the above approximations, we obtain Given that m 2 < m 1 , for similar cross sections r 2 r 1 , and hence we can neglect f c r 1 in the second line of Eq. (4.11) and solve the decoupled system wherer Therefore, the DM asymmetry ratios decrease exponentially with the annihilation crosssection (i.e., DM becomes more asymmetric as annihilations are increased), in a similar way as for 1DM scenarios [5]. Furthermore, in presence of significant conversions, the asymmetry ratio for the heavier component also decreases exponentially with the conversion cross section. This is analogous to the symmetric case (see Eq. (4.5)) where the yield for the heavier components was inversely proportional to the sum of annihilation and conversion cross sections whereas that for the lighter one it was only inversely proportional to its annihilation cross section. In this case, however, the crucial difference is that there is an exponential dependence.

Numerical analysis
In Fig. 3, we plot the evolution of the asymmetric ratios and components after solving Eq. (4.8) and using Eq. (4.7). Note that Y + 1 ∼ Y + 2 ∼ η B , due to our choice of η 1 = η 2 = η B and large annihilation cross-sections, so that the symmetric component is significantly depleted. In general, one obtains Y + i ∼ η i . From the left plot in Fig. 3, it can be seen that the heavier component, χ 1 , is more asymmetric than the lighter one, χ 2 , without conversions. Including conversions significantly reduces the symmetric component of the heaviest state. This can be seen in the right plot where Hence, we observe that while the heavier component is mostly asymmetric, the lighter one is mostly symmetric. Notice that the analytic solutions derived in Sec. 4.2.1 are contrasted with the numerical solutions in Fig. 3 (dotted lines) and are found to be in agreement.  Here, unlike the symmetric case, one cannot fit the asymmetry ratios as power laws in the DM mass. For the case without conversions, we find r i ∼ a e −b m i is a good fit whereas when conversions are included, we can fit the individual asymmetric ratios as the product of a power law and an exponential: r i ∼ a m b i e −c m i , where a , b , c depend on the value of cross sections. Notice that the DM mass enters subtly via x f i = m i /T i in the expression for Φ. We find that for similar values of the cross sections, the asymmetry ratios for the lighter species are larger than that of the heavier ones. This implies that, for a given initial asymmetry η, the heavier species are more asymmetric than the lighter ones, and they become even more asymmetric in the presence of conversions, thus they can be considered to be fully asymmetric. From the right plot of Fig. 4, it can be seen that the individual abundances coincide for the heavier particles with and without conversions and scale with the DM mass m i . Since the heavier species can be considered to be fully asymmetric, their abundances are independent of conversion and annihilation cross sections. This is valid if the cross-sections are larger than a certain value. For example, in this case, annihilation cross-sections larger than O(0.1) pb lead to r < 0.1.
As we can see, the presence of conversions reduces the total abundance also in the partially-asymmetric case. For initial asymmetries close to the baryonic one, i.e., η i η B , we find that only in the case of purely asymmetric DM states with masses above ∼ 10 GeV, is the effect of conversions irrelevant. Furthermore, the direct proportionality of abundance with mass is a good approximation for even lower masses in the presence of conversions, which tend to make the species more asymmetric, as illustrated by the black dashed curve. We also see that the abundance is at or below current Planck values for DM masses around 5 GeV, as expected in asymmetric DM models (see Section 3) due to the initial asymmetry being equal to the baryonic one.

Total relic abundance
We have seen that it is natural that if, DM has different states, the individual symmetric and asymmetric components can contribute significantly to the relic abundance for reasonable values of the cross sections. We follow closely the notation in Ref. [5].
For the case with initial asymmetries, the total DM abundance will consist of both symmetric and asymmetric components and the DM is partially asymmetric, with an energy density given by where r ∞,i depends on Φ i (∞, m i ). The first term corresponds to the asymmetric components (present for r ∞,i → 0), and the second one to the symmetric ones (which scale as the inverse power of the sum of annihilation and conversion cross sections for r ∞,i → 1).
Reproducing the relic abundance yields a constraint between the DM masses, the asymmetries and the annihilation and conversion cross sections through this relation. Notice that in some cases the initial asymmetries may be similar, η i η, i.e., if they come from the same operator, or if conversions/annihilation are large, and we can factorize η from the first term.
In the completely asymmetric case, i.e. r ∞,i = 0, the total DM abundance just depends on the masses and initial asymmetries, (4.15) Figure 5. Relic abundance as the function of DM mass for asymmetric (solid), symmetric (dashed) and partially-asymmetric (dot-dashed) cases. The red and blue curves correspond to two different choices of σ i = 2 pb, 4 pb respectively. The asymmetric case is linearly dependent on mass and η, which is chosen to be equal to η B . The symmetric case is almost independent of mass and depends only on the values of σ i . The partially-asymmetric case becomes independent of σ i for large mass values, e.g., it starts exhibiting a completely asymmetric behaviour.
In the symmetric limit r ∞,i → 1, the expansion of the second term gives which is independent of the initial asymmetries and very mildly dependent on the DM mass via Φ i (∞, m i ). In presence of conversions, σ i has to be summed over the conversion cross sections as well, such that total cross-section becomes , one can also re-express the fractional asymmetry for heavier components in presence of conversions as [5] where λ i,t = (π/45) 1/2 M pl m 1 (σ i + σ c ), and σ WIMP and Φ WIMP correspond to a symmetric DM particle of the same mass, m i . It should be noted that the freeze-out temperature (x f i ) entering the expression for Φ in Eq. (4.13) is modified when conversions are included, hence we use Φ i,t . Naively, the relic abundance in the partially-asymmetric case (Eq. (4.14)) might be expected to be the sum of Eqs. (4.15) and (4.16); however, that's not the case as the symmetric component also depends on η [5]. Only for scenarios where the DM components are completely symmetric or asymmetric, the total abundance can be given by their sum for i = j. In Fig. 5, we plot the dependence of the relic abundance on masses 7 and cross sections to show the differences between the symmetric, asymmetric and partially-asymmetric cases. It can be seen that the curve corresponding to partially asymmetric is not equivalent to the sum of the other two, and starts tracking the asymmetric curve for larger masses. The mass at which it starts tracking the asymmetric case depends on the value of σ i and η i .

Effective Parameterisation
In the analysis for symmetric and asymmetric components performed above, we have considered all cross sections to be independent of DM mass. However, these are model dependent and in particular, also scale with the masses of the particles involved. For instance, considering only s-wave annihilations we can parameterise the DM annihilation cross section at freeze-out as where m i is the DM mass, λ i the coupling characterising the DM-SM interaction and Λ is the mass scale of the heavy mediator involved. We consider two limiting cases for our analysis, corresponding to light and heavy mediator respectively, with a i and b i being model-dependent numerical parameters. In a similar manner, we parameterise the conversion cross section as where λ ij is the interaction strength between the DM components with masses m i and m j . As before, we consider the limiting cases Hence, in order to describe the dynamics of N -component DM, we require N DM masses, N DM-SM couplings and N (N −1)/2 DM-DM couplings i.e. N (N +3)/2 parameters in total. In the following, we analyse the results in a model independent way, treating a, b, c, d as free parameters and fixing them to be equal for all i, j. We then provide an example model in Appendix A.

Light mediators
First, let us consider the case of light mediators, where the cross-sections are independent of the mediator mass, see first line of Eqs. (5.2) and (5.4), and depend on the new parameters a i and c ij , taking N = 10 and varying the DM masses from 300 GeV to 1200 GeV. In the absence of conversions (c ij = 0), all Boltzmann Equations are decoupled and we roughly obtain the individual yields (abundances) scaling as ∼ m i (∼ m 2 i ). This is expected since Y i scales as ∼ m −1 when we treat all cross sections to be similar and independent of mass. The asymptotic abundances are plotted as a function of mass in the left plot of Fig. 6, for different choices of a. In this case, the abundances follow a power-law of the form where κ is to be determined from the fits. Moreover, we find that the scenario becomes interesting in presence of conversions with an additional dependence on the parameter c (we take all c ij ≡ c). The total conversion cross sections of a given DM particle into the lighter states σ conv, i v = j>i σ conv, ij v are of similar size. This is due to the fact that while heavier DM particles have more DM states to convert into, their conversion cross sections are suppressed by their mass, whereas lighter states have less DM states to convert into, but their cross sections are much larger. These effects partially compensate each other. As a result of this, the power law dependence on mass for abundances varies depending on the strength of conversions with respect to annihilations, unlike the case of mass-independent cross sections, where we saw Ω i h 2 ∝ m −1 i in presence of conversions. We illustrate this in the right plot of Fig. 6. As can be seen, for c a, the abundances do not change much with the mass of the DM component. The abundances in this case can be expressed in the form where κ and δ are numbers determined from the fit. If the conversions are weak, i.e. c a, the first term in the denominator of Eq. (5.6) dominates and Ω i h 2 ∝ m 2 i . On the other hand, for significant conversions i.e. c ∼ a or c a, the second term in the denominator dominates and we find Ω i h 2 does not depend on m i . Therefore with conversions, Ω i h 2 ∼ κ m β with 0 β 2.
Note that for extremely large conversions, the fits are determined by ignoring the lightest and next-to-lightest DM particles as they do not follow the same mass-suppression as heavier ones. However, for moderate to weak conversions, we find that it is possible to perform a fit for all the DM components.

Heavy mediators
In the case of heavy mediators and/or light DM, see second line of Eqs.  The case of no conversions is shown in the left plot of Fig. 7 for different choices of the parameter b. We have taken the heaviest DM component mass to be 20 GeV and the lightest 2 GeV. We also fix the mediator scale to Λ = 100 GeV. The asymptotic relic abundances follow a power law of the following form where κ is determined from the numerical fits. In presence of conversions, the dependence on mass will vary depending on the strength of conversions with respect to annihilations, similar to the light mediators case. However, due to the direct proportionality to mass, the total conversion cross sections of each species are not similar unlike in the previous case, and scale as ∼ m 3 i . Therefore, we can express the abundances as where κ , δ are determined from the numerical fits. When conversions dominate, i.e. d b, we get Ω i h 2 ∝ m −3 i , whereas for weak conversions, we obtain roughly Ω i h 2 ∝ m −2 i . Hence, for the heavy mediators case, we have Ω i h 2 ∼ κ m γ with −3 γ −2, see the right plot of Fig. 7.

Asymmetric components
As discussed above, in the case of completely asymmetric dark matter, the relic abundance is independent of the annihilation and conversion cross sections, whereas there is a slight dependence on these cross sections in the case of partially-asymmetric dark matter.
From the analytic expressions derived above, we know that the asymmetry ratio is exponentially suppressed by the annihilation cross section and the initial asymmetry η i . In the following, we discuss the effect this has on the DM relic abundance for the light and heavy mediator cases.

Light mediators
Let us first consider the case without conversions. When the cross sections are inversely proportional to the square of the DM masses, the lighter species will have a larger cross section and their asymmetry ratio gets exponentially suppressed compared to the heavier ones. The trends vary depending on the value of in Eq. (3.2) and a in Eq. (5.2), which enter in Eq. (4.12). In Fig. 8, we show the fractional asymmetries versus mass for different choices of η (e.g. , see Eq. (3.2)). In the left plot of Fig. 8, it can be seen that the solid curves ( = 0.1) are larger than the dashed ( = 1) ones for the same value of a, due to the suppression by η. At the same time, curves corresponding to larger values of a have r 1 and can be treated as completely asymmetric, while the curves corresponding to tiny values of a are close to one, and are mostly symmetric. This is reflected in the relic abundance of the species, where we see that for larger values of a, the curves coincide and are a linear function of m (i.e., the solid red and the dashed blue curves), whereas for smaller values of a, the abundances scale as ∼ m 2 i , similar to what we observed for the symmetric components.

Heavy mediators
In the case of heavy mediators, as heavier components have larger cross sections, they are expected to be highly asymmetric in nature compared to the lighter ones. In Fig. 9, we show the asymmetry ratios and abundances for different choices of b (annihilations) and .
We fix the mediator mass scale at 100 GeV and vary the DM masses in the range 2 − 20 GeV. Comparing the asymmetry ratios with the light mediator case, c.f. Fig. 8, we see a reversed behaviour and now heavier particles in the set-up tend to be more asymmetric than the lighter ones. The suppression due to is evident in the left plot. Depending on the strength of annihilations, the abundance of heavier particles that are completely asymmetric is proportional to their mass whereas the abundance of the lighter ones roughly scales as ∼ m −2 , similar to what we had seen in the symmetric case. Notice also how, the smaller the value of , the larger the value of the DM mass where the abundance changes behaviour from mainly being symmetric-dominated to being asymmetric-dominated.

Effect of conversions
Here we discuss the effect of including conversions for both cases above. Since conversions deplete the symmetric components, the heavier particles which can convert to the lighter ones become more asymmetric, compared to the case of no conversions. However, the heavier states may be more symmetric (asymmetric) than the lighter ones in the case of light (heavy) mediators (c.f the left panels in Figs. 8 and 9). Therefore, in presence of strong conversions, the abundance of the heavier species scales proportional to their mass and the initial asymmetry. For weak conversions, the behavior is similar to the case of only annihilations. Depending on the strength of conversions versus annihilations, one can categorise the multi-component asymmetric DM set-up as shown in Table 1. The only case where the DM can be treated as partially asymmetric is when both annihilations and conversions are weak; all other combinations lead to a completely asymmetric behavior.

Annihilation/Conversion Weak Strong
Weak Partially asymmetric Asymmetric Strong Asymmetric Asymmetric Table 1. Behavior of multi-component DM (except for the lightest component), depending on interaction strength. 'Asymmetric' denotes that the DM is completely asymmetric.

Mass-dependence: power-law scaling
For the parameterisations of annihilation and conversion cross-sections used we observe that the multi-component models exhibit a power-law behavior. The scaling of the relic abundance with DM mass is summarised in Table 2. The power laws are a good fit for the heavier components in the set-up and do not apply to the lightest component. 8
Scaling of the relic abundance with the DM mass for the heavier components in a multi-component DM with symmetric or asymmetric components, with and without conversions, for different types of dependence of the cross sections on the DM mass. The left/right exponents indicate the extreme limits that are obtained for strong/weak conversions (and annihilations) in the symmetric (partially-asymmetric) case. In the cases with an * , left/right exponents indicate weak/strong interactions.
We show the case of treating all cross sections as mass-independent (σ = f (m i )), as well as the cases of light and heavy mediators. In the mass-independent symmetric case without conversions, all of the particles contribute with a similar abundance to the total one. Similarly, when the cross sections are inversely proportional to the DM mass, e.g. for the light mediator case, strong conversions may lead to several stable components with similar abundances in the symmetric case, e.g. the relic abundance becomes independent of the DM mass. On the other hand, in models with heavy mediators, we expect the lighter components of the dark sector to dominate the relic abundance in presence of significant conversions.
For the partially-asymmetric case, if the annihilations and conversions are strong, the heavier components are highly asymmetric and the abundance depends linearly on mass, whereas if both the processes are weak, we recover the dependence found in the symmetric case with weak conversions. Therefore, depending on the values of interaction cross-sections, it may be possible to have several partially-asymmetric particles with similar abundances in the heavy mediator case or in the mass-independent case with weak annihilations and conversions.

Abundance of the lightest state
The power laws provided above are a good fit to describe the behavior of heavier components that undergo conversions. However, this is not true for the lightest component that is being produced by all the heavier components. In this section, we comment on an interesting effect that the lightest component undergoes in presence of conversions, which to our best of knowledge has not been discussed elsewhere. It turns out to be useful to define R = σ conv, 12 v / σ ann, 1 v as the ratio of the conversion and the annihilation thermally-averaged cross-sections. In Fig. 10, we plot the abundance of the heaviest and the lightest component versus R, for N = 2, in the massindependent case. In general, one would expect the abundance of the lightest (heaviest) component to be directly (inversely) proportional to the strength of conversions (i.e. from annihilations of the heavy states into it). While the abundance of heaviest does decrease with conversions, for the lightest component, we find that the abundance does not scale directly with the conversion rate. We observe the presence of a peak in the abundance when R 1, after which the abundance drops but is always (slightly) larger than the case without conversions (R = 0). We observe that the peak appears at R 1, a phenomenon we term as 'annihilation-conversion resonance'. The presence of this peak is inherent to conversions and can be attributed to the interplay of freeze-out dynamics for all the components involved. A word of caution, however, is in order: for R 1, the computation of the relic abundance becomes more complicated, with the possible presence of significant reheating in the dark sector [39].
A suitable fit to the profile of the peak of the lighter component in Fig. 10 (right panel) is given by where α and γ are constants to be determined from the fit, and Ω l h 2 (R) [Ω l h 2 (0)] is the abundance of the lightest component in the case of conversions [no conversions]. For the heavier components (Fig. 10, left panel), the fit goes as α /(1 + R), where α is a constant. Such peaks in the lightest component abundance can be seen in all the scenarios discussed above, whether the DM is symmetric/partially asymmetric or heavier/lighter than the mediator involved, for N ≥ 2, as long as the mass splitting between the DM components fulfills ∆ ij ≥ 0.05. Note that the mass splitting ∆ ij of both DM species plays a central role in the above discussion as a consequence of Eq. (5.4), since for small ∆ ij no peak in the abundance is observed. In Fig. 11, we show the case of mass-dependent cross-sections, using the parameterization of Eqs. (5.2) and (5.4). We consider the light mediator case with different values of a (annihilations) and c (conversions).

Phenomenological Implications
In general, the phenomenological implications are model dependent, and mostly driven by the interaction strength with the SM (annihilations), and to a lesser extent by conversions. Assuming similar values for the different components, the following generic features are expected.
In multi-component DM scenarios, if there are no conversions, typically the component with the smallest annihilation cross section dominates the abundance. Given that indirect detection (ID) signals are proportional to n 2 σ ann,i v ∝ 1/(m 2 σ ann,i v ), for equal massindependent annihilation cross sections, we expect the lightest component to dominate. In presence of conversions, the abundance of the heavier components gets reduced and thus the ID signals should be suppressed compared to the case with no conversions. For example, in 2DM case, the ID rates for the heavier component χ 1 and the lighter component χ 2 are The scaling of indirect detection rates with DM mass for the different scenarios are shown in Table 3. We observe that in the light mediator case with only annihilations, or in presence of weak conversions, the ID rates of the different DM components are of similar order.

Cross-section Indirect Detection
Model σ Only Ann. Ann.+Conv. In the partially-asymmetric case with 1DM, there is an exponential reduction of the asymmetry ratio (and therefore on the indirect signals) with the annihilation cross section compared to the vanilla symmetric case, as pointed in Ref. [5]. In the case of partiallyasymmetric multi-component DM, we find an exponential dependence of the asymmetry ratios on the conversion cross section for the heavier species, making them more asymmetric, whereas the lighter ones may remain mostly symmetric. Therefore, with strong conversions, the heavier species become even more asymmetric, e.g. r i 1, and their indirect detection signals are highly suppressed compared to a WIMP candidate of the same mass [5], The ID rates for lighter components are also suppressed but not as much as the heavier ones R ID asym R ID sym = σ ann, l σ sym × r l,∞ × 2 1 + r l,∞ 2 1 .
On the other hand, if the conversions are weaker and the annihilation rates are similar to the thermal one, the symmetric part of the components may be similar, as already discussed above. One may wonder whether strong conversions may lead to a locally-enhanced population of the lighter states compared to the heavier ones in high DM-density environments, such as the galactic centre or dwarf galaxies, enhancing the previous rates. Schematically, assuming that the local energy density ρ local in the considered body scales as the global one Ω global , as is the case for cold DM [42,43] If the lighter states become more populated, their annihilation channels, say into SM states α, will dominate. This would be interesting phenomenologically, given the stronger ID limits on light DM. For example, in the case of 2DM with masses m 1,2 (m 1 > m 2 ) and number densities n 1,2 , conversions will lead to an increase in the number density of the lighter component. One can estimate ∆n 2 2 n 2 1,i σ conv, 12 v t , (6.5) where t is the relevant time scale from the formation of the local body, and n 1,i is the initial local number density of the heavier component at the time where the body is formed. The change in the lighter species may be estimated as where P ≡ ρ 1 local /ρ local = Ω 1 global /Ω global is the initial fraction of the energy density in the heavy state. Therefore, the current local number densities of the states become n 1,f n 1,i − ∆n 2 , n 2,f n 2,i + ∆n 2 . (6.7) Though in this case the change is negligible, it could be significant for local bodies with much higher dark matter densities, for lighter DM, and for larger conversions. A detailed analysis of the evolution of those systems would be needed in that case. Multi-component DM could also give signatures in direct detection, specially kinks in the energy spectra of averaged rates [44,45], and partial cancellations and a non-sinusoidal behavior in time-dependent rates [46] (being able to reproduce DAMA Phase 2 results). The implications on direct detection rate for reproducing the relic abundance in thermal-freezeout and asymmetric scenarios have been discussed in Ref. [47]. In general, the presence of large conversions allows to reproduce the relic abundance with smaller DM-SM couplings, hence relaxing the strong bounds coming from DD null-results.
Astrophysical signals from capture in the sun and other celestial bodies (e.g. neutron stars) can constrain the set-up. For example, the asymmetric component may yield constraints from creating black-holes, etc. The symmetric part may yield the typical signals, for instance annihilation of DM components into neutrinos, and/or anomalous heating of the celestial body. The conversions χ 1 χ 1 → χ 2 χ 2 with m 1 m 2 (inelastic scattering) can yield a larger population of the lowest mass state in the Sun (assuming thermalisation is reached) [48]. The opposite reaction χ 2 χ 2 → χ 1 χ 1 is kinematically forbidden in typical bodies for non-negligible mass splittings.
Similarly, the conversions χ 1 χ 1 → χ 2 χ 2 with m 1 m 2 (inelastic scattering) can explain small-scale structure problems in dwarfs if large enough (also bounds from clusters) [49] and again yield a larger population of the lowest mass state. The opposite reaction χ 2 χ 2 → χ 1 χ 1 is typically kinematically forbidden also in clusters for non-negligible mass splittings.
At colliders, even if strong conversions make the DM to be mainly composed of the lighter state, the heavier states, if kinematically accessible, may also be produced. For LHC studies of multi-component scenarios see for example Ref. [50,51].

Conclusions
Guided by the visible sector, one may argue that the most natural dark sector is multicomponent and partially-asymmetric. In this work, we have studied the relic abundance of several stable particles in the dark sector with symmetric and asymmetric components, with special emphasis on the role that conversions have in the final relic abundance. We have found that the total abundance gets reduced by conversions in all cases, generically helping to hide the dark sector from experiment.
In many scenarios, we have found that several DM components may have similar abundances, contributing on equal footing to the relic abundance. These are the cases of: i) symmetric or partially-asymmetric DM with mass-independent cross sections and small conversions, or light-mediators with strong conversions; ii) heavy-mediator with weak conversions. We expect the dependence on the DM mass shown in Table 2 to be valid for a large class of models, and to be rather universal as long as the masses and cross sections are at the weak scale.
In the symmetric case, for mass-independent cross sections, the abundance of heavier species is suppressed due to significant conversions and they are mostly asymmetric, while the lightest ones dominate the energy density. Furthermore, for large conversions and low masses, the local density of lighter states may be enhanced with respect to the global one in very high-density environments. These effects specially affect ID signals, with lighter states being more populated and yielding the largest SM final state products, while heavier ones remain mostly asymmetric and less populated.
We conclude that, in multi-component scenarios, it is natural to have symmetric and asymmetric components in presence of conversions. These can help to explain the absence of any positive signal up to now, and furthermore provide distinctive experimental signatures in current and future direct and indirect detection experiments, as well as at colliders, where cosmologically subdominant components may be produced.

A Model Example
Here we provide an example model to realize multi-component DM and demonstrate the behavior of the cross section with DM mass, in line with the discussion in Section 5. One of the simplest models where this can be seen is the scalar singlet DM model, where the DM couples to the standard model via Higgs portal [52,53]. Most of the parameter space of the real scalar singlet model is ruled out, leading to DM masses either in the funnel region of Higgs resonance or m DM > 1 TeV [54,55]. In the following, we use it as an illustration and extend it to N components. The parameter space for the two component model has been studied in detail in Ref. [11].
Consider the standard model extended by N scalar singlets φ i , all of them stabilised by Z 2 × Z 2 . . . × Z N 2 symmetry and coupling to the SM via a Higgs portal coupling λ i . The part of the Lagrangian concerning DM interactions is given by where λ s i is the quartic self coupling, m 0 i is the DM bare mass, λ ij is the DM-DM coupling and H is the SM Higgs doublet, which acquires a vacuum expectation value after electroweak symmetry breaking and is parameterised by H = (0, (υ + h)/ √ 2) T , where υ = 246 GeV, while φ i = 0. Hence, the mass of the DM candidates are m 2 i = m 2 0 i + λ i υ 2 2 and we choose m i > m j for j > i. The Higgs portal couplings (λ i ) opens annihilation channels for DM states, φ i φ i → f f, W W, ZZ, hh. Furthermore, they also lead to conversions φ i φ i → φ j φ j along with the direct interactions λ ij , which dominate for larger DM masses.
In the large DM mass regime (light mediator case), m i m h , the annihilation cross section (s-channel) for the model can be roughly approximated as According to our adopted convention, this implies a i ≡ λ 2 i /16π. For conversions, the interaction rate is given by Taking the threshold value s = 4m 2 i , in the heavy mass regime, the direct interaction term λ ij dominates over λ i,j , and we can approximate Therefore, c ij ≡ λ 2 ij /32π.
In the low DM mass regime (heavy mediator case), the annihilation and conversion cross sections can be approximated by The approximation for the conversion cross section above holds in the case where the DM mass is given by m 2 i ∼ λ i υ 2 . Thus, b i ≡ λ 2 i /8π and d ij ≡ λ 2 j /8π in this case.