BBN constraints on universally-coupled ultralight scalar dark matter

Ultralight scalar dark matter can interact with all massive Standard Model particles through a universal coupling. Such a coupling modifies the Standard Model particle masses and affects the dynamics of Big Bang Nucleosynthesis. We model the cosmological evolution of the dark matter, taking into account the modifications of the scalar mass by the environment as well as the full dynamics of Big Bang Nucleosynthesis. We find that precision measurements of the helium-4 abundance set stringent constraints on the available parameter space, and that these constraints are strongly affected by both the dark matter environmental mass and the dynamics of the neutron freeze-out. Furthermore, we perform the analysis in both the Einstein and Jordan frames, the latter of which allows us to implement the model into numerical Big Bang Nucleosynthesis codes and analyze additional light elements. The numerical analysis shows that the constraint from helium-4 dominates over deuterium, and that the effect on lithium is insufficient to solve the lithium problem. Comparing to several other probes, we find that Big Bang Nucleosynthesis sets the strongest constraints for the majority of the parameter space.


Introduction
The nature of dark matter (DM) is still unknown, despite many decades of dedicated research. Proposed DM candidates span a mass range from ∼ 10 −22 eV to the Planck mass for elementary particles and further up to several solar masses for composite objects, such as primordial black holes. A well-motivated class of DM candidates is found at the lower end of this mass range. Here, the DM is represented by an ultra-light -of mass well below 1 eV -boson with extremely weak couplings to the Standard Model (SM) fields. Representatives of this class are the axion, which was proposed as a solution to the strong CP problem of QCD [1][2][3][4], string theory moduli [5,6] and dark photons [7,8]. Such an ultra-light DM (ULDM) candidate is characterized by large particle occupation numbers in phase space and can be described as a classical field.
The weakness of the coupling between ULDM and the SM makes detecting ULDM a challenging prospect. Many efforts have been undertaken to identify observations and experiments that are sensitive to ULDM, ranging from cosmological and astrophysical probes [9] to terrestrial experiments [10]. One promising avenue for probing ULDM is to look for its imprints in the abundance of the primordial elements. The formation of the primordial elements, such as helium, deuterium, and lithium, is governed by Big-Bang Nucleosynthesis (BBN), which describes an epoch of primordial nuclear transformations in the expanding Universe. Measurements of the primordial abundances of helium-4 ( 4 He) and deuterium (D) are very precise and agree well with the theoretical predictions within the SM, whereas the inferred abundance of lithium is smaller by about a factor of 3 than the theoretical one [11,12]. The latter discrepancy may be due to a systematics in the observations or some kind of new physics. Putting this lithium problem aside, even a mild modification to the conditions in the early Universe would generically lead to observable deviations in the abundances of 4 He and D making BBN a powerful probe for beyond the SM physics [13,14]. In particular, it has been extensively used to constrain the variations of fundamental constants [15][16][17][18][19] and modified theories of gravity [20][21][22][23]; a summary of many of these studies can be found in [24].
In this work, we investigate the constraints imposed by BBN on the interaction of a real scalar ULDM field φ with the SM. We focus on the case of a universal coupling, meaning that φ interacts with the SM fields through an effective metric, g µν = g µν 1 + 2α(φ) . (1.1) This case is interesting since it preserves the weak equivalence principle (WEP) in the SM sector and thereby avoids many laboratory bounds. This type of coupling naturally appears in scalar-tensor theories of gravity [25]. Further, we will assume that the function α(φ) is even, so that its Taylor expansion starts with the quadratic term, where Λ is a scale of the underlying UV-physics. At the fundamental level, this property can be enforced by requiring the symmetry of the theory under the reflection φ → −φ.
As a consequence, the linear term in the coupling is absent and this scenario avoids the bounds coming from the tests of r −2 fall-off of the gravitational force [26] as well as from the bending of light by the Sun [27]. Consequently, this type of coupling is only weakly constrained by the present-day observations. The prospects to probe it using binary pulsar timing have been discussed in [28,29]. The amplitude of φ is diluted by the expansion of the Universe and hence had a much larger value at earlier epochs. Correspondingly, the coupling α(φ) was also larger at early times than it is now, with a more pronounced effect on the dynamics of the SM sector. This makes BBN a sensitive probe of the above scenario.
Setups similar to ours have been previously considered in the literature. Refs. [20,21] studied the effects on the BBN in scalar-tensor gravity theories. However, the scalar field in these works was assumed to be massless and hence could not play the role of DM. The non-zero mass was included in [30] and more recently in [31]. Our work complements these previous studies in three key aspects: 1. We fully take into account the back-reaction from the SM particles on DM, which modifies the DM mass and leaves a substantial impact on its evolution, 2. We use a kinetic description of neutron freeze-out, instead of the instantaneous approximation, which we find to be inadequate, 3. We perform a numerical BBN analysis using the AlterBBN package [32], which allows us to calculate the primordial abundances for all the light elements. This third point is facilitated by considering the theory in the Jordan frame, where the effects of φ reduce to a change in the expansion rate of the Universe. We also clarify the relation between the calculations in the Jordan and Einstein frames.
The paper is organized as follows: We define the model in section 2, where we describe our assumptions and give a preliminary discussion of the BBN sensitivity. In section 3, we quantify the cosmological evolution of the scalar field starting with the assumption that φ makes up the DM today. The detailed analysis of the BBN constraints is performed in section 4. We first analytically estimate the effects of φ on the 4 He abundance in section 4.1 and then present our numerical constraints in section 4.2. These are compared to other existing bounds on the model in section 5. We conclude in section 6.

