Anatomy of astrophysical echoes from axion dark matter

If the dark matter in the Universe is made of μeV axion-like particles (ALPs), then a rich phenomenology can emerge in connection to their stimulated decay into two photons. We discuss the ALP stimulated decay induced by electromagnetic radiation from Galactic radio sources. Three signatures, made by two echoes and one collinear emission, are associated with the decay, and can be simultaneously detected, offering a unique opportunity for a clear ALP identification. We derive the formalism associated with such signatures starting from first principles, and providing the relevant equations to be applied to study the ALP phenomenology. We then focus on the case of Galactic pulsars as stimulating sources and derive forecasts for future observations, which will be complementary to helioscopes and haloscopes results.


Introduction
QCD axions are hypothetical particles introduced to solve the so-called strong CP problem [1][2][3][4], namely to explain the extremely small value of the neutron's electric dipole moment.This pseudo-scalar Nambu-Goldstone boson, associated with the Peccei-Quinn symmetry breaking, can be produced in the early phase of the Universe, before or after inflation [5][6][7].In this way, the axion can also be a good candidate for cold dark matter (DM), for values of the Peccei-Quinn scale, f a , of the order of 10 10 − 10 12 GeV.Such values of f a , in turn, correspond to an axion mass, m a , in the rage 10 −6 − 10 −3 eV.QCD axion represents a top-priority target for DM, and many detection strategies in the laboratory and in the sky have been designed in order to find it [8], exploiting its very weakly couplings with Standard Model (SM) particles and, in particular, with photons.
It was later realized that axion-like particles (ALPs) are also present in many beyond-thestandard-model theories, such as string theory models [9].These "cousins" of the QCD axion can be (pseudo-)scalar particles, with masses as low as zeV, with very weak couplings to the SM.As QCD axions do, they couple to photons through the Lagrangian term L ∝ g aγ aF µν F µν .Moreover, they can also be good DM candidates in some portions of the parameter space.
Thanks to their coupling with photons, besides the well-known possibility for an ALP to convert into a photon in the presence of external magnetic fields [10], an ALP has a finite probability to produce two photons from its decay.Moreover, in the presence of an ambient radiation field, a stimulated enhancement of the decay rate occurs.In the radio band and for certain astrophysical environments, this can amplify the photon flux by several orders of magnitude.The observational prospects concerning the ALP stimulated decay inside a photon emitting source, such as dwarf spheroidal galaxies, the Galactic Center, and galaxy clusters, were studied in [11,12].The two photons arising from the decay are emitted backto-back, i.e. one with the same direction as the stimulating photon and one in the opposite direction (or, more precisely, approximately in the opposite direction since the ALP velocity is non-relativistic, but not zero).This means that stimulated ALP decay generates "echoes", i.e. photons traveling in the opposite direction with respect to the flux of the stimulating source.The production and detection of the ALP echo from an artificial beam was proposed in [13] with in-depth calculation in [14] and a possible concrete experiment discussed in [15] (see also [16,17]).It was soon realized that also "natural" beams, from the emission of astrophysical sources, could induce echoes in the radio sky [18].In particular, the potential of the echo signatures in our Galaxy was explored for supernova remnants, which are particularly promising since they could have been significantly brighter in the past [19,20].
The main aim of this work is to present a full and complete derivation of the astrophysical echoes from ALP DM, starting from first principles, and determining the dependence of the signal on the DM properties (such as the spatial and velocity distributions) and the source characteristics (such as distance, age, and motion).We will discuss three signatures, made by two echoes and one collinear emission, associated with the decay, which can be simultaneously detected, offering a unique opportunity for a clear ALP identification.
The paper is organized as follows.In Section 2, we describe the derivation of the echo signal from a point source, with additional details of the computation reported in Appendices A-C.Section 3 describes the observational properties of the echoes and their dependence on the DM halo and source characteristics.In Section 4 and Appendix D, we treat the case of collinear emission.Then we focus on the case of Galactic pulsars as stimulating sources, deriving observational prospects in Section 5. Section 6 concludes.

