Addressing Astrophysical and Cosmological Problems With Secretly Asymmetric Dark Matter

We present a simple model of dark matter that can address astrophysical and cosmological puzzles across a wide range of scales. The model is an application of the Secretly Asymmetric Dark Matter mechanism, where several flavors of dark matter are fully asymmetric despite an exact dark matter number symmetry $U(1)_{\chi}$. The dark matter relic abundance arises from these asymmetries, generated in the early universe through right-handed neutrino decays. The $U(1)_{\chi}$ is gauged by a massless dark photon, and asymmetries with opposite signs in the different DM flavors result in the formation of bound states. Dark acoustic oscillations in the early universe lead to a suppression in the matter power spectrum for addressing the $\sigma_8$ problem. The dark photon as a relativistic degree of freedom contributes to $\Delta N_{eff}$, easing the tension between the CMB and late-time $H_{0}$ measurements. Finally, scattering between the bound states after structure formation leads to a flattening of the dark matter distribution at the cores of haloes.


Introduction
Secretly Asymmetric Dark Matter (SADM) has recently been proposed [1] as a mechanism to generate the dark matter (DM) relic abundance through an asymmetry in the dark sector despite an exact (global or gauged) DM number symmetry U (1) χ . In the implementation of this idea in ref. [1], the asymmetry is first generated in the visible sector through highscale leptogenesis [2] 1 and then transferred to a dark sector containing three DM flavors χ i , a flavorless mediator φ, and a dark photon γ D via a coupling to the right-handed leptons of the Standard Model (SM). The mediator in that case is charged both under SM hypercharge and under U (1) χ , and therefore leads to mixing between the photon and the dark photon at the one-loop level. This is phenomenologically only acceptable if the dark photon has a nonzero mass [5], and therefore the experimental signatures of the model manifest themselves at short distances, but in the context of long distance physics it falls under the collisionless cold DM category.
In this work we introduce a different implementation of the SADM mechanism, where a massless dark photon is phenomenologically allowed, and leads to a rich set of astrophysical signatures, despite the rather minimal high-energy setup. We argue that this setup is capable of addressing open questions in astrophysics both at the large (Hubble-scale) and small (galaxy-scale) scales simultaneously. The model we propose serves as a short-distance implementation of two mechanisms that have recently been proposed to address the same open questions, namely Dark Acoustic Oscillations [6][7][8][9][10] and Hidden Hydrogen DM [12,13]. By focusing on a particular region of parameter space of the model, we show that the H 0 discrepancy [14], the σ 8 discrepancy [15][16][17][18][19] and the mass deficit problem [15,20,21] in dwarf galaxies and galaxy clusters can all be addressed. There have been many efforts in the literature to address these discrepancies using non-standard DM self interactions [22][23][24][25][26][27][28][29][30][31]. The SADM model presented here extends these efforts by providing a minimal shortdistance setup that nevertheless has rich enough dynamics to address multiple issues across a wide range of scales.
The core feature of the SADM mechanism is that due to the unbroken DM number symmetry, the total dark charge of the universe is zero at all times, but there can be several DM flavors that have individual asymmetries. Due to the crucial role of there being multiple DM flavors, SADM can be considered as a special case of the Flavored Dark Matter scenario [32][33][34][35][36][37][38]. Unlike the model studied in ref. [1] however, where the heavier DM flavors decay to the lightest one before Big Bang Nucleosynthesis (BBN) and result in a symmetric distribution of the lightest DM flavor at late times, in the model we study in this work all DM flavors are extremely long-lived and are still present today. This means that flavors with opposite signs of asymmetry (and with opposite dark charges) must coexist, along with a massless dark photon that mediates long-range interactions, and therefore multiple bound states can form, each one behaving as atomic DM [22,39,40]. These bound states, along with unbound DM particles and the massless dark photon then give rise to rich dynamics across a range of distance scales. The dark photon as an additional relativistic degree of freedom helps ease the tension between CMB-based [41] and low redshift [14,42,43] measurements of H 0 , dark acoustic oscillations in the early universe that arise from the scattering of free and bound states leads to a suppression of the matter power spectrum at small scales, and the scattering between bound states after structure formation leads to a flattening of the DM density distribution inside haloes.
The layout of this paper is as follows: In section 2 we introduce the short-distance model, and we consider its phenomenology in the context of particle physics observables and constraints. In section 3 we then study the cosmological features of the model, and the astrophysical signatures that it gives rise to at the large and small scales. We conclude in section 4 and we consider directions for future research.

