Axion minivoids and implications for direct detection

In the scenario in which QCD axion dark matter is produced after inflation, the Universe is populated by large inhomogeneities on very small scales. Eventually, these fluctuations will collapse gravitationally to form dense axion miniclusters that trap up to $\sim$75% of the dark matter within asteroid-mass clumps. Axion miniclusters are physically tiny however, so haloscope experiments searching for axions directly on Earth are much more likely to be probing ``minivoids'' -- the space in between miniclusters. This scenario seems like it ought to spell doom for haloscopes, but while these minivoids might be underdense, they are not totally devoid of axions. Using Schr\"odinger-Poisson and N-body simulations to evolve from realistic initial field configurations, we quantify the extent to which the local ambient dark matter density is suppressed in the post-inflationary scenario. We find that a typical experimental measurement will sample an axion density that is only around 10% of the expected galactic dark matter density. Our results also have implications for experimental campaigns lasting longer than a few years, as well as broadband haloscopes that have sensitivity to transient signatures. We show that for a $\mathcal{O}$(year)-long integration times, the measured dark matter density should be expected to vary by 20--30%.

In the scenario in which QCD axion dark matter is produced after inflation, the Universe is populated by large inhomogeneities on very small scales. Eventually, these fluctuations will collapse gravitationally to form dense axion miniclusters that trap up to ∼75% of the dark matter within asteroid-mass clumps. Axion miniclusters are physically tiny however, so haloscope experiments searching for axions directly on Earth are much more likely to be probing "minivoids"-the space in between miniclusters. This scenario seems like it ought to spell doom for haloscopes, but while these minivoids might be underdense, they are not totally devoid of axions. Using Schrödinger-Poisson and N-body simulations to evolve from realistic initial field configurations, we quantify the extent to which the local ambient dark matter density is suppressed in the post-inflationary scenario. We find that a typical experimental measurement will sample an axion density that is only around 10% of the expected galactic dark matter density. Our results are taken as conservative estimates and have implications for experimental campaigns lasting longer than a few years, as well as broadband haloscopes that have sensitivity to transient signatures. We show that for a O(year)-long integration times, the measured dark matter density should be expected to vary by 20-30%.

I. INTRODUCTION
Axions have rapidly been accumulating interest in both the theoretical [1,2] and experimental [3][4][5] particle physics communities recently. A particularly wellmotivated subset of axion models is the quantumchromodynamics (QCD) axion of Peccei & Quinn [6,7], which solves the puzzle of the missing neutron electric dipole moment [6][7][8][9][10][11][12][13][14]. Simple arguments show that the QCD axion is also a rather minimal explanation for the abundance of dark matter in the Universe [15][16][17][18]. As a result, a flourishing campaign of new experiments now seeks to test whether these particles really are what constitute the invisible halos that envelop galaxies like our own Milky Way.
The primary challenge in axion direct detection is both one of breadth-the available mass range spans potentially m a ∈ [∼ 10 −11 , O(1)] eV [48][49][50]-but also depth in that the couplings to Standard Model particles may be vanishingly small if the axion was proa g.pierobon@unsw.edu.au duced at a very high energy scale. Therefore, the field requires both resonance-based experiments that can search deeply while accessing tiny couplings, but also broadband-sensitive experiments that can search with lower sensitivity over wider mass windows.
Since the parameter space is still wide open, it is useful to think about possible theoretical biases we might have towards one axion model over another. In this regard, a target that is certainly worthy of investigation is the 10 −5 -10 −3 eV axion mass window predicted under the so-called "post-inflationary" scenario. This is the scenario in which the Peccei-Quinn (PQ) symmetry-of which the axion is the associated Goldstone boson-is broken after the end of inflation. The resulting dark matter abundance turns out to be tied precisely to the energy scale of PQ breaking, and thus to the axion mass. Because the symmetry breaking occurs after inflation, large field gradients and topological defects emerge in our Universe. These highly non-linear features are generally not amenable to any analytical treatment, so dedicated cosmological field simulations are required to establish the precise relationship between the axion model parameters and the resulting dark matter abundance. To this end, the first simulation-backed predictions of the QCD axion abundance and mass have been made in recent years [51][52][53][54][55][56][57][58][59], although many unresolved theoretical issues remain.
Independently of the precise value of the axion mass, one generic consequence of the post-inflationary scenario that has been known for some time is the production of dark matter substructures in a galactic halo [60][61][62][63][64][65][66][67]. The large-amplitude yet ultra-small-scale density perturbations left over by the axion field dynamics around the QCD era lead to the formation arXiv:2212.00560v2 of gravitationally-bound axion clumps already during radiation-domination. Previous studies have shown that around 75% of all axions become bound into these socalled "miniclusters" [68], which have masses similar to those of asteroids, i.e., M ∈ [10 −13 , 10 −9 ] M , and radii of the order of an astronomical unit (AU∼ 5 µpc). In contrast with the average dark matter density in the solar neighbourhood (ρ DM ∼ 0.01 M pc −3 ≈ 0.4 GeV cm −3 ) inferred from stellar dynamics on scales 100 pc [69][70][71], these miniclusters are orders of magnitude denser than the average expectation. But this statement also implies that the vast majority of the space outside of miniclusters must be relatively empty.
Importantly, axion miniclusters remain stable during the formation of the first galaxy-sized halos and are only tidally disrupted much later by close encounters with comparatively point-like objects like stars [72][73][74][75][76]. Thus, as a firm prediction of the post-inflationary axion dark matter scenario, the late-time consequences of axion miniclusters cannot be ignored in the context of direct or indirect searches.
The existence of significant axion dark matter substructures in galaxy-sized halos has both positive and negative implications for efforts towards its detection. On the one hand, the presence of small bound clumps potentially opens up new opportunities for indirect detection, such as the search for transient radio signals from axion miniclusters colliding with the intense magnetospheres of neutron stars [77,78]. Another proposed indirect detection technique is gravitational microlensing [62,79,80], although axion miniclusters might be too diffuse to generate a detectable signal [81]. 1 On the other hand, the chances of a direct detection of axion dark matter are significantly reduced if most of the axions are bound in miniclusters. A simple back-of-the-envelope calculation of the encounter rate of an experiment in a galaxy full of miniclusters suggests that we could only expect to be lucky enough to pass through one such substructure once every 100,000 years [3].
While the later tidal disruption of miniclusters may slightly refill the space inside the galaxy with axions and facilitate their detection [83,84], we cannot eliminate the possibility that our experiments sit in a "minivoid" well away from centres of the bound axion clumps. In such a worst-case scenario, it is quite possible that there are simply not enough axions in the space we occupy for terrestrial haloscope searches to work as intended. Our aim in this work, therefore, is to quantify how bad it can get: how much of the space in a galaxy has a deficit of axion dark matter and how big this deficit is. To do so, we shall perform simulations of the axion minicluster formation and evolution. We begin from realistic initial conditions left over from the PQ and QCD phase transitions that account for the decay of topological defects, and the formation and collapse of temporary soliton-like field configurations called axitons [53,55] that lay down the fluctuations from which the miniclusters form. Our simulations then follow the gravitational dynamics of the density fields across matter-radiation equality.
The formation and evolution of axion miniclusters have previously been studied using semi-analytic techniques such as the Press-Schechter [66,67,80] and the Peak-Patch [81,85] formalisms, as well as in full numerical N-body simulations [68,75,86]. Our work extends upon and complements the N-body studies of [68,75,86], and uses a mixture of Schrödinger-Poisson and N-body methods to simulate the formation, growth, and merger of halos made of miniclusters. We furthermore quantify our results not only in terms of the properties of the axion miniclusters, but also the statistics of the axions that are not gravitationally bound in substructures. A visualisation of the axion density field produced in our simulations at two times during the evolution is shown in Fig. 1.
The rest of the article is organised chronologically in cosmic time. We begin in Sec. II with a description of the results of the early-universe lattice simulations that provide us with the requisite initial conditions. We switch to simulations of the axion's gravitational dynamics in Sec. III and study the properties of the resulting axion minivoids in Sec. IV. We discuss the consequences of these results for haloscopes in Sec. V and, finally, conclude in Sec. VI.

