Boosting Asymmetric Charged DM via Thermalization

We consider a dark sector scenario with two dark matter species with opposite dark $U(1)$ charges and an asymmetric population comprising some fraction of the dark matter abundance. A new mechanism for boosting dark matter is introduced, arising from the large mass hierarchy between the two particles. In the galaxy, the two species thermalize efficiently through dark Rutherford scattering greatly boosting the lighter dark matter particle, far above the virial and escape velocities in the galaxy, while the dark charge prevents it from escaping. We study the consequences of this scenario for direct-detection experiments, assuming a kinetic mixing between the dark photon and the photon. If the charged dark sector makes up 5% of the total DM mass in our galaxy and the mass ratio is between $10^3-10^4$, we find that current and future experiments may probe the boosted light dark matter for masses down to 100 keV, in a hitherto unexplored parameter range.


Introduction
Astrophysical and cosmological observations of galaxy rotation curves [1], colliding galaxy clusters [2] and the Cosmic Microwave Background (CMB) [3] all indicate that a significant portion of matter in the Universe is not made of baryonic matter, the particles of the Standard Model (SM). Thus far, the missing Dark Matter (DM) has only been observed indirectly by its gravitational influence on baryons, indicating that its non-gravitational interactions with the SM are, at most, very weak. Its underlying nature remains mostly unknown; we do not know its mass, whether it's made up of a single or of several different DM species, or how it interacts with itself or with baryonic matter.
Results from various measurements and simulations point towards several noteworthy astrophysical properties of DM: it forms a halo around our galaxy whose profile is usually parameterized as NFW [4,5] (other parameterizations are Burkert [4,6,7] and Einasto [8][9][10]) with a local density of ρ DM 0.3 GeV/cm 3 [11] at the Earth's galactic location. The DM velocity distribution is usually taken as a Maxwellian distribution, truncated by the galactic escape velocity v esc ∼ 10 −3 c [12], which produces an approximately constant velocity rotation curve as seen in observations.
In the last few decades the leading DM candidate has been the Weakly Interacting Massive Particle (WIMP), a single weak-scale particle weakly-coupled to the SM. For O(TeV) scale masses and couplings of the same order as the weak force, the correct relic abundance is reproduced in what has been colloquially referred to as the "WIMP miracle". Such particles are predicted by several well-motivated theories beyond the Standard Model (BSM), most notably in supersymmetry (see for example [13] for a review on WIMPs and the minimal supersymmetric Standard Model). On the experimental front -state of the art DM direct-detection experiments are taking data, with the main effort focused on searching for the nuclear recoil from a rare scattering of DM against some target material. The current nuclear-recoil direct-detection experiments are XENON1T [14,15], PICO [16], CRESST [17], DarkSide-50 [18], CDMSlite [19], PandaX-4T [20] and LZ [21], and more are planned for the near future, most notably XENONnT [22], SuperCDMS [23] and DAMIC-M [24]. This detection method is sensitive to heavier DM, in the WIMP range, while other experiments and analyses are based on electron recoil sensitive to lighter DM masses: XENON10 [25], SENSEI [26], DAMIC [27], DarkSide-50 [28], CDMSHVeV [29] and EDELWEISS [30]. For even lighter masses, a detection strategy utilizing the wave nature of DM [31] is employed in CAST [32], ADMX [33,34], CASPEr [35], MADMAX [36], IAXO [37] and ABRACADABRA [38].
The immense experimental effort of searching for WIMPs has resulted in no discovery, and a significant fraction of the WIMP parameter space has already been excluded. At this stage, it is therefore important to consider broader ideas for DM beyond the standard paradigms. One such idea is the boosted DM scenario proposed in [39], wherein the DM is boosted to velocities significantly higher than its virial velocity ∼ 10 −3 c . In [39] a second subdominant and lighter component of DM is introduced. The massive species annihilate to the lighter species, transferring their rest-mass energy into kinetic energy, and as a result the lighter species are produced at high velocities. The higher momentum of the lighter species allows them to pass the detection threshold in direct-detection, enabling these experiments to probe mass ranges typically outside their reach. Other mechanisms of boosting DM were studied in [40][41][42][43][44][45][46].
In this work we consider a new type of boosted dark matter, based on thermalization instead of the annihilation in [39]. Here, the dark sector will have a dark U (1) D gauge group and two fermions of opposite charges, a heavy X and a light . We assume an asymmetry in the dark sector, so that the fermion anti-particles are not present while the total charge is still zero because the fermion particles have opposite charge (similarly to the electron and the proton in the SM). The two components couple through dark-electric interactions and can therefore thermalize so that by the equipartition theorem, the light will be boosted to have similar kinetic energy as the heavy X. If the mass hierarchy of the light and heavy particles is large enough, the light DM may reach speeds far above the escape velocity while the dark-electric interaction with the heavy component will keep it bound to our galaxy. The cosmological history and possible signatures of a similar setting was analyzed in [47,48] within the context of the Mirror Twin Higgs model [49].
Our dark sector interacts with the SM through a small kinetic mixing of the dark and SM photons. The kinetic mixing naturally induces an electric milli-charge on the DM [50] while its dark EM charges are ∼ O(1) (see [51][52][53][54][55][56][57][58] for examples of other studies of millicharged DM). As in [39] the lighter but faster DM can pass the energy detection threshold of terrestrial experiments, and we will see that future SENSEI runs will effectively probe milli-charge values of boosted at a range of 100 keV-10 MeV unconstrained by any current data.