A real scalar with universal couplings
We begin by adding to the SM a real scalar field φ with the Lagrangian We assume that any possible self-interaction of φ can be neglected, whereas its couplings to the SM fields preserve the WEP. This implies that φ is coupled to the SM through the effective metric (eq. 1.1). We will further assume the coupling to be weak, |α(φ)| 1 at the values of the field that we are going to consider. This condition sets the domain of validity of all the equations derived below.
As the SM Lagrangian contains fermions, one needs to specify the form of the effective vielbein, rather than just the effective metric,ē a µ = e a µ 1 + α(φ) . In this way we arrive at the total Lagrangian of our setup, where M pl is the Planck mass, R is the Ricci scalar, and L SM is the SM Lagrangian and we have collectively denoted the SM fields with ψ. Note that we have assumed the standard Einstein-Hilbert term for the gravitational action. Expanding to the linear order in α we can write the leading interaction term as where Θ SM [ψ] is the trace of the energy-momentum tensor (EMT) of the SM sector. Let us discuss the contributions of various fields to eq. 2.3. For fundamental fermions (quarks and leptons), using their equations of motion, we obtain, Similarly, for massive vector fields (W and Z bosons) we have, In both cases the interaction amounts to a rescaling of the particle mass by a φ-dependent factor, The situation is more subtle for composite particles, such as hadrons, whose mass is dominated by non-perturbative QCD contributions. Nevertheless, we presently argue that the mass scaling (eq. 2.6) applies to them as well. Indeed, the non-perturbative contributions to the hadron masses are determined by the QCD scale Λ QCD . The latter can be related to the physics at high energies using the renormalization group running of the strong coupling constant α s . Including the mass thresholds from charm, bottom and top quarks, one has [17], where m c,b,t are the masses of the charm, bottom, and top quarks, respectively, and M UV > m t is the high-energy scale where the strong coupling is normalized to a given value. If the φ interactions are to preserve the WEP, the scale M UV must vary according to with α s M UV (1 + α) held fixed. This will be the case, for example, in Grand Unified Theories where the cutoff scale M UV corresponds to a physical particle mass. Then Λ QCD scales in the same way as the fundamental particle masses, implying that the φ-dependence of all hadronic quantities is dictated by their mass dimensions. This applies not only to the hadron masses, but also to nuclear binding energies, decay widths, and cross sections. See appendix A for more discussion of this point. In principle, one could envision scenarios where the UV scale does not obey eq. 2.8. However, this will inevitably lead to a violation of WEP in the low-energy physics and would make the predictions of the theory dependent on the details of the UV completion. In this paper we avoid such complications by restricting to the WEP preserving case.
We can get another perspective on the effects due to φ by transforming the metric from the Einstein to Jordan frame, g µν →ḡ µν . Keeping the leading terms in α we obtain from eq. 2.2, Note that the second term in the first line, though quadratic in ∂ µ α, is enhanced by the Planck mass and cannot be neglected in general. We observe that φ has decoupled from the SM, which now interacts covariantly only with the metricḡ µν . Quantum corrections within the matter sector do not change this covariant form and thus do not spoil WEP. The latter can still be violated by quantum corrections involving metric perturbations in the loops. These are, however, suppressed by the inverse Planck mass. The non-minimal coupling of φ to the Jordan metricḡ µν implies that we are dealing with a variant of scalar-tensor gravity. The field φ affects BBN by modifying the expansion rate of the universe at the BBN epoch. The modification stems from the φ-dependence of the Planck mass, as well as the EMT of φ. Focusing for simplicity on the first effect -the variation of the Planck mass -we can get a rough idea of the BBN reach in constraining the model parameters. The difference between the Planck mass at the time of BBN and today is bounded at the level of 3% [19]. The coupling α(φ) nowadays is very small due to the dilution of the amplitude of φ by the expansion of the universe. Thus, we expect BBN to exclude the parameter region where the coupling |α(φ)| exceeds roughly 0.015 at the BBN epoch. On the other hand, BBN will remain insensitive to the presence of φ if |α(φ)| 0.01 at t BBN . We will see below that the presence of the φ EMT, as well as its time evolution, make the true story somewhat more complicated. In particular, the constraining power can be reduced in some parameter regions due to cancellation between several competing effects.
We now further specify the model by imposing the condition that the coupling α is an even function of φ, so that its Taylor expansion starts with a quadratic term (see eq. 1.2). The higher-order terms are assumed to be negligible as long as φ 2 /Λ 2 1. This implies the validity of the form (eq. 1.2) all the way back through BBN for the relevant range of parameters. This has an important implication for the dynamics of φ in the early Universe. From the expression (eq. 2.3) for the interaction Lagrangian in the Einstein frame we see that the presence of SM matter with energy density ρ SM and pressure p SM induces a contribution into the mass term of φ leading to an effective time-dependent mass, In the next section we study the evolution of the ULDM field taking this effective mass into account. Prior to concluding this section, let us comment on the issue of quantum corrections to the φ-mass and its self-interaction due to its coupling to the SM. Working in the Einstein frame and assuming, as before, the preservation of WEP we obtain that the result of integrating out the SM fields should have the form, where ρ vac SM is the SM vacuum energy density. It is well-known that a naive estimate of loop contributions would yield ρ vac SM ∼ m 4 t . Such a large contribution would completely destroy our scenario. However, we also know that the total vacuum energy density ρ vac tot is very small, implying a delicate cancellation between different contributions into it, including the bare (unrenormalized) value. The mechanism ensuring this cancellation is still unknown, which constitutes the famous cosmological constant problem. We do not attempt to add anything to its solution and just speculate that the cancellation of vacuum energy can happen separately within the SM sector (or its extension universally coupled to φ) implying ρ vac SM ρ vac tot . In this case the effect of quantum corrections (eq. 2.13) is completely negligible.

