{'fraction_non_alphanumeric': 0.03879049950716988, 'fraction_numerical': 0.02165272964293663, 'mean_word_length': 4.239733444398167, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Radio signatures from encounters between Neutron Stars and QCD-Axion Minihalos around Primordial Black Holes

Probing the QCD axion dark matter (DM) hypothesis is extremely challenging as the axion interacts very weakly with Standard Model particles. We propose a new avenue to test the QCD axion DM via transient radio signatures coming from encounters between neutron stars (NSs) and axion minihalos around primordial black holes (PBHs). We consider a general QCD axion scenario in which the PQ symmetry breaking occurs before (or during) inflation coexisting with a small fraction of DM in the form of PBHs. The PBHs will unavoidably acquire around them axion minihalos with the typical length scale of parsecs. The axion density in the minihalos may be much higher than the local DM density, and the presence of these compact objects in the Milky Way today provides a novel chance for testing the axion DM hypothesis. We study the evolution of the minihalo mass distribution in the Galaxy accounting for tidal forces and estimate the encounter rate between NSs and the dressed PBHs. We find that the encounters give rise to transient line-like emission of radio frequency photons produced by the resonant axion-photon conversion in the NS magnetosphere and the characteristic signal could be detectable with the sensitivity of current and prospective radio telescopes. It would be important to investigate in detail search strategies for such signals which would provide a novel pathway for QCD axion detection.


Introduction
In this paper, we propose for the first time a new avenue for detecting the QCD axion via transient radio signatures coming from the encounters between neutron stars (NSs) and axion dark halos around primordial black holes (PBHs). Crossing of the NS magnetosphere through the axion minihalo leads to monochromatic radio signals due to axion-photon resonant conversion [1]. We show that the signal should be detectable on the Earth under suitable conditions.
Even though a wide range of observations coming from, for example, gravitational lensing [2], hot gas in clusters, the cosmic microwave background radiation [3,4] and the hierarchical structure formation of the Universe [5], are well explained by the inclusion of cold dark matter, its origin in the particle physics context remains unknown [6,7]. The strong CP problem [8][9][10] and unification with gravity in the context of string theory make the QCD axion [11] one of the strongest DM candidates. In last years, the quest for the axion dark matter has strongly attracted the attention of the physics community leading to a diverse search program including, for example, haloscopes [12][13][14][15], helioscopes [16] and indirect axion searches [17][18][19][20][21][22][23][24]. Interesting enough, most of these searches are based on the axion-photon coupling through the interacting Lagrangian ∆L ∼ g aγγ φ E · B, where φ is the axion field, g aγγ is the axion-photon coupling constant, and E and B are the electromagnetic field components, respectively.
The axion searches currently place very little constraints on the QCD axion DM and it is therefore extremely important to study other indirect ways to test this theoretically well motivated setup. In this paper, we show that the presence of a small abundance of primordial black holes, which form a subdominant co-dark matter component, enhances the probability for detecting the QCD axion via indirect searches. It should be noted that the LIGO-Virgo gravitational wave detection [25,26] and recent NANOGrav results [27][28][29] may indeed suggest the existence of a small fraction of DM in PBHs in the O(1 − 10)M mass range. PBHs are one of the oldest DM candidate [30][31][32][33] and they have been studied in several setups with other DM candidates such as the axion [34,35] and WIMPs [36][37][38].
We consider a general QCD axion model in which the PQ symmetry is broken before or during inflation. Our results may apply also in the case where the PQ symmetry is broken after inflation, provided that the domain wall number N DW = 1 so that the well known domain wall problem is avoided. When the Universe temperature drops to the QCD scale, the QCD instanton effect lifts the axion potential so that the axion acquires mass. Then, the axion field starts to roll down to one of the N DW degenerate minima so that domain walls are formed keeping away the different vacua. If N DW = 1, axionic strings formed during the PQ symmetry breaking attach to each domain wall forming disk-like objects which quickly collapse under their surface tension [39]. This phenomenon leads to an upper bound for the axion decay constant to avoid an overabundance of DM such that F a (4.6-7.2) × 10 10 GeV [40] . For N DW > 1, the string-wall network is long-lived and the QCD axion is excluded in the standard scenario 1 . If the PQ symmetry is broken before or during inflation, the axion field acquires quantum fluctuations proportional to the Hubble parameter during inflation leading to the so-called isocurvature perturbation problem [41,42]. However, several possible solutions to this issue have been already presented in the literature (e.g., see Ref. [43]). Since inflation expands a tiny region, where the axion field is homogeneous, to a region larger than the current Hubble horizon, the entire observable Universe holds the same value for the axion field. Thus, there are no production of domain walls as before and the axion abundance is given by the initial misalignment angle θ i . For θ i ∼ 1, the upper bound over the axion decay constant to avoid the overclosure of the Universe read as F a 10 12 GeV [44][45][46] For typical QCD axion models, we have the Kim-Shifman-Vainshtein-Zakharov (KSVZ) [51] and the Dine-Fischler-Srednicki-Zhitnitsky (DFSZ) [52,53] models. In the KSVZ model the number of degenarate minima in the QCD potential is the number of heavy quarks which carry U(1) PQ charge, so that N DW = 1 may be realized. In the DFSZ model, the domain wall number equals or is twice the number of flavours of quarks carrying U(1) PQ charge, namely N DW = 3 or 6. Apart from the phenomenological point of view, there is no general theoretical reason to choose one model over the other.
The QCD axion mass is related to the axion decay constant as [9] m a = f π m π √ m u m d where f π and m π are the pion decay constant and mass, respectively, and m u and m d are the up and down quark masses, respectively. Even though the axion interactions with Standard Model particles scale as ∼ f π /F a , there are several non-negligible model-dependent factors. For the case of our interest, the axion-photon coupling constant is given in terms of the electromagnetic (E) and color (N) anomalies of the axial current associated with the axion as follows [54] g aγγ = α em 2πF a where α em is the fine structure constant. The dependence of the axion-photon coupling on the axion decay constant and the ratio E/N makes g aγγ varies several orders of magnitude.
In the present work, we will consider as benchmark models the DFSZ and KSVZ models in which E/N is equal to 8/3 and 0, respectively, so that g DFSZ aγγ 9 × 10 −14 GeV −1 10 10 GeV F a , (1.3) |g KSVZ aγγ | 2 × 10 −13 GeV −1 10 10 GeV F a . (1.4) We will take the axion mass to be a free parameter but within the range defined by the upper bound over the axion decay constant of the classical QCD axion window and the largest axion mass at which the resonant axion-photon conversion takes place in the NS magnetosphere. In the parameter space of our interest, we will see this mass range is approximately equivalent to 10 10 GeV F a 10 12 GeV in terms of the axion decay constant. The setup we are studying, in which the axion is the dominant DM component and PBHs constitute a subdominant DM component, leads to the presence of PBHs surrounded by axion minihalos in the Milky Way today. These astrophysical objects are sometimes called dressed PBHs in the literature and we will also use this terminology here. According to the theory of spherical gravitational collapse [55], any overdensity in the DM distribution will seed the growth of a minihalo. Thus, the smooth axion DM background will be accumulated around PBHs mainly during the matter-dominated era, so that dressed PBHs will end up in galactic halos. Generally speaking, probing the axion DM hypothesis is notably hard since the QCD axion is predicted to interact very weakly with Standard Model particles. That is why the presence of dressed PBHs, where the axion density may reach several orders of magnitude larger than its local DM value, offers an unique opportunity for axion DM detection. The resonant axion-photon conversion in a cosmological setup was discussed for the first time by one of us in Ref. [56]. Two decades later, the same resonant phenomenon was studied in Ref. [1] in the neutron star environment as a novel way of testing the axion hypothesis. This work was extended later in Refs. [20,21], where the sensitivity of current and prospective telescopes with respect to the associated radio signals was analyzed in detail. A threedimensional computation of the photon flux, considering the axion phase-space distribution and the particular configuration of the NS magnetosphere, was performed in Ref. [57]. In addition, the axion-photon conversion has been also studied beyond the DM background scenario by considering axion DM substructures such as axion stars 3 [18,65,66], axion spikes around intermediate mass black holes [67], and axion miniclusters [24]. Note that the axion miniclusters [68][69][70][71] may be formed only if the PQ symmetry is broken after inflation and if the QCD axion potential has just one minimum, according to numerical simulations performed in Ref. [40] as explained before. Our interest is to analyze the photon emission from axion resonant conversion by encounters between NSs and dressed PBHs in the Galactic halo. Considering a general QCD scenario in which the PQ symmetry is broken before (or even during) inflation, if a small fraction of DM in PBHs exists, as suggested by LIGO-Virgo observations [25,26] and NANOGrav results [27][28][29], then the formation of dressed PBHs will take place. Since their typical length scale is the order of parsecs and the axion density within their minihalos reaches values much higher than the local DM density, these astrophysical compact objects provide a unique chance for testing once and for all the axion DM hypothesis.
The present paper is organized as follows. In Sec. 2, we briefly review the DM minihalo formation around PBHs when they constitute a small initial fraction of the total DM content. We make special emphasis on the DM radial profile of these minihalos and propose a reasonable initial mass distribution for dressed PBHs immediately after being incorporated to galactic halos. In Sec. 3, we estimate the effects of tidal forces existing in the Milky Way Galaxy acting on dressed PBHs. In Sec. 4, we estimate the encounter rate between NSs and dressed PBHs in terms of their location in the galaxy. In Sec. 5, we estimate the transient radio signal during a NS-dressed PBH encounter due to the axion-photon resonant conversion in the NS magnetosphere. We also estimate the sensitivity of current and prospective radio telescopes for the QCD axion detection. Finally, in Sec. 6, we discuss and summarize the main ideas of the paper.

Axion Dark Minihalos around PBHs
Primordial black holes which are formed with masses larger than ∼ 10 15 gram do not evaporate but begin to acquire a dark matter halo by gravitationally capturing the surrounding axion DM from the smooth background.
Since PBHs are local overdensities in the axion DM distribution, they inevitable seed the growth of spherically symmetric minihalos. The spherical gravitational collapse (or secondary infall) theory tells us that at about the time of radiation-matter equality, t ∼ t eq , self-similar dark halos begins to grow around PBHs, during the matter dominated era proportional to the cosmological time as t 2/3 [55]. Assuming the absence of tidal forces over PBHs, negligible peculiar velocities and an initial dark matter background in the Hubble flow, the virialized minihalo mass and radius scale with the redshift as [72,73] M halo (z) = 3 1000 where the minihalo radius is about one third of the turnaround radius or radial distance at which the DM shell was able to break free from the Hubble flow. Here we need to mention that Eq. (2.1) should be considered as an optimistic upper bound to the mass growth of minihalos. Indeed, the assumption of isolated and stationary PBHs made here begins to break down around z ∼ (30 − 10), e.g. at the time of first galaxies formation. The density profile of the internal structure of these minihalos can be readily derived from Eqs. ( which is valid for distances at which M halo (r) > M PBH . This profile agrees up to a numerical factor of two with the more recent work performed in Ref. [74]. The ρ halo ∼ r −9/4 density profile was confirmed by the first N-body simulations of a universe containing a subdominant fraction of PBHs in a background of DM particles performed in Ref. [37]. One crucial question to be answered is the total final mass of the induced minihalo around the central PBH. Taking into account the presence of the cosmological constant, there is a last bound shell, namely, a last shell of DM which may be accreted by the black hole. All shells with a radius larger than that particular shell will unavoidably continue expanding in the Hubble flow. Based on that, the maximum minihalo mass is estimated to be M halo ≈ 1500 M PBH in Ref. [72] by calculating the radial infall of a DM shell under the effect of the cosmological constant. However, the minihalo mass probably will not reach this final mass since dressed PBHs begin to interact among them and/or run out of matter to accrete as minihalos grow. This scenario breaks the assumptions of isolated and stationary PBHs embedded in a homogeneous background made in the secondary infall theory. The final induced minihalo mass is calculated in Ref. [74] by assuming that the induced minihalo stops the accretion process in the nonlinear regime when density perturbations around dressed PBHs are the order of the growing minihalo mass. Under this consideration, the final average mass for minihalos is reported to be about M halo ∼ (10 1.5 − 10 2.5 )M PBH for the range M PBH ∼ (10 −8 − 10 2 ) M (see Eq. (35) in Ref. [74]). At around z ∼ 6, dressed PBHs will be incorporated to galactic halos and the massradius relation of the minihalos is readily derived from Eqs. (2.1) and (2.2) as We do not expect a monochromatic mass function for the minihalo mass in terms of the central PBH mass but certain spread of masses around an average value. We model the initial mass function (IMF) before any disruption effects by a Gaussian distribution with an average minihalo mass M halo,0 and variance σ halo,0 . We assume the distribution is constrained between minimum and maximum minihalo mass values M min halo,0 and M max halo,0 so that (2.5) The normalization constant C halo,0 is given by (2.6) As a benchmark model, we take M halo,0 = 150M PBH , M min halo,0 = 30M PBH , M max halo,0 = 270M PBH with σ halo,0 = 24M PBH such that M min halo,0 ≤ M halo,0 ≤ M max halo,0 correspond to ±5σ halo,0 . In the expressions above, we have introduced the symbol 0 to make it explicit that we refer to quantities before disruption effects. Later in the text we will suspend this for brevity.
As galaxies evolve, further dressed PBHs and axion DM background are incorporated to galactic halos. After that dressed PBHs will undergo different levels of disruption depending on their orbital radii around the galactic center. Some of these dressed PBHs will undergo total disruption and some of them will undergo a partial mass loss so that the initial mass distribution given in Eq. (2.5) will be modified. We discuss the mass function after disruption in Sec. 3.4.

Disruption of Dressed PBHs in the Milky Way
After dressed PBHs are incorporated to the galactic halos at around z ∼ 6 [75], their dark minihalos will undergo different levels of disruption depending on their minihalo masses, central PBH mass and radial orbit in the Galaxy. Generally speaking, main sources of disruption for dark matter substructures in the Milky Way are: (1) global tides coming from the mean-field potential of the galaxy, (2) high speed encounters with stars and (3) tidal shocking during disk crossing.
Large dark matter substructures with masses about 10 7 M − 10 9 M within a Galactocentric radius of ∼ 30 kpc would undergo a depletion of their abundance by a factor of (2-3) due to disk shocking [76]. Axion miniclusters would undergo a tiny depletion of their abundances at the solar neighborhood, being the main source of disruption high speed encounters with stars. Taking a typical axion minicluster mass and radius of ∼ 10 −12 M and ∼ 10 7 km, respectively, and assuming a power-law density profile (∼ r −1.8 ), a depletion of about (2−5)% was reported in Refs. [77,78]. Recently, a Montecarlo simulation was performed in Ref. [79], where an axion minicluster mass range of ∼ (10 −19 − 10 −5 ) M is assumed. At the solar position, they reported survival probabilities of 99% and 46% for power-law (∼ r −9/4 ) and Navarro-Frenk-White [80] density profiles, respectively.
For the case of our interest, the cuspy density profile shown by dark minihalos around PBHs would offer them certain level of protection against disruption in the Milky Way. Analytical and numerical estimates performed by some of us in Ref. [81] show that typical dressed PBHs with M halo 100 M PBH lose about 9% of their minihalo masses, when circular orbits at the local neighborhood are assumed. Here, we extend our previous analysis performed in Ref. [81] to estimate the disruption map of dressed PBHs in the Milky Way. We will focus on disruption coming from the mean field potential of the Milky Way and disk shocking when dressed PBHs cross the galactic disk. Disruption coming from high speed encounters with stars is disfavored for the PBH mass regime that interests us the most. As was discussed in Sec. 2.2 of Ref. [81], the critical impact parameter b c at which the gained internal energy of the minihalo after the encounter, ∆E(b), is equal to the binding energy of the dressed PBH, where V rel is the relative velocity between both astrophysical objects and M is the star mass.
Take as a Benchmark model M = M and M PBH = 30M . Only tiny impact parameters in terms of the minihalo radius are able to immediately disrupt the minihalo or lead to a significant minihalo mass loss. This kind of encounter is highly disfavored so that it is very unlikely that minihalos undergo one-off disruption. Thus, the total probability for disruption, N total , is mostly given by multiple star-dressed PBH encounters holding impact parameters such that b > b c . The needed time for total disruption is longer than the Milky Way age in the parameter space of our interest even for large stellar densities, n . The probability for disruption from multiple encounters during a time t is estimated as [81] Suppose that dressed PBHs undergo circular orbits close to the Galactic center at a radius of 100 pc. Thus, the needed time for minihalo disruption (N total = 1) is estimated to be t 200 Gyr (21 pc −3 /n ) for minihalos with masses 30 M PBH M halo 270 M PBH . We approximate the total stellar density at 100 pc by the bulge contribution using the spherical version of the bulge profile described in Ref. [82] (see below Eq. (3.6) and Table I in Sec. 3.1).

Milky Way Mass Model
There is a rich literature with respect to mass models of the Milky Way. The first example dates from the sixties in Ref. [83]. Later, further models were put forward in Refs. [84][85][86]. Here, we use the Milky Way model reported in Ref. [82] (M16), which complies with theoretical modelling and fits the observational constraints. This mass model is composed of a spherically symmetric DM halo profile, a two-component (thin/thick) galactic disk in cylindrical coordinates, a two-component (H I and molecular) gas disk and a bulge profile.
The M16 model assumes a Navarro, Frenk, and White (NFW) profile [87] for the DM halo, which is parameterized as where ρ s and r s are the scale density and the scale radius, respectively, and r is the Galactocentric radius. This kind of profile is a common approximation to the density profile obtained in dark-matter-only cosmological simulations. The two-component stellar disk is modelled in cylindrical coordinates (z, r cyl ) with an exponential decay as we departure from the Galactic plane. The density of both disk components is modelled as where for each component z d and r cyl,d are the scale height and scale length, respectively, and Σ 0 is the central surface density such that the total disk mass reads as M d = 2πΣ 0 r 2 cyl,d . In addition to the stellar disks, the M16 model includes a two-component gas disk. This inclusion is necessary to make the number of stars with high vertical velocity in the local neighborhood dynamically consistent with the observed stars far from the Galactic plane. The gas disk density profile reads as where parameters Σ 0 , z d , r cyl, d are similar to those from the stellar disk case. The parameter r cyl, m is the scale length associated with the central hole of the configuration. The density of the gas disks falls off as exp(−z/z d ) as z → ∞ in the same way as those for the stellar cases.
Lastly, the M16 model uses an axisymmetric bulge profile so that in cylindrical coordinates we have r = r 2 cyl + (z/q) 2 , where q = 0.5 is the axial ratio. We take for simplicity the spherical version of this profile so that the density profile reads (3.6) Table 1: Best-fit model for the Milky Way Galaxy performed in Ref. [82] including the baryonic component parameters. Values of the parameters, including those associated with the baryonic components, for the best-fit model in Ref. [82] (M16) are listed in Table I. For the bulge profile, we report the parameter values (r 0 , α b , r b , ρ 0,b ) associated with the axysymmetric model but we approximate q ≈ 1.

Global tides from the Milky Way
Unless dressed PBHs are orbiting close to the Galactic Center, the high density of dark minihalos around PBHs offers them protection against the tidal stripping. The overall Galactic potential acts on dressed PBHs stripping the region of the minihalo which is beyond the tidal radius, r tidal . The stripped fraction of particles within the tidal radius is only a second-order correction corresponding to the highspeed particle population close to the escape speed [88].
Applying the distant-tide approximation, the tidal radius is calculated as [89] where R is the radial distance of the dressed PBH from the Galactic center under a circular orbit (or the perigalactic distance under eccentric orbits [90]) and M MW (R) is the mass of the Milky Way enclosed within the radius R using the M16 model. The enclosed mass depend on the Galactic mass profile composed of the DM halo, the stellar and gas disks and the bulge. Using M halo (R halo ) M PBH and the spherically symmetric mass profile of the minihalo, we have [81] which depends on the minihalo mass in units of the central PBH mass rather than the M PBH itself. Figure 1 shows a contour plot for the ratio r tidal /R halo in the parameter space (M halo /M PBH , R). The lighter the minihalo mass in units of the central PBH mass, the stronger its resistance against tidal stripping. For M halo = 100 M PBH ( M halo = 30 M PBH ), we have r tidal /R halo 1 when dressed PBHs orbit around the Galactic Center at a closest radius R 2.5 kpc (R 0.1 kpc).

Gravitational Field of the Galactic disk
As dressed PBHs cross the Galactic plane, they undergo a compressive force coming from the gravitational field of the disk. This force pinches them along the normal to the disk plane so that the cumulative effect of successive disk crossings may eventually disrupt minihalos. This source of tidal disruption known as disk shocking was first proposed in Ref. [91] in the context of globular clusters.
Modelling the Galactic disk as an infinite slab when dressed PBHs are just about to cross the disk plane and under the approximation of circular orbits for axion particles in minihalos, the gained energy of a minihalo per unit of axion mass in a single crossing reads as (see Sec. III.D.2 in Ref. [88] and Sec. 2.3 in Ref. [81]) where ρ d (r cyl , 0) is the density profile of the thin disk at the Galactic plane and A(η) = (1 + η 2 ) −3/2 is the adiabatic correction which depends on the adiabatic parameter η. The adiabatic parameter is defined as η(r, R) = ω(r)τ cross (R), where ω(r) is the orbital frequency of the axion particles at distance r from the central PBH and τ cross (R) is the effective crossing time of dressed PBHs under circular orbits with a Galactocentric radius R. The impulse approximation holds when the axion orbital time within the minihalo is much longer than the disk crossing time, namely when η(r) → 0 so that A(η) → 1. By contrast, when η(r) > 1 the efficiency of the disruption due to disk shocking starts to damp out. The angular momentum conservation associated with an axion particle orbiting several times during a single disk cross acts as a protection against disruption. We may estimate the axion orbital frequency by taking the inner dispersion velocity at radius r as ω(r) = v DM (r) 2 /r. Using the isothermal approximation where the square of each Cartesian component of the velocity dispersion is equal to the half of the square of the circular velocity, we have The effective crossing time is estimated as function of the half-height H of the disk and the vertical component of the dressed PBH velocity at a radius R in the Galactic frame, V z (R). We approximately have (3.11) Figure 2 shows the contour plot for the adiabatic correction A(η) considering a flat Galactic velocity curve of ∼ 220 km/s in the parameter space of (r/R halo , M halo /M PBH for r/R halo ∼ 1. Thus, the lighter the dressed PBH the more resistant they are against disk shocking. Note that the adiabatic correction depends on the fraction M halo /M PBH rather than the central PBH mass as shown in Eq. (3.10), when the axion orbital radius is expressed in terms of the minihalo radius. We assume dressed PBHs undergo circular orbits around the Galactic center during the age of the Milky Way, T MW ∼ 10 Gy, so that the total number of disk crossing is calculated as where M MW (R) is the total Galactic mass within a radius R using the M16 model. Before the first disk crossing, we take the minihalo radius to be R halo (N cross = 0) = min(r tidal , R halo ), where r tidal is the tidal radius considering global tides from the Milky Way as calculated in Sec. 3.2. The associated minihalo mass M halo (N cross = 0) is read from Eq. (2.4).
Define the gravitational potential of the minihalo at radius r as so that after the first disk crossing we may recalculate the new tidal radius by comparing the gained energy of the minihalo per unit of axion mass, Eq. (3.9), with the change in the gravitational potential as [81,88] where r is the minihalo radius before the disk crossing. The above equation is solved iteratively until the total number of disk crossings, N cross (R). The final tidal radius is related with the final tidal mass of the minihalo via Eq. (2.4). Our analysis assume that between successive disk crossings dressed PBHs are able to virialize after a partial loss of particles. The virialization time is given by the free-fall time, t free-fall = 3π/(32G N ρ halo ) where ρ halo ∼ 3M halo /(4πR 3 halo ). Since we are considering circular orbits around the Galactic Center with a Galactocentric radius R, we require that the half of the orbital period is at least the order of the free-fall time, so that For dressed PBHs with minihalo masses ranging as 30 M PBH M halo 60M PBH , this condition is satisfied until radii very close to the Galactic center, R 100 pc. However, as the minihalo mass increases the virialization assumption begins to break done for orbits with larger radii introducing a certain degree of inaccuracy in our results. For example, for a typical minihalo mass of M PBH ∼ 100 M PBH , we require R 450 pc. Here we point out that the realistic picture includes eccentric orbits with much larger orbital periods relaxing Eq. (3.15). Figure 3 shows the contour plot for the survival mass fraction of minihalos in the Milky Way in the parameter space (M halo /M PBH , R), where we are assuming dressed PBHs undergo circular orbits of radius R around the Galactic center. The survival mass fraction f surv is calculated by taking the ratio between the final minihalo mass after N cross crossings and the initial minihalo mass before undergoing global tides and disk shocking disruption. In the whole range of minihalo masses studied by us, the mass loss is negligible when dressed PBHs undergo circular orbits around the Galactic center at a radius R 13 kpc. As we approach the Galactic center, the mass loss begins to increase and the survival mass fraction decreases down to f surv ∼ (0.05 − 0.1) for 135 M halo /M PBH 270 at R ∼ 100 pc. We see that the heavier the minihalo in units of the central PBH mass or the closer the minihalo orbital radius, the larger is the mass loss due to disruption in the Milky Way.

Minihalo Mass Function after Disruption
Since dressed PBHs undergo loss of mass due to tidal forces and disk shocking in the Milky Way Galaxy, the disruption effects cause the mass function to deviate from its initial form Eq. (2.5). Figure 3 shows that there is a sizeable loss of mass from minihalos as the galactocentric radius decreases. This loss of mass will produce a redistribution of the spectrum of masses for minihalos according to their orbital radius. If minihalos are not totally destroyed, which holds for the galactocentric radii considered in Fig. 3, the total number of minihalos is conserved so that we need to re-normalize the mass distribution after considering the effects of tidal forces from Milky Way. Figure (4) shows the mass function after disruption. For a galactocentric radius R 13 kpc the disruption effects are negligible and the mass function practically coincides with the IMF (solid black line). As the galactocentric radius becomes smaller, the mass function begins to narrow and move towards lighter masses. At the solar position (red dashed line) and at 6.5 kpc (blue dashed line), most of minihalos lose about Figure 4: Mass function dP/dM halo at different galactocentric radii assuming circular orbits for dressed PBHs and the initial mass function (before disruption) given in Eq (2.5).
∼ 11% and ∼ 33% of their initial masses, respectively. The mass spectrum redistribution leads to a significant deviation from the IMF which will further impact the encounter rate between NSs and dressed PBHs.

Number Density for Neutron Stars
There have been several efforts to model the NS distribution in the Milky Way (see e.g. Refs. [92][93][94]). Neutron stars are mostly born during the core-collapse of massive stars. We may obtain a rough estimate of the total number of NSs produced in the Milky Way through the present-day core-collapse supernovae rate [95]. The total number of NSs born is about 10 8 − 10 9 , so that they represent a sizeable fraction of the Galactic stellar content. By performing Monte Carlo simulations to model NS orbits under suitable assumptions with respect to the Galactic potential, distribution of progenitors and birth velocities, Ref. [96] reported that about 80 percent of NSs are in bound orbits and the remaining ones are unbound due to natal kicks. Following Ref. [97], where the space and velocity distributions of Galactic isolated old NSs is addressed, we estimate that 40% and 60% of the total NS population were formed in the Galactic disk and bulge, respectively. In particular, we take NS disk = 3.2 × 10 8 and NS bulge = 4.8 × 10 8 , so that NS total = NS disk + NS bulge = 8 × 10 8 .
Here the interesting fact that the Galactic bulge contains less stellar mass than the Galactic disk but more NSs are produced is because the bulge stellar mass function is skewed towards heavier stars than that in the disk [97].
Due to the lack of pulsars observed in the inner regions of the Milky Way, we do not have a direct measure of the isolated NSs which are distributed through the Galactic bulge. While Ref. [22] considered a Hernquist profile for that distribution, Ref. [24] assumed that this distribution tracks the stellar population in the bulge. We will assume that NSs tracks the stellar bulge so that the NS number density in the bulge n NS, b (r) is estimated as where M NS = 1.4 M is the characteristic mass for NSs having a characteristic radius of R NS = 10 km and the spherically symmetric stellar bulge density ρ b (r) is given by Eq. (3.6).
The prefactor in Eq. (4.1) comes from the integration of the bulge stellar density to ensure a total number of NSs in the bulge equal to NS bulge . For the case of the NS distribution in the Galactic disk, we model the NS number density n NS,d as a double exponential profile as follows where r NS = 5 kpc and z NS = 1 kpc as was reported in Ref. [98] where the distribution of pulsars in the Galactic disk is studied in detail.

Encounter Rate
We model the initial distribution of dressed PBHs per mass before any disruption effects as n halo,0 (r, M halo,0 ) = n halo,0 (r) dP 0 (M halo,0 , r) dM halo,0 , where the initial mass function dP 0 (M halo,0 , r)/dM halo,0 is given by Eq. (2.5). From this we obtain the initial mass density of dressed PBHs as ρ halo,0 (r) = (M halo,0 + M PBH )n halo,0 (r), where M halo,0 is the initial average minihalo mass for a PBH of mass M PBH . We assume the initial distribution n halo,0 tracks the DM halo profile in the Galaxy and denote the inital fraction of DM in dressed PBHs before the disruption as f halo,0 ≡ ρ halo,0 (r)/ρ NFW (r). The initial minihalo distribution can then be written as where f PBH = ρ PBH (r)/ρ NFW (r) and we have assumed n PBH = n halo,0 , i.e. all the PBHs have accreted minihalos around them. Including the disruption effects, the minihalo density per mass becomes n halo (r, M halo ) = n halo,0 (r) dP(M halo , r) dM halo , (4.5) where dP(M halo , r)/dM halo is the mass function with the disruption effects accounted for and it depends on r as illustrated in Fig. (4). Using this, we can estimate the encounter rate between NSs and dressed PBHs as follows where n NS (r) = n NS,d (r cyl , z)+n NS,b (r) is the radial number density of neutron stars given by Eqs. (4.1) and (4.2), ... is the average over the NS-dressed PBH relative velocity distribution in the Milky Way halo, and σ eff is the effective cross section of the collision. Since M halo = M halo (R halo ) we have The effective cross section is given by the usual geometrical cross section plus the gravitational focusing enhancement as follows is the mutual escape speed between the minihalo and the NS. We have taken R halo R NS and M halo (M PBH , M NS ) in the rightmost part of Eq. (4.8). The gravitational enhancement, which comes from the deflection of trajectories by the gravitational attraction between the NS and the minihalo, can be written as (4.9) Suppose that dressed PBHs and NSs have velocities dispersion σ PBH and σ NS , respectively, and both follow a Maxwell-Boltzmann distribution so that the probability P (v rel )dv rel that dressed PBHs and NSs have a relative speed in the interval (v rel + dv rel ) reads as (see page 712 in Ref. [89]) where σ rel = σ 2 PBH + σ 2 NS ≈ √ 2× G N M MW (R)/R ≈ √ 2×220 km/s is the relative velocity dispersion assuming a flat Galactic velocity curve and C 0 is the normalization constant so that 4π vesc 0 P (v rel )dv rel = 1, with v esc as a characteristic escape velocity in the Galaxy. In most cases of interest, we have v esc √ 2σ rel leading to C 0 ≈ 1 (e.g., see Fig. 3 in Ref. [99]). Thus, we may estimate the average over the relative velocity of the effective cross section as where the approximation is obtained by pushing the upper limit of the integral to infinity supported by the quick decay of the exponential. (4.12) As shown before, the tidal disruption due to the Milky Way mean field potential and disk shocking both depend on the ratio M halo /M PBH rather than M PBH . Thus the survival mass fraction in Fig. 3 does not depend directly on the PBH mass. However, the encounter rate shows a direct dependence on M PBH through the minihalo radius and the radial number density of dressed PBHs. From Eq. (4.9), we see that the gravitational focusing may be safely neglected for M PBH ∼ (1 − 10)M . On the other hand, R halo ∝ (M PBH /M ) 1/3 and n halo (r) ∝ (M PBH /M ) −1 so that the differential encounter rate runs with the PBH mass as dΓ NS-PBH /dr ∼ (M PBH /M ) −1/3 . The larger the central PBH mass, the smaller the encounter rate due to the decrease of the radial number density of dressed PBHs which dominates over the increase of the cross section. Figure 5 shows the differential encounter rate between NSs and dressed PBHs with respect to the galactocentric radius assuming a Gaussian distribution for the initial mass function of dressed PBHs in Eq. (2.5) and taking into account the mass loss of minihalos due to tidal forces and disk shocking in the Milky Way. We have assumed that dressed PBHs undergo circular orbits around the Milky Way center. The fraction of DM in PBHs is taken to be f PBH = 10 −3 . As shown in Fig. 3, tidal forces acting on dressed PBHs have negligible effects for R 13 kpc so that the quick drop of the differential encounter rate for larger galactocentric radii in Fig. 5 is due to the decrease in the numbers of NSs and dressed PBHs available for collisions. As we approach the innermost parts of the Milky Way there are two factors which compete against each other. From one side the number of NSs and dressed PBHs increases but on the other hand the mass distribution of dressed PBHs moves towards lighter minihalo masses as shown in Fig. 4. Both effects tend to compensate each other at distances of O(100) pc leading to an encounter rate of O(10 −7 ) kpc −1 s −1 . As explained above, the heavier the central PBH mass, the smaller the differential encounter rate at the same galactocentric radii. The blue, red and green solid lines correspond to the cases M PBH = (1, 30, 50)M , respectively, showing total encounter rates for R 100 pc as Γ 100 NS-PBH ≈ (0.17, 0.055, 0.047) day −1 , respectively. These encounter rates are about two orders lower than the rates reported in Ref. [24] for the case of axion miniclusters colliding with NSs, when these dark matter substructures are assumed to have a power-law density profile, ∼ r −9/4 , and constitute the entire DM.

Transient Radio Signal
As the NS goes through the minihalo of a dressed PBH, axion particles fall towards the NS surface and they resonantly convert into photons within its magnetosphere. For the axion mass in the MHz-GHz range, the conversion takes place in a small region around r c > R NS , where r c is called the conversion radius [21]. This particular radius corresponds to the place at which the plasma frequency ω p matches the axion mass m a . For simplicity and for the sake of comparison with Ref. [24] where aligned NSs are considered, we focus on the case of aligned or slightly oblique NSs by assuming that the misalignment angle between the magnetic field and the NS rotating axis is small, θ m 1. In this case, the magnetosphere can be reasonably approximated by the Goldreich-Julian model [100] and the power radiated per unit solid angle at the observation angle θ and to zeroth order in θ m in the WKB and stationary phase approximations reads as [21] dP(θ) dΩ ≈ g 2 aγγ B 2 0 ρ a (r c )πv c 6m a Here B 0 is the magnetic field strength at the NS poles, ρ a (r c ) and v c are the axion density and velocity at a distance r c from the NS, respectively, and P is the NS spin period. Now we consider the individual trajectories of the axion particles close to the NS surface. The maximum impact parameter b a, max , which ensures that an axion particle ends up crossing the conversion radius, is given by At distances of the order of b a, max , the axion particle mostly feels the NS gravitational potential and the gravitational potential from the effective mass located at the center of the dressed PBH may be safely neglected. In addition, we have neglected velocity dispersion of the axion particles within the minihalo because, using the isothermal approximation, we have v DM (r) 2 = 3G N m halo (r) 2r The above analysis only breaks down if the impact parameter between the NS and the center of mass of the dressed PBH is extremely small but such encounters are statistically disfavored. Suppose that the NS crosses the minihalo with an impact paramater b as indicated in Fig. 7 (left). The crossing time is estimated to be T cross = (2/v rel ) R 2 halo − b 2 . Substituting the expressions for v c and r c into Eq. (5.1), we obtain Here we have taken ρ a (r c ) ≈ ρ halo ( R halo ) given by Eq. (2.3) for r c R halo R halo and The allowed angles depend on the axion mass and the NS spin period and the magnetic field at poles. Figure 6 (left) shows the ratio r c /R NS in terms of the polar angle for several axion masses with all other parameters fixed. As the axion mass increases, the angular regions corresponding to the absence of resonant axion-photon conversion slowly begin to grow until they eventually extend over the full angular space, i.e. Θ = {}.
We take the axion mass as a free parameter within the range m a ∈ [m a,lower , m a,upper ] constrained as follows. The lower limit m a,lower is determined by the upper bound on the axion decay constant, so that m a,lower (F a = 10 12 GeV) GHz. The upper limit m a,upper is set either by (i) the maximum mass for which the conversion radius is larger than the NS radius, Eq. (5.2), or (ii) the maximum frequency detectable by a given radio telescope, depending on which one gives a stronger bound.
The photon flux is peaked around the central frequency ν peak = m a /(2π). Assuming energy conservation, the signal bandwidth is taken to be proportional to the initial DM dispersion as ∆ν ∼ ν peak v 2 DM in Ref. [21]. For the case of axion DM background, this bandwidth is typically of the order of kHz, but in our setup where the initial DM dispersion is set by Eq. (5.7), this would lead to a narrower spectral line. In Ref. [24], where the resonant signal is studied in the context of NS-axion minicluster encounters, the signal bandwidth is also predicted to be very narrow due to the small internal velocity dispersion of these astrophysical objects. These authors argue that the signal bandwidth is expected to be widened by other mechanisms and fix the bandwidth of the signal to 1 kHz, which equals the resolution of current and prospective radio telescopes [101][102][103]. In the context of oblique NSs and when axions are non-relativistic at the conversion surface, authors in Ref. [104] point out that the broadening of the spectral line should be dominated by the relative motion between the resonant surface and the observer instead of the initial axion DM dispersion. In this case, the bandwidth of the signal would read as (see Eq. (84)  where is a geometrical factor which depends of the NS properties. For an oblique NS, the intersection of the conversion surface with a plane perpendicular to the NS rotation axis forms an ellipse with eccentricity rather than a circle. As an example, consider the particular case of the isolated NS RX J0806.4-4123 located at about 250 pc distance from the Earth, which has a period P ≈ 11.37 s and a magnetic field at poles B 0 ≈ 2.5 × 10 13 G [105]. For axion masses in the range 300 MHz m a 50 GHz, the signal bandwidth ∆ν ≈ kHz for (0.01− 0.03) and ∆ν ≈ 100 kHz for (0.1 − 0.3). For the sake of comparison with Ref. [24], we will first set the signal bandwidth equal to 1 kHz and analyze the parameter space for detection in this case. In [24] the value 1 kHz was taken to estimate the bandwidth broadening due to astrophysical effects not specified to further detail. We follow the same phenomenological approach here but point out in addition that requiring the Doppler broadening effect (5.13) to be below 1 kHz level requires small eccentricity, ∼ 10 −3 for the parameters shown in our figures. Later, we will also consider larger bandwidths dominated by the Doppler effect Eq. (5.13) to compare our results with those in Ref. [104].
The spectral flux density is calculated as S(θ) = (dP(θ)/dΩ)/(d 2 ∆ν), where d is the distance from the encounter to the Earth and ∆ν is signal bandwidth as explained above. Due to the physical extension of the minihalo, we have a transient radio signal as NS crosses it. We will present our results in terms of the mean spectral density S(θ) = (1/∆t) S(θ)dt, where the time lapse is given by ∆t. From Eq. (5.8), we see that the spectral flux has a strong dependence on the polar angle, so that S(θ) ∈ [S θ,min , S θ,max ], where S θ,min ≡ S(θ = π/2) and S θ,max ≡ lim θ→θ 1 S(θ) = lim θ→θ 4 S(θ). The signal is very sharply peaked around the maxima, as demonstrated in Fig. 6 (right) for a given parameter set, and tracking a maximum for an extended time period with a sufficiently high angular precision may be technically challenging. For the sake of completeness, we will present our results using different angles including θ = π/2. As the NS crosses the dressed PBH, the axion density around it will change according to the density profile of the minihalo, Eq. (2.2). Considering an impact parameter ranging as 0.01R halo b R halo , the typical crossing time is O(10 11 s) for the minihalo mass range of interest. Figure 7 (right) shows the contour plot of the crossing time in the parameter space of (M halo /M PBH , b/R halo ) at the local neighborhood. The smaller the impact parameter or the larger the minihalo mass in units of the central PBH mass, the larger the crossing time as expected.

Estimates for radio telescope sensitivity
We estimate the current and projected sensitivity from radio emission during a NS-dressed PBH encounter. The minimum detectable flux for a given radio telescope reads as [22] S min = SNR min SEFD η s √ 2∆B∆t obs , (5.14) ∼ 220 µJy where SNR min is the minimum signal-to-noise ratio, SEFD ≡ 2k B T sys /A e is the system equivalent flux density (where T sys and A e are the system temperature and the effective area, respectively), ∆B is the bandwidth, η s is the system efficiency, and ∆t obs is the observation time.
The Expanded Very Large array (EVLA) located in New Mexico comprises 27 independent antennas and covers frequencies from 1 GHz to 50 GHz. The antenna SEFD is O(100) Jy with the minimum in the X-band (central frequency: 10 GHz) and maximum in the Q-band (central frequency: 45 GHz). The antenna SEFD in each band is reported in Table 1 of  Table 7 of Ref. [101]. From now on, we assume that the optimized bandwidth in Eq. (5.15) matches the signal bandwidth so that ∆ν = ∆B. Regarding the NSs properties, a Monte Carlo-based population synthesis is used by Faucher and Kaspi [92] to model the birth properties and time evolution of pulsars in the Parkes and Swinburne Multibeam surveys. While the magnetic field at the NS pole follows a log-normal distribution ranging as ∼ 10 12 − 10 14 G, the spin period follows a normal distribution ranging as ∼ (0.04 − 10) s (see Figure 6 in Ref. [92]). As discussed in Refs. [1,21], potentially good targets for axion detection via resonant conversion are the so-called Magnificent Seven, which corresponds to a group of ∼ 7 nearby isolated NSs holding strong magnetic fields (∼ 10 13 G), long spin periods (∼ 5 s), short distances from the Earth (less than 500 pc), and absence of non-thermal emission and observed radio emission [105,107,108]. From Eq. (5.9), we see that such NS properties would enhance the flux during NS-dressed PBH encounters. In addition, while the lack of non-thermal emission allows us to assume the validity of the Goldreich-Julian model within the NS magnetosphere, the lack of radio pulsarlike emissions allows us to estimate the projected sensitivity to the axion-photon coupling neglecting the NS background radiation [21]. Within this group, we find the isolated NS RX J086.4-4123 which we mentioned before [105]. We are assuming that the NS is crossing the minihalo zone at which the flux is maximum for an impact parameter b = 0.01R halo , i.e. in the middle of the observing period ∆t/2 the distanceR halo is minimised, see Fig. 7 (left). The zones for detection for the SKA1-mid array [101] and EVLA [102] are shown in the shaded orange and blue bands, respectively. Each group of 10 flux curves is obtained by varying the NS spin period from 1 s (lowermost curve) to 10 s (uppermost curve). We have shown the particular cases for the minimum mean flux associated with the KSVZ (green curves) and DFSZ (purple curves) axion models and the mean flux at θ = 0.93 rad for the KSVZ model (pink curves). In addition, the solid (dashed) black line shows the estimated mean spectral flux density of the particular NS RX J0806.4-4123, assuming the KSVZ model and a polar angle θ = 0.93 rad (0 rad). Generally speaking, for sufficiently large axion mass, we see that the signal associated with both axion models should be detectable by both SKA-1 mid and EVLA telescopes when we consider ∆B = 1 kHz. In the bottom panel of Fig. 8, the effect of the impact parameter is analyzed by taking otherwise the same parameters as above but using just the KSVZ model and fixing the axion mass to 30 GHz. The group of 10 blue (7 orange) curves are obtained by taking B 0 = 10 14 G and θ = 0.6 rad ( B 0 = 3 × 10 13 G and θ = 0 rad) and varying the NS spin period from 1 s to 10 s (1 s to 7 s). In addition, the solid and dashed red lines show the particular case of RX J0806.4-4123, assuming the KSVZ model, m a = 10 GHz, and θ = π/2 and 0.8 rad, respectively. The green (brown) shaded region corresponds to the detectability zone for EVLA Band Ka (X), which is sensitive to the spectral line associated with m a = 30 GHz (10 GHz). Generally speaking, the smaller the impact parameter, the larger the flux as the NS crosses inner regions of the minihalo with higher densities. The axion detection is favored for b ∼ O(0.01)R halo for the predefined observation time and bandwidth. Figure 9 shows the projected sensitivity to the axion-photon coupling constant as a function of the axion mass for the isolated NS RX J0806.4-4123 during one encounter with a dressed PBH with M halo = 10 2 M PBH . We use a single-dish Arecibo-like telescope with a typical SEFD = 2 Jy in the frequency range from 300 MHz to 10 GHz, η s = 0.5, and SNR min = 1. We are estimating the corresponding mean spectral flux density as S(θ) ≈ S(θ) for ∆t = ∆t obs = 100 hr and ρ a (r c ) ≈ ρ a (r), where r = 10 −3 R halo is the radial distance of the NS from the central PBH. The gray and green shaded bands correspond to the exclusion regions from the CAST solar axion experiment [16] and the ADMX haloscope collaboration [14], respectively. The yellow shaded region is the QCD axion parameter space, where we have indicated with a dashed and solid orange lines the particular KSVZ and DFSZ models, respectively. For ∆B = 1 kHz, the red, brown, and black solid lines show the projected sensitivity for RX J0806.4-4123 assuming a polar angle θ = (π/2, 0.85, 0.98θ 1 ) rad, respectively, where θ 1 = θ 1 (m a , P, B 0 ) is calculated using Eq. (5.12). For ∆B = ∆ν calculated using Eq. (5.13), the solid and dashed blue lines show the projected sensitivity for RX J0806.4-4123 assuming 2 = 0.1 and a polar angle θ = π/2 and 0.9 rad, respectively. We see that most part of the sensitivity curves are within the QCD axion parameter for different angles of observation and 2 ≤ 0.1. We assume the NS is crossing the zone within the minihalo at which the flux is maximum for b = 0.01R halo . Each group of 10 curves (pink, green, purple) is obtained for P from 1 s (lowermost curve) to 10 s (uppermost curve). The shaded blue (orange) band is the detectability zone for EVLA (SKA1-mid). (Bottom) S for the KSVZ model and m a = 30 GHz as function of b/R halo . The group of 10 blue (7 orange) curves refers to B 0 = 10 14 G, θ = 0.6 rad (B 0 = 3 × 10 13 G, θ = 0 rad) and are obtained for P from 1 s to 10 s (to 7 s). The lowermost (uppermost) curve refer to the minimum (maximum) spin period in each group. The green (brown) shaded region is the detectability zone for EVLA Band Ka (X). We have added in both panels results for RX J0806.4-4123 using different θ angles and axion masses. : Projected sensitivity to g aγγ as a function of the axion mass m a for ∆t obs = 100 hr, η s = 0.5, SNR min = 1, and SEFD = 2 Jy for an Arecibo-like telescope. The shaded yellow band is the QCD axion parameter space, while the dashed (solid) orange line is the particular KSVZ (DFSZ) model. Exclusion regions from CAST [16] and ADMX [14] are shown by the gray and green shaded bands, respectively. The red, brown, and black solid lines correspond to the projected sensitivity for the isolated NS RX J0806.4-4123 using θ = (π/2, 0.85, 0.98θ 1 ) rad, respectively, and ∆B = ∆ν = 1 kHz. The solid and dashed blue lines are calculated using ∆ν from Eq. (5.13) with 2 = 0.1 and θ = π/2 and 0.9 rad, respectively. For all cases, we are assuming a NS passing through a dressed PBH (M halo = 10 2 M PBH ) at a distance r = 10 −3 R halo from its central PBH.

Discussion and Conclusion
In this paper, we have discussed for the first time a novel way to detect the QCD axion by means of transient radio signatures from resonant axion-photon conversion during encounters of NSs and axion minihalos around PBHs in the Milky Way. In the scenario where the PQ symmetry is broken before (or even during) the inflation and a small fraction of DM is composed of PBHs, they will unavoidably acquire minihalos from the axion background. Thus, dressed PBHs will end up in the Milky Way halo at z ∼ 6 to undergo later different levels of disruption. Mainly motivated by the LIGO-Virgo gravitational waves detection [26,109] and the recent NANOGrav results [27], as a Benchmark model, we take f PBH = 10 −3 as the initial fraction of (naked) PBHs with masses M PBH = O(1 − 10) M . For such range in PBH masses, the main sources of disruption are tidal forces due to the mean field potential of the Milky Way and the gravitational potential during disk crossing. Under the simplification of circular orbits around the Galactic center, we estimate the mass loss of minihalos depending on their Galactocentric radii. Assuming an initial Gaussian minihalo mass distribution, we conclude that the mass loss is negligible for radii R 13 kpc, but as we approach the Galactic center the loss increases reaching up to ∼ 90% of the original mass for M halo ∼ (135 − 270)M PBH at R ∼ (0.1 − 2) kpc. The heavier the minihalo in units of the central PBH and the smaller its orbital radius, the larger the mass loss due to disruption in the Milky Way.
Taking into account the disruption effects that act on the dressed PBHs, we estimate the differential encounter rate between NSs and dressed PBHs as function of the Galactocentric radius. In the local neighborhood, we find dΓ NS-PBH /dr ∼ 2×10 −7 kpc −1 s −1 (M PBH /M ) −1/3 . As we approach to the inner parts of the Milky Way the number of NSs and dressed PBHs increases but the mass distribution of dressed PBHs also moves towards lighter minihalo masses (and correspondingly smaller minihalo diameters). The two effects tend to compensate each other at distances below leading to a differential encounter rate of ∼ 10 −7 kpc −1 s −1 . At such small distances our assumption of re-viralization of the minihalos between disk crossings starts to break down which is expected to lead to certain overestimation of the differential encounter rate. This should not affect our final results, however, as the small distance regime gives only a subleading contribution to the total encounter rate. The encounter rate scales with the PBH mass as Γ NS-PBH ∼ (M PBH /M ) −1/3 , so that the larger the central PBH mass the smaller is the encounter rate between NSs and dressed PBHs. The minihalo radius grows as ∼ (M PBH /M ) 1/3 but the number density of dressed PBHs in the Milky Way halo decreases as ∼ (M PBH /M ) −1 . Therefore, when M PBH is increased, the decrease of the number of dressed PBHs available for encounters dominates over the increase of the collision cross section. As the total encounter rate, integrated up to R ≥ 100 pc, we find O(10 −2 ) day −1 in the PBH mass range of our interest. Our estimates for the rate of NS-dressed PBH encounters could be improved by including eccentric orbits around the Galactic Center in the study and/or performing numerical simulations to statistically track the evolution of the minihalo initial mass distribution in the Milky Way halo, we leave these tasks for a future work.
By estimating the mean spectral flux density produced by the NS-dressed PBH interaction and comparing it with the sensitivity of current and prospective radio telescopes, we show that the transient radio signals should be detectable on the Earth under suitable conditions. Because of the power radiated per unit solid angle scales as dP/dΩ(θ = π/2) ∼ (ρ halo )(m a ) 4/3 (B 0 ) 5/6 (P ) 7/6 (g aγγ ) 2 , the signal is the stronger the smaller the impact parameter, and the larger the magnetic field, spin period, axion-photon coupling constant and axion mass.
Interesting targets for possible detection would be nearby isolated NSs having long spin periods, strong magnetic fields, and absence of non-thermal emission and pulsar activity. We use as a particular example, the isolated NS RX J0806.4-4123, which is part of the Magnificent Seven group [105,107,108]. Based on the sensitivity of current and prospective radio telescopes, we show in Figs. 8 and 9 that RX J0806.4-4123 could be detectable on the Earth if it undergoes a close encounter with a typical dressed PBH. Since the mean spectral flux density is proportional to the axion density at the conversion radius and the minihalo density quickly increases as we approach the inner shells, the minimum observation time required for detection significantly decreases if we assume smaller impact parameters. In Fig. (8) we show that the transient radio signal coming from a RX J0806.4-4123-dressed PBH encounter could be detectable on the Earth assuming EVLA and SKA1-mid sensitivities, b ∼ O(0.01)R halo , a kHz-bandwidth and 1 year of observation. In Figure 9 we show that the projected sensitivity to the axion-photon coupling constant, assuming an Arecibo-like telescope and 100 hours of observation time, is mostly within the QCD axion parameter space, including the particular KSVZ and DFSZ models, for a bandwidth of the order of (1 − 10 2 ) kHz. Note that our sensitivity results for RX J0806.4-4123 with a bandwidth ∼ 10 2 kHz are more than 3 orders of magnitude stronger than those of Ref. [104], which is currently one of the state-of-the-art estimates. This is because the axion minihalo density at a distance r = 10 −3 R halo from the central PBH is about 8 orders of magnitude higher than the local DM density and the sensitivity of the axion-photon coupling constant scales as g aγγ ∼ 1/ ρ a (r c ).
Motivated by the enhancement in the signal via the enhancement in the axion DM background, there has been much interest in the literature for studying the signal coming from axion stars located close to the Galactic center, such as the magnetar SGR J1745-2900. This astrophysical object is located at about 0.1 pc from the Galactic center and holds a strong magnetic field at poles B 0 ≈ 1.6×10 14 G, and a long spind period P ∼ 3.76 s [110][111][112][113]. However, such kind of an environment is not the ideal in our setup since the strong tidal forces coming from the Milky Way's central supermassive black hole would remove most of the DM shells from minihalos.
Extragalactic encounters of NSs with dressed PBHs may potentially also be detectable on the Earth if we consider close encounters which lead to a final merger. In such a scenario, we could have a multi-messenger probe of the QCD axion via the observation of both the axion photon conversion signal and gravitational waves from PBH-NS binary inspirals, as studied in detail in Ref. [67] for the case of intermediate mass black holes with axion DM spikes.
Even though locations of most of the NSs are not know, we consider that a good starting point for searches should be along the Milky Way disk and towards the Galactic center (but sufficiently far away to avoid sizeable minihalo disruption) based on the NSs sub-population on the disk and bulge. On the other hand, resolving the signal from radio foregrounds is expected to be a general challenge for such searches involving large fields of view but analysis of this topic goes beyond the scope of the current work.
Recently, a study of radio signals coming from the axion-photon conversion during NS-axion minicluster encounters was performed in Ref. [24]. According to the sensitivity of current and prospective radio telescopes, they showed that the QCD-axion should be detectable assuming an axion mass of ∼ 5 GHz and a power-law density profile for axion miniclusters with large mean density, ∼ (10 4 −10 8 )M /pc 3 . If axion miniclusters are assumed to have a NFW density profile, the associated emitted flux would not be large enough to be detectable (for a discussion about the axion minicluster density profile in the context of Nbody simulations see, for example, Ref. [114]). If the PQ symmetry breaking occurs before or during the inflation the axion fluctuation is too small to form axion miniclusters, but if the PQ symmetry is broken after the inflation, large density fluctuation of the axion will be produced. However, it seems uncertain how large axion fluctuations are produced after the PQ symmetry breaking. Our scenario on the other hand assumes the presence of PBHs which is certainly uncertain, although the LIGO-Virgo observations and recent NANO-Grav results may provide some suggestive hints of their existence.
To conclude, our results give rise to the striking prediction of a transient characteristic line-like emission which would need to be confirmed via dedicated search strategies. A more careful analysis about the final bandwidth of the radio signal and model for the conversion process happening in the NS magnetosphere should be done. We leave these tasks for future work.