Dark Matter Interferometry

The next generation of ultralight dark matter (DM) direct detection experiments, which could confirm sub-eV bosons as the dominant source of DM, will feature multiple detectors operating at various terrestrial locations. As a result of the wave-like nature of ultralight DM, spatially separated detectors will each measure a unique DM phase. When the separation between experiments is comparable to the DM coherence length, the spatially-varying phase contains information beyond that which is accessible at a single detector. We introduce a formalism to extract this information, which performs interferometry directly on the DM wave. In particular, we develop a likelihood-based framework that combines data from multiple experiments to constrain directional information about the DM phase space distribution. We show that the signal in multiple detectors is subject to a daily modulation effect unique to wave-like DM. Leveraging daily modulation, we illustrate that within days of an initial discovery multiple detectors acting in unison could localize directional parameters of the DM velocity distribution such as the direction of the solar velocity to sub-degree accuracy, or the direction of a putative cold DM stream to the sub-arcminute level. We outline how to optimize the locations of multiple detectors with either resonant cavity (such as ADMX or HAYSTAC) or quasistatic (such as ABRACADABRA or DM-Radio) readouts to have maximal sensitivity to the full 3-dimensional DM velocity distribution.


I. INTRODUCTION
Cold, bosonic dark matter (DM) candidates with masses much smaller than the eV scale have macroscopic occupation numbers and may be described in the solar vicinity by classical fields. Two well-studied DM candidates in this category, which we broadly refer to as wave-like DM, are axions [1][2][3][4][5][6][7][8][9][10][11] and dark photons [12][13][14]. Wave-like DM candidates require distinctive experimental techniques for discovery that take advantage of their spatial and temporal coherence (see, e.g., [15]). The spatial coherence length of the DM waves, λ c , and the coherence time τ are given by 1 [16] where v 0 parameterizes the DM velocity dispersion,v is the mean velocity, and m DM is the DM mass. In the solar neighborhood we expectv ∼ v 0 ∼ 10 −3 for the bulk of the DM, in natural units, such that the coherence length is around 10 3 times the Compton wavelength, and the coherence time is around 10 6 times the oscillation period for the DM wave. In this work we show that multiple phasesensitive wave-like DM detectors separated by distances of order λ c may join their data -through a process we refer to as "DM interferometry" -to measure properties of the DM phase-space distribution that are inaccessible to single experiments operating in isolation. 1 We briefly review both concepts in App. A.
Many axion and dark-photon detection strategies already leverage the axion coherence time as a "quality factor" that amplifies the DM signal in the experiment. For example, axion haloscopes [17][18][19][20][21][22] use a resonant cavity with a strong static magnetic field to convert axion DM into electromagnetic cavity modes, which build up coherently over the DM coherence time; in this setup the DM Compton wavelength is of order the size of the experiment. Experiments operating in the quasistatic regime (where the DM Compton wavelength is much larger than the experiment) -including searches for the axion-photon coupling [23][24][25][26][27][28], axion interactions with nuclear spins [29], or dark photons [30] -aim to detect a time-varying magnetic flux through a pickup loop, which can build up coherently in a lumped-element circuit [31,32].
In this paper, we explore the phenomenology of spatial coherence for wave-like DM by exploiting spatiallyseparated detectors that probe the same DM field. It is straightforward to understand why multiple detectors offer unique insights for wave-like DM. Generically, the wave-like DM field may be written as a(x, t) = a 0 cos(ωt − k · x + φ), where ω is the oscillation frequency, k is the wave vector, φ is a random phase, and a 0 is the amplitude. 2 If the DM wave is traveling in the directionk with speed v 1, then ω ≈ m DM (1 + v 2 /2) Two detectors separated by a vector x12, however, are sensitive to the speed distribution modulated by the k · x12 phase of the DM wave, replacing f (v) with functions F c,s 12 (v) as defined in (3). As the figures demonstrate, the modified speed distributions exhibit daily modulation and carry additional information about the velocity distribution f (v) that would be invisible to a single detector. For this example we take mDM = 25.2 µeV [33], near the window where the HAYSTAC collaboration is searching for axion DM. Taking the Standard Halo Model ansatz for f (v) in (50), we place one detector at a latitude and longitude of (41 • N, 73 • W), and a second instrument ∼ 20 m to the North, corresponding to d ∼ 2λc. A curve is shown for every ten minutes starting from midnight on January 1st of 2020. Note that as F c,s 12 (v) are functions of mDMd, qualitatively similar effects exist for e.g. mDM ∼ 10 −9 eV, in the mass range probed by ABRACADABRA and DM-Radio, for d ∼ 500 km. and k ≈ m DM vk. For a single detector we may always choose coordinates such that x = 0. This means that a single detector is only sensitive to the speed through ω and is not sensitive to the direction of the DM velocity. 3 By contrast, two experiments located at positions x 1 and x 2 will be sensitive to phase factors k · x 1 and k · x 2 . Only one of these can be removed by a coordinate choice, leaving a residual k · x 12 , with x 12 = x 1 − x 2 , which manifestly probes the velocity rather than the speed. The interferometry proposed in this work is directly at the level of the DM field: the effect arises due to the phase difference wave-like DM exhibits between spatially separated points. 4 Indeed, due to the nonzero velocity dispersion v 0 , DM waves are coherent up to distances of order λ c ; as we will show, phase-sensitive data combined from two experiments exhibits maximal modulation when d ≡ |x 12 As we will demonstrate, this opens up striking new signatures, such as a unique daily modulation signal applicable only to wave-like DM with multiple detectors, because the direction of x 12 rotates over a sidereal day with respect to the DM field. 5 The main result of this paper is that interference effects between a pair of detectors separated by a distance x 12 are fully characterized by the modified speed distributions with f (v) the DM velocity distribution. Examples of these distributions at various times throughout the day are shown in Fig. 1 for optimally-separated detectors. If the goal is simply to enhance the total signal reach, we should maximize the constructive interference and take d λ c , in which case F s 12 (v) = 0 and F c 12 (v) = f (v), with f (v) the DM speed distribution. The observation that there is an enhanced sensitivity for an array of detectors located within the DM coherence length has been made previously in Ref. [45]; this is also the basis for the multiplexed cavity setup proposed by ADMX [46]. However, if the goal is to extract information about the full 3-dimensional DM phase space distribution, which encodes e.g. the boost of the Solar System with respect to the Galactic Center as well as possible DM substructure (including the Sagittarius stream [47] and the Gaia Sausage [48,49]), we should take d ∼ λ c as in (2). Ref. [45] points out the possibility of observing this daily modulation effect for experiments separated by distances of order λ c ; here we extend this analysis by focusing on constraining directional parameters in the phase space distribution. We will show that the sensitivity to the phase space information that may be extracted from multiple detectors is comparable to the sensitivity to the initial discovery, since the interference effects have an O(1) effect on the data when d ∼ λ c . As such, in principle these unique signatures could be used to immediately verify a putative axion signal. More optimistically, DM interferometry would allow for the detailed mapping of the local DM phase space distribution after an initial detection.
For concreteness, we focus in this work on the case of axion DM coupled to electromagnetic signals, but our results would apply equally well to scalar and vector DM as long as the readout is proportional to the DM field. Similarly, for simplicity we will present most results for the case of two experiments, but our formalism holds for any number N ≥ 2 of experiments, and we will provide our key results for a general N also. Our results also apply equally well to resonant-type experiments and to broadband-type experiments (such as ABRACADABRA-10 cm [24,25]), so long as the resonant experiments are able to preserve the phase of the data, as opposed to e.g. recording the power directly. One advantage of resonant experiments for wave-like DM, in addition to generically having enhanced sensitivity [31,32], is that putative signal candidates may immediately yield detailed and high-significance studies, since the signal-to-noise ratio rapidly grows with measurement time when frequency-scanning is no longer necessary.
We organize the remaining discussion as follows. In Sec. II, we sketch a derivation for the statistics of the correlated Fourier-transformed data from multiple experiments. A more extensive derivation and discussion is presented in App. B, with some useful orthogonality relations summarized in App. C. In Sec. III, we construct a likelihood function for the axion signal as observed at N experiments, following the formalism of [16]; a practical data-stacking procedure is outlined in App. D. In Sec. IV we use several simplified toy examples to illustrate analytic estimates of uncertainties on parameters of the velocity distribution using the Asimov data seta technique where the asymptotic properties of the data are assumed in order to replace Monte Carlo simulations with analytic estimates (see Sec. III) -and demonstrate that uncertainties on directional parameters in several simple examples with two detectors are minimized for d ∼ 2λ c . We also highlight the important distinction between the de Broglie wavelength and the coherence length for cold but boosted DM substructure. Furthermore, we show that there is a rotational symmetry of the likelihood which can lead to degenerate best-fit parameters for N = 2 experiments. In Sec. V, we extend the likelihood analysis to include daily modulation from the changing detector orientation throughout the day. We also perform analyses of simulated data sets to demonstrate how the likelihood may be implemented in practice to constrain the morphology of the DM phase-space distribution. Using the realistic examples of the Standard Halo Model (SHM) velocity distribution and the Sagittarius Stream, we show how daily modulation breaks the symmetry discussed in Sec. IV and use this to perform parameter estimation using the effect. We conclude in Sec. VI with some practical implications for current and upcoming axion experiments. In App. A we provide a brief review of the coherence length and time.

