Full-shape BOSS constraints on dark matter interacting with dark radiation and lifting the S8 tension

In this work we derive constraints on interacting dark matter-dark radiation models from a full-shape analysis of BOSS-DR12 galaxy clustering data, combined with Planck legacy cosmic microwave background (CMB) and baryon acoustic oscillation (BAO) measurements. We consider a set of models parameterized within the effective theory of structure formation (ETHOS), quantifying the lifting of the S8 tension in view of KiDS weak-lensing results. The most favorable scenarios point to a fraction f ∼ 10-100% of interacting dark matter as well as a dark radiation temperature that is smaller by a factor ξ ∼ 0.1-0.15 compared to the CMB, leading to a reduction of the tension to the ∼ 1σ level. The temperature dependence of the interaction rate favored by relaxing the S8 tension is realized for a weakly coupled unbroken non-Abelian SU(N) gauge interaction in the dark sector. To map our results onto this SU(N) model, we compute higher-order corrections due to Debye screening. We find a lower bound αd ≡ gd 2/(4π) ≳ 10-8 (10-9) for dark matter mass 1000 (1) GeV for relaxing the S8 tension, consistent with upper bounds from galaxy ellipticities and compatible with self-interactions relevant for small-scale structure formation.


Introduction
Our knowledge about the constituents, structure and evolution of the Universe has increased considerably by combining precision data covering a broad range of length and time scales.Although the ΛCDM model still agrees with many datasets, precise measurements have unveiled tensions in the values of cosmological parameters estimated by distinct probes, possibly indicating unaccounted systematic effects or alternatively hints for deviations from ΛCDM.First, the Hubble constant H 0 determined from the local distance ladder [1] and inferred indirectly within ΛCDM from the cosmic microwave background (CMB) anisotropies [2] as well as the scale of baryon acoustic oscillations (BAO) [3] disagree at a level of around 4σ (see [4,5] for reviews and [6,7] for summaries of proposed solutions).Second, the value of S 8 = σ 8 (Ω m /0.3) 0.5 , with σ 8 being the amplitude of matter clustering on a scale of 8 Mpc/h and Ω m the matter density parameter, measured from a variety of large-scale structure (LSS) probes (including weak lensing [8][9][10][11], galaxy cluster counts [12,13], and lensing combined with galaxy clustering [14]) are systematically below the value inferred from CMB measurements within ΛCDM, with a significance of up to 3σ (see also [15,16]).Third, the so-called small-scale crisis poses an intriguing puzzle [17,18].
The aforementioned tensions can serve as a guideline to address another long-standing open question in cosmology: the nature of dark matter (DM).For example, the small-scale puzzles can be addressed in ultra-light [19] or self-interacting [20] DM models.On the other hand, the S 8 tension points to scenarios where the matter power spectrum on scales k ≳ 0.1h/Mpc is somewhat suppressed compared to larger scales of order k ∼ 10 −2 h/Mpc probed via CMB anisotropies.A priori, such a suppression can naturally arise in scenarios that deviate from the minimal paradigm of cold and collisionless DM.Examples include a nonzero free-streaming length for warm dark matter [21], a macroscopic de-Broglie wavelength such as for fuzzy dark matter [22], dark matter decay into invisible daughter particles [23], or certain self-interacting dark matter scenarios [24].Nevertheless, specific models imply particular patterns imprinted on the CMB and matter power spectra along with possible changes in the background evolution, such that it is non-trivial to identify viable and wellmotivated models successfully addressing the S 8 tension.
A mechanism for suppressing power that is known to operate in nature is the Compton scattering of baryons (i.e.visible matter) with CMB photons in the early Universe.It is a plausible possibility that an analogous mechanism is responsible for a power suppression in (part of) the DM component.For example, DM can be gauged under a symmetry group akin to the SM.If the symmetry is unbroken and the gauge coupling sufficiently small, the gauge fields are light or massless, and we end up with a dark sector with two components [25]: an interacting dark matter (IDM) coupled with a new component of ultrarelativistic dark radiation (DR).
The scenario of DM-DR interactions can be realized microscopically by a variety of setups, leading in general to different implications for the power spectrum and structure formation.On large scales, the effects of DM-DR models on structure formation can be parameterized through the effective theory of structure formation (ETHOS) [26].ETHOS provides a (Boltzmann) framework to follow the evolution of DM and DR perturbations in presence of DM-DR interactions as well as DR self-interactions.DM-DR models have often been discussed in the context of small-scale puzzles (see e.g.[27]) and also considered with respect to the H 0 and S 8 tensions [28,29].As far as the H 0 tension is concerned, they behave similar to models with extra radiation, meaning that they cannot solve the H 0 tension in a minimal setup (see however [30] using a two-component DR model and [31] assuming a phase transition in the DR sector that addresses H 0 ).However, the interaction of DM with DR allows to reduce S 8 relative to ΛCDM and therefore provides a possibility to address the S 8 tension [28].
In this work we update constraints on DM-DR models within the ETHOS framework considering a range of scenarios depending on (i) the temperature dependence of the interaction rate, (ii) the fraction of DM that is interacting, and (iii) the case of free-streaming or efficiently self-interacting DR (also known as dark fluid).The main new ingredient of our analysis are the monopole, quadrupole and hexadecapole galaxy clustering power spectra in redshift space provided by DR12 of the BOSS galaxy survey [32,33], that we combine with CMB Planck 2018 legacy data [2,34] and complementary BAO-scale measurements.We use the CLASS-PT framework [35] to take the full shape information of the BOSS galaxy clustering data into account [36,37].This implementation is based on a perturbative approach to non-linear corrections of the power spectrum of biased tracers in redshift space [38,39], built on top of the large-scale bias expansion (see [40] for a review) and an effective field theory description of redshift space clustering [41][42][43][44].We analyze the capability of the various scenarios to address the S 8 tension, comparing also to a prior obtained from KiDS-1000 [10].
The features of the DM-DR interaction that are most promising to solve the S 8 tension can be realized by a SU (N ) gauge symmetry within the dark sector, leading to a suitable temperature dependence of the interaction rate as well as DR self-interaction [25].Both of these properties are intimately related to the non-Abelian gauge boson interaction.We therefore interpret our results within the parameters of this model.The interaction rate is sensitive to Debye screening and we compute and include higher order corrections in the mapping between ETHOS and the microscopic parameters of the theory.
The structure of this work is as follows.In Sec. 2 we review the theoretical basis of both the ETHOS framework as well as the perturbative description of the galaxy power spectrum in redshift space within the effective theory approach.Next, in Sec. 3 we present the datasets and statistics used in this work and detail the Markov chain Monte Carlo (MCMC) setup.The main results are presented in Sec. 4. Finally, we discuss in Sec. 5 the mapping between ETHOS and a dark sector described by a SU (N ) gauge theory.We conclude in Sec. 6.We dedicate appendices to complementary material of IDM-DR modeling and a complete calculation of the DM-DR drag opacity for the SU (N ) model, taking the complete leading-order as well as higher-order corrections into account that have been neglected so far.

