Dark Sector Equilibration During Nucleosynthesis

TL;DR

This paper explores how dark sectors equilibrating during nucleosynthesis can alter the universe's expansion rate, affecting light element abundances and the effective number of neutrino species, with implications for cosmological models.

Contribution

It provides a model-independent framework to analyze dark sector equilibration effects on nucleosynthesis and examines specific models where dark sectors interact with neutrinos or photons.

Findings

Significant dark sector effects can be consistent with observed light element abundances.

Timing of dark sector equilibration critically influences nucleosynthesis outcomes.

The impact on $N_{\mathrm{eff}}$ varies depending on the model and epoch of equilibration.

Abstract

Light, weakly-coupled dark sectors may be naturally decoupled in the early universe and enter equilibrium with the Standard Model bath during the epoch of primordial nucleosynthesis. The equilibration and eventual decoupling of dark sector states modifies the expansion rate of the universe, which alters the predicted abundances of the light elements. This effect can be encompassed in a time-varying contribution to $N_{\mathrm{eff}}$, the effective number of neutrino species, such that $N_{\mathrm{eff}}$ during nucleosynthesis differs from its measured value at the time of recombination. We investigate the impact of such variations on the light element abundances with model-independent templates for the time-dependence of $N_{\mathrm{eff}}$ as well as in specific models where a dark sector equilibrates with neutrinos or photons. We find that significant modifications of the expansion rate are consistent with the measured abundances of light nuclei, provided that they occur during specific periods of nucleosynthesis. In constraining concrete models, the relative importance of the cosmic microwave background and primordial nucleosynthesis is highly model-dependent.

Tables

Table 1. Table 1: Notation and various temperature scales discussed throughout this work.

Notation	Definition
$T_i$	temperature of species $i=\mathrm{ds},\nu ,\gamma$
$T$	shorthand for the photon temperature ($T_{\gamma }$)
${\xi }_i$	temperature of species $i$ normalized to the photon temperature
${\xi }_{\mathrm{ds}}^0$	value of ${\xi }_{\mathrm{ds}}$ before DS-SM equilibration
${\xi }_{\nu }^{\text{SM}}$	value of ${\xi }_{\nu }$ in a standard (unmodified) cosmology
$T^{\text{ds eq}}$	temperature of the photon bath at DS-SM equilibration

Table 2. Table 2: Nuclear rates updated relative to AlterBBN v1.4  [ 49 ] . In the text, we compare the predictions using NACRE-II  [ 55 ] and Coc et al. (2015)  [ 58 ] for the reactions in the last three rows due to their importance in estimating the deuterium yield.

Rate	Reference
$p(n,\gamma )\text{D}$	[56]
${}_{}{}^3\mathrm{He}(n,p){}_{}{}^{3}H$	[57]
$\text{D}(\text{D},\gamma ){}_{}{}^4\mathrm{He}$	[55]
${}_{}{}^{3}H(\text{D},n){}_{}{}^4\mathrm{He}$
${}_{}{}^3\mathrm{He}(\text{D},p){}_{}{}^4\mathrm{He}$
$\text{D}({}_{}{}^4\mathrm{He},\gamma ){}_{}{}^6\mathrm{Li}$
${}_{}{}^6\mathrm{Li}(p,\gamma ){}_{}{}^7\mathrm{Be}$
${}_{}{}^7\mathrm{Li}(p,\gamma ){}_{}{}^8\mathrm{Be}$
${}_{}{}^6\mathrm{Li}(p,{}_{}{}^4\mathrm{He}){}_{}{}^3\mathrm{He}$
${}_{}{}^7\mathrm{Li}(p,{}_{}{}^4\mathrm{He}){}_{}{}^4\mathrm{He}$
$\text{D}(\text{D},n){}_{}{}^3\mathrm{He}$	[55, 58]
$\text{D}(p,\gamma ){}_{}{}^3\mathrm{He}$
$\text{D}(\text{D},p){}_{}{}^{3}H$

Equations

$${\Gamma }_{\text{eq}}\sim {\alpha }_{\mathrm{ds}}^{2}T(\text{DM-SM scattering})$$

$${\Gamma }_{\text{eq}}\sim {\alpha }_{\mathrm{ds}}{m}_{\mathrm{ds}}^2/T(\text{DS decays to SM})$$

$$H\sim T^2/M_{\mathrm{Pl}},$$

$$T^{\mathrm{ds}\text{ eq}}$$

$$\mathcal{L}\supset \lambda \phi {\nu }^2+h.c.,$$

$$\mathcal{L}\supset \frac{\phi }{4\Lambda }{\tilde{F}}_{\mu \nu }{F}^{\mu \nu },$$

$$\rho \approx {\rho }_{\gamma }+{\rho }_{e^{\pm }}+{\rho }_{\nu }+{\rho }_{\mathrm{ds}}.$$

$${\rho }_{\nu }+{\rho }_{\mathrm{ds}}=\frac{7}{8}\left(\frac{{\pi }^2}{15}\right)\times N_{\mathrm{eff}}(T)\times (T_{\nu }^{\text{SM}}{)}^4.$$

$${\xi }_i\equiv \frac{T_i}{T}.$$

$$N_{\mathrm{eff}}=3\bigg{[}\left(\frac{\xi_{\nu}}{\xi_{\nu}^{\text{SM}}}\right)^{4}+\frac{g_{*}^{\mathrm{ds}}}{g_{*}^{\nu}}\left(\frac{\xi_{\mathrm{ds}}}{\xi_{\nu}^{\text{SM}}}\right)^{4}\bigg{]}\,,$$

$$({\rho }_{\mathrm{SM}}+{\rho }_{\mathrm{ds}})a^4=\mathrm{constant}\text{(SM–DS equilibration)},$$

$$(s_{\mathrm{SM}}+s_{\mathrm{ds}})a^3=\mathrm{constant}\text{(SM–DS decoupling)},$$

$${\dot{\rho }}_i+3H({\rho }_i+p_i)=g_i\int \frac{d^{3}p}{(2\pi {)}^3}C[f_j],$$

$$g_i\int \frac{d^{3}p}{(2\pi {)}^3}C[f_j]\simeq -m_{\phi }{\Gamma }_{\phi }\left[n_{\phi }(T_{\mathrm{ds}})-n_{\phi }^{\mathrm{eq}}(T_{\mathrm{SM}})\right],$$

$$N_{\mathrm{eff}}\approx 3+\frac{4}{7}{g}_{\ast }^{\mathrm{ds}}{\left(\frac{{\xi }_{\mathrm{ds}}^0}{{\xi }_{\nu }^{\mathrm{SM}}}\right)}^4,$$

$$N_{\mathrm{eff}}\approx 3{\left(1+\frac{g_{\ast }^{\mathrm{ds}}}{g_{\ast }^{\nu }}\right)}^{1/3}\left(1+\frac{g_{\ast }^{\mathrm{ds}}}{g_{\ast }^{\nu }}({\xi }_{\mathrm{ds}}^0{)}^4\right).$$

$$N_{\mathrm{eff}}\approx 3+\frac{4}{7}{g}_{\ast }^{\mathrm{ds}}{\left(\frac{{\xi }_{\mathrm{ds}}^0}{{\xi }_{\nu }^{\mathrm{SM}}}\right)}^4,$$

$$N_{\mathrm{eff}}\approx 3\left(\frac{g_{\ast }^{\gamma }+g_{\ast }^{\mathrm{ds}}}{g_{\ast }^{\gamma }+g_{\ast }^{\mathrm{ds}}({\xi }_{\mathrm{ds}}^0{)}^4}+\frac{g_{\ast }^{\mathrm{ds}}/g_{\ast }^{\nu }}{({\xi }_{\nu }^{\mathrm{SM}}{)}^4}\right).$$

$$g_{\mathrm{critical}}^{\mathrm{ds}}\approx\frac{1}{(\xi_{\mathrm{ds}}^{0})^{4}}~{}\left[g_{*}^{\nu}\,(\xi_{\nu}^{\text{SM}})^{4}\,\Big{(}(\xi_{\mathrm{ds}}^{0})^{4}-1\Big{)}-g_{*}^{\gamma}\right]~{}.$$

$$N_{\mathrm{eff}}\approx 3\left[\Big{(}1+\frac{g_{*}^{\mathrm{ds}}}{g_{*}^{\gamma}}\Big{)}^{1/3}\Big{(}1+\frac{g_{*}^{\mathrm{ds}}}{g_{*}^{\gamma}}\,(\xi_{\mathrm{ds}}^{0})^{4}\Big{)}\right]^{-1}\,,$$

$${\eta }_b={\left({\xi }_{\nu }\right)}^3{\eta }_b^0.$$

$${\eta }_b^0={\xi }_{\nu }(T_{\mathrm{CMB}}{)}^{-3}{\eta }_b^{\mathrm{CMB}},$$

$$\xi_{\nu}\approx\xi_{\nu}^{\text{SM}}\times\begin{cases}1&(T>T^{\mathrm{ds}\text{ eq}})\\ \left(\frac{g_{*}^{\gamma}+g_{*}^{\mathrm{ds}}}{g_{*}^{\gamma}+g_{*}^{\mathrm{ds}}\,(\xi_{\mathrm{ds}}^{0})^{4}}\right)^{1/4}&(m_{\mathrm{ds}}\lesssim T<T^{\mathrm{ds}\text{ eq}})\\ \frac{\big{(}g_{*}^{\gamma}\big{)}^{1/3}}{\big{(}g_{*}^{\gamma}+g_{*}^{\mathrm{ds}}\big{)}^{1/12}\big{(}g_{*}^{\gamma}+g_{*}^{\mathrm{ds}}\,(\xi_{\mathrm{ds}}^{0})^{4}\big{)}^{1/4}}&(T\lesssim m_{\mathrm{ds}})\,.\end{cases}$$

$$Y_p\approx \frac{4n_{\text{He}}}{n_p+n_n}\approx \frac{4(n_n/2)}{n_p+n_n}=\frac{2(n/p)}{1+(n/p)}\approx \frac{1}{4}.$$

$$\frac{n_{\mathrm{D}}}{n_p+n_n}\sim {\eta }_b{\left(\frac{T}{m_p}\right)}^{3/2}{e}^{B_{\mathrm{D}}/T},$$

$$T_{\text{BBN}}\sim B_{\mathrm{D}}/\log (1/{\eta }_b)\sim \mathcal{O}(100)\text{ keV}\times \left(1+\mathcal{O}(10^{-2})\frac{\delta {\eta }_b}{{\eta }_b}\right),$$

$$(Y_p{)}^{\text{obs.}}=0.245\pm 0.003,(\text{D}/\text{H}{)}^{\text{obs.}}=(2.569\pm 0.027)\times 10^{-5}.$$

$$Y_p=0.24633\pm 0.00034,10^{5}D/H=2.57\pm 0.12(10^{10}{\eta }_b^{\text{CMB}}=6.10)$$

$$Y_p=0.24633\pm 0.00034,10^{5}D/H=2.57\pm 0.13(10^{10}{\eta }_b^{\text{CMB}}=6.10\pm 0.04),$$

$$Y_p=0.24642\pm 0.00032,10^{5}D/H=2.446\pm 0.037(10^{10}{\eta }_b^{\text{CMB}}=6.10)$$

Full text

1. 1 Introduction

2. 2 Late Dark Sector Equilibration

3. 2.1 Motivation

4. 2.2 Impact on the Expansion Rate and the Baryon Density

5. 2.2.1 Coupling to the Neutrino Bath

6. 2.2.2 Coupling to the Photon Bath

