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.

Ultralight dark matter fields are composed of a large number of individual modes with randomized phases. Thus, as noted in the main text, due to the central-limit theorem, such virialized ultra light fields (VULFs) belong to a broad class of Gaussian random fields that are encountered in many subfields of physics. The statistical properties of such fields are fully described by the two-point correlation functions in the time domain or power-spectral densities in frequency space. A VULF-specific discussion is presented in Ref. [1].
Here we consider ultralight dark matter searches in the regime where the observation time T is shorter than the coherence time τ c of dark matter field. Because the VULF oscillations are nearly coherent on these short time scales, the VULF-generated signal s(t) in a detector can be described by an oscillatory function s(t) = γξϕ(t) = γ ξ Φ 0 cos(2πf ϕ t + θ) . (1) Here f ϕ is the field frequency (specific to a particular realization of the VULF as the Compton frequency, f c , is Doppler shifted by the relative motion with respect to the detector, f ϕ = f c + m ϕ v 2 /(2h)), θ is an unknown yet fixed phase, ξ is a constant determined by the details of the experiment, and γ is a coupling strength that is, assuming a null result, to be constrained by the measurement. We fix the amplitude of the DM field to be positive, Φ 0 ≥ 0 and the phase θ of the field to lie in the range [0 − 2π). The coupling strength γ can have an arbitrary sign.
In the conventional approach, the DM field amplitude Φ 0 is fixed to its root-mean-square value Φ DM , defined in the main main body of the paper. We refer to such an approach as deterministic. The key issue is that for a given experimental run T ≪ τ c , the value of the VULF amplitude Φ 0 can substantially differ from Φ DM ; the VULF amplitude randomly fluctuates. This is apparent from the simulation presented in Fig. 1 of the main text. Thus Φ 0 has to be treated as a nuisance parameter in the analysis. We refer to this treatment as stochastic. This has important implications for establishing constraints on γ.
The probability densities for the DM field amplitude in the two approaches can be summarized as The stochastic probability distribution for the field amplitude can be recast into the probability distribution for the local energy density ρ = Φ 2 0 f 2 c /2 (where f c is the Compton frequency, not the Doppler shifted f ϕ ) of the dark matter field, in agreement with the distribution used in Ref. [2].
As an analytically treatable case, we assume that the data stream was acquired on a uniform time grid. Considering the oscillating nature of the signal (1), we will work in frequency space, applying a discrete Fourier transform (DFT) to the data. The data set is assumed to include N measurements, taken over the total observation time T . To make the derivation as transparent as possible, we assume that the intrinsic instrument noise is white with variance σ 2 .
The assumption of white noise, however, is hardly necessary, and the derivation remains valid for the more general case of colored noise. The generalization can be accomplished by replacing N σ 2 with the noise power spectral densitỹ S(f p ), where f p is the DFT frequency. We use DFT notation and definitions of Ref. [1] (see appendix of that paper for a review of Bayesian statistics in frequency space).
Suppose we would like to establish constraints on the coupling strength γ at one of the DFT frequencies f p = f ϕ withd p being the DFT component of the original data. A commonly used test statistic is the excess-power statistic [3] and we will use its Fourier amplitude analog, A p = |d p |/ √ N σ 2 , see Sec. Supplementary Note 3 for a more detailed discussion of the likelihood. Note that θ becomes important when sampling time T is short when compared to 1/f c , which we treat in Sec. Supplementary Note 1 C. The likelihood at f p (we now drop the unnecessary p subscript) is then for a given excess-signal amplitude A s = |s p |/ √ N σ 2 withs p = N γξΦ 0 e iθ /2 being the DFT component of the signal (1) at f p = f ϕ and I 0 being the modified Bessel function of the first kind. This likelihood is general and in agreement with those used in other broadband searches [1,4]. In a more general case this likelihood entails a product over the relevant frequency space bins, and over multiple channels of a network [1,4,5].
The likelihood is the starting point of most, if not all, inference frameworks and here we will demonstrate several approaches using Eq. (4). We approach the problem from a Bayesian perspective and derive the posterior distribution for γ in the stochastic and deterministic cases using both a uniform and Jeffreys' prior in Sec. Supplementary Note 1 A. We discuss the subtle difficulties in the use of the uniform prior for this particular problem and the benefits of utilizing objective priors. In Sec. Supplementary Note 1 B we demonstrate a marginal likelihood approach to compare the two cases, where the results are also briefly discussed for the case of the gradient coupling in Sec. Supplementary Note 1 C.
As discussed in Refs. [6][7][8][9][10], the exclusion limit can depend on the chosen analysis method. We discuss several approaches in the subsequent sections. Crucially in our case, Lindley's paradox (disagreement between Bayesian and frequentist limits), is resolved with the use of the Berger-Bernardo reference prior. It is important to emphasize that, regardless of the chosen analysis, accounting for the stochastic behavior generically relaxes constraints on couplings of VULFs to standard model particles and fields.