Evolution of dark matter
To understand the effect of DM on BBN, we calculate the cosmological evolution of the DM field φ from the present epoch back to the time when the temperature of the Universe was T ∼ MeV. This epoch corresponds to the freeze-out of the weak interactions and sets the neutron abundance, which is one of the key ingredients of BBN. The energy density of φ is subdominant until the epoch of matter domination. In addition, we will neglect in this section the φ-induced changes in the evolution of the SM fields. These changes would lead only to O(α 2 ) corrections in the eventual BBN analysis, whereas the leading effects are linear in α. Thus, working in the Einstein frame, we can treat φ as a probe field embedded into a Friedmann-Robertson-Walker (FRW) universe with the scale factor a(t) obeying the standard expansion history. Assuming that the field is spatially homogeneous, we obtain its equation of motion,φ where˙denotes the derivative with respect to time, H ≡ȧ/a is the Hubble rate, and m 2 φ,eff is given by eq. (2.12).
Let us discuss the trace of the SM EMT, which enters into the effective mass. Relativistic particle species (photons, neutrinos) do not contribute into it 1 . For most of the time the SM EMT is dominated by the non-relativistic baryonic matter and reads, where ρ b (t) is the time dependent energy density of baryons, H 0 and Ω b are the present-day Hubble constant and the baryon density fraction, respectively, and we have normalized the scale factor to be unity today, a 0 = 1. However, around the epoch of BBN, Θ SM receives a large additional contribution from the electron-positron plasma. Due to the numerical coincidence between the electron mass and the BBN temperature, electrons and positrons become non-relativistic during BBN, while their number density still greatly exceeds that of protons and neutrons until e + annihilate with e − at somewhat lower temperatures. Using the standard thermodynamic expressions for the energy density and pressure of a Fermi gas we obtain the trace of e + e − EMT, where m e is the electron mass and we have neglected the electron chemical potential. The factor 4 in front counts the number of spin degrees of freedom. The total SM EMT trace is given by and is shown in fig. 1 (left panel) as a function of the scale factor 2 . Note that the EMT trace never exceeds 10% of the SM energy density during the BBN epoch [35], The evolution of φ has three distinct regimes characterized by the dominant term in the equation of motion eq. 3.1: • Induced mass dominance (I): Depending on the model parameters, the field goes through these regimes in various sequences, which may be quite complicated, as shown in fig. 1 (right panel). This is due to the non-trivial time-dependence of the induced mass. 1 We neglect the trace anomaly which is relevant only at temperatures above several MeV when electrons are relativistic. 2 To convert the temperature into the scale factor we used the formula [33] where T0 is the present-day Cosmic Microwave Background temperature and g * S(T ) is the effective number of relativistic degrees of freedom in the entropy density, which we obtain from the microOMEGAs package [34]. In the regime (H), the field is frozen at a constant value. On the other hand, in the (B) regime the field oscillates, where ϕ is a constant phase and the amplitude decreases with time as Φ(t) ∝ a −3/2 (t). In this regime the field φ behaves as DM [9]. In particular, its energy density scales inversely proportional to a 3 , Equating this to the measured average DM density today ρ DM,0 = 1.26 × 10 −6 GeV/cm 3 determines the present-day amplitude To reproduce the success of the ΛCDM cosmology, the field must be in the regime (B) throughout the matter-dominated stage of the history of the universe, which puts the constraint on the model parameters, where a eq 10 −4 is the scale factor at the epoch of matter-radiation equality. The behavior of the field φ in the regime (I) differs qualitatively, depending on the sign of the coupling (1.2).
Negative coupling In this case the induced mass term is tachyonic and leads to an exponential growth of the field, with the approximate solution, In principle, this can serve as a mechanism for DM production. Note, however, that it requires a fine-tuning of the initial φ-value in order not to over-produce DM. We do not discuss a possible origin of such tuning. Instead, we adopt a phenomenological approach and, in order to see what constraints BBN imposes on this type of coupling, evolve the field backward in time from today by matching the oscillations in the (B) regime to the exponential growth (3.10) in the (I) regime. Positive coupling In this case the solution in the regime (I) is oscillating, similar to the regime (B). As the oscillations are much faster than the expansion of the Universe, we can obtain a unified description in these two regimes using the WKB-type expansion. This yields, where the amplitude scales as, Using this scaling and the known amplitude Φ 0 today, we determine the amplitude back in time until the field enters into the (H) regime, where it freezes at a constant. To patch the solutions in (H) and (I)/(B) regimes, we solve eq. (3.1) numerically. As already noted, the field φ may pass between different regimes several times. In that case, the frozen and oscillatory behavior alternate. In our BBN analysis, we use an effective solution where φ 2 /Λ 2 is averaged over fast oscillations while slow oscillations and frozen regions are resolved without averaging. This is obtained by patching the numerical solution with the WKB amplitude for the rapidly oscillating regions. We also take into account the factor 1/2 which appears due to averaging (see appendix B for details). An illustrative example of the φ-evolution for m φ = 10 −20 eV and Λ = 10 17 GeV is shown in fig. 2 (right panel).

Analytic estimates
To understand the modification of BBN predictions due to ULDM, we first consider a simplified picture of BBN that will allow us to capture the main physical effects analytically. The full implementation of ULDM effects in the BBN code AlterBBN [32] will be discussed in the next subsection.