The Setup
For our boosted dark matter model we introduce a dark sector with a U (1) D gauge group kinetically mixed with the SM U (1) EM . The dark photon γ D is assumed to be massless, and the dark sector is comprised of two Dirac fermions, denoted X and , with opposite dark charges ±1. They are singlets under all other gauge groups. The Lagrangian is µν its field strength, and m i the mass of fermion i (i = X, ). g D is the dark elementary charge and is the mixing of the dark and SM photons. Assuming small , we can change the basis of the Lagrangian up to O( ) to decouple the dark and SM photons, which induces an electric charge on the dark fermions: with Q ≡ 2 α D /α QED the effective EM charge of the DM and α D = g 2 D /4π the U (1) D fine-structure constant. With a small , therefore, X and obtain milli-charges. We set α D to α QED from here on out for simplicity so Q = | |.
We will further assume both milli-charged particles (MCPs) have an asymmetric abundance where X and remain, whileX and¯ were annihilated. This asymmetry is similar to the SM baryon asymmetry, and can be produced with any of the many proposed baryogenesis mechanisms (see e.g. [59] for a review). We will assume that the MCPs are thermalized via the dark-electric interactions, and that X is heavy while is light and that their mass ratio is m X /m 1. When this mass ratio is large enough, the temperature of the X and particles in the halo will be larger than the binding energy of (X ) bound states, and hence the dark bound states which were formed during "dark" recombination are re-ionized.
The binding energy of hydrogen-like atoms is B X = m α 2 D /2 while the temperature is where v X ∼ 10 −3 is roughly the virial velocity in the Milky Way, so requiring that T X B X results in m X /m 10 2 . In this work we will consider larger mass ratios m X /m = 10 3 , 10 4 such that the dark plasma can be assumed to be completely ionized.
While there are strong constraints on the magnitude of DM self-interactions [53,54,60,61], the constraints disappear if the self-interacting component is less than 10% of the DM mass [47,55]. We will assume for the rest of the paper that the MCP plasma makes up 5% of the DM by mass. Our model also implicitly contains a standard cold dark matter (CDM) candidate χ, neutral under U (1) D , that accounts for most of the DM in the galaxy. We remain agnostic about its exact nature and the early universe dynamics which produce the abundances of the CDM and the charged dark sector. It should also be noted that, while there are claims that MCPs are evacuated from the galactic disk by supernova shock waves and prevented from returning to it by magnetic fields [62,63], this is not the case for our model since any energy gained by shocks is quickly turned into heat via their self-interactions [64,65].
The constraints and new signatures arising from the cosmological history in a similar scenario were analyzed in [47], where the equivalent X and are the mirror proton and the mirror electron of the Mirror Twin Higgs model [49], and the U (1) D is mirror electromagnetism. A viable cosmological history can be accommodated, with the main constraints stemming from the dark U (1) D contribution to N eff . Additionally, new signatures arise due to the dark acoustic oscillations, which arise similarly to regular baryon acoustic oscillations. The oscillatory evolution of the density perturbation modifies the matter power spectrum on large scales, with the main effect of suppressing structure formation. A comprehensive study of the cosmological signatures of this model was performed in [66].
Another interesting phenomena to consider is plasma instability, which can occur in much of the parameter space of MCPs and U (1) D mediators, including the massless dark photon with strong self-interactions in our model [54,67,68]. However, although the plasma almost certainly exhibits instability, its observational effects are poorly understood due to the highly nonlinear nature of instabilities. Simulations of dark plasma instabilities have been performed previously [69] in the context of galaxy cluster collisions, and may offer a viable explanation for some features of the DM distribution of some observations. Dark plasma instabilities offer an intriguing avenue of exploration for dark photon models, with potentially huge constraining power in much of the parameter space, but would require detailed numerical studies and simulations to make concrete statements about said models. We ignore these complexities here, and therefore we will refrain from deriving robust bounds on our scenario in this work.
We summarize the dark sector particle content in table 1.

