Transient Radio Signatures from Neutron Star Encounters with QCD Axion Miniclusters

The QCD axion is expected to form dense structures known as axion miniclusters if the Peccei-Quinn symmetry is broken after inflation. Miniclusters that have survived until today will interact with neutron stars (NSs) in the Milky Way to produce transient radio signals from axion-photon conversion in the NS magnetosphere. We quantify the properties of these encounters and find that they occur frequently ($\mathcal{O}(1-100)\,\mathrm{day}^{-1}$); last between a day and a few months; are spatially clustered towards the Galactic center; and can reach observable fluxes. These radio transients are within reach of current generation telescopes and therefore offer a promising pathway to discovering QCD axion dark matter.

Introduction -Peccei-Quinn (PQ) theory [1,2] predicts the existence of the QCD axion [3,4], which could simultaneously solve the strong-CP problem and act as a compelling candidate for particle dark matter (DM) [5][6][7][8]. The QCD axion is the pseudo-Nambu-Goldstone boson [9][10][11] of the new global PQ symmetry. Laboratory searches are underway worldwide to directly detect this QCD axion [12][13][14][15]. Astrophysical observations are also a promising avenue for detecting the axion [16][17][18][19][20]. In particular, radio observations could be used to look for emission from the conversion of axions into photons in neutron star (NS) magnetospheres [21][22][23][24][25][26][27]. In this Letter, we propose and characterize a new class of radio source, arising from encounters between NSs and overdense structures known as axion miniclusters (AMCs) [28][29][30][31][32][33]. 1 As axions fall towards an NS, they can resonantly convert into photons within the magnetosphere. This occurs at a radius R c at which the plasma frequency ω p in the magnetosphere equals the axion mass m a [23]. Assuming a Goldreich-Julian model for the NS magnetosphere [34], the power radiated per unit solid angle is derived in the WKB and stationary phase approximations as [21][22][23][24] where ρ a is the axion density at the conversion radius [23], g aγγ is the axion-photon coupling, and we have averaged over viewing angles. We also fix the rotation axis to be aligned with the NS dipole field and set the NS radius R NS = 10 km (see Supplemental Material for * Electronic address: thomas.edwards@fysik.su.se † Electronic address: kavanagh@ifca.unican.es ‡ Electronic address: luca.visinelli@sjtu.edu.cn § Electronic address: c.weniger@uva.nl 1 We use the terms 'miniclusters' and AMCs interchangeably. details of the NS modeling). The power scales with B 2 0where B 0 is the magnetic field strength at the NS poles -emphasizing why NSs are the most promising astrophysical target for these searches, with the highest known magnetic fields in the Universe O(10 10 −10 15 ) G [35]. The flux also scales with ρ a , meaning that regions of large axion density -such as AMCs -can give rise to very bright radio sources.
AMCs are a generic feature of models in which the PQ symmetry is broken after the end of inflation [32]. Large spatial variations of the axion density around the QCD epoch lead to the formation of minicluster 'seeds' [36] which collapse into gravitationally bound AMCs around matter-radiation equality [37]. This evolution has been confirmed by numerical simulations, which show that a significant fraction of DM axions might be contained within such bound structures [36,38]. Despite their low mass (10 −19 M < ∼ M AMC < ∼ 10 −5 M ) and large radius (10 −8 pc < ∼ R < ∼ 10 −2 pc), the density of an AMC can be many orders of magnitude larger than the local DM density [39]. On the other hand, AMCs are significantly more diffuse than stars. In our companion paper [40] (hereafter KEVW20), we show that tidal interactions with stars can have a dramatic effect on the survival of AMCs in the Milky Way (MW). AMCs towards the inner regions of the Galaxy undergo significant stripping and disruption, whereas those further out remain intact. This in turn is reflected in the observational signatures of AMCs.
Here, we build upon the results of KEVW20 to predict the rate, brightness, and sky distributions of encounters between AMCs and NSs. For a typical MW virial velocity of v ∼ 200 km/s ∼ 10 −11 pc/s, the time taken for a NS to pass through an AMC is expected to be O(10 3 − 10 9 ) s, meaning that these interactions would appear as radio transients. As we will show, this is particularly true for the brightest events which last between 10 5 to 10 7 seconds. We consider a Kim-Shifman-Vainshtein-Zakharov (KSVZ)-like QCD axion [41,42] of mass m a = 20 µeV, motivated by recent simulations [38,43]. This corresponds to a radio frequency of f = 4.84 GHz and an axion-photon coupling of g aγγ ≈ 8 × 10 −15 GeV −1 . All code associated with this work is available online at github.com/bradkav/axion-miniclusters [44].
Axion Miniclusters in the Milky Way -Due to the randomness of the initial overdensity fluctuations, miniclusters are born with a wide range of masses and densities. During matter-domination, the minicluster halo mass function (HMF) evolves under hierarchical structure formation, allowing ever heavier AMCs to form. N -body simulations modeling AMC evolution from recombination to z ≈ 99 predict a featureless HMF with a characteristic slope dn/d log M ∼ M −0.7 [33], corroborating semi-analytic studies [45][46][47]. We consider AMC masses between 3.3 × 10 −19 ≤ M AMC /M ≤ 5.1 × 10 −5 [45,46] and assume that these AMCs make up 100% of the DM, tracing the NFW profile of the MW halo. Miniclusters are also characterized by their overdensity parameter δ, which depends upon the random initial conditions of the axion field and its gradient at the onset of axion oscillations [29][30][31]. We take the distribution of δ from recent simulations [38] and map it to the AMC characteristic density ρ AMC as in Ref. [31].
Tidal interactions between AMCs and their local environment can have a significant effect on AMC properties [40,[48][49][50][51]. Stellar encounters prove to be the most important and can easily lead to the total disruption of AMCs. In addition, many successive weak encounters can cause surviving AMCs to lose mass, as well as altering their internal density. In KEVW20, we present Monte Carlo simulations used to assess the effects of these stellar tidal interactions, starting from the initial distribution of masses and densities described above. These simulations allow us to describe the properties of AMCs across the MW today.
The internal density profiles of AMCs are not well understood. For example, Ref. [33] finds that AMCs with masses M AMC > ∼ 10 −13 M have approximately Navarro-Frenk-White (NFW) [52] profiles in their outer regions whereas lighter AMCs are expected to have Power-law (PL) profiles [37]. For each simulation we therefore treat the entire population of AMCs as having a single universal structure given by either an NFW or a PL density profile (see Supplemental Material for examples). For a fixed characteristic density ρ AMC , the mean internal densityρ of our assumed NFW profile is O(10 5 ) times lower than for the PL profile. This leads to quantitative differences in the survival probability and distributions of masses and radii. Using these two density profiles therefore allows us to generously bound the uncertainties coming from the internal AMC structure.
Another source of uncertainty is related to the formation of axion stars (ASs). These are non-relativistic compact objects, described by solutions to the Schrödinger-Poisson equation [29,53], which can potentially form in the centers of miniclusters [54][55][56]. ASs have an inverse relationship between their mass and radius, leading to a potentially problematic scenario for a low-mass AMC in which its radius is smaller than that of the AS in its center. To avoid this issue, we first follow the evolution of all AMCs, described initially by the HMF above, then apply a cut to remove these potentially problematic light AMCs. 2 The remaining miniclusters form our fiducial sample and are used throughout the rest of this work (see Supplemental Material and KEVW20 for further details).
Minicluster -Neutron Star Encounters -The AMC-NS encounter cross section is given by σ(u) = πR 2 1 + 2GM NS /(Ru 2 ) , where R is the minicluster radius and M NS = 1.4 M is the NS mass. This expression includes a gravitational focusing term [57, p. 627] that depends on the relative velocity of the encounter u. We assume that the velocity dispersion of the NS and minicluster populations is σ v (r) = GM encl (r)/r ≈ 200 km/s, depending on the enclosed MW mass M encl (r) within a galactocentric radius r. The encounter velocity then fol- The velocityweighted cross section can then be written: Interactions on the outskirts of large, diffuse AMCs dominate the encounter rate but do not produce a signif-icant increase in the axion density close to the NS. We therefore consider only interactions with impact parameters b < b cut (with b cut ≤ R), such that the peak overdensity during the encounter is at least 10% of the local DM density ρ DM (r). We parametrize this in terms of an effective cross section σu (r) which saturates at the standard cross section σu (r) for dense AMCs. The expected rate of encounters over the entire MW is then [58,59] where n NS (r) is the NS number density at position r and dn AMC (r)/dR is the differential number density of AMCs with radius R, computed in KEVW20. Taking into account the distribution of AMC properties, we find that dΓ/dM AMC ∼ 1/M AMC , rising more steeply than this at low AMC masses, where the gravitational focusing effect becomes important for R < ∼R . We assume a population of 10 9 NSs in the MW, with 60% formed in the bulge and 40% in the disk [60,61], of which 20% have become unbound due to natal kicks [61]. Explicit expressions for n NS (r) are given in the Supplemental Material.
In Fig. 1, we show the integrand dΓ/dr of Eq. (3). At the largest radii (r > ∼ 10 − 20 kpc), encounters are rare due to the falling number of both AMCs and NSs. Near the Galactic center, the densities of both AMCs and NSs instead rise rapidly. However, the encounter rate is suppressed by the low survival probability of AMCs in this dense environment, leading to a plateau. 3 We find that the encounter rate is larger for eccentric than for circular orbits; AMCs on eccentric orbits spend less time at small radii, leading to a larger survival probability. For NFW profiles, there is a comparable contribution from encounters with bulge NSs at small radii and encounters with disk NSs at larger radii (close to the Solar circle r ). Miniclusters with PL profiles are more dense and therefore smaller than those with NFW profiles, leading to an overall decrease in the encounter rate. However, these dense AMCs are also more resistant to disruption in the Galactic center, leading to a greater survival probability at small r. This compensates for their smaller size and means that for PL miniclusters most encounters occur with NSs in the bulge. Over the entire Galaxy, we expect encounter rates of Γ PL = 3.4 day −1 and Γ NFW = 186.7 day −1 . If we had neglected the stellar disruption of AMCs described in KEVW20, these encounter rates would be larger by a factor of 1.4 and 45.4 for PL and NFW profiles respectively.
Signal Estimation -For each choice of minicluster profile, we sample 10 7 encounters to calculate the expected distributions of fluxes, durations, and sky locations. We sample the galactocentric radius of the encounter according to the encounter rate dΓ/dr in Fig. 1. We draw the height of the encounter z cyl above the Galactic plane from the distribution of NSs along the z cyl -axis (assuming that the AMC distribution is spherically symmetric), and we draw the galactocentric azimuth angle uniformly between 0 and 2π. We sample the AMC radius R following dΓ/dR at fixed galactocentric radius and sample the AMC density, given R, from the distributions derived in KEVW20.
Considering the trajectories of individual axions close to the NS, the maximum impact parameter which still crosses the conversion radius R c is: AMCs have radii several orders of magnitude larger than this, so we can consider the NS as tracing the internal AMC density ρ int (R) during the encounter. Such a direct encounter is likely to completely disrupt the AMC. However, the relaxation time for the AMC is much longer than the encounter time and we therefore neglect the evolution of the minicluster during the NS transit.
For each encounter, we estimate the radio flux density: where s is the distance of the encounter from Earth. The signal bandwidth BW is typically set by the axion velocity dispersion far from the NS [23], leading to narrowband line emission. However, because of the small internal velocity dispersion of the AMCs ( < ∼ 1 km/s), this is unlikely to be the main source of the signal bandwidth. We therefore fix the bandwidth of the signal to a larger value, 1 kHz, representative of the resolution of current and planned radio telescopes [63,69,70]. Determining the full directional dependence of the radio emission is highly non-trivial, though there have been a number of recent developments dealing, for example, with non-radial infall of axions [26,71]. Here, we assume for simplicity that the emission is isotropic; Eq. (1) has been averaged over viewing angle and we have fixed R c to its angular average. If future studies show that the radio emission is instead concentrated in a fraction f beam of the unit sphere, then our results can be straightforwardly re-interpreted: the observed rate will be reduced by a factor ∼ f beam and the flux density increased by a corresponding factor ∼ f beam .
We estimate the mean flux density S of each encounter by averaging Eq. (5) over the duration of the encounter T enc = 2 √ R 2 − b 2 /u. Note that the peak flux (where we assume that the observation time equals the duration of the encounter) for the Very Large Array (VLA) [63], SKA1-mid [64,65], and SKA2 [24]. We assume AMCs have Power-law (left) or NFW (right) internal density profiles.
during the encounter is generally comparable to the mean flux. Since each encounter is independent, the number of expected encounters starting within a time step ∆t is Poisson distributed with mean Γ∆t. Our simulations can therefore be combined into a time series for the predicted signal by taking a uniform distribution of start times. In Fig. 2, we show the distributions of T enc and S derived from our full sample of AMC-NS encounters. We also plot a smaller sample of 10 4 individual encounters, colored by the mean internal density of the AMC.
For PL AMCs, the distribution of encounters peaks at flux densities between 10 −6 µJy and 10 2 µJy, with a typical duration of 1-100 days. The distribution also includes a number of bright events O(1 Jy) which should be detectable by current radio telescopes such as the Very Large Array (VLA) [72]. These brightest events come from encounters with dense AMCs ( > ∼ 10 5 M pc −3 ). AMCs with NFW internal density profiles have a density around 10 5 times smaller than their PL counterparts, making high-flux events rarer. However, this is partially compensated by the larger encounter rate between NSs and NFW AMCs. The rate of encounters above a flux of 1 mJy (a sensitivity which has been achieved in recent searches for NS radio emission [26,27,73]) is Γ PL (Ψ > 1 mJy) = 0.04/day and Γ NFW (Ψ > 1 mJy) = 0.007/day for PL and NFW miniclusters respectively. Given the rate and duration of the encounters, we expect at least one bright event in the sky at all times. We also note that the rate of the brightest events ( S > 1 Jy) is relatively insensitive to our assumptions on the AMC density profile, once stellar perturbations are taken into account (see Supplemental Material).
The sky distribution of AMC-NS encounters is shown in Fig. 3 for AMCs with NFW profiles. In this case, the encounters occur predominantly towards the Galactic center although there is also a population of events extending along the disk, to a longitude of | | < ∼ 60 • . This morphology reflects the two populations of NSs in the bulge and disk. In the case of PL AMCs (not shown), encounters are concentrated almost exclusively towards the Galactic center, with 68% of events lying within 7 • of the center. Considering only the brightest events, we find that the distributions become even more concentrated towards the Galactic center, as shown in the bottom panel of Fig. 3.
Discussion and Conclusion -In this Letter, we have characterized the radio signatures of axion-photon conversion from encounters between NSs in the Milky Way (MW) and a population of QCD axion miniclusters (see KEVW20 [40]). These signatures will appear as regular transient radio point sources ( Fig. 1) with timescales varying from days to over a year. Interestingly, these transients will be spatially clustered towards the Galactic center ( Fig. 3) with potentially observable fluxes (Fig. 2). This suggests that radio observations could be used to discover QCD axion DM in the near future.
Within the MW, there are a variety of sources of transient radio emission, especially from the Galactic center. A recent analysis of archival VLA data [74] found a number of potential transients within ∼ 10 −3 degrees of Sgr A * (see also Ref. [75]). Potential explanations of these radio transients are pulsar emission, radio flares from dwarf stars, and outflows from X-ray binaries, all of which emit a broad energy spectrum. AMC-NS encounters could contribute to a population of transients towards the Galactic center. However, our results predict a characteristic line-like emission which would need to be confirmed with dedicated search strategies.
The slope of the AMC halo mass function (HMF) is not well constrained [33,[45][46][47]. To test the dependence of our results on this slope, we re-ran the entire pipeline assuming dn/d log M ∼ M −0.5 , as obtained using the Press-Schechter formalism [45]. Flattening the HMF increases the mean AMC mass and therefore reduces the total number of AMCs in the MW. 4 Fortunately, this is partially counteracted by an increase in the encounter cross-section in Eq. (2) which scales as R 2 ∼ M 2/3 . Overall, we find that the encounter rate has only a mild dependence on the slope; flattening from −0.7 to −0.5 leads to a factor of 5-10 fewer events above 1 mJy. A recent study [76] found an HMF slope that broadly agrees with the results of Ref. [33], but with an overall shift to lower masses. This would primarily result in a decrease to the number of AMCs passing the AS cut and therefore a reduced encounter rate.
The typical AMC density ρ AMC ∼ δ 4 is more strongly affected by the uncertainty in the AMC internal density profile (ρ PL /ρ NFW ∼ 10 5 ) than by a change in δ. We therefore do not expect that small variations in the distribution of δ would affect the detectability of the signal.
The production of axion miniclusters in the early Universe is a robust prediction of the post-inflationary scenario of axion cosmology [28][29][30][31][32]. The fraction of axions bound in these structures remains unclear [33], but is likely to be substantial. As we show in KEVW20, if this fraction is large then direct detection efforts may be ineffective. Our results are therefore complementary to these ongoing direct searches, alongside searches for continuous radio emission from the smooth halo of axions interacting with NSs [21][22][23][24].
Although we have calculated the population-level distribution of signals, much work is still needed to characterize the details of each event. More concretely, the precise signal bandwidth [77] and the modelling of the conversion process in realistic NS magnetospheres remain unclear. Both of these can have dramatic effects on the properties of the final signal and should be addressed in future work. We therefore emphasize that a non-detection cannot be reliably used to set upper limits on the axion parameter space. Nevertheless, this paper characterizes the unique transient nature of these interactions, and shows that current and near future radio telescopes have the sensitivity required to detect QCD axion DM. Finally, we acknowledge the use of the Python scientific computing packages NumPy [78] and SciPy [79], as well as the graphics environment Matplotlib [80]. Considerable effort has been put into modelling the population of neutron stars (NSs) in the Milky Way (MW) [66,67] given the sample of those we can actually observe (see for example those reported by the Australia Telescope National Facility pulsar catalogue [81]). We assume that the spatial distribution of millisecond pulsars in the MW can be used to approximate the corresponding distribution of old NSs, as in Ref. [24].
The MW hosts around 10 9 NSs [61], of which 20% have been unbound due to natal kicks [61]. Of these NSs, 60% are formed in the bulge and 40% in the disk [60,61]. We normalize the spatial distributions in the bulge and in the disk assuming the total numbers N bulge = 4.8 × 10 8 and N disk = 3.2 × 10 8 , respectively. We model the NS spatial distributions in terms of the galactocentric cylindrical coordinates r cyl and z cyl , which describe the radial distance from the axis of symmetry and the height from the Galactic plane respectively. Here, we assume that the spatial distribution of NSs in the bulge tracks the stellar population. We fix this in the companion paper [40] as a truncated Power-law distribution [82,83] n bulge (r cyl , z cyl ) = N bulge 11.1 kpc 3 where we use the parameters from Ref. [84], namely the core density ρ bulge 0 ≈ 99.3 M /pc 3 , r = r 2 cyl + (z cyl /q) 2 with q = 0.5, the bulge cutoff r 0 = 0.075 kpc, the exponent λ = 1.8, and r cut = 2.1 kpc. The numerical factor 11.1 accounts for the integration of the NS density over the bulge volume. Note, that this choice differs from other literature on the subject in which a Hernquist profile is assumed [24].
We use a Lorimer profile to model the distribution of millisecond pulsars in the Galactic disk [68] n disk (r cyl , z cyl ) = N disk C B+2 e −C 4π r 2 σ z Γ(B + 2) with parameters that are obtained from a fit to the population of almost one hundred millisecond pulsars -these are taken from Table III of Ref. [85], namely B = 3.91, C = 7.54, and σ z = 0.76 kpc. We have not incorporated any decay mechanisms for the NS's magnetic field, such as ohmic dissipation [86], ambipolar diffusion [87,88], or Hall drift [89]. We assume that all NSs have a mass of M NS = 1.4 M and radius R NS = 10 km.

B. Neutron Star Magnetosphere
Here, we use the Goldreich-Julian model [34] of the NS magnetosphere, for which the magnetic field along the axis of rotationη is where the radial dependence shows the typical dipole behavior falling as ∝ r −3 . For simplicity, we have assumed that the magnetic field is aligned with the axis of rotation, which are both inclined at an angle θ obs ∈ [−π/2, π/2] with respect to the observer. Each NS in the population is described by a magnetic field strength at the poles B 0 and a period P which are drawn from log-normal distributions, with mean and dispersion given by log 10 (B/G) = 12.65 and σ B = 0.55 for the magnetic field strength [66,67], and log 10 (P/ms) = 2.7 and σ P = 0.34 for the period [68]. Given the angular velocity vector of the NS Ω with absolute value Ω = 2π/P , the charged plasma in the magnetosphere at distance r has a number density [34] n c = 2ΩBη(r, θ obs ) e + relativistic corrections . The plasma frequency can be expressed as ω p = 4πα EM n c /m c where α EM is the fine structure constant and m c is the charge carrier mass. For electrons, we obtain The conversion radius R c is defined as the region for which the plasma frequency equals the axion mass. Using Eq. (S5), which is valid in the electron-dominated region, the conversion radius is given by [23] R c (θ obs ) = 224 km R NS 10 km 3 cos 2 θ obs − 1 B 0 where the resonant conversion only takes place if R c (θ obs ) > R NS .

C. Axion Minicluster Density Profiles
As described in main text, we use two different parameterizations for the internal density profiles of the AMCs. Since we do not know the internal density profiles precisely, these two choices are made to reflect the range of potentially observable radio signatures. An example of both density profiles and their corresponding truncation radii can be seen in Fig. S1.
An AMC with a Power-law (PL) profile is described by [45,90] where Θ (x) is the Heaviside step function. We truncate the PL profile at a radius where we fix ρ s r 9/4 s = ρ(δ)(R PL AMC ) 9/4 /4 [45], to give mean density ρ(δ) and the correct total mass for the AMC. On the other hand, AMCs with NFW density profiles are described by where the function f NFW (c) = ln(1 + c) − c/(1 + c) is defined in terms of a concentration parameter c ≈ 100 [33,47]. The truncation radius is now given by R NFW AMC = c r s .

D. Flux Distributions
Here, we give more details concerning the expected distributions of radio fluxes from AMC-NS encounters. In the left panel of Fig. S2, we plot the cumulative probability distribution of the mean flux density S (that is, the fraction of events above a given value of S ). We show results for AMCs with Power-law (solid blue) and NFW (solid olive) internal density profiles. The typical flux density from an encounter between an NS and a Power-law minicluster is larger because these AMCs are substantially more dense than those with NFW profiles. We show also the results for AMCs which have not undergone perturbations due to stellar encounters (dashed lines). For Power-law miniclusters, these results are very similar to the perturbed case; their higher density also makes them more resistant to disruption. Instead, for NFW miniclusters, the typical flux which we would expect when neglecting perturbations is much smaller than when perturbations are included.
We can see this expressed also in terms of the encounter rate above a given threshold in flux S , as shown in the right panel of Fig. S2. In the NFW case, going from the perturbed to unperturbed distributions, the encounter rate drops by a factor of around 40. However, the rate of very bright encounters actually increases once perturbations are taken into account. As we show in detail in Ref. [40], the survival probability for NFW miniclusters is typically larger than 50% throughout the MW. However, surviving AMCs are stripped of a significant fraction of their mass, typically leaving behind a much more dense remnant AMC. Thus, what would be common encounters with large, diffuse AMCs in the unperturbed case become rarer but brighter encounters with small, dense AMCs once perturbations are accounted for.
Of particular interest is that for the very brightest events (above around 1 Jy) the NS encounter rates for Power-law and NFW miniclusters start to converge, typically to within an order of magnitude. Despite the substantial differences in their sizes and densities, we find that the rate of bright encounters between NSs and AMCs in the MW is somewhat insensitive to the initial density profiles of the AMCs.

E. The Role of Axion Stars
An axion star (AS) [91] is a condensate made of cold axions, described by a solitonic solution of the relativistic Klein-Gordon equation [92,93]. The axions inside the star are usually non-relativistic, so that a description in terms of the Schrödinger-Poisson (SP) equation is often a suitable approximation. Since axions are pseudo-scalar particles, ASs differ from the analogous solutions for 'boson stars' obtained in scalar boson theories [94]. While for scalar bosons a static solitonic solution to the SP exists, for pseudo-scalar axions the solution has to be oscillating periodically in time. An AS is then made up of a self-gravitating, oscillating axion field.
ASs may form in the dense central region of an axion minicluster, where the density is high enough that two-to-two processes enable the cooling of its inner core and lead to the formation of the condensate [29,53]. This process has been observed in recent numerical simulations [54][55][56]. The existence of ASs could be indirectly probed through their