Standard picture
We begin with a recap of the key stages of the standard BBN dynamics following [36]. At early times, neutrons and protons are in thermal equilibrium through the reactions n + ν e ↔ p + e − and n + e + ↔ p +ν e which results in a neutron abundance (i.e. the ratio of the amount of neutrons to the total amount of baryons), where T is the temperature and m np 1.29 MeV is the difference between the neutron and proton masses. At T W 0.8 MeV the rates of the weak reactions become comparable to the Hubble expansion rate and further evolution of the neutron abundance is governed by the kinetic equation, Here λ n→p is the rate of neutron-to-proton conversion, which can be approximated as where G F is the Fermi constant, g A 1.26 is the nucleon axial charge, and the phase space integral has the form, 3 By solving the linear equation (4.2) we find the neutron abundance after the weak freezeout, where we have switched from time to the scale factor as the integration variable for later convenience. The equilibrium neutron abundance X eq n and the rate λ n→p in this formula are understood as functions of the scale factor. Evaluating the integral using the standard thermal history yields, After freeze-out, the number of neutrons continues to slowly decrease due to neutron decay, where the neutron lifetime is Γ −1 n 880 s. The onset of BBN is held back by the efficient dissociation of deuterium, which is often referred to as the deuterium bottleneck. This delays BBN until the universe cools to a temperature which allows deuterium to survive long enough to form heavier elements. This temperature can be estimated from the Saha equation for the equilibrium abundance of deuterium, MeV is the deuterium binding energy and η b = (6.104 ± 0.058) · 10 −10 is the baryon-to-photon ratio, as inferred from the Planck data [12]. The bottleneck opens when the rate of conversion of deuterium into heavier elements becomes comparable to the expansion of the universe. This occurs when its abundance reaches X D ∼ 2 × 10 −5 [36]. At this moment, the deuterium burning is still much slower than the reaction n + p ↔ D + γ maintaining deuterium in equilibrium with neutrons and protons, and the deuterium abundance continues to grow until it reaches its maximal value X max D ∼ 10 −2 . According to eq. (4.8), this happens at temperature T BBN 75 keV. Note that it is only logarithmically sensitive to the the precise value of X max D and corresponds to the time t BBN 230 s used in eq. (4.7). At this epoch the deuterium burning is very fast and its abundance drops, whereas the majority of available neutrons end up in 4 He due to its high binding energy, resulting in the primordial 4 He mass fraction of This estimate agrees well with the precise calculation using numerical codes [12], The standard theory also predicts the primordial abundances of other light elements, such as deuterium and lithium [12], 4 The experimentally measured values of these abundances inferred from astrophysical observations are [37], We see that for 4 He and D the theory and experiment are in very good agreement. In particular, for 4 He the relative difference is where the error is dominated by the experimental uncertainty. This relative difference constrains how much the addition of dark matter is allowed to change the SM prediction.
On the other hand, for Li the theoretical prediction is almost 3 times higher than the observed value -the mismatch is known as the "lithium problem". For constraints on ULDM, we will use only 4 He and D abundances.

Effects of dark matter on Helium-4 abundance
The abundance of 4 He is of particular interest as it depends only on a few factors. We now derive an analytic estimate of the ULDM impact on it using the formulae from the previous subsection. To this end, we consider the relative change where we have used eqs. (4.7), (4.9) and have written the factor describing neutron decay as an integral over the scale factor between the weak freeze-out and BBN. Note that we have allowed the neutron decay width to be time-dependent due to its modulation by ULDM. We first work in the Einstein frame and then present the Jordan frame description.
Einstein frame As discussed in sec. 2 (see also appendix A), the universal coupling of DM to SM affects the masses of particles and other dimensionful quantities, while leaving the dimensionless ratios intact. The Fermi constant has mass dimension −2, whereas the neutron-proton mass difference, the neutron width and the deuterium binding energy all have unit mass dimension. Thus, we have, (4.14) Next, it is shown in appendix C that the SM energy density and the Hubble rate in the Einstein frame remain the same functions of the scale factor, as in the standard cosmology, ρ SM (a) ρ The corrections to these expressions are doubly suppressed by the product of α and the small ratio Σ introduced in eq. (3.5). As the main effect of DM on BBN is of order α, these corrections can be neglected. We now perturb the formula (4.5) for the neutron abundance after freeze-out. The perturbation function α φ(a) enters through the Fermi constant and the neutron-proton mass difference. To simplify the result, we assume that during the relevant epoch temperature scales approximately as T ∝ 1/a. This is equivalent to neglecting the change in the effective number of relativistic degrees of freedom, which is a good approximation as long as T m e . We arrive at, As long asλ n→p > 1, i.e. the reactions are faster than the expansion of the universe, the integrand is strongly suppressed by the exponential factor. This is the case before the weak freeze-out. On the other hand, at low temperatures, the integral is cut off by the hyperbolic cosine in the denominator. Thus, the integral is saturated at temperatures around T ∼ T W . The two terms in the curly brackets can be traced back to the modification of the equilibrium neutron abundance and the reaction rate respectively. Note that they enter with opposite signs. Depending on the behavior of α(a), one or the other will win. This means that both over-and under-production of neutrons is possible in different regions of ULDM parameters. This is in stark contrast with the approximation that the weakinteractions freeze-out instantaneously and ∆X n,W ≥ 0 for all ULDM parameters. From the Saha equation (4.8) we infer that the BBN temperature scales as the deuterium binding energy. Hence, and This gives the perturbation of the upper integration limit in the integral describing the neutron decay in eq. (4.13). The change of the lower limit is irrelevant. Indeed, the ratio of the neutron decay width to the Hubble rate at week freeze-out is (Γ n /H) W ∼ 3 · 10 −3 . Thus, even an order-one change in a W would affect the helium abundance only at the level of a few per mille. Of course, in our case the change is further suppressed by the small coupling α. Thus, we can set a W in eq. (4.13) to be equal to its standard value a (0) W . The variation of the decay width Γ n takes place in the intervening times between the weak freeze-out and the deuterium burning. Thus, the shift ∆Γ n must be evaluated at the appropriate scale factor, ∆Γ n Γ n (a) = α φ(a) . (4.20) Combining everything together, we arrive at the following expression, where ∆X n,W is given by (4.16). In the next subsection we will use this formula to derive the constraints on the ULDM model.
Jordan frame It is instructive to repeat the above calculation in the Jordan frame to highlight some subtleties in the connection between the two frames. Here, the effect of ULDM enters only via the modified Hubble rate (see appendix C), whereas all masses and widths remain unchanged. The perturbation of the neutron abundance upon freeze-out, eq. (4.5), then reads, (4.23) At first sight, this expression differs from the Einstein frame expression (4.16). In particular, it involves the derivative of α with respect to the scale factor while the latter does not. However, it is straightforward to verify that removing the derivative by integration by parts and taking into account the dependence of the temperature on the scale factor T ∝ 1/a brings eq. (4.23) exactly to the form (4.16). While this equivalence is to be expected, it is worth mentioning that an oversimplification of the problem could break it down. In particular, this would be the case if, to estimate the neutron abundance, one used the approximation of an instantaneous freeze-out at the temperature when the weak reaction rate becomes equal to the Hubble rate. Then the Einstein frame description would be sensitive only to the value of α at that moment, whereas the Jordan frame one would also depend on its derivative. Thus, consideration of the full kinetic equation (4.2) is essential for the consistency of the analysis.
Similarly, the contribution accounting for neutron decay in eq. (4.13) in the Jordan frame is, where passing to the second line we have integrated by parts. As discussed above, the first term is small and can be safely neglected. For the derivative of the Hubble rate in the last term we use eq. (C.13) from the appendix C, which yields, Neglecting, as usual, Σ-suppressed contributions, we again reproduce the Einstein-frame result. This serves as a cross-check of our calculation.