Galactic MCP Distribution
In order to study the potential for detection of the boosted MCPs on Earth, we first analyze their distribution in our galaxy. The charged dark sector has a subdominant mass fraction with little effect on the dominant neutral CDM distribution, for which we take a standard DM mass NFW profile ρ NFW (r) in the Milky Way [4]. For the analysis here we assume the X and plasmas are thermalized and in Sec. 4.1 we will find the parameter space wherein this assumption holds. The typical mean velocity of the X particles is of order v X ∼ 10 −3 , the virial velocity induced by the dominant CDM potential. The hierarchical mass ratio between the massive X and light results in a significant boost of the particles, possibly to speeds above the escape velocity. We found in an explicit calculation that the particles are still bound to the galaxy by their U (1) D interaction with the X particles. We model the halo dynamics of the and X plasmas as ideal gases in thermal and hydrostatic equilibria, with the neutral CDM acting as a background gravitational potential, and neglect dissipation processes in the MCP plasma. We further discuss dissipation, and the regions of parameter space where it is important, in Sec. 4.3. Our analysis assumes spherical symmetry, and we leave the study of more complicated magneto-hydrodynamic effects that may break this symmetry for a future analysis.
The pressure P i of each MCP plasma component balances the gravitational and darkelectric forces acting on it. The hydrostatic equilibrium condition of each component and Gauss's law for the dark charge and gravity is described by the following set of equations: (m X n X + m n + ρ NFW )4πr 2 dr , (3.1) where Q tot (M tot ) is the total dark charge (mass) of the halo. The SM electric force contribution is negligible since the MCP's SM charge is suppressed by with respect to the dark charge. The CDM component dominates the halo mass M tot (r) M CDM (r) and hence we neglect the MCP's contribution to the gravitational potential Φ, given therefore by ∇Φ = GM CDM /r 2 .
Due to the long-range dark electric force, the X and s are separated from each other only by a Debye length, which is much smaller than the length scale of the galaxy. Thus the number densities of the light and heavy particles trace each other, effectively setting n X n (≡ n MCP ). This is the neutral plasma approximation, and we can use it to simplify Eq. (3.1) to a single hydrostatic equilibrium condition, where ρ MCP = m X n X + m n m X n MCP is the total MCP mass density. We have selfconsistently checked that Q tot 1 in our solution by extracting Q tot (r) from Eq. (3.1). We now solve Eq. (3.2) to find ρ MCP , using similar methods to those used in [48]. To fully solve the equation, we need to specify the temperature profile of the halo. As we show in Sec. 4.2 thermal conduction is effective in the relevant range of parameters, and the MCP halo is approximately isothermal. Using plasma neutrality and the ideal gas equation-ofstate P i = n i T i we find that the total pressure of the MCP plasma is twice the pressure of the X species, P X + P = 2P X 2ρ MCP θ X (where θ i ≡ T i /m i is the temperature divided by the mass). Then the density profile is given by where Φ is the gravitational potential of the CDM component in a NFW profile. The central density ρ iso MCP (0) is set by requiring that the total fractional MCP mass in the galaxy is f MCP ≡ M MCP /M CDM , which we will take as 5% in this paper. It is trivial to rescale the results of this paper for any other value of f MCP .
To find the temperature θ iso X , we assume that the initial conditions, and in particular, the total initial energy of the MCP halo, are independent of the U (1) D interaction. Therefore, to find this initial energy, we consider the hypothetical scenario where these interactions are turned off. Then the "collisionless" MCPs would be indistinguishable from the CDM component and would follow the same profile shape. Taking the usual assumption that the system is virialized we may use the virial theorem 2K vir + U vir = 0 to relate the average kinetic energy of the particles to their potential energy. For this, we define the virial temperature as The total energy is then In our scenario, the collisions of the MCPs redistribute this initial energy among themselves but the total energy is conserved. The potential energy of the MCP halo is The total kinetic energy of the isothermal MCP halo is 1 We find θ iso X by requiring energy conservation This fixes the MCP's mass density profile, which is shown in Fig. 1. In similar fashion one can also solve the case where heat conduction is inefficient and heat convection results in an adiabatic halo, whose density profile is also shown in Fig. 1. In Appendix A we present the exact analytical solutions calculated in this section as well as explicitly solve the hydrostatic equilibrium condition in Eq. (3.2) for the adiabatic halo. In the following, unless stated otherwise, the MCP number density and temperature are computed in the isothermal halo approximation.

Thermodynamic Considerations
In the previous section we calculated the MCP profile assuming thermal and hydrostatic equilibria. Here we investigate some of the assumptions we have previously made. Firstly, we have assumed the and X plasmas are thermalized, T = T X . This is determined by the thermalization rate between the two species. Secondly, we have taken the MCP halo to be isothermal. This assumption depends on the rate of heat transfer in the halo. Lastly we have neglected dissipation processes through which the MCP plasma is cooled that would drive the collapse of the MCP halo to discs and other possible structures. Such cooling effects can be ignored if their energy loss is small. Below, we examine the validity of all these assumptions, and find the parameter space in which they are justified. This parameter space is presented on the plot in Fig. 6. We further expand on the analysis in this section in Appendix B.

− X Thermalization
We first analyze the rate of thermalization between the and X components. The two plasmas exchange energy through repeated collisions of and X particles until their average kinetic energies equalize. We compare the thermalization timescale t th through these collisions to the lifetime of our galaxy t gal in the parameter range of interest. We follow the results of [70][71][72][73] to compute the thermalization time, identifying the electrons with our light MCP ( ) and the protons with their heavy MCP partners (X).
The energy exchange rate per unit volume of particles due to Rutherford scattering with a non-relativistic X is given by where σ T ≡ 8πα 2 D /3m 2 is the dark Thomson cross-section. ln Λ is the Coulomb logarithm, where Λ ∼ T /ω p with ω p the relativistic plasma frequency ω 2 p = 4π 3 n α D /T , the effective cutoff for the IR divergence of Rutherford scattering.
The relaxation frequency, the inverse of the thermalization time t th , is given by where ε is the average kinetic energy of particles normalized to their temperature, The relaxation frequency is then We first note that the thermalization rate is dependent on n X and hence varies across the halo, with the denser inner regions thermalizing faster. The Solar System's location in the galaxy is in this inner region, and it is therefore possible that the two species are thermalized locally in our vicinity even if in the outer regions their temperatures differ, e.g. as in the blue curve in Fig. 2 on the left. However, the mean-free-path (mfp) of the particles in the halo can be fairly large in this part of our parameter space (see Eq. (B.1) and Fig. 2 on the right), making it harder to analyze the location from whence the particles arrive to the Earth. To simplify the discussion, we make the conservative assumption that the thermalization of and X particles is determined by the rate averaged over the entire galaxy -which is always slower than the local rate.
To find the parameter space where our thermalization assumption is valid, we compare the thermalization timescale t th to the age of our galaxy t gal . We find that for any choice of the ratio m X /m our assumption is valid for low enough values of m and breaks down above some critical mass. Above this mass the two plasmas have not reached thermal equilibrium within the lifetime of the galaxy and have separate temperatures. In this scenario the particles can still be somewhat boosted by the heat transfer between the plasmas, but we leave the analysis of this case for a future work. In Fig. 6 we show these maximal allowed masses as solid purple lines for two choices of the mass ratio m X /m = 10 3 , 10 4 .