7. 3 Nucleosynthesis

8. 3.1 Dependence on $N_{\mathrm{eff}}$ and $\eta_{b}$

9. 3.2 Procedure

10. 3.3 Standard Constraints on $N_{\mathrm{eff}}$ and $\eta_{b}$

11. 4 Equilibration and Decoupling during Nucleosynthesis

12. 4.1 Model-Independent Results

13. 4.2 Model-Specific Results

14. 5 Conclusion

15. A Collision Terms

16. B Constraints on Dark Sector-Photon Equilibration

1 Introduction

Measurements of the light element abundances provide one of the earliest direct tests of cosmology. Broad agreement of the observed abundances with predictions of standard Big Bang nucleosynthesis (BBN) has been used to constrain scenarios with late reheating [1, 2], anti-matter domains [3], long-lived particle decays [4, 5], invisible [6] and electromagnetic energy injection [7, 8, 9], and new light particles in thermal equilibrium with the Standard Model (SM) [10, 11, 12, 13]. For a survey of phenomena that can be tested by nucleosynthesis see, e.g., Ref. [14]. Given the far-reaching implications of the observed primordial abundances of light nuclei, it is important to understand the robustness of the resulting conclusions.

In this work, we focus on possible modifications to the expansion rate at the time of primordial nucleosynthesis. Simple deformations of standard assumptions can leave open significant room for non-standard phenomena [15]. One key assumption that is often made when constraining light particles is that the new degrees of freedom were already in equilibrium with the SM at the onset of nucleosynthesis. Light, weakly coupled dark sectors are typically not in thermal equilibrium with the SM until late times [16, 17, 18, 19, 20, 21]. For instance, if equilibration takes place after neutrino-photon decoupling, the resulting modification to the expansion rate is suppressed [17, 18, 20, 21], allowing for the presence of many new degrees of freedom at the SM temperature. Equilibration during nucleosynthesis is therefore a natural possibility.

The scenarios we consider are schematically shown in Fig. 1. At the onset of nucleosynthesis, the SM bath consists of two decoupled components – photons and neutrinos. The equilibration of a light dark sector (DS) with one of these two SM components has a distinct impact on the effective number of neutrino species, $N_{\mathrm{eff}}$, which parametrizes the energy density of the primordial plasma in non-electromagnetic degrees of freedom. New particles coupled to neutrinos (or those that are completely decoupled) increase $N_{\mathrm{eff}}$ at late times, while those that interact with photons decrease $N_{\mathrm{eff}}$. Dark sector equilibration and decoupling with the SM can give rise to “pulses” in the modified expansion rate (and therefore $N_{\mathrm{eff}}$) that occur during a specific stage of nucleosynthesis; this is illustrated in Fig. 1. Additionally, new species coupled to photons can dilute the baryon-to-photon ratio, $\eta_{b}$, which encodes the baryon asymmetry of the universe.

We discuss the conditions under which an initially decoupled DS attains equilibrium with the SM during nucleosynthesis in Sec. 2. In Sec. 3, we summarize the sensitivity of the predictions of BBN to changes in the expansion rate and the baryon asymmetry. There, we also discuss our numerical methods and treatment of nuclear rate uncertainties. We use this framework to first study the impact of time-dependent modifications to $N_{\mathrm{eff}}$ in a model-independent manner, by considering generic templates for $N_{\mathrm{eff}}$ variations that can be mapped onto specific models in Sec. 4.1. We then analyze a concrete scenario where the leading interactions between the DS and SM are mediated by decays and inverse-decays of a light mediator in Sec. 4.2. In this case, we obtain a realistic time-evolution for $N_{\mathrm{eff}}$ by solving the relevant Boltzmann equations for the DS and SM plasmas (with the relevant collision terms given in Appendix A). We summarize our findings in Sec. 5. A summary of common notation and definitions used throughout this work is provided in Table 1.

2 Late Dark Sector Equilibration

2.1 Motivation

We focus on models in which a DS relativistically enters thermal equilibrium with the SM bath after neutrino-photon decoupling, which occurs at temperatures of $T\sim\text{few}\times\mathrm{MeV}$ [22, 23]. We will refer to this process as late equilibration.111One can also call this late recoupling, in analogy to recombination; the prefix here does not necessarily imply that the DS and the SM were thermally coupled at earlier times. These models are interesting for several reasons. First, this cosmology is a generic feature of dark sectors with light (sub-MeV) feebly-interacting particles, as we argue below. Second, strong constraints on thermalized sub-MeV particles derived from considerations of the cosmic microwave background (CMB) and BBN are significantly relaxed in these scenarios, as we will show in Sec. 2.2. Finally, late equilibration provides an explicit example of a time-dependent modification of the Hubble expansion rate between nucleosynthesis and recombination, thereby highlighting the complementarity of BBN and CMB measurements in constraining the evolution of the early universe.

Late equilibration naturally occurs if a sub-MeV DS interacts with the SM through a light weakly-coupled force carrier. To see this, let us denote the rate for thermal equilibration between the DS and SM as $\Gamma_{\text{eq}}$. At temperatures much greater than the masses of various particles participating in the equilibrating reactions, the rate is schematically of the form

[TABLE]

for DM-SM scattering and

[TABLE]

for decays and inverse-decays of DS force-carriers into the SM, where $T$ is the photon bath temperature, $\alpha_{\mathrm{ds}}$ is a coupling constant, and $m_{\mathrm{ds}}$ is a characteristic mass-scale in the DS. Equilibration occurs when $\Gamma_{\text{eq}}$ is comparable to the Hubble expansion rate,

[TABLE]

where $M_{\mathrm{Pl}}$ is the Planck mass. Since $\Gamma_{\text{eq}}/H$ increases as the universe cools, the DS will eventually enter equilibrium with the SM, provided that $\alpha_{\mathrm{ds}}$ is sufficiently large. Defining $T^{\mathrm{ds}\text{ eq}}$ to be the temperature at which $\Gamma_{\text{eq}}\sim H$, we find that the DS and SM enter equilibrium when the temperature drops below

[TABLE]

for scattering and decays, respectively. This occurs after neutrino-photon decoupling if the couplings are sufficiently small. For instance, decays and inverse-decays equilibrate the two sectors after neutrino-photon decoupling but before the end of nucleosynthesis if the DS coupling falls in the wide range, $\mathcal{O}(10^{-21})\lesssim\alpha_{\mathrm{ds}}\,(m_{\mathrm{ds}}/\text{keV})^{2}\lesssim\mathcal{O}(10^{-15})$.

Once equilibrated, the DS will remain thermally coupled to the SM bath until the temperature drops below some characteristic mass scale (typically the mass of the heaviest particle involved in the reaction), at which point $\Gamma_{\text{eq}}$ becomes mass- or Boltzmann-suppressed and the DS decouples from the SM bath. As a result, each time the temperature drops below a DS mass-threshold, the corresponding species becomes non-relativistic and its comoving entropy density is transferred to the SM, similar to $e^{+}e^{-}$ annihilations in the early universe. Throughout this work, we will refer to such processes as DS-SM decoupling.

In this work, we consider two concrete examples in which a DS couples to neutrinos or photons. These scenarios affect the expansion rate of the universe in distinct ways. Late equilibration between a DS and the neutrino bath has been studied in Refs. [17, 18] and was recently highlighted in Refs. [20, 21] as a way of realizing cosmologies of thermal dark matter below the MeV scale. A generic feature of these models is a light mediator coupled to neutrinos. Such particles are a natural feature of models of neutrino mass generation through spontaneous lepton-number ($L$) violation. The pseudo-Nambu-Goldstone boson of $L$-breaking, the majoron, is light and has renormalizable interactions with the neutrino mass-eigenstates. Moreover, the majoron can be coupled to additional light states that, for example, account for dark matter [20, 21].

Thus, the existence of a majoron can lead to the late equilibration of a DS with neutrinos. In practice, we will consider a simplified model consisting of a light mediator, $\varphi$, interacting with neutrinos,

[TABLE]

where $\nu$ is a neutrino mass-eigenstate in two-component notation. This coupling arises in a gauge-invariant manner after electroweak symmetry breaking mixes $SU(2)_{L}$-charged neutrinos with singlet right-handed neutrinos (with $\varphi$ directly coupled to the latter) [21]. The dominant process leading to energy transfer and equilibration between the SM and the DS is the decay $\varphi\leftrightarrow\nu\nu$.

Late DS equilibration with photons was first considered in Ref. [16] in a model with a millicharged particle and a dark photon. In that case, the dominant equilibrating reaction is a Compton-like scattering process. For simplicity and in analogy to the neutrino case above, we will instead focus on a model where the dominant process is decay; this is realized by an axion-like particle, $\varphi$, with an interaction [19]

[TABLE]

where $F$ is the photon field strength and $\Lambda$ encodes the scale and couplings of ultraviolet physics that generates this operator. As we show in Appendix B, in minimal models the size of $\Lambda$ that is required for equilibration to occur shortly before or during nucleosynthesis is typically ruled out by stellar cooling arguments and the observed SN1987A burst duration. These bounds can potentially be avoided in more baroque scenarios [21]. For simplicity, we will assume that the resolution of these issues has no impact on cosmology, especially on nucleosynthesis.

Since we are interested in the cosmological impact of late equilibration, we only use three parameters ($m_{\text{ds}},g_{\text{ds}},\Gamma_{\text{eq}}$ corresponding to the DS mass, DS internal degrees of freedom, and DS-SM equilibrating rate, respectively) to specify a model in which a DS couples to either neutrinos or photons. In the next section, we will discuss in detail how this may lead to non-standard modifications to the expansion rate of the universe.

2.2 Impact on the Expansion Rate and the Baryon Density

Light DS states can modify the history of the universe in several ways. The main effect considered in this work is on the Hubble expansion rate, which is determined by the total energy density, $H\propto\sqrt{\rho}$. As $\rho$ is dominated by relativistic species, after neutrino-photon decoupling it receives contributions from photons, electrons, neutrinos, and potentially the DS:

[TABLE]

The energy density of each relativistic species depends on its internal degrees of freedom, $g_{i}$, and its temperature, $T_{i}$. We define the effective relativistic degrees of freedom, $g_{*}^{i}$, such that $g_{*}^{i}=g_{i}$ for bosons and $g_{*}^{i}=(7/8)\,g_{i}$ for fermions. The energy density of each relativistic species is then given by $\rho_{i}=(\pi^{2}/30)\,g_{*}^{i}\,T_{i}^{4}$.

It is convenient to parametrize the contribution of a light DS to the energy density in terms of the effective number of neutrino species, $N_{\mathrm{eff}}$, which is defined as

[TABLE]

Above, $T_{\nu}^{\text{SM}}$ is the neutrino temperature at a given time (or, equivalently, at a given photon temperature) in the standard cosmology and in the instantaneous-decoupling approximation. Under this approximation, $N_{\mathrm{eff}}=3$ in the SM.222Incomplete neutrino decoupling at the time of $e^{+}e^{-}$ annihilation leads to neutrino spectral distortions and $N_{\mathrm{eff}}\approx 3.045$ [24, 25, 26]. This subtlety is not accounted for in our treatment of early universe thermodynamics and nucleosynthesis, so we approximate $N_{\mathrm{eff}}\approx 3$ in the SM. In scenarios with extra light degrees of freedom at the same temperature as the SM neutrinos, $N_{\mathrm{eff}}$ simply counts these particles. Below, we consider models in which the DS has a different initial temperature from that of the neutrinos, and the neutrino temperature evolution is modified; in this case, $N_{\mathrm{eff}}$ no longer directly measures the number of particles in the DS. Equation 8 is consistent with the more common definition, $N_{\mathrm{eff}}\propto(\rho_{\nu}+\rho_{\mathrm{ds}})/\rho_{\gamma}$, often used in the context of the CMB. Our adopted definition of $N_{\mathrm{eff}}$ also accounts for non-standard neutrino temperature evolution.