Analytic constraints from Helium-4
By combining the limits on the 4 He abundance (4.12) with the impact on BBN derived in the previous section (4.21), we arrive at the constraints for ultralight scalar DM. We first consider the case of an ULDM field with positive coupling α = φ 2 /Λ 2 . This yields constraints on the DM parameter space shown in the left panel of fig. 3. Also shown are the constraints from a model which neglects the effects due to the DM induced mass and assumes an instantaneous weak freeze-out [30]. Including the DM induced mass dramatically modifies the constraints. They are strengthened at low masses m φ 10 −22 eV, weakened by up to two orders of magnitude at high masses m φ 10 −18 eV, and show non-trivial features at intermediate masses.
To understand this behavior, we first look at the two limiting cases at low and high mass where the ULDM field, and therefore α, has a straightforward time-dependence. For very light masses, the field is in the Hubble-friction (H) dominated phase during BBN 5 where H m φ,eff . Thus, the field is frozen at a constant value and its effect is analogous to changing M pl as shown in eq. 2.11. Later on the field begins oscillating with decreasing amplitude. The induced mass shifts the start of oscillations to an earlier time, which gives more room for the amplitude to decay. This increases the frozen value of the field at BBN and leads to stronger constraints.
On the other extreme, for high masses (m φ 10 −18 eV) the ULDM field is dominated by the induced mass (I-regime) throughout BBN. Here the field is rapidly oscillating and its amplitude is described by the WKB formula (eq. 3.12) both at BBN and at later times. The effective mass in the denominator strongly reduces the φ-amplitude at BBN, compared to the case without induced mass, weakening the constraints. Matching the WKB amplitude to the present-day DM density ρ DM,0 and using that today m φ,eff = m φ , whereas at BBN it is m φ,eff = √ 2Θ SM /Λ, we can find α as a function of the scale factor at BBN, This expression implies that BBN is sensitive to the combination m φ Λ, as indeed seen from fig. 3. Notice that in this regime the field amplitude features a rapid decrease with the scale factor. Thus most of the constraining power comes from the earliest stage of BBN, i.e. the weak freeze-out. The scaling in the intermediate mass range is more complicated due to the non-trivial time dependence of the induced mass. In this parameter region the field exhibits oscillations on the time scale comparable to the expansion rate of the universe. Depending on whether the oscillations happen to be on the peak (in the trough) during the relevant stages of BBN, the constraints are strengthened (weakened). This explains the oscillating features in the exclusion line clearly visible in fig. 3.
The reduction of BBN sensitivity due to the rapid time-dependence of the DM field allows the coupling α to approach order-one values during BBN. This is dangerous, since in our analysis we assume α to be small and expand at linear order in it. As a criterion for the validity of this approximation we impose the requirement that α should be less than 1 at the moment of the weak freeze-out marking the beginning of BBN. The corresponding parameter region is delimited by the black dotted line in fig. 3. We have verified that our constraints always lie in the region where this criterion is satisfied. Still, our tests show that inclusion of non-linear terms in α can shift the constraints on m φ or Λ by a factor of 2.
The dynamics of neutron freeze-out has an important qualitative impact on the predictions of the model. From eq. 4.16 we saw that, depending on the choice of ULDM parameters, neutrons, and therefore 4 He, can be either over-or under-produced. In the H-dominated regime, the modification of the reaction rates dominates and we increase the amount of 4 He. In the I-dominated regime, the modification of the equilibrium neutron abundance dominates and the amount of 4 He decreases. The transition between over-and under-production happens along the horizontal line Λ −1 ≈ 10 −18 GeV −1 in fig. 3. Note that the instantaneous weak freeze-out approximation leads to an over-production of 4 He across all of the parameter space, in stark contrast to the actual predictions.
For a negative coupling α = −φ 2 /Λ 2 DM becomes tachyonic in the induced mass dominated phase. This corresponds to the blue shaded region in the right panel of fig. 3. In this case the initial conditions for the field φ must be strongly fine-tuned to avoid its overproduction. We find this fine-tuning unappealing and do not explore the tachyonic regime further. The ULDM does not produce any observable effects on BBN outside of the tachyonic region. Interestingly, the allowed parameter space contains a band where the DM stays non-tachyonic during cosmological evolution since the start of BBN, but still has a coupling strong enough to develop instability inside extremely compact objects such as neutron stars. This may lead to a phenomenon known as scalarization (see e.g. [38][39][40][41]) during which the compact object spontaneously acquires a scalar charge, and can have observable signatures in e.g. binary pulsar systems. The criterion for acquiring a scalar charge is that 2Θ NS /Λ 2 max{m 2 φ , R −2 }, where Θ NS is the trace of the neutron star EMT and R ∼ 10 km is the neutron star size. Neglecting for an estimate the pressure contribution and taking the neutron star density ρ NS 3 × 10 14 g/cm 3 [42], we obtain that scalarization is possible in the region above the green line in fig. 3. A more detailed investigation of scalarization is justified, but is beyond the scope of this work.

