Decaying Dark Matter and Lyman-α forest constraints

Decaying Cold Dark Matter (DCDM) is a model that is currently under investigation regarding primarily the S 8 tension between cosmic microwave background (CMB) and certain large-scale structure measurements. The decay into one massive and one (or more) massless daughter particle(s) leads to a suppression of the power spectrum in the late universe that depends on the relative mass splitting ϵ = (1 - m 2/M 2)/2 between the mother and massive daughter particle as well as the lifetime τ. In this work we investigate the impact of the BOSS DR14 one-dimensional Lyman-α forest flux power spectrum on the DCDM model using a conservative effective model approach to account for astrophysical uncertainties. Since the suppression of the power spectrum due to decay builds up at low redshift, we find that regions in parameter space that address the S 8 tension can be well compatible with the Lyman-α forest. Nevertheless, for values of the degeneracy parameter ϵ ∼ 0.1-0.5%, for which the power suppression occurs within the scales probed by BOSS Lyman-α data, we find improved constraints compared to previous CMB and galaxy clustering analyses, obtaining τ ≳ 18 Gyrs for small mass splitting. Furthermore, our analysis of the BOSS Lyman-α flux power spectrum allows for values τ ∼ 102 Gyrs, ϵ ∼ 1%, that have been found to be preferred by a combination of Planck and galaxy clustering data with a KiDS prior on S 8, and we even find a hint for a marginal preference within this regime.


Introduction
The standard model of cosmology known as ΛCDM is a very successful model in explaining the large scale structure (LSS) of the universe.Cold dark matter (CDM) sits at the heart of this model, causing the typical hierarchical bottom-up structure formation we observe in the LSS by being non-relativistic during the clustering process.Despite the success of CDM, there are still unresolved issues that are hinting that there may be more [1].This sparks interest in different cosmological models that are able to address these issues.Moreover, comparing the process of structure formation within extended cosmological models with LSS observations allows us to constrain fundamental properties of the two large unknowns, dark energy and dark matter, such as the equation-of-state, or the lifetime.
One of the open questions is the so-called σ 8 tension, where σ 8 is a measure of the amplitude of matter fluctuations at a scale of 8 Mpc/h.More specifically, it is convenient to use the parameter S 8 = σ 8 Ω m /0.3 that also includes the matter density parameter Ω m .The tension arises between early universe cosmological data preferring larger values of S 8 , and local, low redshift measurements tending towards lower values when interpreted within the ΛCDM model, with a typical significance of the order of 2−3σ [1,2].Measurements of the cosmic microwave background (CMB) temperature and polarization anisotropies by Planck yield S 8 = 0.834 ± 0.016 [3] which, in this respect, is in agreement with other CMB data like from ACT [4].On the other hand, weak gravitational lensing surveys provide constraints via cosmic shear, e.g. S 8 = 0.759 +0.024  −0.021 from the Kilo-Degree Survey KiDS-1000 [5] and S 8 = 0.780 +0.030  −0.033 from HSC [6].The combination from shear and galaxy clustering from threeyear data of DES yields S 8 = 0.776 +0.017 −0.017 [7] and a combination of shear, clustering and galaxy abundance from KiDS-1000 S 8 = 0.773 +0.028  −0.030 [8], while galaxy cluster counts from SPT-SZ report S 8 = 0.766 ± 0.025 [9] and eROSITA results favour S 8 = 0.791 +0.028  −0.031 [10].Individual measurements might not produce as large deviations, but show a trend being significantly lower compared to CMB data.To account for this, and leaving aside the possibility of a statistical fluctuation, there either needs to be some unaccounted systematic error (see e.g.[11]) or an alternative to ΛCDM featuring a suppression of the matter power spectrum in the k ∼ 0.1 − 1 h/Mpc regime.One model to achieve such a suppression is the Decaying Cold Dark Matter (DCDM) model.It is based on the hypothesis that dark matter can decay on cosmological timescales into secondary dark sector particles.The decay products are assumed to be effectively stable on cosmological scales and, like the dark matter itself, sufficiently weakly coupled to visible matter to escape (in-)direct detection.However, the kinetic energy released in the decay process counteracts the growth of structures and leads to a suppression of the power spectrum.The model has been mainly investigated in two variants: a decay into massless secondaries that act as dark radiation (DR), or into a massless and a massive daughter.Depending on the mass splitting between mother and massive daughter particles, the latter acts as warm dark matter (WDM) being gradually produced in the decay process in the late universe.For both variants, the evolution in the early universe is identical to ΛCDM, thereby preserving its success in explaining the CMB and LSS on very large scales.Both models were studied regarding the S 8 and also the Hubble tension [12], taking CMB, BAO and recently also galaxy clustering data into account.The decay into massless secondaries was investigated e.g. in [13][14][15][16][17][18][19][20][21][22][23], allowing also the possibility that only a fraction f of the dark matter decays.It was found that the decay into purely DR will most likely not be able to solve cosmological tensions and requires minimum lifetimes of around ∼ 200 Gyrs for f = 1.The latest works [22,23] for this model confirm this even more while also providing tight constraints.For lifetimes shorter than the age of the universe, [22] finds f < 2.16% and for f → 1 a lower bound of τ > 250 Gyrs [22,23].As a further variant, also the decay of warm dark matter mother particles into massless dark radiation has been considered [24], but similar to the decay of CDM into massless daughters, this setup was found to neither solve the H 0 nor S 8 tensions [25].
The model with a decay of CDM into WDM and DR, on which we mainly focus in this work, is apart from the lifetime τ described by the mass splitting parameter involving the mass of the mother (M ) and the massive daughter (m) particle.The results regarding the Hubble tension are similar compared to the massless case, implying that it is probably not to be resolved with DCDM [26][27][28] (see also [29,30] for earlier work).The situation for the S 8 tension is however not so clear.While [28] suggests that this tension can also not be addressed, [31], which includes an improved treatment of perturbations, finds that it actually can be lessened for τ ∼ 55 Gyrs and ∼ 0.7% based on Planck CMB, BAO, RSD and SN Ia data.In two follow-up works a newly developed code for much faster computation of the DCDM power spectra is used.This allows for a more in depth analysis, like in [32], where a mild preference for DCDM is found depending on the priors for S 8 .The latest work [22] also includes full-shape information from BOSS DR12 galaxy clustering, and finds that DCDM can ease the S 8 tension even though it is not performing significantly better than ΛCDM when disregarding KiDS data.However, when including KiDS, DCDM is preferred, and the best-fit model occurs for a lifetime of τ ∼ 120 Gyrs and ∼ 1.2%.
Another possibility to study DCDM is via galaxy and halo properties with more regards towards the small scale issues (see e.g.[33]) like the cusp-core problem of DM halos [34].For DCDM with a massive and a massless daughter particle, [35] connects the model to the observed population of Milky Way satellites.Combining numerical and semi-analytic methods, they find constraints at τ 30 Gyrs for 20 v k 200km/s.Here v k is the socalled kick-velocity which is transferred to the daughter particles during the decay, being related to via v k ∼ c for 0.5.The analysis [36] builds up on this work and uses Milky Way satellite galaxies observed by DES, excluding τ < 18 Gyrs for v k = 20 km/s.This probe is extremely sensitive to the low regime which still affects the halo distribution and substructure due to the low virial velocities in dwarf galaxies.
Even when not considering cosmological tensions, it is still interesting to constrain fundamental properties of dark matter like its lifetime via different complementary probes, regarding how little we know about the actual particle nature, and the fact that very few particles are naturally stable [13].Therefore, the degree to which DCDM is compatible with various cosmological and astrophysical observations is worth studying.
In this work we confront DCDM with measurements of the one-dimensional Lyman-α forest flux power spectrum, using data from BOSS DR14 [37].The Lyman-α forest is an important probe for dark matter models that lead to a modification of the power spectrum on scales k 1h/Mpc, which is typically the case for models addressing the S 8 tension.A pecularity of DCDM is that the suppression of the power spectrum occurs at late redshifts, such that it is different for weak lensing, galaxy clustering and cluster number count observations that are sensitive mainly to z 1 as compared to Lyman-α measurements at z ∼ 2 − 4. Therefore, one expects that a larger amount of power suppression is possible at low redshift as compared to models where the power suppression is imprinted already in the early universe, making DCDM a promising model in view of the S 8 tension and Lyman-α constraints.
The main challenge in any Lyman-α forest analysis is the extraction of the actual matter fluctuations from the measured flux power spectrum, requiring a description of the complex intergalactic medium (IGM).In this work, we make use of an effective model that was already used to analyse BOSS data and extensively validated against hydrodynamical simulations for a variety of dark matter models as well as massive neutrino cosmologies in the past [38,39].It contains a number of free parameters that account for the IGM behavior as well as uncertainties from strongly non-linear scales entering via the line-of-sight projection, while taking advantage of the increased reach of a perturbative treatment of the underlying three-dimensional matter distribution at the relevant redshifts z ∼ 3.This allows us to determine robust constraints on the DCDM parameters from the Lyman-α forest on the relatively large scales measured with a high precision by BOSS.
The possibility to address the S 8 tension with DCDM raises the question about an embedding of this scenario in a more complete particle physics framework.We make a first step in this direction by exploring the generalization from two-to three-body decays, that generically occur in models where the involved particle species are fermions.A small mass splitting can be realized naturally by a pseudo-Dirac fermion pair in that setup.
The structure of this work is as follows: In Sec. 2, we give an overview of the formalism of DCDM, the basic background dynamics and the generated power spectrum.Then, in Sec. 3, we review the data set used in this work as well as the effective model and its input and free parameters.Afterwards, in Sec. 4 we present our results within the DCDM parameter space of lifetime and mass splitting.We also set them in context with earlier works with emphasis on the S 8 tension.In Sec. 5 we comment on an extension from two-to three-body decays.Finally, we conclude in Sec. 6.

