Stochastic fluctuations of bosonic dark matter.

Numerous theories extending beyond the standard model of particle physics predict the existence of bosons that could constitute dark matter. In the standard halo model of galactic dark matter, the velocity distribution of the bosonic dark matter field defines a characteristic coherence time τc. Until recently, laboratory experiments searching for bosonic dark matter fields have been in the regime where the measurement time T significantly exceeds τc, so null results have been interpreted by assuming a bosonic field amplitude Φ0 fixed by the average local dark matter density. Here we show that experiments operating in the T ≪ τc regime do not sample the full distribution of bosonic dark matter field amplitudes and therefore it is incorrect to assume a fixed value of Φ0 when inferring constraints. Instead, in order to interpret laboratory measurements (even in the event of a discovery), it is necessary to account for the stochastic nature of such a virialized ultralight field. The constraints inferred from several previous null experiments searching for ultralight bosonic dark matter were overestimated by factors ranging from 3 to 10 depending on experimental details, model assumptions, and choice of inference framework.

Introduction

It has been nearly ninety years since strong evidence of the missing mass we label today as dark matter (DM) was revealed[1], and its composition remains one of the most important unanswered questions in physics. There have been many DM candidates proposed and a broad class of them, including scalar (dilatons and moduli[2-5]) and pseudoscalar particles (axions and axion-like particles[6-11]), can be treated as an ensemble of identical bosons, with statistical properties of the corresponding fields described by the standard halo model (SHM)[12,13]. In this work, our model of the resulting bosonic field assumes that the local DM is virialized and neglects non-virialized streams of DM[14], Bose–Einstein condensate formation[15-18], and possible small-scale structure such as miniclusters and axion stars[19-21]. To date, it is typical to ignore such DM structure when calculating experimental constraints, and within this isotropic SHM DM model, we demonstrate the general weakening of inferred constraints due to the statistical properties of the virialized ultralight field (VULF)[21-24]. We note that some astrophysical and cosmological simulations can and do resolve these stochastic properties[25,26], however in this paper we discuss their impact on inferences drawn from direct detection experiments.

During the formation of the Milky Way the DM constituents relax into the gravitational potential and obtain, in the galactic reference frame, a velocity distribution with a characteristic dispersion (virial) velocity vvir ≈ 10−3c and a cut-off determined by the galactic escape velocity. Following Refs. [27,28] we refer to such virialized ultralight fields, ϕ(t, ), as VULFs, emphasizing their SHM-governed stochastic nature. Neglecting motion of the DM, the field oscillates at the Compton frequency fc = mϕc2h−1. However, there is broadening due to the SHM velocity distribution according to the dispersion relation for massive nonrelativistic bosons: fϕ = fc + mϕv2(2h)−1. The field modes of different frequency and random phase interfere with one another resulting in a net field exhibiting stochastic behavior. The dephasing of the net field can be characterized by the coherence time [29]. We note that there is some ambiguity in the definition of the coherence time, up to a factor of 2π, and adopt that which was used in the majority of the literature. See the discussion in Supplementary Note 4.

While the stochastic properties of similar fields have been studied before, for example in the contexts of statistical radiophysics, the cosmic microwave background, and stochastic gravitational fields[30], the statistical properties of VULFs have only been explored recently. The 2-point correlation function, , and corresponding frequency-space DM “lineshape” (power spectral density, PSD) were derived in Ref. [28], and rederived in the axion context by the authors of ref. [31]. While refs. [28,31] explicitly discuss data-analysis implications in the regime of the total observation time T being much larger than the coherence time, T ≫ τc, detailed investigation of the regime T ≪ τc, until now, has been lacking (although we note that ref. [31] includes a brief discussion of the change in sensitivity due to coherent averaging for this regime in their Appendix E). Note that a preprint of this paper has been available online since 2019, and multiple experimental groups have already used it to correct their exclusion limits for stochastic fluctuations or noted the effect[32-41].

We focus on this regime, T ≪ τc, characteristic of experiments searching for ultralight (pseudo)scalars with masses ≲ 10−13 eV[42-48] that have field coherence times ≳ 1 day. This mass range is of significant interest as the lower limit on the mass of an ultralight particle extends to 10−22 eV and can be further extended if it does not make up all of the DM[49]. Additionally, there has been recent theoretical motivation for “fuzzy dark matter” in the 10−22–10−21 eV range[23,49-53], and the so-called string “axiverse” extends to 10−33 eV[54]. Similar arguments also apply to dilatons and moduli[55].

Here, we show that for experiments operating in the T ≪ τc regime it is incorrect to assume a fixed value of Φ0 when inferring constraints on the coupling strength of bosonic DM to standard-model particles. The constraints inferred from several previous null experiments searching for ultralight bosonic DM were overestimated by factors ranging from 3 to 10 depending on experimental details, model assumptions, and choice of inference framework.