Prerequisites
In this section we review the theoretical basis of this work.We start with the ETHOS formulation of DM-DR interaction.Then we summarize the large-scale bias expansion and the perturbative description of the power spectrum in redshift space within an effective field theory approach.

The ETHOS parameterization
ETHOS considers a non-relativistic DM species coupled to a DR component, taking the coupled hierarchy of moments of the phase-space distribution function for DR into account.The evolution equations for the energy density contrast δ, velocity divergence θ and higher moments Π l for DR are given by [26] δDR + 4 3 where {α l , β l } are the (l-dependent) angular coefficients for IDM-DR and DR-DR scatterings, respectively.ϕ g and ψ g are the gravitational potentials and σ DR = Π DR,2 /2 is the shear stress.The opacity of each process is given by Γ (sometimes also named κ).For non-relativistic IDM only the density contrast and velocity divergence need to be taken into account, where H = aH is the Hubble parameter rescaled by the scale factor and c IDM is the IDM sound velocity, which is small for non-relativistic IDM and does not contribute on large scales.Non-linear corrections to Eqs. (2.4) and (2.5) are included as presented in Sec.2.2.We can expand the opacities Γ as a power-law in temperature and therefore in redshift, and using energy-momentum conservation we can relate Γ IDM−DR and Γ DR−IDM , The dimensionless functions x IDM (z) and x DR−DR (z) absorb non-trivial redshift dependences, i.e. deviations from a power-law scaling, and are for simplicity set to unity in this work.Furthermore, we use z D = 10 7 as a normalization and truncate the series in the first (dominant) term in n.
The collision terms are then encapsulated in the set of parameters {α l , β l , a n , b n }, such that distinct microscopic models can be mapped onto those coefficients (see Sec. 5 for an example).Regarding DR, we consider two limiting cases: • strongly self-interacting DR (also known as dark fluid ) that leads to an efficient damping of higher moments Π DR,l → 0 for all l ≥ 2 (formally corresponding to Γ DR−DR → ∞), such that the parameters α l , β l , b n are irrelevant, • and DR that does not interact with itself at all (Γ DR−DR → 0, i.e. b n = 0), in which case we set α l = 1 while β l is irrelevant.Since there is no self-interaction, we follow the usual convention to refer to this case as free-streaming, though DR still interacts with IDM.
The most relevant parameters are the amplitude of DM-DR interaction and its temperature dependence.We parameterize the former by and consider n = 0, 2, 4. The case n = 0 corresponds to the SU (N ) interaction discussed in Sec. 5, while n = 2 and n = 4 arise for example for an unbroken Abelian U (1) interaction or a massive vector mediator, respectively [26].In addition, there is one parameter characterizing the density or equivalently temperature of dark radiation.We adopt written in terms of the ratio between the dark radiation temperature T DR and the CMB temperature T CMB .This quantity can also be expressed in terms of extra light species as ∆N eff ≡ ρ DR /ρ 1ν (also called ∆N fluid for the dark fluid case), where ρ 1ν is the energy density of one neutrino flavour.See App.A.2 for more details.The density parameter of DR is given by where ζ = 1 (7/8) for bosons (fermions) and η DR is the DR spin/color degeneracy factor.For concreteness, we take DR to be bosonic and with η DR = 2 unless stated otherwise, but note that our results can be rescaled to other values while keeping ω DR fixed (see Sec. 5).Finally, we consider the possibility that only a fraction of the total DM interacts with DR, with the remaining fraction behaving as usual cold dark matter (CDM).Typically, two interesting limits of those parameters are discussed in the literature [28].First, the weakly interacting (WI) limit [25,45], in which the DM-DR interaction is comparable to or smaller than the Hubble rate around matter-radiation equality (Γ IDM−DR (a eq ) ≲ H(a eq )).Since the power suppression for that case is milder, this kind of model in general allows for f = 1.Second, the strongly coupled limit for which DM and DR form a tightly coupled dark plasma (DP) [28,46], in which Γ IDM−DR (a eq ) ≫ H(a eq ).The precise value of Γ IDM−DR ∝ a dark drops out in the tight coupling limit, but this case is typically only viable when demanding f ≪ 1.We will discuss those scenarios further below.
Finally, we note that an interesting feature of the case n = 0 is that the ratio Γ IDM−DR /H has only a weak time dependence and therefore the interaction can be active while a wide range of scales enters the horizon, leading to a relatively mild scale dependence of the power suppression [45] (see e.g.left panel of Fig. 3 in [28]).The amount of DR, given by the parameter ξ, mainly determines the scale where the suppression in the power spectrum sets in, while the interaction strength a dark and the fraction f determine the amount of suppression.This feature of n = 0 is favorable in view of the S 8 tension, and our results confirm this observation as we shall see below (see also Appendix A).

The large-scale bias expansion
To relate the observed galaxy density contrast δ g to the underlying matter distribution we use the large-scale bias expansion, in which we write the tracer field (in our case the galaxy overdensity) in terms of the most general set of independent operators O [47-49] (2.13) Each operator is accompanied by its respective bias parameter b O . 1 The ϵ fields take into account the stochasticity of structure formation.The set of operators, up to third-order in perturbation theory2 is given by where ϕ v is velocity potential and ϕ g the gravitational potential.G n are the n-order Galileon operators and Γ 3 is the difference between the Galileon for the density and of the velocity potentials [48,49].The coefficient of the last operator can be viewed as a nonlocal bias and in addition absorbs small-scale uncertainties within an effective field theory treatment of dark matter clustering itself [41,54].
Taking the correlations among all pairs of operators into account, we can write the galaxy-galaxy power spectrum as where P lin is the linear matter power spectrum and P 1L the one-loop contribution.To avoid cluttering the text, we have omited the z dependence of all bias parameters.For a full expression for [35].For a generalization of the largescale bias expansion to more tracers, see [55].The two last terms correspond respectively to a term proportional to ⟨∇ 2 δδ⟩ ∼ k 2 P lin and the stochastic contribution proportional to k 0 and k 2 .In redshift space, the most general expression for the galaxy-galaxy power spectrum is constructed using a multipole expansion with respect to the angle between the Fourier wavevector and the line-of-sight direction ẑ.For the full expression for the monopole, quadrupole and hexadecapole used in this work, we refer again to [35].We use the same bias parameters and structure as described in [56].We use CLASS-PT [35] to calculate each term in the power spectrum.CLASS-PT uses the FFT-Log algorithm from [57,58] to boost the integral calculations.It also takes into account the Alcock-Paczynski (AP) distortions [59] from assuming a fiducial cosmology for converting angles and redshift differences in distances transverse and along the line of sight, respectively, and implements the IR-resummation [60] from [61].

