Axion astronomy with microwave cavity experiments

Terrestrial searches for the conversion of dark matter axions or axion-like 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 N-body 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 post-discovery era are able to perform"axion astronomy".


I. INTRODUCTION
Axions are light pseudoscalar particles that appear in the solution of Peccei and Quinn [1,2] to explain the unnatural absence of CP-violation in quantum chromodynamics (QCD). More modern motivation from the landscape of axion-like particles (ALPs) appearing in string theory [3][4][5], inspires the generalisation of axions to light pseudoscalars with a theoretical origin outside of the original Peccei-Quinn solution. Such ALPs can cover an extremely wide range of masses and couplings to the standard model [6]. Axions and ALPs are an attractive cold dark matter candidate and can be produced in the early Universe through a variety of non-thermal mechanisms such as vacuum misalignment or the decay of topological defects [7,8] in ways that are consistent with the known cosmological abundance and phenomenology of dark matter [9,10]. For a recent review of axion cosmology see e.g., Ref. [11]. In the following we use the term 'axion' to refer to both the QCD axion and generic axion-like particles.
Certain axion mass and coupling ranges can be ruled out with various astrophysical observations such as the cooling of white dwarfs [12,13], neutron star interactions [14], the lifetimes of horizontal branch stars in globular clusters [15] and supernovae neutrinos [16,17]. Axions may also be observable in the lab. Experimental * Electronic address: ciaran.ohare@nottingham.ac.uk 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 [18]. 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 [19] (and the planned IAXO [20]) searching for axions produced inside the Sun, and haloscopes such as ADMX searching for Galactic dark matter axions [21]. 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 signal-to-noise level of the experiment. Planned dark matter axion experiments such as QUAX [22][23][24], CULTASK [25] and the layered dielectric haloscope MADMAX [26,27] 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 lab-based experiments such as the Any Light Particle Search [28] using the technique known as light-shiningthrough-a-wall [29]. In haloscopes beyond ADMX, it may be possible to search for lighter axions with broadband readout circuits [30] or LC circuits [31] as well as heavier masses in the meV range with Josephson junctions [32,33]. Other couplings such as those to nuclei [34] can be probed in nuclear magnetic resonance experiments such as CASPEr [35,36] and the coupling to electrons can be constrained using conventional dark matter direct de-tection searches for weakly interacting massive particles (WIMPs) [37]. For a recent review of axion experiments see for example Ref. [38].
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 axion-photon 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 [39]. 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' [40][41][42]. There is also a potentially exploitable O(1) 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 [43].
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 axion-photon 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 [44][45][46] or applying constraints to different halo models [42,47]. 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 well-established methods for extracting astrophysical information in the case of WIMP direct detection experiments e.g., Ref. [48]. 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 ADMX-like experiments might be for doing axion astronomy. We discuss this idea in the context of simple analytic halo models, distributions from N-body simulations, and minicluster streams -a phenomenon unique to axion dark matter [49].
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 N-body data from the Via Lactea II (VL2) [50] 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 axion-photon conversion power inside a microwave cavity. Full details of these calculations can be found in Refs. [39,51,52]. 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 Breit-Wigner function.
The axion to two photon coupling is given by the formula, which includes the fine structure constant α and the axion decay constant f a . The dimensionless coupling g γ is, In which E/N is the ratio of the colour axion anomaly to the electromagnetic axion anomaly and z is the ratio of the up and down quark masses. The value of this constant is model dependent: E/N = −0.97 for the KSVZ model [53,54] and 0.36 for the DFSZ model [55,56] for example. In the interest of model independence and to generalise to ALPs we express the interaction in terms of the coupling g aγγ . The effective Lagrangian for axions coupled to electromagnetism is, where F µν is the electromagnetic field strength tensor, andF µν = 1 2 µνρσ F ρσ its dual. The axion potential V (a) is provided by QCD instanton effects and can be approximated with a simple mass term 1 2 m 2 a a 2 . The axionic charge density and the electromagnetic current density are written as ρ q and j µ . Writing F µνF µν = −4 E · B we then see the axion-photon interaction in terms of electric and magnetic field strengths, This interaction modifies Maxwell's equations to include an additional axion current, However these equations simplify for the setup we consider here. Firstly we assume the axion field has no spatial dependence on laboratory scales (∇a = 0). 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: ρ q = 0 and j µ = 0. This results in the following simple set of equations, Under the above assumptions the equation of motion for the axion field is, Dark matter axions in the Milky Way undergo essentially no interactions, so in a quadratic potential V (a) 1 2 m 2 a a 2 , the field oscillates coherently at the axion mass a(t) = a 0 e imat ≡ a 0 e iωt . However due to thermalisation in the Milky Way the coherence of the oscillations is spoiled slightly by a dispersion in axion velocities: ω = m a (1 + 1 2 v 2 + O(v 4 )). One can account for this by moving to a Fourier description of the field, written as A(ω), where T is some large reference time used to take the averages. The quantity |A(ω)| 2 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, 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 f lab (v) in the laboratory frame (i.e., f lab (v) = f gal (v + v lab )) by temporarily introducing a number density, where dn is the number density of "particles" with speeds between v and v + dv. The constant n 0 is found from integrating dn over all velocities and is used to define the local axion number density n 0 ≡ ρ a /m a . This allows the connection to a classical field oscillating at m a , which should have a 2 (t) = n/m a , to be made [39]. An expression for the axion power spectrum |A(ω)| 2 can be obtained by satisfying Parseval's relation and changing variables from ω to v, we can then substitute for d a 2 (t) /dv using, where this expression clarifies the distinction between a 3-dimensional velocity distribution f (v) and its 1dimensional speed distribution f (v). Hence, the formula for the axion power spectrum on Earth can be written as, The simplest assumption for a dark matter halo is the Standard Halo Model (SHM) which is spherically symmetric with a 1/r 2 density profile and yields a Maxwell-Boltzmann velocity distribution, To simplify the following expressions we use the peak velocity v 0 for the shape of the distribution, however it can also be written in terms of a velocity dispersion, v 2 0 = 2σ 2 v . The speed distribution follows from an integral over all directions, However the power spectrum should reflect the fact that we observe f lab not f gal so we need to compute the speed distribution in the laboratory frame. To do this we make the transformation v → v − v lab (t) which yields for the velocity and speed distributions, and, Since we are now in the moving laboratory frame a time dependence appears in the speed distribution. Finally after changing variables to ω = m a (1 + v 2 /2) we arrive at the axion power spectrum. The axion power spectrum must be 0 for ω < m a which is enforced by requiring that f (v) be real. For use in later examples we also define the velocity distribution for streams f str lab (v) which are spatially and kinematically localised substructure components of the dark matter halo. Their velocity distribution can also be described with a Maxwellian, where the stream is parameterised by its Galactic frame velocity v str ∼ O(100 km s −1 ) and dispersion σ str ∼ O(1 km s −1 ). We will assume that the stream comprises some fraction ρ str /ρ a 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, 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. [57] for a review of these calculations). In the SHM the velocity of the LSR is usually written as v 0 = (0, v 0 , 0) in Galactic co-ordinates, with v 0 = 220 km s −1 . An often quoted value for the peculiar velocity of the Sun from Schoenrich et al. [58] is v pec = (11.1, 12.24, 7.25) km s −1 with roughly 1 km s −1 sized systematic errors.