A. Minicluster seeds
The initial conditions are of critical importance when studying the later formation of miniclusters, and generating a suitable initial configuration for the field requires simulating nonlinear dynamics of the axion well before gravitational effects become relevant. In the postinflationary scenario, the dynamics of the axion field goes through several distinct eras during which important and nontrivial field configurations emerge. Firstly, global axion strings form as a result of the U (1) PQ spontaneous symmetry breaking at very large temperatures T 10 10 GeV. Axion strings then enter a scaling regime and survive until around the QCD phase transition when the axion mass becomes comparable to the Hubble parameter, H(t). At this point, the axion field begins damped oscillations around the minimum of its potential and domain walls form between strings. We call the redshift or time when this happens the characteristic time, t 1 or z 1 , and can be thought of as a measure of when in cosmic evolution the axion becomes dark matter in earnest.
The characteristic times can be calculated by satisfying the condition m a (t 1 ) = H(t 1 ). This redshift turns out to be very high, and depends weakly on the axion mass . (1) The comoving size of the causal volume at t 1 is, accordingly, the characteristic length scale 2 , where n is a parameter that controls the growth of the axion mass with temperature, i.e., m 2 a ∼ T −n . Assuming n 7 (see, e.g., [87]), we obtain L 1 ∼ m −0.18 a . Therefore, given the fact that m a cannot be arbitrarily varied over many orders of magnitude-as there is only a restrictive mass range that give the correct dark matter abundance-the spatial scale L 1 falls between 0.02 pc and 0.05 pc, i.e., it does not vary substantially.
After a few multiples of t 1 , the string-wall system collapses under the wall tension, 3 during which almost the entirety of the energy of the system is converted into nonrelativistic modes of the axion field. The field can then be treated as an inhomogeneous nonrelativistic fluid. Except for a few r ∼ m −1 a localised regions which form quasi-stable solitonic configurations of the axion field known as axitons, the field is largely in the linear regime. Axitons are expected to disappear when the axion mass reaches its present-day value [55].
The topological defects and the subsequent axitons both leave behind large overdensities in the energy density distribution that act as seeds for gravitational structures to form later in the cosmological evolution. In our recent work [58] we performed high-resolution simulations of this early stage as a function of the axion mass growth index, n, and studied the resulting distribution of minicluster seeds. We continue the evolution of this system here, when the field values are small enough that the nonrelativistic approximation is safe to assume.
For concreteness, we focus only on the QCD axion, adopting the value n = 7. Our initial configurations are the final snapshots of two previous simulations with 4096 3 and 8192 3 lattice sites (see Appendix A for further details), and box side lengths of L = 8L 1 and L = 16L 1 respectively. 4