Datasets and statistics
We start this section by exposing the different datasets considered in this work.Next, we explain the statistics used to compare how the models perform compared to ΛCDM and to quantify the S 8 tension in those models.
• Full-shape (FS): Data from BOSS DR12 from different cuts of the sky (NGC and SGC) [32,33,67], divided in two z slices between 0.2 < z < 0.5 (z 1 ) and 0.5 < z < 0.75 (z 3 ) .The likelihoods3 are constructed using the covariance matrix estimated from a set of 2048 mock catalogs from 'MultiDark-Patchy' [68,69] and the window(free) function of [70,71].We use the information from the power spectrum monopole, quadrupole and hexadecapole (up to k max = 0.2h Mpc −1 ), real-space power spectrum proxy from [72,73] (up to k max = 0.4h Mpc −1 )4 , the BAO post-reconstructed spectrum from [74], similarly to what was done in [56] and based on [75].For the priors in the bias and counter-term parameters from the large-scale bias expansion of Sec.2.2 we considered the same setup as described in [56], with each bias and stochastic parameters being fit separately in each dataset.We explicitly checked that none of the bias terms reaches the prior bounds. 5kin done by [23], we always include Planck and BAO data in our analysis.We focus on comparing three main setups: • Planck + BAO + FS and The inclusion of KiDS in the analysis has not the intend to provide a full joint analysis of KiDS with other datasets, since that would require a model-specific KiDS likelihood.Instead, including the KiDS prior serves as an indication of how well the IDM-DR models can account for the low-S 8 trend seen by KiDS.Moreover, combining datasets that disagree is, of course, troublesome, e.g. for the case of Planck and KiDS for ΛCDM.As we will point out, this is not the case for KiDS and Planck + BAO + FS for interacting DM.
Regarding the cosmological parameters, we explore the following parameter space: {ω b , ω cdm , log 10 10 A s , n s , τ reio , H 0 , ξ, a dark } . (3.1) We consider as our benchmark fiducial scenario the dark fluid limit Γ DR−DR → ∞ for DR self-interaction, with f = 1 and n = 0. Variations of the IDM fraction f are discussed in Sec.4.2.We also consider in Sec.4.3 the free-streaming scenario Γ DR−DR → 0 of DR without self-interaction, and considering n = 0, 2, 4.
We have used MontePython [78,79] to run the Markov chain Monte Carlo.MontePython connects to CLASS [80] to solve the Boltzmann system of equations.As previously mentioned, the non-linear power spectra were calculated using CLASS-PT. 7We used the Gelman-Rubin criteria [82] for the convergence of the chains (R − 1 < 0.02).For a dark we sampled the chains in log coordinates, using the prior log 10 (a dark ) < 20, similarly to [83]. 8Note that an upper bound for a dark is needed, otherwise for ξ → 0, an arbitrarily large value of a dark is allowed.For ξ we set a conservative prior ξ < 1, inside the CMB bounds for N eff .
One last comment is in order.Within IDM-DR models, the interaction between both components is only relevant during the early stages of our Universe.This means that for the latest stages of evolution, IDM behaves as CDM on large scales with its interaction with DR being strongly suppressed.Therefore the perturbative computation of power spectra based on perfect fluid equations complemented with an effective stress tensor is still valid and the interactions do not have to be modeled at the latest (non-linear) stages of cosmic expansion.We refer the reader to Appendix A.1 for more details.

Statistics
In order to quantify how well a model performs in fitting the data compared to ΛCDM, we use the χ 2 difference.The χ 2 difference for a model M with respect to ΛCDM is given by where χ 2 min is the best-fit value. 9From ∆χ 2 we can define the Akaike Information Criterium where N is the number of degrees of freedom of the model.The AIC is a way to penalize models that have too many free parameters using Occam's razor criteria.
In order to quantify a tension in a dataset within a model (and therefore with the same number of free parameters), we use the 'difference in maximum a posterior' as [86] We can quantify how much a dataset is in tension with other datasets by taking the square root of Q M,data DMAP .The tension of a model M with KiDS can be quantified by Notice that the statistics used here to quantify both the tension with respect to KiDS data and the preference of a model over ΛCDM resemble the statistics used to discuss how different models address the H 0 tension with respect to SH0ES in [6].

Results
In this section we discuss results for different models within the ETHOS framework.We consider three different scenarios, already pointed out in Sec.3.1: Planck + BAO +RSD, Planck + BAO + FS and Planck + BAO + FS + KiDS.We stress again that the posteriors when including KiDS should be taken with a grain of salt: the KiDS result is implemented as a Gaussian around S 8 measured for ΛCDM and not as a full model-dependent likelihood.
Here we show its effect for comparison with other models and to figure out directions in parameter space that are preferred when the low S 8 value reported by KiDS is included.We start in Sec.4.1 discussing the fiducial case, in which all of DM interacts with DR and the latter behaves as a dark fluid.Then, we dedicate Sec.4.2 for studying scenarios with different IDM fractions and Sec.4.3 for the case of free-streaming DR.The main results are summarized in Tab. 1.

