Dark Sector Equilibration During Nucleosynthesis

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.


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 farreaching 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 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 eff at late times, while those that interact with photons decrease N eff . Dark sector equilibration and decoupling with the SM can give rise to "pulses" in the modified expansion rate (and therefore N 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, η 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 eff in a model-independent manner, by considering generic templates for N eff variations that can be mapped onto specific models in  In the left panel, N eff does not change during equilibration due to energy conservation; it is only modified when the cold DS particles become non-relativistic and decouple, heating the neutrinos above the SM expectation (indicated by the gray dashed line). In the middle panel, N eff increases when photons equilibrate with a cold DS. When these new particles decouple, they reheat the photon bath, leading to a decrease in N eff . It is also possible that a DS that is initially hotter than the photon bath leads to a decrease in N eff during equilibration and decoupling, as shown in the right panel. For the same number of DS degrees of freedom, the maximum deviation to N eff is larger for photon equilibration. See Sec. 2 for more details.
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 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.

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 ∼ few × MeV [22,23].
We will refer to this process as late equilibration. 1 These models are interesting for several reasons. First, this cosmology is a generic feature of dark sectors with light (sub-MeV) feeblyinteracting 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 Γ eq . At temperatures much greater than the masses of various particles participating in the equilibrating reactions, the rate is schematically of the form Γ eq ∼ α 2 ds T (DM-SM scattering) for DM-SM scattering and Γ eq ∼ α ds m 2 ds /T (DS decays to SM) (2) for decays and inverse-decays of DS force-carriers into the SM, where T is the photon bath temperature, α ds is a coupling constant, and m ds is a characteristic mass-scale in the DS. Equilibration occurs when Γ eq is comparable to the Hubble expansion rate, where M Pl is the Planck mass. Since Γ eq /H increases as the universe cools, the DS will eventually enter equilibrium with the SM, provided that α ds is sufficiently large. Defining T ds eq to be the temperature at which Γ eq ∼ H, we find that the DS and SM enter equilibrium when the temperature drops below T ds eq ∼ α ds M Pl (scattering) T ds eq ∼ α ds m 2 ds M Pl 1/3 (decays) , 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, O(10 −21 ) α ds (m ds /keV) 2 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 Γ 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 massthreshold, 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, ϕ, interacting with neutrinos, where ν 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 ϕ directly coupled to the latter) [21]. The dominant process leading to energy transfer and equilibration between the SM and the DS is the decay ϕ ↔ νν.
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 T ds eq temperature of the photon bath at DS-SM equilibration Table 1: Notation and various temperature scales discussed throughout this work.
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, ϕ, with an interaction [19] L where F is the photon field strength and Λ encodes the scale and couplings of ultraviolet physics that generates this operator. As we show in Appendix B, in minimal models the size of Λ 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 ds , g ds , Γ 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.

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 ∝ √ ρ. As ρ is dominated by relativistic species, after neutrino-photon decoupling it receives contributions from photons, electrons, neutrinos, and potentially the DS: 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 ρ i = (π 2 /30) g i * T 4 i .
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 eff , which is defined as Above, T 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 eff = 3 in the SM. 2 In scenarios with extra light degrees of freedom at the same temperature as the SM neutrinos, N 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 eff no longer directly measures the number of particles in the DS. Equation 8 is consistent with the more common definition, N eff ∝ (ρ ν + ρ ds )/ρ γ , often used in the context of the CMB. Our adopted definition of N 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 eff . For example, in the SM, e ± decoupling increases N eff to ≈ 3.045 [24][25][26]. Throughout this work, we use the photon temperature, T ≡ T γ , as a standard reference for cosmological epochs. For later convenience, we define ξ i as the relative temperature of species i compared to that of the photon bath [27], We use ξ SM ν to denote the value of ξ ν under the assumptions of a standard cosmology and instantaneous decoupling, i.e., ξ SM ν (T ) = 1 for T m e and ξ SM ν (T ) = (4/11) 1/3 0.7 for T m e [28]. Equation 8 can be converted into an explicit expression for N eff , where g ν * = 3 × 2 × (7/8) = 21/4 and g ds * = g ds * (T ) is the (temperature-dependent) effective number of relativistic degrees of freedom in the DS. In general, ξ ν , ξ SM ν , ξ ds , and g 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+ν 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 (ρ SM + ρ ds )a 4 = constant (SM-DS equilibration), where a is the scale factor. The temperature evolution of the SM + 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 SM+DS bath through, e.g., decays or annihilations and eventually decouple. This process occurs in equilibrium and therefore conserves entropy, leading to (s SM + s ds )a 3 = constant (SM-DS decoupling), which can be used to estimate the SM + DS temperature before and after decoupling. Equations 11 and 12 will be used to qualitatively understand the temperature evolution of N 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ρ where i = γ, ν, 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, ϕ, into SM particles (ϕ → SM) with the decay rate, Γ ϕ . The corresponding collision term is then given approximately by where the superscript "eq" denotes an equilibrium distribution with n eq ϕ (T i ) ∼ T 3 i 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 ϕ (T ds ) = n eq ϕ (T ds ). Pauliblocking 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.

