Dark radiation from the primordial thermal bath in momentum space

Motivated by the stunning projections for future CMB surveys, we evaluate the amount of dark radiation produced in the early Universe by two-body decays or binary scatterings with thermal bath particles via a rigorous analysis in momentum space. We track the evolution of the dark radiation phase space distribution, and we use the asymptotic solution to evaluate the amount of additional relativistic energy density parameterized in terms of an effective number of additional neutrino species ΔN eff. Our approach allows for studying light particles that never reach equilibrium across cosmic history, and to scrutinize the physics of the decoupling when they thermalize instead. We incorporate quantum statistical effects for all the particles involved in the production processes, and we account for the energy exchanged between the visible and invisible sectors. Non-instantaneous decoupling is responsible for spectral distortions in the final distributions, and we quantify how they translate into the corresponding value for ΔN eff. Finally, we undertake a comprehensive comparison between our exact results and approximated methods commonly employed in the existing literature. Remarkably, we find that the difference can be larger than the experimental sensitivity of future observations, justifying the need for a rigorous analysis in momentum space.

1 Introduction Motivated extensions of the standard model (SM) often incorporate new light and feebly interacting degrees of freedom [1][2][3].Sometimes, these elusive particles are even stable.Whether their stability follows from the symmetry structure of the underlying theory (e.g., new global symmetries) or their lifetime is just much longer than the age of the Universe, it is important to quantify accurately and precisely their relic abundance as well as the associated detectable signals.Explicit realizations include frameworks where these new particles serve as dark matter candidates [4][5][6] with mass bound to be above the keV range if they are produced via processes involving particles from the primordial bath [7][8][9][10][11].
In this work, we focus on new physics scenarios with weakly-coupled particles that never contribute to the energy density of non-relativistic matter.They are instead always in the relativistic regime and therefore provide an additional radiation component in the early Universe.We probe the energy density stored in relativistic degrees of freedom at two key moments in cosmic history.The first one is Big Bang Nucleosynthesis (BBN), which is not a sudden process, and it is sensitive to the radiation energy density when the Universe was one second old and the temperature T BBN ∼ 1 MeV.Another crucial probe is the last scattering surface, namely the moment when protons and electrons form neutral hydrogen, and photons could travel undisturbed afterward.This happens much later than BBN although the Universe is still relatively young: its age is approximately 380,000 years, and the photon temperature T CMB ∼ 0.3 eV.The photons from the last scattering surface that we detect today are dubbed the Cosmic Microwave Background (CMB) radiation.These two cosmic probes provide two independent measurements of the radiation energy density at early times.Incidentally, the measured number of effective relativistic degrees of freedom does not have to be the same since things may change as the Universe evolves between T BBN and T CMB .
The only SM relativistic degrees of freedom at both BBN and CMB are photons and neutrinos.However, theories motivated from the top down often feature the presence of one (or more) light particle X that can be present in the relativistic regime and with a significant cosmic abundance.Historically, the presence of this dark radiation is quantified in terms of an effective number of additional neutrino species.More explicitly, the energy density of relativistic particles at the time of CMB formation can be written as The first term in the bracket accounts for both photon polarizations, whereas the second one is proportional to the statistical Fermi-Dirac factor of 7/8 and to the fourth power of the neutrino to photon temperature ratio T ν /T γ = (4/11) 1/3 .This relation defines the dimensionless factor N eff that quantifies the number of effective neutrino species.If we stick to the SM alone, one may naively conclude that this quantity is just the number of fermion generations (in this case three).In fact, a careful study of neutrino decoupling leads to a slightly larger number [12][13][14][15][16].
The recent analysis in Ref. [17] reports the value N SM eff = 3.043.How well do we measure the radiation energy density in the early Universe?Or, equivalently, how well do we measure N eff ?BBN and CMB have similar sensitivities at the moment.The analysis of light element abundances, where the effective number of neutrino species is usually denoted by N ν , puts the constraint N BBN ν = 2.889 ± 0.229 [18,19].From the CMB side, the most constraining result comes from the Planck satellite, N Planck eff = 2.99 ± 0.17 [20].Other ground-based experiments, such as ACT [21] and SPT-3G [22], probe N eff with a lower sensitivity and are consistent with Planck.The measured values at BBN and CMB are consistent with each other, and they are both consistent with the SM prediction.This does not have to be the case necessarily because extra dark radiation components can alter the SM prediction, and things may change between BBN and CMB formation.As a result, the presence of dark radiation in the early Universe is severely limited.
Thus N eff is a powerful probe of SM extensions containing light degrees of freedom [23], and this is the first motivation for our work.The second one is supported by the expected improvements for the future.At the moment, the best sensitivity available is the one due to Planck and it results in σ(N eff ) ≃ 0.17.In the near future, the first significant advancement will be due to the Simons Observatory [24] with σ(N eff ) ≃ 0.05.Later on, CMB-S4 [25-27] will reach σ(N eff ) ≃ 0.03.There are also futuristic proposals [28] with the target σ(N eff ) ≃ 0.014.This astonishing program will make N eff an even more robust new physics probe.Even in the absence of new physics, the increased sensitivity will make the difference between the naive expectation for N SM eff and its actual value detectable.This known SM case is an unambiguous illustration of how important it is to provide not only accurate but also precise predictions.
Our main goal is to develop the framework needed to provide robust predictions for the amount of dark radiation within SM extensions.We focus on a specific injection source: decays and/or scatterings involving degrees of freedom belonging to the primordial thermal bath.Our study is based upon two rather mild assumptions: there was a thermal bath before BBN, and the dark radiation candidate X has feeble interactions with the thermal bath itself (e.g., SM particles).The first one, which is not supported by observations, is a reasonable extrapolation of what we probe at BBN [29].Although our analysis assumes that the thermal bath dominates the energy budget, it is straightforward to extend our results to cases where the thermal bath is a sub-dominant component but it is still responsible for the production of dark sector particles [30][31][32][33][34][35][36][37][38].The second one is rather common in UV complete models.
This thermal production is in some sense unavoidable, and the way it works is rather simple from a qualitative point of view.Processes from the bath dump X particles into the early Universe, and if they happen often enough they may even lead to the thermalization of dark radiation.Whether this happens or not, at some point the Universe is so cold and diluted that nothing can happen to X particles anymore, they just free-stream uninterrupted.A rapid way to estimate the amount of dark radiation is via the instantaneous decoupling approximation.This quick method leads to answers that are definitely wrong by an amount larger than future experimental sensitivities, but also larger than the current ones if one considers production around the time of the QCD crossover where the number of SM effective degrees of freedom changes rapidly.Alternative methods based on ordinary differential equations that allow to tracking of the X's number or energy densities improve the accuracy of the predictions.However, as we discuss in this paper, they are still based upon some assumptions that are not always justified.And the well-known example of SM neutrinos illustrates how important is to investigate the decoupling epoch carefully.
We present here a general method to investigate the thermal production of dark radiation in momentum space based on an integro-differential Boltzmann equation that allows us to track the X's distribution in the phase space.A study based entirely on a phase space analysis has multiple benefits.We list and explain below the main novelties of our approach.
• Non-Thermalized Relics.We always consider initial conditions at early times with no additional dark radiation.If the interaction strength is not enough, these particles may never thermalize in the early Universe but still be around today with a detectable abundance.Our methodology accounts carefully for parameter space regions where X never reach thermalization.This has to be contrasted with other methods conventionally adopted in the literature.The instantaneous decoupling approximation cannot work in this case for obvious reasons.Even if we adopt improved procedures based on integrating the Boltzmann equation in momentum space to track a given moment of the distribution function (e.g., the number density), there are always some assumptions about reaching a thermal profile.
• Decoupling epoch.Even if thermalization is achieved, the decoupling is neither instantaneous nor momentum-independent. Methods based on ordinary differential equations go beyond the instantaneous decoupling, but they still need to assume that the phase space profile of the distribution function remains thermal after decoupling.As we will show, different momenta decouple at different times and residual spectral distortions are somewhat unavoidable.These effects are not dramatic and do not alter the predicted amount of dark radiation by orders of magnitude, but they still lead to effects larger than the CMB-S4 sensitivity and therefore need to be accounted for.
• Quantum statistical effects.There is no reason that forces us to employ classical statistical distributions in our framework for bath particles or dark radiation.We discuss throughout this work what difference it makes to treat classically or quantummechanically the different actors.Here, we emphasize how the assumption of classical statistics is never needed in our case.As we explain in the paper, this assumption is not needed if one employs the instantaneous decoupling approximation either.However, if we adopt the procedure where ordinary differential equations describe the evolution of number or energy densities, one is forced to employ the Maxwell-Boltzmann statistics for the dark radiation to reduce the original integro-differential equation to ordinary ones.Bath particles can always be treated quantum-mechanically instead.
• Feedback on the thermal bath.The Boltzmann equation in momentum space describes accurately how dark radiation particles populate the associated phase space, but it is only part of the full story.Production processes also subtract energy from the thermal bath and it is not possible to study the two systems independently.We complete our framework with an additional Boltzmann equation accounting for the evolution of the energy density of the thermal bath.Besides its evolution following from the geometry of the Hubble expansion, the right-hand side of the Boltzmann equation contains a collision term that quantifies the energy exchanged with the dark radiation sector.It is enough for our purposes to track the energy density of the primordial bath without going to the phase space in this case; SM interactions are quite efficient at early times and ensure thermalization.
We structure the presentation of our results as follows.Sec. 2 contains the general formalism to track the phase space evolution of any thermally produced dark radiation particle.This work focuses on a broad class of production channels, the ones where only one dark radiation particle is produced in the final state of a generic process involving bath particles.We dub this scenario single production, and the Boltzmann equation analysis simplifies substantially as explained in Sec. 3. Furthermore, we develop the formalism for two explicit production channels: two-body decays and binary collisions.Technical details about phase space integrals appearing in the collision operators of the Boltzmann equation are relegated to App. A. Our results are collected in Sec. 4 where we solve for the phase space distribution and integrate over the particle momenta to compute the effective number of neutrino species as a function of dark sector masses and couplings.We perform the analysis for both classical and quantum statistics, and we highlight the differences when there are any.An undeniable advantage of a phase space analysis is the opportunity to track the evolution of each momentum bin, and in particular, analyze how different momenta decouple at different epochs in the cosmic history.This is the origin of spectral distortions we find in our solutions, and ultimately the reason why our predictions differ slightly from the ones obtained via standard methods.Thermalization and decoupling in each momentum bin are scrutinized in the App.B. We put our study into the context of the current literature in Sec. 5 where we compare our results with others derived via approximate methods.These standard procedures, which we review in App.C, are the instantaneous decoupling approximation, and Boltzmann equations for the number density and the energy density of the dark radiation candidates.We explain in the paper all the assumptions needed to recover these approximate methods from the complete phase space framework.Remarkably, we find parameter space regions where our results differ from conventional methods by more than the sensitivity of future experiments.Finally, we summarize our findings and discuss future applications of our general framework in Sec. 6.
B 1 (P 1 ) The most general process for thermal production of dark radiation considered in this work: n bath degrees of freedom (B i ) collide and lead to a final state with m bath and ℓ dark radiation (X) particles, respectively.We also include production via decays (n = 1).