B. The Schrödinger-Poisson system
In the linear regime, the dynamics of the axion field can be described as the slow free-streaming of almost-frozen relic waves. Gravitational effects, on the other hand, can be accounted for in the weak-field limit, where it is convenient to work in the Newtonian gauge defined by the line element ds 2 is the scale factor in a Friedmann-Lemaître-Robertson-Walker background, and Φ N 1 can be identified with the Newtonian gravitational potential.
Given that the now nonrelativistic axion field is rapidly oscillating with its mass m a as the leading order frequency, we can employ a WKB approximation to write the axion field a(x) as where ψ denotes a slowly varying complex scalar field. In the limit ∂ t ψ m a ψ, the evolution of ψ is governed by the Schrödinger-Poisson equations, where δ a ≡ ρ a / ρ a − 1 is the axion overdensity, ρ a = m a |ψ| 2 is the energy density, and · · · denotes an average over the simulation volume. The set of equations (4) and (5) describes any scalar field in the classical regime, i.e., when the occupation number is very large. Their solutions exhibit wave-like effects such as interference patterns and the formation of soliton solutions, as have been observed in several numerical studies [100][101][102][103][104][105][106][107][108][109] in different cosmological scenarios. These features are expected to show up on scales smaller than the so-called "Jeans wavelength" [110], i.e., or the de Broglie wavelength λ dB = 2π/(m a v a ). For axion masses m a 10 −5 eV, λ J and λ dB are comparable to, or smaller than, the discretisation scale of a typical simulation (L = 8L 1 box size and ∼ 4000 3 points; see, e.g., Fig. 7 in Appendix A). In contrast, the Schrödinger-Vlasov correspondence [111] stipulates that the evolution of ρ a on scales larger than λ J cannot be distinguished from that of collisionless pressureless matter such as standard particle cold dark matter. Therefore, standard N-body simulations can be used to evolve the axion field, provided that wave-like effects appear only on scales smaller than a simulation's spatial resolution [68]. We solve the system of equations (4) and (5) beginning at an initial redshift of z = 5 × 10 11 , until the numerical solution of the Schrödinger-Poisson system becomes unreliable at z ∼ 10 6 (discussed further in Appendix A) and the first miniclusters start to form. From this point onwards, we need to switch to N-body simulations in order to capture the gravitational collapse of axion overdensities; this will be discussed in the next section. Because in addition to free-streaming, we include the effects of linear gravity between the early Universe simulations and the beginning of the N-body simulations in our modelling, this work represents an improvement upon the numerical study of minicluster formation of Ref. [68]. Additional details on the numerical evaluation can be found in Appendix A.

III. N-BODY SIMULATIONS
The largest overdensities are expected to collapse and decouple from the Hubble flow to form gravitationally bound structures at redshifts z ∼ 10 6 . As briefly justified above, from this point onwards we switch to Nbody simulations and continue the numerical evolution of the axion field under gravitational interactions using the N-body code GADGET-4 [112]. We use box sizes corresponding to L/L 1 = 8, 16 and compare our simulation results with those of Ref. [68] that used a larger box size of L/L 1 = 24. Even though the initial conditions of the latter have been evaluated without the Schrödinger-Poisson evolution, we will show that our main results on the minivoid distribution and statistics are largely independent of the choice of initialisation method. An indepth analysis of the choice of initial conditions and their implications on the structure of axion miniclusters will be discussed in a forthcoming publication. Figure 1 shows two visualisations of the resulting density field from our L/L 1 = 8 simulation.

A. Initial conditions
In order to create the initial conditions for the N-body evolution, the axion field values on the lattice need to be converted into particle configurations. One way to do this, explored in Ref. [68], is to create particles of the same mass value and vary the local number of particles placed in the simulation according to the local density contrast. Specifically, a number of particles equal to floor(ρ i / ρ ) are created on every grid point i and each particle is randomly displaced by an amount drawn from a Gaussian distribution with a standard deviation equal to a quarter of the grid spacing. This method is particularly efficient in describing the small-scale structure and the density profiles of the miniclusters, at the cost of a lower resolution in the initially underdense regions that eventually become minivoids.
To better sample the underdense regions, we explore in this work an alternative initialisation method. Here, we place particles in the initial snapshot homogeneously, i.e., one particle per grid cell, but each particle can now carry a different mass that reflects the axion density on the grid cell. Specifically, we scale masses according to m i = δ a,i m av , where m av is the average particle mass. In what follows, we refer to the first mapping method as the "same mass ICs" and the second method as the "different mass ICs". As in Ref. [68], our lattice grids must be down-sampled because of limitations on the total particle number tractable by available computational resources; currently, we can simulate up to 1024 3 particles, which shrinks the lattice grids by a factor of 4-8 per spatial dimension when switching to the N-body method.
We note in passing that we have also implemented another refinement in the initialisation procedure by including initial particle velocities, which had been neglected in Ref. [68]. To do this, we computed the gradient of the phase of ψ, i.e., v i = ∇ arg ψ i /m a , for the ith particle. However, in short, we did not observe many substantial differences between including and not including the velocity information-at least not in the context of the minivoid analysis we shall present shortly.

B. Final simulation time
Having implemented the initial particle realisation, their subsequent evolution under gravitational interaction as tracked by GADGET-4 simply follows wellestablished and widely used N-body techniques. The only remaining aspect that needs to be discussed is the final time imposed upon the simulation by the finite box size.
Because we use boxes of several different sizes, the final simulation redshifts also vary slightly. In the L = 24L 1 simulation box of Ref. [68], the N-body simulation becomes unreliable past z = 99. This can be seen from the dimensionless power spectrum ∆ 2 (z), which becomes O(1) on the largest simulated scales k ∼ 1/L at these late times. The simulated power spectrum on these large scales also begins to deviate from the linearly-evolved power spectrum, given by is the linear growth factor, z i the initial redshift, and z eq is the redshift at matter-radiation equality. Note also that for the smaller simulation box sizes, L/L 1 = 8, 16, used in this work, the final simulation redshift must necessarily be higher, because the largest simulated scales in these boxes are smaller and hence become nonlinear at earlier times.
Here, we choose the final simulation redshift z f by demanding that the ratio between the simulated power spectrum ∆ 2 (z f ) and the linearly-evolved one at k 15 pc −1 should not exceed some threshold p thr , i.e., The choice of p thr = 0.2 gives us the following z f values: If we were to evolve past z f , some small-scale properties such as the shapes of the miniclusters would unlikely be affected. However, statistical statements about the large-scale properties of the system would not be trustworthy. Given that we are quantifying the voids which take up the bulk of the simulation volume, the final redshift is something we must enforce. In the following, we will refer to the box size with its length in units of parsecs. Further details on the N-body simulations can be found in Appendix A.