Echo signal from a point source
We want to describe the following process: a photon emitted by a source with momentum ⃗ q stimulates the decay of a DM ALP1 with momentum ⃗ p. Two photons are produced in the decay: one with momentum ⃗ q, identical to the incoming photon, and one with momentum ⃗ k, the echo photon.In this Section, our task is to derive an expression for the flux density of the echo as seen from Earth.We will treat the collinear emission in Section 4.
The Boltzmann equation for the axion number density n a tells us how many axions decays  occur per unit time: (2.1) Here λ |M 0 | 2 is the matrix element squared, summed over the polarizations of the two final state photons, and f a and f γ are the phase space distributions of the axions and of the photons emitted by the source.Notice that the enhancement factor from phase space density of photons with momentum ⃗ k does not appear in Eq. (2.1) because we are considering a configuration in which this state is not initially highly occupied.
Let's focus on a point source located at the origin of our coordinate system emitting radially and isotropically.We borrow part of our notation from Ref. [20].We consider an ALP decay happening at a location ⃗ x ds that we choose to lie on the xz-plane.At ⃗ x ds , the photons from the source have a unique direction with θ s determined by ⃗ x ds .We take the Earth to lie on the z-axis.The geometry is shown in Figures 1 and 2. We consider two possible configurations: 1) the back-light echo coming from the direction opposite to the source, shown in Figure 1 and 2) the front-light echo coming from behind the source, shown in Figure 2.
For the axion phase space distribution, we choose a Maxwell-Boltzmann where ⟨⃗ p ⟩ is the average DM momentum in Earth's rest frame and √ 3δp(⃗ x ) is the total momentum dispersion at a location ⃗ x assuming no anisotropy.The axion number density is given by From Eq. (2.1), with Eqs.(2.2) and (2.3), through the derivation described in Appendix A, we obtain with where we have used the fact that for non-relativistic axions ⃗ p = m a ⃗ v.In the equations above, Ω k is the direction pointing from the decay location to the detector, ⟨v ∥ ⟩ and ⟨v ⊥ ⟩ are the components of the axion average velocity parallel and perpendicular to the line of sight (los) in the xz-plane, ⟨v t ⟩ is the average DM velocity perpendicular to the xz-plane, ε = 2ω k /m a − 1, and θ d is the angle between the los and the direction from the source to the decay location.In Eq. (2.5), n a , f γ , ⟨⃗ v ⟩ and θ d are all functions of the decay location ⃗ x ds .Explicit expressions for the components of ⟨⃗ v ⟩ are given in Appendix C.
The interpretation of Eq. (2.5) is the following.At the decay location ⃗ x ds , the source's photons arrive from a certain direction θ s .Any ALP present at ⃗ x ds can be stimulated to decay, but only if the ALP in question has velocity components v ∥ = ϵ, v ⊥ = ω k θ d /m a and v t = 0 will the echo photon travel in the direction of Earth with frequency ω k .The likelihood of such a decay depends then on how many axion particles are available with the required momentum.
For a fixed observation angle θ i ≪ 1, the deflection angle θ d takes the following form where the upper sign corresponds to the back-light echo, the lower to the front-light echo.
For the back-light echo, x s < x ds , so that a small θ i implies a small θ d .For the front-light echo, instead, x s /x ds can be large if the decay happens close to the source.Thus, in this case, a small θ i does not necessarily imply a small θ d .However, for non-relativistic axions, decays with large θ d are suppressed by the axion phase space distribution.So, in both cases, we can approximate The decay rate is the largest at a frequency and is suppressed exponentially for ω − ω * > m a δv.Since ⟨⃗ v ⟩ is small, we can approximate ω k ∼ m a /2 in Eq. (2.5), keeping corrections of order ⟨⃗ v ⟩ only in the exponents.
The rate of change of the axion number density can be related to the echo flux density observed from Earth.For each axion decay, two photons, each of energy ∼ m a /2 are produced.Only one of those travels toward Earth.The power emitted towards us per unit volume, unit frequency, and unit solid angle is then equal to m a /2 |d ṅa /(dν k dΩ k )|.Integrating along the los, we obtain the spectral intensity where we assumed absorption effects to be negligible.Next, we need to relate the source's photon phase space density to its flux density observed from Earth.Using Eq. (2.2), we obtain where ρ γ is the photon energy density.The flux density through a surface area perpendicular to the direction Earth-source is where S (0) ν is the source's flux density as observed from Earth, Ω i is the solid angle in the direction of the los (θ i , ϕ i ), and B describe the observing beam.In Eq. (2.12), we have not taken into account two possible effects: the motion of the source with respect to Earth and the time dependence of the source's luminosity.If we introduce these effects, Eq. (2.12) generalizes to ) where xs (t) is the distance from Earth to the actual position of the source at time t (not to the position where the source appears to be), xds (t) is the distance from the source to the decay location at time t, and L ν is the luminosity of the source.At time t em = t − x d − xds (t em ), the source emitted a photon that subsequently traveled a distance xds (t em ) to the decay location, stimulated the decay of an ALP which in turn created a photon that traveled a distance x d to Earth, arriving at a time t.The time t em can be expressed as a function of x d and the velocity of the source ⃗ v s as described in Appendix B. The function h is evaluated at a fixed location, determined by θ i and x d independent of t em .However, it develops a time dependence through θ d .This will become clear in Section 3.2.2.

Angular distribution of the echo on the sky
To study the angular distribution of the signal in the sky, we isolate from S ν the following dimensionless quantity where ρ ⊙ and δv ⊙ are the DM density and velocity dispersion at the Sun's location.