A. Bayesian approach
Bayesian inference of some parameter of interest, y, requires prior information, p(y), to derive the respective posterior distribution, p(y D), given the data D. If there is a nuisance parameter z, marginalization of the joint posterior is performed to remove it using Given the likelihood L(D y, z), application of Bayes theorem yields the joint posterior where the denominator is a normalization constant. With the joint prior p(y, z) = p(y)p(z y) and Eq. (5), it is equivalent to applying Bayes' theorem directly to the marginal likelihood with the marginal likelihood defined as For our particular problem, the parameter of interest is the coupling strength γ, the nuisance parameter is Φ 0 , and the observed data are A. We will compare the resulting limits on γ using the two distributions for Φ 0 shown in Eq. (2). The shape of the resulting posterior distribution p(γ A) has immediate implications for placing constraints, also in the event of a DM detection, as the stochastic nature of the VULF affects uncertainty in the value of γ. Here, we assume that the dark matter signal is well below the detection threshold, an assumption consistent with current experimental results. Then the constraint at the 95% credible interval (the Bayesian analog to the confidence interval), γ 95% , can be determined from The correction factor obtained will be a ratio of the constraints determined from the deterministic and stochastic Φ 0 distributions, γ stoch 95% /γ det 95% . For the remainder of this section we derive constraints in terms of the excess signal amplitude A s where the results can be converted into those of γ. We note that working in terms of Fourier amplitude implies our constraint is actually on γ , if we worked in the complex Fourier domain the coupling could be negative and a factor of two would be introduced to Eq. (9).
Applying Bayes theorem, Eq. (6), to our likelihood, Eq. (4) the posterior distribution for the excess signal amplitude, A s is where the exact shape of the posterior, and resulting limits, depend strongly on the choice of prior p(A s ). The requirement of this prior is the main departure from frequentist or purely likelihood based approaches. While a uniform prior can seem a natural and intuitive choice for an unknown model parameter (all values are equally likely), it can lead to a number of issues that can be avoided with the use of objective priors [11,12]. Objective Bayesian statistics derive prior distributions using formal rules instead of basing them off of degree of belief. Importantly, these objective priors can ensure invariance of the inference results under a change of variables (e.g. between power and amplitude), see Ref. [13] for a review. We demonstrate this lack of invariance and a resulting improper posterior when using a uniform prior after a change of variables to power in Sec. Supplementary Note 1 A 2.