Formalism
The DCDM model we study comprises collisionless cold dark matter particles that are unstable and decay into two components, DCDM → WDM + DR . (2.1) One is a massive daughter acting as a warm dark matter (WDM) component whereas the other is massless dark radiation (DR, see Sec. 5 for an extension to three-body decays).This model can be described by introducing two new parameters Γ and .The first one, Γ, is the decay width of the CDM mother particle which we usually replace by the decay time τ = Γ −1 .It determines when the decay sets in.The second parameter defined in (1.1) is related to the mass splitting between the mother and massive daughter particle, and characterizes the ratio of energy transformed into DR and WDM.The parameter is also related to the amount of energy that is transformed from rest mass into kinetic energy, and only depends on the mass ratio of mother and daughter and not on the absolute mass values.In the case of m → 0 corresponding to → 0.5, only dark radiation is produced by a decay into two massless daughters.In the opposite case of m → M corresponding to → 0, the daughter particle has almost the same mass as the mother particle and therefore the energy transferred to DR vanishes.Effectively, the decay becomes irrelevant for = 0, independently of the lifetime.Therefore, in both the limits of τ → ∞ as well as → 0 one recovers ΛCDM.
Following [32], we work in synchronous gauge which is comoving with the mother particle.The physical energy-momentum four vectors then take the form P dcdm = (M, 0), P wdm = ( m 2 + p 2 , p) and P dr = (p, − p).Therefore, the physical momentum p ≡ | p| of the daughter particles is fixed by energy and momentum conservation to The density parameter of DCDM today at t 0 is given by an initial density times an exponential factor describing the decay Ω 0 dcdm = Ω ini dcdm e −Γt 0 . (2.3) We can characterize the DCDM energy content by ω ini dcdm = Ω ini dcdm h 2 which is the density of DCDM today if no decay had taken place.At an arbitrary time the density is given by where a −3 gives the additional expansion factor.It can alternatively be written as the time dependent number density Ndcdm times the rest energy of the mother particle.Next, we consider the Boltzmann equations relating the total time derivative of the phase-space distribution function to the collision term governed by the decay.At the homogeneous background level, they are given by ḟdcdm (q, τ ) = −aΓ fdcdm (q, τ ) and ḟwdm (q, τ ) = ḟdr (q, τ ) = aΓ Ndcdm 4πq 2 δ(q − ap 2-body ) , where q = ap is the comoving momentum.The collision term for DCDM only depends on Γ due to the exponential decay.The scale factor a arises from switching to conformal time with dτ = dt/a, and we use the notation where a dot denotes d/dτ .For WDM and DR, the collision term has the opposite sign and is proportional to the number density of DCDM.
The factor 1/4πq 2 takes into account the spherical symmetry and the comoving momentum q is fixed by the delta function to the momentum transferred to the daughter particles ap 2-body = aM due to two-body kinematics.The phase-space-distribution is related to the mean energy density ρ and pressure P by the integrals Here E ≡ m 2 a 2 + q 2 is the comoving energy which is related to the physical energy by E = aE.Making use of these definitions we can transform the Boltzmann equations to ρdcdm = −3H ρdcdm − aΓρ dcdm , with H = aH, Hubble rate H, and equation-of-state parameter ω = Pwdm /ρ wdm for WDM (and assuming an equation-of-state parameter equal to zero for DCDM as usual, which During matter domination the ratio decreases due to a shift from matter to radiation which lowers the energy density.To account for this change while keeping a fixed angular diameter distance to the last scattering surface the ratio has to increase for low redshifts.We note that for small values 0.5 relevant for the S 8 tension the modification of H(z) is negligibly small as compared to the case = 0.3 shown here for illustration.The main change occurs at the perturbation level in that case, see Fig. 3. amounts to neglecting its velocity dispersion).Here, we see the impact of DCDM in the second term on the right-hand side, reducing or adding to the energy densities depending on Γ and .
The resulting evolution of the energy density parameters is shown in Fig. 1 for τ = 40 Gyrs and = 0.006, versus redshift z.Additionally, the cold dark matter density for ΛCDM is shown in blue.For large redshifts, DCDM clearly converges towards ΛCDM as expected.Only for redshifts long after recombination does the decay cause a drop in the DCDM density while WDM and DR are produced in return.Note that the densities approach a plateau for low redshift due to the logarithmic scale and since z → 0 corresponds to the limit t → t 0 .The lifetime τ determines when the deviation from ΛCDM sets in, while controls the relative size of the WDM and DR densities.
The modified background quantities also result in a different Hubble expansion rate as compared to ΛCDM, with the difference being more pronounced the larger (corresponding to a larger fraction of DR) and the smaller τ .The evolution of H(z)/H(z) ΛCDM for a large value = 0.3 is shown in Fig. 2 for various lifetimes, including rather extreme values for illustration.Here we adjusted the dark energy density within DCDM such that the angle under which the first CMB peak appears is identical to the reference ΛCDM model in all cases, i.e. all models feature an identical angular diameter distance at times t τ .The evolution at high redshift is identical to ΛCDM.Once the decay starts to set in, some amount of matter is replaced by radiation, which redshifts faster and subsequently contributes less to the energy content.This leads to a decrease of H(z)/H(z) ΛCDM , lasting until dark energy becomes relevant (marked by the dotted vertical lines).To keep the angular diameter distance to the last scattering surface constant in all models, the dark energy content needs to be larger for shorter τ .This leads to an increase of H(z)/H(z) ΛCDM at low redshift that over-compensates the earlier suppression.While this increase has been considered as a possibility to address the H 0 tension in the past, the increased dark energy content is not consistent with a combination of CMB, BAO and SN Ia data, as discussed previously.For 0.5, the regime that is mostly relevant for addressing the S 8 tension, the modification of H(z) is almost negligible, since most of the energy density is transferred to the WDM component.The increase in Ω Λ additionally leads to an increased late Integrated-Sachs-Wolfe effect in the CMB anisotropy spectrum.Again, for realistic values of τ and this effect is very small, similar to the effect on the Hubble expansion rate.A deviation from ΛCDM occurs instead at the perturbation level, leading to a suppression in the power spectrum at late times.