Fiducial (dark fluid,
We start by presenting the results for the fiducial scenario, in which DR is described as a dark fluid, all DM interacts with DR (f = 1), and the temperature dependence of the interaction rate is comparable to that of the Hubble rate (n = 0).
In the right panel of Fig. 1, we show the results for the posterior of S 8 for ΛCDM (dashed) and the fiducial model (solid).The difference between the red dashed line and the gray band, representing the nominal 1σ KiDS bounds, clearly shows the tension between Planck + BAO and weak-lensing data within ΛCDM.Combining Planck and BAO with FS does not change the posterior to smaller S 8 values (blue dashed line), indicating that Planck data dominates the posterior.
For the fiducial IDM model, we see that Planck + BAO + RSD (red solid) already prefer smaller values of S 8 if compared to ΛCDM.Not only does the width of the posterior increase,  which is expected from adding more free parameters, but also a skew in the direction of smaller S 8 arises (note that the same does not happen on the right side of the posterior).The inclusion of FS leads to an even more skewed tail and a small shift towards low S 8 , compatible with [56] 10 .We see that Planck+BAO+FS is not incompatible with KiDS for IDM and in that case it makes sense to combine both datasets.The S 8 tension is alleviated from 2.9σ for ΛCDM to 1.5σ for the fiducial IDM model (see Tab. 1).In Sec.4.2 we will see that models with a reduced fraction of IDM are more effective in reducing the tension.Moreover, one may wonder how degenerate are the bias parameters of the FS analysis with a shift in S 8 .We explicitly checked that the shifts in S 8 are not being trivially absorbed by b 1 or other bias parameters.
In the left panel of Fig. 1 we display the combined posteriors for the fiducial IDM cosmology.We can see, e.g., in the S 8 vs. ξ panel a clear feature: an allowed direction in parameter space that can substantially reduce S 8 .When combining with KiDS data (light green), this is the region that is mostly preferred, favoring at 1σ a DR temperature ξ = 0.117 +0.026 −0.016 Furthermore, for Planck + BAO + FS data, ∆χ 2 = −5.5 for IDM compared to ΛCDM (see Tab. 1), showing slight preference for IDM with f = 1 according to the AIC criterion.However, the inclusion of KiDS data leads to ∆χ 2 = −11.7,indicating a strong preference for IDM compared to ΛCDM.
From Fig. 1 we can see that the blue (Planck+BAO+FS) posteriors favor somewhat larger values of a dark and ξ compared to Planck + BAO + RSD (red), and become narrower for ω cdm and H 0 .This means that the FS analysis does provide relevant additional information within the IDM model.In addition to providing constraints for some of the cosmological parameters, FS data strongly shifts the values of IDM parameters away from the ΛCDM limit (larger IDM temperature ξ and stronger interaction).
The contributions of each likelihood to the total χ 2 are displayed in Tab. 2. We highlight that, when including Planck and FS data, the analysis of IDM for the fiducial case improves χ 2 relative to ΛCDM, but within the expected amount given the extra free parameters.This expected reduction is subtracted off for the ∆AIC statistics (see Tab. 1), which does not strongly discriminate between models when leaving out weak lensing data, but slightly prefers IDM.Note that when including weak lensing (KiDS) in the likelihood for ΛCDM, the value of χ 2 increases by 8.5, most of this contribution coming from KiDS itself (and also another part coming from FS), exposing the S 8 tension within ΛCDM.In the case of IDM with f = 1, χ 2 only changes by 2.2 when including KiDS, with Planck high ℓ being the main responsible for this shift.FS χ 2 improves by about 2 when including KIDS.We can deduce that IDM accommodates KiDS and FS data way better than ΛCDM, with a small deterioration coming from the Planck high ℓ fit.This indicates that an analysis including Atacama Cosmology Telescope (ACT) data and South Pole Telescope (SPT) data is an interesting direction of investigation in the future.An important additional probe of scenarios suppressing the power spectrum comes from the Lyman-α forest.In this work we focus on the BOSS FS analysis, and leave a quantitative analysis including Lyman-α data to the future.However, we note that for the scenario with n = 0 considered here, the suppression of the power spectrum is rather flat, since the DM-DR interaction rate has a similar time dependence as the Hubble rate, such that the interaction is active while a wide range of scales enter the horizon [28] (see Appendix A.1).This is also the reason why Lyman-α data are less relevant for this scenario.Indeed, in [29] it was found that Planck data are more constraining than HIRES and MIKE Lyman-α observations [88] for IDM.A similar conclusion was reached in [89] based on BOSS Lyman-α obervations [90,91] and using a conservative analysis regarding astrophysical uncertainties.The suppression of the power spectrum is also rather gradual if only a fraction of DM interacts with DR [28,46].In the next section we investigate this possibility.
Before moving on, however, one more comment with respect to H 0 is noteworthy.The H 0 posteriors for ΛCDM and IDM are shown in the right panel of Fig. 1.We see that the inclusion of FS data improves the constraints on H 0 for IDM when compared to Planck + BAO + RSD.The peak position does not shift towards higher H 0 (in the SH0ES direction [1]), i.e. the H 0 tension is not solved by IDM, in agreement with former observations [6] (though it could be mitigated when including lower bound priors on N eff [83], but with the price of making the fit worse).In our case we find H 0 = 68.11+0.36  −0.37 (Planck + BAO + FS, for IDM) .The result for H 0 is similar for different values of f or n, considered in the following sections of this work.

Mixed cold and interacting DM
The scenario in which only part of DM interacts with DR is motivated for instance by a WIMP dark sector (the PAcDM scenario of [46]) with two matter fields χ 1 and χ 2 , in which only χ 2 couples to DR [92].χ 1 , in that case, remains collisionless and dominates the DM abundance.For that scenario, it is the suppression of χ 2 clustering that leads to smaller S 8 .
-  We start by discussing the case in which only 10% of DM interacts with DR.Notice that the ratio of IDM and CDM densities for that case is of similar magnitude as that of baryons and DM.From Tab. 1, we can see that the tension in S 8 is substantially alleviated to 1.3σ.Also ∆χ 2 = −9.6 (with ∆AIC = −5.6)when including KiDS data, showing a preference of this model over ΛCDM.Fig. 2 shows the combined posteriors for f = 0.1.Compared to the case in which 100% of DM interacts with DR, the 10% case provides a good trade between reducing S 8 via DM-DR interaction and at the same time not inducing too much suppression and being in agreement with Planck and BAO data.Note that the 10% scenario allows for strong DM-DR interaction (i.e.large a dark ) if compared to the results of Fig. 1 when including KiDS.We can connect that to the dark plasma and weakly interacting scenarios pointed out in Sec.2.1 and [28]: when IDM composes 100% of DM, a relatively weak DM-DR interaction strength is slightly preferred.In that case the rate is comparable to the Hubble rate while the relevant modes enter the horizon.When only 10% of DM interacts, the interaction can be much stronger, leading to a tight coupling of IDM with DR, corresponding to the dark plasma scenario.Since our parameterization encompasses both limiting cases, the difference in the posteriors for a dark in Fig. 1 and in Fig. 2 can be clearly associated with the weakly interacting and dark plasma limits, respectively.
Moreover, note however that when only the FS likelihood is included, there is no preference for 10% IDM (∆AIC = 1.3).Tab. 2 also shows the contribution of each likelihood to the total χ 2 of the 10% model.Note that when including KiDS, the FS likelihood strongly favors 10% IDM, with some deterioration in χ 2 from Planck data.The preferred value of the DR temperature when including KiDS is ξ = 0.140 +0.026 −0.022 (Planck + BAO + FS + KiDS, f = 0.1) . (4.3) In the left panel of Fig. 3, we show the MCMC results for the case in which IDM composes 1% of DM.In this scenario, the suppression in the power spectrum is mild and not enough to shift S 8 towards the KiDS reference value.The tension, in that case, is only mildly alleviated to 2.3σ and ∆χ 2 = −4.7 when including KiDS, with no preference for this model.The right panel of Fig. 3 shows the case in which the IDM fraction f is also a free parameter.We see that in that case the S 8 tension is alleviated to 1.9σ.When including FS, ∆χ 2 = −5.2 and ∆AIC = 0.8 showing no preference 11 .The inclusion of KiDS information induces a sharp drop of ∆χ 2 = −10.2and ∆AIC = −4.2,favoring then IDM over CDM.Notice that the case in which f = 0, that resembles the case of CDM plus extra DR, is not preferred over other parameter points when KiDS is included.There is even a slight preference for non-zero values for f , with a bump at f ∼ 10%, being the case already discussed above that produces good agreement with KiDS data.

Free-streaming DR
We now move to the free-streaming scenario for DR.In that case, we assume that the selfinteraction rate of DR in (2.7) is negligible and the entire hierarchy of moments for the Boltzmann system is solved up to l = 17 in CLASS.The free-streaming case is especially relevant for example in the case of DR being sterile neutrinos or the case in which DM is charged under an Abelian U (1) symmetry.In the first scenario, the scaling of the DM-DR interaction with temperature would lead to n = 4.For the second case, it would lead to n = 2 [26].We also show the n = 0 free-streaming case.
The MCMC results for the free-streaming case with n = 0, 2, 4 are shown in Fig. 4. We see that only the n = 0 scenario is able to reach smaller values of S 8 .The cases n = 2, 4 have a strong temperature dependence, leading to a sharp transition from the case in which Γ/H ≫ 1 to the case Γ/H ≪ 1 (see App. A.1).As a consequence, structure formation is strongly suppressed on scales that enter the horizon before this transition.In contrast, the DM-DR interaction for n = 0 is almost constant during radiation domination, leading to a gradual suppression of the power spectrum.Figure 4: Posteriors for the free-streaming scenario when taking n = 0, 2, 4 (left, middle, right).
Correspondingly, we find that the parameter region compatible with low values of S 8 exists only for n = 0.In that scenario, we find the smallest tension regarding KIDS data, with only 0.8σ.However, as remarked before, the Q DMAP statistic relies on finding minima, for which MCMC is not optimal, and therefore the precise number should be taken with a grain of salt.The MCMC sampling required more time to converge according to the Gelman-Rubin criterion for this scenario, and several chains were necessary.Still, the results from Tab. 1 indicate a clear pattern.IDM lifts the tension between KiDS and other probes to the ∼ 1σ level, as long as n = 0 and f ≳ 0.1.
Interestingly, for n = 2, 4, we find that adding the FS information significantly improves the upper limit on the interaction strength a dark , even for rather small values of ξ.This can be traced back to the scale-dependence of the power suppression for n = 2, 4, and the sensitivity of FS data to smaller scales compared to Planck.Furthermore, we stress an important difference between n = 0 and n = 2, 4: in the first case, there is an allowed lower S 8 region that drives the preferred ξ to higher values when including FS, whereas in the second case FS analysis tightens the constraints on IDM parameter space.We also note that as opposed to n = 0, Lyman-α data are expected to provide relevant additional information for n = 2, 4 [83], but a combined analysis of Lyman-α with Planck and FS data is beyond the scope of this work.

Dark sector with non-Abelian gauge interaction
In this section we provide an example for a microscopic DM-DR interaction model that falls into the category of models studied above.After briefly reviewing the model, we present a computation of the DM-DR interaction rate taking Debye-screening into account, including a non-analytic correction at next-to-leading order in the coupling expansion.We then map the previously derived constraints on the model parameter space.Furthermore, we discuss implications for small-scale structure formation as well as constraints from requiring the gauge interaction to be weakly coupled until today.
The non-Abelian IDM-DR model [45] is described by an unbroken SU (N ) dark gauge symmetry and a massive Dirac fermion dark matter multiplet χ in the fundamental representation, where is the dark gauge field.The free parameters are, apart from N , the dark gauge coupling g d (or equivalently α d = g 2 d /(4π)) and the DM mass m χ .DR is described by massless dark gauge bosons with η DR = 2(N 2 − 1) degrees of freedom, while η χ = 2N for DM.Here we assume no interactions of DM or DR with the Standard Model that are relevant during horizon entry of the perturbation modes probed by CMB and LSS, and during structure formation.
The model can be complemented by an additional massive field that transforms trivially under the gauge symmetry, providing a non-interacting contribution to the total dark matter density (called χ 1 in Sec.4.2).This option is necessary for scenarios in which only a fraction f of the DM is interacting.Apart from the value of f and the property of being cold and collisionless, no microscopic information about this component is required.
The non-Abelian dark gauge boson self-interaction ensures that DR behaves as a dark fluid, described by its density contrast δ DR and velocity divergence θ DR , while all higher multipole moments are driven to zero and taken to be vanishing.Furthermore the temperature dependence of the DM-DR interaction rate maps onto the n = 0 scenario (see below for details), such that the non-Abelian model coincides with the settings considered in Sec.4.1 and Sec.4.2 for all or a part of DM being interacting, respectively.We also assume the gauge coupling to be sufficiently weak such that the confinement scale of the gauge interaction is much smaller than the DR temperature today, as discussed below.

DM-DR interaction rate
The interaction between DM and DR is analogous to Compton scattering of baryons and photons, with the difference of an additional t-channel contribution to the scattering amplitude (see Fig. 5) involving the non-Abelian three-gluon vertex.DM-DR interactions that affect the power spectrum take place in the non-relativistic limit for the DM particle, T DR ≪ m χ . 12n this limit, the non-Abelian interaction leads to an important difference compared to the Abelian (Thomson) scattering process.For small scattering angles the momentum transfer in the t-channel goes to zero, leading to a strong enhancement of the non-Abelian contributions.Indeed, formally, the cross-section and thereby the interaction rate would diverge logarithmically when integrating over all angles.This divergence is regulated by Debye-screening in the non-Abelian plasma, described by the Debye mass [93]

DM-DR Interaction in ETHOS Formalism
where C A = N is related to the gauge group and N f is the number of light fermionic degrees of freedom, that is zero for the relevant temperature range and within the minimal IDM-DR model considered here, but that we include below for generality.More precisely, the t-channel propagator needs to be replaced by the hard thermal loop (HTL) resummed propagator [93], taking the Dyson resummation of self-energy diagrams shown in Fig. 6 into account.The HTL approximation captures the leading dependence in the large temperature limit, i.e. for a momentum transfer that is much smaller than T DR , and taking only the dominant contribution from (ultra-)relativistic particles in the self-energy loop into account.The HTL approximation amounts to implementing the following replacement of the gluon propagator [94] where L and T refer to longitudinal and transverse projections of the gluon self-energy Π µν given e.g. in [95], corresponding to the electric and magnetic interactions, respectively.Up to corrections that are relatively suppressed by higher powers of T DR /m χ only the longitudinal part contributes, and the momentum transfer is dominated by its spatial part q ≃ (0, ⃗ q), corresponding to the Coulomb limit.Furthermore, the Debye-screening is relevant for |⃗ q| ≪ T DR .Altogether this means that to obtain the leading dependence on T DR /m χ , we can approximate Π L → m 2 D , while neglecting the transverse contribution.An explicit computation taking the Debye-screening into account yields for the interaction rate (see App. B) 13Γ IDM−DR = −a π 18 with α d = g 2 d /(4π).The logarithmic term arises from Debye-screening and agrees with previous results [26].Here we provide also the constant term as well as the first-order correction with coefficients c 0 = 1 + ln 6 2N + N f + ln(4π) − 24 ln (A) , and with A ≃ 1.28243 being the Glaisher-Kinkelin constant.The non-analytic correction linear in g d originates from averaging over the thermal Bose-Einstein distribution of dark gauge bosons, which is strongly enhanced for small momenta.A similar effect is well-known in the computation of various transport rates at finite temperature [93].While the logarithmic and linear terms receive contributions only from the t-channel process, also the interference of s-and u-channel with the t-channel amplitudes contribute to the constant term c 0 .On the other hand, the Thomson-like contributions (being the square of the s-and u-channel as well as their interference) are suppressed by a relative factor T 2 DR /m 2 χ .We refer to App.B for details of the computation.

Mapping of cosmological constraints on the model parameter space
The DM-DR interaction rate (5.4) can be mapped onto the ETHOS parameterization (2.8).Neglecting contributions that are suppressed in the non-relativistic limit T DR ≪ m χ leads to a n = 0 for n > 0 and a dark ≡ a 0 14 where ξ N = T DR /T γ (see below), T γ,0 is the present CMB temperature, and Ω γ h 2 its density parameter.In the second line we used natural units to convert the rate into the dimensions used in the previous sections.
The DR energy density is related to its temperature according to ρ DR = π 2 30 η DR ξ 4 N T 4 γ .Within the parameterization used for the MCMC analysis, the DR temperature enters explicitly only via ρ DR .However, for the input parameter ξ used in the previous sections we assumed a fiducial value η fid DR = 2, while for the dark SU (N ) model η DR = 2(N 2 − 1).Nevertheless, we can map both cases by requiring that they correspond to the same DR energy density, which implies that ξ N = ξ(N 2 − 1) −1/4 .Here ξ N is related to the actual temperature of the SU (N ) dark plasma and ξ is the MCMC parameter used in the sections above.Alternatively, the DR energy density can be parameterized by (see App A) (5.7) If the dark sector was in thermal equilibrium with the thermal bath (of the visible sector involving SM particles) above some decoupling temperature T dec , its temperature is determined by entropy conservation to be where s vis = π 2 45 g vis * s T 3 and s d = π 2 45 g d * s T 3 d are the entropy densities of the visible and dark sector, respectively, and the second line assumes T dec > m χ and uses that T ≪ m e , T DR ≪ m χ during the epoch when perturbation modes relevant for CMB and LSS enter the horizon.For T dec < m χ the only change is that the second factor in the expression after the second equality sign is absent.This can be compared to the required amount of DR to address the S 8 tension, of order ξ N S 8 ∼ 0.1(10ξ MCMC )(N 2 −1) −1/4 , see Sec. 4. Thus, if the dark sector was in thermal contact in the very early Universe, this would require additional heavy particle species within the thermal bath to be present such that g vis * s (T dec ) ≫ 106.75.Alternatively, the perhaps more plausible option is that the dark sector was never in thermal contact, allowing for in principle arbitrarily small values of ξ N .
The SU (N ) gauge interaction becomes stronger at low energies due to the running coupling determined by the β function where C A = N and we assumed no light degrees of freedom apart from the gauge bosons for the relevant regime T DR < m χ .The condition that the confinement scale Λ c is below some maximal value Λ max c can be expressed as an upper bound on α d evaluated at some reference scale µ ref via where we used the position of the Landau pole obtained from the one-loop beta function as a proxy of the confinement scale.In practice, we evaluate the reference scale at the typical energy scale of the IDM-DR scattering taken to be the DR temperature when the relevant modes k ≲ 1h/Mpc for CMB and LSS start to enter the horizon, Apart from the confinement constraint on α d , we also include a bound on the long-range interaction strength derived from the observed ellipticity of the gravitational potential of the galaxy NGC720 [97], where we included a factor C F /2 = (N 2 − 1)/(4N ) relative to the Abelian case obtained from the color factor for χχ → χχ scattering and averaging over the dark SU (N ) degrees of freedom of one of the incoming χ, while summing over all others, as appropriate for the scaling of the scattering rate of a given χ particle.In addition, we estimate the change of the bound if only a fraction f of DM is interacting by rescaling it according to the decrease of the scattering rate.The magnitude of the self-scattering cross-section can be estimated as [97] where v is the relative velocity.Self-interaction cross-sections of order cm 2 /g are considered in the context of small-scale puzzles, specifically the core-cusp problem [20].We use the result (5.6) for the IDM-DR rate to map the Planck+BAO+FS constraints obtained in Sec.4.1 and Sec.4.2 onto the parameter space of the dark fine-structure constant α d = g 2 d /(4π) versus the DR temperature ξ N .The results are shown in Fig. 7, assuming m χ = 1TeV and N = 3 (left panels) or m χ = 1GeV and N = 2 (right panels) for illustration.case.Furthermore, for f = 100% there is a tendency to prefer smaller values of α d for points with S 8 in the KiDS range, although large values are also possible.For f = 10%, the interaction strength is only bounded from below by α d ≳ 10 −8 (10 −9 ) for m χ = 1000(1) GeV when requiring S 8 compatible with KiDS, in accordance with the limiting case of a tightly coupled dark plasma.
In addition, the mapping allows us to compare complementary constraints from confinement and the impact of long-range interactions on elliptical galaxies, as discussed above.The corresponding exclusion regions are shown as hatched and orange shaded areas in Fig. 7, respectively.We see that ellipticity bounds are more constraining for light dark matter masses, and confinement for heavy masses in the TeV regime.Nevertheless, these constraints are compatible with points in parameter space that are allowed by cosmological Planck+BAO+FS data and have S 8 values favored by KiDS.In addition, the allowed parameter space is compatible with a self-interaction cross-section of order 1 cm 2 /g, being relevant in the context of small-scale structure puzzles.
As anticipated above, the DR temperature favored by KiDS is significantly smaller than the one that would be expected if the dark and visible sectors would have been in thermal equilibrium in the very early Universe.This points to a dark sector that was never in thermal equilibrium with the Standard Model.A possible production mechanism of dark sector particles, in that case, is via freeze-in [98].In a minimal setup, freeze-in can occur via gravitational interactions only [99], naturally leading to a dark sector temperature that is significantly smaller than the SM, depending on the reheating temperature.The further evolution of the dark sector depends on whether the dark gauge coupling is strong enough to establish thermal equilibrium [100].If that is the case, the annihilation χχ → A d A d can lead to a further depletion of the density of χ particles, and its freeze-out within the dark sector determines the DM abundance [97].A similar evolution occurs if the dark and visible sectors interact via higher-dimensional operators suppressed by some large energy scale.This interaction would leave the phenomenology of structure formation discussed here unchanged, but could set the initial DR temperature and DM abundance via UV freeze-in [101] and subsequent dark sector thermalization and freeze-out [102].In all these scenarios the DR temperature and DM abundance are related to the reheating temperature, which is in turn related to the scale of inflation and thereby to the tensor-to-scalar ratio r.In particular, gravitational production typically requires a high reheating/inflation scale, leading to a lower bound on r, potentially testable with upcoming CMB B-mode polarization measurements [103].We leave a further exploration of DM and DR production mechanisms to future work.

Conclusion
In this work we have derived constraints on dark matter interacting with dark radiation based on the full shape (FS) information of the redshift-space galaxy clustering data from BOSS-DR12, combined with BAO and Planck legacy temperature, polarization, and lensing data.In addition, we have quantified the extent to which the S 8 tension can be relaxed due to DM-DR interactions, by taking the value of S 8 measured by KiDS into account.We have considered a range of scenarios using the ETHOS parameterization, that differ in the temperature dependence of the interactions rate (n = 0, 2, 4), the assumptions on DR selfinteractions (dark fluid vs. free-streaming), and the fraction f of DM that interacts with DR (f = 100%, 10%, 1% or free).
We find that, for IDM-DR models, FS data add some information when compared to the case of using only measurements of the clustering amplitude and growth rate f σ 8 extracted from redshift space distortions.Furthermore, we find that Planck and FS data are compatible with low values of S 8 preferred by KiDS for models with approximately constant Γ/H (n = 0).We use several statistical indicators to assess the ability to reduce the S 8 tension.According to the Q DMAP statistic the tension between Planck+FS and KiDS can be reduced from 2.9σ for ΛCDM to a level of ∼ 1σ for an IDM-DR model with n = 0, and fraction f ≳ 10% of interacting DM, and a DR temperature ξ ∼ 0.1 about an order of magnitude below the CMB temperature.
The class of scenarios favored by relaxing the S 8 tension can be mapped on a microscopic model with a dark sector that interacts via an unbroken, weakly coupled SU (N ) gauge symmetry that was never in thermal equilibrium with the Standard Model.We compute the relevant interaction rate taking higher order corrections from Debye screening into account, and map the constraints on the parameter space of the SU (N ) model.We find that a solution of the S 8 tension requires a dark fine-structure constant α d ≳ 10 −8 (10 −9 ) for DM mass m χ = 1000(1) GeV.This lower bound is compatible with upper bounds from the asphericity of elliptical galaxies, and allows for self-interaction cross-sections relevant for small-scale structure puzzles.It would be interesting to work out a production mechanism and thermal history of a dark sector with properties hinted at by the S 8 measurements from Planck, BOSS and KiDS, for example along the lines of gravitational or UV-dominated freeze-in, which is left for future work.In addition, it will be interesting to consider complementary probes of large-scale structure, such as the cross-correlation of weak lensing and galaxy clustering or galaxy cluster number counts.the various ETHOS scenarios.

A.1 The scaling of the IDM-DR coupling
In order to understand when the IDM-DR coupling is relevant, we can start from (2.8) and write (assuming that only a single term in the n sum is relevant) where the superscript 0 denotes quantities today and M Pl is the Planck mass.We find then Here we assumed for simplicity that DR is bosonic (ζ = 1) and has two polarizations (η DR = 2).Substituting typical values for the constants, we obtain with where the dependence of P on ξ is subdominant since typically ξ 4 ≪ 1.
We display Γ IDM−DR /H in Fig. 8 for n = 0, 2, 4, ξ = 0.1 and a n = 1000Mpc −1 .We can see that the choices n = 2, 4 lead to a sharp decrease in the interaction rate with time (i.e. in the direction of lower redshift).In practice, this means that only modes that enter the horizon before the ratio falls below unity are affected by the interaction, leading to a sharp k-dependent suppression of the power spectrum on small scales.This is the reason why these scenarios cannot account for the S 8 tension, which requires a more shallow suppression of the power spectrum.In addition, it implies that model parameters where the suppression occurs on scales that are smaller than those probed by CMB and LSS data cannot be tested in that way, except for the usual impact of the extra DR density.In contrast, the n = 0 scenario exhibits an interaction rate that evolves in time in the same way as the Hubble rate, during radiation domination.That implies that the interaction is relevant while a wide range of scales enter the horizon during the radiation era, and correspondingly leads to a rather flat suppression of the matter power spectrum.as the ratio between the DR energy density and the energy density of one neutrino family.
In the last line we inserted ζ = 1 and η DR = 2 as example.

B Computational details of DM-DR interaction
This appendix is dedicated to presenting some of the computational details of the DM-DR interaction rate within the SU (N ) model.As already stated, the non-Abelian DM-DR interaction (γ d (p 1 )+χ(p 2 ) → γd p 3 )+χ(p 4 )) has three distinct contributions shown by Feynman diagrams in Fig. 5, where P i is the i th four-momentum and p i is the i th three-momentum.The first ingredient to compute is the squared amplitude, averaged over initial states and summed over final states.After evaluating the color numbers and the traces of gamma matrices using FORM [104], the result of each contribution is with s, u and t being the Mandelstam variables.Adding everything together and evaluating at t = 2p 2 1 (μ − 1) and s = m 2 χ + 2p 1 m χ (with μ = p1 • p3 , pi being unit vectors, and p 1 = |p 1 |), we get the total squared amplitude Note that we used here the corrected propagator of the gluon, meaning t → t − m 2 D in the denominator (see the discussion in Sec.5.1 for a definition of the Debye mass m D ).The first term in (B.7) comes from the t−channel alone, the second term is equivalent to the Thomson scattering in the Abelian case, and the last term is the sum of what is left from the s/u−channels and their interference with the t−channel.The approximation sign comes from neglecting terms proportional to m χ p 3 1 (μ − 1) in the numerator of the t−channel contribution, being the first expression on the right-hand side.They are suppressed in the non-relativistic limit p 1 ≲ T DR ≪ m χ .