Equilibration or decoupling of any species after neutrino-photon decoupling typically results in time-evolution of $N_{\mathrm{eff}}$. For example, in the SM, $e^{\pm}$ decoupling increases $N_{\mathrm{eff}}$ to $\approx 3.045$ [24, 25, 26]. Throughout this work, we use the photon temperature, $T\equiv T_{\gamma}$, as a standard reference for cosmological epochs. For later convenience, we define $\xi_{i}$ as the relative temperature of species $i$ compared to that of the photon bath [27],

[TABLE]

We use $\xi_{\nu}^{\mathrm{SM}}$ to denote the value of $\xi_{\nu}$ under the assumptions of a standard cosmology and instantaneous decoupling, i.e., $\xi_{\nu}^{\mathrm{SM}}(T)=1$ for $T\gg m_{e}$ and $\xi_{\nu}^{\mathrm{SM}}(T)=(4/11)^{1/3}\simeq 0.7$ for $T\ll m_{e}$ [28]. Equation 8 can be converted into an explicit expression for $N_{\mathrm{eff}}$,

[TABLE]

where $g_{*}^{\nu}=3\times 2\times(7/8)=21/4$ and $g_{*}^{\mathrm{ds}}=g_{*}^{\mathrm{ds}}(T)$ is the (temperature-dependent) effective number of relativistic degrees of freedom in the DS. In general, $\xi_{\nu}$, $\xi_{\nu}^{\text{SM}}$, $\xi_{\mathrm{ds}}$, and $g_{*}^{\mathrm{ds}}$ all evolve with time (or temperature) as the universe expands and cools.

Dark sector equilibration is governed by the conservation of comoving energy. Consider, for example, a DS that equilibrates with neutrinos after they have themselves decoupled from photons. The energy of the closed DS+$\nu$ system is conserved up to the work done by the bath in driving the expansion of the universe. As long as the participating DS and SM species are relativistic, energy conservation in an expanding universe takes the form

[TABLE]

where $a$ is the scale factor. The temperature evolution of the $\mathrm{SM}+\text{DS}$ bath can be estimated by equating the value of the left-hand side of Eq. 11 immediately before and after equilibration (assuming equilibration is instantaneous). As this process is irreversible, entropy is not conserved.

As the universe continues to expand and cool, the temperature eventually drops below the masses of various DS degrees of freedom. As these DS particles become non-relativistic, they deposit their entropy into the $\mathrm{SM}+\text{DS}$ bath through, e.g., decays or annihilations and eventually decouple. This process occurs in equilibrium and therefore conserves entropy, leading to

[TABLE]

which can be used to estimate the $\mathrm{SM}+\text{DS}$ temperature before and after decoupling. Equations 11 and 12 will be used to qualitatively understand the temperature evolution of $N_{\mathrm{eff}}$ for the scenarios considered below in Secs. 2.2.1 and 2.2.2.

Equilibration and decoupling are not instantaneous processes. The quantitative evolution of energy and entropy densities in the various sectors can be obtained from Boltzmann equations such as

[TABLE]

where $i=\gamma,\nu,\mathrm{ds}$ and $C$ is the appropriate collision term. As discussed in Sec. 2.1, late equilibration of a feebly-interacting DS naturally occurs, for example, when the dominant process mediating DS-SM interactions is the decay of a light mediator, $\varphi$, into SM particles ($\varphi\to\mathrm{SM}$) with the decay rate, $\Gamma_{\varphi}$. The corresponding collision term is then given approximately by

[TABLE]

where the superscript “eq” denotes an equilibrium distribution with $n_{\varphi}^{\mathrm{eq}}(T_{i})\sim T_{i}^{3}$ in the relativistic limit. We will also assume the existence of number-changing processes in the DS that drive the DS (pseudo)-chemical potentials to zero; this simplification reduces the number of Boltzmann equations that need to be solved by allowing us to set $n_{\varphi}(T_{\mathrm{ds}})=n_{\varphi}^{\mathrm{eq}}(T_{\mathrm{ds}})$. Pauli-blocking and Bose-enhancement factors modify the precise form of the collision term in Eq. 14, but their effects are always to establish thermal and chemical equilibrium between the DS and SM baths. We include quantum-statistical effects in numerical calculations in the model-specific results of Sec. 4.2; the corresponding collision terms are calculated in Appendix A.

In the next two sections, we apply these results to particular scenarios where a light DS equilibrates with either neutrinos or photons.

2.2.1 Coupling to the Neutrino Bath

We use energy (Eq. 11) and entropy (Eq. 12) conservation to understand the evolution of $N_{\mathrm{eff}}$ throughout DS equilibration and decoupling. In the analytic expressions below, we take $T\lesssim m_{e}$ for simplicity, such that $g_{*}^{\gamma}=g_{\gamma}=2$ and $\xi_{\nu}^{\mathrm{SM}}=(4/11)^{1/3}$. Prior to equilibration, the DS contributes to the energy density like “dark radiation” such that

[TABLE]

where $\xi_{\mathrm{ds}}^{0}$ is the DS-to-photon temperature ratio prior to equilibration. Although we take $\xi_{\mathrm{ds}}^{0}$ to be a free parameter, one expects $\xi_{\mathrm{ds}}^{0}<\xi_{\nu}^{\mathrm{SM}}$ if the DS has been diluted by entropy dumps (from SM particles decoupling or otherwise – see, e.g., [29, 30]) or if the DS was not populated efficiently during reheating [31]. It is also possible that the DS is initially hotter than the SM, but this scenario is generically in conflict with considerations of BBN and the CMB regardless of when it attains equilibrium with neutrinos. This is easily seen from Eq. 15 since for $\xi_{\mathrm{ds}}^{0}/\xi_{\nu}^{\mathrm{SM}}>1$ we have $\Delta N_{\mathrm{eff}}\gtrsim 0.6$. The subsequent dynamics of equilibration and decoupling from neutrinos only increase $\Delta N_{\mathrm{eff}}$ as described below, worsening the tension.333This is not obviously the case in the photon-equilibration scenario discussed below, since heating of the photons relative to neutrinos lowers $N_{\mathrm{eff}}$. For this reason, we focus on dark sectors that were initially colder than the SM neutrinos. Note that this dark radiation contribution in Eq. 15 rapidly becomes negligible as $\xi_{\mathrm{ds}}^{0}$ is lowered.

During DS-SM equilibration, the total comoving energy density is conserved (see Eq. 11). Since $N_{\mathrm{eff}}$ is a measure of the total $\text{DS}+\nu$ energy density, it remains unchanged during this epoch. Therefore, immediately after DS-SM equilibration, the form of $N_{\mathrm{eff}}$ is the same as in Eq. 15. Note, however, that the neutrino and DS temperatures evolve non-trivially in order to conserve the comoving energy density (see, e.g., Fig. 1 of Ref. [20]).

As the universe cools, the temperature of the DS-$\nu$ bath may drop below a mass threshold of a DS particle. At this point, processes such as annihilations ($\text{DS}~{}\text{DS}\to\nu\,\nu$) or decays ($\text{DS}\to\nu\,\nu$) efficiently deplete the thermal abundance of this DS particle, transferring its entropy into neutrinos. This slightly increases the neutrino energy density relative to photons. To see this, let us assume that all $g_{*}^{\mathrm{ds}}$ degrees of freedom in the DS equilibrate with neutrinos after neutrino-photon decoupling and later simultaneously (and instantaneously) decouple. In that case, entropy conservation (Eq. 12) implies that

[TABLE]

For $\xi_{\mathrm{ds}}^{0}<1$ and $g_{*}^{\mathrm{ds}}\geq 1$, we therefore have $\Delta N_{\mathrm{eff}}\gtrsim 0.18$. This result should be compared to the one obtained using the standard assumption that the DS species was in equilibrium with the SM before neutrino-photon decoupling. In this case, $\xi_{\mathrm{ds}}^{0}=1$ and $\Delta N_{\mathrm{eff}}$ is generically in conflict with considerations of BBN and the CMB even for $g_{*}^{\mathrm{ds}}=1$ [13, 32] .

The evolution of $N_{\mathrm{eff}}$ for a DS equilibrating with neutrinos is shown schematically in the left panel of Fig. 1. Several detailed examples are also shown in Fig. 2, for various representative values of DS degrees of freedom ($g_{*}^{\mathrm{ds}}$) at a common mass scale ($m_{\mathrm{ds}}$); these were calculated using the Boltzmann equation in Eq. 13. For each curve, the decay rate for DS $\leftrightarrow\nu\nu$ is chosen such that equilibration occurs at $T_{\nu}\approx 10\,m_{\mathrm{ds}}$, fixing the initial DS-to-photon temperature ratio to be $\xi_{\mathrm{ds}}^{0}=0.3$. Figure 2 illustrates that before and during equilibration, $N_{\mathrm{eff}}$ is very close to 3. The DS mass scale, $m_{\mathrm{ds}}$, sets the decoupling temperature. After DS-$\nu$ decoupling, $N_{\mathrm{eff}}$ is given by Eq. 16. If the various DS particles had parametrically different masses, these curves would contain additional “steps” corresponding to the different decoupling events as the temperature crosses these mass thresholds. The value of $N_{\mathrm{eff}}$ before and after each decoupling event can be estimated from repeated applications of comoving entropy conservation.

2.2.2 Coupling to the Photon Bath

Late equilibration of a DS with photons gives rise to a qualitatively different behavior of $N_{\mathrm{eff}}$ compared to the neutrino case described above for two reasons. First, after a new species equilibrates with the photon bath, a given photon temperature (which is effectively used as a “clock” that determines the key cosmological epochs) corresponds to a larger energy density (due to the increased degrees of freedom). This means that $N_{\mathrm{eff}}$ changes during DS-SM equilibrium even though the comoving energy density is conserved during this process. Second, changes to the photon density modify the baryon-to-photon ratio, $\eta_{b}=(n_{b}-n_{\bar{b}})/n_{\gamma}$, which also alters the outcome of primordial nucleosynthesis.

As in the neutrino case, the DS contributes to $N_{\mathrm{eff}}$ as dark radiation prior to equilibration such that

[TABLE]

where, as before, $\xi_{\mathrm{ds}}^{0}$ is the initial DS-to-photon temperature ratio. DS-$\gamma$ equilibration fixes $\xi_{\mathrm{ds}}=1$ and therefore modifies $N_{\mathrm{eff}}$ to

[TABLE]

This equation exhibits two effects that modify $N_{\mathrm{eff}}$ upon DS-$\gamma$ equilibration. The first term in Eq. 18 corresponds to $\xi_{\nu}^{4}$ in the general form of Eq. 10. Depending on the initial temperature of the DS (encapsulated in $\xi_{\mathrm{ds}}^{0}$), DS-$\gamma$ equilibration either deposits energy into or draws energy from the $\gamma$ bath. As a result, the neutrino temperature (relative to photons) is lowered or increased, respectively. The second term in Eq. 18 is the direct contribution from the thermalized DS component and simply counts the additional relativistic degrees of freedom.