Coupling to the Neutrino Bath
We use energy (Eq. 11) and entropy (Eq. 12) conservation to understand the evolution of N eff throughout DS equilibration and decoupling. In the analytic expressions below, we take T m e for simplicity, such that g γ * = g γ = 2 and ξ SM ν = (4/11) 1/3 . Prior to equilibration, the DS contributes to the energy density like "dark radiation" such that where ξ 0 ds is the DS-to-photon temperature ratio prior to equilibration. Although we take ξ 0 ds to be a free parameter, one expects ξ 0 ds < ξ 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 ξ 0 ds /ξ SM ν > 1 we have ∆N eff 0.6. The subsequent dynamics of equilibration and decoupling from neutrinos only increase ∆N eff as described below, worsening the tension. 3 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 ξ 0 ds is lowered.
During DS-SM equilibration, the total comoving energy density is conserved (see Eq. 11). Since N eff is a measure of the total DS + ν energy density, it remains unchanged during this epoch. Therefore, immediately after DS-SM equilibration, the form of N 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-ν bath may drop below a mass threshold of a DS particle. At this point, processes such as annihilations (DS DS → ν ν) or decays (DS → ν ν) 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 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 For ξ 0 ds < 1 and g ds * ≥ 1, we therefore have ∆N eff 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, ξ 0 ds = 1 and ∆N eff is generically in conflict with considerations of BBN and the CMB even for g ds * = 1 [13,32] .
The evolution of N 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 ds * ) at a common mass scale (m ds ); these were calculated using the Boltzmann equation in Eq. 13. For each curve, the decay rate for DS ↔ νν is chosen such that equilibration occurs at T ν ≈ 10 m ds , fixing the initial DS-to-photon temperature ratio to be ξ 0 ds = 0.3. Figure 2 illustrates that before and during equilibration, N eff is very close to 3. The DS mass scale, m ds , sets the decoupling temperature. After DS-ν decoupling, N 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 eff before and after each decoupling event can be estimated from repeated applications of comoving entropy conservation.