Idealized case
As a first step, we want to establish a basis onto which to add the various effects relevant to the distribution on the sky of the echo intensity.To this end, we consider the academic case of an infinite DM halo with constant energy density ρ a (⃗ x ds ) = ρ ⊙ , velocity dispersion δv(⃗ x ds ) = δv ⊙ , and null DM and source velocity ⟨⃗ v ⟩ = 0. We also assume constant source luminosity L ν (t) and velocity ⃗ v s = 0. We fix the energy such that we are at resonance ω = ω * .Under these assumptions, g(θ i , ϕ i ) reduces to A first observation is that, if in addition to requiring that the DM is at rest on average, we also require it to be perfectly cold, we obtain g id ∝ δ(θ i ).This is to be expected because, in the absence of average velocity and velocity dispersion, the echo travels back to the location where the source was when the photons that stimulated the decay were emitted.For a source at rest emitting radially, the only way to receive the echo radiation is to be between the source and the decay location.
When velocity dispersion is present, the intensity on a spherical shell of radius x d centered at Earth drops exponentially as θ i increases beyond xs x ds θ i ≳ 2δv.For fixed x d , we thus have the typical size of the intensity patch where we have approximated x ds = x d ± x s .Again, the upper sign holds for the back-light echo, the lower sign for the front-light echo.For given x d , the intensity of the front-light echo is more compact than that of the back-light echo, because the distance to the source is smaller.We also notice that, for decays happening on shells of increasing x d , the size of the patch becomes larger, as shown in Figure 3. Correspondingly, for fixed θ i , the back-light echo emission is peaked at x d = 0 if θ i ≲ 2δv, while for larger observation angles θ i > 2 √ 2δv, the decay needs to happen at larger distances in order to have the same θ d .In this case, the largest contribution is from decays happening at and the flux is suppressed by the distance.For the front-light echo, the flux from a given direction θ i is always suppressed at distances x d < x s (θ i /(2 √ 2δv) + 1).In both cases, the integrand decreases as x −2 d at large distances.Carrying on with our idealized case, we can integrate over x d analytically in the approximation x ds ≈ x d ± x s .The lower integration limit is Earth's radius for the back-light echo, while for the front-light echo we take r min = x s .Sending R → ∞, the result is For the back-light echo, we can neglect r min compared to x s .Eq. (3.4) reduces to Notice that the result is finite for For the front-light echo we can set the error function in Eq. (3.4) equal to 1, obtaining The divergence at θ i = 0 is due to the assumption of a point source.The source's flux density diverges at the source location.In the case of the front-light echo, we can receive photons produced in decays happening very close to the source (x ds ≈ 0), which cause the divergence at θ i = 0, while this is not the case for the back-light echo.In practice, the finite size of the source will regulate the divergence.
The idealized distribution on the sky is shown in Figure 4.The intensity is almost constant for θ i ≲ δv for the back-light echo, while for large θ i it drops as θ −1 i .The typical size of the intensity patch on the sky is 10 ′ , taking δv ∼ 10 −3 .For the front-light echo, the intensity decreases as θ −1 i everywhere.
Note also that the echo emission can be disentangled from the direct radiation from the source in the front-light echo case by using polarization.Indeed the echo polarization is orthogonal to the one of the source's emission.

Realistic case
Now that we have identified the main features of the distribution on the sky of the echo intensity, we proceed to turn on one by one the effects coming from the characteristics of the DM halo and of the source.Whenever an effect is not taken into consideration, the corresponding quantity in Eq. (3.1) is set to its default value.The default values are: ρ(r) = ρ ⊙ , x s = 0.1 kpc, δv(r) = δv ⊙ , ⟨v ⊥ ⟩ = 0, the age of the source τ age = ∞, ⃗ v s = 0 and L = const.In all plots of the next two subsections, we set the ω = ω * and we vary the observation angle θ i leaving the Galactic longitude of the los direction unchanged.Thus, when we talk about negative (positive) θ i , we mean that we are observing directions of lower (larger) Galactic latitude than the source.
We illustrate the effects on four hypothetical sources with properties listed in Table 1.ρ(r)

Effects of the dark matter halo
We consider a Navarro-Frenk-White halo, with parameters ρ s = 0.41 GeV cm −3 and r s = 20 kpc [41], yielding a DM density at the Sun's location ρ ⊙ = 0.5 GeV cm −3 .We calculate δv(r) integrating the spherically symmetric Jeans equation assuming no velocity anisotropy (see for example Ref. [42]).We obtain δv ⊙ = 144 km/s.We assume that the DM is on average at rest in a non-rotating frame attached to the Galaxy.Then the average velocity of DM in Earth's rest frame ⟨⃗ v ⟩ is mainly given by −⃗ v ⊙,ϕ , where ⃗ v ⊙,ϕ is the velocity of the Sun in the tangential direction in the Galactic plane with magnitude 245.6 km/s [43] in the direction of Galactic coordinates (l, b) = (90 • , 0 • ).We neglect the Sun's velocity in the other directions as well as the motion of Earth with respect to the Sun.
We study the effect on the echo of the following quantities: 1. ρ(r), 2. ρ(r) in combination with the distance of the source x s , 3. ρ(r) in combination with x s and the velocity dispersion δv(r), 4. the DM velocity relative to the Sun ⟨⃗ v ⟩.
In Figures 5 and 6, we show each of the effects above in a separate panel, for the back-light and front-light echo, respectively.The thick black line in all panels represents the idealized case g id (θ i ).
On the top left panel of each figure, we see the effect of the DM density profile.The lines for sources A, B, and C are on top of each other because the DM density along the los is the same in these three cases.We observe a minor suppression of the echo flux compared to the idealized case.On the other hand, for source D, the back-light echo flux is more significantly reduced because the los points away from the Galactic center.As expected, we see an increased flux for source D in the case of the front-light echo.
The top right panels show the combined effect of the density profile and the distance to the source.Notice, that if we let x s vary among sources while keeping ρ = ρ ⊙ fixed, there would be no effect on the g (see Eqs. (3.5), (3.7)).The lines for sources A and C are unchanged from the top left panels.For sourced B and D, whose distances are larger than the default value, we see a reduction of the flux for the back-light echo due to the lower DM density along the los.Conversely, for the front-light echo, we notice a substantial increase in flux for source D which is now placed closer to the Galactic center.This is due to two facts: the DM density is higher where x ds ∼ 0, and the contribution from the high DM density at the Galactic center is not suppressed by a large value of x ds .
The bottom left panels display the effect of the variation of the velocity dispersion along the los in combination with the source distance and the variation of the DM energy density.In all cases, except for the front-light echo of source D, the effect is negligible.The velocity dispersion goes to zero at large Galactic radii, causing a reduction of the echo flux from these locations.However, in this case, the flux is already suppressed by the distance and by the low DM density.We see instead a non-negligible effect on the front-light echo of source D, for which the los passes close to the Galactic center, where the velocity dispersion is small.For observation angles small enough to have θ d < δv close to the Galactic center, we see an enhancement of the flux due to the normalization of h(ω, ⃗ x ds ), i.e. there are more axions available to decay that will produce photons with the desired direction.If the observation angle is larger than that, the suppression from the exponential wins over the enhancement from the normalization factor and we observe a reduction of the flux.
The bottom right panel of each figure shows the effect of the average DM velocity.For source A, ⟨⃗ v ⟩ is almost completely along the los.Hence it does not affect the distribution of the echo on the sky, but rather it will shift its frequency.For source B, instead, we have a sizeable value of ⟨⃗ v ⊥ ⟩.According to Eq.(2.6), the echo flux is the largest when ω k θ d /m a = ⟨v ⊥ ⟩ or For the back-light echo, the emission is peaked at a location where θ i has the same sign as ⟨v ⊥ ⟩.Now we have to be careful to interpret this correctly.From Eq. (A.12), we see that ⟨v ⊥ ⟩ has the opposite sign as ⟨v x ⟩ (cos θ ⃗ k < 0).Since θ i grows in the positive x-direction, the peak of the back-light echo emission comes from the direction opposite to that of ⟨v x ⟩.
In other words, the peak of the back-light echo emission is shifted in the direction opposite to that of the average DM velocity.When we look in the direction opposite to source B, the component of the average DM velocity perpendicular to the los points in the direction of decreasing Galactic longitude, i.e. ⟨v x ⟩ < 0, causing the emission to peak at positive θ i .The contrary is true for the front-light echo, for which the largest emission comes from the direction where the DM velocity points.However, in this case, we don't see a shift of the peak, but only an asymmetric distribution of the flux.