The Model, and particle physics considerations
In this section we will introduce our SADM model and consider its high energy aspects. We assume heavy right handed neutrinos N i exist and that they have Majorana masses as well as a Yukawa coupling to the SM leptons where H is the SM Higgs doublet. This well-studied extension of the SM generates both light neutrino masses through the see-saw mechanism and it also allows the creation of a net B − L number as the origin of the matter-antimatter asymmetry in the SM sector through CP violating phases in the coupling matrix Y N ij . These features are explored in depth in reviews of leptogenesis [3,4].
To this setup we add three flavors of Dirac fermions (χ, χ c ) i and a scalar φ with the couplings (in the χ-mass basis) L ⊃ M 2 φ |φ| 2 + m χ,i χ c i χ i + λ ij N i χ j φ + h.c.. (2.2) The λ ij coupling matrix contains physical phases and is a source of CP-violation in the dark sector. Purely for reasons of simplicity we choose the N to only couple to the left-handed χ but not the right-handed χ c . A more general setup with both types of couplings and with independent coupling matrices λ ij andλ ij is possible, and this would not significantly change the conclusions of our paper, as long as the couplings are still CP-violating. We will soon introduce a convention for assigning the flavor labels i. Note that an exact U (1) χ symmetry exists, under which all χ have charge +1, while all χ c and φ have charge -1. We will take this symmetry to be gauged by the dark photon, with a fine structure constant α d . The mediator φ, being a scalar, will be taken to be heavy, though we will assume that it is lighter than the right-handed neutrinos. Once the right-handed neutrinos become non-relativistic and the interaction of equation 2.2 drops out of equilibrium, the SM and dark sectors decouple from one another. As the heavy neutrinos N decay, they will then generate asymmetries both in the SM leptons and in the dark sector (in φ and the χ i ); however as mentioned before U (1) χ is never broken and therefore the total dark charge of the universe is zero at all times. With the interaction of equation 2.2 out of equilibrium, the asymmetries generated in the three χ flavors cannot wash each other out. As the temperature drops further, the φ particles become nonrelativistic and the symmetric part of the φ particle distribution annihilates to dark photons. The asymmetric φ particles then decay as φ → χ ± H ∓ , transferring their dark charge to the χ flavors, such that after the φ decays the asymmetry in the dark sector will reside entirely in the three DM flavors, in such a way that the total dark charge of the universe remains zero. Note that the heaviness of both φ and N means that the χ are very challenging to access in collider, direct detection, or indirect detection experiments. Our model does however have astrophysical signatures, which are the focus of this paper and which we will study in the next section.
Assuming generic phases in the entries of the Y N ij and the λ ij matrices, the ratio of the typical magnitudes of the elements of Y N ij and the typical magnitudes of the elements of λ ij will determine the branching ratio of the decaying N into the visible vs the dark sector, and thereby the ratio of the overall lepton number (more precisely, B − L number) and the flavor-by-flavor χ asymmetries that are generated through heavy neutrino decays.
Therefore, for similar magnitudes in the entries of the Y N ij and the λ ij matrices, the fact that Ω DM ∼ 5Ω B suggests the GeV scale for the mass of the heaviest χ flavor, but if the entries of one coupling matrix are larger than those of the other, then this mass can also vary in either direction. As we will show in section 3, the region of parameter space that is of interest for us has the mass of the heaviest DM flavor ∼ O(10) GeV, the mass of the lightest DM flavor ∼ O(1) MeV, and a value for α d of order several percent. In this parameter region, the symmetric parts of the χ distributions annihilate efficiently, leaving behind only the asymmetric part generated during the N decays.
The asymmetry generated in the different χ flavors will differ from one another. In fact, due to the total dark charge of the universe being zero, there will always be one χ flavor with one sign of the asymmetry, and two flavors with the opposite sign, after the φ particles have all decayed. We will adopt the following naming convention for the three χ flavors: The χ flavor which is unique in its sign of the asymmetry will be labeled χ 1 , whereas the two χ flavors with the other sign of the asymmetry will be labeled χ 2 and χ 3 , such that m χ 2 > m χ 3 . As we will see in the next section, the most interesting phenomenology is obtained when χ 3 and not χ 1 is the lightest DM flavor. Thus we will take m χ 1 > m χ 3 as well, while either one of χ 1 or χ 2 could be the heaviest flavor. The possible bound states are then formed between χ 1 and χ 2 , labeled H 12 , and between χ 1 and χ 3 , labeled H 13 . Since the total dark charge of the universe is zero at all times, the asymmetries of χ 2 and χ 3 add up in magnitude to that of χ 1 .
The overall size of the entries of the λ ij coupling matrix determines the branching fraction of N to the dark sector, and therefore the overall size of the asymmetries of all χ flavors, while the specific flavor and phase structure in the coupling matrix determines the branching fractions among the DM flavors. The three comoving asymmetries ∆Y i of the χ flavors (after φ decays) satisfy 3 i=1 ∆Y i = 0. In terms of the physical number densities of the three χ particles (or antiparticles, depending on the sign of the asymmetry), this can also be written as n 1 = n 2 + n 3 . Thus the n i are not independent, however the ratio n 3 /n 1 can in principle have any value in the interval [0, 1].
We can estimate the probability distribution for this ratio using Monte Carlo methods with a prior where each entry of the 3 × 3 matrix λ ij is randomly chosen over the unit complex circle. To calculate the asymmetry generated in each flavor from a given matrix λ ij , we follow standard techniques from leptogenesis [3]. We assume the heavy neutrino masses are hierarchical with m N 1 m N 2 , m N 3 (we use m N 1 = 1.0 × 10 16 GeV as a benchmark value). N 1 decays are fast compared to the Hubble scale, and therefore the asymmetry production takes place in the strong washout regime. Ratios of the ∆Y i generated in the N 1 decays can be obtained form the ratios of the CP asymmetry factors for each flavor i where the matrix [m] ij is given by Figure 2. Decay mode for a heavier DM flavor to decay to a lighter one.
For the branching ratios of the subsequent φ decays, only the tree level contributions with an off-shell N 1 are used, since any additional CP violation introduced at this stage is subdominant. Note also that an overall rescaling of all λ ij entries does not affect the probability distribution of n 3 /n 1 , and therefore we are not committing ourselves to a particular size of the λ ij couplings by generating random numbers over the complex unit circle. For a given choice of m N 1 and m χ i , such a rescaling can be chosen such that the correct overall ρ DM and ρ B are obtained. The result for the n 3 /n 1 probability distribution is shown in figure 1. As we will see in the next section, the region of greatest phenomenological interest corresponds to n 3 /n 1 < ∼ 0.1, which does not require significant fine-tuning. Even without a compressed spectrum among the DM flavors, all χ are extremely longlived, and are still undecayed today. The leading decay mode arises from the one-loop diagram shown in figure 2. After integrating out N and φ, this diagram gives rise to a flavor off-diagonal dipole term for χ at low energies, the coefficient of which can be estimated as (we use λ to stand in for any O(1) entries in the λ ij coupling matrix): The lifetime is many orders of magnitude larger than the age of the universe for the parameter region that will be used in this paper. The leading tree-level decay diagrams are even further suppressed than this decay mode, and are therefore irrelevant. As we will see in the next section, the most interesting phenomenological features of the model will be obtained when the mass of the heaviest flavor (χ 1 or χ 2 ) is O(10) GeV, the mass of the next heaviest flavor is O(0.1 − 1) GeV, and the mass of the lightest flavor (χ 3 ) is O(1) MeV. This means that H 12 will be strongly bound and difficult to ionize, serving as a neutral DM component, while H 13 may be easy to ionize, such that a nonnegligible fraction of χ 3 (and the corresponding amount of χ 1 ) can be unbound at certain times during the early universe, delaying the decoupling of the dominant DM component from the dark radiation, and leading to interesting dynamical effects. For a massless dark photon, it needs to be checked that the mixing with the SM photon is acceptably small. More precisely, in the basis where γ D only interacts with χ but not the SM fermions, the millicharge acquired by χ under electromagnetism needs to be small. The γ-γ D mixing is induced at the three-loop level: see figure 3. Matching on to an effective operator F μν D F μν , one can estimate Let us compare this to the constraints for millicharged DM (figure 1 in ref. [44], with m ∼ O(10) GeV, the mass of the heaviest DM flavor), keeping in mind that some of the constraints will be relaxed, as these were derived with the assumption that all DM particles are millicharged while we will focus in this work exclusively on the possibility that all but a small fraction of the χ particles are in bound states, and therefore have no net millicharge. In particular, the leading constraints from galactic and cluster magnetic fields probe the DM at small momenta, or at distance scales larger than the size of the bound states, thus these constraints cannot resolve the millicharges of the constituents. While the momentum scale that is probed at LUX is still too small to resolve the constituents of H 12 , it may be sufficiently large to resolve the constituents of H 13 . However, in the parameter space of interest to us, H 13 is the subdominant DM component, which relaxes the constraints below the generic values of in our model. Furthermore, for a significant fraction of the parameter region of interest to us (especially when m χ 2 > m χ 1 ), the mass of H 13 is below the sensitivity of LUX. It is however interesting to note that future direct detection experiments with a low mass threshold may be able to test the scenario presented here. Finally, for the parameter region of interest, the value of in our model is also below the bounds from supernovae [5].
As a final consideration about the particle physics nature of our model, we want to address the potential concern that flavor oscillations between the DM flavors may cause a washout of the asymmetries, as the particles of χ 1 oscillate into the antiparticles of χ 2 or χ 3 , or the other way around. However, the large mass gap between the χ flavors in the parameter region of interest in this paper, and the extremely small intrinsic widths of these states makes oscillations so small that they are negligible for all intents and purposes.