Both of these effects lead to non-trivial evolution of $N_{\mathrm{eff}}$ that qualitatively depends on the size of the initial DS temperature ratio, $\xi_{\mathrm{ds}}^{0}$, and the number of DS degrees of freedom, $g_{*}^{\mathrm{ds}}$. Equation 17 implies that $N_{\mathrm{eff}}\gtrsim 3$ before equilibration. However, from Eq. 18 we see that after DS-$\gamma$ recoupling, $N_{\mathrm{eff}}\gtrsim 3$ if $g_{*}^{\mathrm{ds}}\gtrsim g_{\mathrm{critical}}^{\mathrm{ds}}$, or $N_{\mathrm{eff}}\lesssim 3$ if $g_{*}^{\mathrm{ds}}\lesssim g_{\mathrm{critical}}^{\mathrm{ds}}$, where $g_{\mathrm{critical}}^{\mathrm{ds}}$ (the critical number of degrees of freedom in the DS) is given by

[TABLE]

For an initially cold DS ($\xi_{\mathrm{ds}}^{0}<1$), $g_{\mathrm{critical}}^{\mathrm{ds}}$ is unphysical ($g_{\mathrm{critical}}^{\mathrm{ds}}<0<g_{*}^{\mathrm{ds}}$) and hence, $N_{\mathrm{eff}}$ strictly increases during equilibration; this reflects the cooling of the photon bath relative to the neutrinos.

In principle, it is possible to realize a decrease of $N_{\mathrm{eff}}$ in a concrete model. We must demand that $g_{\mathrm{critical}}^{\mathrm{ds}}\gtrsim g_{*}^{\mathrm{ds}}>1$, which is only possible if $\xi_{\mathrm{ds}}^{0}\gtrsim 1.75$. In this case, Eq. 17 implies that $N_{\mathrm{eff}}\gtrsim 24$ before equilibration. Although such non-trivial behavior of $N_{\mathrm{eff}}$ is interesting in this fine-tuned regime, this scenario is usually constrained since it necessarily leads to large deviations in the expansion rate at earlier times which will impact the decoupling of electroweak interactions.

When the temperature of the DS-photon bath drops below a DS mass threshold, processes such as annihilations or decays transfer entropy to the photons, heating them relative to neutrinos and decreasing $\xi_{\nu}$. As a result, $N_{\mathrm{eff}}$ decreases to

[TABLE]

which is strictly less than the SM expectation of $N_{\mathrm{eff}}\approx 3$. The evolution of $N_{\mathrm{eff}}$ throughout DS-$\gamma$ equilibration and decoupling is illustrated schematically in the middle and right panels of Fig. 1, for an initially cold and hot DS, respectively. In the left panel of Fig. 3 we show several examples of $N_{\mathrm{eff}}(T)$ calculated using the Boltzmann equation in Eq. 13.

In addition to its effect on the expansion rate, DS-$\gamma$ equilibration modifies the predictions of BBN by changing the evolution of $\eta_{b}$. Suppose that before DS-$\gamma$ equilibration and $e^{+}e^{-}$ annihilations, the baryon-to-photon ratio is given by $\eta_{b}^{0}$. Since $(n_{b}-n_{\bar{b}})/s_{\nu}$ is constant after neutrinos have decoupled, $\eta_{b}$ at a later time is given by

[TABLE]

This parameter is independently determined by observations of the CMB to be $\eta_{b}^{\mathrm{CMB}}\approx 6.1\times 10^{-10}$ at the time of recombination [32], fixing its value at late times. For a given cosmological history, which determines $\xi_{\nu}(T)$, the initial value of $\eta_{b}^{0}$ must be chosen such that the late-time value (Eq. 21) is compatible with CMB measurements, i.e.,

[TABLE]

where $T_{\text{CMB}}\sim 0.1\text{ eV}$. In the case of DS-$\gamma$ equilibration, the evolution of the neutrino-to-photon temperature ratio, $\xi_{\nu}$, throughout equilibration and decoupling is given by

[TABLE]

Therefore, $\xi_{\nu}(T_{\text{CMB}})$ in Eq. 22 depends on when the DS equilibrates and decouples with respect to recombination. The evolution of the baryon asymmetry (assuming the initial value given by Eq. 22) is shown in the right panel of Fig. 3 for several choices of $g_{*}^{\mathrm{ds}}$ and $m_{\mathrm{ds}}$.

3 Nucleosynthesis

The per-nucleon binding energies of nuclei are at the MeV scale. When the universe cools below this temperature, the formation of nuclei becomes energetically favorable. However, the large abundance of photons compared to baryons ensures that these bound states are not synthesized in significant numbers until much later, when the baryon plasma is even more dilute. As a result, only D, 3He, 4He, and 7Li are formed in appreciable quantities. Helium-4 and D have the largest primordial abundances, which have been measured with several-percent precision. These observations are consistent with standard predictions of BBN, and in the following we use the abundances of these two elements to constrain beyond SM physics. In this section, we explain the effect of variations of $N_{\mathrm{eff}}$ and $\eta_{b}$ on ${}^{4}{\rm He}$ and D yields; we also discuss the numerical treatment of nucleosynthesis for our modified cosmologies and derivation of model constraints.

We do not consider measurements of ${}^{3}{\rm He}$ and ${}^{7}{\rm Li}$ in obtaining bounds on light dark sectors. There is no consensus on the measured ${}^{3}{\rm He}$ primordial abundance due to limited observational regions and uncertainties in theoretical modeling [33, 34]. The measured ${}^{7}{\rm Li}$ abundance is about a factor of three lower than the SM prediction, a long-standing discrepancy known as the “Lithium Problem” [35]. The ${}^{7}{\rm Li}$ abundance can be significantly modified after BBN in halo stars [36]. Alternatively, the discrepancy can be due to physics beyond the SM [37, 38]. Some of the scenarios considered in this work alleviate this tension, but at the price of making ${}^{4}{\rm He}$ and D yields inconsistent with observations. Throughout this work, we therefore assume that the resolution of the lithium problem is astrophysical in nature.

3.1 Dependence on $N_{\mathrm{eff}}$ and $\eta_{b}$

In Sec. 2.2, we examined modifications to $N_{\mathrm{eff}}$ and $\eta_{b}$ that arise as a result of a DS equilibrating and decoupling from either the neutrino or photon bath at temperatures below an MeV. How these (or any related) scenarios impact primordial nucleosynthesis can be understood in terms of a few key epochs that control ${}^{4}{\rm He}$ and D yields: neutron-proton freeze-out, the deuterium bottleneck, and the freeze-out of deuterium burning processes – see, e.g., Refs. [39, 40, 41]. We now briefly review the essential physics for each of these epochs.

In the early universe, neutrons and protons interconvert through electroweak interactions such as $n\,e^{+}\leftrightarrow p\,\bar{\nu}_{e}$. At temperatures smaller than the neutron-proton mass-splitting ($T\lesssim\Delta m_{np}\approx 1.3\text{ MeV}$), the neutron-to-proton number density ratio, $n/p$, becomes smaller than unity, i.e., $n/p\approx\text{exp}(-\Delta m_{np}/T)\lesssim 1$. Eventually, neutron-proton conversion processes freeze-out at $T\approx 0.8\text{ MeV}$ at which point $n/p\approx\text{exp}(-1.3/0.8)\approx 1/5$.444In reality, neutron-proton freeze-out is not instantaneous, but occurs over a range of temperatures. For more details, see, e.g. Ref. [42]. The decays of free neutrons further reduces the neutron density to $n/p\approx 1/7$ before the onset of nucleosynthesis. Nearly all of these neutrons are eventually bound into ${}^{4}{\rm He}$ because of its large binding energy and the corresponding enhanced equilibrium number density. The resulting ${}^{4}{\rm He}$ mass-fraction, $Y_{p}$, is therefore given by

[TABLE]

Although the per-nucleon binding energies of light nuclei are $\mathcal{O}(\text{MeV})$, nucleosynthesis does not occur until much smaller temperatures near $T\sim 100\text{ keV}$. There are two reasons for this: the large entropy of the universe and the small deuterium binding energy. The first fact implies that even if, e.g., ${}^{4}{\rm He}$ was in chemical equilibrium with $n$ and $p$ its abundance would be negligible until $T\sim 300\;\mathrm{keV}$ [28]. The synthesis of these heavier nuclei is delayed even further because their production rates, e.g., via \textrm{D}\;\textrm{D}\to p\;\mbox{{}^{3}{\rm H}} (the first step in making ${}^{4}{\rm He}$), depend on the deuterium abundance. These reaction rates are slow compared to Hubble expansion until $n_{\textrm{D}}/{(n_{p}+n_{n})}\sim\mathcal{O}(10^{-4})$ [41]. This rate-limiting step is often referred to as the “deuterium bottleneck.” The temperature at which the D-burning reactions become important can be estimated from the relative equilibrium abundance,

[TABLE]

where $B_{\rm D}\sim 2\text{ MeV}$ is the deuterium binding energy. This ratio is initially $\ll 1$ because of the large entropy of the plasma (small baryon-per-photon number, $\eta_{b}\sim 10^{-10}$). This abundance becomes large enough to support the production of heavier nuclei only when the temperature falls below $\sim 100\;\mathrm{keV}$. At this point, nucleosynthesis commences, and most of the deuterium is eventually processed into ${}^{4}{\rm He}$. However, deuterium-burning processes (such as \textrm{D}\;p\to\mbox{{}^{3}{\rm He}}\;\gamma) remain effective until they freeze out at temperatures of $T\sim\text{few}\times 10\text{ keV}$ and the deuterium abundance has been depleted down to trace amounts of $\textrm{D}/\textrm{H}\sim 10^{-5}-10^{-4}$.

We show how theoretical predictions of the neutron-to-proton ratio ($n/p$), the ${}^{4}{\rm He}$ mass-fraction ($Y_{p}$), and the deuterium fraction (D/H) are affected by modifications to the expansion rate in Fig. 4 and the baryon density in Fig. 5, encapsulated in $N_{\mathrm{eff}}$ and $\eta_{b}$, respectively. Values of $N_{\mathrm{eff}}$ that are greater or smaller than the SM prediction of $N_{\mathrm{eff}}\approx 3$, correspond to an increased or decreased expansion rate during nucleosynthesis. For $N_{\mathrm{eff}}\gtrsim 3$, neutron-proton interconverting processes decouple at earlier times (compared to the SM prediction), leading to an increase in the predicted value of $n/p$ and hence a larger ${}^{4}{\rm He}$ yield (see Eq. 24). Also for $N_{\mathrm{eff}}\gtrsim 3$, D-burning reactions freeze out earlier in the more rapidly expanding universe, which results in a larger predicted deuterium abundance. Conversely, $Y_{p}$ and D/H decrease for $N_{\mathrm{eff}}\lesssim 3$.

The ${}^{4}{\rm He}$ mass-fraction is only logarithmically sensitive to the baryon density, $\eta_{b}$. This effect arises due to how $\eta_{b}$ controls the deuterium bottleneck. Defining the temperature at which nucleosynthesis commences as $T_{\text{BBN}}\sim\mathcal{O}(100)\text{ keV}$ through the criterion $n_{\rm D}/(n_{p}+n_{n})\sim\mathcal{O}(10^{-4})$ [41], we see from Eq. 25 that

[TABLE]

where in the second equality, we have considered small variations of $\eta_{b}$ around the CMB-preferred value, i.e., $\eta_{b}\sim\mathcal{O}(10^{-10})+\delta\eta_{b}$. Hence, for much larger $\eta_{b}$, nucleosynthesis begins at earlier times, and less unbound neutrons decay before ending up in ${}^{4}{\rm He}$. As a result, $Y_{p}$ increases for larger $\eta_{b}$. The onset of nucleosynthesis corresponds to times much shorter than the neutron lifetime, such that $Y_{p}$ is only power-law (not exponentially) dependent on $T_{\text{BBN}}$. Unlike ${}^{4}{\rm He}$, deuterium is much more sensitive to deviations in the baryon density, since $\eta_{b}$ directly controls the strength of D-burning rates. In particular, larger $\eta_{b}$ corresponds to more efficient D-burning and a smaller trace abundance at late times.

