Stepped partially acoustic dark matter: likelihood analysis and cosmological tensions

We generalize the recently proposed Stepped Partially Acoustic Dark Matter (SPartAcous) model by including additional massless degrees of freedom in the dark radiation sector. We fit SPartAcous and its generalization against cosmological precision data from the cosmic microwave background, baryon acoustic oscillations, large-scale structure, supernovae type Ia, and Cepheid variables. We find that SPartAcous significantly reduces the H 0 tension but does not provide any meaningful improvement of the S 8 tension, while the generalized model succeeds in addressing both tensions, and provides a better fit than ΛCDM and other dark sector models proposed to address the same tensions. In the generalized model, H 0 can be raised to 71.4 km/s/Mpc (the 95% upper limit), reducing the tension, if the fitted data does not include the direct measurement from the SH0ES collaboration, and to 73.7 km/s/Mpc (95% upper limit) if it does. A version of CLASS that has been modified to analyze this model is publicly available at https://github.com/ManuelBuenAbad/class_spartacous.


Introduction
Over the last decade, cosmology has entered a golden age in which cosmological observables have been measured to an unprecedented level of precision.While ΛCDM, the standard model of cosmology, has been able to successfully describe a broad range of measurements over this period, in recent years there has been increasing tension between some of the most precise experimental results.The greatest source of this tension arises from the differences in the various measurements of H 0 , the current expansion rate of the universe.Indirect measurements involving a fit of ΛCDM to the cosmic microwave background (CMB) [1][2][3] and large scale structure (LSS) data [4][5][6] favor lower values of H 0 ≲ 68 km/s/Mpc [1][2][3][4][5][6][7][8][9][10][11][12][13][14] than more direct methods based on the so-called cosmic ladder of standard candles such as Type Ia Supernovae (SN) [15][16][17][18][19][20][21][22][23][24][25][26], which prefer H 0 ≳ 70 km/s/Mpc (see also Ref. [27] for H 0 measurements using strongly lensed quasars which also find large values for H 0 in tension with CMB, although there are potential systematics that could be impacting these measurements [28]).Comparing the most precise measurements in each of these two categories, namely the ΛCDM fit to Planck CMB data [1] which gives H 0 = 67.36 ± 0.54 km/s/Mpc, and the supernovae measurements made by the SH0ES collaboration calibrated to Cepheid variable stars [15] which yield H 0 = 73.04 ± 1.04 km/s/Mpc, the tension has reached 5σ significance [29][30][31].Another long standing source of tension, albeit more modest in significance than that of H 0 , involves the amplitude of the matter power spectrum at relatively small scales, conventionally expressed in terms of the S 8 parameter.Direct measurements of this parameter, which is defined in terms of the matter energy density fraction Ω m and the variance σ 2  8 of matter overdensities at 8 Mpc/h as S 8 ≡ σ 8 Ω m /0, 3, have also consistently been in 2 − 3 σ tension with the value inferred from the ΛCDM fit to the CMB. 1 These tensions could be the first evidence of the need to depart from the ΛCDM paradigm, and have motivated a considerable effort in the study of extensions of the standard model of cosmology that can accommodate them.For a sample of proposals that aim to solve the H 0 tension, see Refs. .A more comprehensive list may be found in Refs.[30,31,56,57] and references therein.Proposals to solve the S 8 tension include, for example, Refs.[58][59][60][61][62][63][64][65][66][67][68][69][70][71][72][73].Unfortunately, many of the most promising proposals to address the H 0 tension lead to an increase in S 8 , making the latter tension more significant.It is therefore important to consider models that can simultaneously address both tensions.Efforts in this direction include Refs.[74][75][76][77][78][79][80][81].
In a recent publication we put forward a new joint solution to the H 0 and S 8 tensions, the Stepped Partially Acoustic Dark Matter model, ("SPartAcous") [77].This scenario naturally combines two mechanisms that have been proposed to alleviate the cosmological tensions in interacting dark sector models, namely a dark radiation (DR) bath with a massthreshold [49] that leads to a step-like feature in the fractional energy density in radiation, and a subcomponent of dark matter which is kinetically coupled to this DR [63,64].The change in the energy density in radiation can address the Hubble tension by decreasing the sound horizon [82][83][84].At the same time the interactions between dark matter and DR give rise to Dark Acoustic Oscillations (DAOs) that suppress structure at small scales, and can thereby help resolve the S 8 tension [77].The presence of the mass threshold distinguishes SPartAcous from Partially Acoustic Dark Matter ("PAcDM") [63,64], and leads to very different predictions for the matter power spectrum.The reason for the difference is that the interacting dark matter subcomponent decouples from the DR once the temperature falls below the mass threshold, so that while there is little deviation from ΛCDM on larger length scales, the amplitude at the scales relevant to S 8 is reduced.This latter feature distinguishes SPartAcous from most other models that address the Hubble tension by increasing the energy density in radiation, since they typically enhance the matter power spectrum at scales relevant for S 8 , worsening the S 8 tension.
In this work we quantitatively investigate how well SPartAcous fits a wide range of cosmological data compared to ΛCDM, and to what extent it can solve the H 0 and S 8 tensions.We also propose a simple generalization of the model, labelled SPartAcous+, in which we enlarge the number of massless states in the dark sector, thereby altering the size of the step.As with SPartAcous, we evaluate how well SPartAcous+ fits the data and addresses the tensions.We have implemented these models in a modified version of the Boltzmann code CLASS, which we make publicly available at github.com/ManuelBuenAbad/classspartacous.Our implementation includes a generalization of the tight-coupling approximation [85][86][87] for a fluid that has a mass threshold, which speeds up the code.This is described in detail in Appendix B. Our implementation also includes a correction to the superhorizon initial conditions of the adiabatic perturbations due to the time-dependence of the equation of state, which was missed in earlier work on dark sector models with a mass threshold.
This paper is organized as follows.In Sec. 2 we review the original SPartAcous model and present its generalization, SPartAcous+.Sec. 3 sets forth the implementation of the SPartAcous and SPartAcous+ models in CLASS, lists the data we used to fit the models, and describes the results.Our conclusions are summarized in Sec. 4. Details about our CLASS code are expounded upon in Appendix A (equations and initial conditions), and Appendix B (approximations), while further numerical results and tables are included in Appendix C.