IV. MINIVOIDS
While most of the mass of dark matter is contained in miniclusters and minihalos, these gravitationally bound objects make up only a tiny fraction of the total volume of the simulation. The "minivoids" that span the space between the miniclusters are where we are mostly likely to be. Thus, as argued in Sec. I, it is of critical importance for direct detection that we quantify the energy density of axions in the voids. To our knowledge, an extensive study on void statistics for the axion post-inflationary scenario has not previously appeared in the literature.

A. Bound fraction
Before we get to the minivoids, let us briefly remark on where the majority of axions actually find themselves. After matter-radiation equality, the evolution of the density field essentially consists of mergers of miniclusters, of typical mass The left panel of Fig. 2 shows f b as a function of time as dashed and dotted lines. For the L = 0.6 pc simulation (dashed blue line), we recover the result of Ref. [68] that the fraction of axions bound in miniclusters reaches a value of f b ∼ 0.75 by around z = 99 and begins to plateau. Comparing this with the bound fraction in our L = 0.2 pc simulations, where we distinguish between runs using the "same mass ICs" and the "different mass ICs" technique described in Sec III A, it is evident from Fig. 2 that the results across these different simulations are broadly similar. In particular, at z = 999 the bound fraction is f b ∼ 0.6, with a 5% difference between the different simulations. This reassures us of consistency between the original "same mass ICs" initialisation procedure of Ref. [68] and the "different mass ICs" approach proposed in this work.

B. Finding minivoids
Given that ∼25% of the axions are defined as being outside of miniclusters, we might already expect na'ıvely that the typical density at a given point outside of miniclusters in the simulation to be ∼ 25% of the average density over the whole simulation volume. However, this ∼ 25% suppression is likely to be an overestimate, as it assumes that the unbound axions are evenly distributed in the box, which is certainly not the case. We therefore analyse the structures that primarily fill the simulation volume directly to obtain a more precise quantification of the suppression of the average axion density.
We find void regions on the grid by first building the density contrast field from the particle snapshots. A Gaussian smoothing filter is then applied on the snapshot, using a range of filter radius between R min = 5∆ x , where ∆ x denotes the mean particle separation distance (or the grid spatial unit ∆ x = L/N ), and R max ∼ L 1 . This choice of R min ensures that our results are not contaminated by discretisation effects on small scales; on the other hand, for R max , we do not expect voids to have comoving radii larger than ∼ 1.5L 1 38 mpc. Void cells are tagged if their average density contrast is smaller than a pre-determined threshold δ thr a , and minivoids are identified as spherical regions that do not overlap with other minivoids previously found. We test thresholds in the range δ thr a ∈ [−0.4, −0.75], with δ thr a = −0.7 as the fiducial value. 6 We catalogue the voids by their radii, in order to study their density profiles and statistical distributions.
As depicted in Fig. 2, we estimate the fraction f v in volume occupied by minivoids, as well as their typical energy density ρ v , by gathering all the minivoids found in the previous steps. One might think that both f v and ρ v depend on the choice of density threshold δ thr a . This is so at the beginning of the simulation, but the dependence on our choice of δ thr a clearly dies away after matter-radiation equality and towards the end of our simulations; see Fig. 8 in Appendix B. In addition, our estimates of f v differ from those of Ref. [68] by less than 5% differences, and the minivoid energy density estimates are in excellent (sub-percent) agreement. In particular, we noticed in the L = 0.2 pc simulation that these results are completely independent of the initialisation mapping method and whether or not particle velocities have been incorporated in the initial conditions. We find convergence to the following values: ρ v / ρ 0.075, f v 0.81, (z = 99). Of the remaining ∼ 20% of the total volume, miniclusters occupy only up to 1%; the rest is in the form of slightly underdense or average-density regions (ρ ∼ ρ ). More details can be found in Appendix B.

C. Minivoid structure
Now we turn to the internal structure of voids by computing their density profiles. The profiles are constructed by integrating over spherical shells around the centre of all voids of a given size, from a minimum radius r min = ∆ x to the maximum radius r max = 10R v , where R v denotes the void radius. The upper panel of Fig. 3 shows the void density profiles of all voids with comoving radii R v = 5.5, 7, 10, and 14 mpc identified in the L = 0.2 pc simulation at redshift z = 999. We find an approximately constant profile within the void radius, and a slow decrease at the discretisation level, r 3∆ x . Immediately beyond r ∼ R v , the void density increases steeply and overshoots the average value ρ , before dropping back down to the average again at large radii well outside the void.
We also compare in the lower panel of Fig. 3 the density profiles of voids with the smallest radius, R v = 5.5 mpc, from the L = 0.2 pc simulation at redshift z = 999, to those of a comparable population from the L = 0.4 pc simulation at z = 499, which have R v = 11 mpc. This close-up reveals spiked regions in the void profiles at radii r < R v that are even more pronounced for larger R v .
We attribute these spikes to small miniclusters that have formed at the earliest times but have not merged into larger structures. Notwithstanding the highly nontrivial substructure on ∼mpc scales, the average density profiles of the minivoids-which incorporates full particle information from the simulation-agree with the result obtained earlier from the void finding procedure, which uses smoothed energy densities on the grid. This is seen by comparing the right panel of Fig. 2 with the bottom panel of Fig. 3 at radii r 3∆ x , where ρ v / ρ 0.1 at z = 999 and ρ v / ρ 0.075 at z 499. Similarly to the halo mass function for the miniclusters, we can describe the minivoid statistics by computing the minivoid size function, defined as the comoving number density of minivoids per differential radius interval dn/dR. Figure 4 shows the distribution of minivoids from the simulation with the largest statistics, namely, the one from Ref. [68]. We find that relatively large void regions start to form even before matter-radiation equality, and by z ∼ 1000 the shape of the minivoid size function does not change anymore. Between R v = 5∆ x and R v = L 1 , the void size function is almost scale-invariant and can be described by a power law, where the parameter values a = 1.43 and b = 3.11 provide a good fit to the simulation results at z = 99, as indicated by the dashed line in Fig. 4. Note that our smaller L = 0.2, 0.4 pc simulations exhibit the same void size functions and a comparable redshift dependence, indicating that these results are generally convergent. See Appendix B for further details.