Numerical constraints from Helium-4 and Deuterium
The calculation of the abundances of the light elements beyond 4 He requires the study of a complex network of nuclear interactions. Traditionally, this problem is handled with numerical codes such as AlterBBN [32], PRIMAT [11] and PArthENoPE [43]. Implementing our model into such a code in the Einstein frame involves continuously changing fermion masses and reaction rates, which would take a substantial effort. Instead of undertaking such an effort, we make use of the Jordan frame where the SM is left unchanged but instead evolves in a universe with modified expansion rate (eq. 4.22). In this numerical study we restrict to the case of positive DM coupling, α > 0.
We use the numerical code AlterBBN [32], which parameterizes changes to the evolution of the scale factor (and Hubble parameter) by introducing a dark density component in the Friedman equation. The dark density required to reproduce the Hubble rate (eq. 4. 22) has the form, Note that, since α can be rapidly decreasing or oscillating, the derivative term in brackets can be negative, pushing ρ dark < 0. Using the standard thermal history of the universe, we convert ρ dark into a function of temperature. Then we import it into AlterBBN v2.2 through the option Init dark density table, such that the evolution is modified to correspond to that implied by our model. In this way one can make numerical estimates of not only 4 He but also D, 3 He, 6 Li, 7 Li, and 8 Be. Unfortunately, experimental data suitable for comparison with these predictions only exist for the elements 4 He, D, and the combined lithium abundance [37]. One needs to proceed with caution in defining ρ dark . We find that in some regions of parameter space the DM amplitude may become large (φ 2 /Λ 2 ∼ O(1)) in the time range referenced by AlterBBN. Here ρ dark , calculated according to (eq. 4.27), can exceed the standard model density. This is inconsistent because it signals the breakdown of the linear expansion in α used in our analysis. To avoid this problem, we implement a cut which limits the magnitude of dark density to half of the SM density, i.e. we require that the dark density is bounded by |ρ dark | < 0.5ρ SM . This cut, which is only necessary in the region near the breakdown of the model, can be understood as a conservative measure, which discards constraints arising from α ∼ O(1). We only report constraints where α < 1 at the time of the weak freeze-out. However, AlterBBN begins integrating at earlier times, where α may be larger than at weak freeze-out. Thus, the cut is used for regions where α is small at freeze-out but grows to O(1) towards the earliest times probed by AlterBBN.
The numerical 4 He constraint agrees well with the analytical counterpart in the vast majority of parameter space, as seen in fig. 4. The minor discrepancies in the oscillating region can be traced back to the cuts, which are required in the Jordan frame analysis as described above, but are absent in the analytical treatment. The cuts have a particularly large effect in the oscillating region, because this region features a pronounced degree of cancellation, such that larger values of α are probed. This also implies that the oscillating region is more sensitive to non-linearities in α and model-specific higher-order terms in the Taylor expansion of the function α(φ). These are not captured in our analysis, and so the constraints in this region should be interpreted with care. There is no analytical counterpart for the numerical D exclusion, which is shown in the same figure. Note that the 4 He constraint is dominant for most of parameter space.
With regards to lithium, it is well known that predictions of standard BBN theory do not match the observed abundances. This is also present in our model for the regions of parameter space which are compatible with the observed 4 He and D abundances. Still, it is worth noting that the DM modification tends to slightly reduce lithium production, but this change is vastly insufficient to resolve the lithium problem.  Figure 4: The numerical constraints on 4 He (black, solid) and D (black, dashed) derived using AlterBBN as compared to the constraints from the analytical approximation of the preceding section (red and orange regions). Again, red corresponds to underproduction of 4 He and orange corresponds to its overproduction. Deuterium tends to be overproduced. The black dotted line delimits the region of parameters where the coupling becomes large at the scale factor a W = 10 −9.6 corresponding to the weak freeze-out, α(a W ) > 1. To the left of this line higher-order terms in the Taylor expansion of the function α(φ) become important. The right panel zooms in to the region where the constraints exhibit oscillations. Here, the constraints should be taken with caution as they are sensitive to non-linear effects in α.