Cosmology and astrophysical signatures
Having discussed the short-distance physics aspects of our model, we now turn our attention to the role it plays in the dynamics of the early universe, as well as in late-time astrophysical processes. Although the dark sector is hidden from direct experimental probes via nongravitational interactions, the dynamics within the SADM sector can significantly change the gravitational potential in the early universe and it leaves visible signals on both the large and small scale structures. We will show that there is a region of parameter space where observational anomalies on scales of dwarf galaxies (∼ kpc size) [15,20], galaxy clusters (∼ Mpc size) [21], the σ 8 problem (∼ 10 Mpc size) [15][16][17][18][19], and the tension in the measurements of H 0 from the CMB [41] and from low redshift measurements [45] may all be addressed.
The input parameters of the model that are relevant at large distances can be listed as the comoving asymmetries of the χ flavors (of which only two are independent, and which can therefore be parameterized in terms of n 1 and n 3 /n 1 once the symmetric part of the χ distributions annihilates away), the masses of the three χ flavors, and α d . In order to calculate how the temperature of the dark sector is related to that of the SM sector, consider energies above m N where the interaction of equation 2.2 keeps the two sectors in equilibrium. As the temperature drops below m N , the visible and dark sectors decouple. At this time, the dark sector contains the three Dirac fermions χ 1,2,3 , the complex scalar φ, and the dark photon γ D . φ is only somewhat lighter than the N , and for the parameter region of interest for us, the mass of the lightest flavor is O(1) MeV. Therefore, at the time of structure formation, only γ D is still relativistic.
If there are no other degrees of freedom in the visible sector other than those of the SM up to energies of m N , then by using the ratio of the number of relativistic degrees of freedom in the two sectors at the decoupling scale and at the scale of structure formation, one obtains ∆N ef f = 0.75. However, since the two sectors decouple at a very high temperature scale, one needs to take into account the possibility that the visible sector contains additional degrees of freedom (e.g. connected to the solution of the electroweak naturalness problem, grand unification, etc.) that release additional entropy in the visible sector as the universe cools down, and thereby further suppress ∆N ef f . The value 0.75 should therefore be considered as an upper limit, with the actual value depending on other possible extensions of the SM. Since a value of 0.75 is outside the 2σ contour for reconciling the H 0 discrepancy (see e.g., [41]), we will adopt a benchmark value of ∆N ef f = 0.60 for the remainder of the paper 2 , which can be obtained for example if there are six Dirac fermions beyond the SM in the visible sector. Adding even more degrees of freedom to the visible sector would further reduce ∆N ef f . It should also be noted that at the time of BBN, ∆N ef f is even smaller (∆N ef f = 0.42) since χ 3 has not yet become nonrelativistic at that time for m χ 3 < ∼ MeV, which is the preferred value for addressing the σ 8 problem as we will see later. Thus our model is compatible with BBN constraints [47].
There are two recombination processes during the early universe, of χ 1 with χ 2 into H 12 , and of the remaining χ 1 with χ 3 into H 13 . H 12 forms earlier due to its larger binding energy, and the remaining χ 1 and χ 3 particles scatter with each other and remain in thermal equilibrium at this time. Similar to the proton-hydrogen scattering in the SM sector, the scattering cross section between the non-relativistic χ 1 and the H 12 bound state is larger than the geometrical size of H 12 [48], sufficient for keeping them in thermal equilibrium with each other. While χ 1 -χ 3 scattering is also efficient, the entire dark sector is then in thermal equilibrium with the dark photon, and the resulting DM oscillations delay the formation of large scale structure. The oscillation stops as H 13 recombines, which leaves too few free (dark-)charged particles to sustain the oscillations. These dark acoustic oscillations generate a small but visible damping of the matter power spectrum and may provide a solution to the σ 8 problem. A more complete parameter fitting procedure including also the CMB and BAO data is necessary to confirm this claim in full detail, however in this work we will take a simpler approach and we will calculate the size of power spectrum suppression to argue that the claim is plausible, leaving a more detailed analysis to future work.
We will also show that at later times, during the formation of DM halos, the scattering between H 12 bound states leads to thermal equilibrium and provides a solution to the core/cusp problem of dwarf galaxies. The non-trivial velocity dependence due to the inelastic scattering H 12 H 12 → H 12 H 12 γ D through the dark hyperfine transition also gives the right cross section for cores to form in relaxed cluster halos. As we will see, achieving this will favor the parameter region where the heaviest DM flavor mass is ∼ O(10) GeV, where the value of α d is a few percent, and n 3 /n 1 ∼ 0.1. Furthermore, we will see that achieving the correct amount of damping of large scale structure favors the MeV range for m χ 3 . Thus there is a region of the parameter space of our model where all parameters are technically natural and where all three structure formation problems (dwarf, cluster and σ 8 ) and the H 0 tension could be addressed. Below, we explore each of these aspects of our model in detail.

