Axion astronomy with microwave cavity experiments
Abstract
Terrestrial searches for the conversion of dark matter axions or axionlike particles into photons inside magnetic fields are sensitive to the phase space structure of the local Milky Way halo. We simulate signals in a hypothetical future experiment based on the Axion Dark Matter eXperiment (ADMX) that could be performed once the axion has been detected and a frequency range containing the axion mass has been identified. We develop a statistical analysis to extract astrophysical parameters, such as the halo velocity dispersion and laboratory velocity, from such data and find that with only a few days integration time a level of precision can be reached matching that of astronomical observations. For longer experiments lasting up to a year in duration we find that exploiting the modulation of the power spectrum in time allows accurate measurements of the Solar peculiar velocity with an accuracy that would improve upon astronomical observations. We also simulate signals based on results from Nbody simulations and find that finer substructure in the form of tidal streams would show up prominently in future data, even if only a subdominant contribution to the local dark matter distribution. In these cases it would be possible to reconstruct all the properties of a dark matter stream using the time and frequency dependence of the signal. Finally we consider the detection prospects for a network of streams from tidally disrupted axion miniclusters. These features appear much more prominently in the resolved spectrum than suggested by calculations based on a scan over a range of resonant frequencies, making the detection of axion minicluster streams more viable than previously thought. These results confirm that haloscope experiments in a postdiscovery era are able to perform “axion astronomy”.
pacs:
14.80.Va; 95.35.+dI Introduction
Axions are light pseudoscalar particles that appear in the solution of Peccei and Quinn Peccei and Quinn (1977); Kim and Carosi (2010) to explain the unnatural absence of CPviolation in quantum chromodynamics (QCD). More modern motivation from the landscape of axionlike particles (ALPs) appearing in string theory Svrcek and Witten (2006); Arvanitaki et al. (2010); Cicoli et al. (2012), inspires the generalisation of axions to light pseudoscalars with a theoretical origin outside of the original PecceiQuinn solution. Such ALPs can cover an extremely wide range of masses and couplings to the standard model Arias et al. (2012). Axions and ALPs are an attractive cold dark matter candidate and can be produced in the early Universe through a variety of nonthermal mechanisms such as vacuum misalignment or the decay of topological defects Davis (1986); Hiramatsu et al. (2012) in ways that are consistent with the known cosmological abundance and phenomenology of dark matter Ipser and Sikivie (1983); Wantz and Shellard (2010). For a recent review of axion cosmology see e.g., Ref. Marsh (2016). In the following we use the term ‘axion’ to refer to both the QCD axion and generic axionlike particles.
Certain axion mass and coupling ranges can be ruled out with various astrophysical observations such as the cooling of white dwarfs Raffelt (1986); Isern et al. (2008), neutron star interactions Berenji et al. (2016), the lifetimes of horizontal branch stars in globular clusters Raffelt (1999) and supernovae neutrinos Burrows et al. (1989); Mayle et al. (1988). Axions may also be observable in the lab. Experimental tests for axions predominantly rely on their coupling to electromagnetism resulting in the Primakoff conversion of axions into photons inside strong magnetic fields. Cavity resonators can exploit this effect if the resonant frequency is chosen to match the axion mass Sikivie (1983). As the axion mass is unknown experiments must be designed to scan over a range of resonant frequencies corresponding to a range of axion masses. Experiments include the helioscope CAST Zioutas et al. (2005) (and the planned IAXO Armengaud et al. (2014)) searching for axions produced inside the Sun, and haloscopes such as ADMX searching for Galactic dark matter axions Asztalos et al. (2010). These experiments operate over a narrow range of frequencies and hence make constraints on the axion mass in small bands, where the smallest accessible photon coupling is controlled by the signaltonoise level of the experiment. Planned dark matter axion experiments such as QUAX Barbieri et al. (2016); Crescini et al. (2017); Ruoso et al. (2016), CULTASK Chung (2016) and the layered dielectric haloscope MADMAX Caldwel et al. (2016); Millar et al. (2016) are being designed to probe ranges of axion masses inaccessible due to the technical restrictions of the ADMX design. As well as axion ‘observation’ experiments there exist solely labbased experiments such as the Any Light Particle Search BÃ¤hre et al. (2013) using the technique known as lightshiningthroughawall Van Bibber et al. (1987). In haloscopes beyond ADMX, it may be possible to search for lighter axions with broadband readout circuits Kahn et al. (2016) or LC circuits Sikivie et al. (2014) as well as heavier masses in the meV range with Josephson junctions Beck (2015, 2013). Other couplings such as those to nuclei Arvanitaki and Geraci (2014) can be probed in nuclear magnetic resonance experiments such as CASPEr Budker et al. (2014); Graham and Rajendran (2013) and the coupling to electrons can be constrained using conventional dark matter direct detection searches for weakly interacting massive particles (WIMPs) Aprile et al. (2014). For a recent review of axion experiments see for example Ref. Graham et al. (2015).
The goal of a haloscope experiment is to tune the frequency of a cavity mode to the axion mass resulting in the resonant enhancement of the axionphoton conversion. ADMX achieves this with the use of movable tuning rods placed inside the cavity itself to modulate the resonant frequency over a range of several hundred MHz. Although usually unimportant when scanning over a relatively large range of resonant frequencies, the velocity distribution of axions in the halo would cause a small frequency spread in the resonance Krauss et al. (1985). Furthermore there are also 0.5% and 1% modulations in time due to the relative velocity of the Earth and Sun with respect to the halo dark matter ‘wind’ Turner (1990); Ling et al. (2004); Vergados and Semertzidis (2017). There is also a potentially exploitable modulation dependent on cavity orientation with respect to the incoming axion direction for axion masses (and haloscope volumes) experiencing a loss of coherence over the cavity dimensions Irastorza and Garcia (2012).
Here we consider a scenario in which the axion mass has been determined down to a level of precision dictated by the scanning approach of ADMX and a dedicated high spectral resolution experiment is then performed at a single resonant frequency. In such a situation, the shape of the spectrum of axionphoton conversion will be accessible. This power spectrum is related to the velocity distribution of the local dark matter halo, hence the precise functional form and parameters which arise from astrophysics are important. Past axion searches with ADMX have incorporated some of these astrophysical uncertainties, for example by searching for discrete flows of axions Duffy et al. (2006); Hoskins et al. (2016, 2011) or applying constraints to different halo models Sloan et al. (2016); Vergados and Semertzidis (2017). In the situation we consider here however it will be possible to perform “axion astronomy” in the sense that a measurement can be made directly of the axion power spectrum to learn about the structure of dark matter in the Galaxy. For this reason we develop an analysis that shares similarity with wellestablished methods for extracting astrophysical information in the case of WIMP direct detection experiments e.g., Ref. Strigari and Trotta (2009). Since axion haloscopes effectively observe the axions directly  as opposed to WIMP direct detection experiments which observe the WIMP flux convolved with a stochastic scattering process  the prospects for sensitive measurements of the dark matter halo are much greater. Here, we show how powerful future ADMXlike experiments might be for doing axion astronomy. We discuss this idea in the context of simple analytic halo models, distributions from Nbody simulations, and minicluster streams  a phenomenon unique to axion dark matter Hogan and Rees (1988).
To begin in Sec. II we review some of the basic theory for axions and ALPs, as well as the laboratory frame speed distribution relevant for calculating signals inside a haloscope experiment. In Sec. III we outline the steps in calculating the expected power inside a magnetic cavity resonating at a given frequency. Then in Sec. IV we describe our mock experiment and explain the procedure used to extract astrophysical information from the simulated data. The first results applying these methods to the reconstruction basic sets of input parameters are presented in Sec. V and then extended to Nbody data from the Via Lactea II (VL2) Diemand et al. (2007) simulation in Sec. VI. Finally we extend this simulation to tidal streams from disrupted axion miniclusters in Sec. VII, before summarising in Sec. VIII.
Ii Axions and ALPs
First we outline some of the essential steps in calculating the resonantly enhanced axionphoton conversion power inside a microwave cavity. Full details of these calculations can be found in Refs. Hong et al. (2014); McAllister et al. (2016); Krauss et al. (1985). Importantly we wish to make the connection to realistic halo velocity distributions so we depart from an often used approximation that the axion power spectrum can be described with a BreitWigner function.
The axion to two photon coupling is given by the formula,
(1) 
which includes the fine structure constant and the axion decay constant . The dimensionless coupling is,
(2) 
In which is the ratio of the colour axion anomaly to the electromagnetic axion anomaly and is the ratio of the up and down quark masses. The value of this constant is model dependent: for the KSVZ model Kim (1979); Shifman et al. (1980) and for the DFSZ model Dine et al. (1981); Zhitnitsky (1980) for example. In the interest of model independence and to generalise to ALPs we express the interaction in terms of the coupling .
The effective Lagrangian for axions coupled to electromagnetism is,
(3)  
where is the electromagnetic field strength tensor, and its dual. The axion potential is provided by QCD instanton effects and can be approximated with a simple mass term . The axionic charge density and the electromagnetic current density are written as and . Writing we then see the axionphoton interaction in terms of electric and magnetic field strengths,
(4) 
This interaction modifies Maxwell’s equations to include an additional axion current,
(5)  
(6)  
(7)  
(8) 
However these equations simplify for the setup we consider here. Firstly we assume the axion field has no spatial dependence on laboratory scales (). We can do this because the size of ADMX is around the 1 meter scale and is well below the de Broglie wavelength of the axion field for the mass ranges we consider (100 m). This allows us to assume that there is no spatial dependence in the axion field over the dimensions of the cavity and hence no additional modulations due to the changing orientation of the cavity with respect to the axion wind. We also assume that there is no axionic charge and no electromagnetic current inside the cavity: and . This results in the following simple set of equations,
(9)  
(10)  
(11)  
(12) 
Under the above assumptions the equation of motion for the axion field is,
(13) 
Dark matter axions in the Milky Way undergo essentially no interactions, so in a quadratic potential , the field oscillates coherently at the axion mass . However due to thermalisation in the Milky Way the coherence of the oscillations is spoiled slightly by a dispersion in axion velocities: . One can account for this by moving to a Fourier description of the field, written as ,
(14)  
(15) 
where is some large reference time used to take the averages. The quantity is referred to as the axion power spectrum. The rms of the axion field squared is connected to the axion power spectrum by the Parseval relation,
(16) 
The convention in dark matter detection literature is to use a velocity distribution to describe the kinematics of the local halo. In this context we must blur the distinction between the interpretation of axionic dark matter as a classically oscillating field and as a collection of particles. The velocity distribution is related to the axion power spectrum in the following way. First we write down the distribution of axion velocities in the laboratory frame (i.e., by temporarily introducing a number density,
(17) 
where is the number density of “particles” with speeds between and . The constant is found from integrating over all velocities and is used to define the local axion number density . This allows the connection to a classical field oscillating at , which should have , to be made Krauss et al. (1985).
An expression for the axion power spectrum can be obtained by satisfying Parseval’s relation and changing variables from to ,
(18) 
we can then substitute for using,
(19)  
(20) 
where this expression clarifies the distinction between a 3dimensional velocity distribution and its 1dimensional speed distribution . Hence, the formula for the axion power spectrum on Earth can be written as,
(21) 
The simplest assumption for a dark matter halo is the Standard Halo Model (SHM) which is spherically symmetric with a density profile and yields a MaxwellBoltzmann velocity distribution,
(22) 
To simplify the following expressions we use the peak velocity for the shape of the distribution, however it can also be written in terms of a velocity dispersion, . The speed distribution follows from an integral over all directions,
(23) 
However the power spectrum should reflect the fact that we observe not so we need to compute the speed distribution in the laboratory frame. To do this we make the transformation which yields for the velocity and speed distributions,
(24) 
and,
(25) 
Since we are now in the moving laboratory frame a time dependence appears in the speed distribution. Finally after changing variables to we arrive at the axion power spectrum. The axion power spectrum must be 0 for which is enforced by requiring that be real.
For use in later examples we also define the velocity distribution for streams which are spatially and kinematically localised substructure components of the dark matter halo. Their velocity distribution can also be described with a Maxwellian,
(26) 
where the stream is parameterised by its Galactic frame velocity km s and dispersion km s. We will assume that the stream comprises some fraction of the local dark matter density.
The description of the lab velocity is well known in the context of WIMP direct detection but is not usually considered for axion detection. Since we are reliant on its precise description we will briefly elaborate on its contents. The lab velocity is written,
(27) 
containing respectively, the bulk rotation velocity of the stellar disk (the local standard of rest, LSR), the peculiar velocity of the Sun with respect to the LSR, the orbital speed of the Earth around the Sun and the rotation speed of the Earth about its axis. The latter two velocities are responsible for the annual and diurnal modulations respectively and are known theoretically with effectively perfect precision (see the Appendix of Ref. Mayet et al. (2016) for a review of these calculations). In the SHM the velocity of the LSR is usually written as in Galactic coordinates, with km s. An often quoted value for the peculiar velocity of the Sun from Schoenrich et al. Schoenrich et al. (2010) is km s with roughly 1 km s sized systematic errors.
Iii Resonance power
We model a microwave cavity experiment with a static uniform magnetic field maintained inside a cylindrical cavity of radius and length , with radial, azimuthal and vertical coordinates labelled respectively. The magnetic field is generated by a solenoid with current density in the direction. The electric and magnetic fields we write as,
(28)  
(29) 
where is the Heaviside step function, is the current and is the number of wire turns in the solenoid per unit length. For convenience we use the magnitude of the magnetic field in the following expressions.
In the cylindrical cavity design the important cavity mode orientations are the TM modes which have transverse magnetic fields (in the direction and hence have associated electric fields in the direction). It is useful to write these induced fields in terms of their Fourier components,
In this case, Ampère’s law from Maxwell’s equations reduces to,
(30) 
Solving this equation inside and outside the cavity and matching boundary conditions leads one to a solution for the Fourier components of the axion generated magnetic and electric fields. The solutions are resonances at particular frequencies corresponding to the zeros of a Bessel function (although we will only be interested in the lowest resonance which we label ). Following the derivation of Ref. Hong et al. (2014), the axion power is calculated by evaluating the following integral over the volume of the cavity ,
(31) 
where is the energy stored in the electric and magnetic fields inside the cavity. This expression introduces the quality factor which is a number that quantifies how well the cavity stores energy and depends on the material properties of the cavity wall. Evaluating the above formula with the solution for the Fourier components of the axion electric and magnetic fields (which are expressed in terms of ) one arrives at,
(32) 
where is the th zero of the 0th Bessel function of the first kind. We have also defined , which is a Lorentzian that describes the loss in power off resonance,
(33) 
Usually the haloscope power is written in terms of a cavity form factor. For the transverse magnetic field considered here (TM) this is written ^{1}^{1}1Other mode orientations, the transverse electric (TE) and transverse electromagnetic (TEM) modes both have no axial electric field meaning they have negligible form factors.. We are principally interested in the TM mode which has . ADMX can tune the TM mode from roughly 500 MHz to 900 MHz Asztalos et al. (2010). In general the electric field of the TM mode can be written Jackson (1998),
(34) 
In which, is the time dependent component of the field, is a Bessel function, is the th root of , is the cavity radius and is the cavity length. Modes with and have very small form factors.
Our simulation is based upon the calculation of Eq. (32) so for our purposes it would be sufficient to stop here. But in the interest of comparison with previous calculations we will calculate the power on resonance. To do this we simply set and use a BreitWigner approximation for the axion power spectrum with an analogous factor (this permits an analytic evaluation of the integral in Eq. (32)). We also introduce the axion density by writing . Resulting ultimately in,
(35) 
where we have restored the factors of , and for completeness. If the quality factor of the resonant cavity is very high (i.e., the cavity is very good at storing energy and the dissipation is very slow) then the axion conversion power is limited by the spread in axion kinetic energy. The factor arises from the integral of two BreitWigner functions and indicates how the total power received on resonance is set by the wider of the two power spectra.
Inputting typical values for the experimental parameters we arrive at a total power of the order W as is usually quoted,
(36)  
Iv Mock experiment
Our simulation is an approximation of the current ADMX setup. We list a set of benchmark experimental parameters in Table 1. The magnetic field strength, quality factor and noise temperature are roughly in line with what is currently achievable. For calculating the time dependence we also include the latitude and longitude of the experiment.
Axion: 
3.4671 eV  
10 GeV  
Experiment:  8 T  
70,000  
220 l  
0.2 s  
10 days  
2 years  
4 K  
Latitude  
Longitude  
Maxwellian halo:  0.3 GeV cm  
km s  
km s  
Stream:  km s  
20 km s  
0.05 
In this section we will consider a hypothetical scenario in which the axion has been discovered after a successful low resolution scan over a wider mass range. Once the resonance has been found then an experiment can be performed at a single frequency. The running time of the experiment needs to be long enough to ensure that the signaltonoise ratio is high but for our purposes also needs to be comprised of long timestream samples to obtain high frequency resolution in the resulting spectrum.
For now we pick a benchmark set of particle parameters that lie in the QCD axion band: MHz (= 3.4671 eV) and GeV. This choice evades existing constraints but is easily within the reach of ADMX given a long enough running time at the correct frequency. We use only a single particle benchmark in this study as we are placing the focus on the underlying astrophysical parameters. This is justified however because many of the conclusions are either independent of the choice in mass and coupling (provided the running time and resonant frequency are suitably adjusted) or have dependencies that are simple to explain from the scaling of the axion power. We discuss how one might extend our conclusions to other axion mass and coupling ranges in the Summary Sec. VIII.
The sensitivity of a haloscope experiment is limited by the strength of the axion conversion power compared to the noise level. There are two main sources of background noise in resonant cavity experiments: the signal amplifier and the cavity walls. The cavity walls produce thermal blackbody photons or Johnson noise whereas the amplifiers produce electrical noise which depends on the precise technology, however both can be modelled as white noise Daw (1998); Hotz (2013); Brubaker et al. (2016). The signaltonoise ratio for a haloscope experiment of duration , is set by the Dicke radiometer equation Dicke et al. (1946)
(37) 
where is the bandwidth of the axion signal and is the noise temperature.
Our mock experiment consists of a long total running time which is divided into separate time integrated bins of length . Inside a given time bin we calculate a power spectrum which would correspond to the average of Fourier transformed timestream samples of duration . The Fourier transform of a given sample is a power spectrum with frequency resolution . The noise we simulate as Johnson white noise which has rms power inside a given frequency bin with an exponential distribution Duffy et al. (2006). The noise power spectrum of the average of individual exponential power spectra corresponding to the Fourier transformed timestream samples then approaches Gaussian white noise in accordance with the central limit theorem. Hence our simulated noise inside the larger time bin is Gaussian white noise with mean value and standard deviation . The full dataset then consists of a total number of time integrated power spectra each of which consists of the axion power spectrum averaged over the time added to the Gaussian white noise.
The major motivation for a long running time, aside from simply reducing noise, is to utilise the annual modulation due to which provides a Galactic perspective to the signal. It has previously been shown that the annual modulation signal allows astrophysical parameters to be measured more accurately using WIMP direct detection data, as it breaks degeneracies Savage et al. (2006). Below we show that this is also the case for axion searches. We test our simulation by first generating a mock dataset and then attempting to reconstruct the input particle and astrophysical parameters with a maximum likelihood analysis. Two examples of such data are displayed in Fig. 1 corresponding to two halo models, a smooth isotropic Maxwellian distribution and a pure stream (with parameter values listed in Table 1). The annual modulation of the signal is indicated by the purple line labelled .
We base our likelihood on a statistic which measures the offset between the observed value of power , and the expected power (signal + rms noise) in each bin, where and label frequency and time bins respectively,
(38) 
where the sums run over frequency bins and time bins. The error is given by the suppressed rms noise power . We then construct a likelihood based on this statistic. Mathematically the likelihood as a function of a set of parameters given data is,
(39) 
where we assume , and are free parameters. We also use the generic to label a set of astrophysical parameters as we will perform tests with varying numbers of free parameters. We use to incorporate the optional uncertainty in the knowledge of an astrophysical parameter (it can be set to unity if no prior knowledge is assumed about a given parameter). The final term parametrises the likelihood of the noise power which can be measured externally (although we set this to unity unless otherwise stated).
Our likelihood analysis consists of first finding the parameter values that maximise the likelihood of Eq. 39, we interpret this set of parameters as the best fit points. We then construct 68% and 95% confidence regions around these points which are either 1dimensional intervals when we are only interested in the reconstruction of one parameter or 2dimensional contours when we are interested in the reconstruction of two parameters and their correlation. We do this by first profiling over all other parameters other than those of interest and then calculate a likelihood ratio between the maximum likelihood for each point in the remaining likelihood function. The likelihood ratio we define as, , where are the maximum likelihood estimators. According to Wilks’ theorem Wilks (1938) this is asymptotically distributed for parameters. We then find intervals or contours around the best fit points which enclose regions of parameter values with less than a certain critical value. The critical value of is that for which the cumulative distribution of is the desired confidence level. For example for the 68% interval encloses values of and the 95% encloses values of .
V Reconstructing parameters
In this section we use the simulation and analysis methodology described in Sec. IV to attempt to reconstruct sets of input particle and astrophysics parameters. The aim is to quantify how accurately and with what correlations and degeneracies a future ADMXlike haloscope experiment would measure the local axionic dark matter distribution. In the following results we show 1 and 2dimensional 68% and 95% confidence intervals/contours calculated using the profile likelihood, along with best fit parameters values which maximise the likelihood. To explore the likelihood function we use nested sampling algorithms provided by the MultiNest package Feroz and Hobson (2008); Feroz et al. (2009, 2013), and set a tolerance of and use between and live points depending on the number of parameters being reconstructed.
In Fig. 2 we show the reconstructed axion parameters and (left) and the astrophysical parameters and (right). We show three sets of contours which correspond to experiments of different durations: 10 days, half a year and 1 year. The 10 day long experiment corresponds to a single time integrated bin of the 0.5 and 1 year long experiments. The annual modulation signal does not play a large role in constraining these parameters, hence the effect of increasing the experiment duration is to shrink the confidence intervals by a factor . The axion mass and coupling can be measured to a high level of precision even with only 10 days of data taking, however there is some bias in the best fit values since the dataset consists of a single realisation of stochastic noise. The shapes of the contours are roughly one sided for masses due to the fact that the axion power spectrum is only nonzero for . The astrophysical parameters can be measured to a high level of accuracy too. With a 1 year duration the level of precision would reach around the 1 km s level which roughly matches the accuracy of current astronomical observations Schoenrich et al. (2010).
With a full annual modulation signal we can also access the 3dimensional components of . However since and are summed in the Galactic frame we can only measure directly the and components of . The component (i.e., that which lies along the direction of the rotation of the Milky Way) can only be measured in combination with the LSR speed . In Fig. 3 we show the measurement of these parameters for the same three experiment durations of 10 days, 0.5 and 1 year. Since the 10 day duration experiment consists only of a single time integrated bin we have no annual modulation signal and only the reconstruction of the largest component () is possible as this has the greatest influence on the shape of the spectrum. The remaining two components have essentially flat likelihoods as the single time bin spectrum is not sensitive to their values. However for longer durations with modulation in time, the measurement of all three components becomes possible. Even with only half a year of the annual modulation signal we can still make a measurement of the three components of . However, as the signaltonoise is lower the measurement is biased by particular large fluctuations, which in this example leads to the input values lying outside of the 95% contour. With a full year of data however a very accurate measurement can be made with 95% confidence intervals smaller than 5 km s and the true values (indicated by dashed lines and stars) lying within the 95% interval in all cases.
Finally in Fig. 4 we show the 1 and 2 sigma error bars for various parameter measurements as a function of the total experiment duration . We use three experiment durations from 1 day to 1 year and for each we repeat the experiment 30 times with different randomly generated noise in each to demonstrate the sensitivity to the individual data realisation. As shown in Figs. 2 and 3 the short duration experiments as well as setting much weaker measurements are also biased by the particular data causing some reconstructions to lie further than 2 sigma away from their input values. In the case of the axion mass we expect onesided measurements due to the onesided nature of the power spectrum. This is the case for 10 day and 1 year durations, however for the 1 day duration we see multiple experiments reconstruct a mass smaller than the input mass due to large noise fluctuations in bins slightly below the axion mass. Interestingly for the longer duration experiments the constraint on the axion mass reaches a level smaller than a single frequency bin (5 Hz), this is because the shape of the power spectrum and the annual modulation signal also provide additional information about . The size of the error bars for the remaining parameters decrease roughly as and for durations long enough to exploit the annual modulation signal we see a significant decrease in the scatter in the reconstructed values over different realisations of the experiment. This means that a future experiment such as this would be able to make fine measurements of the axion particle parameters in conjunction with astrophysical parameters and with no significant biases.
Vi Nbody data
We can source more realistic examples of dark matter distributions from Nbody simulations of Milky Waylike halos. These might more accurately reflect the inhomogeneities and anisotropies that will likely be present in a real dark matter halo. This is of particular interest for a high resolution axion experiment because, as shown in the previous section, it is far more sensitive to astrophysical parameters than standard axion searches and WIMP direct detection.
We use data from the Via Lactea II (VL2) Diemand et al. (2007) simulation and select 200 analogue Earth locations at a Galactic radius of 8 kpc and calculate a velocity distribution from all particles contained within 1 kpc spheres centred on each of these locations (we also enforce that no spheres overlap). Although there are more recent hydrodynamic simulations which will better reflect a Milky Waylike dark matter distribution, the VL2 data is sufficient for the illustrative examples we show here and will not change the general conclusions.
We display the range of these 200 velocity distributions in Fig. 5 with certain samples labelled which contain a significant substructure component. We label these samples from #1  #4. The substructure appearing prominently here are types of tidal stream which are present in real galaxies due to the orbits of satellite galaxies. As these smaller galaxies orbit close to their host halo both the stellar and dark matter components can be tidally stripped leaving a long trail of material around the galaxy which may intersect the main galactic disk. In our own Milky Way there has been evidence from observations and simulations that a particular tidal stream from the Sagittarius dwarf galaxy could pass very close to the location of the Solar System Purcell et al. (2012). Tidal streams are of particular interest here as they are very kinematically localised. In particular, streams incoming with velocities at an angle with respect to the motion of the Solar System would give rise to unique annual modulation signals.
We calculate the axion conversion power spectrum in the same way as before but we substitute the analytic with a discretised version calculated by binning particle velocities with a bin size roughly corresponding to the frequency resolution of the experiment. Importantly for each time bin at we rotate all particle velocities into the laboratory frame with the time dependent Galactic to laboratory transformation detailed in Appendix A. We must also boost all particle velocities by .
In Fig. 6 we show a selection of four axion conversion power spectra for a range of sample VL2 velocity distributions (the same selection as labelled in Fig. 5). The four examples are selected because they contain significant substructure components in the form of streams. These show up in the power spectra as sinusoidally modulating features in time, some examples such as #2 and #3 having single dominating streams whereas others such as #4 possess multiple streams with different amplitudes and phases.
We can parameterise the frequency dependence of the modulating streams with the function,
(40) 
In principle the three parameters of this function are related to the three Galactic frame components of the stream velocity, although this will not be a onetoone mapping. The frequency of the stream modulation is proportional to the quantity .
We can extract substructure components from the data we have presented here by searching for sinusoidally modulating features that have a period of 1 year (whilst also fitting for the function Eq. (40)). First we can reduce the data by subtracting the time averaged spectrum and then dividing by the standard deviation of the remaining fluctuations. Next we perform a cut of bins with power fluctuations below a certain level of significance leaving a series of points which if the stream component is large enough will retain the sinusoid modulation. The resulting data points for each example are shown in the left hand panel of Fig. 6. These data points can then be fit to a model for the stream. We again use the same Maxwellian form for the stream velocity distribution as in Eq. (26) with power spectrum shown in the lower panel of Fig. 1. Whilst the stream is unlikely to be perfectly described by a Maxwellian, any deviations will be smaller than the error induced by the finite frequency resolution and noise fluctuations.
A given stream is described by its density , dispersion , and three components of velocity making a total of five parameters. Since we have a method for extracting the stream from the data, we can use the data that remains once the stream is removed to fit for the axion, halo and lab velocity parameters and break the degeneracy with the stream parameters. In Fig. 7 we show the reconstructed stream velocities for the four example halo samples displayed in Fig. 6. Note that in all cases all components of the stream velocity can be reconstructed to high accuracy due to the prominence of the annual modulation signal. This is because the three components of the velocity can all be independently measured with the use of the phase, mean frequency and amplitude of the modulation in Eq.(40), although this relationship is nonlinear due to the transformation from the Galactic to the laboratory frame.
Also in Fig. 7 we show the measurement of stream density fraction and dispersion for each sample. Because the density fraction and dispersion are respectively related to the power amplitude and width of the modulating feature, a reconstruction of these parameters is possible in addition to the velocity components. The four samples we have considered here all have relatively large stream contributions which aids the measurement of these parameters. For weaker streams it is likely that longer duration experiments would be required to increase the signaltonoise. Here the lowest density stream that is detectable with our method is set by the power with respect to the level of noise. Furthermore we have not explored the full stream velocity parameter space with these four examples. It is likely that the accuracy of the reconstruction will be dependent on the direction of the stream with respect to the direction of the lab velocity. Additionally with higher signaltonoise examples it should also be possible to reconstruct more than one stream (as in sample #4). We leave these issues however to future work.
Vii Axion miniclusters
There has been sustained interest in small high density bound structures of axions called miniclusters (see e.g., Refs. Hogan and Rees (1988); Kolb and Tkachev (1993, 1994a, 1994b, 1996); Berezinsky et al. (2013); Tinyakov et al. (2016)). Miniclusters are formed in the early Universe from density perturbations in the axion field. Perturbations which can form miniclusters result from various types of nonlinear dynamics involved with axion oscillations such as vacuum misalignment or the decay of axion defects such as strings and domain walls Chang et al. (1999). Previous work has predicted the existence of up to Tinyakov et al. (2016) locally if all of the dark matter was in the form of miniclusters, though a direct encounter would occur less than once every years Kolb and Tkachev (1994b). Through close interactions with stars however axion miniclusters would become tidally disrupted leading to a network of streams wrapping the Milky Way (possibly in addition to a smooth component of the dark matter halo). The miniclusters will pass through the stellar disk many times over the age of the Milky Way ( Gyr). It has been estimated in Ref. Tinyakov et al. (2016) that a small population of disrupted miniclusters would lead to several streams along the path of the Earth through the Galaxy that are large enough to induce an enhancement in the observed total power. The final result of Ref. Tinyakov et al. (2016) is a value for the number of expected stream crossings with a density larger than the local smooth halo density , which is interpreted as an amplification factor. However if the axion minicluster streams are an additional component to the smooth component then the stream density does not need to be larger than the local density to provide an enhancement to the signal. Since the velocity dispersion of the minicluster streams is extremely small compared to the halo ( km s km s), in a high resolution axion experiment all of the minicluster stream axions would convert to photons in a small number of frequency bins. Hence for a minicluster stream to be observable we simply need the total power from the stream to be larger than the power over a few bins.
Individual miniclusters are parameterised by the density contrast in the axion field, which is a number typically of order unity. The distribution of values of found from the simulations of Ref. Kolb and Tkachev (1996) appears to follow a function similar to which we will use as an approximation. The mass of a minicluster is set by the total mass of axions inside the Hubble radius around the time when axion oscillations begin, (which is allowed by lensing bounds Zurek et al. (2007)). Ref. Kolb and Tkachev (1996) states that the distribution of minicluster masses is concentrated tightly around a large fraction of this mass.
Solving for the collapse of a spherical region with density contrast and evolving through cosmic time to the present day gives a range of minicluster densities,
(41) 
We assume that the miniclusters are spherically symmetric with central density and radius . We assume a Maxwellian distribution for the speeds of axions inside a minicluster with a dispersion set by the virial velocity . From Ref. Zurek et al. (2007) we assume the miniclusters have a power law density profile with for but enforce the gradient to turn towards 0 at to give a central density of .
The number of streams expected at the Solar radius results from evolving the initial distribution of axion miniclusters through the age of the Galaxy to today. Each time the minicluster crosses the stellar disk there is a probability that it will encounter a star close enough to become disrupted. Following previous calculations of this type Goerdt et al. (2007); Schneider et al. (2010), Ref. Tinyakov et al. (2016) gives the probability of a particular minicluster being disrupted,
(42) 
where here is the orbital speed of the minicluster, and the number of crossings of the stellar disk the minicluster undergoes. This calculation has already averaged over an isotropic distribution of minicluster trajectories and has been written in terms of the stellar contribution to column density in the direction perpendicular to the disk, Kuijken and Gilmore (1989). Given this, we can just use miniclusters with circular orbits intersecting the Solar position () to evaluate the number of crossings over the age of the Galaxy () to be roughly . Note that the dependence on drops out of Eq. (42). This is because although faster miniclusters cross the stellar disk more frequently (), they are also less likely to encounter a star during a given crossing (). We also note that has no dependence on since and are both proportional to .
A stream can be specified alone by four parameters: the density contrast and mass of the original minicluster, the age of the stream , and the orbital velocity of the minicluster/stream, v. All other parameters can be derived (we indicate dependence on each by parentheses). Once a minicluster is disrupted by a star it will begin to leave a trail of axions along its orbit, the length of which will stretch linearly with time as the cluster orbits the Galaxy . Assuming the stream retains the original radius of the minicluster and is simply deformed from a sphere of radius into a cylinder of length , the density of the axions for a minicluster stream of age is,
(43) 
Reference Tinyakov et al. (2016) calculates the number of expected stream crossings in a 20 year period for two values for the age of the Galaxy and two masses. We extrapolate the final result of this work down to stream densities of as this is around the lowest density stream that would be observable in this case. We estimate that if this extrapolation of Ref. Tinyakov et al. (2016) is valid then, for Gyr and , there could be up to stream crossings in a 20 year period (although the precise number is not important for the illustrative example we present here).
We simulate the signal for minicluster stream crossings by selecting samples from the parameter space . First we select values for from the distribution . We then select v from an isotropic MaxwellBoltzmann distribution. Finally we draw a value of such that the number of stream crossings with follows the function presented in Fig. 2 of Ref. Tinyakov et al. (2016). The length of time taken to cross the stream is then approximately,
(44) 
which is derived from the distance travelled through the stream, , where is the angle between the stream and the path of the Earth. We distribute each of these crossings uniformly over the running time of the experiment. The power spectrum observed during a crossing is enhanced with an additional Maxwellian component (as with the streams the previous section) with relative velocity and dispersion . Also in a given time bin the minicluster stream signal will gain an additional spread in frequency from the change in over the duration of the bin.
To deal with Eq. (44) diverging for stream directions that align with the path of Earth we remove all streams which orbit with relative to the plane of the stellar disk, where is the width of the stellar disk. This is a safe approximation as this is only a small fraction of the streams, and miniclusters that orbit in the plane of the stellar disk will become disrupted much earlier than those orbiting at a large angle and the streams will hence have much lower present day densities.
In Fig. 8 we display a simulated power spectrum observed over a total period of 10 years for a halo consisting of a smooth population of axions and a network of tidal streams from miniclusters. The streams appear as peaks in the power spectrum over a very narrow range of frequencies (as in Sec. VI) but here since minicluster radii are on the scale of km they are shortlived enhancements compared with usual tidal streams which extend over volumes larger than the scale probed by the Galactic orbit of the Solar system.
The total power measured in the form of these shortlived enhancements would provide an estimate of the fraction of local axion dark matter contained in minicluster streams from which the abundance of miniclusters could be inferred. We emphasise however that a detailed theoretical treatment of the disruption of a population of miniclusters is still needed in order to fully explore the prospects for their detection. Our example here shows that even if miniclusters comprise only a very small contribution to the local axion density, they appear much more prominently in a high resolution experiment. In principle one could make use of the methods described in Secs. V and VI to extract information about individual streams such as their radius, age and Galactic frame velocity as well as place constraints on the minicluster population such as their mass spectrum and abundance. This would require isolating the signal from miniclusters from both the noise and the background axion power spectrum. A possible strategy could be to use the observations during periods without any minicluster enhancement to make accurate measurements of the underlying parameters (as in Sec. V) to then subtract the background spectrum thus isolating the stream crossing events.
A further complication that we have not discussed in this work is the presence of any shortlived environmental peaks which may appear in real resonant cavity experiments and could mimic a positive axion signal. These would usually be dealt with by performing a repeat experiment in the frequency range at which the peak was observed. However in the case of minicluster streams which are themselves shortlived this check would not necessarily be successful if the timescales for the environmental peak and the stream crossing were comparable. However a careful treatment of the frequency modulation of the peak over time may in some cases be enough to distinguish a Galactic signal from a labbased one. We leave a more detailed study of axion minicluster streams and implications for experiments to a future work.
Viii Summary
We have performed a simulation of a hypothetical high resolution ADMXlike experiment following a successful detection of an axion dark matter signal. Our focus here has been on extracting astrophysical information and performing axion astronomy. Our main conclusions are as follows:

The measurement of the axionphoton conversion power spectrum enables the accurate reconstruction of both axion particle parameters in conjunction with the underlying astrophysical parameters.

With the use of the annual modulation signal one can make accurate measurements of the components of the Solar peculiar velocity. With an experimental duration longer than a year the accuracy can reach below 1 km s, which would improve upon the measurement from local astronomical observations.

Substructure such as tidal streams appearing in simulations of Milky Waylike halos show up prominently in the resolved axion power spectrum and can hence be measured to levels of sensitivity not possible in the direct detection of WIMPs. The annual modulation signal plays an important role here too as the precise shape of the modulating stream allows the reconstruction of its properties: the Galactic frame velocity, density and dispersion. This in principle would allow axion haloscopes to trace the formation and accretion history of the Milky Way.

We have simulated an approximation to the expected signal from a population of streams from disrupted axion miniclusters. We have extrapolated a result for the calculation of the expected number of stream crossings from Ref. Tinyakov et al. (2016). In an experiment that resolves the axion spectrum the signal from minicluster streams would appear much more prominently in the data and could be isolated to place constraints on their mass spectrum or abundance.
The issues we have discussed here are relatively unstudied in the context of axion detection. Hence there are a number of areas in which this study might be extended. We have shown that measuring the axion power spectrum allows accurate reconstruction of underlying parameters and although we have only considered simple models here, in principle the same should be true of other models for the dark matter velocity distribution such as those containing anisotropy parameters or additional substructure such as debris flows Kuhlen et al. (2012), dark disks Schaller et al. (2016) and caustic rings Duffy and Sikivie (2008). What remains to be seen however is the extent to which the correct selection of a particular model is possible with data of this kind. This is an important consideration for WIMP direct detection experiments with very low statistics, multiple competing experiments and degeneracy between assumptions about the underlying velocity distribution. These issues have given rise to a number of approaches for making astrophysics independent limits and measurements Frandsen et al. (2012); Fox et al. (2011); Gondolo and Gelmini (2012); Del Nobile et al. (2013); Fox et al. (2014); Feldstein and Kahlhoefer (2014); Kahlhoefer and Wild (2016); Gelmini et al. (2016) and developing general parameterisations for the velocity distribution Peter (2011); Kavanagh and Green (2013); Kavanagh (2014); Kavanagh and O’Hare (2016). In the case of axions however, because the power spectrum could be measured to an arbitrary level of precision given sufficient duration it may not be necessary to develop any such astrophysics independent methods, however this would require a separate investigation.
We have used only one axion benchmark mass and coupling, since our focus is on measuring the underlying astrophysical parameters. However our conclusions can be simply extended to other values by considering Eq.(36). Since the total axion power is proportional to one can extend any of the reconstructions to smaller couplings by scaling up the experiment duration, , by the same factor squared. Although it should be noted that haloscopes can reach smaller couplings by both reducing noise as well as simply extending the duration of the experiment and both of these approaches are necessary for improving constraints on the axion coupling. Since the total power is proportional to , our conclusions still hold for smaller (larger) values of the axion mass if is increased (decreased) by the same factor. The reverse argument goes for values of the local density since the power is linearly proportional to . However we must take care in extending these results to axion masses much larger or smaller than eV) since a given experimental design is only able to probe masses in a small range. There are several reasons for this. First, it is the frequency range of the experiment that dictates the range of axion masses that can be probed. ADMX is suited to masses eV and has currently set constraints between Asztalos et al. (2010); Hoskins et al. (2011). Larger masses require adjustments to the cavity and amplification technology Slocum et al. (2015); Baker et al. (2012). The Yale Wright Laboratory experiment of Refs. Brubaker et al. (2016); Al Kenany et al. (2016) for example operates between 5  25 GHz (corresponding to 20  100 eV) and is the first to set limits for eV over a 100 MHz range. A number of experimental challenges are present in designing experiments for different mass windows. For higher resonant frequencies the effective volume of the cavity falls off quickly as meaning the cavity geometry must be revised to preserve form factors and thus maintain the sensitivity of the experiment. There are also limitations on the frequency ranges for which the SQUID amplification technology is useful meaning new techniques must be developed such as Josephson parametric amplifiers Al Kenany et al. (2016) for the GHz range. For masses towards 40  400 eV the dielectric disk setup of MADMAX Caldwel et al. (2016); Millar et al. (2016) has been designed and avoids the restriction placed on resonators brought about by the dependence on the cavity volume. Smaller masses eV may also be accessible with nuclear magnetic resonancebased experiments such as CASPEr Budker et al. (2014); Graham and Rajendran (2013) or alternative designs with resonant and broadband readout circuits Kahn et al. (2016), and LC circuits Sikivie et al. (2014).
Ultimately the prospects for axion astronomy will depend on the success of one of the aforementioned search strategies, at which point the development of the optimum technology to measure dark matter axionphoton conversion can begin. In addition to the annual modulation signal, which we have shown to be powerful for making more accurate measurements of some astrophysical parameters, it may also be beneficial to search for possible direction dependent methods (e.g., Refs. Irastorza and Garcia (2012); Horns et al. (2013); Jaeckel and Knirck (2016)) as the angular signature of a dark matter signal has been shown to encode much astrophysical information in the context of WIMPs Lee and Peter (2012); O’Hare and Green (2014). However in any of these possible scenarios the methods developed in this study will be a valuable step in progressing toward axion astronomy.
Acknowledgements.
The authors are grateful for useful correspondence with Igor Irastorza and Benjamin Brubaker. C.A.J.O. is supported by a United Kingdom Science and Technology Facilities Council (STFC) studentship. A.M.G. acknowledges support from STFC Grant No. ST/L000393/1.Appendix A Galactic to Laboratory transformation
Here we briefly summarise the coordinate transformation used to rotate particle velocities from the Galactic system into the rest frame of the laboratory.
The Galactic coordinate system is defined such that points towards the Galactic center, points in the plane of the Galaxy towards the direction of Galactic rotation and points towards the Galactic North pole. We define the Laboratory coordinate system which point towards the North, West and zenith respectively. To move between these coordinate systems we also need an intermediate step in the geocentric equatorial frame , where and point towards the celestial equator with right ascensions of 0 and 90 respectively and points to the celestial north pole.
We transform vectors (e.g., VL2 particle velocities) from the Galactic to the laboratory frame with the following transformation,
(45) 
where the transformation from the Galactic to equatorial system is encoded in the matrix Bozorgnia et al. (2012),
(46)  
and the transformation to the laboratory frame is given by,
In which we have used for the Earth latitude of the laboratory. We use for the local apparent sidereal time expressed in degrees which is related to the Julian day (JD) by the following,
(48)  
where is the longitude of the laboratory location. We also must convert the Julian day to Universal Time (UT) using,
(49) 
References
 Peccei and Quinn (1977) R. D. Peccei and H. R. Quinn, Phys. Rev. Lett. 38, 1440 (1977).
 Kim and Carosi (2010) J. E. Kim and G. Carosi, Rev. Mod. Phys. 82, 557 (2010), arXiv:0807.3125 [hepph] .
 Svrcek and Witten (2006) P. Svrcek and E. Witten, JHEP 06, 051 (2006), arXiv:hepth/0605206 [hepth] .
 Arvanitaki et al. (2010) A. Arvanitaki, S. Dimopoulos, S. Dubovsky, N. Kaloper, and J. MarchRussell, Phys. Rev. D81, 123530 (2010), arXiv:0905.4720 [hepth] .
 Cicoli et al. (2012) M. Cicoli, M. Goodsell, and A. Ringwald, JHEP 10, 146 (2012), arXiv:1206.0819 [hepth] .
 Arias et al. (2012) P. Arias, D. Cadamuro, M. Goodsell, J. Jaeckel, J. Redondo, and A. Ringwald, JCAP 1206, 013 (2012), arXiv:1201.5902 [hepph] .
 Davis (1986) R. L. Davis, Phys. Lett. B180, 225 (1986).
 Hiramatsu et al. (2012) T. Hiramatsu, M. Kawasaki, K. Saikawa, and T. Sekiguchi, Phys. Rev. D85, 105020 (2012), [Erratum: Phys. Rev.D86,089902(2012)], arXiv:1202.5851 [hepph] .
 Ipser and Sikivie (1983) J. Ipser and P. Sikivie, Phys. Rev. Lett. 50, 925 (1983).
 Wantz and Shellard (2010) O. Wantz and E. P. S. Shellard, Phys. Rev. D82, 123508 (2010), arXiv:0910.1066 [astroph.CO] .
 Marsh (2016) D. J. E. Marsh, Phys. Rept. 643, 1 (2016), arXiv:1510.07633 [astroph.CO] .
 Raffelt (1986) G. G. Raffelt, Phys. Lett. B166, 402 (1986).
 Isern et al. (2008) J. Isern, E. GarciaBerro, S. Torres, and S. Catalan, Astrophys. J. 682, L109 (2008), arXiv:0806.2807 [astroph] .
 Berenji et al. (2016) B. Berenji, J. Gaskins, and M. Meyer, Phys. Rev. D93, 045019 (2016), arXiv:1602.00091 [astroph.HE] .
 Raffelt (1999) G. G. Raffelt, Ann. Rev. Nucl. Part. Sci. 49, 163 (1999), arXiv:hepph/9903472 [hepph] .
 Burrows et al. (1989) A. Burrows, M. S. Turner, and R. P. Brinkmann, Phys. Rev. D39, 1020 (1989).
 Mayle et al. (1988) R. Mayle, J. R. Wilson, J. R. Ellis, K. A. Olive, D. N. Schramm, and G. Steigman, Phys. Lett. B203, 188 (1988).
 Sikivie (1983) P. Sikivie, 11th International Symposium on Lepton and Photon Interactions at High Energies Ithaca, New York, August 49, 1983, Phys. Rev. Lett. 51, 1415 (1983), [Erratum: Phys. Rev. Lett.52,695(1984)].
 Zioutas et al. (2005) K. Zioutas et al. (CAST), Phys. Rev. Lett. 94, 121301 (2005), arXiv:hepex/0411033 [hepex] .
 Armengaud et al. (2014) E. Armengaud et al., JINST 9, T05002 (2014), arXiv:1401.3233 [physics.insdet] .
 Asztalos et al. (2010) S. J. Asztalos et al. (ADMX), Phys. Rev. Lett. 104, 041301 (2010), arXiv:0910.5914 [astroph.CO] .
 Barbieri et al. (2016) R. Barbieri, C. Braggio, G. Carugno, C. S. Gallo, A. Lombardi, A. Ortolan, R. Pengo, G. Ruoso, and C. C. Speake, (2016), arXiv:1606.02201 [hepph] .
 Crescini et al. (2017) N. Crescini, C. Braggio, G. Carugno, P. Falferi, A. Ortolan, and G. Ruoso, Nucl. Instrum. Meth. A842, 109 (2017), arXiv:1606.04751 [physics.insdet] .
 Ruoso et al. (2016) G. Ruoso, A. Lombardi, A. Ortolan, R. Pengo, C. Braggio, G. Carugno, C. S. Gallo, and C. C. Speake, Proceedings, 14th International Conference on Topics in Astroparticle and Underground Physics (TAUP 2015): Torino, Italy, September 711, 2015, J. Phys. Conf. Ser. 718, 042051 (2016), arXiv:1511.09461 [hepph] .
 Chung (2016) W. Chung, Proceedings, 15th Hellenic School and Workshops on Elementary Particle Physics and Gravity (CORFU2015): Corfu, Greece, September 125, 2015, PoS CORFU2015, 047 (2016).
 Caldwel et al. (2016) A. Caldwel, G. Dvali, B. Majorovits, A. Millar, G. Raffelt, J. Redondo, O. Reimann, F. Simon, and F. Steffen (MADMAX Working Group), (2016), arXiv:1611.05865 [physics.insdet] .
 Millar et al. (2016) A. J. Millar, G. G. Raffelt, J. Redondo, and F. D. Steffen, (2016), arXiv:1612.07057 [hepph] .
 BÃ¤hre et al. (2013) R. BÃ¤hre et al., JINST 8, T09001 (2013), arXiv:1302.5647 [physics.insdet] .
 Van Bibber et al. (1987) K. Van Bibber, N. R. Dagdeviren, S. E. Koonin, A. Kerman, and H. N. Nelson, Phys. Rev. Lett. 59, 759 (1987).
 Kahn et al. (2016) Y. Kahn, B. R. Safdi, and J. Thaler, Phys. Rev. Lett. 117, 141801 (2016), arXiv:1602.01086 [hepph] .
 Sikivie et al. (2014) P. Sikivie, N. Sullivan, and D. B. Tanner, Phys. Rev. Lett. 112, 131301 (2014), arXiv:1310.8545 [hepph] .
 Beck (2015) C. Beck, Phys. Dark Univ. 78, 6 (2015), arXiv:1403.5676 [hepph] .
 Beck (2013) C. Beck, Phys. Rev. Lett. 111, 231801 (2013), arXiv:1309.3790 [hepph] .
 Arvanitaki and Geraci (2014) A. Arvanitaki and A. A. Geraci, Phys. Rev. Lett. 113, 161801 (2014), arXiv:1403.1290 [hepph] .
 Budker et al. (2014) D. Budker, P. W. Graham, M. Ledbetter, S. Rajendran, and A. Sushkov, Phys. Rev. X4, 021030 (2014), arXiv:1306.6089 [hepph] .
 Graham and Rajendran (2013) P. W. Graham and S. Rajendran, Phys. Rev. D88, 035023 (2013), arXiv:1306.6088 [hepph] .
 Aprile et al. (2014) E. Aprile et al. (XENON100), Phys. Rev. D90, 062009 (2014), arXiv:1404.1455 [astroph.CO] .
 Graham et al. (2015) P. W. Graham, I. G. Irastorza, S. K. Lamoreaux, A. Lindner, and K. A. van Bibber, Ann. Rev. Nucl. Part. Sci. 65, 485 (2015), arXiv:1602.00039 [hepex] .
 Krauss et al. (1985) L. Krauss, J. Moody, F. Wilczek, and D. E. Morris, Phys. Rev. Lett. 55, 1797 (1985).
 Turner (1990) M. S. Turner, Phys. Rev. D42, 3572 (1990).
 Ling et al. (2004) F.S. Ling, P. Sikivie, and S. Wick, Phys. Rev. D70, 123503 (2004), arXiv:astroph/0405231 [astroph] .
 Vergados and Semertzidis (2017) J. D. Vergados and Y. Semertzidis, Nucl. Phys. B915, 10 (2017), arXiv:1601.04765 [astroph.GA] .
 Irastorza and Garcia (2012) I. G. Irastorza and J. A. Garcia, JCAP 1210, 022 (2012), arXiv:1207.6129 [physics.insdet] .
 Duffy et al. (2006) L. D. Duffy, P. Sikivie, D. B. Tanner, S. J. Asztalos, C. Hagmann, D. Kinion, L. J. Rosenberg, K. van Bibber, D. B. Yu, and R. F. Bradley (ADMX), Phys. Rev. D74, 012006 (2006), arXiv:astroph/0603108 [astroph] .
 Hoskins et al. (2016) J. Hoskins et al., Phys. Rev. D94, 082001 (2016).
 Hoskins et al. (2011) J. Hoskins et al., Phys. Rev. D84, 121302 (2011), arXiv:1109.4128 [astroph.CO] .
 Sloan et al. (2016) J. V. Sloan et al., Phys. Dark Univ. 14, 95 (2016).
 Strigari and Trotta (2009) L. E. Strigari and R. Trotta, JCAP 0911, 019 (2009), arXiv:0906.5361 [astroph.HE] .
 Hogan and Rees (1988) C. J. Hogan and M. J. Rees, Phys. Lett. B205, 228 (1988).
 Diemand et al. (2007) J. Diemand, M. Kuhlen, and P. Madau, Astrophys. J. 667, 859 (2007), arXiv:astroph/0703337 [astroph] .
 Hong et al. (2014) J. Hong, J. E. Kim, S. Nam, and Y. Semertzidis, (2014), arXiv:1403.1576 [hepph] .
 McAllister et al. (2016) B. T. McAllister, S. R. Parker, and M. E. Tobar, Phys. Rev. Lett. 116, 161804 (2016), [Erratum: Phys. Rev. Lett.117,no.15,159901(2016)], arXiv:1607.01928 [hepph] .
 Kim (1979) J. E. Kim, Phys. Rev. Lett. 43, 103 (1979).
 Shifman et al. (1980) M. A. Shifman, A. I. Vainshtein, and V. I. Zakharov, Nucl. Phys. B165, 45 (1980).
 Dine et al. (1981) M. Dine, W. Fischler, and M. Srednicki, Phys. Lett. B104, 199 (1981).
 Zhitnitsky (1980) A. R. Zhitnitsky, Sov. J. Nucl. Phys. 31, 260 (1980), [Yad. Fiz.31,497(1980)].
 Mayet et al. (2016) F. Mayet et al., Phys. Rept. 627, 1 (2016), arXiv:1602.03781 [astroph.CO] .
 Schoenrich et al. (2010) R. Schoenrich, J. Binney, and W. Dehnen, Mon. Not. Roy. Astron. Soc. 403, 1829 (2010), arXiv:0912.3693 [astroph.GA] .
 Jackson (1998) J. D. Jackson, Classical Electrodynamics (Wiley, 1998).
 Daw (1998) E. J. Daw, A search for halo axions, Ph.D. thesis, MIT (1998).
 Hotz (2013) M. T. Hotz, A SQUIDBased RF Cavity Search for Dark Matter Axions, Ph.D. thesis, Washon U. (2013).
 Brubaker et al. (2016) B. M. Brubaker et al., (2016), arXiv:1610.02580 [astroph.CO] .
 Dicke et al. (1946) R. H. Dicke, R. Beringer, R. L. Kyhl, and A. B. Vane, Phys. Rev. 70, 340 (1946).
 Savage et al. (2006) C. Savage, K. Freese, and P. Gondolo, Phys. Rev. D74, 043531 (2006), arXiv:astroph/0607121 [astroph] .
 Wilks (1938) S. S. Wilks, Annals Math. Statist. 9, 60 (1938).
 Feroz and Hobson (2008) F. Feroz and M. P. Hobson, Mon. Not. Roy. Astron. Soc. 384, 449 (2008), arXiv:0704.3704 [astroph] .
 Feroz et al. (2009) F. Feroz, M. P. Hobson, and M. Bridges, Mon. Not. Roy. Astron. Soc. 398, 1601 (2009), arXiv:0809.3437 [astroph] .
 Feroz et al. (2013) F. Feroz, M. P. Hobson, E. Cameron, and A. N. Pettitt, (2013), arXiv:1306.2144 [astroph.IM] .
 Purcell et al. (2012) C. W. Purcell, A. R. Zentner, and M.Y. Wang, JCAP 1208, 027 (2012), arXiv:1203.6617 [astroph.GA] .
 Kolb and Tkachev (1993) E. W. Kolb and I. I. Tkachev, Phys. Rev. Lett. 71, 3051 (1993), arXiv:hepph/9303313 [hepph] .
 Kolb and Tkachev (1994a) E. W. Kolb and I. I. Tkachev, Phys. Rev. D49, 5040 (1994a), arXiv:astroph/9311037 [astroph] .
 Kolb and Tkachev (1994b) E. W. Kolb and I. I. Tkachev, Phys. Rev. D50, 769 (1994b), arXiv:astroph/9403011 [astroph] .
 Kolb and Tkachev (1996) E. W. Kolb and I. I. Tkachev, Astrophys. J. 460, L25 (1996), arXiv:astroph/9510043 [astroph] .
 Berezinsky et al. (2013) V. S. Berezinsky, V. I. Dokuchaev, and Yu. N. Eroshenko, JCAP 1311, 059 (2013), arXiv:1308.6742 [astroph.CO] .
 Tinyakov et al. (2016) P. Tinyakov, I. Tkachev, and K. Zioutas, JCAP 1601, 035 (2016), arXiv:1512.02884 [astroph.CO] .
 Chang et al. (1999) S. Chang, C. Hagmann, and P. Sikivie, Phys. Rev. D59, 023505 (1999), arXiv:hepph/9807374 [hepph] .
 Zurek et al. (2007) K. M. Zurek, C. J. Hogan, and T. R. Quinn, Phys. Rev. D75, 043511 (2007), arXiv:astroph/0607341 [astroph] .
 Goerdt et al. (2007) T. Goerdt, O. Y. Gnedin, B. Moore, J. Diemand, and J. Stadel, Mon. Not. Roy. Astron. Soc. 375, 191 (2007), arXiv:astroph/0608495 [astroph] .
 Schneider et al. (2010) A. Schneider, L. Krauss, and B. Moore, Phys. Rev. D82, 063525 (2010), arXiv:1004.5432 [astroph.GA] .
 Kuijken and Gilmore (1989) K. Kuijken and G. Gilmore, Mon. Not. Roy. Astron. Soc. 239, 605 (1989).
 Kuhlen et al. (2012) M. Kuhlen, M. Lisanti, and D. N. Spergel, Phys. Rev. D86, 063505 (2012), arXiv:1202.0007 [astroph.GA] .
 Schaller et al. (2016) M. Schaller, C. S. Frenk, A. Fattahi, J. F. Navarro, K. A. Oman, and T. Sawala, Mon. Not. Roy. Astron. Soc. 461, L56 (2016), arXiv:1605.02770 [astroph.GA] .
 Duffy and Sikivie (2008) L. D. Duffy and P. Sikivie, Phys. Rev. D78, 063508 (2008), arXiv:0805.4556 [astroph] .
 Frandsen et al. (2012) M. T. Frandsen, F. Kahlhoefer, C. McCabe, S. Sarkar, and K. SchmidtHoberg, JCAP 1201, 024 (2012), arXiv:1111.0292 [hepph] .
 Fox et al. (2011) P. J. Fox, G. D. Kribs, and T. M. P. Tait, Phys. Rev. D83, 034007 (2011), arXiv:1011.1910 [hepph] .
 Gondolo and Gelmini (2012) P. Gondolo and G. B. Gelmini, JCAP 1212, 015 (2012), arXiv:1202.6359 [hepph] .
 Del Nobile et al. (2013) E. Del Nobile, G. B. Gelmini, P. Gondolo, and J.H. Huh, JCAP 1310, 026 (2013), arXiv:1304.6183 [hepph] .
 Fox et al. (2014) P. J. Fox, Y. Kahn, and M. McCullough, JCAP 1410, 076 (2014), arXiv:1403.6830 [hepph] .
 Feldstein and Kahlhoefer (2014) B. Feldstein and F. Kahlhoefer, JCAP 1408, 065 (2014), arXiv:1403.4606 [hepph] .
 Kahlhoefer and Wild (2016) F. Kahlhoefer and S. Wild, (2016), arXiv:1607.04418 [hepph] .
 Gelmini et al. (2016) G. B. Gelmini, J.H. Huh, and S. J. Witte, (2016), arXiv:1607.02445 [hepph] .
 Peter (2011) A. H. G. Peter, Phys. Rev. D83, 125029 (2011), arXiv:1103.5145 [astroph.CO] .
 Kavanagh and Green (2013) B. J. Kavanagh and A. M. Green, Phys. Rev. Lett. 111, 031302 (2013), arXiv:1303.6868 [astroph.CO] .
 Kavanagh (2014) B. J. Kavanagh, Phys. Rev. D89, 085026 (2014), arXiv:1312.1852 [astroph.CO] .
 Kavanagh and O’Hare (2016) B. J. Kavanagh and C. A. O’Hare, Phys. Rev. D94, 123009 (2016), arXiv:1609.08630 [astroph.CO] .
 Slocum et al. (2015) P. L. Slocum, O. K. Baker, J. L. Hirshfield, Y. Jiang, A. T. Malagon, A. J. Martin, S. Shchelkunov, and A. Szymkowiak, Nucl. Instrum. Meth. A770, 76 (2015), arXiv:1410.1807 [physics.insdet] .
 Baker et al. (2012) O. K. Baker, M. Betz, F. Caspers, J. Jaeckel, A. Lindner, A. Ringwald, Y. Semertzidis, P. Sikivie, and K. Zioutas, Phys. Rev. D85, 035018 (2012), arXiv:1110.2180 [physics.insdet] .
 Al Kenany et al. (2016) S. Al Kenany et al., (2016), arXiv:1611.07123 [physics.insdet] .
 Horns et al. (2013) D. Horns, J. Jaeckel, A. Lindner, A. Lobanov, J. Redondo, and A. Ringwald, JCAP 1304, 016 (2013), arXiv:1212.2970 [hepph] .
 Jaeckel and Knirck (2016) J. Jaeckel and S. Knirck, JCAP 1601, 005 (2016), arXiv:1509.00371 [hepph] .
 Lee and Peter (2012) S. K. Lee and A. H. G. Peter, JCAP 1204, 029 (2012), arXiv:1202.5035 [astroph.CO] .
 O’Hare and Green (2014) C. A. J. O’Hare and A. M. Green, Phys. Rev. D90, 123511 (2014), arXiv:1410.2749 [astroph.CO] .
 Bozorgnia et al. (2012) N. Bozorgnia, G. B. Gelmini, and P. Gondolo, JCAP 1206, 037 (2012), arXiv:1111.6361 [astroph.CO] .