3.2 Procedure

The primordial ${}^{4}{\rm He}$ and D densities are inferred from direct observations of various astrophysical sites. For instance, the primordial ${}^{4}{\rm He}$ abundance is measured from observations of recombination emission lines originating from low-metallicity HII regions of galaxies [43, 44, 45]. The observed deuterium abundance constitutes a lower bound on its primordial value since it is easily destroyed in stellar cycles. It is directly observed from, e.g., isotope-shifted $\text{Ly}\alpha$ features in the absorption spectra of distant quasars [46, 47]. In this work, we adopt the recommended values of Ref. [48] for the observed abundances,

[TABLE]

We calculate ${}^{4}{\rm He}$ and D yields by modifying version 1.4 of the publicly available code AlterBBN [49]. Similar to many public BBN programs (see, e.g., Refs. [50, 51]), AlterBBN follows the structure and techniques of the seminal Kawano code [52]. These codes take cosmological parameters and tabulated nuclear reaction rates as inputs and estimate the primordial abundances by solving a set of differential equations governing BBN. For more details, see Refs. [52, 53, 54].

The default version of AlterBBN computes primordial abundances assuming a standard cosmological history. In Sec. 2, we discussed how light dark sectors that equilibrate during nucleosynthesis lead to temperature-dependent modifications to the expansion rate and the baryon density. We have modified AlterBBN in order to calculate the primordial nuclei yields for such general temperature-dependent forms of $N_{\mathrm{eff}}(T)$ and $\eta_{b}(T)$.555Throughout this work, we assume that the temperature of neutrino-photon decoupling ($\sim\text{few}\times\text{MeV}$) is unchanged.

The particular values of nuclear reaction rates and their experimental errors have a large impact on the predicted primordial abundances and their theoretical uncertainties. The default rates used in AlterBBN are listed in Table 1 of Ref. [49]. We have updated the nuclear reactions that are relevant to ${}^{4}{\rm He}$ and D with newer calculations, largely following the choices of Ref. [13]. The modified rates are summarized in Table 2. We mostly utilize results from the NACRE collaboration, which compiles and evaluates the latest updated cross sections, as given in the NACRE-II compilation [55]. However, two important reactions are not provided in the NACRE-II compilation: $p\;n\to\textrm{D}\;\gamma$ and \mbox{{}^{3}{\rm He}}\;n\to p\;\mbox{{}^{3}{\rm H}}. For these, we use calculations from Refs. [56, 57], respectively. We note that various reactions rates have been updated since the NACRE-II release. In particular, Ref. [58] updated three important reactions of D burning. Not only do they shift the central values of these rates by $\sim 10$%, they also report much smaller errors. Both are crucial for constraining new physics scenarios with observations of deuterium. In our baseline analysis, we follow Ref. [58] for these three rates because the methodology used to calculate the central values and their spreads is comprehensively documented, with special attention paid to experimental systematic errors and theoretically motivated fitting functions. In Sections 3.3 and 4.2, we compare results for different choices of D-burning rates. In the models we study in this work, the neutrino temperature is modified compared to the SM, which changes the $p\leftrightarrow n$ conversion rate. We have checked that varying the neutrino temperature by a factor of a few negligibly affects the conversion rate and we therefore ignore this effect.

Theoretical uncertainties of the predicted ${}^{4}{\rm He}$ and D yields stem from uncertainties in various nuclear reaction rates and the neutron lifetime. We determine the corresponding error bars from a Monte Carlo procedure similar to that outlined in Ref. [13]. We estimate the theoretical uncertainties for ${}^{4}{\rm He}$ and D yields by lognormal-sampling nuclear rates [59, 60] and computing the yields $10^{4}$ times for each value of $N_{\mathrm{eff}}$ and $\eta_{b}$. The resulting yield distribution is well-described by a correlated Gaussian likelihood. When we utilize the D-burning rates from the NACRE-II compilation (see the last three rows in Table 2), we obtain uncertainties similar to those of Ref. [13]:

[TABLE]

where we have taken the SM expectation of $N_{\mathrm{eff}}\approx 3$, in the first line we have fixed the late-time baryon density to the central value inferred from observations of the CMB, and in the second line we have additionally incorporated the spread around the central value, $\eta_{b}^{\text{CMB}}=(6.10\pm 0.04)\times 10^{-10}$ [13, 61]. In the first and second lines of Eq. 28, we find a correlation coefficient of $\approx-0.22$ and $\approx-0.26$, respectively. The anti-correlation is driven primarily by the sizable spread in \textrm{D}\;\textrm{D}\to n\;\mbox{{}^{3}{\rm He}} and \textrm{D}\;\textrm{D}\to p\;\mbox{{}^{3}{\rm H}} (larger values of these rates reduce D/H and enhance ${}^{4}{\rm He}$ abundances [13]).

Compared to the NACRE-II compilation, the uncertainties of important D-burning processes (last three rows of Table 2) are much smaller in Ref. [58]. If we instead adopt these rates from Ref. [58], we find that for $N_{\mathrm{eff}}\approx 3$,

[TABLE]

where $\eta_{b}$ is fixed as in Eq. 28. We find a correlation coefficient of $\approx 0.006$ and $\approx-0.10$ in the first and second lines, respectively; the correlation is mostly due to the $\eta_{b}$ sensitivity of the two yields, as discussed in Sec. 3.1. The central value for the predicted deuterium abundance in Eq. 29 is slightly smaller than that of Eq. 28. This shift is due to the larger D-burning rates of Ref. [58].

The theoretical uncertainty of the ${}^{4}{\rm He}$ abundance is negligible compared to the observational one. This is because the ${}^{4}{\rm He}$ abundance is mainly sensitive to the neutron-to-proton ratio. Therefore, the dominant theoretical uncertainty for $Y_{p}$ comes from variations of the neutron lifetime, which we take from Ref. [48]:

[TABLE]

There is a well-established tension between bottle and beam measurements of the neutron lifetime [62]. The value quoted above is an average that is dominantly determined by bottle experiment measurements. We note that even if the true value of $\tau_{n}$ is closer to that inferred from beam-based measurements ($\sim 888$ s), the corresponding shift to $Y_{p}$ and $\textrm{D}/\textrm{H}$ would be at the sub-percent level.

We find that the fractional theoretical uncertainties, $\sigma(X)/X$ for $X=Y_{p}$ and $\textrm{D}/\textrm{H}$, are to a good approximation independent of the functional form of $N_{\mathrm{eff}}(T)$ and $\eta_{b}(T)$. This can be understood as follows. Reference [63] has shown that the nuclear rate uncertainties in the final yields are well-described using linear error propagation. As a result, the sensitivity of yields of element $X$ to the rates, $\Gamma_{i}$, is encompassed by a set of constants, $\alpha_{i}$ (the logarithmic derivatives of the yields with respect to the rates), such that

[TABLE]

In the nucleosynthesis Boltzmann system, these rates enter as $\Gamma_{i}/H$, which suggests that the effect of variations of $H$ (via $N_{\mathrm{eff}}$) can be likewise linearized, at least for small deformations of $H$, leading to $X\propto{N_{\mathrm{eff}}}^{\alpha^{\prime}}$, for some other constant $\alpha^{\prime}$. The resulting set of logarithmic derivatives, $\{\alpha_{i},\alpha^{\prime}\}$, is presented in Ref. [13]. A similar argument can be made for modifications to $\eta_{b}$. Using this linearized form of the yields, it is simple to show that fractional errors from nuclear rate uncertainties, $\sigma(X)/X$, are constant as $N_{\mathrm{eff}}$ and $\eta_{b}$ are varied. This is expected to hold for small perturbations to these quantities (as long as the linearization is valid); fortunately, the measured and predicted abundances of the light elements are so precise that only small deviations from standard values are allowed. While these arguments are straightforward for constant shifts to $N_{\mathrm{eff}}$ and $\eta_{b}$, they are more difficult to make for arbitrary time-variations of these quantities. Therefore, we have explicitly checked the constancy of $\sigma(Y_{p})/Y_{p}$ and $\sigma(\textrm{D}/\textrm{H})/(\textrm{D}/\textrm{H})$ for time-dependent variations of $N_{\mathrm{eff}}$ and $\eta_{b}$, as considered in the following sections, using the Monte Carlo approach described above. This justifies our use of $\sigma(Y_{p})/Y_{p}$ and $\sigma(\textrm{D}/\textrm{H})/(\textrm{D}/\textrm{H})$ from Eqs. 28 and 29 throughout this work. Although we adopt Eq. 29 for our baseline analysis, in Secs. 3.3 and 4.2 we show how these two rate choices affect the constraints for constant shifts to $N_{\mathrm{eff}}$ and $\eta_{b}$ and in concrete particle physics models, respectively.

3.3 Standard Constraints on $N_{\mathrm{eff}}$ and $\eta_{b}$

Before we discuss our main results in Sec. 4, we first present a simple estimate for the standard bounds on $N_{\mathrm{eff}}$ and $\eta_{b}$, assuming that they take constant values throughout the epoch of primordial nucleosynthesis. These results exemplify our methods that we will apply again in Sec. 4.

We utilize the methodology outlined in Sec. 3.2 in order to calculate ${}^{4}{\rm He}$ and D yields, $(Y_{p})^{\text{theory}}$ and $(\textrm{D}/\textrm{H})^{\text{theory}}$, respectively. We then compare these predictions to the observed abundances in Eq. 27 through the $\chi^{2}$ test statistic,

[TABLE]

where we have defined

[TABLE]

Above, $\sigma_{Y_{p},\textrm{D}/\textrm{H}}^{\text{obs.}}$ and $\sigma_{Y_{p},\textrm{D}/\textrm{H}}^{\text{theory}}$ are the observational and theoretical uncertainties for $Y_{p}$ and $\textrm{D}/\textrm{H}$, as detailed in Sec. 3.2. Fixing $N_{\mathrm{eff}}$ and $\eta_{b}$ to their values expected in a standard cosmology, we find that

[TABLE]

This gives $\chi^{2}\approx 5.7$, assuming a standard cosmological history.

The above definition of $\chi^{2}$ does not account for correlated uncertainties. This is justified in our baseline analysis, since we adopt the nuclear rates of Ref. [58], in which case the correlation coefficients are much smaller than unity (see the discussion below Eq. 29). However, this is not the case for the NACRE-II rates of Ref. [55] (Eq. 28). In most regions of parameter space that we consider throughout this work, it is typically the case that $\Delta Y_{p}\gg\Delta\textrm{D}/H$ or $\Delta Y_{p}\ll\Delta\textrm{D}/H$, and, hence, we do not expect correlated uncertainties to significantly modify our results. To check this, we have explicitly rerun our analysis accounting for these correlations. We find that they have negligible impact in the majority of the parameter space. Hence, when displaying our results, we will ignore such correlations between $Y_{p}$ and $\textrm{D}/\textrm{H}$ even when adopting the NACRE-II rates.