Dark recombination(s)
The formation of the H 12 and H 13 bound states plays an important role in structure formation. In particular dark acoustic oscillations end when there are no longer sufficiently many unbound χ particles left to sustain them. Similar to the recombination of hydrogen in the SM, for both H 12 and H 13 , recombination proceeds through the formation of the excited states (n ≥ 2), with a subsequent decay into the ground state either through a two photon emission or through a Lyman-α transition, where the dark photon becomes redshifted before ionizing other bound states. A brief note on notation: we have introduced n i in the previous section to denote the physical densities in each flavor. Since in this section we will also need to keep track of just the density of unbound particles in each flavor, we will introduce the notation n f i , where f stands for "free". With this definition n 1 = n f 1 + n H 12 + n H 13 , n 2 = n f 2 + n H 12 and n 3 = n f 3 + n H 13 . We also define the ionization fraction in each flavor as X i ≡ n f i /n i .
< l a t e x i t s h a 1 _ b a s e 6 4 = " 3 5

< l a t e x i t s h a 1 _ b a s e 6 4 = " 3 5 9 A p Q n b t l x x s Z U d w g L j W j Z l E N w = " > A A A B 8 3 i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L J a C p 5 p V Q S 9 C 0 Y v H C v Y D m h A 2 2 0 2 7 d L M J u x u h h P 4 N L x 4 U 8 e q f 8 e a / c d v m o N U H A 4 / 3 Z p i Z F 6 a C a + O 6 X 0 5 p Z X V t f a O 8 W d n a 3 t n d q + 4 f d H S S K c r a N B G J 6 o V E M 8 E l a x t u B O u l i p E 4 F K w b j m 9 n f v e R K c 0 T + W A m K f N j M p Q 8 4 p Q Y K 3 k y O D + V A b 7 G r l c P q j W 3 4 c 6 B / h J c k B o U a A X V T 2 + Q 0 C x m 0 l B B t O 5 j N z V + T p T h V L B p x c s 0 S w k d k y H r W y p J z L S f z 2 + e o r p V B i h K l C 1 p 0 F z 9 O Z G T W O t J H N r O m J i R X v Z m 4 n 9 e P z P R l Z 9 z m W a G S b p Y F G U C m Q T N A k A D r h g 1 Y m I J o Y r b W x E d E U W o s T F V b A h 4 + e W / p H P W w G 4 D 3 1 / U m j d F H G U 4 g m M 4 A Q y X 0 I Q 7 a E E b K K T w B C / w 6 m T O s / P m v C 9 a S 0 4 x c w i / 4 H x 8 A 9 p E k D 4 = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " 3 5 9 A p Q n b t l x x s Z U d w g L j W j Z l E N w = " > A A A B 8 3 i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L J a C p 5 p V Q S 9 C 0 Y v H C v Y D m h A 2 2 0 2 7 d L M J u x u h h P 4 N L x 4 U 8 e q f 8 e a / c d v m o N U H A 4 / 3 Z p i Z F 6 a C a + O 6 X 0 5 p Z X V t f a O 8 W d n a 3 t n d q + 4 f d H S S K c r a N B G J 6 o V E M 8 E l a x t u B O u l i p E 4 F K w b j m 9 n f v e R K c 0 T + W A m K f N j M p Q 8 4 p Q Y K 3 k y O D + V A b 7 G r l c P q j W 3 4 c 6 B / h J c k B o U a A X V T 2 + Q 0 C x m 0 l B B t O 5 j N z V + T p T h V L B p x c s 0 S w k d k y H r W y p J z L S f z 2 + e o r p V B i h K l C 1 p 0 F z 9 O Z G T W O t J H N r O m J i R X v Z m 4 n 9 e P z P R l Z 9 z m W a G S b p Y F G U C m Q T N A k A D r h g 1 Y m I J o Y r b W x E d E U W o s T F V b A h 4 + e W / p H P W w G 4 D 3 1 / U m j d F H G U 4 g m M 4 A Q y X 0 I Q 7 a E E b K K T w B C / w 6 m T O s / P m v C 9 a S 0 4 x c w i / 4 H x 8 A 9 p E k D 4 = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " 3 5 9 A p Q n b t l x x s Z U d w g L j W j Z l E N w = " > A A A B 8 3 i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L J a C p 5 p V Q S 9 C 0 Y v H C v Y D m h A 2 2 0 2 7 d L M J u x u h h P 4 N L x 4 U 8 e q f 8 e a / c d v m o N U H A 4 / 3 Z p i Z F 6 a C a + O 6 X 0 5 p Z X V t f a O 8 W d n a 3 t n d q + 4 f d H S S K c r a N B G J 6 o V E M 8 E l a x t u B O u l i p E 4 F K w b j m 9 n f v e R K c 0 T + W A m K f N j M p Q 8 4 p Q Y K 3 k y O D + V A b 7 G r l c P q j W 3 4 c 6 B / h J c k B o U a A X V T 2 + Q 0 C x m 0 l B B t O 5 j N z V + T p T h V L B p x c s 0 S w k d k y H r W y p J z L S f z 2 + e o r p V B i h K l C 1 p 0 F z 9 O Z G T W O t J H N r O m J i R X v Z m 4 n 9 e P z P R l Z 9 z m W a G S b p Y F G U C m Q T N A k A D r h g 1 Y m I J o Y r b W x E d E U W o s T F V b A h 4 + e W / p H P W w G 4 D 3 1 / U m j d F H G U 4 g m M 4 A Q y X 0 I Q 7 a E E b K K T w B C / w 6 m T O s / P m v C 9 a S 0 4 x c w i / 4 H x 8 A 9 p E k D 4 = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " 3 5 9 A p Q n b t l x x s Z U d w g L j W j Z l E N w = " > A A A B 8 3 i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L J a C p 5 p V Q S 9 C 0 Y v H C v Y D m h A 2 2 0 2 7 d L M J u x u h h P 4 N L x 4 U 8 e q f 8 e a / c d v m o N U H A 4 / 3 Z p i Z F 6 a C a + O 6 X 0 5 p Z X V t f a O 8 W d n a 3 t n d q + 4 f d H S S K c r a N B G J 6 o V E M 8 E l a x t u B O u l i p E 4 F K w b j m 9 n f v e R K c 0 T + W A m K f N j M p Q 8 4 p Q Y K 3 k y O D + V A b 7 G r l c P q j W 3 4 c 6 B / h J c k B o U a A X V T 2 + Q 0 C x m 0 l B B t O 5 j N z V + T p T h V L B p x c s 0 S w k d k y H r W y p J z L S f z 2 + e o r p V B i h K l C 1 p 0 F z 9 O Z G T W O t J H N r O m J i R X v Z m 4 n 9 e P z P R l Z 9 z m W a G S b p Y F G U C m Q T N A k A D r h g 1 Y m I J o Y r b W x E d E U W o s T F V b A h 4 + e W / p H P W w G 4 D 3 1 / U m j d F H G U 4 g m M 4 A Q y X 0 I Q 7 a E E b K K T w B C / w 6 m T O s / P m v C 9 a S 0 4 x c w i / 4 H x 8 A 9 p E k D 4 = < / l a t e x i t >
N ef f = 0.6  D c s u x r b g 8 I p 0 k q K 9 O I i F Y + T G g 6 R h l s y q m y + Y p U q l U q u f Q f 2 e Q 5 N T y z L N K i y n S g G k a L n 5 d 3 s Q 4 I g T X 2 G G p O y X z V A 5 M R K K Y k a S n B 1 J E i I 8 Q S P S 1 9 R H n E g n n q + V w C O t D O A w E P r 6 C s 7 V 7 x 0 x 4 l J O u a c r O V J j + d u b i X 9 5 / U g N 6 0 5 M / T B S x M e L Q c O I Q R X A W U Z w Q A X B i k 0 1 Q V h Q / V e I x 0 g g r H S S O R 3 C 1 6 b w f 9 K x S m X N r 6 u F 5 n k a R x Y c g E N w D M q g B p r g E r R A G 2 B w B D c s u x r b g 8 I p 0 k q K 9 O I i F Y + T G g 6 R h l s y q m y + Y p U q l U q u f Q f 2 e Q 5 N T y z L N K i y n S g G k a L n 5 d 3 s Q 4 I g T X 2 G G p O y X z V A 5 M R K K Y k a S n B 1 J E i I 8 Q S P S 1 9 R H n E g n n q + V w C O t D O A w E P r 6 C s 7 V 7 x 0 x 4 l J O u a c r O V J j + d u b i X 9 5 / U g N 6 0 5 M / T B S x M e L Q c O I Q R X A W U Z w Q A X B i k 0 1 Q V h Q / V e I x 0 g g r H S S O R 3 C 1 6 b w f 9 K x S m X N r 6 u F 5 n k a R x Y c g E N w D M q g B p r g E r R A G 2 B w B D c s u x r b g 8 I p 0 k q K 9 O I i F Y + T G g 6 R h l s y q m y + Y p U q l U q u f Q f 2 e Q 5 N T y z L N K i y n S g G k a L n 5 d 3 s Q 4 I g T X 2 G G p O y X z V A 5 M R K K Y k a S n B 1 J E i I 8 Q S P S 1 9 R H n E g n n q + V w C O t D O A w E P r 6 C s 7 V 7 x 0 x 4 l J O u a c r O V J j + d u b i X 9 5 / U g N 6 0 5 M / T B S x M e L Q c O I Q R X A W U Z w Q A X B i k 0 1 Q V h Q / V e I x 0 g g r H S S O R 3 C 1 6 b w f 9 K x S m X N r 6 u F 5 n k a R x Y c g E N w D M q g B p r g E r R A G 2 B w B  Note that the dark charge neutrality of the universe enforces both n 1 = n 2 + n 3 and also n f 1 = n f 2 + n f 3 . The Boltzmann equation for X 2,3 , can be written as [13]