D. Minivoid fluctuations
Apart from the statistics of minivoids and their profiles, the energy distribution within the simulation volume also provides further insights. We are interested to estimate the degree of variance in the typically-observed dark matter density, in order to guide experiments that will be measuring over long timescales. However, our smallest discretisation scale ∆ x = 0.4 mpc (for L = 0.2 pc, N = 512) is still larger than the distance d year 0.2 mpc travelled by the Solar System around the Milky Way halo in one year at a speed of v ∼ 220 km/s. This implies that, with our current simulation resolution, we cannot resolve the density variations for observation times under one to two years. Nevertheless, we have analysed O(10 5 ) randomly selected trajectories in the final snapshot of the L = 0.2 pc simulation over a distance of 11∆ x 4.5 mpc. For each trajectory, we estimate the density variation as a function of the displacement ∆r from the starting point. Figure 5 shows two measures of the density variation relevant for different types of experiments relying on measurements over long periods of time and on individual short-term measurements.  The first measure of the density variation, shown in the left panel of Fig. 5, is what we call the global density variation, and measures how much the density contrast varies along a trajectory, relative to the average dark matter density over the whole simulation volume ρ , ∆δ(r) = max ρ(r) − min ρ(r) ρ .
As can be seen in Fig. 5, while the typical density sampled by an experiment at any one time is about 10% of the average, the number can in fact be expected to vary between ∼5% and ∼15% over an O(1) year timescale. For longer times, i.e., > 10 years, the variation can grow even more. We will refer to this result again in the next section. A varying density could potentially challenge the interpretation of multiple experiments attempting to observe and then test a putative axion signal over several short measurements at different points in time.
The second measure, shown in the right panel of Fig. 5, is a similar quantity which we simply call the local density variation, and quantifies how the typical size of the density fluctuations changes relative to the average density along a particular trajectory, Here, n r denotes the number of grid points making up the trajectory r, and we use ρ(r) traj to make clear that this is the mean density along that trajectory. We define this measure as a standard deviation-rather than as a range as in Eq. (15)-because we use this in the context of single long measurements, though this choice is arbitrary. By this measure, the right panel of Fig. 5 shows that the measured density can be expected to vary by 20-30% initially for an observing time of O(1) year, and by almost 100% if the measurement is carried out over even longer timescales. This measure is consistent with the previous one, but the interpretation is subtly different, in that we must imagine here a continuous measurement of the axion field lasting long periods of time. This variation would affect the reconstruction of the axion's properties by an amount that goes roughly as the square root of the typical scale shown by e.g. the white line in the right panel of Fig. 5. We discuss the implications of this result further in Sec. V.

E. Minivoid survival
With current numerical methods we are not able to extend the simulations after z ∼ 100. Indeed, it is unlikely that any N-body-based simulation that actually resolves individual miniclusters could ever be evolved to the present day on their own, due to the extreme computational demand placed on the resolution and box size. However, we can attempt to qualitatively describe what happens at later stages in the evolution.
To simplify the discussion, minivoids should be thought of as simply the space not occupied by miniclusters rather than actual entities in and of themselves, so what happens to them is necessarily dictated by whatever happens to the miniclusters. The only processes we can think of that can efficiently disrupt or affect the miniclusters are the tidal interactions with dense objects like stars. If, say, a minicluster had enough energy injected in a close encounter for the axions to become totally unbound, then the axion minicluster would grow into a long stream, spilling its mass into a volume given roughly by ∼ R 2 σt, where R is the radius of the minicluster, σ is its velocity dispersion, and t the time since disruption. This is probably best thought of as the formation of a distinct type of object, a stream (see, e.g., [83,84,114]), but effectively what has happened is that the minivoids have partially been re-filled by axions. On the other hand, one of the key observations from our study is that the miniclusters and minivoids have stopped growing by the end of our simulations, with the void the densities approaching a constant value (see e.g. Fig. 2 on the right). Our expectation, as we will argue, is therefore that the stellar and tidal interactions, would likely only increase the fraction of axions outside of miniclusters. An understanding of how the phase space distribution becomes re-filled, and the role of the minivoids and their substructure at z = 0 is beyond our current scope and is left for future works.  15), which describes how much the density varies along a trajectory relative to the large-scale mean density ρ . The right panel shows a related quantity which we call the local density variation, σρ, defined in Eq. (16), which is the typical size of density fluctuations relative to the mean along that particular trajectory. In both cases, the solid white line corresponds to the median over 10 5 random trajectories along 6 mpc, while the dashed lines enclose the 16th to 84th percentiles of the distribution. The distribution itself is normalised at each grid point and indicated by the purple shading.