Figure 6 illustrates how modifications to the baryon density ($\eta_{b}$) and the expansion rate ($N_{\mathrm{eff}}$) are constrained from measurements of the primordial densities of helium-4 and deuterium. The parameters, $\Delta\eta_{b}$ and $\Delta N_{\mathrm{eff}}$, correspond to a shift of the late-time (CMB era) value of $\eta_{b}$ and a time-independent modification of $N_{\mathrm{eff}}$ away from the SM expectation, respectively. The evolution of $\eta_{b}(T)$ during nucleosynthesis for a given $\Delta\eta_{b}$ is evaluated using entropy conservation, as in Eqs. 21 and Eq. 22, assuming the standard form of $\xi_{\nu}$. The SM is defined by $N_{\mathrm{eff}}\approx 3$, with the baryon-to-photon ratio fixed at late times to the CMB-preferred value, $\eta_{b}\approx 6\times 10^{-10}$.

In our calculations of the primordial nuclei abundances, the best-fit point is defined as the value of $\Delta N_{\mathrm{eff}}$ and $\Delta\eta_{b}$ that minimizes $\chi^{2}$, as defined in Eq. 32. In Fig. 6, this occurs at $(\Delta N_{\mathrm{eff}},10^{10}\Delta\eta_{b})=(-0.04,-0.2)$ with a $\chi^{2}$ of $\chi_{\text{min}}^{2}=0.04$. In Fig. 6, we show $1\sigma$ (dark blue) and $2\sigma$ (light blue) regions around this best-fit point, corresponding to $\Delta\chi^{2}\equiv\chi^{2}-\chi_{\text{min}}^{2}=2.30$ and $6.18$, respectively, using the nuclear rates of Ref. [58]. Standard cosmology is consistent with the observed abundances within $\lesssim 2\sigma$. Along the dot-dashed and dotted gray contours of Fig. 6, $|\Delta Y_{p}|=2$ and $|\Delta\textrm{D}/\textrm{H}|=2$, respectively. As discussed in Sec. 3.1, $Y_{p}$ is dominantly sensitive to $N_{\mathrm{eff}}$ with only a logarithmic dependence on $\eta_{b}$. On the other hand, the effect of $\eta_{b}$ on the predicted deuterium abundance is much larger, but is degenerate with $N_{\mathrm{eff}}$, since modifications to the expansion rate alter the freeze-out temperature of deuterium-burning, which can always be compensated by increasing or decreasing the burning rates with larger or smaller values of $\eta_{b}$.

Compared to Ref. [13], our baseline analysis shows a mild preference for smaller values of $\eta_{b}$ and slightly smaller uncertainties because of our different choices for various D-burning rates (see the discussion in Sec. 3.2). In particular, the rates listed in the last three rows of Table 2 are larger and have smaller uncertainties compared to those used in Ref. [13], as reflected in Eqs. 28 and 29. In order to illustrate this point explicitly, we perform the same analysis, but instead adopt the D-burning rates of the NACRE-II compilation [55]. The corresponding $2\sigma$-favored region is shown by the dark green contour of Fig. 6. This result is similar to that of Ref. [13], with small discrepancies due to our different choices for the observed primordial abundances in Eq. 27. We have explicitly checked that our results are consistent with those in Ref. [13] when we adopt the nuclear rates of the NACRE-II compilation and the inferred yields of primordial nuclei noted in Ref. [13]. Also shown in Fig. 6 are regions consistent with recent Planck measurements of the CMB [32].

4 Equilibration and Decoupling during Nucleosynthesis

In Sec. 3, we reviewed how the expansion rate ($N_{\mathrm{eff}}$) and the baryonic density ($\eta_{b}$) control the outcome of primordial nucleosynthesis and how time-independent deviations of these quantities away from their SM expectations lead to changes in the predicted abundances of helium-4 and deuterium. In most studies of BBN, constant shifts to $N_{\mathrm{eff}}$ and $\eta_{b}$ are constrained in this manner.

Light and feebly-coupled dark sectors (DS) that enter equilibrium with the SM bath below the temperature of neutrino-photon decoupling naturally lead to deviations in $N_{\mathrm{eff}}$ and $\eta_{b}$ that are effectively temperature- or time-dependent, as shown in Figs. 1-3. In this case, modifications to primordial nucleosynthesis may occur in specific time/temperature intervals, and adapting bounds from previous studies is not straightforward. In this section, we discuss this more general scenario, in which $N_{\mathrm{eff}}$ and $\eta_{b}$ evolve non-trivially as the universe adiabatically cools below critical temperatures and mass-thresholds of the DS. In Sec. 4.1, we present constraints from considerations of BBN in a model-independent manner, for specific simplified forms of $\Delta N_{\mathrm{eff}}(T)$ (where $T$ is the temperature of the photon bath) and fixing $\Delta\eta_{b}=0$. In Sec. 4.2, we consider a few concrete models that predict non-standard temperature evolution of $\Delta N_{\mathrm{eff}}(T)$ and $\Delta\eta_{b}(T)$ and discuss how the predicted nuclear abundances are modified in relevant regions of parameter space.

4.1 Model-Independent Results

In this section, we only focus on model-independent modifications to the expansion rate and adopt a few representative “temperature waveforms” for the functional form of $N_{\mathrm{eff}}(T)$. The first possibility that we consider is that $N_{\mathrm{eff}}$ tracks the SM expectation of $N_{\mathrm{eff}}\approx 3$ until some later temperature, $T_{\text{on}}$, at which point a deviation turns on, i.e., $\Delta N_{\mathrm{eff}}\neq 0$ for $T\lesssim T_{\text{on}}$. This is shown in the left panel of Fig. 7. As discussed in Sec. 2.2.1, this scenario is motivated by a sub-MeV DS that recouples and later decouples with the neutrino bath after neutrino-photon decoupling.

We parametrize this temperature-behavior as

[TABLE]

where the temperature-dependence is encapsulated in the Heaviside step-function, $\Theta$. We also consider the possibility that $\Delta N_{\mathrm{eff}}\neq 0$ at the beginning of nucleosynthesis, but that this deviation turns off at some critical temperature, $T_{\text{off}}$, i.e.,

[TABLE]

as shown in the right panel of Fig. 7. In general, dark sectors that recouple and later decouple with the SM bath during nucleosynthesis can result in more intricate behaviors of $\Delta N_{\mathrm{eff}}(T)$. We therefore also investigate the generalization of Eqs. 35 and 36,

[TABLE]

This corresponds to a sudden pulse-like modification to $N_{\mathrm{eff}}$ within the temperature range of $T_{\text{off}}\lesssim T\lesssim T_{\text{on}}$.

If these deviations to $N_{\mathrm{eff}}$ occur in the temperature range $\mathcal{O}(10)\text{ keV}\lesssim T_{\text{on,off}}\lesssim\mathcal{O}(1)\text{ MeV}$, then important epochs of nucleosynthesis are potentially modified. Figures 8 and 9 show the fractional shift (compared to $\Delta N_{\mathrm{eff}}=0$) to the helium and deuterium abundances for these various temperature evolutions of $\Delta N_{\mathrm{eff}}(T)$, fixing the amplitude to $\Delta N_{\mathrm{eff}}=0.5$. The final abundance of ${}^{4}{\rm He}$ is modified in the standard way (see Sec. 3.1) if $\Delta N_{\mathrm{eff}}\neq 0$ during neutron-proton freeze-out ($100\text{ keV}\lesssim T\lesssim\text{ MeV}$), while the deuterium abundance is most sensitive to the expansion rate near the freeze-out of D-burning process ($50\text{ keV}\lesssim T\lesssim 100\text{ keV}$). However, it is important to note that drastic modifications to the neutron-to-proton ratio do also affect deuterium yields since its final abundance ultimately depends upon the neutrons that do not end up bound into helium nuclei. This effect can be seen by the feature in the green contours near $T\sim\mathrm{MeV}$ in Figs. 8 and 9, which shows that the deuterium abundance is indeed slightly affected by modifications to the expansion rate that occur well before the deuterium bottleneck is overcome.

In the left and right panels of Fig. 10, we highlight regions of parameter space that are inconsistent with considerations of BBN in the $\Delta N_{\mathrm{eff}}-T_{\text{on}}$ and $\Delta N_{\mathrm{eff}}-T_{\text{off}}$ plane, respectively, assuming a step-like temperature evolution of $\Delta N_{\mathrm{eff}}(T)$, as parametrized by Eqs. 35 and 36. Also shown as blue and green contours are parameters for which $|\Delta Y_{p}|=3$ and $|\Delta\textrm{D}/\textrm{H}|=3$, respectively (see Eq. 33). For this analysis, we adopt the same prescription as described in Sec. 3.3 and exclude parameters that predict $\Delta\chi^{2}\gtrsim 6.18$. Throughout, we have fixed the baryon density, $\eta_{b}$, such that it agrees with the measured value at the time of recombination, as given by Eq. 22, and tracks the standard evolution at earlier times. If the cosmological expansion rate is modified during $50\text{ keV}\lesssim T\lesssim\text{MeV}$, then our calculated bound on $N_{\mathrm{eff}}$ is similar to the standard constraint presented in Fig. 6, i.e., $|\Delta N_{\mathrm{eff}}|\lesssim 0.5$. However, if modifications to the expansion rate occur when $T\lesssim 50\text{ keV}$ or $T\gtrsim\text{MeV}$, then deviations as large as $|\Delta N_{\mathrm{eff}}|\gtrsim\mathcal{O}(1)$ are consistent with the measured abundances of helium-4 and deuterium.

As discussed in Sec. 3.1, measurements of primordial deuterium are largely sensitive to $N_{\mathrm{eff}}$ deviations that occur slightly before the end of nucleosynthesis. This is illustrated in the left panel of Fig. 10, which shows that scenarios in which $\Delta N_{\mathrm{eff}}\neq 0$ only for $T\lesssim\text{few}\times\mathcal{O}(100)\text{ keV}$ are predominantly constrained by measurements of the deuterium abundance. Also apparent in Fig. 10 is the asymmetrical importance of $Y_{p}$ and $\textrm{D}/\textrm{H}$ for constraining $\Delta N_{\mathrm{eff}}>0$ or $\Delta N_{\mathrm{eff}}<0$. This can be understood from the fact that the SM ($\Delta N_{\mathrm{eff}}=\Delta\eta_{b}=0$) predicts a slight overabundance in helium and a more significant underabundance in deuterium (see Eq. 34). As a result, most of the constraining power for $\Delta N_{\mathrm{eff}}<0$ comes from measurements of $\textrm{D}/\textrm{H}$, since such cosmologies lead to an even more significant underabundance of deuterium, as explained in Sec. 3.1.

Constraints on pulse-like evolutions of $\Delta N_{\mathrm{eff}}$ are shown in Fig. 11 (see Eq. 37). In each panel of Fig. 11, we vary the fractional width of the $N_{\mathrm{eff}}$ pulse, $(T_{\text{on}}-T_{\text{off}})/T_{\text{on}}$, as well as the temperature at which the pulse turns on, $T_{\text{on}}$, for different choices of the pulse height, $\Delta N_{\mathrm{eff}}=\pm 1$, $\pm 2$. Regions of parameter space that are inconsistent with the predictions of BBN are shown in solid gray. This analysis is nearly identical to ones shown earlier in this work, but since we are now varying parameters in the three-dimensional parameter space spanned by $T_{\text{on}}$, $T_{\text{off}}$, and $N_{\mathrm{eff}}$, we exclude parameters at the $2\sigma$ level if $\Delta\chi^{2}\gtrsim 8.03$. Figure 11 illustrates that significant modifications to the expansion rate are allowed to occur during nucleosynthesis, provided that these deviations happen after neutron-proton freeze-out and before deuterium burning, i.e., when the temperature of the photon bath is $100\text{ keV}\lesssim T\lesssim\mathrm{MeV}$. As in Fig. 10, measurements of the deuterium abundance are typically more constraining for $\Delta N_{\mathrm{eff}}<0$. On the other hand, for $\Delta N_{\mathrm{eff}}>0$, the relative importance of helium and deuterium strongly depends on when the deviation to $N_{\mathrm{eff}}$ takes place, which is most clearly illustrated in the bottom-left panel of Fig. 11. As discussed above and shown in Fig. 9, deuterium yields are dominantly sensitive to the cosmological expansion rate near the end of nucleosynthesis, with a subdominant sensitivity to earlier epochs near neutron-proton freeze-out. The structure of the exclusion at lower values of $(T_{\text{on}}-T_{\text{off}})/T_{\text{on}}$ (corresponding to narrower pulses) visible in the left column of Fig. 11 can be understood by comparison with Fig. 9.