The DCDM power spectrum
To obtain the total matter power spectrum we solve linear perturbation equations for the mother and daughter particles involved in the decay, coupled to metric perturbations in the usual way [40].For the decaying cold dark matter species it is sufficient to solve the standard continuity equation for its density contrast δ dcdm , given the synchronous gauge choice adopted here.For the daughter particles, the perturbation equations are similar in form to the usual perturbation equations for massless or massive neutrinos, for the dark radiation and warm dark matter component, respectively.In particular, they take the form of a coupled hierarchy for the multipole moments of the perturbed phase-space distribution function.The only difference to the case of neutrinos occurs for the equation of the monopole, that contains an additional source term on the right-hand side due to the decay, given by −aΓ f dcdm δ dcdm [32].Additionally, the homogeneous parts of the distribution functions are time-dependent, as described by (2.5).
In addition to the coupled set of equations for the multipole moments, we consider a fluid approximation for the warm dark matter component as proposed in [32].It amounts to a coupled set of equations for the density contrast and velocity divergence of the massive daughter particle, being the usual continuity and Euler equations.The Euler equation is complemented with an effective pressure term, with sound velocity that depends on the wavenumber k.We adopt the choice given in equation (38) in [32], that was found to reproduce the result of the full coupled hierarchy to a high accuracy, while being computationally much less expensive.
To solve the perturbation equations we use the modified CLASS code presented in [32] by Abellán, Murgia and Poulin. 1 It includes an implementation of the DCDM model in CLASS based on the coupled hierarchy of multipole moments as well as the fluid approximation mentioned above.We checked the agreement of both computational methods, and the dependence on precision parameters like momentum bins or the largest considered multipole moment in case of the computation based on the full hierarchy.We found that in the latter case numerical random oscillations, that may occur on smaller scales depending on the precision settings, are smoothed out for the fluid approximation.Since we are interested in the amount of power suppression in that regime and such oscillations might lead to unphysical artifacts for finite precision settings we found the much faster fluid approximation to be more suitable for our purpose.To generate the linear power spectrum for the ΛCDM model we use the standard CLASS code [41,42].cosmological parameters fixed to Planck 2018 best-fit values.From top to bottom we increase the parameter, increasing the fraction of DR produced in the decay and decreasing the mass of the WDM daughter particle.The parameters and τ clearly have two different effects.For 0.5 the main effect of this parameter is to control the amount of rest mass that is transformed to kinetic energy in the decay, setting the free-streaming scale k fs of the WDM daughter particle.The lighter the particle as compared to its DCDM mother, the larger its kinetic energy, and the smaller is k fs .This imprints a suppression in the power spectrum since structure is washed out for k k fs .Thus, the position of the onset of suppression is determined by with lower values converging to CDM and therefore shifting the suppression more and more to the right.On the other hand, the decay time is responsible for the amount of WDM produced at a given time and thus controls the magnitude of the suppression.Note that the short oscillatory effect on large scales is caused by the switch to the fluid approximation and does not affect scales relevant for our analysis.These are marked by the pink region showing the k ranges that are included in the BOSS Lyman-α data.
Depending on the suppression can start in three regimes.For low values it occurs on scales smaller than those probed by BOSS, implying that the power spectrum does not deviate strongly from ΛCDM (top panel in Fig. 3).On the contrary, for rather high values the suppression occurs already well before the BOSS region, such that the power spectrum is suppressed by an almost constant factor without scale dependence (bottom panel in Fig. 3).As we will discuss below, this shift can approximately be absorbed by astrophysical nuisance parameters within the effective model used here.Therefore, we do not expect a large difference to ΛCDM either, as long as the suppression is not too extreme.The most interesting values are in the range between these limits for which BOSS can directly probe the shape of the power suppression and is therefore most sensitive (middle panel in Fig. 3).If there actually were a preference for DCDM by Lyman-α data we would expect an improvement of the fit over ΛCDM to occur in this parameter region.
Finally, we stress that Lyman-α data are sensitive to the power spectrum at high redshifts (specifically z = 3.0 − 4.2 for our analysis, see below).Due to the exponential decay law, this makes a significant difference in the amount of suppression in the power spectrum as compared to low redshifts.In Fig. 4 we show a DCDM spectrum for τ = 20 Gyrs and = 0.004 at redshifts ranging from z = 0.1 to z = 1000, all normalized to the ΛCDM spectrum at the respective redshifts.For large z, there is no difference because the decay has not set in yet, whereas at the lowest z a large suppression occurs since the amount of WDM has increased and it had more time to wash out structures.The thick lines indicate the relevant redshifts for our analysis at z = 3.0 and z = 4.2, which already show way less suppression of around 10 − 20% compared to the 70% for z = 0.1.
As we are interested in the power spectrum on small scales we have to consider non-linear effects.While our treatment of the complex IGM physics entering the Lyman-α flux power spectrum is discussed in the next section, we first discuss our treatment of non-linearities related to three-dimensional matter clustering.Since we focus on the redshift range z = 3.0 − 4.2, we can take advantage of the much milder non-linearities as compared to z = 0, or equivalently the larger value of the non-linear scale.Indeed, three-dimensional matter clustering is still within the weakly non-linear regime for the redshifts and scales of the BOSS Lyman-α data.Following [38,39], we therefore calculate the one-loop corrections from cosmological perturbation theory and use them as input in our analysis.As in [22], we resort to the conventional EdS approximation for computing non-linear kernels, while it would be interesting to take the exact time dependence for DCDM into account following e.g. the strategy developed in [43].In Fig. 5 we show the relative size of the one-loop corrections P δδ /P linear for redshifts z = 0, 3.0 and 4.2 and τ = 40 Gyrs, = 0.006 (solid) compared to ΛCDM (dashed).We notice two effects: First, as expected, the corrections are much larger at z = 0 since the non-linearities had more time to grow.In contrast, they are sufficiently small within the relevant range of scales and for the earlier redshifts we are interested in, justifying the perturbative treatment at the level of the matter power spectrum (see e.g.[39] for a discussion of the impact of two-loop corrections).Second, for DCDM the corrections are overall smaller compared to ΛCDM within the weakly non-linear regime because of the general suppression of growth, and the quadratic dependence on the linear input spectrum at one-loop order.
3 Lyman-α forest data and model