II. THE STATISTICS OF MULTIPLE DETECTORS
In this section we describe the statistics of an axion DM signal collected by two or more spatially-separated detectors. In particular, while we expect background sources to be generally uncorrelated between detectors, the axion will induce non-trivial cross-correlations indicative of DM interferometry. These correlations will be the source of the additional information available to two or more experiments that we will extract using a likelihood formalism introduced in Sec. III.
We imagine that a given detector, located at a position x, is sensitive to the axion through a time-varying signal Φ proportional to the axion field, The flux Φ is generated by the axion effective current, J a ∼ ∂ t a ∼ m a a, which is the origin of m a in the expression. Accordingly, Φ ∼ κ i J a , revealing κ i as characterizing the individual experimental response to the axion field. In the notation of Ref. [16], we take κ i = A i /ρ DM . The dimensionful constant A i is characteristic of the individual experimental response to the axion field; for instance, in the case of ABRACADABRA [23][24][25], the magnetic flux induced in the pickup loop at the center of the detector is related to the axion field by where g aγγ is the axion-photon coupling, B 0 is the toroidal magnetic field strength, and V B is an effective magnetic field volume associated with the detector. In addition to the experimental factors, A i has been defined to include ρ DM g 2 aγγ , which determines the mean power in the axion field. We assume for simplicity that the detector response A i is purely real and does not include phase delays. Similar expressions are available for other detectors [16]. For our discussion, all that is required is a measurement linear in the axion field, in order to ensure direct access to the axion phase. Measurements intrinsically proportional to a 2 , such as the power in the cavity of an axion haloscope, cannot be directly ported to our formalism. Nevertheless, interferometry can still be performed by these resonant cavity experiments, as long as the phase information is extracted. This may be achieved for example by reading out electromagnetic signals with phase-sensitive amplifiers (e.g., [21,50]).
Ultimately, we envision a set of measurements Φ i of the same axion field, made by N detectors at different spatial locations x i . The correlations between these data sets will arise due to the statistics of the underlying axion field, as we will describe in the following subsections, leaving the full derivation to App. B.

A. Construction of the Axion Field
It is useful to recall the underlying statistics in the axion field that result from its finite velocity dispersion and wave-like nature. In [16] it was shown that we may represent the axion field as seen by a single detector as Here, the sum over j indicates a sum over subsets of particles with speeds in the interval v to v + ∆v. The phase is controlled by ω j = m a 1 + v 2 j /2 and a random contribution φ j ∈ [0, 2π), and further f (v) is the DM speed distribution in the laboratory frame. In Ref. [16] the continuum limit for speeds ∆v → 0 is taken; below (and in App. B) we will generalize to the continuum limit for velocities. In addition to the random phase, the random nature of the axion field is captured in the random variate α j drawn from the Rayleigh distribution p[α] = α e −α 2 /2 . While (5) represents the axion field constructed from the discretized frequency modes specified by the local DM velocity distribution, a more fundamental approach can be understood by considering the local DM field made up of N a axion particles (or wave packets), as detailed in [16]. The enormous occupation numbers characteristic of wave-like DM will then allow us to eventually convert this sum to an integral by taking the N a → ∞ limit; in detail, we should have n DM λ 3 dB 1, where n DM is the DM number density, which is satisfied locally for m a 1 eV. We note that the above construction also assumes DM is a non-interacting wave, which means that selfinteractions should be negligible.
The axion field described in (5) is appropriate for a single detector, but to reveal the effects of DM interferometry we need to extend the description to include the spatial dependence of the DM wave. In particular, the phase will also include a contribution k·x, with k = m a v for a non-relativistic wave. As k depends on the velocity, and not speed, we need to extend the above sum to three independent components, v abc = v ax + v bŷ + v cẑ , where the indexes a, b, c are integers. We may then write where ω abc depends on v abc = |v abc |, and α abc and φ abc are Rayleigh and uniform random variables, respectively, as in (5). Here, (∆v) 3 is a discretization of the 3dimensional velocity, generalizing ∆v for speeds; we will take the continuum limit in App. B. In (6) we have written the axion field in a convenient form for revealing DM interferometry. To reiterate the point, if we measure the axion field at a single location, we can always choose our coordinates such that x = 0. In this case, the velocity information in k is lost, and we are now only sensitive to the speed v = |v| through ω. This collapses f (v) → f (v), and (6) to (5); information about the phase space is lost. However, if we measure the axion field at two locations, an irreducible k dependence remains, and the full velocity information is imprinted in the multi-detector covariance matrix.
We implicitly assume throughout this work that the non-interacting plane-wave superposition in (6) applies for all x. Corrections to this picture should arise from e.g. the gravitational field of the Earth, which would slightly bend the DM trajectories between detectors. However, the DM velocities we consider in this work are much larger than the Earth's escape velocity, and also the detector separations are typically much smaller than the radius of the Earth, so we are justified in neglecting this effect.

B. The Multi-Detector Covariance Matrix
We will now outline how the statistics of the axion field, as described above, lead to a non-trivial covariance matrix in the data collected by N experiments. In this section we will simply state the key results, leaving a derivation to App. B.
We begin by considering the minimal case of N = 1. Suppose that a single experiment takes a time-series of N measurements {Φ n (x) = Φ(x, n∆t)}, with n = 0, 1, . . . , N −1, collected over a time T , so that ∆t = T /N . In order to isolate a signal oscillating at a particular frequency, as we expect the axion to do, we calculate the discrete Fourier transform The transform is indexed by an integer k = 0, 1, . . . , N − 1, which is related to the angular frequency, ω = 2πk/T . We will switch back and forth between talking about frequency ω and wave-number k as convenient. It is convenient to partition the Fourier transform into appropri-ately normalized real and imaginary parts as follows, We can then write the power spectral density (PSD) as, We will assume throughout that T is long enough such that the signal is sufficiently well resolved, i.e. the bandwidth of the Fourier transform 2π/T is much smaller than the width of the signal in frequency space. When specifying Fourier components by frequency as opposed to wave-number we use notation as in S k ΦΦ (x) → S ΦΦ (x, ω). As shown in [16], both R(x, ω) and I(x, ω) are normally distributed with zero mean and variance given by where we have defined v ω = 2ω/m a − 2 as the axion velocity corresponding to a frequency ω, and the speed distribution is defined as This implies, for example, that S ΦΦ (x) is an exponentially distributed quantity, with mean We can understand the velocity dependence by looking back to (5); the signal measured by a detector Φ is proportional to the time-dependent axion axion field, and since a is proportional to f (v), we obtain a power spectrum S ΦΦ proportional to f (v). In any real experiment there will also be background. However, as long as the background is normally distributed in the time domain -as expected for, for example, thermal noise, SQUID flux noise, or Josephson parametric amplifier noise -then both R(x, ω) and I(x, ω) remain normally distributed but with variance where λ B (ω) encapsulates the variance of the potentially frequency-dependent noise from the background sources only.
Note that R k (x) and I k (x) are uncorrelated; in particular, the 2 × 2 covariance matrix for these two quantities is simply This implies that for a single detector, all information about the signal is contained in the PSD S ΦΦ (x). Further, as shown in (14), the location x never enters for N = 1. Even if we chose our coordinates such that k·x = 0, the overall phase remains unphysical as it would vanish when computing the modulus squared in (9). Now let us extend the discussion to the case of interest: data collected by N experiments at positions x i , with i = 1, 2, . . . , N . For each data set, we calculate the real and imaginary parts of the Fourier transform as above. The information collected by all detectors can then be organized into the following 2N dimensional data vector, T . (15) Correlations between the real and imaginary part for any given detector will be identical to the N = 1 case discussed above. However, DM interferometry will reveal itself through non-trivial correlations amongst the different detectors. 6 Indeed, as justified in App. B, d k will be a 2N -dimensional Gaussian random variable with zero mean and a symmetric (2N × 2N )-dimensional covariance matrix given by can be safely neglected. For two detectors with |x 12 | ∼ 50 m and v ∼ 200 km/s, this implies ∆τ 10 −10 s. Typical atomic clocks have timing error of 10 −9 s/day, so the required ∆τ can be achieved by synchronizing the two experiments to an atomic clock over data-taking intervals of about 2.5 hours, which is sufficient for the daily modulation analysis in Sec. V.
Here λ B,i (ω) is the background observed by the i th experiment, and its contribution is purely diagonal. The axion signal, however, induces off-diagonal correlations, which we quantify in terms of with By translation invariance, the entries of the correlation matrix only depend on the relative distances x ij ≡ x i − x j . These expressions then simplify for the correlations amongst a single detector, as F c ii (v) = f (v) and F s ii (v) = 0. But for i = j, the expressions in (18) contain a modulated version of the full velocity distribution, allowing us to extract non-trivial directional information about the velocity distribution f (v) with multiple detectors separated by distances of order the de Broglie wavelength, where the integrand in (18) exhibits maximal variation. We note that the formalism we have developed assumes that the velocity distribution is stationary, or at least varies slowly on timescales compared to the axion coherence time. In Sec. V we will develop a formalism to take into account the daily modulation of f (v) through a joint likelihood over multiple data-taking intervals.