Additional constraints
The quadratic coupling of a new scalar particle to the SM is subject to many other astrophysical constraints, ranging from the anomalous cooling of supernovae, modifications to the timing of binary pulsar systems, and black hole superradiance, as well as laboratory constraints coming from atomic clocks and torsion pendulums. The constraints from torsion balance experiments and atomic clocks probe WEP violation, and thus do not apply to our model. In what follows, we briefly summarize the strongest constraints on the quadratic coupling of a new scalar particle.
Supernova cooling: New light particles with masses below the average temperature of the core of a supernova, T SN core 30 MeV, can be produced within the supernova core and free-stream out, thus removing energy from the supernova. This new channel for energy loss can alter the supernova neutrino luminosity. Although a proper calculation of the free-streaming process is involved and requires complicated simulations, one can approximate an upper limit on the instantaneous luminosity of new particles produced in the supernova based on the observation of SN1987a [44]: L new ≥ L ν 3 × 10 52 erg/s. If the instantaneous luminosity of new particles exceeds this value when the core reaches peak density, ρ c 3×10 14 g/cm 3 and temperature T c 30 MeV, then the energy spectrum of the neutrino burst from SN1987a becomes inconsistent with observations. This is known as the "Raffelt criterion". 6 The dominant production process for the scalar is the bremsstrahlung off neutrons nn → nnφφ due to the high density of neutrons in the proto-neutron star. The rate of φ production can be estimated as [46] Γ nn→nnφφ σ nn × where m n is the neutron mass, and where σ nn 25 mb [47] is the elastic nucleon crosssection. The Raffelt criterion can be rewritten as Γ nn→nnφφ 10 −14 MeV 5 , which corresponds to Λ 14.5 TeV.
Fifth-force experiments: Fifth-force searches traditionally set some of the strongest constraints on ultra-light scalar particles [48,49]. The quadratic coupling of the scalar produces a fifth-force at leading order through the exchange of a pair of φ between two fermions. This generates a potential of the form in contrast to the usual 1/r Yukawa potential from linear-couplings. As a consequence, the constraints from fifth-force experiments on the quadratic coupling are much weaker than those on the linear coupling [48,49]. Limits on the deviation from the typical 1/r behavior of the gravitational potential require Λ 2 TeV for m φ 10 −4 eV [46]. Note, however, that the presence of a non-zero φ background induces an effective linear coupling between φ and matter. It would be interesting to reinterpret existing constraints in the context of this effect. 7 Galaxy formation and Ly-α: A lower limit on the scalar DM mass comes from the cosmological structure formation. The latter becomes affected if the de Broglie wavelength of the scalar is larger of comparable to the relevant astrophysical scales, such as the size of dwarf galaxies R ∼ 1 kpc. Requiring that modifications of the structure formation are compatible with observations leads to the bound m φ 10 −22 eV [50][51][52], assuming that the scalar composes the entirety of DM density. A more stringent lower bound comes from the measurement of the Ly-α forest, which requires m φ 10 −21 eV [53,54]. This is supported by complementary analyses of the galactic rotation curves [55] and the subhalo mass function [56]. If the detection of a global 21-cm absorption signal [57] is confirmed, it will set a similar bounds [58][59][60]. In addition, the fluctuations from light scalar DM can heat up the cores of galaxies and modify their dynamics. Based on the properties of the star cluster in Eridanus II, ref. [61] has inferred the bound m φ 10 −19 eV. All these limits can be relaxed if the scalar is only a sub-component of the dark matter.
Pulsar Timing and Stochastic Gravitational Waves: Binary pulsars are systems whose dynamics are measured with exquisite precision. As a result, they are highly sensitive to any new physics that changes their dynamics [62,63]. As discussed, the addition of an ultralight scalar field with universal couplings will perturb the masses of the SM particles, and thus change the masses of the stars in the binary, resulting in a change in the orbital period. The effect is especially pronounced if the frequency of the perturbation matches the harmonics of the binary orbital frequency. In this situation, the changes in orbital parameters are resonantly amplified and can lead to measurable effects [28,29]. Current data give strong bounds on the coupling in the resonant bands of masses; future observations with the Square Kilometer Array telescope will significantly strengthen the constraints and extends the coverage. Complementary are provided by Pulsar Timing Arrays (PTA) [64] as well as searches for stochastic gravitational waves by the Cassini (CAS) space mission [65]. These constraints are weaker than the BBN constraints derived in this work.
Black Hole Superradiance: Ultralight scalar fields can form gravitationally bound states with black holes if their Compton wavelength is of comparable size to the Schwarzschild radius of the black hole. The scalar fields extract angular momentum from the black hole through a process known as superradiance, which cause rapidly spinning black holes to spin down [6,66,67]. 8 Observations of old, near-extremal black holes can therefore exclude the existence of weakly-interacting scalars in the mass ranges of m φ ∈ [10 −18.2 , 10 −17.6 ], [10 −16.7 , 10 −16.1 ], and [10 −13 , 10 −10.8 ] eV [69]. Measurements of M87 * by the Event Horizon Telescope [70] further constrain scalar masses m φ ∈ [2.9×10 −21 , 4.6×10 −21 ] eV [71]. These bounds assume that the self-interactions of the scalar are negligible. Large self-interactions and other non-linear effects can render superradiance ineffective, thus allowing for the scalar to evade the aforementioned bounds [66,72,73].
A summary of all the constraints described in this section, along with the BBN constraints from this work, are collected in fig. 5.