The Model
In this section, we first quickly review the SPartAcous model proposed in Ref. [77] and then present a simple generalization, SPartAcous+.The original SPartAcous model adds three additional fields to ΛCDM: a new massless (Abelian) gauge field A, and two new fields charged under this gauge group, a light vector-like fermion ψ and a heavy scalar χ.The heavy scalar χ constitutes a subcomponent of dark matter, which we will label interacting dark matter (iDM) due to its interactions with DR.The primary component of dark matter is assumed to be standard collisionless cold dark matter (CDM) (see Ref. [63] for a way to implement both CDM and iDM components within a single theoretical framework).We take the fermion ψ to be sufficiently light that it behaves as DR at early times and only becomes non-relativistic during the CMB epoch.Once the dark sector temperature drops below the mass of ψ, it annihilates away into gauge bosons.Due to its entropy being transferred to the remaining radiation, we obtain a step-like increase in the energy density in DR at that time.As is conventional, we parametrize the energy density in radiation in units of N eff , the effective number of neutrinos.The part of the Lagrangian describing the model is given by, where V µν is the field strength associated with the new gauge boson A, and is the covariant derivative (with g d being the associated gauge coupling).The dark sector is decoupled from the Standard Model plasma (at the times of interest), and therefore has its own temperature T d .This temperature is directly related to the contribution of the DR to ∆N eff at late times, which we denote by ∆N IR eff .The scales we are most concerned with are the ones that enter the horizon not much earlier than matter-radiation equality, at which time χ was already non-relativistic with a cosmic abundance set at a much earlier time.We will express the abundance of χ in terms of its fractional contribution to the total energy density in dark matter, where ρ CDM corresponds to the energy density in the cold (non-interacting) dark matter component (CDM).At temperatures above m ψ , the DR and iDM components behave as a single fluid due to the tight coupling between ψ and χ mediated by the gauge interactions, shown in the first Feynman diagram of Fig. 1, and due to Compton scattering between ψ and A, shown in the second diagram.We will take m χ to be sufficiently heavy that Compton scattering between A and χ is not efficient at the temperatures relevant for the CMB.This ensures that once the dark sector's temperature falls below m ψ , the dynamics changes in two significant ways.Firstly, the energy density in DR increases due to the annihilation of ψ, creating a step in N eff , in a manner analogous to that of the Wess Zumino Dark Radiation ("WZDR") model [49].Secondly, the interactions between DR and iDM decouple due to the exponential suppression in the number density of ψ so that, from this point onward, the iDM becomes collisionless and evolves identically to CDM.This behavior stands in sharp contrast to that of WZDR+, a generalization of WZDR first proposed in Ref. [76], in which all of the dark matter interacts very weakly with the DR, eventually decoupling from it in a much slower manner.The redshift z t at which T d = m ψ characterizes the time of the transition, and it will be used in our parameter scans in place of m ψ .The dark gauge coupling g d is, in principle, another important parameter of the model.However, as we will discuss in Sec. 3, because the couplings considered are sufficiently large as to make the DR and iDM behave as a single tightly-coupled fluid, the precise value of the coupling is largely irrelevant and it only plays a role in determining exactly when ψ freezes out and χ decouples from the DR.Since both of those changes occur due to the number density of ψ becoming exponentially suppressed, the dependence of the decoupling temperature on g d is only logarithmic, leading to minimal dependence of the physics on the value of g d .Even after ψ has frozen out, the DR is prevented from free streaming by a self-interaction among the gauge bosons of the Euler-Heisenberg form, which is generated at loop level when ψ is integrated out [77].

SPartAcous+
A simple way to generalize the SPartAcous model is to allow for different values in the size of the step in ∆N eff at the mass threshold, while retaining the decoupling of DR and iDM.The step size, i.e. the fractional change in ∆N eff , is related to the fractional change in the number of relativistic degrees of freedom, r g , as where the superscript UV (IR) corresponds to quantities above (below) the mass threshold, and g * represents the effective number of relativistic degrees of freedom in the dark sector interacting fluid.For the original SPartAcous model, g * UV = 11/2 and g * IR = 2, and so ∆N eff increases by approximately 40% at the step.As will be shown in the next section, with such a large step, the data does not allow for a non-negligible fraction of the dark matter energy density to be iDM, and the best fit approaches a WZDR limit in which there is effectively no iDM.This motivates us to consider a generalization of our original model with a reduced step size.A decrease in the size of the step can be realized by increasing the number of massless particles in the interacting dark fluid.
In generalizing the SPartAcous model, we will introduce an additional Abelian gauge field A ′ and N df new vector-like massless fermions ξ i , charged under A ′ , which together constitute an extra DR component.ψ is also charged under this new U (1) gauge symmetry.On the other hand, χ only couples to A, and the ξ i only couple to A ′ . 2 This ensures that after the mass threshold, the interacting dark matter component decouples from the dark radiation.
The new terms added to the Lagrangian are therefore: where V ′ µν is the field strength associated with the additional gauge boson A ′ , whose gauge coupling is g ′ d .In addition, the covariant derivative acting on ψ in Eq. (2.1) is modified to include A ′ .We only require that α ′ d be large enough for the new massless particles ξ to be in equilibrium with A ′ at the time when the smallest length scales of interest enter the horizon.This condition is easily satisfied over a broad range for α ′ d [77]: where we have used the fact that T d is not significantly different from T .
Together the gauge interactions mediated by A and A ′ ensure that the entire dark sector behaves like a tightly coupled fluid down until the mass threshold.This alters the values of g * UV/IR in Eq. (2.2).Because we are adding vector-like fermions, g * UV/IR change by multiples of 7/2, and therefore the fractional change in the number of degrees of freedom as a function of N df is given by At the same time, because the new fermions ξ i are not charged under A, they do not scatter efficiently with iDM.This ensures that the model retains the important feature that the interactions between iDM and DR decouple shortly below the mass threshold.The presence of ψ, a state charged under both gauge groups, will nevertheless lead to a kinetic mixing between A and A ′ of size ϵ ∼ g d g ′ d /(16π 2 ), thereby inducing a coupling between the new massless fermions and the iDM at loop level.However, in the range of α d and α ′ d of interest, this interaction is too small to keep the iDM in equilibrium with the DR, and can be safely neglected.To be more quantitative, comparing the interaction rate Qξχ to H leads to the criterion (2.6) Eqs. (2.4) and (2.6) show that as long as 10 −12 < α ′ d < 10 −4 (for these benchmark values of α d and m χ ), the DR subcomponent consisting of A ′ and ξ i will remain self-interacting but decouple from the iDM. 3