< l a t e x i t s h a 1 _ b a s e 6 4 = " B S Y A T f d a l l E m a 7 t 9 Z X o 2 w c 8 2 q v I = " > A A A C A H i c b Z D L S s N A F I Z P 6 q 3 W W 9 S F C z e D R X A h J R F B N 0 L R h S 4 r 2 A s 0 I U y m k 3 b o T B J m J k I J 2 f g q b l w o 4 t b H c O f b O L 0 s t P W H g Y / / n M O Z 8 4 c p Z 0 o 7 z r d V W l p e W V 0 r r 1 c 2 N r e 2 d + z d v Z Z K M k l o k y Q 8 k Z 0 Q K 8 p Z T J u a a U 4 7 q a R Y h J y 2 w + H N u N 5 + p F K x J H 7 Q o 5 T 6 A v d j F j G C t b E C + 0 A E u U c G L H C L K 9 c 7 z T 0 p 0 C 1 t F Y F d d W r O R G g R 3 B l U Y a Z G Y H 9 5 v Y R k g s a a c K x U 1 3 V S 7 e d Y a k Y 4 L S p e p m i K y R D 3 a d d g j A V V f j 4 5 o E D H x u m h K J H m x R p N 3 N 8 T O R Z K j U R o O g X W A z V f G 5 v / 1 b q Z j i 7 9 n M V p p m l M p o u i j C O d o H E a q M c k J Z q P D G A i m f k r I g M s M d E m s 4 o J w Z 0 / e R F a Z z X X 8 P 1 5 t X 4 9 i 6 M M h 3 A E J + D C B d T h D h r Q B A I F P M M r v F l P 1 o v 1 b n 1 M W 0 v W b G Y f / s j 6 / A G a 8 p X A < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " B S Y A T f d a l l E m a 7 t 9 Z X o 2 w c 8 2 q v I = " > A A A C A H i c b Z D L S s N A F I Z P 6 q 3 W W 9 S F C z e D R X A h J R F B N 0 L R h S 4 r 2 A s 0 I U y m k 3 b o T B J m J k I J 2 f g q b l w o 4 t b H c O f b O L 0 s t P W H g Y / / n M O Z 8 4 c p Z 0 o 7 z r d V W l p e W V 0 r r 1 c 2 N r e 2 d + z d v Z Z K M k l o k y Q 8 k Z 0 Q K 8 p Z T J u a a U 4 7 q a R Y h J y 2 w + H N u N 5 + p F K x J H 7 Q o 5 T 6 A v d j F j G C t b E C + 0 A E u U c G L H C L K 9 c 7 z T 0 p 0 C 1 t F Y F d d W r O R G g R 3 B l U Y a Z G Y H 9 5 v Y R k g s a a c K x U 1 3 V S 7 e d Y a k Y 4 L S p e p m i K y R D 3 a d d g j A V V f j 4 5 o E D H x u m h K J H m x R p N 3 N 8 T O R Z K j U R o O g X W A z V f G 5 v / 1 b q Z j i 7 9 n M V p p m l M p o u i j C O d o H E a q M c k J Z q P D G A i m f k r I g M s M d E m s 4 o J w Z 0 / e R F a Z z X X 8 P 1 5 t X 4 9 i 6 M M h 3 A E J + D C B d T h D h r Q B A I F P M M r v F l P 1 o v 1 b n 1 M W 0 v W b G Y f / s j 6 / A G a 8 p X A < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " B S Y A T f d a l l E m a 7 t 9 Z X o 2 w c 8 2 q v I = " > A A A C A H i c b Z D L S s N A F I Z P 6 q 3 W W 9 S F C z e D R X A h J R F B N 0 L R h S 4 r 2 A s 0 I U y m k 3 b o T B J m J k I J 2 f g q b l w o 4 t b H c O f b O L 0 s t P W H g Y / / n M O Z 8 4 c p Z 0 o 7 z r d V W l p e W V 0 r r 1 c 2 N r e 2 d + z d v Z Z K M k l o k y Q 8 k Z 0 Q K 8 p Z T J u a a U 4 7 q a R Y h J y 2 w + H N u N 5 + p F K x J H 7 Q o 5 T 6 A v d j F j G C t b E C + 0 A E u U c G L H C L K 9 c 7 z T 0 p 0 C 1 t F Y F d d W r O R G g R 3 B l U Y a Z G Y H 9 5 v Y R k g s a a c K x U 1 3 V S 7 e d Y a k Y 4 L S p e p m i K y R D 3 a d d g j A V V f j 4 5 o E D H x u m h K J H m x R p N 3 N 8 T O R Z K j U R o O g X W A z V f G 5 v / 1 b q Z j i 7 9 n M V p p m l M p o u i j C O d o H E a q M c k J Z q P D G A i m f k r I g M s M d E m s 4 o J w Z 0 / e R F a Z z X X 8 P 1 5 t X 4 9 i 6 M M h 3 A E J + D C B d T h D h r Q B A I F P M M r v F l P 1 o v 1 b n 1 M W 0 v W b G Y f / s j 6 / A G a 8 p X A < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " B S Y A T f d a l l E m a 7 t 9 Z X o 2 w c 8 2 q v I = " > A A A C A H i c b Z D L S s N A F I Z P 6 q 3 W W 9 S F C z e D R X A h J R F B N 0 L R h S 4 r 2 A s 0 I U y m k 3 b o T B J m J k I J 2 f g q b l w o 4 t b H c O f b O L 0 s t P W H g Y / / n M O Z 8 4 c p Z 0 o 7 z r d V W l p e W V 0 r r 1 c 2 N r e 2 d + z d v Z Z K M k l o k y Q 8 k Z 0 Q K 8 p Z T J u a a U 4 7 q a R Y h J y 2 w + H N u N 5 + p F K x J H 7 Q o 5 T 6 A v d j F j G C t b E C + 0 A E u U c G L H C L K 9 c 7 z T 0 p 0 C 1 t F Y F d d W r O R G g R 3 B l U Y a Z G Y H 9 5 v Y R k g s a a c K x U 1 3 V S 7 e d Y a k Y 4 L S p e p m i K y R D 3 a d d g j A V V f j 4 5 o E D H x u m h K J H m x R p N 3 N 8 T O R Z K j U R o O g X W A z V f G 5 v / 1 b q Z j i 7 9 n M V p p m l M p o u i j C O d o H E a q M c k J Z q P D G A i m f k r I g M s M d E m s 4 o J w Z 0 / e R F a Z z X X 8 P 1 5 t X 4 9 i 6 M M h 3 A E J + D C B d T h D h r Q B
Here, T d stands for the temperature of the dark photon bath, and 2,3 stands for the binding energy of H 12 and H 13 , respectively, given by µ α 2 d /2, with µ the reduced mass of the bound state in question. β 2,3 is the recombination rate, which relates to the ionization rate of the excited (n = 2) state α (2) 2,3 = 9.78 The factor C 2,3 takes into account Peebles' correction to the process, and its value can be approximated as with H being the Hubble rate. Λ 2γ D (2,3) stands for the two photon decay rates of H 12 and H 13 , which we estimate by rescaling the SM result Λ 2γ = 8.227 sec −1 by the ratio of the relevant binding energies. The Lyman alpha production rate is β 2,3 = β 2,3 e 3 2,3 /4T d . When calculating the ionization fractions, we calculate the Hubble expansion rate, keeping in mind that the dark and visible sectors have different temperatures, using the benchmark value ∆N ef f = 0.60 as discussed above.
In Fig. 4, we plot X 3 during H 13 recombination, for a few representative values of m χ 3 and α d , for n 3 /n 1 = 10%. The plots on the left and right correspond to the case where χ 1 is the heaviest flavor with m χ 1 = 25 GeV, and the case where it is the second heaviest flavor with m χ 1 = 1 GeV, respectively. The relic ionization fraction becomes larger either for a smaller value of α d or a larger value of m χ 3 . The residual ionization fraction can be approximated as [10] Dark baryon acoustic oscillations end after H 13 recombination, and the recombination time scale determines the suppression of the power spectrum. Since the H 12 recombination happens at a much earlier time, it does not have a strong effect on the power spectrum.