V. IMPLICATIONS FOR HALOSCOPES
Having quantified the extent to which the typical local density of axions in a galactic halo is suppressed due to the presence of minivoids, we now wish to evaluate the implications for efforts to search for the dark-matter axion in laboratory experiments. For illustrative purposes, we shall focus on experiments utilising the axion's coupling to the photon, g aγ , as these are the subject of the most experimental activity. In principle, though, our results apply to all direct searches for axions via any other coupling, as long as the post-inflationary scenario is true and the resulting miniclusters are present in our universe-this typically will be the case for axions in the mass range ∼ 10µeV-meV. 7 Hence, experiments like CASPEr-electric [20], as well as the proposals of Refs. [21,22], which purport QCD axion sensitivity in the same mass range via the electric-dipole-moment coupling may also be impacted.
Essentially, all direct searches for axions are based on a model for the local behaviour of the field in terms of coherent oscillations, 7 Axion models with domain wall number greater than 1 may have qualitatively different distributions of miniclusters and minivoids, which is deserving of a dedicated study where φ is an arbitrary phase. We ignore all spatial dependence here as it is unimportant for this discussion. The frequency of the oscillations is given by the kinetic energy in the field ω ≈ m a (1 + v 2 /2), where v is the dark matter speed in units of c. Since the dark matter speed is only a correction of v ∼ 10 −3 and varies with a dispersion, σ v , of approximately the same scale, this implies the field oscillates coherently over m a /m a σ 2 v ∼ 10 6 cycles. The velocity will be drawn from the local velocity distribution, however, so when the field is observed for timescales longer than 10 6 /m a , any resulting signal that is tied to those oscillations will be distributed in frequency according to the same distribution up to a change of variables-this, in haloscope jargon, is called the lineshape.
The majority of axion experiments (and all axion experiments using the photon coupling) couple linearly to the axion field, e.g., L ∼ g aγ aFF . This means that signals will scale as a 2 g 2 aγ . Ignoring O(1) factors specific to individual experiments, the generic scaling of electromagnetic signal power on the quantities of interest is where f (ω i ) is the axion lineshape binned at some arbitrary set of frequency bins ω i . The amplitude α(ω i ) is a stochastic quantity drawn from a Rayleigh distribution [114] that is present if the phases in Eq. (17) are randomly drawn in every coherence time (this may not be the case in the vicinity of miniclusters; however, the FIG. 6. Current constraints on the axion-photon coupling in the µeV-meV mass window, for two different assumptions about the ambient density of axions. In the lighter shades, we show the standard assumption that axions make up 100% of the measured galactic dark matter distribution in the solar neighbourhood of ρDM = 0.45 GeV cm −3 . The darker shade bounds our estimate of the lowest assumption for the typical ambient density of dark matter in the worst-case scenario that the local density of axions is primarily inside of miniclusters, yet experiments probe the space between miniclusters. The red haloscope bounds are taken from Refs. , CAST from Ref. [139,140], globular cluster stellar cooling bound from Ref. [141], a search for axions produced in pulsar polar cap cascades [142], and the neutron star dark matter bound from Ref. [143]. Since the last constraint also relies on the assumption that the axions make up galactic dark matter, we have also applied the rescaling factor as with the haloscopes. Limit data and plots available at Ref. [144]. conservative approach is to assume that it is).
From here we see that experimental sensitivity to g aγ scales as 1/ √ ρ a , implying a factor of ∼3 suppression in limits when comparing typical minivoid densities to the baseline homogeneous assumption that the observable axion density is exactly equal to the astronomically inferred dark matter density in the solar neighbourhood. Estimates of the latter vary in the literaturein fact, there is a long history of these types of measurements [69]-but most recent analyses using Gaia data place the local density of dark matter in the range ρ DM = [0.2, 0.7] GeV cm −3 [71]. The most important type of measurement for our purposes are the local techniques which aim to solve the Jeans equation for tracer populations of stars in the solar neighbourhood. Then after subtracting the baryonic contribution to the gravitational potential, one can infer the (sub-dominant) contri-bution coming from dark matter. These also usually give values in the same range, but importantly rely upon stellar tracer populations that span at least a few ∼100 pc. This is far, far above the scale probed by experiments, which is less than a milliparsec. The results of the previous section suggest that the typical density sampled by an experiment at a given instant is around 10% of the average density in the box which we use as a proxy for the distribution of axions in a similar-sized volume within our galaxy. The convention in the axion community has been to adopt the value ρ DM = 0.45 GeV cm −3 and to assume 100% of dark matter is in the form of detectable axions in the ambient dark matter in the solar neighbourhood. In Fig. 6 we show the status of experimental limits on g aγ in the approximate mass window relevant to the post-inflationary scenario. The standard published limits under the dark matter density convention mentioned above are shown in lighter-shaded colours. We show the extent to which this sensitivity is reduced in the worst-case scenario in which the volume in our simulation represents the typical density inside a galaxy at z = 0. A ∼10% reduction in the typical density corresponds to the factor of ∼3 suppression in sensitivity. Interestingly-or perhaps worryingly-this is comparable to the ratios in the couplings of the two common axion model benchmarks, the KSVZ [11,12] and DFSZ [13,14] models: KSVZ/DFSZ ∼ 1.92/0.75 ∼ 2.56. Hence, we can make the statement that an experiment ruling out the DFSZ model over some mass range will still have ruled out KSVZ even when most of the axions are bound up in miniclusters.
We generally expect from the plateauing growth of the miniclusters at the latest redshifts that the voids should not lose much more density in axions. Because of this fact, as well as the fact that ambient density ought to be partially refilled by disrupted miniclusters, we may expect our estimate of the typical density to be a conservative lower bound. It could prove to be higher once the galaxy has formed in full, and as the miniclusters are disrupted with their contents virialised into the main host. However, previous studies have shown that something around 60-70% of miniclusters orbiting around the solar neighbourhood would remain intact [74,75].
This suppression in density is perhaps most impactful for resonance-based experiments like ADMX [118], CAPP [124], HAYSTAC [129], TASEH [137,138], GrA-Hal [127], RADES [134], ORGAN [130,131] and MAD-MAX [29], that rely on only short total integration times at any one frequency. The suppressed typical density reduces the available signal these experiments can observe, and hence should be taken seriously when considering how much of the axion parameter space in the postinflationary mass range they have really constrained.
In contrast, for experiments that do not rely on resonances, and instead operate in a broadband mode over a wide range of frequencies, our aim was to study the possible implications for (i) individual measurements at different times spanning a relatively long period, and (ii) continuous measurements for a long period of time. In the first context (left panel of Fig. 5.), the most relevant implication of this result is for multiple different experiments each testing the veracity of a putative signal at different points in time. Naïvely, we would expect that if a signal were real it would remain at the same strength at every given measurement. Here, we see that measurements can be expected to vary in strength over several-year periods. Furthermore, to reconstruct the value of the axion's coupling, one has to break the degeneracy g 2 aγ ρ a somehow, usually by assuming a value of ρ a . We have shown in this work that a 10% suppression compared to the measured local dark matter density is expected, although this suppression could vary between, say, 5% and 30% on a timescale of only a couple of years.
In the second context-i.e., the right panel of Fig. 5we computed not the range in density suppressions, but rather the standard deviation in densities along many different trajectories. As in this case, we take the density relative to the mean density at that specific location, it is more relevant for experiments conducting single continuous measurements. This is not a relevant question for many current haloscopes, however for upcoming broadband haloscopes, it is. The best example is BREAD [32] which does not attempt any resonant enhancement of the axion signal, but rather has a broadband acceptance of photons over a wide band of frequencies. To compensate for its lower sensitivity at any given mass, it instead looks for signals over much longer integration times. BREAD is additionally relevant as in its proposed guise of "Giga-BREAD" in which a coaxial horn antenna will be used as the detector, it will have sensitivity to axions squarely in the post-inflationary mass range. Projections published by the collaboration assume total data-taking times ranging from 10 to 1000 days. Taking the longer of those durations, we can inspect the right panel of Fig. 5 to see that a signal variation of ∼ 30% is expected within that time.