Results
We have modified the CMB code CLASS [87][88][89][90] to include the SPartAcous and SPartAcous+ models described in the last section.We have implemented a number of approximations, analogous to the ones described in [87], in order to sufficiently speed up the code to allow an efficient exploration of the parameter space.We describe these approximations in Appendix B. Our code is publicly available at github.com/ManuelBuenAbad/classspartacous.We employ this modified version of CLASS, combined with the MCMC sampler MontePython [91,92] to investigate how well the model fits different combinations of data sets and to find the allowed parameter regions for the associated cosmological parameters.We use the Metropolis-Hastings algorithm in MontePython, and the resulting MCMC chains are considered to have converged if the Gelman-Rubin (GR) criterion R < 1.01 is satisfied, where R is the GR statistic [93].
In addition to the standard ΛCDM parameters {ω b , ω dm , θ s , ln 10 10 A s , n s , τ reio }, we include 3 new parameters that capture the important effects of the new models: the DR contribution to ∆N eff at late times, ∆N IR eff ; the fraction of the dark matter energy density in the iDM component, f χ = ρ χ /ρ dm ; and the redshift of the transition, z t , defined as the value of the redshift at which T d (z t ) = m ψ . 4Of course, within ΛCDM, these new parameters are simply ∆N IR eff = 0 and f χ = 0, with ω dm = ω CDM .For both SPartAcous and SPartAcous+ we use the flat priors ∆N IR eff ≥ 0, 0 ≤ f χ ≤ 1, and 4.0 ≤ log 10 z t ≤ 5.5 5 .We also use flat priors on the remaining ΛCDM parameters.In all of our scans we include three active neutrinos, one with a mass of 0.06 eV and the other two massless.
In principle, there are additional parameters, such as m χ and α d ; however, as discussed in [77], these only enter the relevant equations through the combination α 2 d /m χ , and thus we can keep m χ = 1 TeV fixed without any loss of generality.We will be interested in the regime in which the iDM-DR system begins its life as a single, strongly-interacting fluid.The precise redshift at which the ψ-χ scattering freezes-out only has a mild logarithmic dependence on α d , and therefore the precise value of the coupling g d is unimportant (we neglect small corrections to the tight coupling approximation).We adopt the benchmark value α d = 10 −3 in order to In this case, the dark gluons would play the role of the extra DR, coupled to ψ via the gauge interaction.After the temperature drops below m ψ , the dark gluons would still behave as an interacting fluid due to their self-interactions.We will not explore this possibility further in this paper; however, for the same numerical value of the step size, which is a function of rg = 7/(4N dc ) in this case, the cosmological history would be identical to the SPartAcous+ model. 4Throughout this paper we make the assumption that the DR is not present (or has a low enough temperature such that its energy density can be neglected) during BBN (Tγ ∼ 1 MeV), and that it is populated (or heated up) after this point.This assumption allows for larger values of ∆N eff to fix the H0 problem, without at the same time increasing the primordial Helium fraction YHe beyond what is permitted by BBN observations.This scenario can be easily realized in models where the DR has interactions with other energy density components, whose energy density gets transferred to the DR after the BBN era.In our code implementation of the SPartAcous model, we simply set ∆N UV eff = 0 in the routine of the thermodynamics.cmodule that computes YHe. 5 This relatively narrow prior on log zt allows our scans to ignore the parameter space where the mass threshold takes place either very early (in which case the model becomes indistinguishable from that of selfinteracting DR [94]) or very late (where the model becomes identical to PAcDM, [63,64,95]).
increase scanning speed and convergence.As with g d , the precise value of the new coupling g ′ d is also unimportant as long as it is within the range shown in Eq. (2.4).
For the generalized model SPartAcous+, there is a discrete choice for the extra number of fermion flavors we add to the dark fluid, N df .We focus on the scenarios with N df ≤ 3 since the impact of adding more flavors becomes smaller as N df increases, see Eq. (2.5).While this means that we could, in principle, scan over an additional parameter r g which would correspond to SPartAcous and SPartAcous+ for special discrete values of r g , this would lead to most of the points being scanned over representing non-physical realizations of the model.We therefore choose to treat each realization of the model separately.In the next section, we present our numerical results only for N df = 3, which provides the best fit to the data. 6We will make this choice explicit by referring to the model with N df = 3 as SPartAcous+3.