Effects of the source's characteristics
In this Section, we study how the properties of the source affect the echo emission.In particular, we take into account the following ones: 1. the age of the source τ age , 2. τ age in combination with the distance of the source x s , 3. τ age in combination with the time dependence of the flux, 4. the motion of the source with respect to Earth.
Results are shown in Figures 7 and 8.The top left panels show the effect of a finite a τ age .The age of the source translates into a reduced upper integration limit in Eq. (3.1).The time it takes a photon to travel the distance x d + x ds cannot exceed τ age .The upper integration limit becomes then R ∼ 1 2 (τ age ∓ x s ) . (3.9) On the top left panel of Figures 7 and 8, we see the expected reduction of the flux for the younger sources.For the older sources, the upper integration limit cuts out locations from -12 -  where the flux is already suppressed by a large value of x ds , and the effect is negligible.
On the top right panel, we add a varying distance for the sources.For the back-light echo, a larger distance implies a reduced upper integration limit, resulting in a reduced flux from sources B and D compared to the top left panel.For the front-light echo, a large x s means a larger R, but also a larger θ d for fixed θ i .For this reason, we see a dimming of the flux compared to the top left panel even if the upper integration limit is larger.
For the time dependence of the sources' luminosity, we assume a power law where t is the time since the birth of the source.A positive α means that the source was brighter in the past.The bottom left panel shows how the flux increases for α = 1 and τ 0 = 10 3 yr.As expected, the effect is larger for the older sources.
Finally, the bottom right panel shows the aberration effect due to the motion of the source.
Assuming that the source moves with a constant velocity ⃗ v s , its position at time t is where now t = 0 is the time today.Let's focus on a given decay location determined by the observation angle θ i and the distance along the los x d .The motion of the source deforms and rotates the triangle that has Earth, the source, and the decay location as vertices.From ⃗ v s , we can obtain expressions for the time-dependent Earth-source distance xs (t), source-decay location distance xds (t) and the direction from which the source's photon reach the decay location θ s (t), ϕ s (t).Explicit expressions are given in Appendix B.
Let us first comment on the effect on the front-light echo shown in Figure 8.We have chosen the velocity of the sources such that source A is moving towards us, meaning that it was further away in the past.In this case, the only effect is a mild increase in x s /x ds and θ s compared to the case of a source at rest.We find that this effect is negligible.
Source B is moving in the positive x-direction of Figure 2, meaning it was further down the x-axis in the past.If we start with a positive observation angle, at any time in the past the triangle was more open, because the angle between the los to the decay location and the los to the source was larger than today.This leads to a suppression of the flux.On the other hand, if we start with θ i < 0, the triangle closes as we go back in time, becoming a line when the source was along the los.This causes an increase in the flux.
To understand if this effect is important, we estimate at which distance along the los the decay must happen for us to be receiving today the echo radiation stimulated by photons emitted when the source was along the line of sight.The enhancement of the flux will be larger if the los crossing happened recently so that the decay location is close and well within the halo.
Setting v s,y = 0 we see from Figures 1 and 2, that the source was along the los when θ s (t) = θ i .Starting from a positive θ i , this happened at a time t em = ∓x s θ i /v s,x , (see Eq. (B.2))This tells us that the zero crossing of θ i has occurred in the past for the front-light echo if θ i and v s,x have opposite signs, the other way around for the back-light echo.The distance of the decay stimulated by photons emitted at t em is which may well be within the DM halo.
Going back to Figure 8, for source C, ⃗ v s is perpendicular to the plane xz-plane.This causes the triangle to open up, causing an overall reduction of the flux.Finally, source D is moving both away from us and in the positive x-direction.The result is similar to that of source B, the effect being reduced because of a smaller magnitude of v s,x .
Let's now turn to the back-light echo in Figure 7.For sources B and D, we see that in this case, the flux is larger for positive θ i .This is expected because now, due to the different geometry, the source crossed the los in the past if θ i > 0.

