Self-interacting dark matter without prejudice

The existence of dark matter particles that carry phenomenologically relevant self-interaction cross sections mediated by light dark sector states is considered to be severely constrained through a combination of experimental and observational data. The conclusion is based on the assumption of specific dark matter production mechanisms such as thermal freeze-out together with an extrapolation of a standard cosmological history beyond the epoch of primordial nucleosynthesis. In this work, we drop these assumptions and examine the scenario from the perspective of the current firm knowledge we have: results from direct and indirect dark matter searches and cosmological and astrophysical observations, without additional assumptions on dark matter genesis or the thermal state of the very early universe. We show that even in the minimal set-up, where dark matter particles self-interact via a kinetically mixed vector mediator, a significant amount of parameter space remains allowed. Interestingly, however, these parameter regions imply a meta-stable, light mediator, which in turn calls for modified search strategies.

Self-interactions among DM particles -provided that they are frequent enough-lead to a heat transfer that initially decreases the density contrast in the center of DM halos. As DM particles evacuate from the central regions, cuspy density profiles turn into cored ones, overall reducing the halo mass concentration and potentially explaining the 'mass-deficit' problem. Concretely, it is found that the required self-scattering cross section over DM mass, σ χχ /m χ , needs to be larger than 0.1−2 cm 2 /g at the scale of dwarf galaxies [39][40][41][42][43][44][45][46]. On the flip side, the non-observation of an offset between the mass distributions of DM and hot baryonic gas in the Bullet Cluster constrains the same ratio to σ χχ /m χ < 1.25 cm 2 /g at 68% CL [47][48][49], i.e., approximately 1 barn for a 1 GeV DM mass particle. This tension is further exacerbated by recent observations of cluster collisions, im- * nicolas.bernal@uan.edu.co † xiaoyong.chu@oeaw.ac.at ‡ suchita.kulkarni@oeaw.ac.at § josef.pradler@oeaw.ac.at plying σ χχ /m χ < 0.2 − 0.5 cm 2 /g [50][51][52]. Among particle physics models which can produce such large self-interaction cross sections are the scenarios where self-interactions are mediated by light states [37,[53][54][55][56]. Such scenarios have been tightly constrained by linking the experimental observations with the assumption that DM particles stay in thermal contact with the SM bath before they freeze out (e.g. [57][58][59][60][61][62]). However, this needs not to be the case, for instance if the DM abundance was generated through freeze-in [63][64][65], if the reheating temperature is lower than the DM mass [66], or if the thermal history of the universe was non-standard [67].
The purpose of this work is to sever this link and to be agnostic about the production mechanism for DM in the early universe. In fact, many alternative mechanisms to generate the observed DM abundance beyond thermal freeze-out mechanism exist, and even modifications to the latter are entirely possible. In fact, DM production in the context of non-standard cosmologies has recently gained increasing interest, see e.g. . Such scenarios enlarge the parameter space relative to standard freeze-out case, allowing, for example, for smaller annihilation cross-sections and/or higher DM masses.
In light of this, it is the purpose of this work to explore the self-interacting dark matter (SIDM) scenario under the following assumptions: 1 1. in lieu of firm knowledge of the thermal state of the early universe, i.e. for T 1 MeV, the DM production mechanism remains unspecified. However, 2. we require that at T 1 eV the DM abundance has matched onto the CMB-observed value, while acknowledging the possibility that additional restrictions in the window 1 eV < T < 1 MeV may apply. Based on these assumptions, and requiring that the selfinteraction cross section is of suitable strength, we obtain the parameter space that is allowed by current direct DM searches, and by astrophysical and cosmological observations once we couple the dark sector to the standard model (SM). It will be shown that there is a large viable range in the latter coupling, but generically such that it points towards a long-lived force mediator. This in turn, calls for specific search-strategies for indirect DM detection and provides a target for what has been termed the 'lifetime-frontier' in the search of new physics.
The paper is organized as follows: in Sec. II we introduce a simple, prototypical SIDM model featuring a fermionic DM candidate χ and a light dark photon mediator V . Imposing the required self-scattering cross section allows us in Sec. III to place constraints on the kinetic mixing parameter in the plane of DM and mediator mass, the main results of this paper. We conclude in Sec. IV.

II. SELF-INTERACTING DARK MATTER
For concreteness, in this work we consider the scenario of fermionic DM χ that is coupled to the SM through a dark vector portal Here, V µν and F µν Y are the field strengths of the dark vector V and SM hypercharge, respectively. The kinetic mixing with strength Y induces the interaction between charged SM fermions f and V via q f e(f γ µ f )V µ , where ≡ Y cos θ W and q f is the charge in units of the electromagnetic coupling e. As we will focus on mediator masses well below the electroweak scale, we may neglect its suppressed mixing with the SM Z boson. Instead of the DM gauge coupling g χ we will characterize the strength of DM interaction by a dark fine structure constant α χ ≡ g 2 χ /(4π). SIDM may alleviate small-scale structure problems if the self-interaction cross section per particle mass, σ χχ /m χ is in the ballpark of σ χχ /m χ = 1 cm 2 /g.
In this work, we use this as a requirement and fix the DM self-interaction to the value Eq. (2) at a relative DM velocity v = 10 km/s, typical for dwarf galaxies. As is well known, the DM self interaction cross section σ χχ has resonant structure [56], and multiple solutions to Eq. (2) in α χ exist. In the following, and throughout this work, we therefore take the minimum value of α χ that satisfies the equation above for each parameter set of (m χ , m V ) in this way. This choice generates the most conservative bounds on the parameter space, as solutions with larger α χ lead to stronger interactions between DM and SM particles. 2 As an example, Fig. 1 shows σ χχ /m χ as a function of α χ , for m χ = 10 GeV and m V = 0.1 GeV, for a relative DM velocity v = 10 km/s. The peak structure is characteristic of the quantum resonant regime where the scattering cross section has a non-trivial velocity dependence. For this benchmark, the corresponding solution adopted is α χ 1.3 × 10 −2 (dotted vertical line). Additionally, Fig. 2 depicts contour levels for the minimal values of α χ satisfying Eq. (2).
The value of α χ needs to further satisfy several conditions. First, the cluster bounds on SIDM constrain the cross section to σ χχ /m χ 0.5 cm 2 /g at v = 10 3 km/s. 3 Compatibility with Eq. (2) therefore requires the selfscattering cross section to be velocity dependent. We are hence in the light-mediator regime m χ > m V where the typical momentum transfer in the scattering exceeds the mediator mass. Second, we impose that the observed DM density in the centers of dwarf-sized halos  nihilation within their lifetime (∼ 10 10 yr), leading to at v = 10 km/s. This condition also guarantees that DM annihilation decouples at T 1 eV, so that the CMB observations of the DM abundance are not affected. Finally, we require α χ 10 as the perturbativity condition. In practice, this excludes DM masses above the electroweak scale for m V above tens of MeV. On the left side of the plane, the DM mass is at or below the GeV-scale. With decreasing m χ , i.e. with increasing occupation number of DM particles, the constraint on late-time annihilation ('stability') and self-scattering in clusters ('cluster bound') become more severe, implying a maximum mediator mass. Among the two, the cluster bound dominates, as the self-scattering cross section must be guaranteed to remain velocity-dependent, setting the most stringent limit on the value of m V . On the right side, the only constraint is from perturbativity. The sharp onset of the limit at m χ 250 GeV for m V 10 MeV is because the resonance self-scattering strength is insufficient to reach Eq. (2) with perturbative values of α χ .

III. EXPERIMENTAL CONSTRAINTS ON KINETIC MIXING
We now discuss the interaction of DM with the visible sector, induced by the Lagrangian term e(f γ µ f )V µ . In this section we show that if thermal freeze-out for DM production is not assumed, experimental bounds can be alleviated, and upper bounds on the kinetic mixing parameter can instead be extracted.

A. Scattering between DM and SM particles
For portal interactions induced by the mediator V , the spin-independent differential scattering cross section between DM and protons (i = p) or electrons (i = e) can be written in terms of where v is the DM velocity, and E R the kinetic recoil energy of target i of mass m i . Figure 3 summarizes the resulting direct detection limits from both DM-nucleus (left panel) and DM-electron (right panel) scatterings, which are discussed in detail below. Nucleon recoil direct detection. We constrain the size of the DM-proton scattering cross section σ p by taking into account results from four different experiments. Among these are the low threshold experiments CDMSlite with an exposure of 70 kg days, which constrains DM masses from 1.5 GeV onwards [92] and CRESST-II experiment with 52 kg live days, which constrains DM mass above 0.5 GeV [91]. For detectors with higher threshold, the LUX experiment with exposure of 3.35 × 10 4 kg day sets limits for masses above 7 GeV [93], while the latest XENON1T results from 1 ton year exposure limits DM mass above 6 GeV [94].
It should also be noted that the direct detection limits depend on the mediator mass. If m V 10 −3 m χ , the momentum exchange in the dark vector propagator is resolved. In the case of DM-proton collisions, the effect of the light mediator mass has been implemented by using the public tool DDCalc [99,100]. Also note that in the translation of general constraints on the DM-nucleon cross section σ n , the relation to σ p is σ n = (Z/A) 2 σ p where A (Z) are the atomic mass (charge) number of the target nucleus.
Electron recoil direct detection. For DM-electron scattering, we use the latest results from the SENSEI experiment [97] together with data from the XENON experiment. We further use the limits derived from reanalyzing XENON10 data [95,96], and the latest official XENON1T limits [98]. The limits are obtained using the so-called 'S2-only' analysis, where electroluminescence produces a secondary scintillation. The three limits cover complementary DM mass ranges. At every point in parameter space, we use the most stringent upper limit S t a b i l i t y P e r t u r b a t i v i t y Perturbativity  [91][92][93][94] in the left panel and DM-electron scattering [95][96][97][98] in the right panel.
on scattering cross section given the mass of the DM. The DM-electron scattering cross section is conventionally expressed as the cross section on a free electron at a reference momentum transfer where µ χe is the DM-electron reduced mass [101]. The heavy mediator limit m V q ref applies to most of the considered parameter region. For m V below 10 keV, a DM form factor should be taken into account in the standard calculation of recoil events [102]. Here, we instead conservatively assume the bound on at m V = 10 keV applies to smaller m V as well. DM-SM scattering in astrophysics. Beside the direct detection bounds, DM-SM interactions are especially constrained from astrophysical considerations. These include cosmic-ray attenuation by scattering with DM [103], CMB and large-scale structures modified by momentum-transfer with DM (e.g. [104][105][106][107][108][109]), among other probes. 4 Most of them are typically weaker than the bounds from nucleon (electron) direct detection experiments for 4 We do not consider the limit m V → 0, for which the resulting long-range force can be constrained differently; see e.g. [110][111][112].
GeV (MeV) scale DM. One example is the limit derived with direct detection and neutrino experiments, utilizing solar-or cosmic ray-upscattered DM [113][114][115][116][117]. For velocity-independent scatterings, the upper bound on σ χN can be as stringent as 10 −31 cm 2 for MeV DM [115,117], but becomes much weaker in presence of a light mediator [118]. This is because the upscattering process with cosmic rays is dominated by energyexchange much larger than that in direct detection experiments, and thus does not get enhanced even if the mediator particle is much lighter than the reduced mass of the system. Taking m χ = 10 MeV and m V = 10 keV, we have calculated that the upscattered DM flux spectrum is approximately taking into account electron, proton, and helium cosmic rays, which results in 10 −2 using the electron-recoil data from XENON1T; nuclear recoil data does not lead to a useful constraint.

B. Energy injection from dark particles
In the chosen scenario, the effects associated with energy injection are induced by the produced abundance of V particles and their subsequent decay with lifetime τ V . The dark vectors may e.g. be produced via freeze-in from the SM sector [119,120,126] or via DM annihilation [127]. In accordance with our assumptions made in  The various labels summarize principal detection strategies, astrophysical constraints from stellar energy loss ('stellar') [119], and diffuse γ-ray background ('diffuse') [120], from cooling of the proto-neutron star of SN1987A ('SN') [121], from beam-dump and collider experiments, see e.g. [122,123] and references therein. Hatched regions show recent exclusions derived from gamma-ray signatures from SN [124] (blue) and energy-transfer argument inside SN [125] (pink). In addition, contours of mediator lifetime τ V equal to 1 sec (BBN), 10 13 sec (CMB), as well as t U (universe) are shown. The vertical dashed gray line corresponds to m V = 2 m e . Additional cosmological constraints from BBN and CMB for m V > 1 MeV are not shown, as they rely on a V -abundance created at T > 1 MeV [126].
the introduction, we neglect any population of V that may have emerged prior to T = 1 MeV; we however do take into account an abundance of V that arises from (residual)χχ annihilation and from direct production of V particles from the SM bath below a photon temperature of 1 MeV.
Potential cosmological constraints on late V decay are governed by the V lifetime. They are located in the [m V , ] parameter plane along the contours of constant lifetime, shown in Fig. 4. Solid, dashed and dotted lines correspond to τ V = t U 4 × 10 17 sec (age of the universe), τ V = t CMB = 1.2 × 10 13 sec, and τ V = t BBN = 1 sec, respectively. For masses lighter than m V = 2 m e ∼ 1 MeV (vertical gray line), the mediator can only decay via (loop-induced) processes into three photons or two neutrinos. 5 The relevant bounds on 5 See [128] for a more precise calculation of the width for a dark photon decaying to three photons for m V ∼ 2 m e . dark photon production from the SM sector reproduced in Fig. 4 are taken from [119,[121][122][123][124][125]. Bounds from a freeze-in production of V at T > 1 MeV derived in [126] are not applied, in accordance with our assumptions.
The mediator abundance produced from DM annihilation for T ≤ 1 MeV is governed by the non-relativistic annihilation cross section of χχ → V V , where S(v) is the Sommerfeld enhancement factor [56,129], expressed as where κ = 1.6, consistent with our scattering cross section formula adopted from [56]. This enhancement only plays a role for m χ 10 GeV, and we have checked that a more careful treatment of the enhancement [130] does not affect our result qualitatively. Moreover, although radiative bound state formation can happen for m χ 1 TeV and m V ≤ 1 4 α 2 χ m χ [131], this process only improves the bounds mildly [132,133], and will not be further considered here. Equation (8) then allows to study the effects of energy injection from DM annihilation at different epochs, and to estimate the corresponding constraints.
Energy injection during the BBN epoch. Energy injection from dark particles after T 1 MeV may affect the primordial abundances of light elements, such as 4 He and deuterium, see e.g. [134] and references therein. If a population of χ-particles is already present at T = 1 MeV, their annihilation χχ → V V leads to an accumulation of V -particles, Y V | >1 sec = 2 1 sec dt σ an v Y 2 χ s, which -when they decay-inject energy into the primordial plasma; Y i ≡ n i /s where s is the entropy density. However, as is well known, BBN sensitivity falls short to probe a thermal annihilation cross section of 1 pb. In addition, the corresponding bound on s-wave (or velocityenhanced) DM annihilation is much weaker than that from CMB observations; see e.g. [135,136]. In a similar vein, the freeze-in of V -particles post T = 1 MeV will lead to a very weak constraint on that is already excluded otherwise. As a result, we will not further consider the effect of energy injection during BBN in our final bounds.
Energy injection during the CMB epoch. Here we interpret the Planck constraints on DM annihilation from [137][138][139] as upper bounds on the total electromagnetic (EM) energy injection at the time of recombination. 6 Similarly to the BBN case above, one may  start the investigation from the redshift-dependent production rate of the V abundance from DM annihilation, dY V (z)/dz, following it from an arbitrary earlier redshift, z I , down to the CMB time. In the calculation of the accumulated V abundance, we may then safely neglect mediator decays prior to the CMB epoch, since we are interested in parameter regions where the V lifetime inflight is considerably larger than t CMB . In this way, the EM energy density injected from Vdecays at z CMB is found from the redshift integral, where ρ χ,0 is the DM energy density at present,  [126]. In turn, for m V 4 keV, the mediator dominantly decays into neutrinos (Br EM 0), and therefore the bounds on energy injection vanish. The inequality sign in (10) is owed to the fact, that we have may also modify the CMB predictions, its contribution diminishes with time [142], and will not affect our bounds qualitatively. This can also be seen in Figs. 7 and 9 of [143]. taken the smallest, i.e. the relic value of the DM abundance for sourcing V through annihilation.
In the limit m χ m V and t(E z,dep ) t CMB , one finds that the injected EM energy density, E tot , is well approximated by the V abundance produced at the CMB time. This holds true even without any Sommerfeld enhancement of the annihilation cross section. Hence, in practice, the bounds on the case of mediator-decay with a finite lifetime can be obtained by rescaling from the existing constraints on WIMP DM annihilation. Therefore, we impose as requirement, where t * = t CMB corresponds to the characteristic time of the cosmological epoch, and t(m χ ) = (m χ /m V )τ V . For concreteness, we set the DM velocity to 10 −6 c and use for σ an v * the limiting value derived from Planck data [138,139]. As long as σ an v Br EM ≥ σ an v * , Eq. (11) implies an upper bound on , which is shown in the left panel of Fig. 5. Diffuse gamma-ray background. Similar to the CMB bound, we use Eq. (11) with t * = t U to constrain late-time DM annihilation through the isotropic extragalactic gamma-ray bounds obtained by Fermi-LAT for m χ ≥10 GeV [140,141]. Furthermore, it is conservatively assumed that the bound on DM annihilation cross section σ an v is independent of m χ below 10 GeV, while it is expected to become even stronger with the deceasing m χ (e.g. see Fig. 14 of [144]). We also set the SE factor to unity, to obtain the upper bounds on shown in the right panel of Fig. 5. We note in passing that a cosmologically long-lived mediator would strengthen the bound because of a correspondingly smaller optical depth of the gamma ray signal. Neglecting this effect yields a conservative limit. Finally, the diffuse gamma-ray bound on from the freeze-in of V below T = 1 MeV is taken from [120], and incorporated in Fig. 4 ('Diffuse').
Indirect searches from the galactic center. At last, for the heavy DM mass regime, where the isotropic extra-galactic gamma-ray bound weakens, measurements from of the galactic center by H.E.S.S. [145] become relevant. In the parameter region where σ an v Br EM is larger than the limiting value of the H.E.S.S. bound, we require c t(m χ ) ≥ 1 kpc to put an upper bound on . For larger in-flight decay lengths, the morphology of the spectrum is distorted [146], and thus the above bound cannot be applied directly. Overall, the ensuing limits are only competitive in the region of high mass m χ 100 GeV, and do not improve the diffuse gamma-ray bound above noticeably.
The derived upper limits on are shown in Fig. 5, where the left and right panels correspond to CMB and isotropic gamma-rays observations, respectively. Below the gray line the inequality in Eq. (11) is always satisfied and therefore no limits can be extracted. For m V below the horizontal dashed gray line (m V = 2 m e ), the mediator can only decay via loop-induced processes into three photons or two neutrinos. Both decay channels are rather suppressed, and therefore bounds coming from electromagnetic energy injection weaken significantly.

C. Potentially complementary observables
There exist other observables that can be used to constrain the parameter region of this scenario. While they currently do not lead to stronger bounds than the ones discussed above, they may become relevant in the future, or provide competitive bounds in other SIDM scenarios.
First, DM particles heavier than several GeV can be captured by the Sun and annihilate into light mediators, which escape and finally decay. Depending on the assumptions, the bounds from Fermi-LAT and HAWC vary between 10 −42 cm 2 and 10 −46 cm 2 on spindependent DM-nucleon scattering cross section for m χ ≥ 4 GeV [147][148][149][150][151]. 7 However, they only apply to a very narrow parameter region in m V , about a few MeV. This is due to the fact that for the relevant values of DMnucleon scattering cross section, the mediator either decays inside the Sun (m V 1 MeV) or becomes very long-lived (m V ≤ 1 MeV). 7 If SIDM captured by the Sun via self-scattering is taken into account, an enhancement of 10−10 3 may be achieved, potentially leading to stronger bounds [150,[152][153][154].
It is also possible to directly search for a boosted (meta-)stable dark mediator. Similarly to the boosted DM scenario from DM annihilation [155] or decay [156], a long-lived mediator produced relativistically by DM annihilation could induce a signal in a DM or neutrino detector. This possibility is mostly interesting for sub-MeV mediators with ensuing cosmological long lifetime. For thermal DM annihilation cross sections (∼ pb), [155] estimates the bound as which, in our case, is predominantly driven by the inelastic channel e+V → e+γ. Simple dimensional analysis of the relevant cross sections then reveals that this bound is evaded for ≤ 10 −4 − 10 −5 . In addition, if any of the dark particles is light enough, they can be produced in supernovae. Since this thermally produced population is much hotter than the halo DM, it could be observed with direct detection experiments [157]. Upper bounds on dark couplings can also be derived from dissipation in dark halos. Cooling of SIDM via inelastic scattering (and radiative bound state formation in a limited parameter subspace) may happen. After BBN, such processes are typically strongly suppressed due to the low DM relative velocities. At low redshift, it may affect the halo dynamics, but the bound turns out to be weak [158,159]. This can be understood from the SIDM requirement that each DM particle only scatters elastically O(1) times during the whole lifetime of dwarf halos, and the frequency of inelastic scattering is naturally much smaller. The enormous range in as seen from the color legend is owed to the fact that the direct detection and indirect/cosmological searches apply to largely complementary regions in parameter space with significantly different sensitivity to the value of . The limits are further improved once the bounds independent of the DM abundance (Fig. 4) are taken into account. This is especially true in the region of m V < 1 MeV and is summarized in right panel of Fig. 6.

D. Discussions on combined constraints
Our summary figures show that the minimal mediator lifetime that is allowed is greater than 1 sec, found in the region m V ∼ 10 MeV and ∼ 10 −10 . Furthermore, taking into account the newly-derived bound from [125] would require the mediator to have a lifetime longer than O(10 5 ) sec. Together with the BBN observations, it suggests that a thermalization between the dark and visible sectors at T ∼ MeV is very unlikely.
At last, it is worth pointing out that if DM is asym- metric, the above bounds from energy injection can be relaxed. However, a larger DM annihilation cross section may be needed to erase the symmetric DM component in early universe. This adds a layer of complication when the premise is to study SIDM with the least set of assumption about the early universe DM history, as spelled out in the introduction. We leave an exploration of this possibility for future work.

IV. CONCLUSIONS
In this work we analyze the parameter space for selfinteracting dark matter (DM) without assuming any specific mechanism for relic density generation and the thermalization between the dark and visible sectors. The motivation for it is that, currently, no firm conclusions can be drawn about the thermal state and the particle content of the early universe above a photon temperature of several MeV -despite a number of circumstantial hints of what the sequence of events may have been.
For concreteness, we analyze a symmetric fermionic DM candidate χ, where the self-interaction is mediated by a dark photon V that kinetically mixes with the SM counterpart with strength . Throughout the work we require a self-scattering cross section over DM mass of 1 cm 2 /g as the fiducial value that is able to address the astrophysical small-scale problems of ΛCDM. For chosen DM and mediator masses, this fixes the dark gauge coupling and completely determines the hidden sector parameters. The link to SM is controlled by the kinetic mixing parameter and the velocity dependence that is introduced in this model through Sommerfeld enhancement in DM annihilation makes it also interesting from the astrophysical point of view.
Despite relaxing the assumptions on the early universe history, a number of strong constraints remain applicable to this model. Instead of ruling out the model parameter space, these constraints imply upper limits on the value of and hence lower limits on the lifetime of V . The limits on are derived from direct and indirect detection as well as from cosmology. The direct detection of DM through nuclear or electron recoil signatures offer complementary limits on . While DM-electron scattering strongly constrains the low mass region (∼ m χ < 1 GeV), the nuclear recoil signature is sensitive to heavier DM (∼ m χ > 1 GeV). The constraints on the kinetic mixing from DM-nucleus scattering are comparable to those from beam dump experiments at the laboratory. The indirect searches and the CMB limits on the other hand are sensitive to a wide range of symmetric DM and mediator mass combinations and set some of the most stringent constraints for heavier mediators. While literature results are derived under the assumption of prompt energy injection, here we find that these limits need to be rescaled in order to take into account the mediator's finite (boosted) lifetime. For heavier mediators, the limits on arising from indirect searches usually supersede any laboratory limits and imply mediators to be (meta-)stable at a cosmological time scale.
Our analysis suggests a long-lived mediator, and thus the need for analyzing displaced vertices in the sky and small-scale anisotropies in cosmic rays from DM annihilation [146,160], instead of exclusively focusing on astrophysical objects that are DM-dense. For instance, a MeV-mass mediator with a lifetime of 10 10 sec ( ∼ 10 −14 ), produced in GeV-mass DM annihilation in the Galactic center, preferentially decays outside our Galaxy. Given that we have been agnostic about the DMgeneration mechanism, the 'model-independent' bounds on that we derive call for being further interpreted and tested in SIDM scenarios with spelled-out relic density mechanisms beyond the standard case. Such scenarios include feebly-interacting DM or DM in the context of a modified early universe thermal history.

ACKNOWLEDGMENTS
This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agree-ments 674896 and 690575. NB is partially supported by Spanish MINECO under Grant FPA2017-84543-P, by Universidad Antonio Nariño grants 2018204, 2019101 and 2019248, and by the 'Joint Excellence in Science and Humanities' (JESH) program of the Austrian Academy of Sciences. NB thanks the Institut Pascal at Université Paris-Saclay with the support of the P2I and SPU research departments and the P2IO Laboratory of Excellence (program 'Investissements d'avenir' ANR-11-IDEX-0003-01 Paris-Saclay and ANR-10-LABX-0038), as well as the IPhT. XC and JP are supported by the 'New Frontiers' program of the Austrian Academy of Sciences. SK is supported by Elise-Richter grant project number V592-N27. The authors thank the Erwin Schrödinger International Institute for hospitality during this work. This research made use of IPython [161], Matplotlib [162] and SciPy [163].