Thermal Conduction
To determine whether our assumption of an isothermal MCP halo is valid, we estimate the rate of thermal conduction in the galaxy. Heat diffuses from one region of space to another in a series of rapid particle collisions, exchanging energy in the process until they reach kinetic equilibrium. The heat flux density q is proportional to the temperature gradient q = −κ∇T where κ is the thermal conductivity coefficient. The effective thermal conductivity coefficient in plasma is given by [74] The internal heat energy density of the MCP plasma is Q = 5 2 n X T + 5 2 n T = 5n MCP T . 2 Energy conservation then implies the continuity equation ∂ t Q + ∇ · q = 0, from which the heat equation follows: As the temperature and number density change in time so does the total pressure, P X + P = 2n MCP T = 2 5 Q. The pressure fluctuations propagate at the speed of sound in the halo v s ∼ P/ρ ∼ √ θ X ∼ v X (since X dominate the MCP density), and pass through the entire halo ∼ R halo in about 0.5 Gyr t gal . Therefore, we demand hydrostatic it is where ∂T /∂t vanishes for the initial adiabatic halo. As the halo evolves the temperature profile changes with it, ∂T /∂t will no longer be zero and heat will begin to conduct there.
equilibrium at any given time with a quasi-static evolution due to the thermal conduction. We can then differentiate Eq. (3.2) with respect to time, In the last equality we plugged in Eq. (4.6). We solve Eqs. (4.6) and (4.7) to obtain ∂T /∂t , and find the thermal conduction timescale is given by As we now show, the MCP halo density and temperature are well approximated by an isothermal halo in our parameter range of interest, and the profile starts to deviate considerably from an isothermal distribution only when cooling also becomes important. To see this, we first note that in the absence of heat conduction, convection results in an adiabatic halo. To check whether thermalization occurs, we take this to be the halo's initial condition, and check if the rate of heat transfer according to Eq. (4.8) is fast enough to thermalize the halo within t gal . This rate is position-dependent and, in particular, faster in the inner parts of the halo (see Fig. 3). For low m masses we find that heat transfer is inefficient throughout and the halo remains completely adiabatic, and for high m masses heat transfer is efficient throughout and the entire halo is isothermal.
In the intermediate m mass range this distinction blurs, and we find that only beyond some critical radius r c heat transfer is inefficient (see e.g. the blue curve in Fig. 3). To estimate the effect of the outer halo, where heat conduction is inefficient, we crudely separate the halo into two -the isothermal inner halo and the outer regions of the halo that are still in the original adiabatic configuration. The crossover radius is taken to be r c where the thermalization time coincides with the age of the galaxy. The inner halo temperature  Figure 4. The mass-density ρ MCP (left) and the temperature θ X (right) profiles of the MCP of the fully-isothermal, fully-adiabatic and "mixed" halos. We assume r c = 40 kpc for the "mixed" halo.
and density distribution is solved using the hydrostatic equilibrium Eq. (3.2) as before (see Appendix A.3 for more details) and plotted together with the outer adiabatic halo profile in Fig. 4.
We will treat the plasma as isothermal if r c > 40 kpc, since for r c = 40 kpc the resulting local density and temperature in the mixed-halo coincide with the fully isothermal calculation to within less than 50% error, which corresponds to the right of the dashed purple lines in Fig. 6. We note that this crossover regime coincides with the regime where cooling becomes important (see in Sec. 4.3), and so requires a full analysis of both effects which is left to a future study. We will therefore treat the halo as isothermal in our parameter range.

Cooling
The MCP plasma is cooled via two main processes, bremsstrahlung emission and inverse Compton scattering with dark CMB photons -relics from the early universe whose temperature we will take as a free parameter. Bremsstrahlung is the process of particles accelerating due to collisions with X particles and radiating away some of their initial energy. Inverse Compton scatterings are elastic + γ D → + γ D collisions, where the particles transfer energy to the dark CMB photons. In both processes the kinetic energy of particles is lost to the dark photons, cooling the MCP plasma. We investigate in what regions of the parameter space MCP cooling is important.
The rate of energy loss per unit volume due to bremsstrahlung is [75] P brems = n X n 16 3 whereḡ ff is the thermally averaged free-free Gaunt factor, which encodes quantum corrections to this formula. In our thermalized scenario the dark photon energy is roughly ω ≈ m v 2 ∼ T , so the classical limit whereḡ ff = 1 is sufficiently accurate [76]. The formula above is computed in the non-relativistic limit, and relativistic corrections are of order θ 3/2 [77]. In the parameter space of interest θ 3 · 10 −3 , and these corrections are subdominant.
Since the inner regions of the halo are denser than its outer regions, they will lose energy and cool down faster. Heat conduction then will kick in and thermalize the halo, 3 and the effects of cooling will be felt throughout the isothermal region of the halo. To neglect its effects on the local MCPs we wish to observe, the energy lost to cooling by the isothermal region of the halo must be small. We therefore can neglect the effects of cooling if the total energy lost within the lifetime of the galaxy is less than some fraction of this energy, chosen here to be 1/3.
The conduction process is effective in the inner region within r c > 40 kpc (see Sec. 4.2). The total energy of this region is E iso (r c ) = U iso (r c )+K iso (r c ) given in Eq. (A.5) in App. A.3. The total energy loss rate of this region of the MCP halo due to bremsstrahlung is then 10) and we estimate that within the lifetime of the galaxy a total of E cool ≈ ( dE/dt ) brems t gal was lost. Imposing the condition |E cool /E iso (r c )| 1/3 we find that cooling effects on the MCPs can be neglected for masses above (4.11) We plot the minimal m values allowed by cooing in Fig. 6 as solid blue lines for two values of the mass ratio m X /m = 10 3 , 10 4 . The MCP plasma is also cooled via inverse Compton scattering. Its timescale is controlled by the current dark CMB temperature T 0 D and the redshift z = 2 at which the Milky Way and its halo were starting to form. The energy loss rate per unit volume by inverse Compton scattering is [75] 12) and the total energy loss rate is The ratio of energy loss rate via bremsstrahlung to inverse Compton scattering is thus (4.14) For T 0 D ≤ 0.4T CMB more energy is lost to bremsstrahlung than to inverse Compton scatterings for all mass ratios of interest. Larger values of T 0 D > 0.5T CMB are excluded by ∆N eff constraints [3]. We thus ignore Compton scattering in the following discussions.
The dark sector in our scenario is similar to the SM plasma, which we know cools efficiently to form the galactic disk which is the Milky Way. Indeed, for the SM parameters m = m e , m X = m p we see from Eq. (4.11) that cooling is important in the SM. Additionally, cooling due to compton scattering with the CMB photons is significant, while dark compton scattering with the dark CMB is constrained by N eff . It is important to note that as the SM plasma cools down to lower temperatures, more cooling channels become available via atomic processes (e.g. recombination, collisional ionization and collisional excitation). In the dark sector we only consider bremsstrahlung which dominates at high temperatures compared with the Rydberg energy [75]. As we are only interested in determining whether cooling is important, it's enough to check this at the initial high temperatures.
We further note, that these calculations assume the energy lost to the dark photons isn't reabsorbed back into the plasma in photon-particle collisions. We verify this by computing the optical depth of a dark photon travelling from the center of the halo to its edge (in a straight line for simplicity), which is largest for γ D − scatterings, When τ γ D 1, the photons escape the galaxy and carry away energy from the plasma, so it is cooled efficiently. Otherwise, the MCP plasma still loses energy via the dark photons, but only from the outer layer of the galaxy, resulting in much more complicated cooling dynamics which we leave for a future study. We determined that reabsorption of dark photons in the halo becomes relevant for masses satisfying τ γ D ≥ 1. For our two choices of the ratio m X /m = 10 3 , 10 4 , we find that γ D escape the halo for all masses above {31, 14} keV, as shown by the dashed blue lines in Fig. 6.