Data
To infer matter distributions in our universe the Lyman-α forest is a powerful tool for cosmological measurements in a relatively high redshift and small scale regime [44,45].Measurements of the Lyman-α forest are useful in two regards.First, they can help studying the complex photo-ionized hot intergalactic medium.Second, the one-dimensional flux power spectrum is a powerful tool to constrain the underlying matter power spectrum on comparably small scales.High-resolution data, like those measured by HIRES (High Resolution Echelle Spectograph) from the KECK observatory or by MIKE [46,47] are sensitive to small fluctuations and provide a measurement up to very large wavenumbers 0.008 − 0.08 s/km (around 1 − 10h/Mpc).These scales fall in the strongly non-linear regime and are largely affected by Jeans suppression from the IGM pressure.A modelling based on hydrodynamical simulations is indispensable for exploiting these data sets.Mid-resolution data like those provided by SDSS/BOSS [37,48] are sensitive to larger scales, i.e. smaller wavenumbers k ∼ 0.001 − 0.02 s/km (around 0.1 − 2h/Mpc), but compensate for the lower resolution by a much larger number of quasar spectra and correspondingly smaller error bars on the power spectrum.In addition, for the measured redshifts z = 2.2 − 4.6 these scales correspond to the weakly non-linear regime, and are separated by around an order of magnitude from the baryonic Jeans scale.Therefore, BOSS data are potentially amenable to an effective field theory description for which the impact of the complex IGM can be encapsulated in a number of astrophysical nuisance parameters.This is indeed the strategy followed in this work, based on an effective model [38,39] that has been validated with hydrodynamical simulations [49] for ΛCDM and massive neutrino cosmologies as well as a set of IGM parameters, and shown to provide a valid description of BOSS Lyman-α data while being able to absorb IGM uncertainties (temperature, adiabatic index, reionization history) in the effective parameters.We review this model below.
We use the Lyman-α forest flux power spectrum from [37] based on BOSS DR14.We restrict ourselves to redshifts z = 3.0−4.2following the approach in [38], since lower redshifts are more sensitive to non-linearities and higher redshifts to reionization [50][51][52].We show the BOSS data in Fig. 6 where we already included the effective model result for ΛCDM in the linear (dashed) and one-loop (solid) case to give an idea for the later result.The non-linear case clearly fits the data better which is also reflected in the lower χ 2 value with ∆χ 2 = −13.4.