General framework for a phase space analysis
We develop in this Section the general formalism to investigate the thermal production of dark radiation in momentum space.The primordial bath is assumed to be a thermalized gas of weakly-coupled particles B i with temperature T filling the early Universe.Our framework is applicable if such gas is made only of SM particles as well as for extensions with additional particles in thermal equilibrium.We consider the most general process producing an arbitrary number of dark radiation particles X from initial and final states with n and m bath degrees of freedom, respectively.The process is illustrated in Fig. 1 and it explicitly reads We put in parenthesis the associated four-momenta and we use the following parameterization: bath and dark radiation particles have four-momenta P µ i = (E i , ⃗ p i ) and K µ r = (ω r , ⃗ k r ), respectively.Furthermore, we denote the magnitude of the spatial momentum with the same letter without the vector symbol (i.e., |⃗ Our goal is to track the evolution of the phase space distribution (PSD) function f X , and there are two distinct reasons why it changes.On the one hand, there is the dilution effect due to the Hubble expansion, as described by the Friedmann-Lemaître-Robertson-Walker (FLRW) geometry.On the other hand, processes such as the one in Fig. 1 create and destroy X particles and this is an effect controlled by the underlying dynamics.
In the remaining of this Section, we write down the collision operator for the process in Eq. (2.1) and the Boltzmann equation describing the PSD evolution with time.However, the evolution of f X cannot be determined independently on the one of the thermal bath because of the energy exchanged between the two sectors.For this reason, we derive another Boltzmann equation governing the evolution of the thermal bath that accounts for the feedback of dark radiation particles, and we provide a system of Boltzmann equations that can be solved once we couple them with the Friedmann equation quantifying the Hubble rate.Finally, we provide a recipe to compute the additional number of effective neutrino species from the numerical solution of the Boltzmann system.

Dark radiation evolution in momentum space
Without any loss of generality, we focus on the X particle with four-momentum K 1 in Fig. 1.Homogeneity and isotropy of the FLRW geometry ensure that the PSD can only depend on the magnitude of the spatial momentum k 1 and the cosmic time t.Thus the Boltzmann equation governing the evolution of f The Liouville operator L acting on the PSD takes care of the geometry of the expanding universe whereas the collision operator C encodes the dynamics mediating the process in Fig. 1.For an FLRW background, the equation reads The collision operator accounting for Eq.(2.1) explicitly reads (see, e.g., App.B of Ref. [10]) (2.3) with P in = n i=1 P i and P fin = n+m j=n+1 P j + ℓ r=1 K r .The Lorentz invariant relativistic phase space factors are defined as dΠ i = g i d 3 p i /[(2π) 3 2E i ] and dK r = g X d 3 k r /[(2π) 3 2ω r ] with g i and g X the number of internal degrees of freedom for B i and X, respectively.The transition amplitudes for the production (M → ) and destruction (M ← ) are squared and then averaged over initial and final states with the appropriate symmetry factors for identical particles.So far, we have kept the matrix element for the direct and inverse processes distinct.However, if the interactions mediating the process in Eq. (2.1) conserve CP or equivalently T, we have the identity From now on, we assume CP conservation.
Bath particles follow equilibrium PSDs that can be either Bose-Einstein (BE) or Fermi-Dirac (FD).We will discuss throughout this paper when quantum statistical effects can be neglected and the Maxwell-Boltzmann (MB) statistics can serve as a valid approximation.Chemical potentials are negligible at high temperatures, and the explicit expressions for the distribution functions of bath particles in thermal equilibrium read It is possible to simplify the collision operator in Eq. (2.3) by plugging the explicit expressions given in Eq. (2.4) for the PSDs of bath particle.Upon using the identity where in the last step we impose the conservation of energy, we find 1 This is the most general form for the collision operator once we account for CP-conserving processes producing an arbitrary number of dark radiation particles.

Feedback on the thermal bath
Even without dark radiation, the bath energy density ρ B is not constant in time but it gets diluted by the expansion.Processes such as the one in Fig. 1 also impact how ρ B evolves, and we quantify the energy exchanged between visible and dark sectors.Strictly speaking, we would need to couple Eq. (2.2) with the analogous ones for each B i to account for the dark radiation feedback on the thermal bath evolution.In practice, bath particles are always in equilibrium and it is sufficient to focus on their energy densities that are collectively summed to give the total bath energy density.The starting point to determine the ρ B evolution is the integration of both sides of Eq. (2.2) over the phase space The first equality follows from the identification of the dark radiation energy density ρ X as a function of time by integrating the PSD times the energy over the phase space Energy conservation imposes that the quantity appearing on the left-hand side of Eq. (2.7) has to be, with the opposite sign, equal to the right-hand side of the equation accounting for the red-shift of the bath energy density.Thus the bath evolves according to Energy density ρ B and pressure p B are related via the equation of state p B = w B ρ B .

Final system of Boltzmann equations
The evolution of the dark radiation PSD f X and the bath energy density ρ B is controlled by the Boltzmann equations in Eqs.(2.2) and (2.9), respectively.A generic collision operator C appears on the right-hand side of both equations and for the process in Eq. (2.1) and illustrated in Fig. 1 we can use the simplified form provided in Eq. (2.6) if CP is conserved.Both equations are differential in the cosmic time t and they feature PSDs that depend on the bath temperature T .Thus we need a relation connecting time and temperature, and this 1 One can derive the same result by invoking the detailed balance principle that for the CP conserving case is provided by the Friedmann equation accounting for the Hubble rate H ≡ ȧ/a, where a(t) is the FLRW scale factor.To summarize, the full evolution is described by the system (2.10) The last equation of the system assumes that there is no other contribution to the energy density besides the ones from the thermal bath and the dark radiation.In other words, we are assuming that the thermal bath itself is dominating the energy budget at early times as is the case for a standard cosmological history that extrapolates the BBN snapshot of the Universe.Thus dark radiation production happens during a radiation-dominated epoch.
As long as collisions happen at a significant rate, the evolution is nontrivial and the system above has to be solved numerically.However, once number-changing processes stop being effective the phase space evolution is just due to free-streaming and the dark radiation energy density defined in Eq. (2.8) scales as ρ X ∝ a −4 .This is also the scaling of the photon energy density below the e ± threshold (at temperatures below the electron mass).Thus once we are well below the MeV scale, the ratio ρ X /ρ γ reaches a constant value, and this ratio is crucial to quantify the amount of additional radiation at the time of CMB formation.As it is conventionally done in the literature, we isolate the SM contribution to N eff and we split the total amount by defining N eff = N SM eff + ∆N eff .We can identify the contribution to ∆N eff by direct comparison between the two terms on the opposite sides of the last equality of Eq. (1.1), and it explicitly reads (2.11)

A concrete scenario: single production
We study quantitatively a scenario that we call single production: dark radiation is produced via processes like the one in Fig. 1 with only one X particle in the final state (i.e., ℓ = 1).This is often the leading production channel in concrete models due to a severe price to pay in terms of a small coupling for each final state X.
A prominent class of this kind of candidates are axion-like-particles (ALPs) with dimension 5 interactions to standard model fields suppressed by a high scale Λ. Interactions with gauge fields and fermions are both allowed whereas couplings to the Higgs currents can be redefined away.More concretely, for an ALP field ϕ we have the interactions The first class of operators contains the ALP coupled to standard model gauge fields F µν i with the sum running over the three gauge groups.We follow the standard conventions adopted in the literature and we normalize these couplings proportionally to the fine structure constant α i of the associated gauge group.The second contribution contains derivative couplings for the ALPs to vector and axial-vector currents of standard model fermions.If Λ is much higher than the temperature under consideration, single production is clearly the leading production channel.Flavor-diagonal operators (i.e., i = j) mediate only production via binary scatterings.If flavor-violating interactions are allowed, two-body decays are also possible.The QCD axion [39][40][41][42] is a case motivated from the top-down by the strong CP problem.We suppress the index labeling X particles and we have the four-momentum K µ = (ω, ⃗ k) with the relativistic dispersion relation k = ω.The Boltzmann equation for f X simplifies in this case since the PSD appears in the collision operator only via the combination where we identify the equilibrium PSD f eq X (k) for X without the chemical potential.This identity holds for both Bose-Einstein and Fermi-Dirac statistics, and it is valid also when we neglect completely quantum effects and adopt Maxwell-Boltzmann statistics for X.The Boltzmann equation describing the PSD's evolution takes the simple form where we introduce the collision term We stress the difference between the collision operator C[f X (k 1 , t)] and the collision term C(k, t) defined in Eqs.(2.6) and (3.4), respectively; the former is a functional of the PSD f X , the latter is a regular function of the physical momentum k and the cosmic time t.The result in Eq. (3.3) is the most general Boltzmann equation describing the single production of dark radiation, and there is no assumption whatsoever about the functional form of the unknown PSD f X (e.g., we do not have to assume kinetic equilibrium).The quantity C(k, t) has the units of energy and it can be interpreted as the collision rate for a given momentum bin.We focus now on two specific production channels.

Two-body decays
The first channel we investigate is X production via two-body decays B 1 → B 2 + X, and the general collision term in Eq. (3.4) takes the following form The partial decay width for this process results in The bath particles B 1 and B 2 have masses m 1 and m 2 , respectively, and the degeneracy factors g 2 and g X are due to our conventions with the squared matrix element averaged also over final states.Finite B 2 mass effects are captured by the phase space factor The squared matrix element is constant (monochromatic final states), and we can take it outside of the integral and trade with the partial decay width.We spell out the remaining integrations over the phase space in App.A.1, and we write the collision term as follows Here, we employ the bath temperature T as a time variable.The dimensionless function D(k, T ) is defined by Eq. (A.11), and the explicit expressions for different statistics are shown in Tab. 1.
We compare the collision rate normalized as (g X /g 1 )C 1→2 (to remove the g i factors) with the expansion rate.The latter is provided by the Friedmann equation, and we evaluate it for a thermal bath of SM particles with the equation of state parameter w B = w SM as provided by Ref. [43] (see also Refs.[44,45]).Fig. 2 shows the temperature dependence of the ratio (g X /g 1 )C 1→2 (k, T )/H(T ) for m 1 = 1 GeV (left panel) and m 1 = 1 TeV (right panel).For both cases, we fix the dimensionless ratio Γ 1 /m 1 = 10 −14 and m 2 = 0.The different colors correspond to different choices of physical momenta, k/T = {0.1, 1, 10}, and for each color (i.e., for each chosen physical momentum) we show five different lines corresponding to the cases listed in Tab. 1.We fix the values of the physical momenta in units of the temperature T because after decoupling they change just due to the cosmological red-shift, k ∝ a −1 .Likewise, the bath temperature cools down in agreement with entropy conservation and, up to corrections proportional to the variation of the bath degrees of freedom, it scales as T ∝ a −1 .Thus the ratio k/T is approximately constant after decoupling and it can be thought of as a given comoving momentum.We notice from Fig. 2 how different comoving momenta decouple (i.e., C 1→2 ≲ H) at different temperatures.We also observe how statistical effects may alter the temperature at which the dark radiation thermalizes with the bath, while the decoupling temperature is basically independent of the statistics.

Binary scatterings
The second example we consider is single production of dark radiation from binary scatterings.Once we identify a specific collision (e.g., B 1 B 2 → B 3 X), we have to include all processes that are allowed by crossing symmetry and sum over them to get the collision rate (3.8) (g X /g (g X /g Collision terms for two-body decays divided by the Hubble parameter as a function of the bath temperature T .We fix the B 1 mass to 1 GeV (left) and 1 TeV (right), and the final state products to be massless.The decay width is set to Γ 1 /m 1 = 10 −14 .We consider three different physical momenta k: 0.1 T (purple); T (magenta); 10 T (yellow).The different line styles correspond to the choices for the particle statistics explained in the legend.C 2→2 /H/(g The notation is the same as the one adopted in Fig. 2. Each contribution can be obtained from the general result in Eq. (3.3), and we find (3.9) The phase space integrals can be computed with the techniques presented in App.A.2 where we develop a formalism valid for arbitrary interactions, spectra, and statistics.For the sake of illustration, we focus here on a simplified framework where we set the (dimensionless) scattering matrix elements to a constant value independent on the momenta, and this quantity is the same for each permutation: We choose a mass spectrum with only one overall mass scale: B 1 and B 2 have the same mass If we adopt MB statistics for all bath particles, the collision term can be computed analytically up to an integral over the Mandelstam variable s (see Ref. [10]), and it reads The factor of 2 in the second term is due to the fact that Fig. 3 shows the comparison between the normalized collision rate C 2→2 /(g 1 g 2 g 3 ) and the Hubble expansion rate H.We choose two different values for the overall mass scale, m = 1 GeV (left panel) and m = 1 TeV (right panel), and we set the squared matrix element to the constant value |M| 2 = 10 −11 .Lines in the plot colored differently correspond to different physical momenta k in units of the temperature T : k/T = {0.1, 1, 10}.The reason for this normalization of the momenta is the same as the one for the decays.For each color, we show four different cases corresponding to different choices for the statistics of the bath particles.We notice how different comoving momenta decouple (i.e., C 2→2 ≲ H) at different temperatures and how statistical effects may significantly alter the decoupling temperature, unlike in the decay case.Compared to the decay case, the collision term for scatterings is a shallower function of temperature.Thus if some momenta reach equilibrium, they remain coupled for a longer time than in the case of decays, making it easier to retain kinetic equilibrium.

Results
In the previous section, we have introduced the single production scenario and shown how the Boltzmann equation for the dark radiation PSD reduces to the simple form in Eq. (3.3).In particular, the collision term C(k, t) as defined in Eq. (3.4) is a regular function of the physical momentum k and the cosmic time t.We have derived the explicit collision term expression for two specific production processes, two-body decays and binary scatterings, and shown their temperature dependence in Figs. 2 and 3. We are now ready to feed the Boltzmann equation with these collision terms, solve it, and use the output to quantify ∆N eff .
First, we rewrite the Boltzmann system in Eq. (2.10) for dark radiation single production We find it convenient to employ the FLRW scale factor a(t) as the "time variable", and therefore we trade the cosmic time t with a(t).As a reference, we define a I to be the scale factor at some high-temperature scale T I well before production processes become efficient.
For the scenarios considered in this work, physical results do not depend on this choice if T I is sufficiently high: production via decays is always IR dominated, and the same is true for production via collisions if the underlying interactions are renormalizable (see, e.g., Ref. [4]).
We employ as the evolution variable the dimensionless ratio A ≡ a/a I .After the processes accounted for by the collision term C(k, t) stop being efficient, X particles red-shift with the expansion and their physical momenta decrease as k ∝ a −1 .Thus the use of comoving momenta is particularly advantageous since they scale out the effect of the Hubble expansion, and for this reason, we introduce Likewise, we define the comoving energy densities of the thermal bath R B and of the dark radiation R X .We find it convenient to identify the following dimensionless quantities Strictly speaking, the former is not always constant as the Universe expands since we have mass threshold effects that alter the SM equation of state.The latter can be computed at each value of A (i.e., at each moment in time) once f X is known since we just need to integrate Eq. (2.8) over the momenta and scale out the effect of the expansion where we take advantage of the isotropy of the Universe to perform the integral over the solid angle in polar coordinates.
The Boltzmann system in terms of the new variables reads If one looks at the first equation of the system, the one for the evolution of f X , it seems that different Fourier modes are not coupled to each other and the evolution of each momentum bin can be treated independently.In other words, one may erroneously conclude that this is not an integro-differential equation because the PSD for X does not appear in the collision term.In fact, the different Fourier modes for X particles are still coupled.The equation for the evolution of f X contains on the right-hand side the Hubble parameter that is set by the energy content of the universe; this is quantified by the Friedmann equation, the third one in the Boltzmann system.If one looks at Eq. (4.4), it is manifest how the X's contribution to the energy density results from the integral of the PSD over all the Fourier modes.Furthermore, the right-hand side of the second equation of the system also contains an integral over all the modes.The Boltzmann equation is inevitably still integro-differential.This has to be contrasted with the case of dark matter freeze-in where dark matter particles never get close to thermalization, their contribution to the energy budget is negligible and therefore one can take the standard cosmological history result for the Hubble parameter [10].
In what follows, we solve the Boltzmann system numerically2 .We bin the comoving momentum in the range q ∈ [0.005, 20] with N q = 64 logarithmically-spaced bins and we evaluate the N q functions of the temperature C(q, T ) for two-body decays and binary scatterings.Thus we deal with N q + 1 Boltzmann equations, N q for the unknown function f X (q, A) of the associated momentum momentum bin, and a last one for R B .At each "time step" the PSD is interpolated in momenta and used to compute integrated quantities such as R X , and the right-hand side of the second equation in the system (4.5).The temperature of the thermal bath T at each value of A can be found by solving the equation We consider a SM thermal bath with thermodynamic variables as provided by Ref. [43].Initial conditions describe the system when A = 1 (i.e., a = a I ).In the beginning, we assume that there is no dark radiation, f X (q, A = 1) = 0.The generalization to the case where there is some initial X population (e.g., from inflaton decays) is straightforward.On the contrary, the initial bath energy density is non-vanishing and follows from the relation in Eq. (4.6), R B (A = 1) = π 2 30 g * ρ (T I ).We integrate the Boltzmann system up to a value of the scale factor A F corresponding to a temperature T F well below the electron mass and above the recombination temperature.This is necessary because the effective number of relativistic degrees of freedom does not change anymore afterward, g * ρ (T ≪ T F ) = const, and the ratio between the X's and bath's energy densities is constant.Thus we reach our final goal of evaluating the contribution to ∆N eff that results in In the denominator of the last fraction, we identify the energy density of the two photon polarizations since this is the quantity that enters in Eq. (2.11).We find it convenient to identify a dark radiation temperature.One possible definition for this quantity would be to solve the differential equation for the PSD, derive the energy density ρ X as defined in Eq. (2.8), and under the assumption that f X ≃ f eq X define a dark radiation temperature for the dark radiation inverting the relation in Eq. (C.3).This clearly works if the interaction strength is large enough to bring X particles to equilibrium, but it does not work in the most general case.A more general definition, which we choose in our work, does not rely upon the assumption of thermalization.We define the dark radiation temperature from the width of the PSD computing its second moment as follows The overall normalization coefficients are understood by considering the limit where interactions are strong enough to ensure thermalization.In this regime, X particles are in equilibrium with the bath and the PSD reaches the expressions If we apply the definition in Eq. (4.8), we find that the dark radiation temperature is equal to the one of the thermal bath, T X = T .Thus our general definition reproduces correctly the case when thermalization is achieved.
The next two sub-sections contain numerical solutions of the system in Eq. (4.5) for the two production mechanisms considered in this work: two-body decays and binary scatterings.First, we neglect quantum statistical effects for all the particles involved in the production and treat them all with MB statistics.This simplification allows us to rely on analytical expressions for the collision terms, and we provide predictions for ∆N eff obtained from a precise and quasi-analytical phase space analysis.In the following sub-section, we include quantum statistical distributions for all the degrees of freedom participating in the production processes.We consider both bosonic and fermionic dark radiation, and we compare the outcome of this rigorous quantum statistical analysis with approximated results where we neglect either the statistics of all the particles or just the statistics of the bath particles.

Classical Statistics
We consider the decay B 1 → B 2 +X with all particles obeying the MB statistics, and we focus on the simple case g 1 = g 2 = g X = 1.The mass of the final state particle B 2 is assumed to be negligible, and the only mass scale in the game is m 1 .Another key quantity is the decay width Γ 1 for this production process, and we trade it for the dimensionless ratio Γ 1 /m 1 .
Fig. 4 shows the numerical solutions for the PSD.The two different panels correspond to different values for the mass of the decaying particle, m 1 = 1 GeV (left) and m 1 = 1 TeV (right).These two values exemplify a generic case since they allow us to study the production of dark radiation around the QCD crossover where the number of relativistic degrees of freedom changes abruptly, and at very high temperatures where this quantity is constant, respectively.For each panel, the three different rows have different values for the dimensionless ratio Γ 1 /m 1 as given in the legenda.These three choices are meant to represent three distinct physical scenarios: freeze-in of dark radiation (blue lines), thermal equilibrium barely reached (orange lines), and dark radiation strongly coupled to the thermal bath (red lines).
The lines in our plots always show the PSD f X (q, A T ) as a function of the comoving momentum q and at fixed values of the scale factor A T that is understood as the scale factor value when the bath temperature was T .The PSDs appearing on the vertical axis are normalized in a way that we find convenient.The multiplicative factor (g X /2π 2 )q 3 is what is needed to identify the integrand whose integral leads to the dimensionless comoving energy density R X (see Eq. (4.4)).The quantity R X appearing in the denominator is helpful if one wants to compare the PSD at different moments; initially, before the production processes do their job, the value of f X is rather small, and it grows later on.Thus we are comparing the shape of the PSD at different times and not the overall normalization.The choice of the variable on the horizontal axis, where the dark radiation temperature T X (A) given in Eq. (4.8), is also convenient because we want to investigate when we achieve thermalization.We do it by comparing our results with the reference thermal distribution f eq X (q, A) = exp[−T I q/(A T X (A))], and this definition corresponds to f eq X = exp[−k/T X ] with k the physical momentum related to the comoving momentum by Eq. (4.2).Within these conventions, the equilibrium distribution is fixed in time and can be effectively used as a reference to quantify deviation from thermalization.
Each plot in Fig. 4 shows the reference thermal distribution as a dotted black line, and different 'snapshots" of the numerical solution for the PSD at moments of the expansion history when the bath temperature was T = {100, 1, 0.01} × m 1 .These three values are chosen to represent moments when the production of dark radiation has not started yet, is most effective, and has terminated, respectively.The early-time PSD (dot-dashed line), at T = 100 m 1 , appears to be much wider than the equilibrium one but with a peak at lower comoving momenta.The PSD when T = m 1 (dashed lines) is usually the closest possible that it can get to the equilibrium one throughout the cosmic evolution as the interaction rate is maximal and, in the strongly coupled cases, it matches exactly the equilibrium expression.Finally, the late-time PSD (solid line) at low temperatures, T = 0.01 m 1 , is the free-streaming one: X particle production has ceased at this point and thus afterward the (normalized) PSD is frozen in time.Already from this analysis, we can see noticeable differences between the final PSD and the equilibrium one, flagging the presence of spectral distortions.These distortions  are due to the fact that the interaction turns off at different times (temperatures), for the different momentum bins, as already pointed out in Fig. 2. We discuss in further detail this effect and its relation with the kinetic equilibrium assumption in App.B. The discussion for production via the scattering B 1 +B 2 → B 3 +X (and permutations) is very similar to the decay case so we outline here the main peculiarities.As before, we employ the MB statistics for all particles and we take g 1 = g 2 = g 3 = g X = 1.We take the same mass spectrum chosen in the previous section, with B 1 and B 2 having the same mass m and B 3 massless.In Fig. 5, we show the normalized PSD obtained numerically for m = 1 GeV (left) and m = 1 TeV (right).As already done in the previous section, we assume that the squared matrix element |M| 2 is constant and we present results for the three different values 10 −18 (upper row, blue lines), 10 −14 (middle row, orange lines), and 10 −10 (lower row, red lines).These three choices are again representative of three different scenarios where thermalization is never, barely, and completely achieved, respectively.The different snapshots are chosen at the same reference temperatures T = {100, 1, 0.01} × m, and the thermal distribution is again presented as a black dotted line.This figure shows in which regimes and at what temperatures thermal equilibrium is achieved.The behavior is qualitatively very similar to the decays, but overall the distributions are much closer to the reference thermal one, signaling how in the scattering case, where three processes are involved in producing dark radiation, spectral distortions are more difficult to produce.However, we emphasize how these results are derived under the assumption that the squared matrix element is constant and it is not obvious that it still holds if the underlying dynamics leads to transition amplitudes with a Figure 6: Predicted values of ∆N eff for X production via decays (left) and scatterings (right) obtained with the Maxwell-Boltzmann (MB) statistics for all particles involved in the production process.For the decay case, we show how ∆N eff depends on the dimensionless ratio Γ 1 /m 1 for different choices of the decaying particle mass.For the scattering case we show the dependence is on the constant squared matrix element |M| 2 for the collision.
non-trivial momentum dependence.A well-known example is the QCD axion because its Nambu-Goldstone nature implies coupling proportional to the derivative of the axion field itself, and the various momentum bins feel interaction strengths that are rather different.Regardless of the production mechanism, once the PSD is obtained we can integrate it over the momenta to find the energy density and the resulting ∆N eff .In the left panel of Fig. 6, we show the result of ∆N eff as a function of the parameter Γ 1 /m 1 for different choices of the decaying particle mass m 1 .Notice how the rapid changes of relativistic degrees of freedom in the plasma is responsible for a huge variation in the value of ∆N eff from m 1 = 100 MeV to m 1 = 5 GeV.Likewise, the right panel of Fig. 6 shows the prediction for ∆N eff as a function of the constant scattering amplitude for different values of the mass scale m (chosen to be the mass of B 1 and B 2 whereas B 3 is assumed to be massless).

Quantum statistics
We relax the assumption of classical statistics and we treat all particles involved in the production process, bath degrees of freedom B i and dark radiation X, quantum mechanically.The first thing to decide is the spin of the dark radiation particle X, and we consider two cases: spin-0 (bosonic) dark radiation produced via the decay of a Dirac fermion; spin-1/2 (fermionic) dark radiation described by a Weyl field and produced via the decay of a spin-0 bath particle.The spin of the final state bath particle is determined accordingly.We label the first case FD(4)→FD( 4)+BE (1), where the number in brackets denotes the number of internal degrees of freedom (i.e., g 1 = g 2 = 4, g X = 1), and the second case BE(1)→FD( 2)+FD(2) (i.e., g 1 = 1, g 2 = 2, g X = 2).As done before, we consider the decaying particle to be the only massive particle.The collision term for this case is provided by Eq. (3.7), and the explicit expressions for the function D(k, T ) for the different statistics are shown in Tab. 1. Fig. 7 shows the prediction for ∆N eff for the two different mass values m 1 = 1 GeV (left) and m 1 = 1 TeV (right).We compare the full numerical solutions obtained by employing the proper quantum statistical distributions (solid lines) with the approximate analysis done in the previous section where all particles are treated with MB statistics (dotted lines), and with a different approximate analysis where we adopt the MB distribution for all bath particles but we treat X particles with the correct quantum statistics (solid dotted lines).For each plot, we show also the error of each approximate analysis with respect to the exact solution.We notice how treating all particles with MB statistics works up to a 10% error over the whole range of coupling strengths g 1 Γ 1 /m 1 with the absolute error becoming detectable by futuristic ∆N eff probes ∼ 0.01.Instead, the analysis where we treat only the dark radiation quantum mechanically suffers from the same 10% error at small couplings, but agrees very well with the exact result at larger values (g 1 Γ 1 /m 1 > 10 −18 for m 1 = 1 GeV, and g 1 Γ 1 /m 1 > 10 −15 for m 1 = 1 TeV).Thus neglecting the bath particle statistics in the case of decays does not impact the calculation of ∆N eff as long as the coupling is sufficiently strong.Nevertheless, the error made using this approximation for smaller couplings is not detectable by future ∆N eff probes.The underlying physics responsible for these results is illustrated in Fig. 2 where we can observe how different statistics change the temperature at which the interaction starts to be relevant (i.e., the transition from C/H < 1 to C/H > 1), but have no impact on the decoupling epoch (i.e., the transition from C/H > 1 to C/H < 1).Thus decoupling for a given momentum k happens at the same temperature regardless of the statistics (see the discussion in App.B for further details, and in particular Fig. 12).
We consider now bosonic and fermionic dark radiation produced via collisions.The former is assumed to be a spin-less particle produced via collisions of other spin-0 degrees of freedom, and the latter a is Weyl fermion produced via collisions of spin-0 bosons and a Weyl fermion.For both cases, only the bosons B 1 and B 2 are massive and have equal mass m.Fig. 8 compares the results obtained taking into account the quantum statistics of all the particles involved (solid colored lines) with the approximation of treating all particles with the MB statistics, and the improved method with the correct statistics just for the X particle.The error at small couplings (g 1 g 2 g 3 |M| 2 < 10 −15 for m = 1 GeV and < 10 −12 for m = 1 TeV) is at least 10% up to even 100% for all cases.For larger couplings, the improved method is overall better than treating all particles with the MB statistics.However, differently from the decay case, the scattering rate is the sum of three different processes and the total scattering rate is much more sensitive to the statistics with respect to decay.This can be understood referring to Fig. 3, and in more detail in Fig. 12 of App.B: considering different statistics and same momentum, we see that the decoupling temperature (i.e., C 2→2 /H drops below one after being larger than one) depends on the statistics.Thus assuming MB statistics in the collision term C 2→2 is not an approximation as good as in the decay case even at large couplings.

Comparison with Approximate Methods
We take one step back and consider approximate methods to estimate ∆N eff that rely upon ordinary differential equations.In this section, we always adopt MB statistics for both bath particles and dark radiation.The main goal here is to compare the errors that we make when we adopt these approximate methods, and choosing statistical distributions that offer analytical collision rates helps our case.Furthermore, we have seen how the full quantum interaction rates are well approximated by the MB statistics.A careful review of the approximate methods is provided in App. C. Before comparing with the full results obtained via a momentum space analysis, we summarize them briefly here.
• Instantaneous Decoupling.The underlying assumption is that X particles thermalize and eventually decouple.The departure from thermal equilibrium is a sudden phenomenon at a specific temperature T D that is determined via a comparison between the MB(1) → MB(2), FD(2) MB(1) → MB(2), FD(2)  We consider bosonic (upper panels) and fermionic (lower panels) dark radiation.Solid lines show the exact results where we employ quantum statistical distributions for all particles, and we also report results from approximate analysis: all quantum statistical effects neglected (thick dotted black lines), and only X's quantum statistics accounted for (thick dot-dashed black lines).In the lower part of each panel we quantify the errors made by the approximate results.
interaction and the expansion rates.At temperatures lower than T D , dark radiation particles free stream, and one can evaluate their energy density at recombination by employing entropy conservation.This approach is outlined in App.C.1.Once we determine T D , the resulting contribution to ∆N eff results in ∆N eff ≃ 0.027 g X ξ ρ X 106.75/gB * s (T D ) 4/3 , where the dimensionless coefficient ξ ρ X accounts for the statistical properties of X and it is given explicitly in Eq. (C.3) for the different cases.Within this framework, one can account for quantum statistical effects for both bath and dark radiation particles since the interaction rates are evaluated for equilibrium distribution.
• Tracking the number density.This method allows us to deal with the details of the decoupling epoch, and to account for situations where dark radiation never achieves thermalization.The mathematical tool is an ordinary differential equation for the X's MB(1), MB(1) → MB(2), FD(2) MB(1), MB(1) → MB(2), FD(2)  number density n X that is derived from the integro-differential Boltzmann equation in momentum space.We review carefully this derivation in App.C.2.Along the way, we need to assume that assume that ϕ particles are in kinetic equilibrium and described by MB statistics.No assumption is necessary about the statistics of bath particles.The Boltzmann equation takes the final form given in Eq. (C.16).Furthermore, we need to convert the solution for n X into a correspondent energy density ρ X in order to evaluate ∆N eff .A further assumption that is required for this last step is to assume a thermal shape in chemical equilibrium for the dark radiation's PSD.
• Tracking the energy density.The assumptions of kinetic equilibrium and MB statistics for the dark radiation are still unavoidable if one wants to deal with an ordinary differential equation.However, tracking ρ X would have the indisputable benefit that it does not require any conversion to evaluate ∆N eff .In particular, we do not need to assume a thermal shape in chemical equilibrium for the dark radiation PSD anymore.We review this method in App.C.3 where we also explain how to account consistently for the presence of X's in the energy and entropy densities as done in the main text.Methods analogous to this one were used by Refs.[13,[46][47][48][49].
Among these three methods, the last one where we track the energy density has its obvious advantages but it also comes with its own limitations.In particular, it is not possible to go beyond X's kinetic equilibrium and to include quantum statistical effects for X without dealing with the problem in momentum space.
Before presenting the explicit comparison, we find it useful to introduce another approximate procedure that still relies upon working in momentum space.We remind the reader that the general system to solve for an accurate momentum space analysis is the one given in Eq. (2.10).This fourth approximate method provides a simplified version of it.
• Phase space without feedback on the bath.The approximation for his case counts on the fact that the contribution of the dark radiation to the total energy density is negligible, and more precisely of a percent order with respect to the thermal bath for even the largest couplings we consider here.Thus it is reasonable to consider a standard cosmological history with the Hubble parameter given by H = (π g B ⋆ρ (T ) T 2 )/(3 √ 10M Pl ).Consequently and consistently with this approximation the evolution of the bath energy density is determined by the above relation.Thus the system collapses to a single equation for the dark radiation PSD, namely Eq. (3.3), and we can employ entropy conservation to use the temperature as the evolution variable Here, we define the comoving momentum q = (k/T )(g * s (T )/g * s (m)) 1/3 as in Ref. [10].
Fig. 9 shows the comparison for the case of the two-body decay for the two usual benchmarks m 1 = 1 GeV (left) and m 1 = 1 TeV (right).The limitations of the instantaneous decoupling approximation are clear since different comoving momenta decouple at different temperatures (see Fig. 12 in App.B).Moreover, if the coupling is not strong enough to achieve thermalization, the method is not applicable at all.
Small couplings are troublesome also for the n X method.For the m 1 = 1 GeV case, the error can be sizable with the most striking differences appearing in the low coupling region where dark radiation never thermalizes.In general, we notice that to reproduce the exact prediction for ∆N eff with the n X method we need higher couplings.It is interesting to investigate further the origin of this result.When we convert the final number density n X to an associate energy density and compute ∆N eff , we have to assume an equilibrium distribution for the dark radiation.This assumption is clearly not completely justified at small couplings, and this leads to an underestimate of the dark radiation energy density.The details of this conversion are described in App.C.2, and the final relation that allows us to perform this operation can be found in Eq. (C.20).Thus the conversion is done using the temperature of the dark radiation T ρ X found from the energy density, and this has to be contrasted with the T f X defined from the PSD and defined in Eq. (4.8).Therefore the error between the n X method and the phase space solution at small couplings can be estimated by the ratio between these two temperatures, ∆N err eff (n between these temperatures is elaborated in App.B and shown explicitly in Fig. 14).This is illustrated by the orange lines in the lower plots in Fig. 9. Instead, the error at larger coupling is due to the fact that we are solving for an integrated quantity of f X , hence implicitly assuming kinetic equilibrium.Integrated methods lead to an underestimate of the decoupling temperature with respect to phase space solution in the case of decays (see Figure 9: Comparison of ∆N eff predictions for production via two-body decays.All particles are treated with MB statistics, and we consider the two benchmarks values m 1 = 1 GeV (left) and m 1 = 1 TeV (right) for the decaying particle mass.We vary the value of Γ 1 /m 1 over the horizontal axis, and we report the prediction for ∆N eff on the vertical one.The different horizontal lines identify current bounds as well as the sensitivities of future CMB probes.Exact results from a momentum space analysis are shown as black solid lines.The approximate methods we compare with are: instantaneous decoupling (blue squares), tracking the X number density (orange up-sided triangles), and tracking the X energy density (green down-sided triangles).We also show the results of a phase space solution without including the dark radiation feedback on the thermal bath (red crosses).The lower panels show the absolute errors on ∆N eff from the exact solution; the orange lines show the estimate of the error due to the conversion between n X and ρ X .Fig. 12).Moreover, even if this effect is sub-leading, we are not accounting for the feedback of dark radiation into the SM plasma.We will further explain these two aspects below.
The ∆N eff found via the ρ X method agrees very well with the result from the PSD computation at low couplings, but it shows again a sizable (larger than the future experimental sensitivity) difference at couplings large enough that thermal equilibrium is reached during the evolution.Specifically, the ρ X solution overestimates the amount of dark radiation.This again can be explained from the fact that kinetic equilibrium is not attained during the decoupling phase, in which the interaction rate drops below the Hubble parameter.Fig. 13 of App.B shows the factor (f X /f eq X )/(R X /R eq X ) as a function of the momentum at different epochs.Well after decoupling, for T ≪ m 1 , this quantity depends on the momentum bin also for large couplings.From a superficial examination of this figure, one could conclude that since this ratio is much larger for small couplings the ∆N eff computed from ρ X and f X should roughly agree for large coupling and disagree strongly for small couplings.This would lead to the wrong conclusion that for couplings for which thermality is attained at least once in the thermal history of the dark radiation solving for the energy density will give the same ∆N eff as solving for the PSD, while the result will differ for smaller couplings.Fig. 9 shows exactly the opposite behavior.The reason for this resides in the different importance of the back-reaction terms (the one proportional to f X /f eq X in the equations for the PSD and the one proportional to R X /R eq X in the equation for R X ) at different values of the coupling Γ 1 /m 1 .Indeed for Γ 1 /m 1 ≪ 10 −16 we are in the "freeze-in"-like case (R X ≪ R eq X , f X ≪ f eq X ) for which we can safely neglect the back-reaction term in the Boltzmann equations (since X particles never reach equilibrium) which, respectively become Instead, for Γ 1 /m 1 > 10 −16 the equations depend on back-reaction terms Thus if dark radiation particles thermalize and kinetic equilibrium after decoupling is not attained (as shown by Fig. 13 at large momenta), the backreaction in the energy density equation is stronger than the backreaction in each momentum bin and one explains the suppression in the energy density obtained from the integration of the PSD solution with respect to the solution for ρ X .This discussion shows that, while the solution of the Boltzmann equation is accurate as long as X particles are in thermal equilibrium, the relic energy density of dark radiation particles cannot be determined correctly by solving for the energy density because spectral distortions determine a difference in the right-hand sides of the Boltzmann equations, that albeit very small, are exponentially enhanced in the solution of the differential equation.Instead, if dark radiation is produced by a pure freeze-in, i.e.R X ≪ R eq X at all times, then these spectral distortions have no consequence on the right-hand sides of the Boltzmann equations, thus the solution for f X and R X coincide at the percent level.
The above discussion advocates for a complete solution in phase space to avoid overestimating the amount of dark radiation, in particular in the strongly coupled regime.Inspecting closely Fig. 9, we notice that the most interesting region is the transition between smallcoupling and large coupling, around Γ 1 /m 1 ∼ 10 −18 as the error between the approximate methods and the exact is large, experimentally detectable and the value of ∆N eff is not excluded by Planck [20], but will be probed by future searches.
The case of binary scatterings is overall analogous to the two-body decay case but with a few peculiarities visible in Fig. 10.Solving for integrated quantities, both n X and ρ X prove to be a much better approximation than in the case of decays.This is because in the scattering case, the kinetic equilibrium is more easily assured, as there are three contributing processes.Because of this, a momentum which decoupled in a process can still be coupled with another.This results in an overall flatter C/H when compared to decays, something we knew already from the comparison between Fig. 2 and 3. Thus when momenta reach equilibrium they remain coupled for a longer time, ensuring kinetic equilibrium on a higher level.Furthermore, different comoving momenta decouple at roughly very similar temperatures and these temperatures are close to the decoupling temperature estimated through Γ X = H, as clear from Fig. 12, signaling that the relic distribution is closer to the equilibrium one.
It is important to emphasize that for this study we have only considered theories where the scattering matrix element is constant and independent of the external momenta.This is not the case in several explicit models, and the momentum dependence in the matrix element will alter the result in Fig. 3.In particular, the plot for C/H will exhibit less flat lines, and larger spectral distortions are likely.Thus it is not obvious that the remarkably good agreement with approximate methods would persist with non-constant matrix elements.

Conclusions
The vast amount of present and cosmological upcoming data offers a unique opportunity to probe beyond SM physics.Several motivated frameworks feature the presence of light and weakly-coupled particles.The strong motivation for their existence is behind the astonishing experimental effort hunting for them.Terrestrial searches produce solid bounds that do not rely upon any assumption about the cosmological history, but they also come with their own limitations.They require coupling strengths large enough that these particles are indeed produced, and they are able to constrain their interactions with light degrees of freedom.The QCD axion is an emblematic case since terrestrial searches are able to probe directly only axion couplings with photons, electrons, and nucleons.Cosmological surveys provide a complementary tool to test these frameworks because of the high temperatures and densities achieved through the expansion history.Weakly-coupled particles can still be produced copiously at early times.Moreover, the primordial thermal bath contains all SM (and possibly beyond the SM) degrees of freedom in thermal equilibrium, and this provides us with a way to test the interactions of these new feebly interacting states with heavier SM particles.
In this work, we focused on new hypothetical degrees of freedom that are coupled to SM particles and light enough that they provide an additional contribution to the radiation part of the energy budget at early times.We quantified their abundance in terms of an effective number of additional neutrino species, ∆N eff , as it is conventionally done in the literature.The fraction to the energy budget due to radiation can be measured at two key moments in cosmic history: at the time of BBN when the Universe is approximately one second old, and at the last scattering surface when the Universe is much older (380,000 years) but still relatively young compared to its present age.The remarkable agreement of SM predictions with the observations at both BBN and CMB was the first motivation for our work since ∆N eff is a very powerful new physics probe.At the same time, the astonishing sensitivity of current observations and the even more impressive future projections motivated our refinement of the calculation for ∆N eff .With only a few exceptions, this quantity is obtained via the solution of ordinary differential equations that derive from integrations over the phase space of the Boltzmann equation in momentum space.Or, even less accurately, the abundance of new relativistic degrees of freedom is estimated by assuming an instantaneous decoupling.We took one step back to the momentum space Boltzmann equation, and we developed a general framework to compute ∆N eff without any approximation.In particular, we never assumed chemical or kinetic equilibrium for dark radiation at any time, and we tracked the complete distribution in phase space at any moment of the cosmological history.This allowed us to deal also with small couplings such that dark radiation never thermalizes with the primordial bath, and to track the details of decoupling when it thermalizes instead.An additional benefit of our framework was that we never had to employ Maxwell-Boltzmann statistical distributions and we could fully incorporate quantum statistical effects.Finally, we coupled the Boltzmann equation in momentum space with another one describing the evolution of the bath energy density, and we accounted for the energy exchanged between the visible and invisible sectors.
For the sake of concreteness, we focused on situations where bath processes producing dark radiation are two-body decays or binary collisions, and we assumed that only one dark radiation particle is produced in the final state.The decay case allowed us to work in a model-independent fashion and solve the Boltzmann equation in momentum space for a wide range of masses and couplings.On the contrary, binary collisions require the knowledge of the scattering matrix elements and these are model-dependent.For the purpose of this study, we consider theories where matrix elements are independent of the momentum.For both cases, we produced numerical solutions for the phase space distribution of dark radiation and used them to investigate for what couplings the assumption of chemical or kinetic equilibrium is justified and to investigate the details of the decoupling process.And we used the output of the numerical integration of the Boltzmann system to evaluate ∆N eff .Crucially, we found spectral distortions in the final distribution that have an impact on the resulting ∆N eff detectable in the future.We compared these exact results with approximate momentum-space analysis employing classical statistics, and also with the output of approximated methods.Our last appendix contains material where these methods are reviewed.Figs.7 and 8 show when quantum corrections impact the final result by an amount larger than future experimental sensitivity.Likewise, Figs. 9 and 10 illustrate the error made by approximate analysis.We observe that for decays this error could be larger than the sensitivity of future CMB surveys, and therefore our analysis in momentum space is necessary if one wants accurate predictions.From the same figures, one may conclude that approximate methods based on integrated Boltzmann equations work rather well for scattering.This is indeed what we find from our numerical solutions.However, once we apply our formalism to realistic models the approximation of constant scattering matrix elements is likely to be incorrect, and consequently the collision terms would have a stronger momentum dependence that can ultimately induce spectral distortions.This is something that we defer to future work.
Our work can be extended along several other lines.For example, we could go beyond a single production and include an arbitrary number of dark radiation particles in the final state.Several models have double production as the leading process and single production is absent.As discussed in the paper, the resulting Boltzmann equation would be more complicated since the collision operator appearing on the right-hand side of the Boltzmann equation in momentum space would remain an operator and not a regular function as was the case for single production.We always considered in this study the production during a radiation-dominated epoch with the primordial bath dominating the energy budget.However, it is straightforward to extend our study to cases when the thermal bath constitutes a sub-dominant component as long as we eventually recover the RD era before the BBN epoch [50]; all it takes is using the appropriate expression for the Hubble expansion rate.
Finally, we mention a natural application of our framework to a concrete particle physics framework: light and weakly-coupled axion-like particles (ALPs).A leading candidate within this class is the QCD axion given its strong theoretical motivation from the strong CP problem.Early studies estimated the abundance of relativistic axions via the instantaneous decoupling approximation (see, e.g., Ref. [51]).Predictions for ∆N eff obtained by tracking the number density of the relativistic axion population and then converting it to a corresponding energy density (thus assuming chemical equilibrium) followed up by Refs.[52][53][54][55][56][57][58][59].Recently, axion production via pion scattering was analyzed in momentum space by Refs.[60,61] that reported spectral distortions in the final phase space distribution.This confirms what was written above about the agreement between our exact solutions and the ones derived via approximate methods since a non-trivial momentum dependence in the transition amplitudes is likely to alter how decoupling happens for different momenta.It is worth exploring spectral distortions in the axion distribution function for production via derivative couplings to SM fermions.Furthermore, we could observe the cosmological effects of a non-vanishing QCD axion mass [62][63][64][65][66][67][68].With the phase space distribution in hands, it is possible to reconsider existing analyses [69][70][71][72][73] and update the axion mass bound for production channels different than pion scatterings.We leave these developments of our general formalism to future studies.

Finally, we use the well-known relation
and after using the relation p 2 dp 2 = E 2 dE 2 we find The last integration over d cos θ gives a non-vanishing result as long as −1 ≤ cos θ ⋆ ≤ 1, and this imposes a condition on the domain for the dE 2 integral where we report also the limiting expressions for the massless X case (for which ω = k).
We conclude with the application of this general result to the scenario of thermal production via two-body decays introduced in Sec.3.1.In the main text of our paper, we write the collision term as in Eq. (3.7).The identification of the function D(k, T ) only requires a direct comparison with Eq. (A.1), and we find

A.2 Binary scatterings
The phase space integrals for scatterings pose new challenges since the matrix elements do not need to be constant.This is in contrast with the two-body decays just discussed.Nevertheless, it is possible to perform an analogous derivation that does not require information about the underlying process until the last step.We start by defining the functional ), this corresponds to Eq. (3.9) for the specific process B 1 B 2 → B 3 X (i.e., the first term on the right-hand side of Eq. (3.8)).The other contributions to the total collision rate can be obtained by permutations.As discussed before, our analysis is valid for a generic F which is a function of scalars products between momenta and moduli of spatial momenta p i , k.In particular, we do not make any assumption in this appendix about the functional dependence of the squared matrix element.We start again from the manifestly Lorentz invariant phase space measure and we integrate over dΠ 3 with the delta function that fixes P 3 = P 1 + P 2 − K.We introduce the polar coordinate system illustrated in Fig. 11 for the other integrals.In such a frame, the vector ⃗ k is fixed and we measure polar angles with respect to its direction The angular dependence of the integrand in Eq. (A.17) (and in particular of the quantity F) can only depend on the scalar products of spatial momenta given in Eq. (A.15) and therefore only on the cosine of the relative azimuthal angle ϕ.For this reason, the two different roots ϕ ⋆1 and ϕ ⋆2 will give an indentical contribution.We consider one of them and we multiply the full expression by an overall factor of 2. The delta function becomes where we introduce the following function The Heaviside theta function ensures that the argument of the square root is positive, and this condition is equivalent to imposing that the zeroes of the argument of the delta function gives cosines in Eq. (A. 19) that have absolute values equal or less than one.The final expression for the phase space integral results in The remaining four-dimensional integral can be evaluated through Monte Carlo techniques.We reduce the integration domain and make it finite via the change of variables x i = exp(−p i /Λ i ) with dp i = −Λ i /x i dx i .We choose Λ i ∼ 5k to achieve good convergence, and this leads to an acceptance rate > 50%.If F contains thermal distribution functions with temperature T an excellent choice is Λ ∼ 10T .The result is given by where V = 1 2 × 2 2 = 4 is the integration variables volume and ⟨•⟩ is the average over N samples of (x 1 , x 2 , c 1 , c 2 ).We checked that our numerical results agree with the analytical ones for the explicit cases considered in Ref. [10] at the 0.1 % level with 10 7 samples.The general expression in Eq. (A.23) is quite useful for a Monte Carlo evaluation of the collision rate.Its validity is rather general since there is no assumption about the explicit for of the function F besides the obvious fact that it is only dependent on scalar products between the four-momenta involved and the sizes of the spatial momenta.For binary collisions, and if we want to account for all the quantum effects, the specific expression is ).We conclude this Appendix with an analytical derivation of the functional in Eq. (A.23) when we employ MB statistics for bath particles.This approximation turns out to be often useful.Thus we evaluate now the analytical expression for the functional Here, we express the squared matrix element in terms of the Mandelstam variables It is convenient to start from the definition in Eq. (A.12), and rewrite the phase space measure as follows (see, e.g., Ref. [10]) The decoupling temperature for each momentum bin for decays (left panels) and scatterings (right panels).Top panels consider different statistics for the parameter space points in Fig. 2 (left) and Fig. 3 (left).Bottom panels focus on the MB statistics but with different couplings.We show the decoupling temperature (solid lines) and compare it with the one obtained via the instantaneous decoupling approximation (dot-dashed lines).
and feature results for different choices of the statistical distributions.For the decay case, the choice of the statistical distribution has no impact on the decoupling temperature meaning that attaining equilibrium will lead to the same outcome for all the statistics and differences can just arise at smaller couplings.Differences arise instead for binary scatterings.The results in the lower panels are produced for MB statistics, and compare the decoupling temperature for three different interaction strengths (solid lines) with the result obtained via the instantaneous decoupling approximation (dotdashed lines) reviewed in App.C.1.We notice how the mismatch is rather severe for decays where the instantaneous decoupling approximation underestimate typical decoupling temperatures by a factor of 2. This in turn implies that the approximate method overestimates the energy density of dark radiation (i.e., ∆N eff ).This difference is milder for scatterings.Furthermore, the functional dependence of the decoupling temperature is also different for decays and scattering.While in the former case the decoupling temperature decreases steeply for k < T , and flattens for large momenta, in the scattering case it is only weakly dependent on k/T , slightly increasing over the relevant range.
• Spectral Distortions and Kinetic Equilibrium.The different decoupling tem- T I /(AT X ) × q  : Numerical solutions for the ratio (f X /f eq X )/(R X /R eq X ) for two-body decays (top panels) and scatterings (bottom panels).We present the results at the same temperatures chosen in Fig. 4 and Fig. 5.The horizontal lines mark the constant value of 1 that is expected if thermal equilibrium is achieved.Spectral distortions are clearly visible.perature for the various momentum bins is at the hearth of the PSD spectral distortions.An effective way to visualize them is to plot (f X /f eq X )/(R X /R eq X ) as a function of normalized comoving momenta T I /(AT X ) × q at different temperatures.This is also useful to check whether or not kinetic equilibrium is a good assumption.If the ratio (f X /f eq X )/(R X /R eq X ) is constant and therefore independent on the momentum, then kinetic equilibrium is assured.If the ratio is also one, it means that we also have chemical equilibrium.We show this quantity in Fig. 13 for two-body decays (upper panels) and scatterings (lower panels) for the usual benchmark mass values adopted in this work.The different sub-plots for each panel are produced for different coupling values.For all cases, we notice non negligible deviations from kinetic equilibrium at large momenta.For small couplings, dark radiation do not ever thermalize and the deviation from kinetic equilibrium are rather large.Instead, for large interaction strengths, the fact that different momenta decouple at different times generates small spectral distortions.Comparing upper and lower panels, differences between the decay and scattering cases are clear.Spectral distortions for scatterings at strong coupling are milder and concentrated at very high comoving momenta, but three processes contribute to the production so it is easier for X particles to retain kinetic equilibrium if they thermalize.This can be observed directly in Fig. 12 which shows that the dependence of decoupling temperature on comoving momentum is much weaker than in the decay case.
• Dark Radiation Temperature.Another important quantity to track as the Universe expands is the dark radiation temperature.In particular, we want to compare it with the bath temperature T until times well after the era of electron-positron annihilations.We defined in the main text the dark radiation temperature as the width of the PSD with the appropriate normalization.An alternative definition would be to identify the temperature as the fourth root of the energy density, a procedure that clearly works well only if we are close to thermal equilibrium.We present results here for both.We define T f X as the dark radiation temperature obtained via Eq.(4.8).And we also define the temperature T ρ X satisfying the relation ρ X ∝ (T ρ X ) 4 with the correct numerical factors given in Eq. (C.3).In general, the former definition is the trustworthy one, especially in the case of small couplings where equilibrium is never achieved.We show in Fig. 14 the ratio T α X /T with α = f (solid lines) and α = ρ (dashed lines).Notice that initially the temperature ratio is zero for T ρ X since f X (A = 1) = 0, while it is finite for T f X .Notice also how this ratio evolves accordingly to the number of degrees of freedom in the thermal bath after dark radiation production stops being efficient.If thermalization is reached (as in the case for red and orange lines), the difference between the final values of T X /T arises from non-instantaneous decoupling of the different momenta.This is the reason why these differences are milder in the scattering cases.

C Approximate Methods
In this last appendix, we collect and review three approximate procedures that are alternative methods to estimate ∆N eff .First, we present an approach that relies upon the instantaneous decoupling approximation.The second method tracks the evolution of the dark radiation number density n X and it converts its freeze-out value to the corresponding energy density.This last step, which is needed to identify ∆N eff , requires the further assumption of an exact thermal distribution.We address this limitation thanks to a third improved approximate Ratio dark radiation and thermal bath temperatures vs scale factor for decays (top panels) and scatterings (bottom panels).We consider T X defined both from the width of the PSD (solid lines) and the fourth root of the energy density (dashed lines).procedure that instead keeps track of the dark radiation energy density ρ X .For this last method, we also account for the energy exchanged between the bath and the dark radiation which is not always accounted for in this kind of analysis.
Before presenting the methods, we find it useful to collect some well-known results on equilibrium thermodynamics for the case when the dark radiation PSD has a thermal shape with a temperature T X .When the dark radiation is in chemical equilibrium we have where we write the distribution in terms of the dark radiation temperature T X that is a function of the cosmic time t.We neglect the mass of X and therefore the dark radiation has the dispersion relation ω = k.In this appendix, we employ for practicality the symbol T X for the X's temperature appearing in the equilibrium PSD.This is generally different from the bath temperature T , and it is also different from T X defined by Eq. (4.8) and adopted in the main text.It is straightforward to evaluate the number density Likewise, the energy density results in Throught this review of approximate methods, we also discuss the kinetic equilibrium regime.If X particles thermalize in the early Universe, kinetic equilibrium is well justified.However, this results breaks down in the weak coupling regime.Nevertheless, the production source of dark radiation is from decays and scatterings of particles in thermal equilibrium, and the final PSD of X is aware of the underlying bath temperature.Things would be drastically different from other production mechanisms such as dark radiation from inflaton decays where final state particles are produced monochromatic (unless elastic processes can efficiently redistribute particle momenta and induce a thermal profile for the PSD).The dark radiation PSD in the kinetic equilibrium regime takes the form where µ X is the X's chemical potential.Finally, it is useful for the subsequent discussion to the number densities in the kinetic equilibrium regime where in the quantum statistical cases we have the polylogarithm function Li 3 (z).Similar relations hold for the energy densities

C.1 Instantaneous decoupling
This first approach assumes that X particles thermalize with the primordial bath at early times, and subsequently, decouple instantaneously.The temperature when such a decoupling happens, which is dubbed the freeze-out temperature T D , can be estimated through a comparison between two competing effects.On the one hand, collisions favor thermalization.On the other hand, the Hubble expansion dilutes and cools the universe down making it hard for particles to stay in thermal contact.We identify the freeze-out temperature by imposing that at the decoupling moment the collision rate is equal to the Hubble expansion rate The dark radiation energy density at the decoupling time is still the equilibrium one, whereas X particles free-stream afterward and their energy density decreases with the scale factor as ρ X ∝ a −4 .We define the interaction rate as follows3 The above expression is general and valid for any choice of the statistical distributions.The bath particles PSDs are the equilibrium ones in Eq. (2.4), and consistently with our hypothesis we set the PSD for the X to the equilibrium expression in Eq. (C.1).We normalize Γ X with a factor proportional to the inverse equilibrium number density in Eq. (C.2).The Hubble rate is given by the Friedmann equation, The visible contribution to the energy density reads ρ B = (π 2 /30)g B ⋆ρ (T )T 4 whereas the invisible one before decoupling is given by the expression in Eq. (C.3).
Before the decoupling epoch, all the thermodynamic quantities for the dark radiation are the equilibrium ones with T X = T .After decoupling, the use of comoving quantities ease the analysis.The entropy density, defined as s(T ) = (2π 2 /45)g ⋆s (T )T 3 , turns out to be quite useful.The effective number of entropic degrees of freedom is the sum of two contributions, g ⋆s (T ) = g B ⋆s (T )+g X ⋆s (T ).The former is a property of the thermal bath, the latter follows from the fact that f X remains thermal after decoupling but with T X ̸ = T because bath particles become non-relativistic through the expansion, and it explicitly reads We define the dimensionless ratio R X (T ) ≡ ρ X (T )/s(T ) 4/3 .Entropy conservation (s ∝ a −3 ) imposes R X (T < T D ) = R X (T D ), and in particular we have R X (T CMB ) = R X (T D ).
Finally, we evaluate ∆N eff via Eq.(2.11) and we find The result is proportional to the ratio, with the exponent of 4/3, of the total number (i.e., contributions from both the bath and X) of entropic degrees of freedom at the recombination and decoupling epochs.Once we plug the explicit temperature dependence of the X's contribution, as given in Eq. (C.9), we find

C.2 Tracking the number density
Another approximate method commonly employed in the literature is to track the dark radiation number density defined from the PSD as follows n X = g X d 3 k (2π) 3 f X (k, t) . (C.12) We derive the Boltzmann equation describing the time evolution of this quantity.Our starting point is the Boltzmann equation in momentum space in Eq. (2.2), and we multiply both sides by (g X /ω 1 ).For the collision operator, we take the general expression in Eq. (2.6).After we integrate both sides over the three components of the spatial momentum ⃗ k 1 , we find e ωr/T f X (k r ) 1 ± f X (k r ) . (C.13) The above result is still not an ordinary differential equation, and we need some further assumptions to achieve a simpler equation to handle.Given the thermal origin of the produced X particles, it is reasonable to assume kinetic equilibrium and take the PSD given in Eq. (C.4).
If we have f X = f k.eq X then the following relation holds This result is valid for any statistical distribution of the dark radiation degree of freedom, and it allows us to take the term inside the square brackets in Eq. (C.13) outside of the integral since it is independent on the momenta.We find the following Boltzmann equation (1 ± f X (k r )) . (C.15) The term inside the square bracket can be related to the X's number density.Consistently with our assumption of kinetic equilibrium, namely that the PSD distribution has the functional form given in Eq. (C.4), the explicit expression for n X = n k.eq is provided by Eq. (C.5).However, even after we assume kinetic equilibrium the resulting Eq. (C.15) is still not an ordinary differential equation for n X because of the terms containing f X on the right-hand side.Our only option is to approximate 1 ± f X ≈ 1, and this is equivalent to neglect quantum degeneracy effects for X.From now on, we employ MB statistics for the dark radiation.
To summarize, once we assume kinetic equilibrium as well as MB for the X's, we obtain the Boltzmann equation describing the number density evolution dn X dt + 3Hn X = γ X 1 − n X n eq X ℓ . (C.16) Here, the temperature-dependent function γ X accounts for the averaged number of X particles produced per unit of time and volume, and it explicitly reads (1 ± f B j ) .(C.17) This general expression for the production rate accounts for the quantum degeneracy effects of bath particles.This has to be compared with the expression in Eq. (C.8) that accounts also for the quantum degeneracy effects of the dark radiation particle.This is possible because within the instantaneous decoupling approximation we assume that X thermalizes, and we can compute the phase space integral by using its equilibrium distribution.Our next step is to solve Eq. (C.16) and evaluate the X's abundance at the recombination time.As it is often done in the literature, it is convenient to switch to the comoving number density Y X = n X /s to factor out the cosmological dilution due to the Hubble expansion, and to employ the thermal bath temperature as the evolution variable.Upon imposing entropy conservation (i.e., d(sa 3 )/dt = 0), we find the following form of the Boltzmann equation for the comoving number density (C.18) Here, Y eq X = n eq X /s is the comoving equilibrium number density.It is enough to solve the Boltzmann equation up to the point when production stops being effective (e.g., couplings get too small and/or bath particles participating in the production become non-relativistic and disappear from the bath) and Y X reaches a constant asymptotic value.
We still do not know how to account for the contribution of X particles to the energy and entropy densities.The former affects the Boltzmann equation once we express the Hubble parameter H via the Friedmann equation whereas the latter appear in the Boltzmann equation both explicitly and implicitly via the functions s and Y eq X .Typically, this extra contribution is neglected or it is considered to be maximal g * α (T ) ≃ g B * α (T ) Neglecting X g B * α (T ) + g X ξ ρ X Thermalized X (α = ρ, s) .

(C.19)
Under this highly common assumption, where one adopts one of the two extreme choices, the equation for Y X can be integrated and we can find its asymptotic value.The last step is converting the asymptotic comoving number density into an asymptotic comoving energy density.A standard assumption is that the X's PSD has the thermal shape in Eq. (C.1) with a generic temperature T X not necessarily equal to the one for the bath.We emphasize how this is a further assumption with respect to the ones already made in this sub-section.The Boltzmann equation in Eq. (C.16) has been derived under the assumption of kinetic equilibrium (as well as MB) for the dark radiation, whereas the distributions in Eq. (C.1) are valid if chemical equilibrium holds.Assuming kinetic equilibrium only would not be sufficient to make this conversion because of the unknown chemical potential µ X entering the expressions in Eqs.(C.5) and (C.6).Thus the only way to proceed is to make the further assumption of chemical equilibrium for this conversion, with number and energy densities given in Eqs.(C.2) and (C.3), and this allows us to find the relation (C.20) While it is true that we assumed MB statistics for X to derive the Boltzmann equation, this relation can be specialized to different cases for the dark radiation statistical properties upon choosing the appropriate ξ factors.We find that using the appropriate statistics instead of the MB improves the error with respect to the full solution.
Finally, we determine the contribution to ∆N eff as defined in Eq. (2.11)The final result for ∆N eff is expressed in terms of the asymptotic solution of the Boltzmann equation Y X (T CMB ).We include only the entropic degrees of freedom from the bath g B * s (T CMB ).The analogous contribution from dark radiation g X * s (T CMB ) is between the two extremes in Eq. (C.19), and it is likely to be suppressed with respect to the factor of g X ξ ρ X because threshold effects in the thermal bath between decoupling and recombination that lead to T X < T (see Eq. (C.9)).Thus taking the "Thermalized X" case is valid to solve the Boltzmann equation but it would not be appropriate at recombination.
We conclude this sub-section with an improvement that includes also the X's contribution to the entropy and energy densities quantified as follows (see Eq. (C.9)) Here, we still assume that the X's PSD is a chemical equilibrium one with an arbitrary temperature T X .We can invert the first relation to find Once we know how the full contribution g * s (T ) depends on the temperature, we can plug it into the expression for g X * ρ (T ) above and account also for the dark radiation contribution to the energy density.We use these improved estimates for the total number of effective degrees of freedom into the Boltzmann equation, and we can find a new differential equation for Y X (T ).Finally, after we solve it and find the asymptotic value of Y X , we can compute

C.3 Tracking the energy density
As we have just seen the number density is not the most suitable quantity to compute ∆N eff .
Here, we derive a Boltzmann equation describing the evolution of the energy density ρ X .We proceed very similarly as in the previous case, up to Eq. (C.15), which in this case reads

Figure 3 :
Figure3: Collision terms for binary scatterings divided by the Hubble parameter as a function of the bath temperature T .The particle B 3 is massless, whereas B 1 and B 2 have the same mass m equal to 1 GeV (left) and 1 TeV (right).We fix the scattering matrix element to |M| 2 = 10 −11 .The notation is the same as the one adopted in Fig.2.

Figure 4 :
Figure4: PSDs as a function of the comoving momentum for production via two-body decays for m 1 = 1 GeV (left) and m 1 = 1 TeV (right).Different rows are for Γ 1 /m 1 = 10 −18 (top, blue), 10 −14 (middle, orange) and 10 −10 (bottom, red).We show the PSD at different scale factor values A T when the bath temperature was T = (100, 1, 0.01) × m 1 (dot-dashed, dashed, solid lines).We compare them with a thermal distribution with temperature T X defined in Eq. (4.8) (dotted black lines).The choices for the normalizations of both quantities on the Cartesian axes are explained in the text.

Figure 7 :
Figure7: Predicted ∆N eff for dark radiation produced via two-body decays as a function of the coupling strength Γ 1 /m 1 for m 1 = 1 GeV (left) and m 1 = 1 TeV (right).We consider bosonic (upper panels) and fermionic (lower panels) dark radiation.Solid lines show the exact results where we employ quantum statistical distributions for all particles, and we also report results from approximate analysis: all quantum statistical effects neglected (thick dotted black lines), and only X's quantum statistics accounted for (thick dot-dashed black lines).In the lower part of each panel we quantify the errors made by the approximate results.

Figure 8 :
Figure 8: Predicted ∆N eff for production via scatterings as a function of the interaction strength g 1 g 2 g 3 |M| 2 for m = 1 GeV (left) and m = 1 TeV (right).Notation as in Fig. 7.

Figure 10 :
Figure10: Comparison of ∆N eff predictions for binary scatterings in the MB approximation as a function of the squared matrix element |M| 2 .We pick the same mass spectrum adopted in the paper with a single overall mass scale m, and we present results for m = 1 GeV (left panel) and m = 1 TeV (right panel).The notation is the same as in Fig.9.

B 1 → B 2 X, m 1 =BB 1 → 12 Γ1B 2 = 10 − 12 |M| 2 = 10 − 10 |M| 2 =Figure 12 :
Figure12: The decoupling temperature for each momentum bin for decays (left panels) and scatterings (right panels).Top panels consider different statistics for the parameter space points in Fig.2(left) and Fig.3(left).Bottom panels focus on the MB statistics but with different couplings.We show the decoupling temperature (solid lines) and compare it with the one obtained via the instantaneous decoupling approximation (dot-dashed lines).

Figure 13
Figure13: Numerical solutions for the ratio (f X /f eq X )/(R X /R eq X ) for two-body decays (top panels) and scatterings (bottom panels).We present the results at the same temperatures chosen in Fig.4and Fig.5.The horizontal lines mark the constant value of 1 that is expected if thermal equilibrium is achieved.Spectral distortions are clearly visible.

Table 1 :
Explicit expressions for the function D(k, T ) appearing in in Eq. (3.7) for different bath particles statistics with E