Direct-detection Prospects
In this section we calculate the experimental signatures of the light MCP in electronscattering based direct-detection experiments. These experiments detect ionized electrons from a target material on Earth, searching for an excess of events above the expected background rate. The excess electrons were scattered off the target atoms by galactic DM. The incoming DM must surpass a momentum threshold to transfer enough energy to the scattered electron for it to overcome the atomic binding energy. It is typically assumed the DM has virial velocity ∼ 10 −3 , which limits the experiment's sensitivity to low DM masses. This is not the case in our scenario where the 's speed is above the virial velocity, which allows low-mass particles to pass the detection threshold. This is a generic prediction of boosted DM -direct-detection experiments are able to probe boosted DM in lower mass values than in the standard virialized DM scenario.
The incoming scatters off a bound electron in the detector and excites it, the kinematics of which we solve in Sec. 5.1. This scattering is governed by the matrix element, which we calculate in Sec. 5.2. We then compute the rate of − e − scattering events expected to be measured in electron-recoil based direct-detection experiments such as XENON10 and SENSEI in Sec. 5.3.

Bound Electron-Relativistic MCP Kinematics
The light MCP travels through the galaxy and happens to hit the target material in our terrestrial experiment, scattering off an electron of a target atom. The MCP transfers energy and momentum to the scattered electron, and may excite the electron from one energy level to another. We call the initial (final) electron state energy level 1 (2), whose identity will depend on the experiment. To derive the cross-section and event rate, we first consider the kinematics of this process. Derivations of the event rate where the DM is relativistic have been carried out in the literature before, which we review here. In our model the milli-charged scatters off the bound e − via photon exchange, as depicted in Fig. 5. We denote p ( ) = E ( ) , p ( ) as the incoming (outgoing) momentum, and β ( ) its respective velocities (with γ ( ) the Lorentz factor). We denote k as the electron momentum and E e,1(2) its total energy minus its mass in level 1 (2). The typical velocity of a bound electron is v e ∼ Z eff α/n (where Z eff is the effective charge felt by the electron and n is its principle quantum number), so we treat the electron nonrelativistically. q is the momentum transferred through the photon. In summary the fourmomenta are where ∆E 1→2 = E e,2 − E e,1 . They are related by momentum conservation

Off-shell Matrix Element Approximation
We turn to the matrix-element calculation using the Lagrangian in Eq. (2.2). The tree-level amplitude of an electron e − scattering with a particle of charge −Qe corresponding to the diagram in Fig. 5 reads The MCP is not polarized, so we average over the MCP spin states. We are blind to the scattered electron's spin, so we sum over its spin states. The square of the matrix element in Eq. (5.3) is 1 4 spins |M| 2 = e 4 Q 2 4q 4 tr / k e + m e γ µ (/ k e + m e )γ ν tr / p + m γ µ / p + m γ ν = 8e 4 Q 2 q 4 k e · p (k e · p ) + k e · p k e · p − k e · k e m 2 − p · p m 2 e + 2m 2 e m 2 .

(5.4)
If the electron were free, both particles would be on-shell and the expression above may be further simplified using Mandelstam variables. In direct-detection experiments the electron is bound and therefore off-shell, so we refrain from this simplification. Instead, we use momentum conservation, Eq. (5.2), to substitute p and k e with q in Eq. (5.4): (5.5) The dot products qp , k e p and k e q appearing in this expression depend on the angles between the 3-momenta. We simplify our calculation by removing this dependence as much as possible: since is on-shell we can use Eq. (5.2) to show that q 2 = 2qp . The electron is non-relativistic and its rest mass dominates the 3-momentum |k e | E e . After these simplifications we find We may disregard the term proportional to k e · q in k e q since its contribution to the cross-section vanishes once we integrate over all q directions. This is because the detectors considered are insensitive to the directionality of the scattered particles and their atomic form factors are spherically symmetric (see Eqs. (5.20) and (5.21) in the next section). After removing this term, we plug in the relevant momenta specified in Eq. (5.1) and derive the effective matrix-element needed in our calculations: (5.7)