Berger-Bernardo method
In this section we demonstrate the Berger-Bernardo reference prior for the case of a separation of the parameters as discussed above, with a parameter of interest and nuisance parameter y and z respectively [13,14]. First, we find the marginal likelihood as defined in Eq. (8) for our particular case using the stochastic This marginal likelihood is then used to calculate the Berger-Bernardo reference prior, which happens to be equivalent (now that there are no more nuisance parameters [14]) to Jeffreys' prior where I(A DM s ) is the Fisher information, giving the prior (please note that this is effectively a prior on γ given the definition of A DM s ) We can now use this prior to solve for our posterior distribution The corresponding 95% credible interval, Eq. (9), is found assuming A is at the experimental detection threshold, see Sec. Supplementary Note 1 B, A dt = ln (1/α) where α = 0.05 is the type-I error (a false negative: i.e., rejecting the hypothesis when it is actually true) For the deterministic case we use the other distribution of Eq. (2) which is equivalent to setting A s → A DM s in the original likelihood, Eq. (4). Then the deterministic posterior distribution is found using the same Berger-Bernardo reference prior (as we note below, this is technically incorrect but has little effect on the results) Please see Fig. 2 in the main text for a comparison of Eq. (15) and Eq. (17).
Numerical integration under the same assumptions as before, A = A dt , yields a constraint at the 95% credible interval of Comparing the deterministic result with Eq. (16) derived in the stochastic treatment, we find We note that this factor is sensitive to the choice of A and would be reduced if A = ⟨A⟩ rather than A dt was chosen. Rigorously, one should calculate the Berger-Bernardo reference prior for the deterministic likelihood (which would be Jeffreys' prior, Eq. (12), using Eq. (4) with A s → A DM s ). However, this prior does not have an analytic form and the deterministic constraint is relatively insensitive to the choice of prior. Using the aforementioned prior, the stochastic Berger-Bernardo prior, or a uniform prior all yield constraints that agree within 6.3% of Eq. (18). For completeness, we note that using a uniform prior, numerical integration can be avoided resulting in the deterministic posterior This shows that the choice between the three priors has minimal effect on the constraint when using the deterministic model for field amplitude. However, as will be clear in the next section, this insensitivity to the chosen prior is not true for the stochastic case.

Stochastic case for uniform and objective priors
While the deterministic case is insensitive to the choice of prior, we will show that this is not true for the stochastic case. Repeating the calculation done in Eq. (15) to calculate the stochastic posterior but with a uniform prior, yields a constraint with numerical integration of and a correction factor of ∼ 13 which is considerably different from the result using the Berger-Bernardo method, Eq. (19). A serious problem with the uniform prior is revealed by analyzing the effect of a change of variables on the resulting constraint. In order for the results to be meaningful and consistent, the constraint should be independent of the choice to work in amplitude or power (as is the case in the frequentist approach discussed in the next section). In particular, using the uniform prior after a change of variables to excess power, P = |d p | 2 /(N σ 2 ), where we again note that this is effectively a prior on γ with P DM s = P s Φ0=ΦDM = N γ 2 ξ 2 Φ 2 DM /(4σ 2 ). Using this prior yields the constraint agreeing within a fraction of a percent with the result obtained for the analysis using amplitude. For this reason, we find that the effect of the stochastic nature of the VULF on experimental constraints is most reasonably estimated in the Bayesian framework by using the objective prior. Furthermore, as we show in the next section, the effect of stochastic fluctuations of the VULF on constraints analyzed using a frequentist approach yields results similar to that obtained in the Bayesian framework employing the objective prior in Sec. Supplementary Note 1 A 1.

B. Frequentist approach using marginal likelihood
A frequentist constraint can be obtained via hypothesis testing and using the Type-I error (or false alarm rate), α, and Type-II error (or false detection rate), 1 − β (see Ref. [15] or more recent editions). Type-I error is the rejection rate of the null hypothesis when it is true while the Type-II error is the acceptance (or non-rejection) of the null hypothesis when an alternative (signal) hypothesis is true. We illustrate a critical (exclusion) region and threshold that matches these two errors in Fig Our background only hypothesis is the likelihood, Eq. (4), with A s = 0 giving a detection threshold at a confidence level CL = 1 − α The goal is to ensure that false exclusion of the true signal value, 1 − β, does not occur more often than the Type-I error, 1 − β ≤ α, for any signal strength above a certain threshold. The threshold signal strength, A th s , at 1 − β = α determines the constraint at a confidence CL (or ≲ α for discontinuous distributions). So, for the signal plus background case which yields a limit at the 95% CL (α = 0.05) in terms of γ of This limit agrees with the Bayesian approach within 5% using the Berger-Bernardo prior and within 1% using the uniform prior.
To treat the stochastic case we repeat the above calculation using the marginal likelihood, Eq. (11). There is some debate behind treating the likelihood as a PDF in this fashion [12], but it is occasionally done in frequentist approaches [16] and produces reasonable results in our case.
Repeating the previous derivation replacing L with L m yields the constraint This result is within 3% of the Bayesian approach with the Berger-Bernado prior, but is ∼ 4 times smaller than that using a uniform prior. This constraint agrees with a straightforward Monte Carlo (MC) approach that also yielded reasonable fits to the distributions used in this section.

C. Gradient coupling and DC sensitivity
For the gradient coupling, the particular sensitivity axis, time of the experiment, and assumptions made mean that we are unable to calculate a uniform correction factor to apply across the currently published constraints. Many experimental constraints have neglected the power present in the galactic-rest frame (or perpendicular to the lab frame motion), typically setting it to zero. Here we show that under some reasonable assumptions the correction factor is roughly the same as that to the scalar coupling. However, for experiments that neglected the rest-frame contribution, we recommend a reanalysis of their data.
The signal for a gradient-coupled experiment can be modeled as with γ the appropriate coupling constant and ⃗ β the unit sensitivity vector. Similar to Eq. (1), for T ≪ τ c , we can describe a component of the gradient as (s(t)) j ≈ γΨ j cos(2πf ϕ t + θ), where Ψ j is the gradient field amplitude (analogous to Φ 0 ) and we have assumed the sensitivity axis and galactic wave vector are parallel.
Using the results of Ref. [17] we can obtain a distribution for the gradient field amplitude and use it to treat the problem in a similar fashion as in the main text or, as we show here, analogously to Sec. Supplementary Note 1 B. For more details and definitions behind the following equations please see Ref. [17]. The distribution for one component of the gradient (see Eq. 4.12 of Ref. [17]) is which is similarly (to the scalar amplitude) Rayleigh distributed. The shape parameter of the distribution for a particular coordinate,Γ j , is defined asΓ with the Fourier amplitudes, A ⃗ k , defined in the lab frame as where k 0 is set to agree with the SHM virial velocity and the lab velocity is chosen to be along the z direction. We see thatΓ is the quantity of interest that determines the relative anisotropy for a detector along the different directions. For either axis perpendicular to the lab's motion,Γ x =Γ y , we get and for the direction parallel to the motioñ We can now restate the problem for the gradient field amplitude, Ψ, in the same way as Eq. (2) above. Along the z-direction, Ψ 0 (we drop the subscript to simplify the equations) has the following deterministic and stochastic distributions (see Eq. (2) above for the scalar case) We assume the deterministic gradient field amplitude is Ψ DM = √ 2ρ DM v lab and that v lab ≈ v vir . The factor of 5/4 is introduced since the deterministic assumption neglects the rest-frame contribution (the power proportional to k 2 0 ). We compare the scalar and gradient field amplitude distributions in Fig. 2. We can now see it is straightforward to calculate the correction factor using L m A A DM s and the approach from the previous section with A DM s defined in the context of the gradient coupling and rescaled by 5/4. This yields a correction factor for the gradient coupling under these simplifying assumptions to be This correction factor largely depends on the initial assumptions made in the experimental analysis and the detector orientation. Additionally, depending on the frequency regime and duration of the experiment, modulation of the gradient field power will occur at the earth's daily rotation and orbital frequency (one can include this in the k lab term above). We have found that the relative scaling between k 2 0 , k 2 lab and the overall factor can vary depending on the chosen field description, for example comparing the extended results of Refs. [1,17,18].

Constraints below frequency resolution (1/T )
In this section we briefly discuss marginalization over the unknown phase θ of the VULF and, in particular, the effect of this marginalization on searches for VULFs with frequencies below 1/T . This problem was the subject of a recent debate in the literature related to the CASPEr-ZULF-Comagnetometer limits originally presented in Ref. [19]. The comment [20] on the published article [19] correctly pointed out that the limits for frequencies below the 1/T frequency resolution should not improve compared to limits for frequencies above the 1/T frequency resolution. This is true and the mistake in Ref. [19] came from not marginalizing over the uniform distribution of phase, θ in Eq. (31) above, and erroneously assuming the average signal amplitude in this limit ∝ 1 2π |cos(θ)|dθ = 2/π, a deterministic approach to θ similar in principle to the erroneous deterministic approaches to VULF amplitudes and velocities [21]. However, after accounting for this unknown phase as discussed both in the present work and in the reply [21], the authors of the comment [20] still disagreed with our result that the constraint levels off for arbitrarily low frequencies.
The resolution of this debate is rooted in how constraints on couplings of VULF dark matter to standard model particles and fields should be interpreted. Even when assuming that VULF dark matter is well described by the SHM, at a particular time and place, the experimentally observable signal is affected by a number of unknown VULF parameters beyond the coupling strength such as the phase, velocity, and amplitude. Even if the duration T of an experiment is too short to broadly sample the distribution of the unknown VULF parameters, the experiment still excludes the existence of VULF dark matter over some finite volume of a multi-dimensional space of the unknown parameters. To condense the results of a VULF dark matter search into a two-dimensional plot of coupling strength vs. frequency, it is necessary to marginalize over the other unknown VULF parameters. While the authors of the comment [20] object to this marginalization, arguing that for particular phases θ the constraints on couplings could be much weaker, we argue that it is reasonable to plot constraints on couplings ruled out over 95% of the multidimensional unknown VULF parameter space. Furthermore, as we note here, the authors of the comment [20] have treated the amplitude and velocity of the VULF dark matter as deterministic in their own experiment [22]. For particular (but unlikely) VULF velocities and amplitudes, the experiment [22] is insensitive to VULF couplings. In order to constrain VULF dark matter couplings in the regime where T is shorter than either the coherence time or period of the VULF, it is necessary to carry out some form of marginalization over the unknown VULF parameters.
The derived constraints on VULF dark matter presented in Ref. [19] take advantage of the daily modulation provided by earth's rotation. Any gradient-coupled DM signal with frequency f ϕ below the 1/T resolution will be upconverted to the 1/day frequency (f E ) given by Earth's rotation. Thus, instead of trying to resolve frequencies below the 1/T resolution, the authors constrain all frequencies below that threshold by looking only at the f E bin that is within the detection bandwidth. The constraint levels off (becomes constant) for all frequencies below ∼ 1/(25T ) when 2πf ϕ T ≈ 0 and the signal amplitude is dominated by the |cos θ| distribution, see Eq. (31). The height of this constant level is determined by the chosen confidence level, as there are a small range of phases that give vanishing signals for arbitrarily large coupling.
This approach can be used for all gradient-coupled ("wind"-type) experiments that can resolve the f E frequency. We additionally note that scalar searches (for which this upconversion does not occur) which are sensitive to offsets (the zero-frequency component) can also draw similar constraints below the 1/T frequency threshold. There is a caveat that any signal present at the f E frequency would not allow determination of the actual f ϕ , only that it is smaller than 1/T . Additionally, differentiating systematic errors from potential VULF dark matter signals is difficult for upconverted frequencies; one approach would be to confirm the expected yearly modulation in the event of a detection.

D. Blueice Monte Carlo
Here we employ a method widely used in astroparticle physics [23] to compare to the Bayesian and frequentist approaches described in the previous sections. In the case encountered here we find that Wilks theorem, the use of a chi-squared distribution for the log likelihood ratio test statistic, does not hold. The reason for this is in part due to the fact that we are not using the profile likelihood but marginalize over the amplitude fluctuation by drawing from the Rayleigh distribution in each MC realization. In a situation like this, one has to resort to doing a so-called Neyman construction, see e.g. discussion in Ref. [16].
The decisive feature of a Neyman construction is that it provides the statistically desired properties (in this case the coverage of the confidence intervals). It uses the distribution of the test statistic as estimated by many MC simulations instead of the asymptotic distribution used in Wilks theorem. For a given signal strength γ, the teststatistic distribution is obtained by running many MC simulations at that signal strength. The critical value is the value of the test statistic at the 90th (or 95th) percentile of the distribution for a test at the 90% (or 95%) confidence level. This procedure is repeated for different signal strengths such that a critical value curve as a function of the signal strength is obtained. The upper limit on the signal strength for a given measurement (or MC simulation) is then given by the value at which the likelihood crosses the critical value curve.
The correction factor is finally not determined from a single experiment but rather from an ensemble of experiments over which the so called sensitivity (median limit in case of no signal) is calculated. The correction factor is then defined analogous to the case of the Bayesian analysis as: where u is a large number of upper limits at a certain confidence level. Initial results from the simulations at the 95% confidence level are similar to those discussed above when including only the field amplitude distribution, see Sec. 6 of Ref. [24] for simulation details. Work on implementing the correct distributions and signal model is ongoing.

Supplementary Note 2. MODIFIED CONSTRAINT PLOT
Here we present a more detailed parameter exclusion plot for the scalar dark matter coupling where we apply a correction factor based on the analysis method used by the respective experiment. We do not do the same for the gradient coupling per our discussion in Sec. Supplementary Note 1 C.
The scalar coupling has only the field amplitude parameter to account for in the correction factor as treated in the main text. We use the result from the marginal likelihood of γ stoch 95% /γ det 95% = 2.7 from Sec. Supplementary Note 1 B for the experiments using a frequentist based analysis and 3.0 from Eq. (19) for the experiment using a Bayesian framework [25]. Details on the data for the dual rubidium and cesium cold atom fountain FO2 at LNE-SYRTE, black lines, can be found in Ref. [25], for the dysprosium spectroscopy, red line, in Ref. [26], clock comparisons, orange lines, in Ref. [27], and from the global network of optical atomic clocks, blue line, in Ref [28]. The gray dashed line, m ϕ ≲ 10 −22 eV, indicates where the assumption that the field with mass m ϕ makes up all the dark matter should be relaxed [29].

Supplementary Note 3. LIKELIHOOD
Here we demonstrate the change of notation from Ref. [1] to that used in this paper, and discuss the change of parameters between complex Fourier components, amplitude, and power. Starting with the signal, Eq. (1), and with d p the DFT component of the original data, the likelihood at f p is given by wheres p = N γξΦ 0 e iθ /2 is the DFT component of the signal (1). See Appendix of Ref. [1] for a review. For the phase, one can marginalize over θ as done in Ref. [30], however the resulting likelihood can only be used for Bayesian approaches, see discussion in Ch. 2 and Appendix A of Ref. [31]. To remain consistent with frequentist based approaches we can instead do a change of variables to excess amplitude A = d p / √ N σ 2 . First we rewrite Eq. (40) in polar coordinates usingd p = |d p | exp −iϕ and dd p = |d p |d|d p |dϕ as and then do a vector to scalar, (|d p |, ϕ) → A, change of variables with excess signal amplitude A s = |s p |/ √ N σ 2 and I 0 is the modified Bessel function of the first kind.
It is also common to work with the power spectral density [3] and excess power, P = |d p | 2 /(N σ 2 ), where another change of variables yields L P P s = e −(P +Ps) I 0 2 P P s , with P s defined accordingly. See Ref. [3] for the original and alternate derivation of Eq. (45).

Supplementary Note 4. COHERENCE TIME
Another topic is the presence of various definitions of coherence time used throughout the literature. The most common definition used is τ c ≡ f c v 2 vir /c 2 −1 and is the one used in this paper. It is a lower bound derived by considering our motion through the spatial gradients of the field [32,33] given by λ, the De Broglie wavelength, where we have used the Compton frequency f c = m ϕ c 2 /h and v vir ≈ 10 −3 c. However, some definitions differ by factors around ∼ 2π [1]. A factor of 2 can occur from using the kinetic energy mv 2 /2 instead of mv 2 as in Ref. [32,33] or in using 10 6 as the quality factor and corresponding definition of coherence time. These differences in coherence time do not have a significant effect on the results in this paper nor previous published constraints. The strongest effect would be the shift of the upper bound on m a of the CASPEr-ZULF-Sidebands constraint. Since the experiment is incoherently averaging past this point the constraint significantly weakens and was not reported. In addition, the fuzzy definition of coherence time also depends on the relative velocity, making it a distributed quantity.
Regarding this upper bound, a potentially more precise T max rather than the coherence time can be derived (where the regime discussed in this paper would go from T ≪ τ c to T ≪ T max ) considering the lineshape derived in Ref. [1]. Using the derived full width at half maximum ∆ω ≈ 2.5/τ c and setting it equal to the DFT frequency step δω = 2π/T gives T max ≈ 2πτ c /2.5, or with the definition of coherence time used in this paper: T max ≈ τ c /2.5. Thus it is safe to assume the data from the experiments shown are all in the T ≪ T max regime, regardless of the definition of coherence time.

Supplementary Note 5. RAYLEIGH DISTRIBUTION
As shown in Eq. (1), Φ 0 is the observed field amplitude, and in the case of total interrogation time being much less than the field's coherence time, T ≪ τ c , the only field amplitude observed during a measurement. One can relate the local dark matter density ρ DM to the field's energy density Φ 2 0 f 2 ϕ /2 = ρ DM as done for example in Ref. [32]. The problem with this is that Φ 0 is actually a Rayleigh distributed value (approximately, see the following subsection), thus the correct expression is We can define Φ DM as the value satisfying Eq. (47) analogous to previous assumptions of Φ 0 = Φ DM in the deterministic approach. Using Eq. (47) we can solve for the required σ of the Rayleigh distribution p(Φ 0 ).
Note that the average value of this distribution ⟨Φ 0 ⟩ = Φ DM √ π/2 is not simply Φ DM . Approximately 63% of given field realizations will have a field amplitude smaller than Φ DM . The distribution compared to Φ DM is shown in Fig. 4.

A. Simulated energy distribution
The energy density distribution, Eq. (3), and amplitude distribution, Eq. (49), are only exact under the random phase model. The distribution depends on a number of parameters such as galaxy mass and position within the galaxy [17]. As shown in Ref. [17] simulations of these distributions indicate a slightly heavier tail with a best fit distribution ρ DM p(ρ) = 0.1e −0.42ρ/ρDM + 0.9e −1.06ρ 2 /ρ 2 DM . We compare the two distributions below, where a slightly heavier tail is evident in the log-log plot. Marginalizing over this distribution instead of the random phase model, following the steps in Sec. Supplementary Note 1 B, only changes the derived constraint by ≈ 3%. Further investigation of how gravitational interactions impact these distributions for our particular conditions and if their characteristic shape can be utilized for inference is ongoing.