Experiments and Methodology
We perform a full likelihood analysis of the SPartAcous and SPartAcous+ models using various precision cosmological datasets.These include measurements of CMB anisotropies as well as galaxy, supernovae, and weak lensing surveys.We arrange these datasets into three distinct categories, identified with the abbreviations D, H and S, based on their mutual compatibility within the context of ΛCDM: • D: our baseline dataset.This includes the following experiments: -Planck: measurements of TT, TE, and EE CMB anisotropies and lensing from Planck 2018 [1] (the likelihoods dubbed 'Planck highl TTTEEE', 'Planck lowl EE', 'Planck lowl TT', and 'Planck lensing' in MontePython).
-BAO: baryon acoustic oscillations (BAO) datasets in the form of measurements of D V /r drag by the Six-degree Field Galaxy Survey (6dFGS) at z = 0.106 [96] and by the Sloan Digital Sky Survey (SDSS) from the MGS galaxy sample at z = 0.15 [97] ('bao smallz 2014'), as well as measurements of D M (z)/r s,drag and H(z)r s,drag at z = 0.38, 0.51, 0.61 from the Baryon Oscillation Spectroscopic Survey (BOSS) DR12 [98] ('bao boss dr12').
• H: the dataset containing late-universe measurements of the Hubble parameter H 0 .While there are several of these measurements [16-23, 25, 27, 28, 56], at different levels of tension with results from ΛCDM fits to D, we include the most precise: -SH0ES: the latest measurements of H 0 using Hubble Space Telescope (HST) observations of Cepheid variables in galaxies hosting 42 SNe Ia [15].The H 0 value quoted by the SH0ES collaboration is ultimately obtained from their derivation of the absolute magnitude M B of SNe Ia.We use their result, M B = −19.253±0.027, to build a Gaussian likelihood.
• S: the dataset with late-universe weak lensing and galaxy clustering measurements of S 8 ≡ σ 8 Ω m /0.3, which are in mild tension with the results from ΛCDM fits to D. We include in this group the following experiments, with which we construct two-sided Gaussian likelihoods: -DES: the Year 3 results of the Dark Energy Survey (DES): S 8 = 0.775 +0.026 −0.024 [32].-KiDS-1000: results from the Kilo-Degree Survey (KiDS-1000): S 8 = 0.766 +0.020 −0.014 [100].
We denote any combination of these datasets by the corresponding abbreviation, such as DH or DHS, following the notation of Ref. [76].
One could in principle include Lyman-α forest experiments, which are sensitive to even smaller scales and could therefore be very useful in distinguishing between various models that suppress the matter power spectrum.This data consists of the Lyman-α absorption lines present in the spectra of distant quasars (quasi-stellar objects, or QSOs), and arise from intergalactic neutral hydrogen lying along the quasars' line of sight.Indeed, since hydrogen traces the matter distribution (for redshifts of 2 ≲ z ≲ 5 and length scales of [101]), Lyman-α data can be used to probe the matter power spectrum.Some datasets commonly used in the literature include QSO spectra samples from SDSS-II [102], the HIRES/MIKE and XQ-100 experiments [103,104], and from the BOSS and the Extended BOSS (eBOSS) collaborations [105][106][107].However, at present the usefulness of these experiments to constrain models beyond ΛCDM that have an impact on the matter power spectrum is somewhat limited.In particular, some of these measurements rely on modeling the matter power spectrum in the context of ΛCDM [102,108], which may not be valid when considering models beyond ΛCDM.Others require the modeling of various astrophysical nuisance parameters in conjunction with hydrodynamical simulations of structure formation [109][110][111][112][113], which requires a prohibitively large amount of computer resources.Significant efforts to make the Lyman-α data friendlier to analysis of non-ΛCDM models have recently been undertaken.In particular, the authors of Refs.[113][114][115] introduced a description of the linear power spectrum in terms of distinct phenomenological "shape parameters" α, β, and γ (labelled the "{α, β, γ}-parametrization" in the literature), and used this parametrization, along with a suite of dedicated N-body simulations, to perform a MCMC analysis of the HIRES/MIKE and XQ-100 experiments and so derive a likelihood based on these datasets.This likelihood can then be easily used within MontePython to analyze any non-ΛCDM model whose linear matter power spectrum can be mapped into the three shape parameters employed in the analysis [67].Despite its great versatility, this parametrization cannot be mapped to the power spectra from our SPartAcous and SPartAcous+ models, or those from the WZDR+ model of Ref. [76], which means that we cannot use this likelihood.Indeed, the {α, β, γ}-parametrization can only be used when the ratio of the linear matter power spectrum of a given model to that of ΛCDM is suppressed down to zero for sufficiently large wavenumber k (small length scales) in a power-law fashion.Within SPartAcous (and SPartAcous+), the fact that f χ < 1 means that this ratio does not vanish completely, whereas in WZDR+ the suppression scales like ∼ ln(k).While efforts to make the Lyman-α datasets applicable to a wider range of models are ongoing (see for example Ref. [116]), no likelihood that can be easily applied to our model has been made publicly available at the time of the writing of this paper.Therefore, throughout the rest of this work, we limit ourselves to the cosmological observables described in the previous paragraphs, with the hope that in the future Lyman-α forest data can be used to improve our analysis.