Scattering Rate Formulae with Relativistic MCP
Having calculated the matrix-element in Eq. (5.7), we use it to compute the scattering rate of MCPs in electron-recoil based experiments. First we find the cross-section of e − + → e − + where the electron is bound. The standard cross-section calculation assumes free asymptotic states and needs to be amended when bound electron sates are considered. We follow the derivation presented in Appendix A of [78] to account for the effects of the atomic physics, and extend it to incorporate relativistic MCPs. Our results are equivalent to the derivation in [43]. Our starting point is the general cross-section formula in Eq. (A.7) of [78]: is the total energy of the initial (final) state of the system. f 1→2 (q) is the atomic form factor, i.e. the wave-function overlap between the initial (energy level 1) and final (energy level 2) states of the electron. It encapsulates the effect of the target material on the scattering event.
The electron is non-relativistic, so its energy is E ( ) e = m e , and we keep the Lorentz factors of the MCP's energy E ( ) = γ ( ) m . We obtain To calculate the rate of events we average over the MCP velocity and multiply by the local MCP number density where g (β ) is 's velocity distribution. We need eliminate the delta function δ(E f − E i ), that is, to impose energy conservation E f = E i . We make use of the kinematics specified in Sec. 5.1: from Eq. (5.1) we find the total energies of the initial and final states to be E i = E e + E = m e + E e,1 + γ m , E f = E e + E = m e + E e,2 + γ m , (5.11) which implies the difference is We extract γ from 3-momentum conservation, where θ is the angle between q and p . The energy difference in Eq. (5.12) is then We find the minimal velocity required to excite the electron by solving the energy conservation equation E i = E f in Eq. (5.14) when the MCP scatters co-linearly, cos θ = 1.
In our model, the lower integration limit for the transferred 3-momentum is |q| min = ∆E 1→2 , which by Eq. (5.15) corresponds to v min (|q| min ) = 1. In contrast, in the standard virialized DM scenario the particle's maximal speed is the escape velocity v esc and so the excitation velocity must be smaller than it, v min ≤ v esc , which imposes the lower integration limit is instead |q| vir min ≈ ∆E 1→2 /v esc ∼ 10 3 ∆E 1→2 [78]. Since in our model we integrate over a larger phase space and smaller momentum transfers q (which dominate Rutherford scatterings, see Eq (5.7)), the cross-section and rate will be larger and will result in better detector sensitivity compared to those derived in the standard virialized DM scenario.
We eliminate the delta function in Eq. (5.10) with the cos θ integral using Eq. (5.14), which results in the scattering rate This is consistent with Eq. (A.14) of [78] when γ = 1.
Our choice of g (β ) is the isotropic Maxwell-Jüttner distribution, This is the relativistic extension of the Maxwell-Boltzmann distribution usually taken for galactic DM. Here, however, it is not truncated by the DM escape velocity, since particles are still bound to the galaxy even with relativistic velocities. Additionally, we do not account for the relative motion of the Solar System and the MCPs which introduces anisotropy in the distribution function. This is justified because v ≈ 8 · 10 −4 , way below the typical speed of particles we consider. With this choice of 's velocity distribution we find the scattering rate master formula where γ min ≡ γ(v min ).
Our discussion so far made no detailed assumptions about the experimental apparatus at hand -we did not specify what the electron's initial and final states are. The effects of the underlying atomic physics of the detector were quantified in the general atomic form factor f 1→2 (q) introduced in Eq. (5.8). Now we apply our result from Eq. (5.19) for the XENON10 and SENSEI experiments. The derivation is identical to the one in [78], so we shall not repeat it in detail here.
The XENON10 experiment uses a liquid xenon target, searching for electrons ionized from N T isolated atoms. The atom is assumed spherically symmetric with filled electron shells. Before the scattering event, the electron is bound to the atom with binding energy E B . The electron ionized by the MCP is treated as a continuum of positive-energy bound states, approximated as a free particle at asymptotically large radii. Its 3-momentum is k e with recoil energy E R = |k e | 2 /2m e . The ionization rate of an isolated atom is given by where |f ion | 2 is the ionization form factor and ∆E 1→2 = E B + E R . This is equivalent to Eq. (A.21) in [78].
The SENSEI experiment uses semiconductor crystals as the target material. The periodic lattice has a continuum of electron energy levels which form a band structure. The valence bands and conduction bands are separated by an energy gap, and as electrons are excited across the bandgap they create mobile electron-hole pairs which can be detected. The MCP hits one of N cell cells in the crystal and deposits E e energy to an electron which is then excited across the bangap. SENSEI measures the number of electron-hole pairs created by the event, with 1-3 detected pairs set as the threshold. The excitation rate in a semiconductor crystal is given by where |f crystal | 2 is the crystal form factor and ∆E 1→2 = E e . This is equivalent to Eq. (A.32) in [78].

