Constraints on Ultralight Scalar Dark Matter with Quadratic Couplings

Ultralight dark matter is a compelling dark matter candidate. In this work, we examine the impact of quadratically-coupled ultralight dark matter on the predictions of Big Bang Nucleosynthesis. The presence of ultralight dark matter can modify the effective values of fundamental constants during Big Bang Nucleosynthesis, modifying the predicted abundances of the primordial elements such as Helium-4. We improve upon the existing literature in two ways: firstly, we take into account the thermal mass acquired by the ultralight dark matter due to its quadratic interactions with the Standard Model bath, which affects the cosmological evolution of the dark matter. Secondly, we treat the weak freeze-out using the full kinetic equations instead of using an instantaneous approximation. Both improvements were shown to impact the Helium-4 prediction in the context of universally-coupled dark matter in previous work. We extend these lessons to more general couplings. We show that with these modifications, Big Bang Nucleosynthesis provides strong constraints of ultralight dark matter with quadratic couplings to the Standard Model for a large range of masses as compared to other probes of this model, such as equivalence principle tests, atomic and nuclear clocks, as well as astrophysical and other cosmological probes.


Introduction
Most of the matter in the Universe is not accounted for in the Standard Model (SM) of particle physics. This missing matter is referred to as dark matter (DM). Proposed candidates for DM range from many solar mass compact objects to ultralight DM with masses as low as about 10 −19 eV [1][2][3]. In the low mass range, the DM is characterized by high occupation numbers, causing DM to behave like a classical field. Familiar examples of such ultralight DM (ULDM) are axion-like-particles [4][5][6] and dilatons (see, e.g., [7]). Dilatons are scalar fields φ present in models with compactified extra dimensions and have leading order couplings that are of the form [8] where O SM are terms in the SM Lagrangian and Λ is the scale of new physics generating the coupling. Scalars with such couplings generate new long-range forces, and such theories are, therefore, heavily constrained by 5th force experiments [9][10][11][12]. In this work, we will consider instead scalar ULDM with quadratic couplings of the form φ 2 Λ 2 O SM , (1.2) which evades the most stringent constraints on dilatons and possesses novel phenomenology distinct from models with linear couplings. Such quadratic couplings can arise if the new physics assumed at the high-energy (UV) scale Λ is subject to symmetries that restrict the linear contributions. The lack of definitive observation of ULDM indicates that the ULDM must have extremely feeble couplings to ordinary matter. Beyond fifth-force searches, past attempts to search for ULDM have included searches for frequency variation of atomic clocks [13] and searches for perturbations in the arms of gravitational wave experiments [14]. Astrophysical probes of ULDM can be derived from Supernova SN1987A energy loss arguments [15,16], from the spin-down of black holes due to superradiance [17][18][19][20][21][22][23], and from pulsar timing arrays [24,25].
In this work, we investigate cosmological probes of ULDM. Cosmology provides an attractive probe of ULDM because the field values of ULDM decrease with the expansion of the Universe, such that the largest field values are realized at early times. The earliest current probe of physics in the early Universe comes from Big Bang Nucleosynthesis (BBN), which therefore provides a powerful constraint of ULDM. This constraint appears because the presence of a background DM field can shift the effective value of fundamental constants and therefore can modify BBN [26]. To adequately model the effect of ultralight DM on BBN, the backreaction from the SM bath on the ULDM must be accounted for. Furthermore, the BBN analysis must go beyond instantaneous approximations of the neutronproton, or weak, freeze-out. Previous work [27] has shown the relevance of these improvements in the context of universally-coupled ULDM. In this work, we investigate DM with more general couplings, correctly accounting for the evolution of the DM and the dynamics of weak freeze-out. We begin in Section 2 with a description of the quadratically-coupled ULDM model and discuss the DM coupling to low-energy degrees of freedom. We also discuss the evolution of ULDM, properly accounting for the backreaction from the SM bath. In Section 3, we provide an overview of BBN and show how ULDM affects BBN, specifically the predictions of the Helium-4 ( 4 He ) abundance. In Section 4 we will discuss the resulting constraints and compare the BBN constraints to other known constraints. We conclude in Section 5.