VI. CONCLUSIONS
In this article, we have presented the results from a set of lattice and N-body simulations for the postinflationary QCD axion. We performed for the first time a realistic simulation of the axion field from the epoch of cosmic strings before the collapse z ∼ 10 16 and the time when the field configuration collapses into halos and forms minivoids, z ∼ O(10 2 ). This was done by combining a three-stage simulation, with lattice methods to solve the axion relativistic field equation and the nonrelativistic Schrödinger-Poisson system, and N-body methods for the collisionless dark matter dynamics.
We have placed a particular emphasis on the axions that occupy minivoid regions, in the space between small miniclusters and larger minicluster halos. As the minivoids occupy the vast majority of the sub-pc-sized simulation volume, the same will be true of a typical place in our galaxy. Direct detection experiments on Earth therefore do not sample the value of the dark matter density of ∼ 0.4 GeV cm −3 -which is obtained through measurements of stellar dynamics on scales well above this-but rather a suppressed value, given by the typical axion density in the voids. The quantitative answer to how great this suppression is, is around 10% of that density. We have tested various differences in the technical configuration of the simulation, including the physics included in the evolution, as well as the initial conditions, and find that this number is relatively robust.
We also made an estimate of the variance in this density, arguing that the density of axions along typical O(1)-year-long trajectories through the Galaxy can be expected to vary quite considerably. While the 10% suppression relative to the global average on galactic scales does not vary much, for experiments like broadband haloscopes making single continuous measurements of a potential axion signal, the density should be expected to vary by several tens of percent.
We should finish by remarking that our results strictly only apply at the final redshift of our simulationstypically z f 100-which is of course well before the present day. Extrapolating our result to z = 0 is beyond the scope of this work; however, we have tried to make statements that are conservative in several ways.
Firstly, it should be said that the fractions of axions in miniclusters and in voids are both plateauing by the end of the simulation, implying that the evolution is approaching a somewhat steady state. Secondly, the next important process for the lives of the miniclusters will be tidal disruption. The case can be made that our estimation of the density suppression and its variation can both be taken as a conservative lower bound. Due to the tidal and stellar disruption of the miniclusters, the axion dark matter distribution could be described by numerous virialised streams of axions, leading to larger typical densities, but also larger variations. Disruption has been studied for applications for indirect signals, e.g., collisions between miniclusters and neutron stars [78]. However further study is needed to quantify how much our results might change after considering disruption in the case of direct detection. It may well be that sufficient disruption occurs for the haloscopes to reclaim much of their earlier assumed sensitivity. In this appendix, we discuss the technical details of our lattice and N-body simulations.