III. A LIKELIHOOD FOR MULTI-DETECTOR AXION DIRECT DETECTION
Having understood the statistics underlying the data collected by multiple detectors, we now outline how to incorporate these lessons into an appropriate likelihood. The likelihood will be a simple generalization of the axion likelihood (generally applicable to wave-like DM) introduced in [16], and we will closely follow their approach. However, unlike in [16], we work explicitly with the data as represented in R k and I k rather than the PSD, as the former notation exposes the full set of multi-detector correlations, as captured by Σ in (16). We will then outline how we can extract information about the parameters of f (v) using this likelihood, exploiting where possible the asymptotic Asimov procedure [51] to determine results analytically. In Sec. IV we will then put the formalism to use in the context of several toy examples designed to highlight where interferometry opens up new avenues, and build intuition for the more realistic scenarios considered in Sec. V.

A. The Multi-Detector Likelihood
As detailed in Sec. II, we imagine we have a data set collected by N experiments, which each perform a time series of N measurements collected at a frequency f = 1/∆t of a quantity Φ ∝ a. The real and imaginary part of the discrete Fourier transform of each experiments data set is constructed according to (8), and then arranged into a single data set d = {d 0 , d 1 , . . . , d N −1 }, with d k as given in (15). We then define a model M with parameter vector θ that has nuisance parameters θ nuis describing the backgrounds in the individual experiments (encapsulated by λ B,i (ω)) and signal parameters θ sig that characterize the axion contribution. For example, θ sig includes g aγγ , m a , and model parameters that describe the DM velocity distribution f (v). Then, as the data set is distributed according to a multivariate Gaussian, the appropriate likelihood is given by, where |Σ k (θ)| is the determinant of the covariance matrix.
The utility of the likelihood function is that it facilitates inferences regarding the signal parameters, θ sig , from the data. The ultimate goal of the axion DM program would be to infer a nonzero value of A, and hence the existence of a coupling between the Standard Model and DM, for example g aγγ . Taking a frequentist approach to that problem, it is useful to define the following test statistic (TS) from the profile likelihood: In each likelihood,θ nuis denotes the value of the nuisance parameters that maximizes the likelihood for the given signal parameters. The TS is then a function of the signal model parameters. In particular, this means that in the first term in (20) the nuisance parameters are uniquely determined at each θ sig point by the values which maximize the log likelihood. The second term in (20) is evaluated on the null model θ sig = 0, which can be achieved by setting the signal strength parameter A to zero. The TS in (20) is convenient for quantifying the significance of a putative signal, and we will use it throughout the following analysis. In App. D we describe a datastacking procedure which reduces the data storage requirements for practical applications of our formalism.

B. Asimov Test Statistic
In order to build intuition for the information accessible to multiple detectors, we will use the Asimov data set [51] to study the asymptotic TS analytically. More precisely, the Asimov analogue of the TS in (20) is the average value taken over data realization, where the expectation value is taken on the data. In order to evaluate the asymptotic TS, it is convenient to separate the model prediction, which enters through Σ, into background and signal contributions: Referring back to (16), recall that the background is purely diagonal: Using this partitioning of the model prediction, we can then express the TS as follows Note that the values of B appearing in this expression are understood as being set to the value required by the profile likelihood technique. In order to evaluate the Asimov form of this expression, we only need to evaluate the average on the first term, as the average is taken over the data. This can be evaluated as follows, and then as the data has mean zero, we know the above expected value is simply given by the true covariance matrix, Here the truth parameters can be considered as, for example, the parameters one would use when generating Monte Carlo to simulate expected experimental results. For instance, to estimate the expected limit, the truth parameters would commonly have A = 0, whereas if we are estimating our sensitivity to features in f (v), we will take A = 0 in the Asimov data. In this work we are interested in the latter case, and therefore we will further assume the background has been fixed to the true value as a result of the profile likelihood technique, where S t is the true signal model and the same B appears in both the Asimov and model predictions. We further assume that the signal is always parametrically smaller than the background, which is the regime we will be in for any realistic experimental setup. 7 Implementing these assumptions, the Asimov form of (24) is 8 In (28) we have a convenient form of the expected TS that is amenable to analytic study. In the present work, our particular interest is the information contained in f (v) that we can only access as a result of DM interferometry. As such, it is convenient to evaluate a form of the Asimov TS, where all parameters except for those that (18), are set to their true values in the presence of a non-zero signal. If we further assume that our frequency resolution is sufficiently fine with respect to the scales over which the signal and background vary, then we can approximate the sum over Fourier components k with an integral over frequencies ω, or equivalently speeds v = 2ω/m a − 2. Under these assumptions, the TS becomes Much of the remainder of this work is devoted to studying the implications of this result.

C. Limiting Cases of Zero and Infinite Separation
We can use (29) to confirm basic asymptotic scalings expected for an analysis performed with DM interferometry. To begin with, the Asimov TS for a single detector with response A and background λ B , recalling F s ii (v) = 0 and F c ii (v) = f (v), is given by This expression agrees with the result in [16], which was derived for a single detector when analyzing the PSD. Importantly, we emphasize once more that (30) is only dependent on the speed distribution, so directional parameters which affect the velocity distribution but not the DM speed distribution are inaccessible. Before moving on to multiple detectors, we note that the expected discovery significance, which we denote TS, is given by the Asimov Θ evaluated at the model parameters that maximize the likelihood, which are the truth parameters.
In order to extract directional parameters we need at least two detectors. To that end, consider our expression in (29) for N = 2. For simplicity, we take A 1 = A 2 = A and λ B,1 = λ B,2 = λ B , in which case the Asimov TS becomes In particular, the discovery TS is given by Through F c 12 and F s 12 , the discovery TS depends on the spatial separation of the two experiments d = |x 12 |. In the limit where the experiments are close with respect to the DM coherence length, i.e. d λ c , then the two experiments see the same phase of the DM wave (k · x does not vary appreciably between them). In this case, we would expect a coherent enhancement in the signal. Defining for future use we see from (18) that for The N 2 = 4 enhancement of the TS represents a coherent enhancement, a point emphasized in [45]. This configuration provides a benchmark for the largest TS we can achieve for a general N = 2 configuration, and therefore will provide a convenient benchmark in the studies that follows. On the other hand, for widely separated detectors with d λ c , the DM fields will add incoherently. For the problem at hand, again returning to (18), we see that the sine and cosine factors will oscillate rapidly, driving the integrals to zero. What remains is, so that the TS now only scales as N , an incoherent enhancement.
The above argument can be readily generalized to N detectors. Typically the signal strength A is proportional to g 2 aγγ , so for N experiments all with pairwise separations d λ c , we expect our sensitivity to g aγγ should scale coherently as N 1/2 . If instead all experiments have d λ c , the scaling is reduced to N 1/4 , and for scenarios outside these two extremes the scaling will be somewhere in between. However, it is precisely this intermediate regime, where neither F reduces to the speed distribution nor vanishes, where we expect to be able to extract additional information about f (v). We turn to the problem of estimating parameters of f (v) in the context of the Asimov data set in the next section.