Numerical Results
In Figure 2, we show the marginalized posterior distributions for H 0 and S 8 for the three models under consideration, ΛCDM, SPartAcous and SPartAcous+3.The left panel shows the results when fitting to the D dataset, while the right panel shows the results when fitting to the DHS data combination.We can immediately see that, even for the fit to D, both SPartAcous and SPartAcous+3 allow for larger values of H 0 compared to ΛCDM.
In particular, for SPartAcous+3, we see that, even though it does not completely solve the tension, the 95% credible region (CR) overlaps with the 1σ region from the most recent SH0ES measurement of H 0 [15].The allowed ranges for S 8 are approximately the same for all three models when fitting to the D dataset.The posteriors for the fit to DHS, on the right side of the figure, show bigger differences between the 3 models.For ΛCDM, we see that the inclusion of the extra data lowers the value of S 8 , while it leads to only mild changes to H 0 , despite the fact that H 0 is a more significant tension.This illustrates that there is significant tension between the different data sets in ΛCDM.For SPartAcous, we see that with the inclusion of the SH0ES prior, it can accommodate much larger H 0 values, but it is not able to simultaneously lower S 8 .In fact, we find that when the extra data is included, the allowed range for f χ , which controls the decrease in the power spectrum at small scales, becomes very small and peaked at zero iDM, as can be seen in Figure 3.As we will see later, this is driven by the inclusion of H. On the other hand, we see that with SPartAcous+3 we can achieve even larger values of H 0 , while simultaneously lowering S 8 .This was the expectation for this class of interacting DR models, in which increasing ∆N eff allows higher H 0 while increasing f χ reduces S 8 .In Table 1, we show the best-fit points for the parameters of the three models when fit to D and DHS, as well as their corresponding χ 2 , broken down by the contributions from each dataset.Note that when fit to D, both SPartAcous and SPartAcous+3 improve the χ 2 , although only by a limited amount considering that these models have 3 extra parameters compared to ΛCDM.However, once we include H and S, we see that both models lead to a significant improvement over ΛCDM, with SPartAcous+3 giving ∆χ 2 = −24.6.Even when taking into account the penalty for having three extra parameters (∆N IR eff , f χ , and log 10 (z t )) by using the Akaike Information Criterion (AIC), the improvement remains substantial.Indeed, from AIC ≡ χ 2 b.f.+ 2 × d.o.f., we find ∆AIC = −18.6 (−11.2) for fits of SPartAcous+3 (SPartAcous) to DHS; see Table 3.By contrast, fits to only D are disfavored by the AIC as can be seen in Table 2.In fact, note that once the H and S data are included in the fit, SPartAcous and SPartAcous+3 provide an improved χ 2 compared to ΛCDM even when computing only the contribution coming from CMB data.The results in Table 1 also show that SPartAcous+3 provides an overall better fit than SPartAcous, and as expected from Figure 2, we see that SPartAcous+3 provides a more significant reduction of both the H 0 and S 8 tensions.As expected, the inclusion of the H prior pushes H 0 to higher values, in both ΛCDM and the SPartAcous and SPartAcous+3 extensions.Some important points to note are that: (i.) SPartAcous and SPartAcous+3 allow for larger values for H 0 than ΛCDM with or without the H prior, (ii.) the overall global fit is better in SPartAcous and SPartAcous+3 than in ΛCDM, even accounting for the extra degrees of freedom, and (iii.) the improvement in χ 2 SH0ES in the SPartAcous and SPartAcous+ models does not come at the cost of a degradation of the χ 2 CMB , which is lower in both models compared to ΛCDM when fitting to DHS.

Model
We present the 95% CR range for the parameters ∆N IR eff , f χ , H 0 , and S 8 in Table 2 for the fit to D and in Table 3 for the fit to DHS.The posteriors for the new parameters, ∆N IR eff , f χ , and log z t , are also shown in Figure 3 (see Appendix C for the triangle plots with the posteriors of all parameters).One can see from Table 2 that both classes of models allow for much larger values of ∆N eff , even without including the SH0ES results, compared to simple extensions of ΛCDM where DR is purely free-streaming [1] or purely interacting [94].These larger values for ∆N eff lead to larger values for H 0 ; in particular, SPartAcous+3 yields a value of H 0 as large as 71.4 km/s/Mpc at 95% CR, compared to at most 68.6 km/s/Mpc for ΛCDM.This table also shows that, in SPartAcous, the maximum allowed f χ is much smaller than in SPartAcous+3, which explains why SPartAcous+3 is more effective in addressing the S 8 tension.In Tables 4 and 5, we show the mean and 1σ allowed ranges for all the parameters in the three models, fit to D and DHS respectively.We can see that when fitting to D, the main difference in the parameters common to all models is that in both SPartAcous and SPartAcous+, the allowed ranges for ω dm and n s are widened and the central value pushed to slightly larger values to compensate for the inclusion of extra dark radiation, as has been the case in other models with additional radiation [49,94].This is expected in order to keep the redshift at matter radiation equality fixed, and also to compensate for the change in Silk damping.In the fit to DHS, due to the pressure from the SH0ES data to increase H 0 , the mean for ∆N eff is significantly increased.This also pushes the mean values of ω dm and n s to larger values.Table 2.A summary of the fits to the dataset D for the ΛCDM, SPartAcous, and SPartAcous+3 models, showing the allowed parameter range at 95% CR for ∆N IR eff , f χ and H 0 .These intervals are defined to be the narrowest interval containing 95% of the integrated posterior density, and have been computed directly from the posterior densities and not using Gaussian fits to the posteriors.cous+3, we see that despite the significant increase in H 0 values, the allowed range for S 8 remains almost unchanged compared to the fit to D. This is in contrast to many of the most prmoising proposals to address the H 0 tension (see e.g.[31]).

Model
In order to quantify whether the models are compatible with the inclusion of both datasets, we will follow the approach used in [7] (see [117] for more in depth discussion) and compute the Q DMAP of a given model fit to the two datasets, where This quantity is generally interpreted as the number of standard deviations (i.e.σ) as a measure of the tension between the two datasets.For the datasets to be compatible within the context of a given model this quantity is expected to be small; a large value indicates that, for the model under study, there is tension between these datasets.In Table 6, we show the χ 2 min for the fits to DH and the corresponding Q DMAP values for the three models under consideration.We see that there is significant improvement in the tension in both SPartAcous and SPartAcous+3 compared to ΛCDM, going from ∼ 5σ tension in ΛCDM to ∼ 2σ in SPartAcous+3.The improvement achieved with SPartAcous+3 is similar to the one recently found in Ref. [76] using WZDR+, although their definition of D does not include the lensing likelihood, and so a direct comparison of the results is not possible.It is also similar to the results obtained by the models ranked highest in the H 0 Olympics [31], showing that SPartAcous and SPartAcous+3 are amongst the most successful proposals to solve the H 0 tension.Table 6.List of the χ 2 values for the best-fit points of the ΛCDM, SPartAcous, and SPartAcous+3 models when fitting to DH, and their corresponding Q DMAP when comparing to the fit to D.