B.1 Interaction rate (DM drag opacity)
The interaction rate in the ETHOS formalism is given by (for more details see Appendix A of [26]) In order to compute Γ IDM−DR , one needs first to compute the coefficients A 0 (p 1 ) and A 1 (p 1 ) using the general formula for the projection of the spin-summed squared matrix element onto the l th Legendre polynomial P l (x) with P 0 (μ) = 1 and P 1 (μ) = μ.Following the result we got for the full squared amplitude in (B.7), we split the computation of the interaction rate into two parts.We start by the t−channel term, indicated with a superscript t , (B.10a) leading to the following contribution to the interaction rate from the square of the t−channel diagram, Here we defined C = a 1 16πm 3 χ η DR 3 g 4 d for convenience.This result can be further split into four parts.The contribution from the first two terms in the integrand of (B.11) can be computed analytically.We therefore put them together and perform the integration by parts to obtain where A is known as the Glaisher-Kinkelin constant given by ln A = 1/12 − ζ ′ R (−1) where ζ R is the Riemann zeta function, such that A ≃ 1.28243.
The third term in the integrand of (B.11) cannot be done analytically.We start to deal with it by making some simplifications In the second step we used the change of variable x = p 1 /T DR and integration by part, and set f (0) DR (x) = (e x − 1) −1 in terms of x.In the last step we define F(g d ) as the full term in the integral over x.
We are interested in an expansion of F(g d ) in the coupling g d .However, when naively Taylor-expanding the integrand in powers of g 2 d , the individual terms would lead to ill-defined integrals.This means F(g d ) does not have an analytic expansion in g 2 d .Inspecting the integrand we see that for large x ≫ 1 it is exponentially suppressed due to the DR distribution function f (0) DR → e −x .However, due to the infrared singularity f (0) DR → 1/x for small x (related to Bose enhancement), the integral is sensitive to the integration region x ≲ g d ≪ 1.To make progress, we use the following method to obtain the asymptotic behavior for small g d , inspired by the method of regions [105].The idea is to split the integral into two parts First, we perform the integration between 0 and a certain cutoff λ, chosen such that g d ≪ λ ≪ 1, allowing us to expand f DR for small x.Second, we perform the integration from this cutoff λ to infinity.For this second part we are allowed to expand the integrand in powers of g 2 d , since g d ≪ x in the λ ≤ x < ∞ region, while the small-x behaviour is important for 0 ≤ x < λ.In the end, of course, F should be independent of the arbitrary parameter λ.
For the first part F λ 1 , to study the properties of the asymptotic expansion, we take the expansion of the Bose-Einstein distribution up to the fifth order into account and we perform the integration from 0 up to the cutoff λ, We could not analytically compute this integral for an arbitrary constant λ.That is why we computed it for different values of the cutoff to determine its λ-dependence behavior.Keeping all terms appearing in (B.15), we summarize the results in Tab. 3. We find that the result can be written as a Taylor series in g d instead of g 2 d , as one might have naively 3 (λ) d expected.This behavior is typical for thermal corrections in presence of infrared singularities due to Bose enhancement [93].From Tab. 3 one can see that up to the sixth order in g d all terms with an odd power are independent of the cutoff.For the even powers, only g 2 d is independent of λ.We furthermore see that every term in the expansion of the Bose-Einstein distribution (corresponding to the rows in Tab. 3) contributes at fourth and sixth order in g d in a λ-dependent way 15 .The lower subscript in the coefficients c  (B.17) Just as in the former part, it was not possible to perform the integral for an arbitrary value of λ, and we performed the integral for different fixed values to obtain the structure of the result.Since the second part is analytic, Taylor expansion in g 2 d and integration do commute in this case.Since the expansion of the integrand starts at the fourth order in g d , we therefore obtain an expansion of F λ 2 (g d ) starting at the same order, and find that they are cutoff dependent.
Their cutoff dependence has to cancel those of the corresponding terms in F λ 1 (g d ), and we see that the structure of the result allows for this to be the case (nevertheless one cannot check the cancellation at order g 4 d explicitly because for F λ 1 (g d ) there is a contribution to g 4 d terms at all orders in the expansion in x).Still, the result supports the finding that contributions up to order g 3 d in F λ 1 (g d ) have to be cutoff independent, as confirmed by our explicit result.In summary, the splitting into two regions allows us to analytically extract an expansion of the function F(g d ) up to order g 3 After all of the previous discussion, we can conclude that the t-channel squared amplitude gets corrections at linear order in g d and also higher orders that one can safely neglect.We can write this part as The fourth term in the integrand of (B.11) is also a bit tricky to deal with, and we use a similar procedure for it as above, iv = Cm  In the second step, we use again the change of variable (x = p 1 /T DR ), insert the Debye mass (5.2), and perform an integration by part.In the final step, we define the function G(g d ), which again cannot be computed analytically.To proceed, we use the same logic as before, which leads to The computation of the remaining terms in the squared amplitude in (B.7) coming from the square of the s/u-channel and interference terms (we denote them by s, u, . . . ) is straightforward A s,u,...The first term is highly suppressed since we are taking the limit where m χ ≫ T DR .The second term gives a further correction to the interaction rate, at the same order as the O(g 0 d ) terms in the expansion used in (B.

Figure 1 :
Figure 1: Left: MCMC result for the fiducial scenario with f = 1 for three different dataset combinations: Planck + BAO + RSD (red), Planck + BAO + FS (blue) and Planck + BAO + FS + KiDS (green) .The KiDS 1σ bounds on S 8 are shown as gray bands.Right: posteriors of S 8 and H 0 for the fiducial model (solid) and for ΛCDM (dashed) for the same dataset combinations.

Figure 2 :
Figure 2: Same as the left panel of Fig. 1 but for IDM with interacting fraction f = 0.1.

Figure 3 :
Figure 3: Posteriors for IDM fraction f = 0.01 (left) and for f as a free parameter (right).

1 Figure 6 :
Figure 6: Contributions to the self-energy of the dark SU (N ) gauge bosons, relevant for Debye screening.Dotted lines represent ghost propagators.

4 Figure 8 : 3 N 3 (
Figure 8: Ratio between Γ IDM−DR and H as a function of redshift z for different n.A.2 Mapping ξ onto ∆N In this section we review the relation between the distinct ways to parametrize the energy density of dark radiation.The total energy density for radiation ρ r = ρ γ + ρ ν = 1 + 7 8 4 11 4/3 1 ) = (e p 1 /T DR − 1) −1 is the unperturbed DR distribution function, n (0) χ is the spatially homogeneous DM number density, and ρ DR = η DR ζπ 2 T 4 DR /30 is the homogeneous part of the DR energy density (with ζ = 1 for bosonic DR and ζ = 7/8 for fermionic DR).

Table 3 :
Contributions to the small-x region 0 ≤ x < λ of the function F λ 1 (g d ), defined in (B.15).Rows: powers of x in the expansion of f (0) DR (x) = (e x − 1) −1 around zero.Columns: powers of g d .

( 1 ) 1 (i( 2 )
i (λ) refers to the term in the expansion of the Bose-Einstein distribution and the upper index refers to the first part of the integral, i.e. the small-x region 0 ≤ x < λ.By construction, the cutoff dependence of these g 2k+4 d terms (k = 0, 1, . . . ) has to cancel out with the cutoff dependence of corresponding terms appearing in the second part of the integral.For now, we can compact the result we obtained for the first part asF λ (λ)g 6 d + O(g 7 d ) .(B.16)For the second part, i.e. the contribution to F(g d ) from the region λ ≤ x < ∞, the expansion is performed in small values of g 2 d (since g d ≪ x in this region)(λ)g 4 d + d(2) (λ)g 6 d + O(g 8 d ) .

2
combine our results in (B.12, B.18, B.21) to get the full contribution to the IDM-DM interaction rate from the square of the t−channel amplitude, N f ) +ln(4π) − 24 ln (A) + O(g 2 d ) .(B.22)

for p 1 ≫
m D .(B.23b)By combining these two results the Debye mass drops out, and the contribution to the IDM-DM interaction rate from the square of the s/u-channel and interference terms reads Γ s,u,... IDM−DR = C 22).It arises from the interference of the t−channel diagram with the s− and u−channel contributions.We can combine it with the previous result from the square of the t−channel contribution in (B.22) to obtain the full interaction rate of the non-Abelian DM-DR model, Γ IDM−DR = −C

Table 1 :
Total χ 2 , ∆χ2, ∆AIC and Q KiDS DMAP for the best-fit of different models.We compare distinct dataset combinations (see Sec. 3.2).When considering the fluid scenario, we fix n = 0.The fiducial and free-streaming scenarios correspond to f = 1.

Table 2 :
Contribution of different likelihoods described in Sec.3.1 to the total χ 2 of the best-fit for the ΛCDM, fiducial (f = 1) and f = 0.1 models.The values of χ 2 are truncated in the first decimal, which leads to KiDS χ 2 being shown as zero and the total sum being different than the sum of the parts for some models.