IV. ASIMOV PARAMETER ESTIMATION
In this section we will use (32) to perform frequentist parameter estimation and show explicitly that additional information about f (v) can be extracted via DM interferometry. For the purpose of simplifying the discussion, we will restrict our attention to the case of two detectors with equal background and detector responses as given in (32). However, the entire discussion can be readily generalized to N arbitrary detectors by using the asymptotic TS expression in (29).
To be specific, imagine we are interested estimating a set of signal parameters α, which are a subset of the full set of signal parameters α ⊂ θ sig related to f (v), and which have true values α t . The Asimov procedure allows us to study our ability to infer these parameters. For example, it is straightforward to confirm that Θ(α) is maximized for α = α t . 9 Beyond the best fit values, we are interested in determining the associated expected uncertainties and correlations between the various parameters, which are encompassed in the covariance matrix between the parameters, which we denote C. An estimator for this covariance matrix is given by the inverse Fisher information evaluated at the maximum likelihood, C −1 = I(α) (againα are the parameters that maximize the likelihood), where where we use (20). Given this relation, asymptotically our estimate for the covariance matrix is given by This expression involves the following shorthand for derivatives of functions then evaluated at their truth values, The expression in the first line of this result lays bare a simple fact: if Θ has no dependence on a particular parameter, for example the incident direction of a DM stream, or orientation of the Sun's motion through the DM halo, then the associated entries of the inverse covariance matrix vanish along with our ability to estimate that parameter. For the case of a single parameter α, we can readily invert the covariance matrix, and the above expression simplifies to where again all parameters are evaluated at their truth values after derivatives. We can already calibrate our basic expectation for parameter estimation from this result. Optimal estimation of α amounts to maximizing the right hand side of the expression; indeed, as expected, increasing the signal strength, A, or the integration time, T , both achieve this. If a parameter can be estimated from the speed distribution f (v) (in other words, ∂ α f (v) = 0), then that parameter may be estimated by a detector configuration with d λ c . However, the true power in the multi detector setup arises for parameters invisible to a single detector, defined by ∂ α f (v) = 0, but where ∂ α F c,s 12 (v) = 0. In generic cases, such parameters are optimally estimated for d ∼ λ c .
Continuing, let us assume that λ B is independent of frequency, in which case (38) becomes expressed in terms of TS 0 as introduced in (34). In particular, this result demonstrates the expected scaling of σ α ∼ (TS 0 ) −1/2 ; the exact details will require a specific f (v) and experimental configuration. In the following subsections we will continue this line of thinking, demonstrating in several toy examples that a second detector can lift degeneracies from the single detector likelihood.
, that carries the imprint of DM interferometry. Here we show the particularly simple example of an isotropic SHM for N = 2 detectors, in which case the expression is given in (40). The result is shown for various choices of the two detector separation d as compared to the axion coherence length λc = (mav0) −1 , with v0 = 220 km/s. The limiting cases of F c For d ∼ λc, however, the profile is modulated with the interference inherent in the cross-spectrum. In this simple case, there is no additional information about the velocity distribution that may be extracted by having multiple detectors.

A. The Minimal N = 2 Example
We begin our exploration of the above parameter estimation formalism with a simple scenario: N = 2 detectors measuring DM drawn from an isotropic velocity in the laboratory frame, 4πv 2 f (v) ≡ f (v). This example is obviously idealized; in reality, the finite boost velocity of the Sun about the Galactic Center implies that even an isotropic velocity distribution in the Galactic frame will become anisotropic in the laboratory frame. Nonetheless, this example will provide basic intuition for the impact of interferometry.
Invoking isotropy to perform the angular integrals, F c,s 12 (v) can be computed as where again d is the distance between the two detectors. Thus for this example, we see explicitly that for As we will see in the examples below, it is the dispersion v 0 rather than the average speedv which determines the crossover between small and large d.
To progress further, we assume a concrete form for f (v): the Maxwell-Boltzmann distribution, where v 0 is the velocity dispersion. Taking v 0 ≈ 220 km/s, this velocity distribution is an approximation to the SHM that is expected to describe the bulk of the local DM, neglecting the finite velocity boost of the Sun relative to the Galactic Center, which breaks the isotropy in the laboratory frame. We will utilize the Maxwellian ansatz repeatedly in this work as an illustrative example. In Fig. 2 which is a manifestation of the nontrivial correlations in the multi-detector spectrum. Note that we have defined, for the Maxwell-Boltzmann distribution, where v 0 is the velocity dispersion parameter that enters into (41); this is a particular realization of (1). Anticipating the more general scenario where the velocity distribution is not isotropic, it is precisely the deviation from f (v) that we will use to extract information about the full velocity profile.
As we have chosen an isotropic f (v), there is no additional information to extract about the velocity distribution in this case. Indeed, the distribution in (41) is defined by a single parameter, v 0 , which we can envision estimating. Evaluating (39) analytically in this case, we find written in terms of a dimensionless distance scale ξ = m a v 0 d = d/λ c , and Dawson's integral F . We find explicitly that σ v0 is minimized for ξ → 0, i.e. d λ c , since ∂ v0 f (v) = 0.