Collinear emission
In this Section, we consider axion decays happening at locations between us and the source.In this case, the echo (with momentum ⃗ k ) travels back toward the source, while the decay photons with the same momentum as the incoming ones (denoted by ⃗ q ), travel toward us.
The Boltzmann equation now reads Following Appendix D, one obtains with The decay rate is the largest for ε = ⟨v ∥ ⟩, where, as before, ⟨v ∥ ⟩ is the average axion velocity along the direction of the momentum of the photon coming towards us, ⃗ q.
Here Ω q is the angle between the direction of the stimulating photon and the direction of the photon emitted in the ALP decay.The delta-function δ(Ω q ) indicates that the collinear emission comes exactly from the direction of the source, and is not smeared over the sky by the axion momentum dispersion as it happens instead for the echo case.
Proceeding like in Section 2, we obtain the flux density With simple geometrical considerations, one can see that, in the limit of small angles, θ q = θ i (x d + x ds )/x ds , from which it follows δ(Ω q ) = δ(Ω i ) x 2 ds /x 2 s .We can thus recast Eq. (4.4) in a simpler form: where, in analogy to Section 2, we also generalized the equation to include the time dependence of the source's flux.
Compared to Eq. (2.13), the factor L ν (t em )/L ν (t − xs (t)) is not present, because the photons produced in the decay travel together with the ones coming directly from the source, i.e. t em = t − xs (t).Note that in the collinear case, if the source is time-dependent on short timescales, e.g., it has pulsations, the same time-dependency occurs also for the echo signal.
The motion of the source also enters, in principle, through a time-dependent δ(Ω i ).However, as long as the source does not move out of the beam during the course of one observation, the proper motion of the source is completely negligible.

Forecast for pulsars
Now, we apply the formalism we developed in the previous Sections to the case of Galactic pulsars as stimulating sources.Pulsars in connection to ALPs were considered, e.g., in [44], but looking at a signal different from the one considered in our work.Ref. [44] studied the ALP conversion in the pulsar magnetic field.Here we derive prospects for detecting the signals described in the previous Sections.We will focus on the front-light echo and collinear emission, which are peaked at the source location for pointlike sources like pulsars, unlike the back-light echo.Since the back-light echo is spread on angular sized of order 10' for both pointlike and non-pointlike stimulating sources, there is no advantage in considering pulsars rather than, for example, supernova remnants.
Our formalism applies to point sources emitting isotropically.The radio emission of pulsars is however not isotropic, but rather concentrated into a rotating beam confined within the polar cap region.The typical opening angle of the beam is [45] ρ = 5.4 Since the echo intensity is spread over angular distances of order tens of arcminutes, and barring complications due to possible changes of direction of the pulsar's rotation axis or magnetic field lines, we can safely consider the time-averaged pulsar emission to be isotropic in the directions that are relevant to us.A complication may arise due to the pulsar's proper motion.If the source was further away from the los in the past, it may be the case that the angle θ s (t) was larger than ρ.In the following, we remove from our integration over the los locations for which this is the case.We find that this effect is completely negligible.
Since the flux density involves the product S (0) ν x s g(Ω i ), we focus on the pulsars with the largest flux times distance product in the Australia Telescope National Facility (ATNF) pulsar catalog [46] and then take into account all the effects described in the previous Section.
Most of the relevant information is available directly from the ATNF catalog, with the exception of the source's velocity along the los and the time dependence of the source luminosity.From the analysis of Section 3.2.2,we know that the velocity along the los only affects the flux in a negligible way.We thus set it to zero.The time dependence of the flux luminosity is discussed in the next sub-Section.

Modeling of the time dependence of the pulsar flux
The flux from a pulsar typically varies on timescales of order 1 second because of the pulsar's rotation.However, as already noted in [20], this time dependence is not present in the echo flux.In fact, the echo is the aggregate of photons produced in decays that happened at different distances along the los.At these locations, the flux from the source peaks at different times, due to the different travel times from the source location.Since we integrate in Eq. (2.13) over distances of order the size of the halo, the travel time is much longer than the pulsar period and thus we expect the echo flux to be constant in time.Thus, in the following, we consider the pulsar's luminosity averaged over timescales much larger than its period, but much smaller than its age.On the other hand, the time dependence of collinear emission is exactly the same as that of the source's emission.
Pulsars are the brightest at their birth, their luminosity decreasing in time following the spin-down power, Ė, evolution.For the time-dependent evolution of the spin-down power in the magnetic dipole model (braking index n = 3), we consider the following equation, e.g.[47] with β = −2.The constant τ 0 is given by where P 0 and Ṗ0 are the initial spin and spin derivative of the pulsar.It is very challenging to infer the initial P 0 and Ṗ0 distributions.For example, given P and Ṗ today, it is not possible to determine P 0 unless an independent measure of the pulsar's age is available [45].
The authors of Ref. [48] analyzed a population of young neutron stars associated with supernova remnants, for which it is possible to estimate the actual age of the pulsar.They found a log-normal distribution for the initial period, with mean log 10 (P 0 /s) = −1.04+0.15 −0.2 and standard deviation is σ P 0 = 0.53 +0.12 −0.08 .In the dipole model, the magnetic field is constant, since it depends only on the product P Ṗ , so that we can assume that the initial magnetic field distribution corresponds to the present one.From [48], the present B-field distribution is best-fitted by a log-normal with mean log 10 (B/G) = 12.44 and standard deviation σ B = 0.44.We use the mean values of these distributions as benchmarks to evaluate the impact of the pulsar spin-down temporal evolution.
Using the mean values of the initial spin and magnetic field distributions in expression for the characteristic magnetic field [45] we obtain τ 0 = 1.79 × 10 4 yr.
The next ingredient is the relation between the spin-down power and the pulsar's luminosity at a given frequency.Recently, Ref. [49] presented the largest, uniform, census of young pulsars deriving for 1170 objects correlations between the radio (pseudo)luminosity and parameters such as age and spin-down power.They found that, in the L-band, log L L = c Ė log Ė with c Ė = 0.15 ± 0.02.Assuming c Ė to be independent, we obtain that the luminosity of a pulsar follows Eq. (3.10) with exponent α = 0.3.We find that with this value of τ 0 and c Ė the effect on the echo of time dependence of the source's flux is negligible.We thus neglect this effect in our forecasts.In doing this we are being conservative, since, for example, the actual value of τ 0 for a given pulsar may very well be smaller than our benchmark value.