Effective model
To use Lyman-α data for cosmological analyses we need to connect the one-dimensional flux power spectrum to the underlying three-dimensional matter power spectrum.This is often accomplished by running a set of hydrodynamical simulations corresponding to different cosmological and IGM parameters, and then constructing an interpolation among them [53][54][55][56][57].A simulation based strategy is indispensable for small-scale Lyman-α data.However, as mentioned above, for larger scale BOSS data it is also possible to obtain a valid description with a suitable effective model.We use the model discussed in [38,39], that is based on a perturbative description combined with effective parameters that encapsulate the impact of the IGM, designed to work in the BOSS regime far above the baryonic Jeans scale and within the weakly non-linear regime of the underlying matter distribution.As described already above, the model has been validated with a set of hydrodynamical simulations [49] and used to extract neutrino mass bounds as well as constraints on self-interacting dark matter models while marginalizing over IGM uncertainties.We review the setup here.
The Lyman-α photon flux is determined by the transmission F which depends on the optical depth τ via F = e −τ .In particular we are interested in the fluctuations in the transmission spectrum where F is the average transmission.Since on BOSS scales the hydrogen clouds usually do not have large pressure gradients compared to the gravitational forces, and since damped Lymanα systems are removed from the analysis, the underlying matter over-or under-densities can be traced by δ F .Therefore the optical depth and hence the transmission depend largely on the matter fluctuations δ.Additionally, we have to account for the peculiar velocity v p and its gradient along the line of sight, generating distortions in the redshifts of the measured flux power spectrum.This can be described by a dependence on the dimensionless quantity In this case, we can compute the three-dimensional flux where we introduced the parameter β(z) = f (z)b F η (z)/b F δ (z), and P lin (k, z) is the linear matter power spectrum.In the following, we omit the redshift arguments for brevity.
To go beyond the linear approximation, we note that η actually depends on the divergence of the three-dimensional velocity field v given by θ = ∇ v/(aHf ) via η = f µ 2 θ.Using this relationship we arrive at the generalized expression with the additional dependence on the density and velocity power spectra P δδ and P θθ as well as the cross-correlation P δθ .The parameter β can be estimated with the Zel'dovich approximation [58], where it only depends on the adiabatic index γ that in turn is related to the reionization history [59].Since we do not want to impose any prejudice regarding the IGM, we model it with two free parameters α bias and β bias allowing a power-law redshift dependence given by choosing z pivot = 3.0.
In addition, there are also a number of other physical effects we need to account for.First of all, the evolution of baryonic matter is not only determined by gravitational forces but is also tied to its innate pressure.Unlike for dark matter, a collapse cannot happen below the Jeans scale k J = aH/c s , which is related to the sound velocity c s = γT /(µ p m p ) with temperature T , adiabatic index γ and the mean particle mass µ p m p of the IGM.More precisely, we have to look at the filtering scale k F which is given by a time average of the Jeans scale k J [44], For larger k > k F , a suppression is modeled by an exponential factor exp(− (k/k F ) 2 ).A typical value is of order 15 − 20h/Mpc, about an order of magnitude larger than the largest wavenumbers probed by BOSS.Therefore, Jeans suppression has only a minor impact on the relevant scales.Secondly, the spectral lines are subject to thermal broadening, leading to a damping of the flux power spectrum along the line of sight.This damping is also enhanced by other effects like redshift space distortions due to velocity dispersion as well as the finite resolution of the measurements [39,60].To account for this, we include an overall exponential suppression factor of exp(− k /k s 2 ) with the suppression scale k s being mainly determined by the thermal broadening, so k s ≈ m p /T .For example, for a temperature for the IGM of T ≈ 10 4 K, this yields k s ≈ 0.11 s/km, corresponding to around 10h/Mpc.Thus, the effect of broadening is also subdominant within the BOSS regime.
The last additional effect we need to account for are absorption features imprinted by other transitions than Lyman-α on the measured spectrum.The dominant effect stems from SiIII absorption, that can be accounted for by an oscillatory factor with wavelength ∆V = 2π/0.0028km/s due to interference effects [48].These effects are well constrained by the BOSS measurements, and can be described by where we also include the oscillation strength f SiIII = 6 • 10 −3 .The underlying transmission function entering κ SiIII has only a very minor impact and we therefore use a fixed value log( F )(z) = −0.0025(1+ z) 3.7 within (3.7) [38,48].We stress that F is fixed only within the factor described by (3.7).This description was found to be sufficient and allowing the parameters ∆V and f SiIII to vary would not lead to relevant differences in the result since they are very well determined by the oscillatory features in the BOSS data [38,39].
To finally arrive at the one-dimensional spectrum, we integrate along the two directions that are not along the line of sight k leading to where µ = k /k.The integration can be performed for each of the power spectra in (3.4) which leads in turn to the three integrals We observe that the powers of k in the integrals are changing with the powers of µ ∝ 1/k and we also have already included the Jeans suppression that also depends on k.As one can see, the integrals are in principle uncapped and potentially sensitive to extremely non-linear scales.In practice, for I 2 and I 4 these scales are not giving a relevant contribution since the integrand is strongly suppressed for large k due to (i) the inverse powers of k and (ii) the relatively mild enhancement or even suppression of the cross and velocity power spectra relative to the linear spectrum.For I 0 however, a relevant contribution from UV-scales arises.Note that for all k values relevant for BOSS the impact of the UV contribution can be absorbed by an additive constant [38,39].To account for the UV contribution, we therefore include an extra additive counterterm I ct for I 0 , which absorbs the uncertainty from the integration over these UV scales.Again, we allow for a redshift dependence with and add thus two additional free parameters.As mentioned above, we use the δδ, δθ and θθ auto or cross-correlation spectra computed at one-loop order in perturbation theory.Lastly, we add an overall amplitude A accounting for the overall bias factor b 2 F δ in (3.4) that we parameterize by  [3].In case of DCDM, ω cdm is replaced with ω ini dcdm .In addition, we use one massive neutrino species with m = 0.06 eV and two massless species.The σ 8 and S 8 values are model dependent and correspond to ΛCDM. and include the thermal broadening.Finally, with all the additional factors and the integrals, we arrive at In practice, we use overall six free parameters to account for the impact of the IGM and absorb uncertainties from UV modes when integrating across the line of sight.The IGM properties are determined by mainly the velocity bias parameters (shortened now to α b and β b ) as well as α F and β F , while the non-linearities are captured by the counterterm parameters α ct , β ct .As discussed above the impact of the precise value of k s and k F is minor, and can moreover be compensated to a large extent by the other free parameters of the model [38].We therefore use fixed values k s = 0.11 s/km and k F = 18h/Mpc (for a check of the dependence on this choice see [38] and below).
The advantage of this model is that it is very agnostic regarding the complex IGM physics which leads to robust results when marginalizing over the free model parameters.This implies that the model leads per design to conservative constraints on cosmological parameters.

Fitting procedure
To determine the compatibility of a given model with the BOSS DR14 Lyman-α data set [37] we follow a frequentist approach in this work, based on the profile likelihood.We first compute the linear power spectrum, then the one-loop density, velocity and cross power spectra as discussed above, and next the integrals I 0,2,4 for each of the 35 k values and the seven redshift bins z = 3.0, 3.2, . . ., 4.2.For each of these 245 data points we evaluate the difference between the theoretical model and the measured value, and compute a χ 2 taking the statistical uncertainties reported by [37] into account.We then minimize χ 2 with respect to the Lyman-α effective model parameters (3.15).To mitigate the impact of local minima the minimization is performed multiple times scanning different parameter range combinations.
Additionally, we checked the robustness with respect to the model assumptions and repeated the fits with different cutoff scales.This not only tests the cutoff independence (by absorbing the cutoff dependence of I 0 into the values of the counterterm parameters) but provides also an extra check for possible minimization errors.In addition, we performed checks on the dependence on the fixed values for k s and k F as described in section 3.2 and found the result to only deviate below 1%.
We compare our DCDM results to a ΛCDM reference model with the parameters specified in Tab. 1 and results shown in Fig. 6.Additionally, we include one massive neutrino species with m = 0.06 eV and two massless species.This yields a baseline value of χ2 = 192.89 in the one-loop ΛCDM fit.As noted before [38], the absolute χ 2 value should be regarded with care, and we only use the relative χ 2 difference for model comparison.
In order to explore the impact of Lyman-α data, for DCDM we fix the cosmological parameters to the same values as in the ΛCDM reference model, postponing a full analysis to future work.While this choice certainly represents a limitation of our analysis that should be kept in mind, we observe that the corresponding results when taking CMB, BAO and SN Ia data into account do not show any strong degeneracies of the DCDM parameters τ, with the other cosmological parameters [22,32].The main reason is that DCDM behaves identical to ΛCDM around recombination, and when fixing the angular diameter distance to the last scattering surface (ensured by using Θ s as input parameter).The parameter ω cdm is replaced by the equivalent ω ini cdm such that the CDM densities also agree around recombination, before the decay sets in, as do the baryon, photon and neutrino densities.It is therefore reasonable to expect that the cosmological parameters take the same values preferred by Planck in DCDM and ΛCDM, respectively.To span a decent amount of parameter space, we choose 23 different lifetimes τ as well as 28 different values. 2 This results in overall 644 different parameter combinations.