III. RESONANCE POWER
We model a microwave cavity experiment with a static uniform magnetic field B 0 maintained inside a cylindrical cavity of radius R and length L, with radial, azimuthal and vertical co-ordinates labelled (r,φ,ẑ) respectively. The magnetic field is generated by a solenoid with current density in theφ-direction. The electric and magnetic fields we write as, where Θ(r) is the Heaviside step function, I is the current and n L is the number of wire turns in the solenoid per unit length. For convenience we use the magnitude of the magnetic field B 0 = n L I in the following expressions.
In the cylindrical cavity design the important cavity mode orientations are the TM 0l0 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, 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 ω 0 ). Following the derivation of Ref. [51], the axion power is calculated by evaluating the following integral over the volume of the cavity V , where U is the energy stored in the electric and magnetic fields inside the cavity. This expression introduces the quality factor Q 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 |A(ω)| 2 ) one arrives at, where χ 0l is the l-th zero of the 0th Bessel function of the first kind. We have also defined T (ω), which is a Lorentzian that describes the loss in power off resonance, Usually the haloscope power is written in terms of a cavity form factor. For the transverse magnetic field considered here (TM 0l0 ) this is written C 0l0 = 4/χ 2 0l 1 . We are principally interested in the TM 010 mode which has C 010 = 0.69. ADMX can tune the TM 010 mode from roughly 500 MHz to 900 MHz [21]. In general the electric field of the TM nlm mode can be written [59], In which, E(t) is the time dependent component of the field, J m is a Bessel function, x ml is the lth root of J m (x) = 0, R is the cavity radius and L is the cavity length. Modes with n = 0 and m = 0 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 ω 0 = ω a m a and use a Breit-Wigner approximation for the axion power spectrum with an analogous Q-factor Q a ∼ ω/∆ω ∼ 10 6 (this permits an analytic evaluation of the integral in Eq. (32)). We also introduce the axion density by writing a 2 (t) = ρ a /m 2 a . Resulting ultimately in, where we have restored the factors of , c and ε 0 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 min(Q, Q a ) arises from the integral of two Breit-Wigner 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 10 −22 W as is usually quoted, 1 Other mode orientations, the transverse electric (TE nlm ) and transverse electromagnetic (TEM nlm ) modes both have no axial electric field meaning they have negligible form factors.  Table I. The purple line in the frequency-time plane shows the evolution of the frequency at which the power is maximised: 2πνmax = ma(1 + v 2 lab /2) and 2πνmax = ma(1 + |v lab − vstr| 2 /2) for the Maxwellian halo and stream respectively.