Target selection
To estimate our projected sensitivity, we choose to integrate over the bandwidth that maximizes the signal-to-noise ratio.Assuming noise scales like √ ∆ν, a top-hat bandwidth and ν = m a (1 + ⟨v ∥ ⟩)/4π, we find that the optimal bandwidth is ∆ν/ν = 2.8δv [19].
For the targets we are considering, the effect on the echo of a position-dependent velocity dispersion is negligible.We fix then ∆ν based on the velocity dispersion at Earth.The  frequency-averaged echo spectral intensity is with where we have neglected terms of order θ 2 d or vθ d .The Gaussian integrated over the chosen bandwidth provides f ∆ν = 0.84.The frequency-averaged flux density can be then obtained from the usual ⟨S ν ⟩ = dΩ i B(Ω i ) los d⟨I ν ⟩.Similarly, for the collinear emission we have: In Figure 9, we show the bandwidth-averaged front-light echo spectral intensity at a frequency of 400 MHz for the four pulsars in the ATNF catalog with the largest product of distance times flux.We assume an axion photon coupling g aγ = 10 −11 GeV −1 .We take the flux from the best fit to observations from the recent catalog of [50].We note that the flux at 400 MHz is slightly lower than the one reported in ATNF.In Figure 10, we show slices of constant Galactic longitude and latitude to ease the comparison among pulsars.
All effects discussed in Section 3.2 are implemented with the exception of the time dependence of the source's luminosity.We have set v s,z = 0, as this quantity is not known and we have established it is not important for the echo.The age of the pulsar used to determine the upper limit of integration is taken to be the characteristic age τ age plus the pulsar distance.
The Vela pulsar J0835-4510 is the brightest source at this frequency with a flux of about 5 Jy.Its front-light echo flux does not experience significant reduction due to any of the effects discussed in Section 3.2.There is some dimming compared to the idealized case due to the young age τ age = 1.13 × 10 4 yr, and, to a lesser extent, to the DM density along the los.The DM velocity is essentially along the los, thus affecting the resonance frequency, but not the distribution of the echo on the sky.The projection on the sky of ⟨⃗ v⟩ is shown by the blue arrow in Figure 9.The small velocity of the source, ⃗ v s ∼ 60 km/s, shown by the magenta arrow, makes the intensity larger for positive values of ∆l, as also visible on the right panel of Figure 10.
Pulsar J0332+5434 has the largest x s S (0) ν , a factor 1.8 larger than Vela.However, its Galactic longitude l = 144 • , and distance x s = 1.695 kpc, imply that the average DM density along the los is low making it lose its "source advantage".The average DM and source velocities point approximately in the same direction.From the discussion of Section 3.2.1,we know that for the front-light echo, the intensity is larger in the direction of ⟨⃗ v⟩, while in Section 3.2.2we have learned that the intensity is larger in the direction of −⃗ v s (where the source was in the past).For J0332+5434 the two effects tend to compensate, resulting in an intensity patch similar to that of Vela.
Next, we have J1644-4559 and J1645-0317.These two pulsars have similar fluxes at 400 MHz of about 0.4 Jy and comparable distances of order 4 kpc.Their Galactic coordinates are (339 • , −0.195 • ) and (14.1 • , 26.1 • ), respectively, implying that the DM density along both lines of sight is comparable, as well as the projection of the DM velocity on the sky.The main difference between the two is due to our description of their proper motion.For J1644-4559, proper motion information is not available and we set its velocity to zero.J1645-0317 on the other hand has a large velocity of about 375 km/s which causes a drastic reduction in the flux in the direction of ⃗ v s .This discussion highlights the importance of our careful treatment of the signal, which we demonstrate to vary on a source-by-source basis.Different targets therefore reveal a very diverse signal phenomenology.