Dark acoustic oscillations and large scale structure
We now turn our attention to the evolution equations for the DM and dark radiation (DR) density perturbations in our model, after the H 12 recombination, but before the H 13 recombination. Thus the relevant matter degrees of freedom are χ 1 (with X 1 ∼ n 3 /n 1 during this epoch), χ 3 (with X 3 ∼ 1) and the bound state H 12 (with n H 12 ∼ n 2 ). Below, we will refer to these degrees of freedom collectively as "Acoustic DM" (AcDM) that undergo acoustic oscillations. We work in the conformal Newtonian gauge [49] where the fields ψ and φ describe scalar perturbations on the background metric. To linear order in the perturbations, we havė where the derivatives are with respect to η, the conformal time. δ ≡ δρ/ρ is the perturbation of the energy density, k is the wave number, and θ ≡ ∂ i v i is the divergence of the comoving 3-velocity. Since the AcDM components are all non-relativistic at this time, one can ignore the sound speed 3 . We consider the parameter region (α d > 0.02, m χ 3 < 10 MeV), where all AcDM components remain coupled during the dark acoustic oscillations. As long as the momentum transfer rate from the dark Thomson scattering γ D χ 3 → γ D χ 3 is comparable to Hubble, the density perturbations oscillate with the dark photon perturbation, and structures cannot grow. The cross section is given by σ T = 8πα 2 d /3m 2 χ 3 , and n f 3 depends on the ionization fraction X 3 obtained by solving Eq. (3.1). The dark photon perturbations, including higher modes in the Legendre polynomials, F γ D , evolve as [49] 10 a n f 3 σ T F γ D 2 , (3.10) Here the F γ D l are related to the spatial variations in the density fluctuations in the dark photons, in particular δ γ D ≡ F γ D 0 , θ γ D ≡ 3 4 kF γ D 1 /4, and σ ≡ 1 2 F γ D 2 where σ is the shear stress. We truncate the Boltzmann hierarchy at order l max = 4, making use of the approximation outlined in Ref. [49] 4 . The equations are similar to those for the SM photon and baryons.
In the calculation we take ψ = φ and ignore the correction from free streaming radiation. This approximation is good since our ∆N ef f is much smaller than the number of light neutrinos. Gravity perturbations are sourced by the density fluctuations as described by the Einstein equation, where the sum is over the SM photon, the dark photon and the AcDM components. For the initial conditions, the modes that enter the horizon before matter-radiation equality satisfy 14) and the modes that enter during the era of matter domination satisfy We set the initial values of the higher modes F γ D ≥2 = 0, since these higher angular modes quickly damp away when the AcDM-γ D scattering is efficient. We neglect the tilt in the primordial spectrum (n s = 1) and take a k-independent value of ψ = 10 −4 . The final results are independent of the precise value of ψ since we are interested in the ratio of the matter power spectra with and without the dark acoustic oscillations. In the numerical study, we choose the values h = 0.67, Ω γ h 2 = 2.47 × 10 −5 , Ω Λ h 2 = 0.69, Ω b h 2 = 2.2 × 10 −2 , and Ω ν = 0.69Ω γ [41]. After solving this set of differential equations, in order to quantify the importance of dark acoustic oscillations, we compare the DM power spectrum of AcDM to that of collisionless DM with an added non-interacting dark photon component, such that the energy density of dark radiation is equal in both scenarios, and further complications in the fitting of cosmological parameters are avoided: where the terms on the right hand side refer only to the nonrelativistic DM component. We show the power spectrum ratio in Fig. 5 for two representative values of α d in the region of interest. When the H 13 recombination takes place, the momentum transfer term in Eq. (3.7) quickly drops blow the Hubble expansion rate. The density perturbations entering the horizon after this point evolve the same way as they would in the ΛCDM scenario, and thus the ratio for small k modes asymptotes to 1. The matter power spectrum receives a suppression for modes that enter the horizon before recombination, thus for a lower H 13 binding energy (blue curve) there is a larger suppression. This helps to explain the results of low red-shift measurements for σ 8 . We estimate the viable parameter region by requiring a 5 − 15% suppression of the power spectrum at k = 0.2 h Mpc −1 (blue band in figure 5 and yellow band in figure 6). Although in this work we only focus on the suppression of the matter power spectrum in the linear regime, the suppression continues past k = 0.2h Mpc −1 . Once the non-linear corrections to the density perturbation are included, the galaxy survey data and Lyman-α observations, which probe the matter power spectrum at even larger k-modes, can be used to further constrain dark acoustic oscillations [11,[50][51][52], and thereby the parameter space of our model.

Scattering between bound states and small scale structure
After halo formation, the scattering cross section between the non-relativistic H 12 bound states in the halo is roughly geometric in size, σ ∼ (μ 12 α d ) −1 . For the region of parameter space we are interested in, where the heaviest DM flavor has a mass of O(10) GeV and α d ∼ 10 −2 , the resulting cross section over mass ratio σ/m H 12 ∼ 0.1 cm 2 /g can be large enough to thermalize the bound states and change the DM density in the inner part of the halo. It was pointed out in Ref. [12] that dark hydrogen with a similar range of mass and couplings can explain the low DM density cores observed in small galaxies. Moreover, if the hyperfine splitting of the bound state normalized by the energy scale E 0 = α 2 d µ 12 is comparable to the kinetic energy of the bound state, the effect of the inelastic hyperfine upscattering is to enhance the cross section of H 12 self-interactions. The velocity dependence in this process makes the self-interaction cross sections in dwarf halos larger than that in cluster halos, giving the correct σ/m ratio to solve the mass deficit problem in galaxy clusters.
In order to show that the same mechanism also works in the region of parameter space that is of interest to us, we adopt the best fit value from ref. [12] E hf = 10 −4 E 0 (3.18) and we limit ourselves to the range of α d in Fig. 3 of [12], where both the dwarf and cluster data can be explained. Since we want the H 12 self-interaction to solve the small scale structure problem, we choose n 3 /n 1 = 0.1 as our benchmark value such that H 12 and not H 13 is the dominant component of DM, while the n 3 number density is not unnaturally small (see figure 1), and also not too small to maintain the thermal equilibrium in the early universe that is responsible for suppressing structure formation and addressing the σ 8 discrepancy.
Even though it constitutes a smaller fraction of the DM energy density, we need to assess whether the H 13 bound state may still play a role in halo formation. In particular, the H 13 bound state has a much larger radius than H 12 due to the smallness of m χ 3 , and therefore scattering with H 13 could potentially change the desired core size of the H 12 halo. However, we find that this is not the case. For the parameter range we consider, the inverse Bohr radius of H 13 (∼ 10 keV) is still too small compared to to the value that would result in sufficient momentum transfer for keeping H 12 atoms in thermal equilibrium both at dwarf galaxies (∼ MeV) and galaxy clusters (∼ 100 MeV). Thus the scattering process most efficient for momentum transfer is off the χ 1 particle (i.e. the "nucleus") inside H 13 , with a cross section comparable to the H 12 self-scattering. Therefore the geometric size of H 13 does not lead to an enhancement, and the H 12 isothermal profile is not significantly affected.
As structures form, the bound states fall into the overdense region and their gravitational potential energy is converted into kinetic energy, resulting in shock-heating to a temperature [53] T gal ≈ 0.86 keV µ 10 GeV (3.19) for a Milky Way sized galaxy with halo mass 10 12 M and radius 110 kpc. Here µ is the total mass of all degrees of freedom that contribute to ρ DM divided by their total number density, given by (3.20) If T gal is higher than the binding energy, bound states can dissociate. While the more tightly bound H 12 does not dissociate in the region of parameter space that is of interest to us, if H 13 "reionizes" in this fashion, the scattering process χ 3 χ 1 → χ 3 χ 1 γ D can lead to efficient cooling through bremsstrahlung. For simplicity, when checking for the ionization of H 13 , we take the initial condition to be n f i = 0. In the parameter region where H 13 reionizes, we then recalculate T gal with n f 1 = n f 3 = n 3 (fully reionized H 13 ) to use in the estimate for the cooling process. An estimation of the cooling time scale through this process is given by [53] t brem ∼ 6 Gyr 0.02 α d (3.21) The emitted dark photon can have a free streaming length much larger than the size of the halo, leading to halo cooling. If t brem is much shorter than the age of the Milky Way (T M W ), a dark disc may form. Recently, results from the GAIA survey [54] have been used to set an upper bound on the fraction of the DM that can be contained in a dark disc at ∼ 1% [55,56]. A detailed study of the cooling process and the merger history of sub-halos is beyond the scope of this paper, therefore as we scan through the parameter space of our model, we will use t brem /T M W as a conservative indicator of whether there is a significant probability of dark disc formation. In Fig. 6, the region where H 13 can be reionized due to shock heating is shown below the red-dotted curve, and the region where the condition for efficient bremsstrahlung cooling is satisfied is shown above the red dashed curve, resulting in the red shaded region where both conditions are satisfied and where a dark disc may form.

Reconciling the large and small scale structure problems
In figure 6, we show our combined results for several representative parameter choices, and as a function of m χ 3 and α d . As mentioned in the previous section, we fix n 2 /n 1 = 0.9. We consider 25 GeV and 45 GeV as the mass of the heaviest DM flavor, considering both possible hierarchies, m χ 1 > m χ 2 and m χ 2 > m χ 1 .
• Addressing the σ 8 discrepancy: We calculate the power spectrum ratio of Eq. (3.16) for k = 0.2h Mpc −1 , which is close to the perturbation mode for the σ 8 measurement.
The contours of the power spectrum suppression depend mainly on the H 13 recombination time scale, thus a constant power spectrum suppression traces α d ∝ (m χ 3 ) − 1 2 for a constant H 13 binding energy, but they do not depend strongly on m χ 1 . The interesting regions for addressing the σ 8 problem are shown by the yellow bands.
• Small scale structure: We fix the ratio of the hyperfine splitting to the ground state energy as in Eq. (3.18). The mass of the intermediate DM flavor (χ 2 for the upper row of plots, and χ 1 for the lower row) is then determined at each point. This is indicated by the frame labels on the right-hand side of each plot. The preferred α d interval for solving both the dwarf and cluster mass deficit problems from ref. [12] is indicated by the blue shaded bands.  (45) GeV. n 2 /n 1 is chosen to be 0.9 in all plots. In the upper row, χ 1 is the heaviest flavor, so its mass is fixed, while the mass of χ 2 varies according to the E hf /E 0 = 10 −4 benchmark we have adopted (as indicated by the frame labels on the right-hand side of the plots). In the lower row, the roles of χ 1 and χ 2 are reversed. As a function of m χ3 and α d , between the yellow shaded contours the matter power is suppressed by 5-15%, which may explain the smaller σ 8 value from late-time measurements compared to the value obtained by Planck. The blue-shaded bands show the preferred α d interval where the mass deficit problem from dwarf galaxies to clusters is addressed. The red shaded region is disfavored by dark disc constraints as it allows for efficient cooling of the DM during halo formation. The preferred parameter space is therefore given by the overlap of the yellow and blue bands that is outside of the red shaded region.
• Constraints from dark disc formation: As explained in the previous section, the formation of a dark disc is possible in the red shaded region, which is therefore disfavored.
In summary, the overlap region between the yellow and blue bands that is outside of the red region gives the most preferred parameter space for addressing structure formation puzzles at different scales. This favors the ranges α d ≈ 0.02-0.04, 20-45 GeV for the mass of the heaviest DM flavor, and the MeV scale for the mass of the lightest flavor. We reiterate that in this study we have taken a relatively simple approach to demonstrate that our model has the potential to solve the relevant structure formation problems; however a more careful study of the cosmological data and the Lyman-α constraints should be performed to fully establish this claim and determine the precise region in the parameter region where all conditions of interest are satisfied.

Conclusions
We have explored the cosmological and astrophysical implications of a model of Secretly Asymmetric Dark Matter, where flavor-by-flavor asymmetries are generated in the dark sector through the decay of heavy right-handed neutrinos, despite an exact gauged dark U (1). As a result, the total dark charge of the universe is always zero, and the DM flavors have opposite signs of the asymmetry, making it possible for bound states to form. When the heaviest dark matter flavor has a mass of O(10) GeV, the intermediate flavor has a mass of O(0.1−1) GeV, and the lightest flavor has a mass of O(1) MeV, with an α d ∼ 10 −2 and n 3 /n 1 ∼ 0.1, this model can address several outstanding puzzles. In particular, the dark photon as an additional degree of freedom helps resolve the discrepancy between CMB-based and low redshift measurements of H 0 , the resulting dark acoustic oscillations help address the σ 8 problem, and scattering between the bound states H 12 after halo formation, with an inelastic component through the hyperfine transition, helps resolve issues at the cluster and dwarf galaxy scales. The model is consistent with constraints from short distance physics such as bounds on millicharged DM.
We want to emphasize that while we have chosen the DM to have three flavors for simplicity, the SADM mechanism works for a larger number of flavors as well, resulting in potentially even richer dynamical phenomena due to the increased number of relevant energy scales. We leave the exploration of this possibility to future work. Note also that as experimental sensitivity to millicharged DM increases, especially as the mass threshold of direct detection experiments is lowered, the parameter region of interest to this study should become testable in the not too distant future. Also, if the mass of the mediator φ is smaller than the benchmark value taken in equation 2.7, the lifetime of the heavier DM flavors can enter the regime where decaying DM signatures might become observable.