IV. MOCK EXPERIMENT
Our simulation is an approximation of the current ADMX setup. We list a set of benchmark experimental parameters in Table I. 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.
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 signal-to-noise 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: ν a = 842.0 MHz (= 3.4671 µeV) and g aγγ = 10 −15 GeV −1 . 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 [60][61][62]. The signal-to-noise ratio for a haloscope experiment of duration τ , is set by the Dicke radiometer equation [63] where ∆ν a is the bandwidth of the axion signal and T S is the noise temperature.
Our mock experiment consists of a long total running time τ tot 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 N Fourier transformed timestream samples of duration ∆τ . The Fourier transform of a given sample is a power spectrum with frequency resolution ∆ν = 1/∆τ . The noise we simulate as Johnson white noise which has rms power P N = k B T S ∆ν inside a given frequency bin with an exponential distribution [44]. The noise power spectrum of the average of N = τ /∆τ individual exponential power spectra corresponding to the N 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 P N and standard deviation P N /N = P N / √ τ ∆ν. The full dataset then consists of a total number of N τ = τ tot /τ 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 v rev (t) 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 [64]. 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 I). The annual modulation of the signal is indicated by the purple line labelled ν max .
We base our likelihood on a χ 2 statistic which measures the offset between the observed value of power P ij obs , and the expected power (signal + rms noise) P ij a + P N in each bin, where i and j label frequency and time bins respectively, where the sums run over N ν = (ν max − ν min )/∆ν frequency bins and N τ = τ tot /τ time bins. The error σ N is given by the suppressed rms noise power P N / √ τ ∆ν. We then construct a likelihood based on this statistic. Mathematically the likelihood as a function of a set of parameters given data D is, where we assume m a , g aγγ and P N 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 L astro 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 L N (P N ) 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 1-dimensional intervals when we are only interested in the reconstruction of one parameter or 2-dimensional 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, λ(θ) = −2(ln L(θ) − L(θ)), whereθ are the maximum likelihood estimators. According to Wilks' theorem [65] this is asymptotically χ 2 k distributed for k 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 χ 2 k is the desired confidence level. For example for k = 1 the 68% interval encloses values of λ < 1 and the 95% encloses values of λ < 4.

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 ADMX-like 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 [66][67][68], and set a tolerance of 10 −3 and use between 2 × 10 3 and 10 4 live points depending on the number of parameters being reconstructed.
In Fig. 2 we show the reconstructed axion parameters m a and g aγγ (left) and the astrophysical parameters v lab and σ v (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 1 year/10 days. 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 m > m a due to the fact that the axion power spectrum is only non-zero for ω > m a . 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 −1 level which roughly matches the accuracy of current astronomical observations [58].
With a full annual modulation signal we can also access the 3-dimensional components of v lab . However since v 0 and v pec are summed in the Galactic frame we can only measure directly the x and z components of v pec . The y 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 v 0 . 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 (v 0 +v y pec ) 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 v lab . However, as the signal-to-noise 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 −1 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 τ tot . 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  Table I. We label particular samples which contain prominent streams with arrows and the sample number.
sigma away from their input values. In the case of the axion mass we expect one-sided measurements due to the one-sided 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 m a . The size of the error bars for the remaining parameters decrease roughly as 1/ √ τ tot 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. N-BODY DATA
We can source more realistic examples of dark matter distributions from N-body simulations of Milky Way-like 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 astrophys-ical parameters than standard axion searches and WIMP direct detection.
We use data from the Via Lactea II (VL2) [50] 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 Way-like 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 [69]. 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 f (ω) 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 t 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 v → v − v lab (t).
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, 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 one-to-one   mapping. The frequency of the stream modulation ν(t) is proportional to the quantity |v str − v lab (t)| 2 . 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 ρ str , dispersion σ str , and three components of velocity v str 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 signal-to-noise. 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 signal-to-noise 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. [49,[70][71][72][73][74][75]). Miniclusters are formed in the early Universe from density perturbations in the axion field. Perturbations which can form miniclusters result from various types of non-linear dynamics involved with axion oscillations such as vacuum misalignment or the decay of axion defects such as strings and domain walls [76]. Previous work has predicted the existence of up to ∼ 10 10 pc −3 [75] locally if all of the dark matter was in the form of miniclusters, though a direct encounter would occur less than once every 10 5 years [72]. 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 (t MW ∼ 12 Gyr). It has been estimated in Ref. [75] 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. [75] is a value for the number of expected stream crossings with a density larger than the local smooth halo density ρ a , 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 (∼ 10 −4 km s −1 10 2 km s −1 ), 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, Φ = δρ a /ρ a which is a number typically of order unity. The distribution of values of Φ found from the simulations of Ref. [73] appears to follow a function similar to f Φ (Φ) ∼ Φ 0.75 e −Φ 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, M 1 ∼ 10 −12 M (which is allowed by lensing bounds [77]). Ref. [73] 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, We assume that the miniclusters are spherically symmetric with central density ρ mc (Φ) and radius R mc (Φ, M ). We assume a Maxwellian distribution for the speeds of axions inside a minicluster with a dispersion set by the virial velocity σ mc (Φ, M ) = GM/R mc (Φ, M ). From Ref. [77] we assume the miniclusters have a power law density profile with ρ ∝ r −9/4 for r < R mc but enforce the gradient to turn towards 0 at r = 0 to give a central density of ρ mc (Φ).
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 [78,79], Ref. [75] gives the probability of a particular minicluster being disrupted, where here v is the orbital speed of the minicluster, and n 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, S ⊥ = 35 M pc −2 [80]. Given this, we can just use miniclusters with circular orbits intersecting the Solar position (r ) to evaluate the number of crossings over the age of the Galaxy (t MW ) to be roughly n ∼ 2 t MW /t orb ∼ v t MW /πr ∼ 100. Note that the dependence on v drops out of Eq. (42). This is because although faster miniclusters cross the stellar disk more frequently (∝ v), they are also less likely to encounter a star during a given crossing (∝ 1/v). We also note that P (Φ) has no dependence on M since R mc (Φ, M ) and σ mc (Φ, M ) are both proportional to M −1/3 . A stream can be specified alone by four parameters: the density contrast Φ and mass M of the original minicluster, the age of the stream t, 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 L ∼ σ mc t. Assuming the stream retains the original radius of the minicluster and is simply deformed from a sphere of radius R mc into a cylinder of length L, the density of the axions for a minicluster stream of age t is, Reference [75] 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 ρ a /N ν ∼ 0.001ρ a as this is around the lowest density stream that would be observable in this case. We estimate that if this extrapolation of Ref. [75] is valid then, for t MW = 12 Gyr and M = 10 −12 M , there could be up to N str-x ∼ 100 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 N str-x minicluster stream crossings by selecting samples from the parameter space {Φ, v, ρ str }. First we select values for Φ from the distribution P (Φ)f Φ (Φ). We then select v from an isotropic Maxwell-Boltzmann distribution. Finally we draw a value of ρ str such that the number of stream crossings with ρ str > ρ a follows the function presented in Fig. 2 of Ref. [75]. The length of time taken to cross the stream is then approximately, which is derived from the distance travelled through the stream, 2R mc / sin θ, 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 v lab − v and dispersion σ mc . Also in a given time bin the minicluster stream signal will gain an additional spread in frequency from the change in v lab (t) 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 tan θ < 1 2 z disk /r relative to the plane of the stellar disk, where z disk ∼ 0.3 kpc 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 10 7 km they are short-lived 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 short-lived 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 short-lived 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 lab-based 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 ADMX-like 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 axion-photon 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 −1 , which would improve upon the measurement from local astronomical observations.
• Substructure such as tidal streams appearing in simulations of Milky Way-like 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. [75]. 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 [81], dark disks [82] and caustic rings [83]. 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 [84][85][86][87][88][89][90][91] and developing general parameterisations for the velocity distribution [92][93][94][95]. 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 g 2 aγγ one can extend any of the reconstructions to smaller couplings by scaling up the experiment duration, τ tot , 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 m −1 a , our conclusions still hold for smaller (larger) values of the axion mass if τ tot is increased (decreased) by the same factor. The reverse argument goes for values of the local density since the power is linearly proportional to ρ a . However we must take care in extending these results to axion masses much larger or smaller than O(µ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 < 10 µeV and has currently set constraints between 1.9 µeV < m a < 3.69 µeV [21,46]. Larger masses require adjustments to the cavity and amplification technology [96,97]. The Yale Wright Laboratory experiment of Refs. [62,98] for example operates between 5 -25 GHz (corresponding to 20 -100 µeV) and is the first to set limits for m a > 20µ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 ν −3 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 [98] for the GHz range. For masses towards 40 -400 µeV the dielectric disk setup of MADMAX [26,27] has been designed and avoids the restriction placed on resonators brought about by the dependence on the cavity volume. Smaller masses 10 −(6−9) eV may also be accessible with nuclear magnetic resonance-based experiments such as CASPEr [35,36] or alternative designs with resonant and broadband readout circuits [30], and LC circuits [31].
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 axion-photon 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. [43,99,100]) as the angular signature of a dark matter signal has been shown to encode much astrophysical information in the context of WIMPs [101,102]. However in any of these possible scenarios the methods developed in this study will be a valuable step in progressing toward axion astronomy.
system is encoded in the matrix [103], In which we have used λ lab for the Earth latitude of the laboratory. We use t • lab for the local apparent sidereal time expressed in degrees which is related to the Julian day (JD) by the following, where φ lab is the longitude of the laboratory location. We also must convert the Julian day to Universal Time (UT) using, UT = 24 JD + 0.5 − JD + 0.5 .