Matter Power Spectrum
Given the significant improvement in the fit with SPartAcous+3, an important question is whether future measurements can distinguish it from other models that address the H 0 and S 8 tensions.Given the large number of solutions that have been proposed, a comprehensive analysis is beyond the scope of this work.We will therefore limit the discussion to a comparison with two other early universe (pre-recombination) solutions of the H 0 tension, Early Dark Energy (EDE) [37,71], WZDR, and its generalizations [49,76].For EDE and for the original WZDR [49] proposal, increasing H 0 correlates with an increase in S 8 , and therefore future measurements of S 8 with increasing precision should clearly distinguish between the two models.In fact, due to the present S 8 tension, these models are already in some tension with more direct measurements of the matter power spectrum.The recently proposed WZDR+ model, which includes very weak interactions of all of dark matter with the DR, was considered in Ref. [76] (see also Ref. [75] for a similar construction, but with an overall worse fit to DHS).As in SPartAcous+3, the interactions between dark matter and DR can suppress the power spectrum at small scales and address the S 8 tension.Nonetheless, even if both SPartAcous+3 and WZDR+ models lead to similar results for S 8 , that is but a single number; their overall impact on the full power spectrum is very different.
In Figure 5, we show the matter power spectrum for our best-fit models.As discussed in Ref. [77], at small length scales (large wavenumber k), the suppression of the power spectrum compared to the model with no iDM (f χ = 0) becomes constant.When compared to the ΛCDM fit to D, as in Figure 5, the large k limit is not quite a constant due to the difference in n s , becoming a mild power law, but on scales ranging between 0.1-1 h/Mpc, it is suppressed by a few percent.One can also see the dark acoustic oscillations that arise due to the iDM-DR interactions [77] in Figure 5.These oscillations cease once ψ exits the bath at its mass threshold, but get imprinted on scales that entered the horizon prior to that (k ∼ 0.1 h/Mpc).In EDE models there is an increase in power at smaller scales (exacerbating the S 8 problem), mostly due to a significant change in n s .On the other hand, in WZDR+, as expected from similar models in which dark matter interacts with DR [59,64], the suppression becomes larger at smaller scales, decreasing with an approximately logarithmic dependence on k.Therefore, precise measurements of the shape of the power spectrum at scales of order the S 8 scale, k ∼ 10 −1 h/Mpc, may be able to distinguish between the SPartAcous+3 and WZDR+ classes of models, as long as they can handle uncertainties due to non-linear effects.We also see from Fig. 5 that in order to differentiate the best fit of SPartAcous+3 from that of ΛCDM, we would need to understand the power spectrum at these small scales at the percent level.However, since in EDE and WZDR+ the departure from ΛCDM continues to grow at smaller scales, it might be possible to constrain such models using less precise probes of the power spectrum at even smaller scales, such as Lyman-α forests (see e.g.[118] for a study with EDE).In addition, all of these models affect the CMB at small scales (large l) in different ways, and therefore future CMB experiments, such as Simons Observatory [119], CMB-S4 [120], and CMB-HD [121][122][123][124] can play an important role in distinguishing between them.

Conclusions
In this work, we have performed an extensive study of a new class of interacting dark sector models, the Stepped Partially Acoustic Dark Matter (SPartAcous) and its generalizations (SPartAcous+), by fitting them to a wide range of cosmological data.In particular, we have explored how these models compare against ΛCDM when attempting to solve the tensions in measurements of H 0 and S 8 .We compared fits of these models to a combination of datasets that consisted of both a baseline set of experiments and the direct measurements of H 0 and S 8 in tension with it, with fits to data that included only this baseline.We found that both SPartAcous and SPartAcous+3 are very effective in reducing the H 0 tension, particularly in the case where data fitted includes the direct measurements as well as the baseline dataset.
The best improvements were found with SPartAcous+3, for which the best fit point to the baseline model, i.e. excluding the direct H 0 measurement from SH0ES, already allows for H 0 as large as 71.4 km/s/Mpc at 95% CR, much larger than what is allowed within ΛCDM.
When we include all the data in the fit, the χ 2 improvement of SPartAcous+3 over ΛCDM is 24.6, which is a very significant improvement, even when taking into account the three extra parameters in the model.Ideally, the H 0 tension would be fully resolved in the fit to just the baseline dataset.Although this is not the case, we find it promising that even without the SH0ES prior, the tension is already reduced, and that the inclusion of the SH0ES prior, which leads to severe tension with ΛCDM fits to the CMB, now allows a sizeable increase in H 0 without degrading the global fits of the SPartAcous and SPartAcous+3 models.
Our results also showed that the original SPartAcous proposal was less effective at simultaneously addressing the H 0 and S 8 tensions than expected.In this model, once ∆N IR eff is increased to accommodate larger H 0 , the fraction of interacting dark matter that is allowed by the data becomes very small, and thus it cannot sufficiently lower S 8 .In SPartAcous+3, where the step change in ∆N eff is smaller, the allowed range for iDM is larger which leads to a more significant reduction in S 8 .The predicted matter power spectrum for this model differs from other proposed solutions of both tensions at small scales, and in combination with future CMB measurements, may allow one to distinguish this model from other proposals.
Note Added: During the final stages of this work, Ref. [125] appeared, in which the authors study the fits to cosmological data of a phenomenological parametrization of dark sector models.One of the class of models they explore is directly related to SPartAcous and SPartA-cous+.Their analysis differs from ours in the data that is being used for the fits, the number of free parameters describing the models, and the choice of priors used for the MCMC.In addition, shortly after v1 of this paper was posted on the arXiv, Ref. [126] appeared.This paper considers an iDM-DR model with a mass threshold and with varying DR step size, both in its weak (à la WZDR+) and strongly coupled limits (à la SPartAcous, as in this present work).In addition, the authors include ACT, SPT, and BOSS full shape datasets.Where their work overlaps with ours, the results are in excellent agreement.