Coupling to the Photon Bath
Late equilibration of a DS with photons gives rise to a qualitatively different behavior of N 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 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, η b = (n b − nb)/n γ , which also alters the outcome of primordial nucleosynthesis.
As in the neutrino case, the DS contributes to N eff as dark radiation prior to equilibration such that where, as before, ξ 0 ds is the initial DS-to-photon temperature ratio. DS-γ equilibration fixes ξ ds = 1 and therefore modifies N eff to , where m ds is the common mass scale of the g ds * DS degrees of freedom. For each curve, the initial DS-to-photon temperature ratio is fixed to ξ 0 ds = 0.3 and the SM-DS interaction rate is fixed to ensure equilibration at T ds eq = 10 m ds . Note that equilibration with neutrinos does not change N eff because the total neutrino and DS energy density is conserved. N eff is only significantly modified when the DS degrees of freedom become non-relativistic and transfer their entropy to neutrinos when T ν < m ds .
This equation exhibits two effects that modify N eff upon DS-γ equilibration. The first term in Eq. 18 corresponds to ξ 4 ν in the general form of Eq. 10. Depending on the initial temperature of the DS (encapsulated in ξ 0 ds ), DS-γ equilibration either deposits energy into or draws energy from the γ 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 eff that qualitatively depends on the size of the initial DS temperature ratio, ξ 0 ds , and the number of DS degrees of freedom, g ds * . Equation 17 implies that N eff 3 before equilibration. However, from Eq. 18 we see that after DS-γ recoupling, N eff 3 if g ds * g ds critical , or N eff 3 if g ds * g ds critical , where g ds critical (the critical number of degrees of freedom in the DS) is given by For an initially cold DS (ξ 0 ds < 1), g ds critical is unphysical (g ds critical < 0 < g ds * ) and hence, N 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 eff in a concrete model. We must demand that g ds critical g ds * > 1, which is only possible if ξ 0 ds 1.75. In this case, Eq. 17 implies that N eff 24 before equilibration. Although such non-trivial behavior of N 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  Fig. 2, but for a light DS that equilibrates with photons after neutrino-photon decoupling. In addition to N eff evolution (left panel), we also show the evolution of η b (right panel). Equilibration of the DS with photons lowers the photon temperature and introduces additional relativistic degrees of freedom into the plasma. Both effects lead to a large increase of N eff at recoupling, while the lower photon temperature increases η b = (n b − nb)/n γ . When the DS particles become non-relativistic (T m ds ), the DS entropy is transferred back to the photons, reheating them slightly compared to neutrinos, resulting in N eff 3. In the right panel, we fixed the initial value of η b such that the baryon abundance matches the measured CMB value (indicated by the black dotted line) at late times. The light-gray solid line shows the evolution of η b in the SM, where the initial decrease at T γ ∼ m e is due to entropy injection from e + e − annihilations.
such as annihilations or decays transfer entropy to the photons, heating them relative to neutrinos and decreasing ξ ν . As a result, N eff decreases to which is strictly less than the SM expectation of N eff ≈ 3. The evolution of N eff throughout DS-γ equilibration and decoupling is illustrated schematically in the middle and right panels of In addition to its effect on the expansion rate, DS-γ equilibration modifies the predictions of BBN by changing the evolution of η b . Suppose that before DS-γ equilibration and e + e − annihilations, the baryon-to-photon ratio is given by η 0 b . Since (n b − nb)/s ν is constant after neutrinos have decoupled, η b at a later time is given by This parameter is independently determined by observations of the CMB to be η CMB b ≈ 6.1 × 10 −10 at the time of recombination [32], fixing its value at late times. For a given cosmological history, which determines ξ ν (T ), the initial value of η 0 b must be chosen such that the late-time value (Eq. 21) is compatible with CMB measurements, i.e., where T CMB ∼ 0.1 eV. In the case of DS-γ equilibration, the evolution of the neutrino-to-photon temperature ratio, ξ ν , throughout equilibration and decoupling is given by Therefore, ξ ν (T 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 ds * and m ds .

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, 3 He, 4 He, and 7 Li 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 eff and η b on 4 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 He and 7 Li in obtaining bounds on light dark sectors. There is no consensus on the measured 3 He primordial abundance due to limited observational regions and uncertainties in theoretical modeling [33,34]. The measured 7 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 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 He and D yields inconsistent with observations. Throughout this work, we therefore assume that the resolution of the lithium problem is astrophysical in nature.

Dependence on N eff and η b
In Sec. 2.2, we examined modifications to N eff and η 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 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 + ↔ pν e . At temperatures smaller than the neutron-proton mass-splitting (T ∆m np ≈ 1.3 MeV), the neutron-to-proton number density ratio, n/p, becomes smaller than unity, i.e., n/p ≈ exp(−∆m np /T ) 1. Eventually, neutron-proton conversion processes freezeout at T ≈ 0.8 MeV at which point n/p ≈ exp(−1.3/0.8) ≈ 1/5. 4 The decays of free neutrons further reduces the neutron density to n/p ≈ 1/7 before the onset of nucleosynthesis. Nearly all of these neutrons are eventually bound into 4 He because of its large binding energy and the corresponding enhanced equilibrium number density. The resulting 4 He mass-fraction, Y p , is therefore given by Y p ≈ 4 n He n p + n n ≈ 4 (n n /2) n p + n n = 2 (n/p) 1 + (n/p) Although the per-nucleon binding energies of light nuclei are O(MeV), nucleosynthesis does not occur until much smaller temperatures near T ∼ 100 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 He was in chemical equilibrium with n and p its abundance would be negligible until T ∼ 300 keV [28]. The synthesis of these heavier nuclei is delayed even further because their production rates, e.g., via D D → p 3 H (the first step in making 4 He), depend on the deuterium abundance. These reaction rates are slow compared to Hubble expansion until n D /(n p + n n ) ∼ 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, where B D ∼ 2 MeV is the deuterium binding energy. This ratio is initially 1 because of the large entropy of the plasma (small baryon-per-photon number, η b ∼ 10 −10 ). This abundance becomes large enough to support the production of heavier nuclei only when the temperature falls below ∼ 100 keV. At this point, nucleosynthesis commences, and most of the deuterium is eventually processed into 4 He. However, deuterium-burning processes (such as D p → 3 He γ) remain effective until they freeze out at temperatures of T ∼ few × 10 keV and the deuterium abundance has been depleted down to trace amounts of D/H ∼ 10 −5 − 10 −4 . details, see, e.g. Ref. [42].  We show how theoretical predictions of the neutron-to-proton ratio (n/p), the 4 He massfraction (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 eff and η b , respectively. Values of N eff that are greater or smaller than the SM prediction of N eff ≈ 3, correspond to an increased or decreased expansion rate during nucleosynthesis. For N eff 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 He yield (see Eq. 24). Also for N eff 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 eff 3.
The 4 He mass-fraction is only logarithmically sensitive to the baryon density, η b . This effect arises due to how η b controls the deuterium bottleneck. Defining the temperature at which nucleosynthesis commences as T BBN ∼ O(100) keV through the criterion n D /(n p + n n ) ∼ O(10 −4 ) [41], we see from Eq. 25 that where in the second equality, we have considered small variations of η b around the CMBpreferred value, i.e., η b ∼ O(10 −10 ) + δη b . Hence, for much larger η b , nucleosynthesis begins at earlier times, and less unbound neutrons decay before ending up in 4 He. As a result, Y p increases for larger η 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 BBN . Unlike 4 He, deuterium is much more sensitive to deviations in the baryon density, since η b directly controls the strength of D-burning rates. In particular, larger η b corresponds to more efficient D-burning and a smaller trace abundance at late times.

Procedure
The primordial 4 He and D densities are inferred from direct observations of various astrophysical sites. For instance, the primordial 4 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 Lyα 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, We calculate 4 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 eff (T ) and η b (T ). 5 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 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 → D γ and 3 He n → p 3 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 ∼ 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 ↔ 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 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 He and D yields by lognormal-sampling nuclear rates [59,60] and computing the yields 10 4 times for each value of N eff and η b . The resulting yield distribution is well-Rate Reference p(n, γ)D [56] 3 He(n, p) 3 H  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.
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.
where we have taken the SM expectation of N eff ≈ 3, in the first line we have fixed the latetime 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, η CMB b = (6.10 ± 0.04) × 10 −10 [13,61]. In the first and second lines of Eq. 28, we find a correlation coefficient of ≈ −0.22 and ≈ −0.26, respectively. The anti-correlation is driven primarily by the sizable spread in D D → n 3 He and D D → p 3 H (larger values of these rates reduce D/H and enhance 4 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 eff ≈ 3, where η b is fixed as in Eq. 28. We find a correlation coefficient of ≈ 0.006 and ≈ −0.10 in the first and second lines, respectively; the correlation is mostly due to the η 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 He abundance is negligible compared to the observational one. This is because the 4 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]: 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 τ n is closer to that inferred from beam-based measurements (∼ 888 s), the corresponding shift to Y p and D/H would be at the sub-percent level.
We find that the fractional theoretical uncertainties, σ(X)/X for X = Y p and D/H, are to a good approximation independent of the functional form of N eff (T ) and η 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, Γ i , is encompassed by a set of constants, α i (the logarithmic derivatives of the yields with respect to the rates), such that In the nucleosynthesis Boltzmann system, these rates enter as Γ i /H, which suggests that the effect of variations of H (via N eff ) can be likewise linearized, at least for small deformations of H, leading to X ∝ N eff α , for some other constant α . The resulting set of logarithmic derivatives, {α i , α }, is presented in Ref. [13]. A similar argument can be made for modifications to η b . Using this linearized form of the yields, it is simple to show that fractional errors from nuclear rate uncertainties, σ(X)/X, are constant as N eff and η 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 eff and η b , they are more difficult to make for arbitrary time-variations of these quantities. Therefore, we have explicitly checked the constancy of σ(Y p )/Y p and σ(D/H)/(D/H) for timedependent variations of N eff and η b , as considered in the following sections, using the Monte Carlo approach described above. This justifies our use of σ(Y p )/Y p and σ(D/H)/(D/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 eff and η b and in concrete particle physics models, respectively.

Standard Constraints on N eff and η b
Before we discuss our main results in Sec. 4, we first present a simple estimate for the standard bounds on N eff and η 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 He and D yields, (Y p ) theory and (D/H) theory , respectively. We then compare these predictions to the observed  [58]). The same analysis using the rates of Ref. [55] is shown in dark green (2σ). Along the dot-dashed and dotted gray contours, |∆Y p | = 2 and |∆D/H| = 2, respectively, using the rates of Ref. [58] (see Eq. 33). We also compare to regions 2σ-favored by recent Planck measurements of the cosmic microwave background (cyan) [32]. abundances in Eq. 27 through the χ 2 test statistic, where we have defined This gives χ 2 ≈ 5.7, assuming a standard cosmological history.
The above definition of χ 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 ∆Y p ∆D/H or ∆Y p ∆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 D/H even when adopting the NACRE-II rates. Figure 6 illustrates how modifications to the baryon density (η b ) and the expansion rate (N eff ) are constrained from measurements of the primordial densities of helium-4 and deuterium.
The parameters, ∆η b and ∆N eff , correspond to a shift of the late-time (CMB era) value of η b and a time-independent modification of N eff away from the SM expectation, respectively. The evolution of η b (T ) during nucleosynthesis for a given ∆η b is evaluated using entropy conservation, as in Eqs. 21 and Eq. 22, assuming the standard form of ξ ν . The SM is defined by N eff ≈ 3, with the baryon-to-photon ratio fixed at late times to the CMB-preferred value, η b ≈ 6 × 10 −10 .
In our calculations of the primordial nuclei abundances, the best-fit point is defined as the value of ∆N eff and ∆η b that minimizes χ 2 , as defined in Eq. 32. In Fig. 6, this occurs at (∆N eff , 10 10 ∆η b ) = (−0.04, −0.2) with a χ 2 of χ 2 min = 0.04. In Fig. 6, we show 1σ (dark blue) and 2σ (light blue) regions around this best-fit point, corresponding to ∆χ 2 ≡ χ 2 − χ 2 min = 2.30 and 6.18, respectively, using the nuclear rates of Ref. [58]. Standard cosmology is consistent with the observed abundances within 2σ. Along the dot-dashed and dotted gray contours of Fig. 6, |∆Y p | = 2 and |∆D/H| = 2, respectively. As discussed in Sec. 3.1, Y p is dominantly sensitive to N eff with only a logarithmic dependence on η b . On the other hand, the effect of η b on the predicted deuterium abundance is much larger, but is degenerate with N 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 η b .
Compared to Ref. [13], our baseline analysis shows a mild preference for smaller values of η 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σ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].

Equilibration and Decoupling during Nucleosynthesis
In Sec. 3, we reviewed how the expansion rate (N eff ) and the baryonic density (η 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 eff and η 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 eff and η 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 eff and η 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 ∆N eff (T ) (where T is the temperature of the photon bath) and fixing ∆η b = 0. In Sec. 4.2, we consider a few concrete models that predict non-standard temperature evolution of ∆N eff (T ) and ∆η b (T ) and discuss how the predicted nuclear abundances are modified in relevant regions Figure 7: Schematic model-independent temperature evolutions for N eff that we consider in Sec. 4.1. The left and right panels correspond to the left and right panels of Fig. 10. In the left panel, N eff tracks the Standard Model expectation at early times (∆N eff = 0). Later, a deviation to N eff (∆N eff = 0) occurs below some critical temperature, T on . A similar scenario is shown in the right panel, in which a deviation turns off after the universe cools below the critical temperature, T off . In Fig. 11, we also consider the possibility that a deviation to N eff turns on and subsequently turns off during the epoch of primordial nucleosynthesis.
of parameter space.

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 eff (T ). The first possibility that we consider is that N eff tracks the SM expectation of N eff ≈ 3 until some later temperature, T on , at which point a deviation turns on, i.e., ∆N eff = 0 for T T 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 where the temperature-dependence is encapsulated in the Heaviside step-function, Θ. We also consider the possibility that ∆N eff = 0 at the beginning of nucleosynthesis, but that this deviation turns off at some critical temperature, T off , i.e., 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 ∆N eff (T ). We therefore also investigate the generalization of Eqs. 35 and 36, This corresponds to a sudden pulse-like modification to N eff within the temperature range of T off T T on .  Figures 8 and 9 show the fractional shift (compared to ∆N eff = 0) to the helium and deuterium abundances for these various temperature evolutions of ∆N eff (T ), fixing the amplitude to ∆N eff = 0.5. The final abundance of 4 He is modified in the standard way (see Sec. 3.1) if ∆N eff = 0 during neutronproton freeze-out (100 keV T MeV), while the deuterium abundance is most sensitive to the expansion rate near the freeze-out of D-burning process (50 keV T 100 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 ∼ 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 ∆N eff − T on and ∆N eff − T off plane, respectively, assuming a step-like temperature evolution of ∆N eff (T ), as parametrized by Eqs. 35 and 36. Also shown as blue and green contours are parameters for which |∆Y p | = 3 and |∆D/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 ∆χ 2 6.18. Throughout, we have fixed the baryon density, η 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 keV T MeV, then our calculated bound on N eff is similar to the standard constraint presented in Fig. 6, i.e., |∆N eff | 0.5. However, if modifications to the expansion rate occur when T 50 keV or T MeV, then deviations as large as |∆N eff | 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 eff  Fig. 8, but for a pulse-like modification to N eff (see Eq. 37). Its impact on the yields of 4 He and D as a function of the transition temperature, T on , is shown for T on /T off = 3 and 10 (solid and dashed colored lines, respectively) and for a fixed pulse amplitude of ∆N eff = 0.5.
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 ∆N eff = 0 only for T few × O(100) keV are predominantly constrained by measurements of the deuterium abundance. Also apparent in Fig. 10 is the asymmetrical importance of Y p and D/H for constraining ∆N eff > 0 or ∆N eff < 0. This can be understood from the fact that the SM (∆N eff = ∆η 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 ∆N eff < 0 comes from measurements of D/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 ∆N eff are shown in Fig. 11 (see Eq. 37). In each panel of Fig. 11, we vary the fractional width of the N eff pulse, (T on − T off )/T on , as well as the temperature at which the pulse turns on, T on , for different choices of the pulse height, ∆N eff = ±1, ±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 on , T off , and N eff , we exclude parameters at the 2σ level if ∆χ 2 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 keV T MeV. As in Fig. 10, measurements of the deuterium abundance are typically more constraining for ∆N eff < 0. On the other hand, for ∆N eff > 0, the relative importance of helium and deuterium strongly depends on when the deviation to N 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 on − T off )/T on (corresponding to narrower pulses) visible in the left column of Fig. 11 can be understood by comparison with Fig. 9.

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, ϕ: ϕ ↔ νν or ϕ ↔ γγ. 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 eff (T ) and η 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 ds , and DS degrees of freedom, g 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 ds ) and the DS effective number of relativistic degrees of freedom (g ds * ) for an initial temperature ratio of ξ 0 ds = 0.3. The 2σ-excluded regions (solid gray) correspond to ∆χ 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 eff (see Eq. 35) occurring at T on ∼ m 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 ds MeV), T on MeV and N eff is modified during every key epoch of nucleosynthesis, which is very strongly constrained for g ds * 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 ∼ 100 keV are typically unconstrained from considerations of BBN, since in this case N 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  Figure 11: Model-independent constraints on the pulse-like temperature evolution of N eff (see Eq. 37). Shown in gray are regions of parameter space that are inconsistent (within 2σ) with observations of the primordial helium-4 and deuterium abundances (calculated using the rates of Ref. [58]). Along the blue dot-dashed and green dotted contours, |∆Y p | = 3 and |∆D/H| = 3, respectively (see Eq. 33). Significant deviations to the expansion rate are allowed to occur during nucleosynthesis, provided that this occurs after neutron-proton freeze-out and before deuterium burning (100 keV T MeV). See Sec. 4.1 for a more detailed discussion of the various temperature-dependent features of these exclusions. the bound on g ds * in Fig. 12 by using T on ∼ m ds /5 and Eq. 16 to translate ∆N eff into g ds * .
The horizontal dotted line in Fig. 12 shows the 2σ upper bound on g ds * derived from Eq. 16 and Planck measurements of N eff at recombination 6 [32]: 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 eff . In contrast, bounds derived from primordial nucleosynthesis depend on precisely when modifications to N 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 eff is SM-like throughout the most important epochs of nucleosynthesis. Hence, bounds derived from BBN are weakened for m ds 100 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 eff ≈ 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 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 η b is fixed to the CMB-preferred value. This parameter point, however, is robustly ruled out by Planck [32], since it corresponds to ∆N eff ∼ 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 ds eq = 3 m 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 3) with the photon bath (assuming that this occurs after neutrino-photon decoupling), as a function of the dark sector degrees of freedom (g ds * ) and the dark sector mass scale (m ds ). We have fixed the equilibration temperature to be T ds eq ≈ 3m ds , such that these dark sector degrees of freedom are relativistic at this time. The 2σ-excluded region is shown in gray. The baryon density, η b , is determined by requiring consistency with the late-time CMB measurement (see Eq. 21). The resulting constraints are shown for two parametrizations of D-burning rates: those from the NACRE-II [55] compilation (dashed line) and Coc et al. (2015) [58] (solid line). 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σ upper bound on g ds * that is derived from Eq. 20 and Planck measurements of N eff at recombination [32]: 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 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, η b . If we were to ignore this effect on η b , then BBN would only constrain a finite interval of m ds (for a fixed value of T ds eq /m 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 ds , respectively), then the predicted yields are not modified relative to the SM case. However, once we account for modifications to η b , the resulting bounds extend to much lower DS masses. The baryon-to-photon ratio is independently measured at recombination, and hence we fix η 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, η 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 η 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 ds * in Fig. 13 can be understood by comparing the maximum deviation of N eff during photon equilibration (Eq. 18) to the constraint on standard constant shifts of ∆N eff 0.5 in Fig. 6: solving for g ds * , one finds that the expected bound should lie near g ds * ∼ 0.1, which is borne out in the complete calculation of Fig. 13. The observed abundances of light nuclei constrain g 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 ds eq < 3m ds ), then the effective value of g 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 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 ds /4 T ds eq < 3m ds . This is because for T m ds /4, the DS-SM equilibration rate no longer increases relative to the Hubble parameter with the expansion of the universe, while for T 3m ds , g ds * ≥ 1 and equilibration is relativistic. Furthermore, even for scenarios in which T ds eq ∼ m ds /4 (corresponding to g ds * 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 m e ), the effective mass of transverse photon excitations is m 2 t ≈ 4παT 2 /9, which falls quickly to m 2 t ≈ 4παn e /m e as T drops below m e [64]. In certain regions of parameter space, this mass kinematically forbids the equilibrating processes ϕ ↔ γγ if m ds = m ϕ < 2 m t (T ) [19]. However, for our fiducial choice of T ds eq = 3m ds , m ϕ > 2m t (T ds eq ) for any m ϕ 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 ds eq = 5m ds , then the decay channel is open only if m ds < 100 keV, significantly reducing the mass range over which ϕ 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 ∆N 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.

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 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 neutrinophoton decoupling impacts the predictions of 4 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 in-vestigated several possibilities for this time-dependence, encoded in the temperature evolution of N eff . We considered both model-independent and concrete particle physics-motivated examples. Both approaches illustrate the seldom-mentioned point that, in general, N eff extracted from the CMB spectrum is different from the one inferred from the light element abundances. Depending on the time-evolution of N 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 ∆N eff 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 D(p, γ) 3 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 3 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 ∼ 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 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].

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 where Λ 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 Λ ∼ 2πM/α em , where α 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, ϕ ↔ γγ. The zero temperature rate is [71] The visible and dark sectors equilibrate when Γ ϕ ×(m ϕ /T ) ∼ H (see Sec. 2.1). The equilibration temperature is then If we conservatively demand that the equilibration and decoupling take place while SM photon number changing processes are efficient at T 0.4 keV [72] we find an approximate requirement on Λ 7 Λ < 5 × 10 5 GeV × m ϕ keV 2 (equilibration at T 0.4 keV).
These couplings are generically in conflict with limits from stellar cooling arguments. For example, for m ϕ < 10 keV, the strongest constraint comes from the shortening of the He burning phase in Horizontal Branch (HB) [73] and massive stars [74], leading to Λ 10 10 GeV (stellar cooling).
For 10 keV m ϕ 10 MeV, the dominant bound comes from observations of the SN1987A. The agreement between the predicted and measured SM neutrino burst duration implies that [75] Λ 10 9 GeV (SN1987A burst duration).
If ϕ decayed to photons once outside of the SN, even more powerful limits from the nonobservation of a γ-ray excess become applicable [76]. We see that generically the bounds on the ϕ coupling to photons are in conflict with late equilibration. While these limits are extremely powerful, they are not entirely model-independent. For example, the ϕ − γ 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 ϕ − γ coupling and self-trapping naturally occur if ϕ has sufficiently strong self-interactions, i.e., a non-trivial potential.