Conclusions
In this work, we have investigated the effect of ultralight scalar DM with universal quadratic couplings to SM fields on predictions of BBN. This type of coupling preserves the weak equivalence principle and does not give rise to a long-range force, thereby evading numerous laboratory tests. We took into account the full dynamics of BBN, as well as DM cosmological evolution. We have found that the precision measurements of the primordial 4 He abundance constrain a large portion of ULDM parameter space as shown in fig. 5.
The quadratic coupling has an important effect of the evolution of DM endowing it with an effective mass proportional to the energy density of the SM environment. This leads to substantial corrections to previous calculations [30], where this effect was neglected. Specifically, for positive coupling the constraints from the full ULDM evolution are about two orders of magnitude weaker for m φ 10 −18 eV and have non-trivial features at m φ 10 −18 eV where the ULDM exhibits oscillations on a time scale comparable to the expansion rate of the universe. In contrast, at very low masses our constraints are stronger. For negative coupling we find the DM has a tachyonic instability in a large portion of the parameter space, requiring an extreme fine-tuning of initial conditions to avoid its overproduction. In addition, we moved beyond the instantaneous approximation for neutron freeze-out and used the full kinetic description. We have found that it qualitatively changes the result for the 4 He abundance, which is predicted to decrease in much of the parameter space; in contrast, the instantaneous approximation erroneously predicts an over-production of 4 He across all of parameter space.
The universal coupling studied in this work allows for treatment in an alternative frame, the Jordan frame, where the ULDM modifies the metric, but does not couple directly to the matter fields. The modified metric manifests itself as a modification to the Hubble parameter, which we implemented into a numerical code AlterBBN v2.2. The numerical analysis allows us to determine not only 4 He abundance, but also that of other light elements, such as deuterium and lithium. We found the limits based on deuterium are subdominant to 4 He throughout the parameter space. Curiously, the ULDM leads to a slight reduction of the Li abundance, although this is vastly insufficient to resolve the "lithium problem".
Our results show that BBN sets the strongest constraint on the ULDM coupling for the majority of the parameter space where 10 −19 m φ 10 −6 eV. However, we note a few caveats. Our constraints are subject to the assumption that SM and ULDM are described by a simple Lagrangian (eq. 2.2) all the way back to BBN. We further assumed that we can Taylor expand the coupling α(φ) keeping only the leading-order term. Violation of these assumptions may modify the constraints. Inclusion of a self-interaction of the DM field or its more complicated coupling to SM can significantly alter the dynamics of the ULDM and lead to rich phenomenology (see e.g. [31,[74][75][76][77][78]). The specific form of these non-linear terms could be determined within a more complete setting describing the UV origin of the universal coupling. Further study along these lines is warranted and we leave this for future work.
assuming α(φ) is small. The neutron-proton mass difference m np receives contributions from the isospin symmetry breaking by the quark masses and from the electromagnetic interaction, where δ iso = 2.05 ± 0.30 MeV and δ em = −0.76 ± 0.30 MeV [79]. Here α em is the fine structure constant. A more recent ab initio lattice calculation finds 10 δ iso = 2.52±0.29 MeV and δ em = −1.00 ± 0.17 MeV [80]. Independently of the precise values, we have where we have used eqs. (2.6), (2.9). The neutron decay rate can be approximated by [17] where g A 1.26 is the nucleon axial charge and Taking the variation of all masses as before we arrive at the fractional change in decay rate The sensitivity of the deuterium binding energy B D can be estimated in a variety of ways [24]. Our results are not particularly sensitive to the exact scheme. Following [17], we take where ∆m q /m q is the variation of the quark masses.

Appendix B Dark matter evolution
In this appendix we describe our patching procedure for the cosmological evolution of the DM. The patching procedure is required because of rapid oscillations, which are present both in the bare mass regime (B) and in the induced mass regime (I), when positive couplings are considered. These oscillations can make calculating full numerical solutions computationally expensive. Alternatively, the WKB-approximation allows us to accurately reproduce the evolution in the regime where the field is rapidly oscillating. Since the WKBapproximation is only valid in the rapidly oscillating regime, we must patch between the full numerical solution and the WKB-approximation. We implement the patching as follows: 1. Compare the total effective mass m φ,eff to the Hubble parameter to determine if and when the Hubble friction becomes relevant. Specifically, we require that the Hubble parameter H is bigger than 0.2 m φ,eff before the numerical solution is applied. This determines when to transition between the full numerical solution and the WKBapproximation. Note that in the case of intermediate fast regimes, triggered by the induced mass, there will be several transitions.
2. Numerically solve the equations of motion from some initial time (log(a) = −12) and until the latest transition time (i.e. the highest value of a) found above.
3. In order to ensure the correct normalization, we need to match the solutions at the peaks of the oscillations. We therefore numerically search for peaks in the oscillations near the desired transition times. We perform the searches away from the slowly varying regions to ensure that we always find peaks in the fast regime.
4. The amplitude of the WKB solution is fixed to the present day DM density. We can therefore fix the amplitude of the numerical solution by matching the amplitude of the numerical solution to the WKB solution at the peaks found above. If an intermediate fast regime is present, then the process is reversed and the amplitude of the intermediate WKB solution is matched to a peak of the numerics.

5.
For the BBN analysis we need the value of φ 2 /Λ 2 averaged over fast oscillations. This is smaller by a factor 2 than the actual envelope of the field. As we normalize our WKB amplitude to track the average of the numerical solution, we have to multiply the numerical solution by a factor 2 when transitioning from the WKB to the Hubble friction dominated regime. To prevent the evolution from being discontinuous we gradually impose the factor of 2 by turning it on or off over a period of a Hubble time.
We illustrate the matching procedure in fig. 6, indicating the various types of solutions: full numeric, WKB-amplitude, and the effective solution. Here, we show the different transition times, indicated by thin black vertical lines, which are found by comparing the total mass and the Hubble parameter. The solutions are matched at the peaks, marked with thick, black vertical lines. The amplification function, which takes the factor of 2 from averaging into account, is also visualized with a gray dotted line, although the function is shown with an artificial amplitude for visualization purposes. The effective solution used for the BBN analysis is shown as a red, solid curve. It results from the patching of the full numeric solution and the WKB-amplitude for the rapidly-oscillating regions. For negative coupling the patching procedure is similar, except that in the induced mass regime (I) the solution is exponentially growing, instead of oscillating. In this case we do not use the WKB approximation in the I regime, but solve the equation exactly.  ULDM perturbs the expansion of the universe during BBN. We keep only the leading order corrections in α(φ). The general expressions can be found in Ref. [21]. This system must be complemented by the equation of state relating p SM to ρ SM and φ, where the dependence on φ appears due to the φ-modulation of the SM parameters (see below).
The last expressions imply the Einstein frame equation of state, Notice that this expression is in general φ-dependent. Only for a linear equation of state p (0) (ρ SM ) = wρ SM , with constant w, the φ-dependence drops out and one recovers the standard relation p SM = wρ SM . Let us now discuss how the ULDM perturbs the expansion of the universe in the two frames during BBN. In the Einstein frame the Hubble rate obeys the Fiedmann eq. (C.1a) where the last two terms represent the DM energy density. We argue that the latter is negligible. The DM density is comparable to that of radiation at the scale factor a eq ∼ 10 −4 . If the DM dynamics is dominated by the bare mass, the ratio ρ φ /ρ SM scales like a/a eq and is less than 10 −4 at the BBN epoch, which is too small to affect BBN in an observable way. In the Hubble friction regime, when the value of φ is frozen, the ratio ρ φ /ρ SM is suppressed even further. It remains to check what happens when the induced mass dominates. In this case we can use the WKB solutions (3.10), (3.11) to obtain This expression is doubly suppressed by the small quantities |α(φ)| and the small ratio Σ introduced in eq. (3.5). We saw in the main text that the leading effect on BBN is of order O α(φ) , so the contributions of the form (C.11) can be neglected. Now we need to determine the effect of the ULDM coupling on the SM energy density ρ SM . This is non-trivial, since the SM energy is not conserved due to the direct DM coupling, see eq. (C.1c). Moreover, the equation of state relating the SM pressure to the energy density is modified, eq. (C.10). To overcome these complications, we use the map to the Jordan frame, where the energy density has the standard dependence on the scale factor (C.3). Using the first of eqs. (C.9), we obtain, ρ SM (a) = (1 + 4α)ρ We see that, though ρ SM changes with respect to the standard cosmology, the deviation is doubly suppressed by α and Σ. Thus, we can omit it within our approximation. Substituting into the Fiedmann equation, we conclude that, up to terms of order O(αΣ), the Hubble rate in the Einstein frame is described by the same function of the scale factor, as in the standard cosmology. In this way we arrive to eqs. (4.15) from the main text. The Hubble rate in the Jordan frame is determined by eq. (C.2a). Using the same arguments as for the φ-density in the Einstein frame, one can show that the two terms in the second line of this equation are negligible. Recalling also that the dependence of the SM density on the Jordan scale factor is standard, we immediately obtain, which is equivalent to eq.