Lattice simulations
As in Refs. [53,58] we use the public code jaxions 8 to simulate axions around the QCD phase transition. In particular, we consider the Klein-Gordon evolution equation for the complex field φ, which has as its solutions radial and massive saxion modes, angular and massless axion modes, and strings. Equation (A1) is used to evolve φ from an initial redshift of z i ∼ 10 16 until the collapse of topological defects that happens around z ∼ 10 13 . From there on, only the angular degree of freedom needs to be considered, and the equation of motion reduces tö where a/f a = arg φ, and m a (T ) ∝ T −n is the temperature-dependent axion mass with power-law index n = 7 [53,87]. We solve Eqs. (A1) and (A2) using a finite difference method and the symplectic integrator RKN4 [149]. The spatial derivatives are calculated from two neighbouring points and accurate to O(∆ 4 x ), where ∆ x = L/N , while the time-step is dynamically adjusted to resolve the fastest modes in the box [53]. To solve Eq. (A1) we implement the so-called fat-string or Press-Ryden-Spergel (PRS) trick [150,151] that keeps the string width constant in comoving coordinates by rescaling the potential  coupling λ → λ/R 2 . A resolution of N 3 = 8192 3 grid points is required for simulations in boxes of side length L > 8L 1 , in order to avoid the unphysical destruction of the string-wall network. For smaller simulation volumes a resolution of N 3 = 4096 3 is adopted. After a few background field oscillations, the axion field a(x) enters the linear regime as the field values become smaller. An adiabatic approximation is applied at z ∼ 5 × 10 11 , well after the axitons are expected to disappear as the axion mass reaches its zero-temperature value.
At this point, we apply the nonrelativistic approximation and map the conformal axion field aτ , where τ denotes the conformal time dτ = dt/a, into the complex conformal field ψ as in Eq. (3). This is done by identifying the real and imaginary parts of ψ, i.e., Note that, although this identification effectively means that we double the degrees of freedom stored in the field, the system of equations, Eqs. (4) and (5), that describes the evolution of ψ is a first-order system. Thus, the computational cost required to solve them is the same as that needed to solve the second-order Eq. (A2). We solve the Schrödinger-Poisson system with a resolution of 1024 3 grid points and similar schemes as described above, i.e., a sympletic time integration and finite difference for the spatial derivatives, also implemented in jaxions. The effects of gravity in the Schrödinger-Poisson system enter via the Madelung transformation in the equation of motion for the phase S = argψ, while the equation of motion for the modulus |ψ| can be recast into the energy density continuity equation. The phase field S will then drive the axion density towards the minima of the potential Φ N , and itself starts to wrap the fundamental domain [0, 2π) according to the depth of the well. Thus, the numerical evaluation will lose its accuracy when 2π phase jumps occur near the discretisation scale. In other words, there is a maximum value that can be reached in the gradient of S (i.e., the velocity field), which limits the evolution of the Schrödinger-Poisson system. This maximum phase gradient can be related to the maximum difference in the values of the gravitational potential that can be read at each time step. The condition (A6) is reached at redshift z * ∼ 2 × 10 6 , after which we continue the evolution of the system using N-body simulations, since with our choice of m a the Jeans wavelength is not resolved anymore, λ J < ∆ x ; see Fig. 7.

a. A note on velocities
As long as the particles are in the single-streaming regime, we can assign a single velocity at each local point by solving the Euler equation and hence v(x) = ∇ arg(ψ) m a .
This task can be performed using the full complex field ψ by exploiting the relation where ψ r,i denote the real and imaginary parts, respectively. Typical velocities in the initial conditions are the order v ∼ 10 −9 in natural units.

N-body simulations
Given the large range of time evolution, an N-body run is costly in terms of computational resources and runtime. We therefore limit the resolution of our simulations in this work to 512 3 dark matter particles and leave higher resolution runs for a future publication. For the simulations with box sizes L = 8L 1 and L = 16L 1 , we use a TreePM method as implemented in GADGET-4, with a PM grid of 512 3 points and a numerical softening length of s 1.95 × 10 −5 pc/h 4 AU/h in comoving units. This value is ∼ 30 and ∼ 60 times smaller than the average particle separation in the two simulation volumes used in this study. We modified  Minivoid fraction, FIG. 9. Minivoid volume fractions fv for "mild" underdensities, defined as regions with ∼ 80% of the average value ρ over the simulation volume. At the final simulation time of the L = 0.6 pc run (i.e., z f = 99), these mildly underdense regions occupy ∼ 20% of the simulation volume.
GADGET-4 to include the effects of radiation that dominates the energy budget at the initial time and perform our simulations using the following cosmological parameters: Ω m = 0.3, Ω Λ = 0.7, Ω r = 8.49 × 10 −5 , and h = 0.71. The average particle mass in our simulations is m av = 1.15 × 10 −16 M /h and m av = 4.4 × 10 −17 M /h, respectively, in our L = 0.2 pc and L = 0.4 pc simulations.
Appendix B: More on minivoids As described in the main text, minivoids are identified based on a pre-determined density threshold δ thr a on the axion energy density built from the particle positions. We  show in this appendix that the resulting void statistics do not depend on the specific choice of δ thr a . Figure 8 shows the minivoid volume fraction f v and energy density ρ v / ρ extracted from the L = 0.2 pc simulation using a range of threshold values −δ thr a = {0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.8}. Evidently, both f v and ρ v / ρ converge to the values given in Eq. (13) for all thresholds δ thr a ≥ −0.7. We find deviations for the choice of δ thr a = −0.8, as minivoids with ∼ 20% energy of the overall average value ρ occupy a sizeable fraction of the box. This behaviour is also observed in the L = 0.4, 0.6 pc simulations.
In addition, we want to stress that, while 80% of the volume is occupied by low-density minivoids, the remaining ∼ 20% of the volume is occupied by regions that are only mildly underdense-usually in the form of filaments connecting minicluster halos-and by regions of average density that limit the boundaries of miniclusters. We have checked this by imposing a lower density threshold in the void finder algorithm, i.e., 1 ≥ ρ v / ρ ≥ δ thr a .
(B1) Figure 9 shows the resulting f v from the L = 0.6 pc simulation for several choices of δ thr a . Finally, Fig. 10 shows the void size function at z = 999 for various choices of minivoid density thresholds for two different simulation volumes. Evidently, they are in good agreement with each other. Note that the void size function does not change substantially between z = 999 and z = 99 (see Fig. 4).