Detector Reach and Predictions
In this section we present the results of our analysis of the direct-detection prospects for our boosted MCP. For this, we compute the expected detection rate of the milli-charged particle in these experiments using Eqs. (5.20) and (5.21) following the methods outlined in [25,26,78]. The projected future reach of SENSEI with the single-pair threshold is  [79], SN1987A (gray) [80], XENON10 and SENSEI constraints on virial milli-charged DM, assuming the same n as in the boosted MCP (light green and green) [25,26], the future projection of SENSEI in the virial case (dotted green) [78] and nuclear-recoil direct-detection constraints on the X partner (pink) [81,82]. See the text for more details.
shown in the dash-dotted red curve in Fig. 6, assuming a MCP detection rate of 3.6 signal events in 1 kg − yr (the optimistic benchmark value used in [78]). We see that a large part of the parameter space accessible by the future runs of SENSEI is not excluded by the current data of XENON10 [25] and SENSEI [26] (3-pair threshold 4 ), which probes the region above the solid orange and red curves, respectively. Our analysis of the MCP halo dynamics is simplified; we did not take into account dissipation, different plasma temperatures, the effects of interstellar magnetic fields or plasma instabilites. We therefore refrain from interpreting the solid orange and red curves as direct bounds, pending a more detailed analysis. In green we show the constraints for a standard virial milli-charged DM for the same experiments [25,26,78], assuming the same number-density as in our model, for comparison. As expected, the direct-detection experiments are able to probe boosted milli-charged DM in a significantly lighter mass range and with better sensitivity compared with the standard virial case.
We move to analyze the constraints on our scenario. In the pink region labeled "NR DD" we plot the nuclear-recoil direct-detection bounds on the heavier MCP X. This species is non-relativistic and is subject to the usual constraints set by these experiments. We translate the bounds from the m X − Q plane to the m − Q plane, keeping the mass ratio fixed at m X /m = 10 4 and m X /m = 10 3 , respectively. The constraints on X were obtained by XENON1T [14], PandaX-II [83], LUX [84], PICO [16], CRESST [17] and CDMSlite [19], derived in [81,82] for light mediator exchange. We also include indirect constraints from BBN/CMB. Charged DM and massless dark photons contribute to the radiation density of the Universe at matter-radiation equality to which the CMB spectrum is sensitive, and the effect is usually parameterized by the effective number of neutrinos N eff . The additional radiation also changes the relic abundance of primordial elements generated by BBN. The bounds on N eff from CMB and BBN are translated into bounds on the MCP mass and charge in [79]. These are shown in blue and light blue contours in Fig 6. Inside stars MCPs are produced in pairs by plasmon decay and contribute to the total energy loss of the star. The MCP cooling channel modifies stellar evolution and shortens the lifetime of stars, thereby changing their expected population. The observed population of red giants, horizontal-branch stars and white dwarfs therefore constrain these additional channels of stellar cooling [79]. A similar argument can also be applied to the supernova SN1987A: the expected number of neutrinos would decrease if the proto-neutron star is cooled by more particles, which can be compared to the observed flux of the supernova [80]. Recently, however, questions were raised with regards to how robust these bounds are [85]. Stellar cooling bounds are shown in light orange and bounds from SN1987A are shown in gray in Fig 6. We now consider the regions of parameter space wherein our assumptions regarding the MCP's distribution and speed are justified. The solid purple line indicates the allowed region in which the X and plasmas are thermalized, as explained in Sec. 4.1. To the right of the purple line the two plasmas have different temperatures and transfer heat slowly. It is plausible that close to the thermalization bound the plasma is heated enough to boost the particles as in the thermalized scenario we considered, but we leave the full analysis for future work. The solid blue line indicates the boundary of the region in which the bremsstrahlung emission rate is sufficiently high and cools the MCP plasma (Sec. 4.3), allowing for a formation of a dark disk within the galaxy, resulting in lower MCP velocities v X ∼ 10 −4 [48,[55][56][57][58]. In the cooled MCP region further investigation of the MCP halo dynamics is needed to draw definitive conclusions about our model. We thus refrain from plotting the detection reach in this region, although we do not rule out the possibility this region may still be within the reach of SENSEI. The dashed purple line indicates the crossover between the isothermal and adiabatic halo as dictated by the thermal conduction timescale (see Sec. 4.2). As we can see, the halo can be treated as isothermal wherever cooling is negligible. Beyond the dashed blue line the dark photon's optical depth is greater than one, so dark photons emitted in cooling processes are reabsorbed into the plasma, further complicating the dissipative dynamics. Additionally, throughout the paper we assumed the MCPs sustain a spherical halo. It is conceivable that due to gravothermal collapse the halo had already collapsed within the lifetime of our galaxy. In Appendix C we argue why the MCP halo will not undergo gravothermal collapse and is stable in our parameter region of interest.
Finally, we have presented the results for two values of m X /m . As m X /m increases, direct-detection experiments are able to probe smaller m masses in the parameter space. This is because the higher temperature T /m ∝ m X /m increases the average velocity of particles. With higher speeds, lighter MCPs are able to overcome the momentum threshold required to excite the bound electron in the detector (Sec. 5). For even higher values of m X /m , SENSEI would only be able to probe our model in an already excluded region due to stellar cooling bounds, while mass values above 10 −2 -10 −1 MeV would be beyond the thermalization limit, where the effects of heat transfer need to be understood in order to give meaningful predictions about our model.