A.2 Initial Conditions
In our work we consider adiabatic perturbations of the various components of the universe in the absence of curvature.These perturbations have initial superhorizon conditions that, in principle, depend on the equation of state of the fluids to which they belong.These initial conditions have been studied in detail in the past (see, for example, Refs.[86,127]).While these are straightforward to derive in the simplest case of a constant equation of state (such as for matter-like or radiation-like fluids), they are significantly more complicated in more general cases.
The DR in SPartAcous and SPartAcous+ is one such example.It undergoes an entropy dump and a step around the redshift z t , which means that w dr and c 2 dr,s deviate from 1/3 around this time.One could in principle simply set the initial conditions at a redshift sufficiently before z t , when w dr = c 2 dr,s = 1/3 still.However, the Hubble timescale associated with such an early time can be so small that the CLASS code can encounter memory problems dealing with the correspondingly large number of time steps.Because of this we instead solve for the superhorizon initial conditions of the adiabatic perturbations at any time by expanding around 1/3 to first order in the deviations δw dr ≡ w dr − 1/3 and δc 2 dr,s ≡ c 2 dr,s − 1/3.These deviations depend on r g , and can reach a change relative to 1/3 of ≲ −20% for r g ≤ 2. We therefore expect that any errors introduced to the initial conditions by dropping contributions of order O((δw dr ) 2 , (δc dr,s ) 2 ) and higher will be ≲ 4%.
Within this prescription the initial conditions, although taking an analytic form that can be easily implemented in CLASS, are still rather cumbersome to write down.Therefore, we instead include a Mathematica notebook (named spartacous initial conditions.nband located in the notebooks/ folder) in our class spartacous code, where we systematically derive the initial conditions of the adiabatic perturbations of all the fluids present in the SPar-tAcous and SPartAcous+ models.These have been added to the perturbations initial conditions routine of the perturbations.cmodule of our class spartacous code.

B Approximation Schemes
In this appendix we describe the different regimes that the DR-iDM system experiences throughout its history, as well as the approximation schemes we employ to facilitate their numerical description within the context of the CLASS code.Initially the DR and iDM fluids have a large momentum-exchange rate, particularly well suited to what is known as the tightcoupling approximation [85][86][87].After the DR decouples from the iDM, and sufficiently deep in the matter-or dark energy-dominated era, the DR behaves as test particles in an external gravitational field, since their average energy density is negligible.This allows for another simplification of the DR evolution equations, called the radiation-streaming approximation [87], which avoids significant computational efforts.
The two subsections that follow deal with the (dark) tight-coupling approximation and the (dark) radiation-streaming approximation of the DR-iDM and DR systems respectively.

B.1 Dark Tight-Coupling Approximation (DTCA)
In the SPartAcous model (and its SPartAcous+ extension) the ψ and χ particles, constituting the DR and iDM, scatter primarily via the t-channel exchange of the massless gauge boson A, under which both are charged.The relevant quantity is the momentum-exchange rate Γ which, as long as the DR temperature T d is larger than the mass m ψ , evolves in time just like the Hubble expansion rate H does during radiation domination [59,77].For the parameter space relevant to the SPartAcous model Γ ≫ H, which means that the DR and the iDM are tightly coupled; their ratio reaches values of the order of Γ/H ∼ 10 11 for α d = 10 −3 and m χ = 1 TeV.
Such a large hierarchy in the timescales relevant for the cosmological evolution of the iDM and DR components is not without its technical difficulties.Indeed, the fact that small differences in the DR and iDM velocity divergences θ dr − θ idm are multiplied by a large rate Γ means that the we are in the presence of a stiff system, a kind of system which is notoriously unstable to numerical methods.
In order to address this problem it is better to think about the DR and the iDM fluids not as separate substances but as a single fluid in equilibrium.Clearly, since they are tightly coupled by the large size of Γ, any deviation from equilibrium will be quickly smoothed out over timescales 1/Γ ≪ 1/H.Mathematically, this means finding the equations for DR and iDM describing their mutual equilibrium, and writing their departures from the same in terms of a series expansion in the small parameter H/Γ.This method, called the tightcoupling approximation, has been successfully applied to the baryon-photon plasma, where Compton scattering tightly couples electrons and photons [86,87].We denote the application of this approximation scheme to our dark sector by the moniker DTCA.
In our CLASS implementation of the SPartAcous model, the perturbations variable dtca on/dtca off controls whether the DTCA is used or not, while the precisions variable dark tight coupling approximation denotes the method by which it is implemented, currently to first-order in H/Γ.We begin our evolution of the DR-iDM system of equations with the DTCA on (unless their interactions are very small, or one or both of those fluids are not present), and turn it off when one of the following three conditions is satisfied: the ratio of Hubble to the momentum-exchange rate becomes large (H/(aΓ) ≥ 0.003), the timescale associated with the relevant mode is also large (k/(aΓ) ≥ 0.0035), or the dark temperature is close to its value at the time of ψ freeze-out (x/x fo ≥ 0.8, with x fo the value of x = m ψ /T d at freeze-out).
Throughout the rest of this section we derive the evolution equations for the perturbations in the iDM and DR fluid in the tightly-coupled limit (Γ ≫ H), starting with their general versions Eqs.(A.12)-(A.15),and following the procedure laid out in Ref. [87], suitably modified to fit our case.

B.1.1 The DTCA Equations
The goal is to derive the equations in the DTCA, when τ c ≪ τ, 1/k.Taking Eq. (A.13) and adding R ×Eq.Defining the useful functions we find that the slip satisfies the exact equation: Note that f /τ is small; f then can serve as an expansion parameter, and the idea is to solve for Θdr-idm perturbatively in f .It can be shown [87] that this perturbative solution is: We will only be interested in the expansion to first order, so we will then only keep ġ in the equation above.Defining the shorthand ∆ ≡ − ) we know that our DTCA equations depend on w dr , c 2 dr,s , Γ, and their derivatives with respect to the scale factor a. However, as shown in Eq. (A.3), Eq. (A.4), and Eq.(A.21), these quantities are more easily expressed in terms of analytic functions of x = m ψ /T d .Therefore, we need to convert all derivatives with respect to a into derivatives with respect to x.This can easily be done by using the chain rule as well as Eq.(A.11) in order to relate x to a.