Results

Model of bosonic dark matter and amplitude distribution

Figure 1 shows a simulated VULF field, illustrating the amplitude modulation present over several coherence times. At short time scales (≪τc), the field coherently oscillates at the Compton frequency, see the inset of Fig. 1, where the amplitude Φ0 is fixed at a single value sampled from its distribution. An unlucky experimentalist could even have near-zero field amplitudes during the course of their measurement.

Fig. 1

Simulated VULF based on the approach in ref. [28] with field value ϕ(t) and time normalized by ΦDM and coherence time τc, respectively.

The inset plot displays the high-resolution coherent oscillation starting at t = 0.

On these short time scales, the DM signal s(t) exhibits a harmonic signature,where γ is the coupling strength to a standard-model field and θ is an unknown phase. Details of the particular experiment are accounted for by the factor ξ. In this regime, the amplitude Φ0 is unknown and yields a time-averaged energy density . However, for times much longer than τc the energy density approaches the ensemble average determined by . This field oscillation amplitude is estimated by assuming that the average energy density in the bosonic field is equal to the local DM energy density ρDM ≈ 0.4 GeV/cm3, and thus .

The oscillation amplitude sampled at a particular time for a duration ≪τc is not simply ΦDM, but rather a random variable whose sampling probability is described by a distribution characterizing the stochastic nature of the VULF. Until recently, most experimental searches have been in the mϕ ≫ 10−13 eV regime with short coherence times τc ≪ 1 day[56-70]. However, for smaller boson masses it becomes impractical to sample over many coherence times: for example, τc ≳ 1 year for mϕ ≲ 10−16 eV. Assuming the value Φ0 = ΦDM neglects the stochastic nature of the bosonic dark matter field[42-48].

The net field ϕ(t) is a sum of different field modes with random phases. The oscillation amplitude, Φ0, results from the interference of these randomly phased oscillating fields. This can be visualized as arising from a random walk in the complex plane, described by a Rayleigh distribution[31]analogous to that of chaotic (thermal) light[71]. This distribution implies that ≈63% of all amplitude realizations will be below the r.m.s. value ΦDM. Equation (2)[31] is typically represented in its exponential form[72] (Supplementary Note 3), and is well sampled in the T ≫ τc regime. However, this stochastic behavior should not be ignored in the opposite limit. Simulations of galactic p(Φ0) distributions for fuzzy dark matter show slightly heavier tails than the random phase model[35], Eq. (2), but these differences have a negligible effect on the results of this paper as we discuss in Supplementary Note 5.

Establishing constraints on coupling strength

We refer to the conventional approach assuming Φ0 = ΦDM as deterministic and approaches that account for the VULF amplitude fluctuations as stochastic. To compare these two approaches we choose a Bayesian framework and calculate the numerical factor affecting coupling constraints, allowing us to illustrate the effect on exclusion plots of previous deterministic constraints[42-48]. It is important to emphasize that different frameworks to interpret experimental data than presented here can change the magnitude of this numerical factor[73-76], see Supplementary Note 1 for a detailed discussion. In any case, accounting for this stochastic nature will generically relax existing constraints as we show below.

We follow the Bayesian framework[77] (see application to VULFs in ref. [28]) to determine constraints on the coupling-strength parameter γ. Bayesian inference requires prior information on the parameter of interest to derive its respective posterior probability distribution, in contrast to purely likelihood-based inference methods. The central quantity of interest in our case is the posterior distribution for possible values of the coupling constant γ, derived from Bayes theorem,The left-hand side of the equation is the posterior distribution for γ, where D represents the data, and the Compton frequency fϕ is a model parameter. is the normalization constant, and the likelihood is the probability of obtaining the data D given that the model and prior information, such as those provided by the SHM, are true. The integral on the right-hand side accounts for (marginalizes over) the unknown VULF amplitude Φ0, which we assume follows the Rayleigh distribution described by Eq. (2). For the choice of prior p(γ, Φ0) we use what is known as an objective prior[78]: the Berger–Bernardo reference prior[79]. Note that this approach is equivalent to starting with the marginal likelihood and using Jefferey’s prior to calculate the posterior[80]. See details in Supplementary Note 1.

Results from Bayesian inference are sensitive to the choice of prior[79], and we find better agreement with frequentist-based approaches when using an objective prior rather than a uniform prior p(γ) = 1 (as shown in Supplementary Note 1). Additionally, the uniform prior yields constraints that are noninvariant under a change of variable.