Exclusion bounds
We find that DCDM is not statistically preferred over ΛCDM by BOSS DR14 Lyman-α data even though the fit partly improves over ΛCDM with a best-fit value of χ 2 = 189.77and ∆χ 2 = −3.1 for τ 40 Gyrs and 0.006.The region in the two-dimensional parameter space spanned by the decay rate Γ = τ −1 and the degeneracy parameter that is excluded at 95%C.L. when compared to the ΛCDM reference model (i.e. has ∆χ 2 > 3.841) is shown in Fig. 7.At the bottom left corner, DCDM converges to ΛCDM since it corresponds to large lifetime and low , meaning colder DM.Similarly, DCDM approaches ΛCDM for both τ → ∞ and → 0, i.e. models close to the bottom and left axis are most ΛCDM-like.This is also observed in the fit with a χ 2 value close to the one of ΛCDM when moving in those regions.
At the right side of Fig. 7, the lifetime is very short, down to 1Gyr.The suppression in the power spectrum is therefore very strong in this regime which leads to it being excluded for all 10 −4 .For even smaller the power suppression only arises for scales below those tested by BOSS Lyman-α data.The bound on the lifetime is strongest for values of ∼ 10 −2 − 10 −3 , where the kdependent suppression of the power spectrum due to DM decay falls within the BOSS range.We find that a lifetime of τ 18 Gyrs is approximately excluded at 95%C.L. in this region.For large values 0.1, corresponding to large mass splitting, our analysis yields again weaker bounds on the lifetime since the power suppression is approximately k-independent within the BOSS window for these scales, such that it can be compensated by a corresponding shift in the opacity bias parameter (or equivalently α F within our model parameterization).One could improve the Lyman-α bound in this case by imposing a prior on α F , for example due to a calibration of the opacity bias with hydrodynamical simulations, similar as was done for the neutrino mass analysis in [38].However, as we see below, this large region is already tightly constrained by CMB and BAO data, and therefore we stick to a conservative setup without imposing strong priors on the model parameters.One peculiar feature is the patch on the lower right where one small island in the parameter space is allowed.This is due to a small fluctuation of the χ 2 value such that it is barely not excluded.
Apart from the excluded region, we also show some contour lines of constant S 8 in Fig. 7, assuming Planck 2018 cosmological parameters as given in Tab. 2. Larger values are generated at the lower left corner and lower values at the upper right, due to an increasing level of suppression from DM decay.In addition, in the upper left corner the large mass splitting and low lifetime correspond to the regime where the change in background evolution for DCDM becomes relevant.This leads to a partial compensation of the power suppression ), around the Ly-α best-fit value (marked by a star).The hatched region is allowed at 95%C.L. by Planck, BAO, BOSS DR12 galaxy clustering full-shape (FS), and Pantheon data (referred to as Planck+BAO+FS), adapted from [22].The gray band shows the KiDS range S 8 = 0.759 +0.024 −0.021 [5].BestFit1 and BestFit2 correspond to the best-fit points for Planck+BAO+FS and Planck+BAO+FS+KiDS reported in [22], respectively, see Tab. 2.
due to a decreased Hubble rate at intermediate redshifts (see Fig. 1), that leads in total to a relative enhancement of the growth rate.This explains the shape of the S 8 contours.The grey band shows the region in parameter space where the value of S 8 within DCDM is compatible with the KiDS range S 8 = 0.759 +0.024 −0.021 [5].Thus, we find that Lyman-α exclusion bounds do allow for low S 8 values within the DCDM model.

Allowed region and comparison with CMB and BAO data
In Fig. 8 we show the region in DCDM parameter space that is allowed at 95%C.L. by BOSS DR14 Lyman-α data, corresponding to the region with ∆χ 2 = χ 2 − χ 2 min ≤ 5.991 around the best-fit point at τ 40 Gyrs and 0.006 (marked by a star).We note that for each point in parameter space we minimize χ 2 with respect to the astrophysical parameters of the effective Lyman-α model, i.e. compute the profile likelihood to obtain the allowed region.For comparison we also show the Bayesian posterior obtained from an analysis of Planck, BAO, Pantheon and full shape BOSS DR12 galaxy clustering data (referred to as Table 2: Parameter values and χ 2 for the best-fit points from Planck+BAO+FS (BestFit1) and Planck+BAO+FS+KiDS (BestFit2) reported in [22], as well as the additional contribution to χ 2 from BOSS DR14 Lyman-α data analysed in this work and the Q DMAP statistic quantifying the compatibility with KiDS, when including also Lyman-α data.Due to the proximity of BestFit2 with the Lyman-α best-fit point, the tension is slightly lowered when including the results of this work (compared to Q DMAP = 1.5σ for DCDM without Lyman-α data).For comparison, Q DMAP 3σ for ΛCDM [22].
Planck+BAO+FS) reported in [22].Although the comparison should be treated with care due to the different statistical method and treatment of remaining cosmological parameters, it is instructive to see that the allowed regions are partially complementary to each other.In particular, there exists a region in parameter space for ∼ 10 −3 that is allowed by Planck+BAO+FS, but not allowed by Lyman-α data.Within this regime, the suppression of the power spectrum occurs on scales probed by the BOSS Lyman-α measurements.
In addition, we show in Fig. 8 the region in parameter space favoured by the value of S 8 measured by KiDS as gray shaded region.We see that there is an overlap of Lyman-α, KiDS and Planck+BAO+FS allowed regions.Furthermore, it is intriguing that the KiDS band overlaps with the best-fit point of the Lyman-α analysis.Moreover, this point (marked with a star) is also close to the best-fit point from Planck+BAO+FS+KiDS data identified in [22], shown by the upward pointing triangle in Fig. 8.This indicates that Lyman-α data are well compatible with the DCDM scenario that is preferred for relaxing the S 8 tension, and even slightly favours it.
To further quantify the ability of DCDM to alleviate the S 8 tension, we consider the two best-fit points obtained from Planck+BAO+FS as well as Planck+BAO+FS+KiDS data [22], referred to at as BestFit1 and BestFit2.In Tab. 2 we show the corresponding model parameters as well as the χ 2 values taken from [22].In addition, we compute the contribution to χ 2 from the BOSS DR14 Lyman-α analysis performed in this work.The slight preference of BestFit2 leads to a small reduction of the Q DMAP statistic quantifying the S 8 tension within DCDM from 1.5σ to 0.9σ when taking also Lyman-α data into account.For comparison, within ΛCDM the S 8 tension quantified in this way is ∼ 3σ.
In summary, we find that BOSS DR14 Lyman-α data allow for values of the DM lifetime and mass splitting of mother and daughter particle that are favourable for resolving the S 8 tension between CMB and LSS data.This can be attributed to the strong redshift dependence of power suppression arising from DM decay, such that low S 8 values at z 1 can be compatible with Lyman-α measurements at z ∼ 3 − 4.