B. The Infinitely-Cold Stream
We now consider our first example of an anisotropic velocity distribution, a DM stream, and show that we can infer the direction of this stream using DM interferometry. In addition to the bulk SHM, it is expected that the local DM velocity distribution could contain non-virialized substructure, such as cold tidal streams [16,34,[52][53][54][55][56][57][58][59][60][61][62]. Streams are characterized by low velocity dispersions but large velocity boosts in the solar frame. Let us suppose that in the laboratory frame the stream is boosted at velocity v str and has velocity dispersion v 0 |v str |. In the limit v 0 → 0, the velocity distribution approaches a delta function, which has an infinite coherence length but a finite de Broglie wavelength. This is clearly an artificial example -it is the maximally anisotropic velocity distribution -but it is one we can evaluate fully analytically. Further, a number of the conclusions that we will reach for the infinitely-cold stream will hold also in more realistic cases. Note that for this example which has no dependence on the direction of the stream, and therefore a single detector cannot infer the direction. As claimed, for this simple scenario, we can compute the exact global TS using (32), and find Θ(θ str , φ str ) = TS 0 cos m a dv str (v str −v t str ) ·x 12 . (44) We consider the TS as a function of the spherical coordinates of our test stream direction, α = {θ str , φ str }, with the aim being to use the TS to infer the true direction of the stream, given by α t = {θ t str , φ t str }. In this case we can also compute TS 0 , as defined in (34), and we obtain 10 Now consider the angle-dependent factor in (44). Without loss of generality, we takex 12 =ẑ and define spherical coordinates with respect tox 12 , so that the argument of the cosine in (44) simplifies to where θ str and θ t str are the usual polar angles in spherical coordinates. Neither azimuthal coordinate φ appears in this expression, and hence the azimuthal angles are also absent in the TS. This implies we cannot infer one of the angular coordinates of v t str from the data. For our particular choice of coordinates, we can infer the parameter θ t str , as we will describe below, but the likelihood has a flat direction in φ str , so that we cannot infer the associated truth value. The degeneracy is physical. In our coordinates the symmetry of the likelihood is represented by an invariance under changes in φ, but more generally the TS is unchanged by rotations about the detector separation axis,x 12 . This can be seen from the dependence of the TS on (v str −v t str ) ·x 12 : any change in the test or truev str that is perpendicular tox 12 has no impact.
This symmetry of the TS under rotations aroundx 12 is in fact not a relic of our idealized example. Our ability to infer the direction of a velocity parameter vector that defines a given f (v) enters through the v · x 12 in F c,s 12 . But as v·x 12 is itself invariant to rotations of the velocity about thex 12 axis, one can show that this flat direction in the likelihood exists generally -indeed we will see it in more realistic cases (a direct analogue is apparent in the symmetry observed in Fig. 4, related to the SHM example discussed below). This symmetry will be broken by a dependence in the likelihood on multiple detector axes that are not parallel, provided either by a third detector or alternatively by daily modulation, where the singlê x 12 (t) will vary throughout the day at different times t.
We will explore this latter example in detail in Sec. Vindeed the optimal detector configuration will be determined by maximally violating this symmetry -but until then the symmetry will represent a basic feature of the physics.
Returning to our specific coordinate system wherê x 12 =ẑ, we may perform parameter estimation on the angle between the stream and detector. From (39), we have Note that the uncertainty on the parameter θ str is minimized for θ t str = π/2, i.e. when the stream is perpendicular to the two-detector axis. On the other hand, if the two vectors are parallel, θ t str = 0 or π, then we see σ θstr diverges. Yet we can still infer θ t str in this case. Indeed, looking to (46) we see that the asymptotic TS depends on θ str ; the likelihood is not globally flat, and we can estimate the angle from contours around the maximum likelihood. Instead, in this case there is a breakdown of the quadratic approximation around the maximum likelihood. If we were to incorporate higher derivatives than in (37), we would confirm that the likelihood is not truly flat at these points. This of course should be contrasted with the true flat direction in the likelihood associated with φ str . Note, however, that as θ t str approaches either 0 or π, becoming parallel to x 12 , the undetermined parameter φ str is less relevant. In the limit where the two vectors are parallel, we can infer the true direction of the stream, in spite of this degeneracy.
There is another interesting feature in (47): the result suggests that we can take d → ∞ to constrain this one direction of the stream to arbitrary precision. This is a manifestation of our assumption that the stream has no velocity dispersion: it remains coherent over arbitrary large distances, allowing for an improved baseline over which we can measure the stream direction. To study this feature further, imagine making this example slightly more realistic by introducing a finite velocity dispersion v 0 , with v 0 v str , such that f (v) has support in a small volume of radius ∼ v 0 around v str . For small enough v 0 we would expect the results of the δ-function stream to hold. Yet there is an important conceptual difference: the coherence length is no longer infinite because the different waves that constitute the local DM field now have speeds that vary by O(v 0 ). Parametrically, the argument of the interferometric terms scale as m a |v||x 12 | ∼ m a d(v str +O(v 0 )), but with the O(v 0 ) term varying between states. If we now take d (m a v 0 ) −1 , then the different waves will add incoherently, suppressing the power. But if we choose d ∼ (m a v 0 ) −1 , a degree of coherence can be maintained, along with the interference pattern carrying the information we seek to extract (see also the orange curve d = 2λ c in Fig. 2). Accordingly, for the optimal separation, the scaling of the sensitivities in (47) is (taking sin 2 θ t str ∼ 1/2 for definiteness) In the more realistic examples we will confirm the conclusion that d ∼ λ c provides the maximum sensitivity. There is another consequence of this choice. Taking d = (m a v 0 ) −1 in (44), the prefactor of the dot product in (46) is v str /v 0 1, by definition of this being a cold stream. Small variations in (v str −v t str ) ·x 12 will induce large variations in the argument of the cosine, implying that the global structure of the TS is highly nontrivial. Although the maximum TS will be attained at the true θ, there will be a pattern of local maxima with comparable TS (this result is depicted in Fig. 5, and persists even with daily modulation as shown in Fig. 7).

C. The (boosted) Standard Halo Model
The bulk DM halo of the Milky Way is expected to be Maxwell-Boltzmann distributed as in (41) in the Galactic frame, except for a possible cut-off around the escape velocity ∼500 km/s [63]. On the other hand, the Sun is boosted with respect to the Galactic frame by [64] v ≈ (11, 232, 7) km/s , in Galactic coordinates, wherex points towards the Galactic Center,ŷ points in the direction of the local rotation of the disk, andẑ points towards the Galactic north pole. Thus in the laboratory frame (neglecting the Earth's motion), the velocity distribution becomes that of the SHM, with a velocity dispersion v 0 ≈ 220 km/s [65,66]. Note in particular that v 0 ∼ |v | ≡ v , so for the SHM λ c ∼ λ dB . The associated speed distribution is As we have emphasized many times already, single detectors are only sensitive to the speed distribution, which only depends on v but not the orientation of the solar velocityv . Thus, a single detector may constrain the model parameters v 0 and v (as shown in [16]), but determining the orientation requires multiple detectors. 11 To determine the expected sensitivity to the direction v we need to compute the derivatives of F c,s 12 (v) that appear in (39): cos θ cos θ t cos(m a vd cos θ) where I 0,1 are both modified Bessel functions (an analogous expression holds for F s 12 ). In computing this result we have again chosen coordinatesx 12 =ẑ, but left the direction ofv arbitrary, defined by (θ t , φ t ). The most important feature of this result is that it exhibits no dependence upon φ t : again, there is a symmetry in the likelihood for rotations aroundx 12 . Beyond this, we can also see that the derivative vanishes when v is parallel to the detector separation (θ t = 0). Accordingly, in this case we will find σ θ diverges, as we did for the stream. But again this is not a global flat direction in this case; the likelihood is just sufficiently flat at the maximum that the first three derivatives vanish.
To proceed beyond these analytic insights, we will compute the remaining results numerically. We define the angle between v and x 12 as θ . To begin with, we take a generic value of θ t = π/4 and consider how well we can infer this angle as a function of detector separation. The results are shown in Fig. 3. Unlike for the δ-function stream, there is now a minimum at a finite value of d, and as argued on general grounds this occurs when d ∼ λ c = (m a v 0 ) −1 . That the uncertainty diverges for d → 0 is consistent with the fact that a single detector cannot infer this direction. In more detail we find the minimum occurs at d ∼ 2λ c , where we obtain σ θ d ≈ 2/ √ TS 0 . For example, if TS 0 = 25, corresponding to a 5σ local significance detection with d = 0, then at the distance d ∼ 2λ c corresponding to minimum uncertainty, the solar velocity direction with respect to the detector axis could be localized to 0.4 rad ∼ 20 • on the sky. We can understand the magnitude of σ θ d at its minimum from (48): the SHM has the form of a stream where v 0 ∼ v str = v , and therefore we would expect σ θ d × √ TS 0 ∼ 2, exactly as observed. However, as we have emphasized already, it is important to keep in mind that our estimate of σ θ is a measure of the expected curvature of the likelihood in the vicinity of the true value and does not capture the global structure of the expected likelihood function. To illustrate these features we fix d = 2λ c for definiteness and illustrate the global map Θ(θ , φ )/TS 0 in Fig. 4, for three different values of θ t . Note that we have divided out the overall significance TS 0 , so that exactly how well we can localize the direction will depend on how significantly the DM signal has been measured. However, the expected global structure of the TS will be a rescaled version of these maps. In each case the truev that we are 11 In principle annual modulation may be used by a single detector to inferv , as discussed in [16,34]. In this example we have set the true orientation to θ t = π/4. With this configuration, we find that the maximum precision is obtained for d ≈ 2λc.
seeking to infer is located in the center of the Mollweide projection maps. The left panel illustrates the scenario withx 12 =v (θ = 0), the center has θ = π/4, while in the right panel two directions are perpendicular and x 12 points between the poles of the map (θ = π/2). In all cases the symmetry of the TS around thex 12 axis is apparent. The only case where this flat direction in the maximum TS is not an obstruction to determining the true direction ofv is when θ d = 0. In that case we are still able to localize the true direction, although we note the likelihood is relatively flat around the maximum (consistent with the second derivative vanishing). In Sec. V we illustrate how daily modulation generically allow us to fully determine both of the angles associated with the direction of v .

D. The Sagittarius Stream
As a final example working with a single staticx 12 , we return to the case of the cold stream with non-vanishing velocity dispersion. We expect many of the conclusions reached in Sec. IV C to hold in this case. In particular, the symmetry around thex 12 axis will remain, but we will see explicitly in this case the non-trivial structure induced in the global likelihood by the ratio of v 0 /v str ∼ λ dB /λ c 1. To make the example concrete, the DM component of the Sagittarius stream may ex-   Θ(θ, φ) for the SHM divided by the co-located detection significance TS0. The detectors are configured so that the displacement vector between them is parallel to the SHM boost velocity, and the Mollweide plot is rotated so that it is centered on the maximum test statistic. (Center ) As on the left, but for a detector configuration where the displacement vector is at a 45 • angle to the North (θ t = π/4) with respect to the SHM boost velocity. (Right) As on the left, but for a detector configuration where the displacement vector is perpendicular to the SHM boost velocity (θ t = π/2). In this configuration the location of the boost velocity can only be localized to a great circle on the celestial sphere.  Fig. 3, but for the Sagittarius (SGR) stream rather than for the SHM. As before, the maximum precision for the inferred value of θstr is achieved at mav0d ≈ 2, although the overall dependence is somewhat softened outside of the extremes at mav0d = 0 and mav0d = 2π. The values of σ θ str × √ TS0 are also considerably smaller than those found in the SHM example, indicating that the angle θstr can be reconstructed with much greater precision for the SGR stream as compared to the SHM. (Right) The Asimov TS Θ(θstr) for the SGR stream rescaled by the co-located detection significance TS0 as a function of θstr for a detector configuration where the true stream direction is θ t str = π/4 (dashed vertical line). We have fixed mav0d = 2. The TS Θ(θstr) is maximized at the true value of θstr, but there is considerable nontrivial global structure with a large number of local minima and maxima in Θ.
tend to the Sun's location, and estimates [67,68] suggest that it could make up ∼5% of the local DM density. However, the DM associated with the stream would be highly collimated in phase space; we follow [16] and model the Sagittarius stream DM velocity distribution by a boosted Maxwellian as in (50), but with v 0 = 10 km/s and v replaced by v str = (0, 93.2, −388) km/s [53,55]. We consider the stream in isolation, as opposed to in conjunction with the bulk SHM DM phase-space distribution, because even though the stream component is sub-dominant in terms of DM density, it still dominates in the narrow region of phase space where the Sagittarius stream has compact support. To simplify the discussion we will simply take v str = 400 km/s, with a direction that we will again specify by its angle with respect to the detector axis (given the degeneracy in rotations about that axis). Note that this example could apply equally well to other putative DM streams, such as the newly discovered S1 stream [60][61][62].
To begin with, in the left panel of Fig. 5 we show the expected uncertainty on the recovered angle between the stream and detector, θ str , as a function of the distance in units of λ c = (m a v 0 ) −1 , for a true value θ t str = π/4. This figure is the stream analogue of what we showed for the SHM in Fig. 3. Once more, following the general discussion in Sec. IV B, the optimal sensitivity is achieved for d ∼ λ c , and from (48), we expect σ θ d × √ TS 0 ∼ 2v 0 /v str ∼ 0.05 at the minimum-uncertainty distance, compatible with what we see in Fig. 5.
However, just like in the case of the SHM it is important to also examine the global properties of the TS in addition to the curvature of the expected TS at the true parameter values. Towards that end, on the right of Fig. 5 we show the expected TS Θ, normalized to TS 0 , as a function of the reconstructed angle between the stream and detector, θ str . For this figure we have fixed the true orientation at θ t str = π/4 along with the separation d = 2λ c . We see that Θ drops off quickly around the true value of θ str = π/4 (vertical dashed), but that there is non-trivial structure with local maxima at larger and smaller θ str values. This is a direct manifestation of the non-trivial interference patterns discussed in Sec. IV B for cold streams: the large ratio v str /v 0 enters into the argument of the trigonometric functions in F c,s 12 (v).

V. DAILY MODULATION
One of the most dramatic signatures of DM interferometry is the unique daily modulation signal available to multiple detectors. This effect, which we describe in the current section, would be a smoking gun signature that an emerging excess has a DM origin, and it also allows two detectors to better determine geometric parameters describing the velocity distribution. The basic idea is simply that for two detectors fixed at generic locations on the surface of the Earth, the separation vector x 12 is rotating in the inertial Galactic frame throughout the day. This is in contrast to the angular parameters entering in the DM velocity distribution, such as the Solar directionv , which should always point in the same Galactic direction, regardless of the orientations of the detectors at any point in time on Earth. The rotation of x 12 with respect to the fixedv implies that we will sample a variety of angles between the two vectors, and therefore vary the modulation of the speed distributions in F c,s 12 (v), as already depicted in Fig. 1. Critically this will lift the flat direction in the maximum likelihood associated with rotations around x 12 that we observed repeatedly in Sec. IV: as the likelihood will now depend on a collection of different vectors x 12 (t) , the symmetry that exists around any one of them will not be preserved in the full TS.
In the rest of this section we divide the discussion into three parts. Firstly, we describe how to construct the likelihood for the generic case of N detectors incorporating daily modulation and describe how it is straightforward to generalize our full formalism to this case. We then focus on the specific case of N = 2 and show, within the Asimov formalism, how the examples of the SHM and Sagittarius stream discussed above are modified in the presence of daily modulation. Finally, we turn to a Monte Carlo simulation of a realistic example and demonstrate how, within a day, a resonant experiment could constrain the direction of the solar velocity vector, v , that controls the SHM to sub-degree accuracy.

A. A Likelihood with Daily Modulation
So far in this work, we have envisioned a set of N experiments collecting measurements of the signal-plusbackground frequency spectra for a duration of time T while the detector separations were fixed with respect to the boost velocity of the DM component under consideration. However, this framework cannot be extended to the case of daily modulation, as the signal prediction will fundamentally be varying over a 24-hour period. In order to properly account for this effect, the data must be collected in time intervals of duration T 24 hours and analyzed with a joint likelihood over all the collected intervals. In detail, if we imagine that we collect M such time intervals, indexed by r = 0, 1, . . . , M − 1, then for each of these we will have a data set d r = {d k,r }, where again k labels the Fourier mode. For each data set d r , we can compute the likelihood as in (19), and the full joint likelihood is the product of these over r. Explicitly, we have Importantly, note that we have also attached an index r to the signal prediction Σ(θ), as we need to account for the variation of the detector separations x ij throughout the day. In a similar fashion the full formalism of Secs. III and IV can be extended to include the varied detector orientation: within a given sub-interval we simply adjust x ij as appropriate, and then we form joint quantities by combining these as in the likelihood above. To provide just a single illustrative example, the Fisher information computed in (36), would become with other expressions similarly generalized.

B. Asimov Examples with Daily Modulation
While the alteration to our formalism imposed by daily modulation is minimal -as exhibited in (53) -the impact on the results can be dramatic. We will demonstrate this with several examples in this section, all within the Asimov formalism. To begin with, we consider using N = 2 detectors in order to determine the direction of  Fig. 4, we construct Mollweide projections of the Asimov test statistic Θ(θ, φ) for the SHM rescaled by the co-located detection significance TS0. However, we now perform a joint likelihood over data collected over a 24-hour period so that the daily modulation of the detector displacement vector produces a time-varying signal, which helps break degeneracies in the reconstructed directional parameters. The Mollweide projection for a configuration in which the detectors are oriented along an East-West (North-South) orientation is shown on the left (right). While the results of obtained in an East-West configuration do not depend on the latitude of the detectors, the North-South configuration results do, so for definiteness, we have taken the detectors to be located in New Haven, CT, the site of the HAYSTAC detector. In both configurations, the SHM boost velocity direction can be localized effectively, although there remains a non-trivial degeneracy in the East-West map between two points on the sphere. v in the SHM. This is the same problem we considered in Sec. IV C, which produced the results shown in Fig. 4, where there is a clear degeneracy associated with rotations around x 12 . We will now see explicitly that daily modulation helps lift this degeneracy. To do so, let us suppose that the DM velocity distribution follows the SHM in (50), with v 0 = 220 km/s and v = 232 km/s. Our goal, as previously, will be to infer the direction ofv . We consider two detectors separated by d = 2λ c = 2/(m a v 0 ), and for definiteness we place one detector at a latitude 41.3 • N and longitude 72.9 • W. In Fig. 6 we show results where a second detector is placed a distance d to the East (left) or North (right) of this detector with data stacked at two-hour intervals over the 24-hour period. 12 For the North-South configuration, we see that the direction can be well-localized: a high significance axion detection in this case would lead to a precise estimation of the direction of v , as we show explicitly in Sec. V C below. This configuration clearly outperforms an East-West configuration, where there remains a degeneracy that has not been fully lifted by the daily modulation. Additionally, the maximum test statistic realized in the North-South configuration would be approximately 10% larger than one realized in an East-West configuration for otherwise identical data collections. 12 Note that since the Earth's rotation is aligned with the East-West direction, results obtained for the East-West configuration are independent of the exact experimental locations, so long as the detector separation is much smaller than the Earth's radius of curvature. For any other configuration, however, the result will generically depend on latitude.
Using the same experimental design, we can also revisit the example of the Sagittarius stream discussed in Sec. IV B. In Fig. 7 we construct the analogue of Fig. 6, but now for the much colder stream. Note that since v 0 for the stream is a factor of ∼ 20 smaller than for the SHM, the optimal detector distance d = 2λ c is a factor of 20 larger than in Fig. 6. Although in both configurations the TS is maximized at the expected location on the sphere, nontrivial structure due to the presence of many local maxima are apparent in both the North-South and East-West configurations. We note that, as in the SHM example, there is only one global maximum for the North-South configuration, located at the true direction of the stream. However, there remains a degeneracy in the East-West configuration.
The degeneracy represented in the Mollweide maps for the SHM and the Sagittarius stream in the East-West configuration is exact. It has its origin in the dimensionality of the space swept out by the detector separation vector x 12 over the course of the day. As studied in Sec. IV, for data taken at fixed x 12 , the test statistic Θ(v) evaluated as a function of the orientation of the boost velocity depends only on the angle betweenv and x 12 . As a result, Θ(v t ) = Θ(v ) wherev t is the true boost direction andv is a velocity obtained by reflectingv t across any plane which contains x 12 . For detectors in an East-West configuration, the Earth's rotation produces a daily modulation of x 12 that is confined to the plane orthogonal to the Earth's rotational velocity vector. As a result, the TS measured at each point in the day, and therefore the sum of such TSs, will be exactly preserved under reflections of the boost velocity across that plane. This means that accounting for daily modulation in the  Fig. 6, but for the Sagittarius stream example. For a fixed axion mass, the physical detector separation d = 2λc is a factor of 20 larger than in Fig. 6 because of the larger coherence length of the stream. While there are many local maxima in both configurations, the North-South orientation produces only a single global maximum, at the true detector localization, while the East-West orientation leads to two degenerate global maxima (one at the true detector location and the other displaced). An animated version of these figures, showing how the localization improves throughout the day as more orientations of x12 are sampled, can be found at github.com/joshwfoster/DM Interferometry.
East-West configuration the directional parameters can only be determined up to a reflection across the plane perpendicular to the Earth's rotation axis. By contrast, for detectors in the North-South configuration, the set of detector separation vectors throughout the day will generically not be co-planar, and thus there is no analogous degeneracy. 13

C. Monte Carlo Example with Daily Modulation
As a realistic demonstration of our ability to perform parameter estimation using the daily modulation effect, we generate a Monte Carlo realization of data in the North-South SHM scenario (as depicted in the right panel of Fig. 6) using A = 38.25 and λ B = 1, both of which we take to be dimensionless without loss of generality. The values of A and m a was chosen to generate a signal of expected 5σ significance during a 100-second collection in a single detector to mimic a realistic resonant scanning strategy in which one of two independently operated detectors detects an excess and both are then used for a 24-hour observation of the excess candidate. We constructed 24 hours of Monte Carlo data for this signal, taking a detector separation of d = 2λ c ; with these parameters, the excess would be expected to appear at TS 0 ≈ 60, 000 after 24 hours. While large, this TS is consistent with the power of a resonant strategy once the axion mass is known.
Using uniform priors on A between [33,43], on v between [212.5, 252.5] km/s, on v 0 between [200, 240] km/s, on λ B between [.999, 1.001] and a uniform prior on the sphere for (θ , φ ), we construct a Bayesian posterior distribution for the model parameters. The results of an analysis performed using Multinest [69][70][71][72] with 2,000 live points are shown in Fig. 8. In particular, we see that the true location of the stream has been located to degree precision. This precision can be understood from (48), which gives the expectation σ θ ∼ 0.5 • , consistent with what is shown in the figure. Let us suppose that the Sagittarius stream, as modeled in this work, comprises 10% of the local DM. In the example above, we would expect that after 24 hours the location of the stream could be localized to ∼ 10 ; interestingly, this represents greater accuracy for stream localization than localization of the bulk SHM even though the stream is a sub-dominant component of the DM.

VI. CONCLUSION
In this work we have demonstrated the power of DM interferometry for wave-like DM. The spatial coherence of the DM field imprints phase correlations on the signals observed at spatially-separated detectors, and these phase correlations are sensitive to parameters in the full 3-dimensional velocity distribution f (v), whereas a single detector is blind to all effects beyond the speed distribution f (v). As a result, the advantages of DM interferometry go beyond a simple coherent enhancement of the signal strength as the number of detectors is increased. By taking advantage of the fact that the correlation matrix of the Fourier-transformed signals at multiple detectors depends on modified speed distributions which contain modulated forms of f (v), we have demonstrated that parameters such as the solar velocity vector may be reliably extracted from two detectors separated by a distance d ∼ λ c . Furthermore, directional parameters of coherent substructure such as DM streams may be estimated at even higher significance, though in that case the optimal separation λ c is parametrically different from the DM de Broglie wavelength λ dB .
Our formalism has immediate practical applications for new and upcoming axion DM experiments. The sensitivity to g aγγ for resonant-cavity axion experiments which use external magnetic fields, like ADMX and HAYSTAC, is typically BV 1/2 , where B is the peak magnetic field strength and V is the magnetic field volume. In order to achieve resonant enhancement, the volume of an individual cavity is fixed to be of order 1/m 3 a , so to achieve greater sensitivity, one must either increase the B-field or construct a multiplexed readout with multiple cavities. Assuming the latter strategy is chosen, our results motivate placing at least one of the cavities at a distance λ c : if a signal is detected, the loss of coherent enhancement of the signal is more than compensated by the ability to localize the boost direction of the DM velocity distribution to within 1 degree with just 24 hours of data.
While there are many challenges to the construction of additional instruments, we emphasize that all of the important phenomenology is captured by a two-detector array. This smoking-gun signature of DM is invisible to a multiplexed setup where all cavities lie inside a single coherence length. A similar analysis applies to experiments in the quasistatic regime like ABRACADABRA and DM-Radio, where the physical volume of the experiment is decoupled from m a . For both types of experiments, our formalism may be applied to determine the optimal detector orientation for localizing the solar velocity to the desired precision (with North-South orientations generally being preferred to East-West). The optimal detector separation corresponds to physically reasonable distances for well-motivated axion masses -O(10) m for the 10 −5 eV mass range of ADMX and HAYSTAC and O(1000) km for the 10 −9 eV targeted by ABRACADABRA/DM-Radio -and as such the coherence length and detector orientation can form an important design parameter for future experiments, in much the same way as L/E determines the design of neutrino oscillation experiments.
The future of axion detection involves readout beyond the standard quantum limit, using tools such as Josephson parametric amplifiers and squeezed states. In this regime, it is important to note that our variables R k and I k are canonically conjugate, and thus cannot be simultaneously measured to arbitrary precision. In future work, we plan to investigate how our formalism must be modified for quantum-limited readouts. As the number of new axion experiments proliferates, this work motivates careful consideration of the spatial configuration of multiplexed detectors. In this appendix we briefly review the concepts of the coherence length and time, as relevant to wavelike DM. We emphasize that in our work both concepts only arise heuristically. Indeed, the coherence length and time are only defined parametrically, and for all quantitative results we instead rely on the likelihood formalism in [16], which produces not only all parametric scalings, but also the required O(1) factors.
Consider first the coherence length λ c ∼ (m DM v 0 ) −1 , the scale over which wavelike DM remains coherent. In discussions of ultralight DM, "coherence length" is often used interchangeably with "de Broglie wavelength." Strictly speaking, though, the de Broglie wavelength λ dB = 2π/(m a v) is a property of particles with fixed velocity v, while the coherence length describes the dephasing of various plane wave components with different velocities. When v 0 ∼ v, these two length scales are com-parable, but there are relevant situations where the two diverge, such as for cold streams, and then the distinction between the coherence and de Broglie wavelengths becomes important.
The coherence time is then the timescale over which a measured signal of ultralight DM will build up coherently. In real space, this is the time it take for a new spatially coherent packet of the DM wave, which has size λ c , to arrive at the instrument. If these packets travel with a mean speed ofv, then the timescale is τ ∼ λ c /v ∼ (m DMv v 0 ) −1 . The same result can be arrived at from a frequency space consideration. The Fourier transform of an experimental data set collected over a time T will have a frequency resolution of ∆ω = 2π/T . If the entire signal fits within a single frequency bin, the result is associated with a single draw from an exponential distribution, as shown in [16]. Once we resolve the signal, however, we obtain multiple draws which will combine incoherently, partially offsetting the benefit of additional integration time. The coherence time is therefore dictated by the width of the signal in frequency space, and then as dω = m DM v dv, we again arrive at τ ∼ (m DMv v 0 ) −1 .
The goal of this appendix is to demonstrate a fact that was used without proof in the main body: the data set d, given in (15), is a random variable drawn from a multivariate normal distribution with zero mean and covariance matrix Σ as given by (16). In order to show this we will start from the known statistics of the axion field, as reviewed in the main body, together with a Gaussian background, and show that the mean and variances of the data sets follow the expected normal distribution. We will further confirm this result with a Monte Carlo realization of the axion field. From here, rather than confirm that all higher moments are also consistent with Gaussianity, we will instead confirm numerically that the distribution is normal. Indeed, the diagonal components of d, which govern the statistics of individual detectors, must be Gaussian as proven in [16].
Let us begin by restating (6) in a simplified notation. We introduce a single multi-index d = abc, and a random variable We now envision collecting a data sensitive to this axion field at each of the N detectors, located at positions x i . Specifically, we imagine collecting N measurements at a frequency f = 1/∆t at each experiment, so that we have at our disposal N × N data points {Φ The second term in this expression captures the background noise. We will assume the noise is Gaussian, which holds for a wide range of sources as described in the main body, and in detail that it satisfies In other words, we assume the noise has zero mean, is uncorrelated between detectors, and has a variance that increases with the measurement frequency f . The variance is controlled by the mean power in the background, λ B,i , and if there are multiple background sources at a single detector, their power can simply be combined. From this data set, we compute the discrete Fourier transform {Φ (i) k } using (7), and then the associated real and imaginary parts, R  (8). These variables are what combine to form the data vector d, and so the goal is to study their statistics. Before proceeding, let's introduce some further notation to keep expressions compact. Firstly, we encapsulate the axion phase into a single term, To capture the trigonometric sums introduced by the Fourier transforms, we write and the equivalent expression for sine is denoted s n,k . Using this, the real and imaginary parts of the data set can be written From these expressions, we can see immediately that R That this holds for the background follows from (B3), and for the axion signal contribution we have The first step follows as the value of α d (and hence f d ) is uncorrelated with φ d (and hence ϕ (i) d,n ), whilst the second utilizes the fact cos ϕ = 0 when the argument ϕ is a random phase. This establishes that d = 0.
Next we consider the covariances. In particular, we will compute R k . The calculation where one or both of the real components is replaced by an imaginary equivalent proceeds similarly, and we will comment on the important differences throughout. In detail, we will compute Note the effect of sending R (j) k is simply to replace c m,k → −s m,k , and similarly for R (i) k . Continuing with the calculation at hand, expanding out the final two lines, we will have expressions involving only the signal, only the background, and cross terms. As the background value is uncorrelated with the signal, the cross terms will be zero via an almost identical argument to the vanishing of the means. Of the remaining terms, consider the background first.
(∆t) 2 T  Figure 9. Monte Carlo validation that the statistics of DM interferometry are as claimed in App. B. In the left figure we confirm that the variances of the real and imaginary signal-only data sets, collected for the N = 2 experiments, is as claimed in (16). This was proven directly in the text, but in the plot we show that the average of 4,000 Monte Carlo simulations provides a consistent prediction for the variances as a function of frequency in the different cases. On the right figure, for the frequency where R (1) R (1) achieves its maximum, we show the distribution of values across the simulations. In detail, we see that the real and imaginary components are normally distributed, and consistent with a mean-zero normal distribution, where the variance is given as on the left, here σ 2 ≈ 25 Wb 2 /Hz. We found that the distributions were consistent with the Gaussian expectation at the level of p > .05 using the D'Agostino and Pearson omnibus normality test [73,74]. In both cases, each Monte Carlo simulation involves a direct construction of the axion field starting from (B15) with Na = 100, 000, taking ma = 2π Hz, and A = 1 Wb 2 . Further, we take the velocity distribution to follow a variant of the SHM in (50), but with v0 = 0.07 and v = (0, 0.08, 0), both in natural units. The (unphysically) large velocity helps simplify the computation of the Fourier transform. The detector separation is x12 = d(0, 1, 0), with d ≈ 4.4λc.
which we can simplify further as, In the final step we see the emergence of the k · x type phase factors that separate F c,s ij (v) defined in (18) from f (v). We have broken the calculation into a number of pieces at this stage, let us begin to put things back together. Combining the different expressions above, we have The final result is the claimed form of R used in the main body, but let us detail the steps in the calculation, working backwards. In the last step we simply recalled the definitions introduced in (18) and (17). The penultimate step simply involved approximating the sum over all velocity components d = abc with an equivalent integral. The only non-trivial manipulation occurred when we evaluated the sums over n and m. These were performed using a set of discrete Fourier transform double orthogonality relations, which for convenience we have collected in App. C. From those relations, we can see that as R involved c n,k c m,k , only the cosine of k · x ij survived. By analogy, if we were evaluating I (where the background contribution vanishes as described above), we have c n,k s m,k , which instead singles out the sine, implying the above result would have c ij (ω) → s ij (ω). The same argument holds for I (i) k R (j) k , up to a sign. Taken together, the above arguments suffice to demonstrate analytically that the variance of the data set is as claimed in the main body. We can also confirm this result numerically. On the left of Fig. 9, we show that a direct construction of the axion field as a sum over N a plane wave components, where v i is drawn from f (v) and φ i is drawn uniformly from [0, 2π), leads to the exact same results. 14 The detailed parameter choices are described in the figure caption, and the curves represent the average over repeating 14 Binning the velocities leads to (6) in the main text, with a Rayleigh-distributed amplitude in each bin.
this procedure 4,000 times. In all cases, there is excellent agreement between this approach and the corresponding theory curves.
On the right of Fig. 9 we confirm a point that we did not demonstrate directly, namely that the individual real and imaginary components are normally distributed. The distribution is shown amongst the 4,000 simulated data sets for the two components measured at two different detectors. In all cases consistency is observed with the predicted Gaussian distribution. We performed a chisquared test to determine the goodness of fit and found p-values greater than 0.05. In detail, the R (1) , I (1) , R (2) , and I (2) data sets shown in Fig. 9, had corresponding p-values of 0.06, 0.12, 0.97, and 0.27.

Appendix C: Orthogonality Relations
In App. B we made use of several unstated orthogonality relations. We collect these in the present appendix. Firstly, the following expressions vanish for any k However, recall that we usually exclude these exceptional k values from our likelihood. The non-zero results above were written in terms of Dirac δ-functions, however this is an approximation. Recall all results are obtained through the discrete Fourier transform, within which the frequency can be interpreted as ω = (2π/T )k, with k = 0, 1, . . . , N − 1. In truth, if we define k d = ω d T /2π , then what appears in the above sums is the Kronecker-delta δ k k d . However, in the spirit of assuming our frequency resolution is sufficient enough to approximate ω as a continuous variable, we take which is the form it appears in (C2) and (C3).

Appendix D: Data Stacking Procedure
In practical situations it is usually neither feasible nor necessary to save the entire time-series data to disk and then construct the Fourier transform of the full data set. The frequency resolution of this complete Fourier transform would be ∆ω = 2π/T , and potentially much smaller than the scale of any expected features induced by the signal due to f (v). As a specific example, the ABRACADABRA-10 cm experiment [24,25] recorded the PSD data over short time periods and then stacked the PSD data over the time subintervals to construct the average PSD data. The advantage of this averaging procedure is that it requires less storage and is easier to deal with computationally, since there are less frequencies involved than would be in the full data set without time sub-binning.
With this in mind, it is useful to understand how we may stack the Fourier transform data over multiple experiments in such a way that we preserve the full power of the likelihood in (19) but that allows us to reduce the data volume needed to be saved to disk. (An optimized procedure for stacking the data from a single experiment is presented in [16].) Let us imagine that we record time-series data in N T equal time subintervals of time ∆T = T /N T , and that in each subinterval the frequency spacing of the ∆N = N/N T Fourier components is sufficient to resolve the axion signal by multiple frequency bins, i.e. we retain sufficient frequency resolution that our signal remains well resolved. We then denote the full data set by d = {d k }, indexed now by both k = 1, . . . , ∆N − 1, denoting the Fourier component, and = 1, . . . , N T , the data subinterval. The appropriate likelihood is then simply the product of the likelihood in (19), but now also over all values of N T . However as N T × ∆N = N , the number of frequency bins in the Fourier transform of the full data, at this stage we have not reduced the size or complexity of the data or likelihood evaluations at all. In order to do so, we can combine the data into the following average data matrix, which can be computed prior to any evaluation likelihood, Here, the indices i and j run over the 2N entries of the data vector in (15), k indexes the discrete Fourier transform, and specifies the appropriate subintervals. In terms of the average data matrix, the likelihood can be written as where we have left the θ dependence of Σ implicit. We can now compare how much data needs to be stored for this stacking procedure compared to the full data set. Again, we have N T subintervals, each with ∆N Fourier components, and for each we have 2N components in our data vector. As (D1) is a real symmetric matrix, we need N (2N +1) components to specify it for each k value. Thus in total, we need to store N ×(2N +1)×∆N entries to disk, although if Σ −1 k has a number of zeros (associated with experiments well within or outside the coherence length λ c ), fewer points may be required. This number should be contrasted with the 2N × N = 2N × ∆N × N T values that would be needed in the absence of a data stacking procedure. Thus, as long as N T N , a significant reduction in the data set can be achieved. For simplicity, in the main body of the paper we assume that no data stacking has been performed, though it is important to keep in mind that all results we derive may also be applied to the stacked data likelihood. An important caveat is that care should be taken when accounting for daily modulation to make sure the data is stacked with other data taken at a similar time of day, otherwise the effect can be washed out.
Finally, we briefly demonstrate using the Asimov procedure that as long as the subintervals retain sufficient frequency resolution that the signal remains well resolved, the stacked and full likelihoods are equally sensitive. If the signal prediction remains unchanged in each subinterval, then the averaged data set defined in (D1) has the following expected value, It is straightforward to then evaluate the equivalent Asimov Θ, and one finds a result enhanced by N T , but with T → ∆T when replacing the sum over Fourier components with an integral over speed. For instance, the equivalent of (29) has T → N T ∆T . Yet as N T ∆T = T , by definition, the test statistic is identical, and therefore the stacking procedure is optimal as claimed.