4.2 Model-Specific Results

In this section, we consider two concrete minimal models in which a DS equilibrates with the neutrino or photon bath after neutrino-photon decoupling via decays and inverse decays of a light bosonic mediator, $\varphi$: $\varphi\leftrightarrow\nu\nu\;\mathrm{or}\;\varphi\leftrightarrow\gamma\gamma$. These processes can be realized by the interactions given in Eqs. 5 and 6 for decays to SM neutrinos and photons, respectively. Throughout this section, we use the notation introduced in Sec. 2. Since equilibration and decoupling between the DS and the SM do not occur instantaneously, we solve the Boltzmann equations to find $N_{\mathrm{eff}}(T)$ and $\eta_{b}(T)$ (see Eqs. 13 and 14). The relevant collision terms are given in Appendix A. Solutions of the Boltzmann equations for different illustrative values of the common DS mass scale, $m_{\mathrm{ds}}$, and DS degrees of freedom, $g_{*}^{\mathrm{ds}}$, are shown in Figs. 2 and 3 for neutrino and photon couplings, respectively.

We evaluate the light element abundances for these various cosmologies and compare them to the observed values, as in Sec. 3. In Fig. 12, we show regions of parameter space for the neutrino-coupled model that are excluded from considerations of BBN, as a function of the DS mass scale ($m_{\mathrm{ds}}$) and the DS effective number of relativistic degrees of freedom ($g_{*}^{\mathrm{ds}}$) for an initial temperature ratio of $\xi^{0}_{\mathrm{ds}}=0.3$. The $2\sigma$-excluded regions (solid gray) correspond to $\Delta\chi^{2}=6.18$. Gray regions outlined by solid and dashed lines are obtained using two different sets of D-burning rates, as discussed in Sec. 3.2. The equilibration and decoupling of a cold DS with the neutrino bath gives rise to a step-like deformation in $N_{\mathrm{eff}}$ (see Eq. 35) occurring at $T_{\mathrm{on}}\sim m_{\mathrm{ds}}/5$, as shown in Fig. 2. Hence, the shape of these exclusions can be understood from the left panel of Fig. 8. For larger DS masses ($m_{\mathrm{ds}}\gg\text{MeV}$), $T_{\text{on}}\gtrsim\text{MeV}$ and $N_{\mathrm{eff}}$ is modified during every key epoch of nucleosynthesis, which is very strongly constrained for $g_{*}^{\mathrm{ds}}\gtrsim\text{few}$. As the DS mass scale is decreased, the expansion history is altered for a correspondingly shorter time/temperature interval during nucleosynthesis. Hence, DS masses lighter than $\sim 100\text{ keV}$ are typically unconstrained from considerations of BBN, since in this case $N_{\mathrm{eff}}$ is modified only after the most important processes of nucleosynthesis have concluded. Precisely the same behavior is seen in the model-independent constraint derived in the previous section and shown in the left panel of Fig. 10. In fact, the model-independent result can be used to approximate the bound on $g_{*}^{\mathrm{ds}}$ in Fig. 12 by using $T_{\mathrm{on}}\sim m_{\mathrm{ds}}/5$ and Eq. 16 to translate $\Delta N_{\mathrm{eff}}$ into $g_{*}^{\mathrm{ds}}$.

The horizontal dotted line in Fig. 12 shows the $2\sigma$ upper bound on $g_{*}^{\mathrm{ds}}$ derived from Eq. 16 and Planck measurements of $N_{\mathrm{eff}}$ at recombination666We use the Planck result from the joint fit for $N_{\mathrm{eff}}$ and $Y_{p}$, which does not assume standard BBN. In this case, marginalizing over $Y_{p}$ in the resulting likelihood gives $N_{\mathrm{eff}}=2.99^{+0.43}_{-0.40}$ at $2\sigma$ for the Planck TT,TE,EE+lowE+lensing+BAO datasets [32]. [32]:

[TABLE]

For this neutrino-coupled scenario, the CMB bound is always more stringent than the one derived from considerations of BBN. There are several reasons for this. First, the CMB is sensitive to the expansion rate at much later times, when the presence of the DS has maximally modified $N_{\mathrm{eff}}$. In contrast, bounds derived from primordial nucleosynthesis depend on precisely when modifications to $N_{\mathrm{eff}}$ commence. For instance, if the mass of a cold neutrino-coupled DS particle is smaller than the temperature at which D-burning processes freeze out, $N_{\mathrm{eff}}$ is SM-like throughout the most important epochs of nucleosynthesis. Hence, bounds derived from BBN are weakened for $m_{\mathrm{ds}}\lesssim 100\text{ keV}$. Second, these constraints also depend on the particular set of D-burning rates adopted in the analysis, as discussed in Sec. 3. When nuclear rates from the NACRE-II compilation [55] are adopted, theoretical uncertainties for deuterium yields are large and matching to the observed nuclei abundances (Eq. 27) leads to a preference for a SM-like expansion rate ($N_{\mathrm{eff}}\approx 3$). Instead, when rates from Coc et al. (2015) [58] are utilized, theoretical uncertainties of deuterium yields are significantly reduced (compare Eqs. 28 and 29) and the larger D-burning rates lead to a corresponding preference for a larger expansion rate (compared to NACRE-II). Both of these effects tend to substantially weaken the BBN constraints on $g_{*}^{\mathrm{ds}}$, compared to those derived from the CMB. The point in parameter space that is preferred from considerations of BBN alone (for the Coc et al. (2015) rates [58]), shown as the red star in Fig. 12, realizes a larger expansion and provides a better fit to the observed D abundance than the SM if $\eta_{b}$ is fixed to the CMB-preferred value. This parameter point, however, is robustly ruled out by Planck [32], since it corresponds to $\Delta N_{\mathrm{eff}}\sim 0.8$ at recombination.

Compared to neutrino-coupled models, the relative importance of the CMB and BBN in constraining light dark sectors is reversed in the photon-equilibration case. This is shown in Fig. 13. Here, we have fixed the DS-photon equilibration temperature such that $T^{\text{ds eq}}=3\,m_{\mathrm{ds}}$, which ensures that the DS degrees of freedom are relativistic at equilibration and therefore contribute fully to the expansion rate. This choice of parameters allows for the largest viable range of DS masses to be considered, as we justify below. The solid and dashed lines indicate BBN constraints that are obtained by using the nuclear rates from Coc et al. (2015) [58] and NACRE-II [55], respectively. The dotted line in Fig. 13 shows the $2\sigma$ upper bound on $g_{*}^{\mathrm{ds}}$ that is derived from Eq. 20 and Planck measurements of $N_{\mathrm{eff}}$ at recombination [32]:

[TABLE]

As shown in Sec. 2.2.2 and Fig. 3, modifications to the expansion rate at late times (relevant for CMB measurements) can be much smaller than those during nucleosynthesis in these models. As a result, the measured abundances of ${}^{4}{\rm He}$ and D provide the leading constraints on such new degrees of freedom. In addition to its effect on the expansion rate, equilibration and decoupling of a DS with the photon bath modifies the evolution of the baryon-to-photon ratio, $\eta_{b}$. If we were to ignore this effect on $\eta_{b}$, then BBN would only constrain a finite interval of $m_{\mathrm{ds}}$ (for a fixed value of $T^{\text{ds eq}}/m_{\mathrm{ds}}$), reflecting the pulse-like modification of the expansion rate, as shown in Fig. 3. If this pulse in the expansion rate occurs well before or well after nucleosynthesis (corresponding to large and small $m_{\mathrm{ds}}$, respectively), then the predicted yields are not modified relative to the SM case. However, once we account for modifications to $\eta_{b}$, the resulting bounds extend to much lower DS masses. The baryon-to-photon ratio is independently measured at recombination, and hence we fix $\eta_{b}$ at earlier times as in Eq. 21. Dark sector-photon equilibration and decoupling results in an irreducible increase to the comoving entropy. As a result, before the onset of nucleosynthesis, $\eta_{b}$ has to be larger than in a standard cosmology, if it is fixed to the CMB-preferred value at later times. Hence, even if DS-photon equilibration and decoupling occur well after nucleosynthesis has concluded, the primordial abundances of nuclei are still modified by the non-standard value of $\eta_{b}$ during nucleosynthesis. As a result, considerations of BBN lead to constraints that extend down to sub-keV DS masses.

The strength of the BBN constraint on $g_{*}^{\mathrm{ds}}$ in Fig. 13 can be understood by comparing the maximum deviation of $N_{\mathrm{eff}}$ during photon equilibration (Eq. 18) to the constraint on standard constant shifts of $\Delta N_{\mathrm{eff}}\lesssim 0.5$ in Fig. 6: solving for $g_{*}^{\mathrm{ds}}$, one finds that the expected bound should lie near $g_{*}^{\mathrm{ds}}\sim 0.1$, which is borne out in the complete calculation of Fig. 13. The observed abundances of light nuclei constrain $g_{*}^{\mathrm{ds}}<1$, thereby excluding even a minimal DS that contains a single scalar degree of freedom that equilibrates with photons while relativistic and after neutrino-photon decoupling. Instead, if equilibration occurs while the DS is semi- or non-relativistic ($T^{\mathrm{ds}\;\mathrm{eq}}<3m_{\mathrm{ds}}$), then the effective value of $g_{*}^{\mathrm{ds}}$ can be smaller than unity even if there exist several light states in the DS. This process of semi-relativistic equilibration leads to a smaller shift in $N_{\mathrm{eff}}$ than what is expected from the estimates in Sec. 2.2 [21]. However, this scenario also requires a coincidence of scales, i.e., $m_{\mathrm{ds}}/4\lesssim T^{\mathrm{ds}\;\mathrm{eq}}<3m_{\mathrm{ds}}$. This is because for $T\lesssim m_{\mathrm{ds}}/4$, the DS-SM equilibration rate no longer increases relative to the Hubble parameter with the expansion of the universe, while for $T\gtrsim 3m_{\mathrm{ds}}$, $g_{*}^{\mathrm{ds}}\geq 1$ and equilibration is relativistic. Furthermore, even for scenarios in which $T^{\mathrm{ds}\;\mathrm{eq}}\sim m_{\mathrm{ds}}/4$ (corresponding to $g_{*}^{\mathrm{ds}}\gtrsim 0.25$ during equilibration), there remains a slight tension with the BBN bounds of Fig. 13.