Three-body decay
So far, we have studied a two-body decay with one massive and one massless daughter particle.However, realistic dark matter models may allow only for decays with more particles in the final state due to selection rules imposed by underlying (approximate) internal and spacetime symmetries.This frequently occurs for unstable particles in the Standard Model, such as for example for muons and neutrons, or more generally in nuclear β decays.Therefore, we explore the changes when considering three-instead of two-body decays within the dark sector in this section (see e.g.[29] for some discussion in this direction).Specifically, we consider the decay DCDM → WDM + DR a + DR b , ( featuring one massive daughter particle (WDM) with E 1 = m 2 + p 2 1 and now two massless particles (DR a and DR b ) with E 2 = p 2 and E 3 = p 3 , respectively.Whereas the momentum and energy in the two-body decay are fixed as described in (2.2), an additional daughter particle leads to a continous momentum and energy distribution.In the following we outline the changes in the formalism to account for this situation, and propose a simplified mapping that allows one to approximately translate results for the two-body decay to more general scenarios in the most relevant limit 0.5 (i.e.m/M → 1).However, let us start by discussing the general three-body setup.
The total decay rate for a three-body decay is given by where one can choose freely over which two out of the three energies one integrates.Here |M 2 | is the matrix element squared, averaged (summed) over initial (final) state degrees of freedom.For concreteness, and in order to be agnostic about the precise origin of the three-body decay, we consider the spectrum that results purely from three-body kinematics, assuming that the matrix element squared for the decay can be approximated by a constant.The decay spectra are given by the differential cross sections dΓ/dE i for i = 1, 2, 3, respectively.To compute it for the WDM particle (i = 1), we use the minimal and maximal energies E 3 for a given E 1 allowed by energy and momentum conservation, where we introduced the normalization factor N ≡ |M 2 |/(128π 3 ).
For the massless daughters (i = 2, 3), the shape of the spectra as dictated by kinematics are identical.For i = 3 it is given by and for i = 2 one has dΓ/dE 2 = dΓ/dE 3 | E 3 →E 2 .The total decay width takes the form We show decay spectra for the massive (WDM) and one of the massless (DR) daughter particles in Fig. 9.For convenience we show the distribution with respect to the comoving momentum (see below) instead of the energy.For the value = 0.499 close to the maximal value 0.5 (bottom right), the massive daughter is almost massless and also behaves as DR.Thus, all distributions converge towards each other in that limit.For the lower value of = 0.4 (lower left), the difference between the distributions becomes visible, and for 0.5 (upper row) the DR spectrum is peaked at half of the maximal possible momentum, while the WDM distribution always has its maximum at the maximally possible momentum.This can be understood by kinematics and phase-space arguments: the DR a particle reaches its maximal energy (and thus maximal absolute momentum) when the WDM and the DR b particles have momentum vectors that are pointing opposite to DR a .Energy and momentum conservation require the momentum p 3 of DR b to approach zero in that limit.However, the number of available final states in phase-space is suppressed by a factor p 2 3 for p 3 → 0, explaining the suppression of the DR spectrum close to the endpoint.For WDM a similar restriction does not exist, such that it has a spectrum that remains finite at the endpoint.
For later use we also define the average of some quantity X (e.g. the energy of one of the decay products or some power of the momentum) by Note that in the last relation there is no summation, but i is fixed and it holds for any choice of i = 1, 2, 3. Due to energy conservation one has E 1 + E 2 + E 3 = M , and since the massless daughters have identical spectra, one has E 2 = E 3 .We hence use the notation E wdm ≡ E 1 and E dr ≡ E 2 = E 3 for the average WDM and DR energies in the threebody decay.Energy conservation then implies E wdm + 2 E dr = M .Using the three-body spectrum from above, one finds We show the dependence of the average energies on by solid lines in the left panel of Fig. 10.For → 0.5 (i.e.m → 0) one has E wdm /M → E dr /M → 1/3, consistent with the fact that all decay products become massless in that limit.For → 0 (i.e.m → M ) one has E wdm /M → 1 − and E dr /M → /2, i.e. as expected most of the energy goes into the rest mass of the massive daughter, with only little kinetic energy left over.In that limit the average energies agree with the fixed values of the energy in the two-body case (dashed lines in the left panel of Fig. 10).
For later use we also give the average of p 2 /3E over the decay spectrum for WDM (relevant for the pressure), For → 0 one has p 2 wdm /3E wdm /M → 2 /5 in the non-relativistic limit, and for → 0.5 it approaches the relativistic limit E wdm /3, which together with E wdm → M/3 yields p 2 wdm /3E wdm /M → 1/9.This average is also shown in the right panel of Fig. 10.For the fixed values of the momentum in the two-body decay one obtains 2 /(3 − 3 )M .We note that p 2 wdm /3E wdm is thus smaller by a factor 3/5 compared to the corresponding quantity in the two-body case, for 0.5.Let us now discuss how the evolution equations are changed when considering three-body decays.The Boltzmann equations for the homogeneous part of the distribution functions read (compare to (2.5) for the two-body decay) ḟdcdm (q, τ ) = −aΓ fdcdm (q, τ ) , where dΓ/dq i = (dΓ/dE i × dE i /dp i × dp i /dq i )| p i =q i /a = (dΓ/dE i × p i /(aE i ))| p i =q i /a .While the equation for DCDM is unchanged, the equations for WDM and the two DR components contain the momentum spectrum.In particular, the expression in the bracket replaces the delta function δ(q −ap 2-body ) that occurs for two-body decays by the respective decay spectra that are normalized such that dq i Γ −1 dΓ/dq i = 1 for i = 1, 2, 3 when integrating over the comoving momentum.They depend only on the three-body kinematics, while the dependence on the lifetime enters via the prefactor containing Γ = τ −1 .Since the two massless daughter particles have identical decay spectra, it is sufficient to consider the total DR distribution given by fdr (q, τ ) ≡ f dr,a (q, τ ) + f dr,b (q, τ ) . (5.10) By multiplying the distribution functions with the energies and integrating over q we obtain evolution equations for the energy densities, These equations are a generalization of (2.7) to three-body decays, and contain the energy averaged over the decay spectrum.Using that E wdm = M − 2 E dr they are formally identical to (2.7) when replacing in these equations by 2 E dr /M (see (5.7)).This may suggest a mapping of the three-to the two-body scenario at the background level.However, this putative mapping would not be adequate at the level of perturbations, that are relevant for the power suppression, and thus the Lyman-α analysis, as we argue below.However, we note that when considering the most interesting limit 0.5 (i.e.m → M ), the precise value of actually becomes irrelevant as far as the background densities are concerned.The reason is that in this limit the DR energy density becomes negligibly small (with ρ dr ∝ ), while the WDM component becomes non-relativistic with ω ∝ 2 → 0. Furthermore, the average energy entering the evolution equation (5.11) for the WDM energy density approaches the limit E wdm → m M , independent of .Therefore, the time-evolution of ρwdm depends only on the lifetime τ but not on for 0.5.Altogether, the WDM evolution equations (5.11) and (2.7) for the two-and three-body cases become identical to each other for 0.5.Given that this is true also for the DCDM component, this implies that for the relevant background energy densities there is no difference between two-and three-body decays for 0.5.On the level of perturbations, the information on the decay spectrum enters in general via a momentum-dependent collision term in the Boltzmann equations for WDM and DR.In this work we do not attempt a full solution of this case, but rather propose a simple prescription to approximately map the two-body results to the three-body case in the most relevant limit 0.5.In particular, we note that in this limit the DR component becomes irrelevant while the WDM part can be well approximated by a fluid description, as discussed above.The suppression of the power spectrum is then encoded in the effective sound velocity following [32], where it was argued that the adiabatic sound velocity complemented by a small correction factor yields results in congruence with the full Boltzmann hierarchy for 0.5 in the two-body case.In the following we assume that the dominant effect of going from twoto three-body decays can be captured by the modification of the adiabatic sound velocity, leaving a more detailed analysis to future work.It is given by where p is the pseudo-pressure [61], being a higher moment of the distribution function.This result is a generalization of the one given in [32] to the three-body case, and contains the average over the decay spectrum of the quantity given in (5.8).In the non-relativistic limit → 0 the contribution from pseudo-pressure becomes suppressed by a relative factor 2 , as does ω.Furthermore, inspecting the evolution equation for the pressure itself we find that for small (5.13)This ratio approaches 3/5 for → 0. Following the arguments from above we can approximately account for this reduction in the sound velocity parameter by re-interpreting the results of our Lyman-α analysis obtained in the two-body case for some given set of parameters and τ as constraints that apply also to the three-body decay for a mass spectrum We use an effective model for the one-dimensional Lyman-α forest flux power spectrum, that is designed to work on scales of BOSS data, far above the baryonic Jeans scale and within the weakly non-linear regime of the underlying matter density field.The effective model contains in total six free parameters that account for uncertainties from the IGM as well as the impact of non-linearities, and has been validated with hydrodynamical simulations for various cosmological models in the past.
We find that for certain values of the mass degeneracy parameter = (1 − m 2 /M 2 )/2, BOSS Lyman-α data yield constraints on the dark matter lifetime τ that are stronger compared to a combination of Planck, BAO, SN Ia and galaxy clustering data.The lower bound reaches τ 18 Gyrs for ∼ 0.1 − 0.5%.Interestingly, the lifetime τ ∼ 10 2 Gyrs and ∼ 1% that is favoured for relaxing the S 8 tension is also marginally preferred by BOSS Lyman-α data, as compared to ΛCDM.The S 8 tension according to KiDS data is around 3σ within ΛCDM, and is reduced to 1.5σ for DCDM.When including also BOSS Lyman-α data it is slightly further reduced to ∼ 1σ.While the hint for this mild preference is intriguing, the main conclusion is that Lyman-α data are compatible with dark matter decay being a possible explanation of the S 8 tension.
Apart from the question whether the S 8 tension is due to systematic effects, it would be interesting to investigate realistic and well-motivated models of dark matter decay.As a first step in this direction, we provide a mapping from the two-body decay scenario to a more general three-body decay that is valid in the limit 0.5 (i.e. for m → M ), taking the different phase space into account.The mapping amounts to a rescaling of the parameter.This allows one to easily translate constraints derived for the two-body case to models with three-body decays.We find that the best-fit Lyman-α scenario corresponds to m/M = 0.994 for two-body decays, and m/M = 0.992 for the three-body case.

Figure 1 :
Figure 1: Background evolution of the DCDM (orange), WDM (green) and DR (red) density parameter for τ = 40 Gyrs and = 0.006 in comparison to a conventional CDM scenario without any decay (blue).For large redshifts (corresponding to t τ ) the decay is irrelevant, while subsequently the DCDM density drops below CDM, and WDM as well as DR are produced.The plateau at low redshift is due to the logarithmic z-axis.

Figure 2 :
Figure2: Hubble rate H(z) depending on redshift normalized to ΛCDM for various lifetimes τ and a large value = 0.3 of the mass splitting.The dotted colored lines correspond to the equality of matter and dark energy and the pink region to the redshifts spanned by the Lyman-α data used in this work.During matter domination the ratio decreases due to a shift from matter to radiation which lowers the energy density.To account for this change while keeping a fixed angular diameter distance to the last scattering surface the ratio has to increase for low redshifts.We note that for small values 0.5 relevant for the S 8 tension the modification of H(z) is negligibly small as compared to the case = 0.3 shown here for illustration.The main change occurs at the perturbation level in that case, see Fig.3.

Figure 3 :
Figure 3: DCDM power spectrum normalized to ΛCDM at z = 3 and for = 0.0003, 0.008, 0.05 from top to bottom, respectively.Each panel shows various lifetimes τ .The pink region indicates the range of BOSS Lyman-α data.Since the parameter determines the kinetic energy of WDM produced in the decay, it controls the onset of suppression in the power spectrum.Larger corresponds to more kinetic energy, shifting the suppression scale to the left.The decay time is in turn responsible for the steepness of the suppression since it controls the amount of WDM at a given time.

Figure 4 :
Figure 4: The DCDM power spectrum at τ = 20 Gyrs and = 0.004 at z = 0.1, 3.0, 4.2, 10, 100, 1000 relative to ΛCDM at the respective redshifts.The relative suppression increases with time due to the ongoing decay.The thicker lines indicate the redshift range between z = 3.0 and z = 4.2 used in our Lyman-α analysis.

Figure 6 :
Figure 6: BOSS DR14 data for the one-dimensional Lyman-α forest flux power spectrum at redshifts z = 3.0 − 4.2 used in this work, with errorbars as well as the linear (dashed) and one-loop (solid) best-fit effective model in the ΛCDM case.
) with velocity v p and comoving coordinate x p along the line of sight.Now, we can expand the optical depth fluctuations at first order as[44] δ τ = b τ δ δ + b τ η η,where we introduced the bias parameters b τ δ and b τ η .Applied to the transmission, we get in turn δ F = b F δ δ + b F η η with the new parameters b F i = log(F )b τ i .In linear approximation the velocity gradient would be proportional to the density contrast via η = f µ 2 δ, where the µ 2 = k 2 /k 2 factor arises from only taking the contribution along the line of sight k , and f = d ln D/d ln a is the growth rate, with linear growth factor D(a).