Observational prospects
For our forecasts, we consider two interferometric observatories of the near future: SKAO for the Southern hemisphere and LOFAR 2.0 for the Northern hemisphere.
A further development of the SKAO, which used to be referred to as SKA-2, is foreseen for the 2030's, with telescope specifications still to be decided.According to [52], the goal is to reach a factor of ×10 improvement in sensitivity with respect to SKA-1.
For very low-frequency emission and/or Northern hemisphere targets, we consider LOFAR Figure 11: Projected sensitivity on the ALP-photon coupling g aγ versus ALP mass m a , for different scenarios, described in Section 5.3.The upper colored cases correspond to observations of single pulsars with either SKA-1 or LOFAR 2.0.More precisely, we show the cases of J0835-4510 (blue), J1644-4559 (green), J1645-0317 (red), and J0332+5434 (orange), including both the front-light echo signal (dashed) and collinear emission (solid).The two black curves report the sensitivity for surveys of all Galactic pulsars with SKA-1 (dashed) and SKA-2 (solid) considering only the collinear case.
The first case shown in Figure 11 concerns the sensitivity forecast for 100 hr observations of single pulsars (colored curves).We consider the same four pulsars described in Section 5.2.The pulsars J0835-4510, J1644-4559, and J1645-0317 are in the Southern hemisphere, under the SKA-1 reach.For the observation of J0332+5434 in the Northern hemisphere we instead use LOFAR 2.0.In the case of the front-light echo (color dashed), we remove the contribution from observation angles smaller than the frequency-dependent size of the broadened pulsar image.This would be given by photons passing again through the pulsar, where they might be absorbed or rescattered.For simplicity, we discard those photons.The pulsar's apparent angular size is given by the broadening due to the scattering of photons in the interstellar medium and it is related to the scattering timescale τ s as θ b = τ s /x s [45].
As it can be seen from Figure 11, in this case, prospects fall into a strongly constrained region of the ALP parameter space, already excluded by the CAST helioscope results.We find that the collinear emission provides the best sensitivity.
On top of the forecast for a single pulsar, we then estimate the projected sensitivity in the case of a survey of all Galactic pulsars.Such survey is planned for SKAO [54] with an already ongoing program along this direction conducted with MeerKAT [55].Providing an accurate sensitivity forecast for the survey case is beyond the goal of this paper since it would require knowing the schedule of observations and performing simulations of pulsar populations.We can however derive a simple but reliable projected sensitivity in the following way.Let's assume Vela (J0835-4510) to be a typical pulsar, in terms of luminosity and other properties, and to be the closest one, which also means the brightest one in terms of flux.We consider a radio luminosity distribution N ∝ L −0.8 for Galactic pulsars, as found in [56].Let us now take a conservative simple 1D distribution along the direction to the Galactic center, or, say in other words, let us conservatively assume N ∝ S (0)−0.8ν (in reality it will be somewhat steeper).This implies that the collective pulsar flux measured in a survey is dS S dN/dS ≃ 4 S Vela .Then we take an average pulsar distance of ⟨d⟩ ≃ 8 kpc and an average DM density between us and the pulsar of ≃ 2.9 ρ ⊙ , obtained volume averaging the Navarro-Frenk-White profile used in Section 3 above in the inner 8 kpc from the Galactic center.Finally, we assume 10 hr of observations for each survey field (which, in general, will contain several pulsars), instead of 100 hr considered for the single pulsar case discussed above.Combining all these arguments we obtain an improvement of a factor of ∼ 10.2 in the projected g aγ bound with respect to the case considering the Vela pulsar only.
In Figure 11, we show the forecast for a survey with SKA-1 assuming this simple rescaling factor (black dashed).We focus on the collinear signal since we saw from the single pulsar case that it is more promising than the front-light echo.
Our ultimate forecast is for a survey with SKA-2 (black solid).With such a strategy, we will be able to touch the more interesting region of the ALPs parameter space, providing complementary bounds to haloscopes.A detailed forecast making use of pulsars' population models and more precise survey details is beyond the scope of the present work.

Conclusions
In this work, we present the in-depth formalism for computing the signals associated with photon emission from the stimulated decay of axions in the Galactic halo.We started from first principles and focused on stimulating point-like sources in our Galaxy.The presence of µeV ALPs in the Milky Way halo would imply the existence of three distinct line signatures: 1. Front-light echo: Stationary emission from approximately the same direction of the stimulating source and with a slightly broadened size and an orthogonal polarization.
2. Back-light echo: Stationary emission from the opposite direction with respect to the stimulating source and with a significantly broadened size and an orthogonal polarization.

Collinear emission:
Possibly pulsed emission (if the stimulating source is pulsed) with the same direction, size, and polarization of the stimulating source.
While only the back-light echo has been discussed in the literature so far, we here highlight the co-existence of these three line signatures with the above properties. 2It is difficult to imagine an astrophysical source providing simultaneously all of them.Therefore, they offer a unique opportunity for a clear ALP identification.
Previous works [11][12][13][14][15][16][18][19][20] have discussed similar or connected signals, and here we aim at providing a comprehensive view, and an in-depth discussion of the rich and diverse phe-nomenology that arises from different targets.In particular, we discussed how the echoes and the collinear emission depend on the ALP DM properties, namely on the DM spatial and velocity distributions, and on the source characteristics, such as distance, age, and motion.These dependencies offer extra handles for the ALP signal identification.
In the final part of the paper, we focused on the case of Galactic pulsars as stimulating sources.We showed that pulsar radio surveys of the next decade can be used to identify the described signatures for µeV ALPs.
To compute the echo flux observed from Earth today, we also need expressions for the distance Earth-source and source-decay location (B.8) Notice x s = xs (−x s ) and x ds = xs (−x s ), when choosing t obs = 0.
The echo flux is related to the phase space density of photons at t em , at the location where the source was back then, as follows where r * is the radius of the source.Notice that in f γ (t, x), x is the distance from the location where the source was at time t.The Doppler shift of the source's radiation as seen from the decay location compared to as seen from Earth is unimportant to the echo if f γ (ω k , x ds ) changes slowly as a function of ω k .