It is important to note that experiments searching for couplings of VULFs to fermion spins (axion “wind” searches) are sensitive to the projection of the field gradient onto the sensitive axis of the experiment. Due to this directional sensitivity, the derived coupling strength strongly depends on specific experimental conditions. However, under some reasonable assumptions discussed in Supplementary Note 1, the correction factor is similar in size to the scalar case considered here. Axion-wind experiments can also utilize the daily modulation of this projection, due to rotation of the Earth, to search for signals with an oscillation period much longer than the measurement time T ≪ 1/fϕ. The unknown initial phase θ of the VULF sets the amplitude of this daily oscillation and also needs to be marginalized over. We discuss these topics in Supplementary Note 1, relevant for the experiments[42-45], and focus solely on stochastic variations of the scalar field amplitude, Φ0, here.

Using the posterior distribution, p(γ∣D, fϕ, ξ), one can set constraints on the coupling strength γ. Such a constraint at the commonly employed 95% confidence level (CL), γ95%, is given by

The posteriors in both the deterministic and stochastic treatments are derived in Supplementary Note 1. In short, the two posteriors differ due to the marginalization over Φ0 for the stochastic case, see the integral of Eq. (3). Assuming white noise of variance σ2 and that the data are in terms of excess amplitude A (observed Fourier amplitude divided by expected noise, an analog to the excess power statistic) we can derive the posterior for excess signal amplitude As. The posteriors areHere , I0(x) is the modified Bessel function of the first kind, and p(As) is effectively the prior on γ. In Fig. 2, we plot the normalized posteriors assuming A at the 95% detection threshold and using Berger–Bernardo reference priors for p(As); we compare other choices of prior in Supplementary Note 1. The derivation relies on the discrete Fourier transform for a uniform sampling grid of N points and the assumptions of the uniform grid and white noise can be relaxed[28].

Fig. 2

Posterior distributions for the coupling strength γ in the deterministic and stochastic treatments, Eqs. (5) and (6), respectively.

Due to the fat-tailed shape of the stochastic posterior one can clearly see the 95% limit is larger with . The assumed value of the data is at the 95% detection threshold (see text).

Examination of Eqs. (5), (6) and Fig. 2 reveals that the fat-tailed stochastic posterior is much broader than the Gaussian-like deterministic posterior. It is clear that for the stochastic posterior, the integration must extend considerably further into the tail, leading to larger values of γ95% and thereby to weaker constraints, . Explicit evaluation of Eq. (4) with the derived posteriors results in a relation between the constraintswhere the numerical value of the correction factor depends on CL and assumed value of A (the factor increases for higher CL and decreases for smaller A).

This correction factor becomes ≈10 when derived using a uniform prior, as shown in Supplementary Note 1. However, the result obtained with the uniform prior is not invariant under a change of variables (e.g., from excess amplitude to power). Additionally, using the objective prior yields better agreement with frequentist-based results of a factor ≈2.7. For the gradient coupling of pseudoscalar particles, the directional sensitivity, deterministic assumptions, and initial phase of the field (when relevant) can further impact this factor as discussed in Supplementary Note 1.

Discussion

Ultralight DM candidates are theoretically well-motivated and an increasing number of experiments are searching for them. Most of the experiments with published constraints thus far are haloscopes, sensitive to the local galactic DM and affected by Eq. (7). However, experiments that measure axions generated from a source, helioscopes, or new-force searches, for example, do not fall under the assumptions made here. We illustrate how the existing constraints have been affected in Fig. 3 and provide a more detailed exclusion plots for dilaton couplings[46-48] in Supplementary Note 2.

Fig. 3

The modified constraint, green and blue lines, based on the stochastic approach compared to previous laboratory constraints, gray line, based on the deterministic approach for the dilaton coupling strength de[41, 46–48].

The green and blue lines illustrate the importance of the choice of prior for a Bayesian approach. Supplementary Fig. 3 provides a detailed exclusion plot.

To interpret the results of an experiment searching for bosonic DM in the regime of measurement times smaller than the coherence time, stochastic properties of the net field must be taken into account. An accurate description accounts for the Rayleigh-distributed amplitude Φ0, where the variation is induced by the random phases of individual virialized fields. Accounting for this stochastic nature yields a correction factor of ≈2.7−10, relaxing existing experimental bosonic DM constraints in this regime. In the event of a bosonic DM discovery, the stochastic properties of the field would result in increased uncertainty in the determination of coupling strength or local average energy density in this regime.

It is important to note that observational knowledge of the energy distribution of DM[81] could constrain the stochastic behavior of the amplitude (or energy density). The smallest features observed so far are on the order of ≈0.1 kpc[82] (corresponding to a coherence length of a boson of mass mϕ ≈ 10−22 eV), however the analysis in ref. [82] performs angular averages which would suppress the stochastic variation discussed in this paper.