Figure 7 :
Figure 7: Excluded parameter space by BOSS DR14 Lyman-α data at 95%C.L. (dark blue and hatched).The red lines indicate contours of fixed S 8 within the DCDM model and the gray band shows the region consistent with the value S 8 = 0.759 +0.024 −0.021 favoured by KIDS [5].

Figure 9 :
Figure 9: Momentum distribution of WDM and one DR particle in the three-body decay DCDM → WDM + DR a + DR b for the three different values = 0.01, 0.4, 0.499 of the mass splitting parameter, corresponding to m/M = 0.04, 0.45, 0.99.The solid grey line shows the maximum momentum a particle in the three-body decay can have, being identical to the fixed momentum in the two-body case.For → 0.5, the distributions converge towards each other since all particles effectively behave as DR.For lower values the distributions differ due to phase-space suppression (see text for details).

Figure 10 :
Figure 10: Left: Average energy E wdm of the WDM (orange) and both DR particles 2 E dr (green) in the three-body decay (solid lines), depending in the mass splitting parameter = (1 − m 2 /M 2 )/2.For comparison we also show the fixed energies for a two-body decay, given by (1 − )M and M , respectively (dashed).Right: same for the average of p 2 wdm /3E wdm for three-body decay (solid) and the two-body case (dashed, given by 2 /(3 − 3 )M ).

Table 1 :
Cosmological parameters for the ΛCDM and DCDM model used in our analysis, fixed to Planck values