C Implementation
To evaluate Eq. (2.12) numerically, we introduce two cartesian coordinate systems: 1) XY Z centered at the Galactic center, the Z axis pointing towards the Galactic north pole (GNP) and the X axis pointing towards the Sun's location; 2) X ′ Y ′ Z ′ obtained by translating XY Z along the X axis by r ⊙ , where r ⊙ is the distance of the Sun from the Galactic center.The system X ′ Y ′ Z ′ is centered at the Sun's location (we neglect the vertical displacement of the Sun from the Galactic plane).A point with Galactic coordinates (l, b) at a distance d from the Sun has cartesian coordinates The galactocentric distance of the point r is then given by To evaluate Eq. (2.12), we need to compute r for all points along the line of sight.So, we need to find an equation for the los in the XY Z coordinate system.Let the source have coordinates (l s , b s ) and distance from the Sun x s .Then, for the front-light echo, we will center our observations in the direction (l f s , b f s ) = (l s , b s ), while for the back-light echo we'll have (l b s , b b s ) = (l s + π, −b s ).The los is displaced on the sky from the direction (l j s , b j s ) (with j = f, b) by an angle θ i > 0 in some direction ϕ i .Let's define ϕ i = 0 when the los is displaced with respect to (l i s , b i s ) in the direction of growing Galactic latitude and unchanged Galactic longitude, and let ϕ i grow in the clockwise direction when we look at (l i s , b i s ) from the Sun's location.The direction of the los is then (l los , b los ) = (l i s − θ i sin ϕ i , b i s + θ i cos ϕ i ).Let X ′ d , Y ′ d , Z ′ d be the coordinates of a point along the los at a distance x s from the Sun, then the equations for the los are and the distance from the Sun of a point along the los is x d = x s t and its galactocentric distance is The next step is to calculate the components of the average DM velocity along the axes of the system xyz defined in Figures 1 and 2. To do so we express the unit vectors x, ŷ, ẑ in the X ′ Y ′ Z ′ system.They are given by From these, we can calculate the components of the average ALP velocity in the xyz system and then obtain ⟨v ∥ ⟩ and ⟨v ⊥ ⟩ through the rotation matrix Eq. (A.13).

D Collinear emission
Starting from the Boltzmann equation for the collinear emission Eq. (4.1) we proceed to integrate over the momentum of the echo photon traveling back toward the source (D.2) We have chosen the direction of ⃗ q to be that of the positive z-axis and as usual we keep corrections of order p only when they are relevant to the frequency or direction of the photons.Plugging in the expression for the matrix element  As expected, the transverse average axion momentum has no effect, and the decay rate is the largest for ϵ = ⟨p ∥ ⟩.

Figure 1 :
Figure 1: Geometry of the back-light echo.

Figure 2 :
Figure 2: Geometry of the front-light echo.

Figure 3 :
Figure 3: Angular distribution of the echo from decays happening at fixed distance x d in the idealized case.We set δv = 10 −3 .

Figure 4 :
Figure 4: Emission profile of the echoes after integrating along the los in the idealized case.We set δv = 10 −3 .

Figure 5 :
Figure 5: Effects of the dark matter halo on the back-light echo.

Figure 6 :
Figure 6: Effects of the dark matter halo on the front-light echo.

Figure 7 :
Figure 7: Effects of the source characteristics on the back-light echo.

Figure 8 :
Figure 8: Effects of the source characteristics on the front-light echo.

Figure 9 :
Figure 9: Bandwidth-averaged front-light echo spectral intensity at 400 MHz for four of the most promising pulsars from the ATNF catalog assuming g aγ = 10 −11 GeV −1 .Blue arrows indicate the direction of the component of the average DM velocity perpendicular to the line of sight, while magenta arrows indicate the source's proper motion.The arrows are drawn in arbitrary units, but their relative lengths are to scale with each other.The thin black lines mark contours of constant ⟨I ν ⟩.

Figure 10 :
Figure 10: Bandwidth-averaged front-light echo spectral intensity at 400 MHz for four of the same pulsars as Figure 9 assuming g aγ = 10 −11 GeV −1 .On the left (right) panel, we vary the Galactic latitude (longitude) of the line of sight, keeping the Galactic longitude (latitude) fixed.

X
′ = −d cos b cos l Y ′ = −d cos b sin l Z ′ = d sin b .(C.1) x = (sin b s cos l s cos ϕ i ± sin l s sin ϕ i , sin b s sin l s cos ϕ i ∓ cos l s sin ϕ i , cos b s cos ϕ i ) ŷ = (− sin b s cos l s sin ϕ i ± sin l s cos ϕ i , − sin b s sin l s sin ϕ i ∓ cos l s cos ϕ i , − cos b s sin ϕ i ) ẑ = ∓(− cos b s cos l s , − cos b s sin l s , sin b s ) .(C.5)

Table 1 :
Properties of hypothetical sources.The columns are the source Galactic longitude l, Galactic latitude b, distance from Earth as seen today x s , age τ age and direction of the source's velocity in Galactic coordinates l ⃗ vs and b ⃗ vs .For all sources the magnitude of ⃗ v s is 200 km/s.
i ) For sources C and D, whose lines are on top of each other, we have |⟨⃗ v ⟩| ≈ |⟨v t ⟩|, resulting in an overall reduction of the flux by a factor exp −⟨v t ⟩ 2 /2δv 2 .