Ultralight Dark Matter With Quadratic Couplings
We begin by adding to the SM a real scalar field, where m φ is the scalar mass. Couplings, such as the ones considered below, will lead to corrections to the bare scalar mass. Generically, these are of the form with Λ UV the scale of the UV completion and d (2) i are dimensionless couplings defined below. Naïvely, m φ , smaller than this correction could be considered unnatural because this contribution would have to be tuned out. However, there exist models in which these contributions are suppressed, e.g., see [28,29]. Here, we seek to describe the phenomenology of such a field in a model-independent way, and we therefore leave any discussion of radiative corrections and protection from these to other work on model implementations.
We consider scenarios in which the scalar field is charged under a symmetry, e.g., Z 2 , such that the leading coupling is quadratic, i.e., we take the interactions of the scalar with the SM to be in the form φ 2 Λ 2 O SM , where O SM is a term in the SM Lagrangian and Λ is some high scale at which new physics appears. For ease of comparison with results in the literature, we will follow the conventions of [8,[30][31][32] and parameterize the interactions of the scalar particle with the SM with the following Lagrangian: Here, M pl = 1.22 × 10 19 GeV is the Planck mass, β 3 is the QCD beta function, and γ mi are the anomalous dimensions of the u and d quarks. The superscript (2) specifies that these are the quadratic (as opposed to linear) couplings of the scalar. An alternative convention for parameterizing these coupling found in [26,[33][34][35][36] has This convention implies that the scale of the new physics mediating these couplings to the SM is around Λ . The conversion from this convention to the one used in this work is given by d Note that the interactions defined in equation 2.4 will give rise to φ-dependency in the fundamental "constants" of the SM. In a background φ, the fundamental constants will shift according to M pl . Here, α 1/137 is the fine-structure constant, Λ QCD the QCD confinement scale, and m f are the fermion masses. The quark mass-couplings are more useful written in terms of the symmetric and antisymmetric combinations: as physical quantities come in these combinations. In particular, the quark mass-couplings enter BBN through the neutron-proton mass difference, which depends on the anti-symmetric combination, and through the neutron axial coupling, which depends on the symmetric combination. The effective solution is a splice of the numerical solution and the WKB approximation. The thermal mass used for the photon evolution was computed in the high-temperature limit as explained in section 2.1.2. The difference between the numerical solution and the effective solution in the non-oscillating regions is due to factors of 2 present in averaging. Time is indicated both in terms of temperature T and in terms of the scale factor a, the latter of which is normalized to a0 = 1 today.

Dark matter evolution
We are interested in the field values of the ULDM at BBN. The present-day field value, φ 0 , is fixed by the current abundance of DM: ρ DM,0 = 1 2 m 2 φ φ 2 0 . To determine the value of the field around BBN, we evolve the scalar field backward in time according to its equations of motion, using the present-day field value φ 0 as a boundary condition,φ where the˙denotes derivatives with respect to time. The effective mass is a combination of the bare mass and the thermal mass induced by the SM bath: m 2 eff = m 2 φ + m 2 ind . During BBN, the SM bath is primarily composed of photons and electrons. Thus, only interactions with photons (d (2) e ) and electrons (d (2) me ) give significant contributions to the thermal mass as other species are not sufficiently abundant. Therefore, the quark (d (2) δm and d (2) m ) and gluon (d (2) g ) couplings do not provide a significant contribution to the thermal mass and the backreaction from these two couplings can be ignored.
To solve equation 2.11, we can examine three distinct regimes of DM evolution: In the B and I regimes, φ is rapidly oscillating, and performing the full computation is numerically expensive. Instead, we model the interactions of the ULDM with matter by taking the average over an oscillation period. We use a WKB-type solution, as in [27], to find the amplitude of the oscillation. This WKB-type approximation gives, where φ amp is the amplitude of the rapidly oscillating φ. In the H regime, the field value is constant due to the Hubble friction term dominating the equation of motion. In the transitions between the regimes, the field is slowly oscillating and we can solve the field evolution numerically and resolve the behavior of φ without the need for averaging. Figure 1 shows an example evolution for m φ = 10 −20 eV and d (2) e = 2000 or d (2) me = 2000 comparing the evolution computed using our method to the pure WKB and numerical solution. In the top panel showing the electron coupling case the field transitions from being Hubble dominated at early times, to being induced mass dominated, and finally to being bare mass dominated. In the induced mass domination regime the field is slowly oscillating when the induced mass is comparable to Hubble. We see the same behavior in the bottom panel, which shows the photon coupling case. We see that at early times the induced mass is comparable to Hubble and the field is slowly oscillating, and at late times it transitions to bare mass domination. In both panels, we see that at late times the field is rapidly oscillating, which is needed to satisfy the requirement that the field is DM today. The behavior shown in figure 1 is not generic and the evolution depends on the mass and couplings. As we see the induced mass is crucial to understand the behavior of the scalar field. In what follows, we discuss the induced thermal mass for the electron and photon couplings.

Electron coupling
There are two ways to calculate the induced mass from the interactions with electrons. The first is to calculate the induced mass using thermal field theory. The leading diagram is (2.13) In the imaginary time formalism described in [37,38], this diagram gives (2.14) where the sum is over the n Matsubara frequencies which each fix the time component of p to p 0 = i(2n + 1)T . The temperature-dependent part of this diagram is given by This calculation can be cross-checked with a second way of calculating the effective mass, which is to note that when electrons are on-shell the interaction pl φ 2 m eψ ψ can be written as where Θ e is the trace of the electrons' contribution to the stress-energy tensor. The effective mass is then given by [27,39] Note that for T m e , equation 2.15 implies that m 2 eff ∝ T 2 . In this regime, m 2 eff ∝ a −2 whereas H 2 ∝ a −4 . Thus, at early times (high temperature), the DM evolution is dominated by Hubble friction and the field is in the H regime. Thus, the mass induced on φ from finite density effects can lead to highly non-trivial evolution of the scalar field. This can be seen in figure 2, where the induced mass m 2 ind ≡ m 2 eff − m 2 φ is compared to both H and the bare mass m φ . In the example given in figure  2, which corresponds to the behavior also seen in figure 1, the field is initially dominated by Hubble friction, then starts oscillating because of the induced mass and the returns to a brief Hubble friction dominated state before finally settling into an oscillating state driven by the late time mass. Figure 2 also shows the mass that can be induced by a photon coupling, which are discussed in the next section.

Photon coupling
For the photon coupling, we calculate the contribution to the induced mass using thermal field theory. gives where similarly to equation 2.14 the sum is over the n Matsubara frequencies which set p 0 = i2nπT . This gives a scaleless integral and so vanishes in dimensional regularization, therefore this diagram does not contribute to the thermal mass. Instead, the leading diagram appears at two loops φ φ e e . (2.21) This gives , (2.22) summing over the n and m Matsubara frequencies which set k 0 = i2nπT and p 0 = i(2m + 1)πT . Evaluating the above expression is substantially complicated due to overlapping UV divergences, but it is relatively simple to compute in the limit where T m e . In this high-temperature limit, the contribution to the thermal mass is given by This high-temperature limit is a good approximation for the thermal mass at temperatures much larger than the electron mass. For lower temperatures, the finite electron mass suppresses the contribution. Therefore, we expect the high-temperature result to overestimate the thermal mass.
In the regime where the field is always oscillating, we expect this thermal mass contribution to somewhat relax the constraints. Therefore, the high-temperature result can be taken as conservative. This can be seen from the WKB-approximation, by which the amplitude of the oscillations can be described as φ ∝ m −1/2 eff a −3/2 . The thermal mass contribution decays in time, which then tends to compensate for the growth in a and therefore slow down the redshift of φ. This reduces the hierarchy between φ BBN and φ 0 . Since φ 0 is fixed by the zero-temperature mass and the observed DM abundance today φ BBN will therefore be smaller in a scenario with thermal mass contribution, as long as WKB is always a good approximation of the evolution. This argument breaks down in the regime in which Hubble friction becomes significant and WKB approximation fails. In this regime, a thermal mass contribution can overcome Hubble friction which increases the strength of the constraint. At low m φ the impact on the constraint is therefore a competition between opposing effects. For our purposes, we will present our results in the two limits using the high-temperature effective mass given in equation 2.23 and in the low-temperature limit where m eff m φ , which provides an envelope within which the true bound will sit.

Ultralight Dark Matter and BBN
We are interested in understanding the effect of the ULDM couplings on the predictions of BBN. We, in section 3.1, begin with a brief review of the standard BBN analysis, following the discussion in [40]. From this, we can calculate an analytic estimate of the 4 He abundance. In section 3.2, we perturb the standard result to obtain constraints.

Standard BBN
The 4 He abundance is determined by the abundance of neutrons at the time that fusion to helium becomes efficient. This in turn is determined by the abundance of neutrons when neutrinos decouple and the time at which deuterium can form. Before weak freeze-out, the weak reactions n + e + ↔ p +ν e and n + ν e ↔ p + e − are efficient and neutrons and protons have equilibrium abundances, where X eq n and X eq p are the equilibrium fractions of baryons in neutrons and protons, respectively, and m np is the proton-neutron mass difference. As the Universe expands and cools, the reaction rate for these processes falls below the expansion rate of the Universe, the Hubble scale, and so the neutron abundance falls out of equilibrium. At this point, the evolution of the neutron abundance X n is governed by the kinetic equation, where a is the scale factor (normalized to a 0 = 1 today), H =ȧ/a is the Hubble parameter, T is the temperature, and λ n→p is the reaction rate of n + (e + or ν) → p + (ν or e − ): Here, G F is the Fermi constant, g An is the weak axial coupling, and J is a phase space factor, given by At this point, the neutron abundance continues to decrease due to neutron decay. The neutron abundance at these times follows with Γ n the neutron inverse lifetime given by where P is a phase space factor given by Integrating equation 3.7, we find the neutron abundance at BBN: After weak freeze-out, neutrons and protons fuse to form deuterium through Initially, the photon abundance is much larger than the baryon abundance, and given that the deuterium binding energy is not very large compared to the temperature, equilibrium strongly favors the left side of equation 3.11. The deuterium abundance is given by the Saha equation where X p,n is the abundance of protons and neutrons, respectively, B D is the deuterium binding energy, and ζ is the Riemann zeta function. Note that the 4 He binding energy is much larger than B D so in equilibrium 4 He would dominate. However, the reactions needed to form helium D + D T + p , He + p , depend on the formation of deuterium and are suppressed by the low abundance of deuterium. This suppression of the deuterium abundance due to the background photons is known as the "deuterium bottleneck". Eventually when T ≈ B D 30 , the deuterium abundance becomes large enough such that the reactions 3.13 become efficient. At this time most of the neutrons become bound into helium and the 4 He abundance is Y p ≈ 2X n,BBN . (3.14) Not long after this, the Universe becomes too cool and diffuse for these nuclear reactions to be efficient. This signals the end of BBN. After this, but before the onset of star formation, the tritium and 7 Be decay to 3 He and 7 Li respectively 1 . Equation 3.14 gives Y p ≈ 0.25 [40]. A detailed numerical calculation including the kinetics of the fusion reactions given in 3.13 gives Y th p = 0.24672 ± 0.00061 using the baryon to photon ratio measured by Planck [42]. This is in excellent agreement with the PDG recommended average of Y ex p = 0.245 ± 0.003 [1].

Effect of varying fundamental constants
We now consider how the presence of ULDM modifies the standard analysis above. We first consider how φ varies the fundamental constants and we then consider how this impacts the neutron fraction at freeze-out and at BBN. Finally, we derive the impact of φ on 4 He .

Variation in fundamental constants relevant for BBN
As we saw in the previous section, the 4 He abundance produced by BBN depends on the deuterium binding energy B D , the proton-neutron mass difference m np , and the neutron axial coupling g A . In order to determine how the presence of φ affects BNN, one must determine how φ shifts m np , B D , and g A . These parameters do not have analytic expressions but can only be calculated from first principles using lattice gauge theory. The shift in m np and B D can be estimated using a combination of lattice and analytic methods as discussed in [43]; we summarize the relevant parts below. The neutron-proton mass difference can be written as where b is some constant number determined as in [43], giving bαΛ QCD ≈ −0.76MeV. Using equations 2.6, 2.7, and 2.8 we get As already stated above, the dimensionless parameter ϕ is defined as ϕ = √ 4πφ M pl . The dependence of B D on the fundamental constants can be parameterized as [43] Determining the dependence of g A on fundamental constants is more involved and requires input from lattice results. Lattice-QCD based attempts to estimate g A from first principles have limited statistics for physical quark masses and rely on calculations at heavier quark masses to fit formulae taken from chiral perturbation theory (ChPT) to extrapolate to physical quark masses. To estimate the dependence of g A on the quark masses we use these fits. Specifically, we use the expression with π = mπ 4πfπ . This expression is from NNLO heavy baryon ChPT. It neglects effects from ∆ resonances [44,45] and is used by lattice groups for their extrapolation to physical quark masses [46,47]. The coefficients g 0 , c 2 , and c 3 are determined from fits to lattice calculations. The dependence of g A on π can be used to infer the dependence on the quark masses using m π ∝ √m and neglecting the small dependence of f π on the quark masses. Putting this all together, we have This depends on the specific values for the coefficients g 0 , c 2 , and c 3 that come from lattice fits. The values obtained from [48] and its supplementary materials are . (3.23) These give ∂ln g A ∂ln π = −0.008 ± 0.023 .

(3.24)
Since the dependence of the 4 He abundance is proportional to d (2) m ∂ln g A ∂ln π , and the interval above contains zero, we are not able to give meaningful constraints on d (2) m from BBN at this time. Improved lattice results may constrain ∂ln g A ∂ln π away from 0, which in turn could give constraints on d m from BBN.
Armed with these expressions that encode the effects of the ULDM on the parameters that enter the calculation for the 4 He abundance, we can now determine the effect of ULDM on the 4 He abundance.

Variation in the neutron abundance after weak freeze-out
We begin by varying 3.6 to find the change in X n,W due to the presence of φ (or equivalently to ϕ). To first order in κ = d (2) i ϕ 2 2 , where i = e, m e , g, or δm, we find The change in ∆X n,W depends on ∆X eq n , dX eq n da , and ∆λ n→p . Let's first examine the variation in X eq n ∆X eq n = − m np 2T (1 + cosh m np /T ) ∆m np m np , (3.28) and compare this to . (3.29) The variation inλ n→p is given by with the variation in λ n→p given by The last two terms in equation 3.32 are proportional to m 2 e T 2 , which is small at relevant times and in principle can be neglected. However, we will include them going forward. Note that ∆G F G F = 0 for the models we are considering. Putting everything together we find that the change in the neutron abundance at weak freeze-out is given by (3.33)

Variation in the neutron abundance at BBN
After weak freeze-out, the neutron abundance is modified by neutron decay.

Variation in the helium abundance
The variation of equation 3.14 gives The remaining piece of our calculation is to determine ∆T BBN T BBN . Recall that the time of BBN is determined by the deuterium abundance reaching a critical value allowing the efficient fusing of deuterium into 4 He . Using equation 3.12, we can find T BBN by examining Note that the left-hand side is exponentially dependent on B D /T BBN . Therefore, any shift in B D must be compensated by a shift in T BBN , indicating that T BBN depends mostly on B D , i.e.
Combining everything, we get (3.40) Using equations 2.8, 3.16, and 3.17 for the variations in the fundamental constants we can now find the change in the 4 He abundance in terms of the field value of φ as a function of redshift. φ as a function of redshift was found using the evolution as described in section 2.1. We place constraints on the DM couplings by requiring that the resulting shift in the 4 He abundance keeps the 4 He abundance in the 95% confidence interval of the observed abundance, i.e. we require that Y ex Yp with σ Yp being the uncertainty of the theory prediction and the experimental uncertainty added in quadrature.

Results
The constraints on quadratically-coupled ULDM obtained from BBN in this work are shown as solid red lines in figure 3. For masses larger than about 10 −14 eV, the evolution is in the bare mass dominated regime both during and after weak freeze-out. For these masses, the constraints scale as d (2) i ∝ m 2 φ regardless of coupling. For the quark and gluon couplings, the field is Hubble frozen before BBN for masses lower than 10 −19 eV. This leads the constraint to scale as d (2) i ∝ m 1/2 φ in this low-mass regime. At intermediate masses, the field is Hubble frozen-in during or between weak freeze-out and BBN and the behavior is more complicated. For the electron coupling, the thermal mass is relevant for the evolution of the dark matter for masses less than 10 −18 eV. This leads to slow oscillations on the timescale of relevant BBN dynamics and the complex behavior shown in figure 3. For the photon coupling, we show two different constraints for two different estimates of the effective mass. The dark red line neglects the thermal mass of φ and the red line uses the thermal mass calculated in the high-temperature limit. We expect that the true constraint, calculated using the exact thermal mass, would be in between these two bounds as explained in section 2.1.2.
The BBN constraints derived in this work improve on the constraints derived in [26] by including the effects of the thermal mass on the evolution of the dark matter and by treating weak freeze-out using the full kinetic description. For the electron coupling, at masses of at least m φ 10 −14 eV the two effects lead to a constraint about two orders of magnitude weaker than the constraint given in [26]. This reduced constraint is explained by the thermal mass which reduces the field value at the BBN era. For the quark coupling, the thermal mass is insignificant and the constraint on this coupling is instead enhanced by about a factor 2 relative to [26] in the m φ 10 −14 eV regime. In all cases, for masses below 10 −14 eV we see nontrivial behavior not seen in [26]. For the photon coupling, our results indicate only the region of parameter space in which a full calculation of the thermal effects is expected to place the constraint. The constraint on the photon coupling given in [26] lies within the region where we expect the constraint to lie. Therefore, we cannot show a change in the photon coupling constraint without a more detailed evaluation of thermal effects, which we leave to future work. The gluon coupling was not treated in [26] and thus the constraints on the gluon coupling derived in this work are completely novel.
Our analysis is valid in the regime where d (2) i ϕ 2 1. In this regime, we can neglect the higherorder interactions of φ with the SM that would be present in a UV completion and treat φ as a small perturbation to the SM. Importantly, we can neglect the effect of the thermal quartic and other higher-order terms on the evolution of the ULDM. Because the amplitude of φ decreases with redshift, the threshold d (2) i ϕ 2 ∼ 1 could have been crossed sometime in the early Universe; if this happens after weak freeze-out then details of the UV completion are needed to determine the impact of φ on BBN. The parameter space in which the condition that d (2) i ϕ 2 1 is violated at any time after weak freeze-out corresponds to the lightly shaded regions above the dashed red lines in figure 3. In these regions, a model-dependent treatment is required for any given UV completion. Figure 3. Constraints on quadratic couplings of ULDM to electrons, quarks, gluons, and photons. Red regions indicate the BBN constraints derived in this paper. Other constraints are given in shades of gray. On the photon coupling plot, the conservative high-temperature approximation of the thermal mass yields the constraint in lighter red while the optimistic zero-temperature (i.e., without backreaction) result yields the constraint in darker red, see section 2.1.2 for details. The dashed red lines indicate where the approximation used in this paper breaks down and a more detailed analysis is needed. The grey regions show parameter space ruled out by prior work from the Eöt-Wash [9] and MICROSCOPE [49] experiments, atomic clock experiments [31,50] as well as constraints from Lyman-α [2], UFDs [3] and superradiance [23,51]. Also shown, with dashed gray lines, are projected constraints from the AION and AEDGE experiments [52,53] as well as projections for nuclear clocks [54]. Following the discussion in section 2, we show the line above which the low values of m φ may be rendered unnatural by large loop corrections, assuming a cutoff of Λ = 10 TeV and Λ = 2 GeV. For the gluon and photon plots, the line for a Λ = 10 TeV cutoff lies outside the range shown.

Other Constraints
Scalar ULDM with the couplings discussed in this paper is subject to several constraints in addition to the BBN bound derived here. These include constraints from searches for weak equivalence principle (WEP) violation, atomic clock experiments directly looking for varying fundamental constants, probes of cosmological structure formation, and astrophysical constraints.

Weak Equivalence Principle Violation Searches
Searches for violation of the WEP look for accelerations on test masses in addition to those expected from gravity. Such accelerations can violate the universality of free fall (UFF) and can be generated by long-range Yukawa forces mediated by φ. Classically, this effect appears because the φ-SM couplings perturb fundamental constants which leads to perturbations of the binding energy and thus the mass of SM matter. In general, the φ-dependence in SM matter will depend on the material composition, so that the resultant force will not respect the WEP. In the presence of background φ-field gradients or oscillations, this φ-dependence in the mass of SM matter lead to apparent WEP violations that are highly constrained by experiment.
Such constraints on apparent WEP violation depend sensitively on the background field configuration, which was studied in [31]. For linear interactions, the φ-SM interactions generate a static field configuration around the Earth itself, the gradient of which leads to WEP violation constraints. For quadratic interactions, no static solution beyond the trivial one exists. Instead, φ gets an effective mass inside the Earth and other massive bodies (see e.g. [55]). This effect screens the φ field 2 , creating gradients in the φ field in the vicinity of the earth.
The strongest constraints from WEP violation searches are set by the Eöt-Wash experiment [9] and the MICROSCOPE experiment [49]. The Eöt-Wash experiment searches for equivalence principle violation by monitoring a torsion balance with masses made of dissimilar materials (beryllium and titanium) and looking for torques generated by differing force per mass on the test bodies. For quadratic interactions, φ acquires an induced mass from the interaction with the Earth, which screens φ. This effect limits the sensitivity of Earth-based experiments. The MICROSCOPE experiment is a space-based experiment that looks for differing gravitational accelerations on test masses in orbit. Because of the reduced screening, the MICROSCOPE experiment provides stronger constraints than the Eöt-Wash experiment across all masses for quadratic interactions. These constraints are shown in figure 3.
This discussion of WEP violation constraints follows the derivation in [31]. In this work, it is assumed that far away from the Earth ϕ approaches with ρ DM the local DM density. However, this assumption neglects the virial velocity of the DM; the virial velocity can be neglected when the wavelength of the DM is larger than the radius of the Earth. However, for m φ 10 −11 the wavelength is comparable to or smaller than the size of the Earth and the WEP constraints derived in [31] and displayed in Fig. 3 may need to be modified. Note that this is the first time that the quark and gluon constraints from WEP searches have been presented for quadratic couplings.

Experiments looking for Varying Fundamental Constants
Precision measurements can search for temporal variation in fundamental constants in labs on Earth and in space. By examining the stability of different atomic frequency standards relative to each other, tight constraints can be put on the temporal variation in the ratios of different atomic transitions (not necessarily of the same species). This puts strong constraints on the coupling of ULDM to photons, quarks and gluons [13,31]. The most sensitive clock experiments are not sensitive to the electron coupling due to approximate cancellations for the frequency comparisons they make. Future improvements in this field may come from the development of nuclear clocks based on the nuclear transition between the ground state and the first excited state of the 229 Th isotope. If such a clock is realized, the current constraints could be improved by many orders of magnitude. For a review of recent developments in the field, e.g., see [54].
Matter-wave interferometry experiments are very sensitive to potential time variation in atomic transition energy [56]. A large number of upcoming experiments will probe interesting parameter space. A recent summary of progress in this field can be found in [54]. Among the most sensitive of such atom interferometers are the Earth-based AION [52] and the space-based AEDGE [53]. The Earth-based experiment is hampered by the screening effects discussed in the previous section. These experiments are also sensitive to the electron and photon couplings.

Structure Formation
Cosmological structures cannot form on length scales smaller than the de Broglie wavelength of the DM. To correctly reproduce structure formation, the DM must therefore be sufficiently heavy. Observations of dwarf Milky Way satellites require m φ 3 × 10 −21 eV [57][58][59]. Similar constraints come from measurements of the subhalo mass function and observations of stellar streams [60,61]. The large de Broglie wavelength also delays structure formation compared to standard cosmology, so complementary constraints arise from observations of small-scale structure at high redshift. For example, the Lyman-α forest flux power spectrum requires m φ 2 × 10 −20 eV [2].
The strongest lower bound on the mass of ULDM comes from observations of ultra-faint dwarf galaxies (UFDs). The wave-like properties of ULDM cause DM density fluctuations that transfer energy to stars through gravitational interactions leading to dynamical heating of dwarf galaxies. Measurements of the velocity dispersion of the UFDs Segue 1 and Segue 2 require m φ 3×10 −19 eV [3]. Similar constraints were claimed from observation of the dwarf galaxy Eridanus II [62]. However, that analysis made strong assumptions about the solitonic core of the DM halo of Eridanus II, which were relaxed in a more detailed analysis [63].

Astrophysics
Astrophysical probes are also sensitive to the presence of interacting ULDM. In particular, pulsar timing arrays can be used to constrain very light ULDM in a manner complementary to the structure formation arguments discussed above [24,25]. These constraints are weaker than the constraints from structure formation and are thus not shown here. At higher masses, constraints arise from superradiance. Bosonic radiation incident on a rotating black hole can be amplified through a process called black hole superradiance in which energy and angular momentum are extracted from the black hole in a manner similar to the Penrose process [64] and results in the spin-down of black holes [65]. The observation of old, rapidly rotating black holes can therefore be used to rule out the existence of light bosons [17-22, 66, 67]. Observations of solar mass black holes can be used to constrain light scalars with masses in the interval [2.7 × 10 −13 , 6.1 × 10 −12 ] eV [23]. Supermassive black holes can also be used to constrain DM mass in the interval [2.9 × 10 −21 , 1.6 × 10 −17 ] eV [20,51,67,68], although the spin measurements are not as robust as for solar mass black holes and less is understood about the accretion disks, whose properties may disrupt the superradiance process [65,69]. Quadratically coupled ULDM will get an effective mass from interaction with the accretion disk in the vicinity of a black hole, potentially modifying the dynamics of superradiance. In the absence of a detailed analysis of this effect, we restrict the superradiance constraints to couplings where the induced mass from the accretion disk is subdominant to the bare mass, where ρ BH is the density of the accretion disc at the radius of the boson cloud, and Q i is the dilatonic charge of the accretion disk as in [31]. Furthermore, we neglect any self-interactions. Such selfinteractions can disrupt superradiance, potentially weakening the bound further [23,70,71].

Conclusions
In this work, we examined the effect of ultralight scalar DM with quadratic couplings on the predicted helium abundance produced by Big Bang nucleosynthesis. Figure 3 shows the constraints derived in this work. In addition, we also show constraints from the Eöt-Wash and MICROSCOPE experiments, atomic clocks experiments, and projected constraints from the AION and AEDGE matter-wave interferometry experiments. This is the first time that the quark and gluon constraints from WEP searches have been presented for quadratic couplings. In this work, we treated weak freezeout in the full kinetic description and accounted for backreaction from the SM on the evolution of the DM. From the impact of ULDM on BBN, we constrain the couplings of ULDM to electrons, quarks, and photons for DM masses ranging from about 10 −19 eV to 10 −4 eV. This updates the result in [26] that treated weak freezeout in the instantaneous approximation and neglected backreaction from the SM on the evolution of the DM. For a significant range of parameter space, we show that the bounds from BBN are the strongest constraints on quadratically coupled ultra-light dark matter. We show that these BBN bounds are significantly modified by the effects of the backreaction and the full kinematic decoupling studied in this work.