B.1.3 Summary of DTCA Equations
For convenience, we list the final form of the equations for iDM and DR here.The DTCA equations for θ idm and θ dr are: where we have found the slip Θdr-idm to first order in τ c /τ :

B.2 Dark Radiation Streaming Approximation (DRSA)
Soon after matter-radiation equality, the energy density in the relativistic particles, including DR, becomes negligible.Because of this, once DR is decoupled from the iDM it effectively behaves as a fluid in an external gravitational field.As it turns out, finding the exact behavior of the DR perturbations in this regime is both unimportant and computationally prohibitive.It is unimportant because DR is made up of dark, unobservable particles that lead to no observable signature of their own (e.g. a dark CMB), and it is computationally expensive due to the fast oscillations (hard to calculate precisely) taking place in its perturbations, which arise due to the large value of kτ .In this limit one can find the non-oscillating component of the DR equations using a dark radiation-streaming approximation, which we abbreviate as DRSA, and which is based on the methods described in Ref. [87]. 7n our modified version of CLASS including the SPartAcous model, the perturbations variable drsa on/drsa off denotes whether the DRSA is turned on or off, while the precisions variable dark radiation streaming approximation refers to the method of its implementation, currently a suitably modified version of the approximation employed for the neutrinos in Ref. [87] and the pre-existing idr (interacting DR) fluid in the newest version of CLASS.In our code, the DRSA is only turned on when all of the following conditions are true: the relevant mode is deep inside the horizon (kτ > 44.0), the universe has evolved sufficiently past matter-radiation equality (τ /τ eq > 6.0), the redshift is significantly past z tr , when the step in the DR took place (τ /τ (z = z tr /30) > 5.0), and the DTCA is turned off.

C Numerical Results
In this section we provide comprehensive triangle plots and parameter values obtained from our MCMC numerical results.We illustrate 1D and 2D posterior distributions for the D, DH and DHS dataset in Figures 6-8 respectively.The corresponding best-fit, mean, and ±1σ values for the fits of these datasets are summarized in Tables 1, 4 and 5. Table 7 contain the best-fit, mean, and ±1σ values of the SPartAcous model applied to all the datasets (D, DH, DHS).For the SPartAcous+3 model, please refer to Table 8 for the respective values.Throughout this study, we make the assumption of neglecting the additional relativistic degrees of freedom at the time of BBN, thereby removing their effect on the Standard Model abundance of primordial helium.
(still valid at times well past the step).

Figure 1 .
Figure 1.Feynman diagrams for the most relevant interactions in the dark sector fluid.The left diagram is responsible for keeping the interacting dark matter component χ in equilibrium with the DR.The diagram on the right is responsible for keeping ψ in equilibrium with the gauge boson A.

Figure 5 .
Figure 5. Left: Matter power spectrum for ΛCDM using the best fit point to D and that for SPartAcous and SPartAcous+3 using the best fit point to the DHS data.Right: Ratio of the matter power spectrum of SPartAcous and SPartAcous+3 using the best fit point to the DHS to that of ΛCDM using the best fit point to D.

Table 1 .
Best-fit values of the parameters of the ΛCDM, SPartAcous, and SPartAcous+3 (N df = 3) models, fitted to datasets D and DHS.

Table 3 .
A summary of the fits to the dataset DHS for the ΛCDM, SPartAcous, and SPartAcous+3 models, showing the allowed parameter range at 95% CR for ∆N IR eff , f χ and H 0 .These intervals are defined to be the narrowest interval containing 95% of the integrated posterior density, and have been computed directly from the posterior densities and not using Gaussian fits to the posteriors.Posterior distributions for the new parameters in SPartAcous (left) and SPartAcous+3 (right) fit to the D and DHS datasets.The inclusion of the SH0ES likelihood pushes N IR to larger values for both models, while f χ becomes significantly constrained in SPartAcous once it tries to accommodate a larger H 0 .3.2.1 Fit to DH and Cosmic ConcordanceCurrently, the H 0 tension is much more significant than S 8 .For this reason, we have also made runs including only the DH data combination with all three models.One can see in

Table 4 .
Mean and ±1σ values of the fits of the ΛCDM, SPartAcous, and SPartAcous+3 models to the dataset D. Note that the ±1σ values of the log 10 (z t ) parameter in SPartAcous are "nan".This is because of the very non-Gaussian nature of its posterior.

Table 5 .
Mean and ±1σ values of the fits of the ΛCDM, SPartAcous, and SPartAcous+3 models to datasets DHS.Figure4that the inclusion of the SH0ES prior pulls H 0 to larger values in both SPartAcous and SPartAcous+3, as expected, while it only leads to a very moderate shift for ΛCDM.In SPartAcous, the allowed range for S 8 gets pushed to larger values compared to the results fit to D shown in Figure2, which is a reflection of the fact that the allowed range for iDM in SPartAcous is very small and cannot efficiently lower S 8 .On the other hand, in SPartA-Posterior distributions for H 0 and S 8 for ΛCDM (green), SPartAcous (blue) and SPar-tAcous+3 (red) fit to the DH dataset.Compared to the fit to the D, note that both SPartAcous and SPartAcous+3 get pulled to larger values for H 0 , while ΛCDM sees almost no change.

Table 7 .
Mean±1σ and best-fit values of the SPartAcous model to all the datasets studied in this paper.Note that the ±1σ values of the log 10 (z t ) parameter in SPartAcous are "nan".This is because of the very non-Gaussian nature of its posterior.

Table 8 .
Mean±1σ and best-fit values of the SPartAcous+3 model to all the datasets studied in this paper.