Up to this point, we have ignored the presence of the photon plasma frequency. For temperatures above the electron mass ($T\gtrsim m_{e}$), the effective mass of transverse photon excitations is $m_{t}^{2}\approx 4\pi\alpha T^{2}/9$, which falls quickly to $m_{t}^{2}\approx 4\pi\alpha n_{e}/m_{e}$ as $T$ drops below $m_{e}$ [64]. In certain regions of parameter space, this mass kinematically forbids the equilibrating processes $\varphi\leftrightarrow\gamma\gamma$ if $m_{\mathrm{ds}}=m_{\varphi}<2\,m_{t}(T)$ [19]. However, for our fiducial choice of $T^{\mathrm{ds}\;\mathrm{eq}}=3m_{\mathrm{ds}}$, $m_{\varphi}>2m_{t}(T^{\mathrm{ds}\;\mathrm{eq}})$ for any $m_{\varphi}$ and the results in Fig. 13 remain unaffected. For other choices of equilibration temperatures and DS mass scales, this can be an important effect. For example, if we demand that $T^{\mathrm{ds}\;\mathrm{eq}}=5m_{\mathrm{ds}}$, then the decay channel is open only if $m_{\mathrm{ds}}<100\;\mathrm{keV}$, significantly reducing the mass range over which $\varphi$ decays contribute to DS-photon equilibration.

Cosmologies in which a neutrino- or photon-coupled DS equilibrates after neutrino-photon decoupling demonstrate the need for both CMB- and BBN-based measurements of the expansion rate. Since these two epochs are widely separated in time, non-standard physics can affect one but not the other. This is particularly clear in the photon-coupled case, since the dominant change in the expansion rate can be localized in time between nucleosynthesis and recombination. It is therefore crucial to constrain $\Delta N_{\mathrm{eff}}$ at as many different cosmological times as possible. In specific models, additional constraints may be relevant. This is especially true for light, photon-coupled mediators, which we discuss in Appendix B.

5 Conclusion

Measurements of the light element abundances provide a direct probe of the universe seconds after the Big Bang. The concordance of the predictions of standard Big Bang nucleosynthesis with observations of ${}^{4}{\rm He}$ and D abundances constrains the existence of new physics that contributes significantly to the energy density of the universe at that time. We have investigated how a light dark sector that comes into equilibrium with neutrinos or photons after neutrino-photon decoupling impacts the predictions of ${}^{4}{\rm He}$ and D. This scenario naturally occurs if the processes mediating energy exchange between the dark sector and the SM become important at late times.

The equilibration and eventual decoupling of new particles gives rise to a time-dependent modification of the expansion rate, unlike the standard case of dark radiation. We have investigated several possibilities for this time-dependence, encoded in the temperature evolution of $N_{\mathrm{eff}}$. We considered both model-independent and concrete particle physics-motivated examples. Both approaches illustrate the seldom-mentioned point that, in general, $N_{\mathrm{eff}}$ extracted from the CMB spectrum is different from the one inferred from the light element abundances. Depending on the time-evolution of $N_{\mathrm{eff}}$, either the CMB or BBN can be the more sensitive probe of new physics. For example, if new degrees of freedom equilibrate with neutrinos after neutrino-photon decoupling, then the largest modification of the expansion rate occurs at late times. As a result, measurements of the CMB are typically more constraining for such processes. In fact, CMB-S4 experiments will decisively probe the relativistic neutrino-equilibration scenario, which results in $\Delta N_{\mathrm{eff}}\gtrsim 0.2$ at the time of recombination [20, 21, 65]. Instead, a dark sector that equilibrates with photons after neutrino-photon decoupling may dramatically alter the expansion rate during BBN without significantly impacting it during recombination. In this case, observations of the light element abundances provide the most important test of such new degrees of freedom. Regardless, CMB-S4 will achieve comparable sensitivity to these models in the coming decade. These examples serve to illustrate the fact that measurements of the expansion rate during different epochs are crucial in testing the viability of alternative cosmologies.

In deriving the impact of a modified expansion history on the light element abundances, we have also highlighted the importance of certain nuclear reactions. Different choices for key D-burning rates significantly change both the central value of the prediction and the theoretical uncertainty, leading to significantly different constraints on models of new physics. This ambiguity will be reduced with, e.g., the upcoming measurement of \textrm{D}(p,\gamma)\mbox{{}^{3}{\rm He}} at the LUNA experiment [66].

Our analysis relied on several simplifying assumptions. First, throughout this work we assumed that neutrinos have instantaneously decoupled from the baryon-photon plasma, such that their entropy evolves independently for $T\lesssim 3\;\mathrm{MeV}$. In reality, the decoupling of electroweak interactions is gradual [42] and neutrino-flavor dependent [23]. It would be interesting to investigate how this non-instantaneous neutrino-decoupling affects constraints on models in which dark sector equilibration or decoupling occurs near $T\sim\mathrm{MeV}$. As a first step, one can use the recent simplified neutrino decoupling method from Ref. [67]. Second, in the calculation of primordial abundances, we took the neutrino momentum distribution to be of the equilibrium type. While this is a standard assumption in public BBN codes, new degrees of freedom will in general change the neutrino spectrum and therefore affect certain reaction rates. While the effect on the light nuclei yields from spectral distortions is expected to be small (e.g., in the SM, non-instantaneous decoupling distorts neutrino distributions by $\lesssim 5\%$ [26]), it may become important as the theoretical and observational uncertainties decrease.

As mentioned above, light element abundances typically provide the strongest cosmological constraints on models where new states equilibrate with photons after neutrino-photon decoupling. Such particles can also be produced in stellar environments, leading to bounds from the observed lifetimes of horizontal branch and massive stars, and from the cooling rate of SN1987A. We show in Appendix B that these bounds typically exclude late photon equilibration cosmology. It is therefore important to understand whether this statement is robust, or if there are simple models that avoid (or at least weaken) the stellar constraints, through, e.g., environmental dependence of the dark sector-photon interactions. This is also interesting in the context of the proposed direct detection experiments that seek to discover sub-MeV dark matter [68, 69].

Acknowledgment

We thank Tongyan Lin, Sam McDermott, David Morrissey, Ken Nollett, and Josef Pradler for valuable discussions. AB and SL are supported by the U.S. Department of Energy under Contract No. DE-AC02-76SF00515. AB and NB thank the Kavli Institute of Theoretical Physics (KITP) where part of this work was completed. The research at KITP was supported in part by the National Science Foundation under Grant No. NSF PHY17-48958. NB also thanks TRIUMF for hospitality during the completion of this work. This manuscript has been authored by Fermi Research Alliance, LLC under Contract No. DE-AC02-07CH11359 with the U.S. Department of Energy, Office of Science, Office of High Energy Physics.

Appendix A Collision Terms

In this section, we evaluate the collision terms responsible for the energy transfer between the visible and dark sectors. We focus on the simple example of two-body decays of a DS scalar, $\varphi$: $\varphi\leftrightarrow\nu\nu$ or $\varphi\leftrightarrow\gamma\gamma$, which equilibrate the DS with neutrinos or with photons, respectively. The energy transfer rate in and out of the $\varphi$ bath is then given by

[TABLE]

where $f_{1,2}$ and $f_{3}$ are distributions of the SM states ($\nu$ or $\gamma$) and of $\varphi$; $f_{1,2}$ is a function of $T_{\mathrm{SM}}$ (temperature of $\nu$ or $\gamma$), while $f_{3}$ depends on the DS temperature $T_{\mathrm{ds}}$. In Eq. 40, the factors $(1\mp f_{i})$ are Pauli-blocking or Bose-enhancement factors, with the $-$ sign relevant for $\varphi\to\nu\nu$ and $+$ for $\varphi\to\gamma\gamma$. The two-body decay matrix element $|\mathcal{M}|^{2}$ is related to the decay rate

[TABLE]

where we took the final state particles to be massless. The $d\Pi_{1}d\Pi_{2}$ phase-space integrals in Eq. 40 can be performed by boosting into the $p_{3}$ rest-frame. The remaining integral over $p_{3}$ must be done numerically. The collision term can then be expressed as

[TABLE]

where $r=T_{\mathrm{SM}}/T_{\mathrm{ds}}$, $x=m_{\varphi}/T_{\mathrm{SM}}$. For $\varphi\rightarrow\nu\nu$ we find

[TABLE]

while the $\gamma\gamma$ final state gives the same expressions but with $\cosh\rightarrow\sinh$ and $\operatorname{sech}\rightarrow\operatorname{csch}$. It is clear that since $h_{1}(1,x)=h_{2}(1,x)$, the collision term in Eq. 42 vanishes as $T_{\mathrm{ds}}\rightarrow T_{\mathrm{SM}}$. The $\gamma\gamma$ collision terms also exhibit Bose-enhancement in the form of a logarithmic singularity that is regulated by $m_{\varphi}$, resulting in a larger collision term and, therefore, faster equilibration for a given $\Gamma_{\varphi}$.

Appendix B Constraints on Dark Sector-Photon Equilibration

In this section, we summarize the constraints on models of late DS-photon equilibration. As we showed in Sec. 2, equilibration after neutrino-photon decoupling requires the presence of a light mediator particle with mass below the MeV scale. Such light states are easily produced in supernovae or in red giants and horizontal branch stars, if kinematically allowed. The broad agreement between the observed burst duration of SN1987A, the lifetimes of massive stars, and the corresponding SM predictions place tight constraints on the production of new particles that accelerate the energy loss in these systems [70]. Here, we illustrate these difficulties in the context of an axion-like particle with the interaction

[TABLE]

where $\Lambda$ is a dimensionful scale related to the mass, coupling, and multiplicity of states of the ultraviolet physics that generates this interaction. For example, if this operator is generated from a loop of heavy EM-charged states, then we expect $\Lambda\sim 2\pi M/\alpha_{\text{em}}$, where $\alpha_{\text{em}}$ is the fine-structure constant and $M$ is the mass of the heavy particles.

Thermalization of the visible and dark sectors is accomplished by decays and inverse decays, $\varphi\leftrightarrow\gamma\gamma$. The zero temperature rate is [71]

[TABLE]

The visible and dark sectors equilibrate when $\Gamma_{\varphi}\times(m_{\varphi}/T)\sim H$ (see Sec. 2.1). The equilibration temperature is then

[TABLE]

If we conservatively demand that the equilibration and decoupling take place while SM photon number changing processes are efficient at $T\gtrsim 0.4\;\mathrm{keV}$ [72] we find an approximate requirement on $\Lambda$ 777Equilibration and decoupling may be viable at lower temperatures, since decoupling occurs in thermal equilibrium and naively should not lead to CMB spectral distortions, as long as the DS particles have no chemical potential.

[TABLE]

These couplings are generically in conflict with limits from stellar cooling arguments. For example, for $m_{\varphi}<10\;\mathrm{keV}$, the strongest constraint comes from the shortening of the He burning phase in Horizontal Branch (HB) [73] and massive stars [74], leading to

[TABLE]

For $10\text{ keV}\lesssim m_{\varphi}\lesssim 10\text{ MeV}$, the dominant bound comes from observations of the SN1987A. The agreement between the predicted and measured SM neutrino burst duration implies that [75]

[TABLE]

If $\varphi$ decayed to photons once outside of the SN, even more powerful limits from the non-observation of a $\gamma$-ray excess become applicable [76]. We see that generically the bounds on the $\varphi$ coupling to photons are in conflict with late equilibration. While these limits are extremely powerful, they are not entirely model-independent. For example, the $\varphi-\gamma$ interaction strength can be environment-dependent, and the effective coupling probed by the stellar cooling arguments above would not be the same as the coupling responsible for equilibration in the early universe [77]. Cooling arguments can also be evaded via self-trapping inside of stars or supernovae [78]. Both environmental dependence of the $\varphi-\gamma$ coupling and self-trapping naturally occur if $\varphi$ has sufficiently strong self-interactions, i.e., a non-trivial potential.