Conclusions
In this work we explored a two-species MCP DM sub-component where a large mass hierarchy boosts the lighter particles due to thermalization with their heavy partners. The dark charge of the halo prevents the light MCP from escaping it despite its high velocity. We evaluated the MCP density in thermal and hydrostatic equilibria in the presence of a dominant neutral CDM component. We found the relevant parameter space where the MCP components are thermalized and dissipative dynamics can be neglected, and determined that the MCP halo can be treated as isothermal in our parameter range. We then analyzed the detection rate of the lighter species in direct-detection experiments XENON10 and SENSEI, concluding that the boosted MCP can be probed in lower masses as anticipated.
The light MCP mass range in our model is within 100 keV to 10 MeV. The upper mass limit of our analysis is set by the thermalization requirement, above which the two MCP plasmas are not in thermal equilibrium. Current experimental constraints access parts of the available parameter space, while large parts remain unexplored and could be probed by future runs of SENSEI. Much of our lower mass parameter space is subject to significant cooling, which should be considered in a future study. A departure from the isothermal halo is also expected in this region.
where x ≡ r/R H , with the density scale and scale radius taken to be ρ H = 1.4×10 7 M kpc −3 and R H = 16.1 kpc in the Milky Way [4]. Using the virial theorem we can find the virial temperature of the CDM in Eq. (3.4) to be The isothermal MCP halo distribution is found by solving the hydrostatic equilibrium condition in Eq. (3.2) along with the equation-of-state P i = n i T and is given by: 3) The central density is found by requiring the condition on the total mass of the MCP halo M MCP = f MCP M CDM , where M CDM is the total mass of the CDM halo and f MCP is taken as a free parameter. Consequently The isothermal temperature θ iso X is fixed by the requirement of energy conservation in Eq. (3.8). E vir is given in Eq. (3.5). The potential and kinetic energies of the isothermal halo were found in Eqs. (3.6) and (3.7), and are computed to be (A.5) By numerically solving Eq. (3.8) we find θ iso X 2.3 × 10 −7 .

A.2 Adiabatic MCP Halo
If heat transfer is inefficient the MCP halo develops a temperature gradient, and is instead adiabatically isolated. For ideal gases the polytropic process equation is P i = C i ρ γ i with adiabatic index γ and some constant C i , and combined with the ideal gas law P i = ρ i θ i this gives the equivalent relation θ i = C i ρ γ−1 i . The gas constituents are particles (X or ) without internal degrees of freedom or "monoatomic" where γ = 5/3. As in the isothermal regime in Sec. 3 the equal number densities (n X = n ) and temperatures (T X = T locally) with the forms of ρ iso , ρ adia and θ adia X given in Eqs. (A.3), (A.7) and (A.6) respectively. Since the halo's original configuration was purely adiabatic, we identify the region r > r c with the original adiabatic solution found in the preceding section, so the integration constantŝ C X and θ adia X (0) are the same as before. The remaining integration constants of the isothermal region, the central density ρ iso (0) and the isothermal temperature θ iso , are found by imposing the same constraints as in the previous two sections: (1)  = U iso (0, r c ) + U adia (r c , R halo ) + K iso (0, r c ) + K adia (r c , R halo ), (A. 13) where the potential and kinetic energies of the halos are given in Eqs. (A.5) and (A.9) and are understood to be computed only within the region indicated by the parenthesis. The isothermal, adiabatic and "mixed" MCP halo mass density and temperature profiles are shown in Fig. 4.

B Thermalization and Conduction Times of the MCP Halo
In Sec. 4 we examined some of the thermodynamic properties of the MCP halo, in particular its thermalization time t th and its conduction time t cond . Here we plot these timescales as a function of the halo radius for a few typical m masses in our parameter space. For comparison, on all the plots we also show the age of the galaxy t gal in blue.
The first plots (Fig. 7) compare the thermalization timescale t th with the age of the galaxy t gal as a function of the halo radius. The condition for thermalization is conservatively chosen to be t avg th t gal , where t avg th is calculated using the average number-density of the entire halo. With this condition, the orange curve shows the threshold value of the mass for which this average thermalization condition is saturated.
Next, we plot t cond vs. t gal as a function of r in Fig. 8. The radius at which t cond diverges is where ∂T /∂t changes sign in Eq. (4.8), which is an artefact of the arbitrary choice of defining the conduction timescale using the initial profiles, while in the dynamical process, the conduction timescale, as defined in Eq. (4.8) will change with time. We compare the two timescales at r = r c , which is chosen to be r c = 40 kpc (see the explanation in Sec. 4.2 below Fig. 4), and the profiles for the threshold mass value are again shown in the orange solid lines.
To complete our discussion, we also consider the mfp of the particles, relevant for the calculation of the thermalization rate in Sec. 4.1. We compute it using the self-collision time t [74] l mfp = v t = 1 0.714 where v = √ 3θ is the average speed of s. l mfp is plotted in Fig. 2 on the right.

C Gravothermal Collapse
Self-interacting dark matter (SIDM) halos can in principle form structure, such as in core formation (heat flow from the outer parts of the halo to the center) and core collapse (heat flow from the center outwards). The latter is known as gravothermal collapse, and one might worry if in the parameter space of interest of our model, the MCP halo we've assumed to exist had already collapsed on galactic timescales. In this section we explain why this is not the case. The key property determining core formation/collapse is the rate of heat transfer in the galaxy. It is characterized by the Knudsen number [51,86] Kn ≡ l mfp R halo , (C.1) where l mfp it the mfp of the DM. Kn 1 corresponds to almost collisionless dark matter, for Kn 1 heat transfer is effective and may lead to core formation or collapse, and for Kn 1 heat conduction is suppressed, inhibiting core formation/collapse. The lighter s are more mobile and transfer heat, so the relevant DM mfp is l mfp = l mfp in Eq. (B.1).
For most m masses considered Kn 1, so no core formation or collapse is expected. We find Kn ∼ 1 only for masses close to the thermalization limit where the thermalization timescale t th t gal . However, for DM with elastic Coulomb-like interactions the timescale of collapse is ∼ 330t th [87], which is longer than the age of the galaxy for these masses. We conclude that in the parameter space of interest gravothermal collapse is either inhibited or hasn't occurred yet.
Our discussion neglects the effects of dissipation, which is justified when the energy loss due to bremsstrahlung is small as explained in Sec. 4.3. When dissipation is taken into account it can accelerate gravothermal collapse [88]. A detailed analysis of halo evolution with dissipative processes is beyond the scope of this work.