The axion dark matter echo: a detailed analysis

It was recently shown that a powerful beam of radio/microwave radiation sent out to space can produce detectable back-scattering via the stimulated decay of ambient axion dark matter. This echo is a faint and narrow signal centered at an angular frequency close to half the axion mass. In this article, we provide a detailed analytical and numerical analysis of this signal, considering the effects of the axion velocity distribution as well as the outgoing beam shape. In agreement with the original proposal, we find that the divergence of the outgoing beam does not affect the echo signal, which is only constrained by the axion velocity distribution. Moreover, our findings are relevant for the optimization of the experimental parameters in order to attain maximal signal to noise or minimal energy consumption.


I. INTRODUCTION
The identity of dark matter is one of the biggest puzzles in modern science. The QCD axion [1][2][3], as well as axion-like particles [4], are leading candidates for its fundamental composition [5][6][7][8]. Under this assumption, several experimental efforts are being carried out. See references [9,10] for recent reviews. Most of these experiments are based on the axion to two photons vertex [11] described by the interaction term gaF µνF µν , where a is the axion field, F µν the photon field strength tensor and g the coupling constant. One significant feature of this interaction is that the axion decay rate can be dramatically enhanced by the presence of a background of photons with energy close to half the axion mass. This axion stimulated decay is a Bose-enhancement effect, which is intrinsical of Bose statistics. The topic has been thoroughly discussed in the recent years [12][13][14][15][16][17][18]. One peculiarity of the stimulated decay is that the produced photons are constrained to propagate along the incident photon path in the axion rest frame. From the decay of a single axion, one produced photon adds to the incident one, while the other is ejected in the opposite direction. This backward radiation was studied in Ref. [19] and baptized as the "axion dark matter echo". Its physical effects were estimated, getting promising results for axion dark matter searches. To be more specific, the authors propose to send out to space a powerful beam of radio or microwave radiation and to detect the echo at a spot located nearby the emitter using current radio astronomy technology. A very similar idea was recently implemented in the case of extragalactic radio sources [20].
In [19], it was claimed that, in the axion rest frame, the echo wave has the same spatial shape as the outgoing beam. The important consequence is that in such a case, the echo returns exactly at the emission spot, implying that the echo signal does not depend on the beam shape. On the other hand, a realistic velocity distribution of the local axion flow makes the echo spread in space, weakening any signal collected by limited size detectors. Despite this, it was found that the signal is strong enough to pursue this idea.
This paper presents an in-depth analysis of the echo signal using the Green's function method. This approach allows a detailed joint examination of the role of the outgoing beam shape and the axion velocity distribution. As claimed in [19], we confirm that a divergent shape of the outgoing beam does not affect the echo signal at all. Concerning the effects of the local axion velocity distribution, our findings agree with those estimated in [19] as well.
This article is structured as follows. In Section II, we first compute the total energy and power of the echo in general terms. Then, we highlight the role of the local axion velocity distribution by a naive estimation of the echo intensity at the emission spot. To do so, we use the particular case of a beam emitted by a parabolic antenna in the far-field zone limit, where analytical expressions are known. It also allows us to confirm the null effect from the divergent shape of the beam. The subsequent Sections are devoted to the computation of the echo intensity as a function of space-time coordinates. In Section III, we use a one-dimensional framework as an illustrative example, while in section IV we establish the mechanism to address the 3-dimensional analysis using the Green's function method in the paraxial limit. In Section V, we put this machinery at work assuming a paraxial Gaussian beam, similar to that of a laser, which is an excellent approximation of the one produced by a parabolic antenna in its far-field zone limit. We give an analytical description of the effects of the axion velocity distribution on the echo intensity. In Section VI, we perform numerical computations to get the intensity map and the power spectrum of the echo in the isothermal sphere and caustic ring models for the dark matter halo. Finally, in Section VII, we estimate the expected sensitivity of the proposed experiment, while in Section VIII we summarize and conclude.

II. MAIN FEATURES OF THE ECHO SIGNAL
The first step of our discussion is to highlight, in a simple manner, the main properties of the echo signal. We will first compute the energy and power stored in the echo wave. Although the method used here has the advantage of not requiring the knowledge of the explicit form of the outgoing beam, it only provides formulas for the total energy and power of the echo, saying nothing about how these are distributed in space.
In the second part, we will partially fill the deficiencies of our first method by describing the behavior of the echo signal at the emission spot x 0. To do so, we take, for the outgoing beam, the explicit form of a beam emitted by a parabolic dish antenna. Then, assuming that most of the echo comes from the far-field zone of the beam, we calculate the echo field and intensity for x → 0, using the Green's function method. Although this analysis is limited, it gives us some idea of what the intensity looks like. In particular, it explicitly shows the role of the axion velocity distribution and reveals the null effect from the beam divergence.
We would also like to remark that even though this Section only covers a rough characterization of the echo, much of the notation defined here will be used in the subsequent Sections.
A. Energy and power of the echo The differential equation for the electromagnetic vector potential A(x) with an axion background field a(x) is given by where the Coulomb gauge ∇ · A = 0 has been chosen. We have also neglected terms containing gradients of a. This is justified by the assumption of a non-relativistic axion background, meaning that the axion momentum p is much smaller than its mass m. From now on we will always ignore this kind of contributions. For the time being, let's just consider one axion momentum mode with energy density ρ( p ). The axion field can be written as where E p = m 2 + p 2 . As E p = m + O(p 2 ), we will use E p ≈ m throughout this article. Let's study the first order perturbative correction that the axion background produces over an incident electromagnetic wave. If the incident wave has a vector potential A (0) , the first order correction A (1) is determined by Let us expand the source field A (0) as well as the correction A (1) in Fourier modes as We omit the complex conjugate symbol c.c. for every expression of the axion and electromagnetic fields. Its contribution will be assumed by default. Above we have assumed an incident field linearly polarized in the directionê. A (1) is polarized necessarily in the directionk ×ê, except for small corrections of the order of p. We have also defined ω(k) ≡ | k|. Keeping only terms relevant to the stimulated axion decay, i.e. the photon momentum modes k and q = p − k, Eq. (2.3) becomes By looking at the source term of the above equation, we see that when m > ω(k), the solution describes a wave which travels backwards with respect to the incident one. We call it the echo wave. The echo wave is excited when ω(q) is equal or very close to m − ω(k), i.e. when ω(k) ≈ m/2. Now we look for a resonant solution using the ansatz A (1) (t, q ) = A(t, q )e −iω(q) t where A(t, q ) varies slowly in time with respect to A (1) (t, q ). Neglecting second derivatives in Eq. (2.6), we get Solving (2.7) with the initial condition A(0, q ) = 0, we obtain In the limit t 1, we can substitute sin( t/2)/ → πδ( ) in Eq. (2.9). Integrating over all axion momentum modes, we obtain the energy of the echo as where U 0 is the energy of the outgoing beam. In our non-relativistic limit, ≈ 2ω(k)−m−p , where p is the component of p in the direction of k. Assuming that the incident beam propagates mostly in a preferred direction, the direction of k and p are fixed, so U = π 4 g 2 t dp ∂ρ ∂p dω ∂U 0 ∂ω δ( (ω, p )) . (2.11) Let's perform the integral (2.11) in two scenarios: when the bandwidth of the beam δω is much bigger than the bandwidth of the axion δp and in the opposite case. For δω δp, the axion distribution can be considered as if it were condensed at a single valuep . We then approximate ∂ρ/∂p = ρ δ(p −p ). We straightforwardly find where ω * = m +p 2 . (2.13) The echo energy grows linearly in time. This follows from the assumption of an infinite extension of the axion background. The emitted beam travels in the axion background extracting energy at a constant rate. If the energy U 0 is provided from a source emitting a constant power P 0 from t = 0 until t = t off , the energy of the beam is U 0 = P 0 t off and the power of the echo is (2.14) For δω δp, we have ∂U 0 /∂ω = U 0 δ(ω −ω). We find as well as where For the DFSZ model of the QCD axion, an input power of 1 kW over a bandwidth of 1 kHz working for one hour 1 , provides a total echo power of P ∼ 10 −18 W m 10 −4 eV 2 , which is in principle not too difficult to detect with current radio astronomy technology. Unfortunately, as discussed in Ref. [19], due to the non-trivial axion velocity distribution, this power spreads over a surface that eventually exceeds the detector's size. This effect is easy to visualize if we consider what happens for a single axion decay. If the decaying axion moves with a velocity not perfectly aligned with the incident photon (outgoing beam), the echo photon is released in a direction also different from the incident photon trajectory. In fact, if the echo photon has energy Ω, momentum conservation implies Ω sin(χ) = mv ⊥ , where χ is the angle formed by the trajectories of both photons and v ⊥ is the component of the axion velocity perpendicular to the incident photon momentum. From energy conservation, Ω is equal to m/2 plus small corrections, therefore sin(χ) 2v ⊥ . As the axion velocities are of the order of 10 −3 or smaller, we can write χ 2v ⊥ . It follows that, for instance assuming the isothermal sphere model for the galactic halo with velocity dispersion of 270 km/s, the echo spreads over approximately 10 6 km after one hour. Even in the caustic ring model of the galactic halo, where the minimal axion transverse velocity is about 5 km/s, and therefore the echo method is more promising, the spread is at least 10 4 km after the same amount of time. The echo signal is thus strongly dependent on the model for the axion phase-space distribution as well as the size of the detection apparatus.

B. Echo of a dish antenna beam
To develop a sense of how the transverse axion velocities affect the signal, we will perform a simple and useful computation of the echo intensity for some particular models of the axion velocity distribution. Moreover, as an incident wave, we will take the beam emitted by a parabolic dish antenna. Although the parabolic antenna case is the most realistic setup we can address in this work, our findings are limited by scenarios where analytical approximations are possible. These approximations will be explained throughout the text.
Given the axion field a(t, x) and outgoing beam magnetic field B (0) (t, x), the echo vector potential field A (1) (t, x) is determined by For instance, using a klystron. 6 We write the axion dark matter field, in a large volume V , as an expansion in momentum modes as where p is the axion mometum in the experiment rest frame and φ p random phases. These phases are not known, and are usually modeled as uniformly distributed random numbers [21][22][23]. In this article, we limit ourselves to considering the ensemble average where ens indicates the average over a large number of draws of the random phases. We leave a discussion of the statistics of the signal to future work. We normalize f a ( p ) as such that the energy density in the volume V , averaged over time and random phases is We consider a parabolic antenna of radius R fed by a plane wave with electric field amplitude E 0 (ω), ω being the wave's frequency. The electric field emitted by the antenna is known analytically in the far-field zone, i.e., for distances of order ωR 2 or larger from the emission spot. Assuming that the center of the antenna is located at x = 0 and that the central component of the beam propagates alongẑ, the antenna's electric field can be written in spherical coordinates as where J 1 (x) is the Bessel function of the first kind of order one. We have taken the electric field to be linearly polarized alongx. For Eq. (2.18) to be a good approximation, we need most of the echo to be produced in the far-field zone. Therefore, for an outgoing beam turned on at t = 0, our approximation is valid for t ωR 2 . We perform integral (2.18) in the limit ωR 1, which is valid for any experimental setup we are considering in this article. We can then make two additional approximations. First, as most of the echo comes from distances of the order or larger than ωR 2 , we can assume | x | R. The field at the emission spot can then be found using the limit x → 0. Second, the beamwidth of the outgoing wave is of the order (ωR) −1 , allowing us to approximate θ 0. If the outgoing beam is emitted continuously from t = 0, the echo vector potential at the emission spot reduces to where we have defined andρ(ϕ ) =x sin(ϕ ) +ŷ cos(ϕ ). After integration, we get where u = t − t and ξ = min(t, mR/p ⊥ ). To obtain the result above we have extended the upper integration limit dθ to infinity. This is justified for ωR 1 and | p ⊥ |r 1, as then the largest contribution to the integral comes for θ π/2. To find the spectral intensity, we average the quantity Ω 2 A (1) (t, 0, Ω) 2 over fast time oscillations and over the random phases φ p . Assuming also that the axion momentum distribution can be separated in a forward and transverse component as 2 we have Consistently with Eq. (2.21), we normalize f a (p z ) and f a (p ⊥ ) as dp z |f a (p z )| 2 = 1 (2.31) and respectively. From now on, we assume for f a (p z ) the Gaussian shape (2.33) In our notation, stands for averaging over the distribution |f a ( p )| 2 , i.e. the expectation value for every functionQ( p ) is Thus p z is the average value for p z and δp z = p 2 z − p z 2 the dispersion. Notice that, with the normalization Eq. (2.21), ifQ does not depend on one of the momentum components, the corresponding integration trivially yields 1.
At resonance, i.e. when ω = (m + p z )/2 ≡ ω * , we have dp z |f a (p z ) where v ⊥ = p ⊥ /m is the axion transverse velocity and Ω * = (m − p z )/2. We can see from the above result that the echo signal depends strongly on the local axion velocity distribution. Indeed, it scales as the inverse of the momentum dispersion in the forward direction and as the inverse of the average transverse velocity component. In this paper, we consider two models for the velocity distribution in the Milky Way ( v = p/m), the isothermal sphere [24], and the caustic ring model [25]. In the isothermal model, |f a ( v )| 2 has a Maxwell-Boltzmann behavior with dispersion δv = 270 km/s and average velocity given by the velocity of the Sun in the galactic rest frame, i.e. | v | 230 km/s. In the caustic ring model, the local dark matter is dominated by a single cold flow with velocity dispersion δv < 70 m/s and with an average velocity of | v | 290 km/s relative to us. More details on both models will be given in Sec. VI. To include both halo models in our analysis, we write the transverse velocity distribution as where for the isothermal model we make δv ⊥ = 2/3 δv in order to get the Maxwell-Boltzmann distribution for |f a ( v)| 2 . For the caustic ring model, as δv ⊥ is usually much smaller than v p , we take the limit δv ⊥ → 0, getting a delta function centered at v ⊥ = v p , as a reasonable approximation. Of course, this approximation breaks down if the outgoing beam points to a direction parallel to the axion flow (or to a direction such that δv ⊥ ≥ v p ), however the direction of the big flow is determined by the position of the IRAS [26], Planck [27] and GAIA [28] triangles with an uncertainty of 0.01 radians. It implies that v p cannot be reduced to values smaller than 5 km/s, clearly much larger than δv ⊥ .
To compute the quantity v −1 ⊥ that appears in Eq. (2.36), we use Eq. (2.37). We explicitly get where I 0 (x) is the modified Bessel function of the first kind of order 0. For v p δv ⊥ , compatible with the caustic ring model, δv ⊥ , which is compatible with the isothermal model when the outgoing beam is pointed mainly in the direction of the axion wind, we have v −1 3 It is true for almost all the axion masses we are considering in this work In this analysis, we have seen how the axion velocity distribution influences the echo signal. Eq. (2.36) shows the effects explicitly. For t < R v −1 ⊥ the echo has not yet spread beyond the emission region, leading to linear growth of the intensity, in agreement to the results of Sec. II A. After t = R v −1 ⊥ , the intensity saturates suddenly to a constant value, due to the transverse velocity effects. This abrupt change in the intensity behavior is a consequence of the approximations employed, the assumption that most of the echo is coming from the far-field zone. If we had also taken into account the near behavior of the outgoing beam, the transition to the saturated regime would have been smooth. Another limitation of our result is that Eq. (2.36) gives us an approximated value for locations nearby the emission spot. A complete characterization of the outgoing beam would allow us to know the local values of the intensity in more detail. Unfortunately, there is no simple analytical expression for the antenna emitted field in the near zone. We then leave this task for future work.
To gain further insights into the intensity as a function of the position and into the transition to the saturated regime, we will take a simplified model for the outgoing beam. The model is based on the paraxial Gaussian beam used in laser physics. This beam features a simple shape in the near-field zone and matches the antenna beam's far-field zone behavior correctly. Given this beam, we will use the paraxial approximation to get the local values of the echo intensity as well as its time evolution. To warm up engines, in Section III, we will show a complete one-dimensional calculation of the echo signal, then in Sections IV and V, we will take care of the three-dimensional analysis.

III. ONE-DIMENSIONAL ANALYSIS
In order to gain a better understanding of how to calculate the echo intensity, in this Section, we are going to perform a one-dimensional analysis before moving on to the more involved three-dimensional case in the next Section. Even though it is just an illustrative example, as in the previous Section, some of the notation defined here will be used throughout the rest of the paper. First, we chooseẑ as the propagation direction of the outgoing beam and set z = 0 at the emission spot. We also assume that the outgoing beam is linearly polarized in thex-direction, which implies that the echo is linearly polarized in thê y-direction. In this setup, Eq. (2.3) can be written as where B (0) is the magnetic field of the outgoing beam. We assume that the outgoing beam is turned on within the time interval 0 < t < t off and, for simplicity, we assume that during this period it is emitted with a constant amplitude 4 . This means that the RHS of Eq.
for t > t off .
We write the axion field in terms of its Fourier expansion as where f a (p z ) is normalized as in Eq. (2.31). The axion energy density in the length L, averaged over time and random phases 5 , is related to a 0 as in Eq. (2.22). The magnetic field of the incident beam is a plane wave given by where the field amplitude B 0 is related to the beam intensity by For every ω, the echo wave is sourced by a term that has frequency Ω, so we write the echo vector potential as As we look for resonant solutions, we assume |∂ z A| ΩA and |∂ t A| ΩA. Neglecting second derivatives of A, we get The general solution of Eq. (3.8) is where G(t, z) is the retarded Green's function associated to the differential operator ∂ t − ∂ z . The explicit calculation of G(t, z) is found in Appendix A. Using Eq. (A5), we have where the integration boundaries are The echo spectral intensity, at z = 0, averaged over fast oscillations and random phases is where t < is the minimum between t and t off .

IV. THREE-DIMENSIONAL ANALYSIS
For our three dimensional analysis, we model the outgoing beam as a transverse wave propagating in the z direction with small transverse corrections. We write the beam magnetic field as whereˆ is a constant polarization vector. The function η( x, ω) defines the spatial shape of the beam. It is normalized as η( 0, ω) = 1, such that B 0 (ω) is related to the beam time-averaged intensity I 0 (ω), at x = 0, by The beam is turned on at t = 0. At a time t > 0, it extends from z = 0 to z = t. After the beam is turned off at t = t off , it extends from z = t − t off to z = t. In other words, the integration domain in the z direction is given by Eq. (3.11), while for the transverse directions, it is determined by η( x, ω).
We describe the local dark matter axion field as a superposition of plane waves, in the same way as in Section II B (see Eqs. (2.19)-(2.22)). In addition, we also assume that f a ( p ) can be factorized in its forward and transverse parts (see Eqs. (2.29), (2.31) and (2.32)). Discarding non-resonant terms, Eq. (2.3) reads where x ⊥ = (x, y, 0), while Ω and are defined in Eqs. (2.25) and (2.26), respectively. We now write the echo field as Notice that Ω can be written as Ω = ω −p z − . As at resonance → 0, the frequencies of the echo are shifted with respect to ω by p z . As we will discuss in Section VI, this effect is relevant for the frequency spectrum of the echo because, although p z is small compared to ω, it may cause a complete separation between the outgoing beam's and echo's frequency ranges. It might have important consequences from the experimental point of view. However, this shift does not affect the amplitude of the echo. Now we plug Eq. (4.4) into Eq. (4.3) and use the fact that the outgoing beam as well as the echo field satisfy the paraxial approximation with i = t, z. For the echo, we will check in Subsection V C that our results are consistent with the paraxial assumption. In this limit, Eq. (4.3) reduces to the paraxial equation with a source term The solution is where As derived in Appendix B, the Green's function for the paraxial equation is The echo spectral intensity averaged over fast time oscillations is given by Although Eq. (4.10) is the most general formula of this work, its value at z = 0 can be reduced notably in the limit δp z t 1. Since our goal is the estimation of the echo amplitude, for which corrections O(p) are irrelevant, we use the approximation Ω = ω − p z − ≈ ω for the following analytical computations. With the axion momentum distribution Eq. (2.33), and after averaging over random phases, the spectral echo intensity, at resonance, can be written as where and The integration domain for z is specified in Eq. (3.11), evaluated at z = 0. For a derivation of Eqs. (4.11) and (4.12) see Appendix C. Notice that Eq. (4.11) is the three dimensional analog of Eq. (3.13). It means that all the 3D effects enter through in T (t, x ⊥ ). Finally, to compute the total power of the signal, we should know the spectral intensity of the outgoing beam in order to integrate over all frequencies that might contribute. We will use for simplicity a Gaussian spectrum of the form where I 0 is the total intensity,ω the central value of the distribution, and δω the dispersion. Thus, the power collected over a surface S (see also Appendix C) can be computed as where S 0 is the effective cross sectional area of the beam at z = 0, P 0 = I 0 S 0 the power of the beam and 13

V. THE ROLE OF THE OUTGOING BEAM'S SHAPE AND THE AXION VE-LOCITY DISTRIBUTION
In this Section, we are going to use all the machinery developed in Section IV to analyze 3D effects assuming a particular model for the outgoing beam shape. We take our model incident beam as the Gaussian beam obtained in laser physics by solving Maxwell equations in the paraxial limit. In cylindrical coordinates, (r, φ, z) it is given by Here w(z) is the radius at which the beam intensity is reduced by 1/e compared to its axial value and R(z) is the radius of curvature of the beam wavefronts. Their mathematical expressions are where R = w(0) and z R = ωR 2 . In this model, the beam decays radially as a Gaussian function and also develops a divergence that becomes noticeable for z > z R . Moreover, in the limit ωR 1 (always satisfied for the parameters' values relevant to this work), we recover the behavior of the dish antenna beam Eq. (2.23) for z ≥ z R . It is worth mentioning that, in Eq. (5.1), we have ignored the Guoy phase ψ(z) = arctan(z/z R ), which for our purposes does not contribute in any way.
The time-averaged power of the beam is P 0 = I 0 πR 2 , then the surface parameter entering Eq. (4.15) must be given by S 0 = πR 2 . Thus, we can also define R as the effective radius of the beam at z = 0.
The purpose of this Section is to provide analytical expressions for T (t, x ⊥ ) defined in Eq. (4.12). According to Eq. (4.11), T contains all the properties of the echo intensity that derive from the transverse axion velocity distribution and the shape of the beam. Therefore, the behavior of T (t, x ⊥ ) is equivalent to the behavior of the echo intensity.
We start our analysis, by computing ξ from Eq. (4.13) using the beam model Eq. (5.1). We get Using Eq. (5.2), it seems that Eq. (4.12) could be very difficult to integrate analytically in both dz and d 2 p ⊥ , however it is manageable in some regimes. Let's analyze these regimes and see what we can learn from them. To do so, we first define the function T (t, x ⊥ , v ⊥ ) such that T (t, x ⊥ ) = T (t, x ⊥ , v ⊥ ) . From Eq. (4.12), it is found to be where |f a ( v ⊥ )| 2 is given in Eq. (2.37). In the following Subsections, we will do so and reveal the transverse velocity effects that are not present in the one-dimensional case. The transverse velocity components cause a spread of the echo in those directions. Not all the power can be collected at the emission spot as in the 1D case. The lateral spread leads to saturation of the intensity on a time scale that we will determine below. On the other hand, as also shown in Section II B, we find null effects coming from the outgoing beam divergence. Our analysis will be split into two cases; small dispersion (δv ⊥ v p ) and large dispersion (δv ⊥ v p ). Before discussing the 3D effects in detail, as a consistency check, notice that when the transverse velocities are not important, i.e. v ⊥ → 0, we get from Eq. (5.3) that T (t, 0) = t for t < t off and T (t, 0) = t off for t > t off . In other words, we recover the result from the one dimensional analysis of Section III.

A. Small dispersion
, then the integration in Eq. (5.4) becomes trivial. We have where φ is the angle formed by x ⊥ and v p . We will analyze Eq. (5.5) in two scenarios; a) t off < R/v p and b) t off > R/v p , each characterized by two regimes; t < t off and t > t off .
In scenario a), we can Taylor expand T (t, x ⊥ ) in terms of the small quantities v p t/R for t < t off , and v p t off /R for t > t off . We find for t > t off .
We see that for t < t off , T (t, x ⊥ ), and therefore the intensity, grows linearly in time, featuring a transverse Gaussian shape with effective radius R. For t > t off , i.e., after the outgoing beam is turned off, the Gaussian keeps its maximum value reached at t = t off and moves rigidly with velocity v p .
In scenario b), T (t, x ⊥ ) behaves exactly as the upper formula in Eq. (5.6) for t < R/v p . Later, for R/v p < t < t off , we observe from Eq. (5.5) that the echo spreads along the direction of v p at a speed v p , forming a "sausage" shaped intensity profile at z = 0. In this regime, we can find the maximum possible intensity compatible with the small dispersion limit. The energy that in scenario a) accumulates at the emission spot, now spreads laterally, leading the echo intensity to saturate at any coordinate x ⊥ when t > r cos(φ) vp . The maximum saturated value is found for locations where φ = 0 and R < r cos(φ) < v p t. The corresponding value for T within this region is twice the value obtained at r = 0. We name it Finally, for t > t off , the "sausage" shape moves rigidly at the speed v p having a length v p t off and a width 2R.

B. Large dispersion
Let's discuss now the case δv ⊥ v p . As in the previous discussion, we split the analysis in two scenarios; a) t off < R/δv ⊥ and b) t off > R/δv ⊥ , and each of them in two time regimes; t < t off and t > t off .
In scenario a), we Taylor expand T (t, x ⊥ , v ⊥ ), given in Eq. (5.3), in terms of the small values v ⊥ t/R for t < t off and v ⊥ t off /R for t > t off . The result is, of course, identical to Eq. (5.6) after performing the substitutions v p → v ⊥ and φ → ϕ. Now, integrating over the velocity distribution, we straightforwardly obtain the upper formula in Eq. (5.6) for t < t off . For t > t off , the integration yields This behavior can be separated into two sub-regimes; t off < t < R/δv ⊥ and t > R/δv ⊥ . In the first one, the intensity maximum stays constant at the value obtained at t = t off , and also, its radial extension R does not change. In the second one, the radial extension becomes δv ⊥ t, i.e., the energy spreads radially at velocity δv ⊥ . On the other hand, the value of the intensity decreases as (δv ⊥ t) −2 . Scenario a) can be summarized as follows for t > R/δv ⊥ .

(5.9)
In scenario b), T (t, x ⊥ ) behaves exactly as the upper formula of Eq. (5.6) for t < R/δv ⊥ . After that, for R/δv ⊥ < t < t off , the energy spreads radially leading the intensity to saturate for all the positions satisfying r < δv ⊥ t. The corresponding saturated value for T (t, x ⊥ ) can be found taking the limit erf(x) → 1 for the second erf function appearing in the upper formula of Eq. (5.3). Thus, Eq. (5.4) gives us The maximum value of T in this saturated regime is found at r = 0 and is We can also find analytically how the system transits to this saturated value. Evaluating Eq. (5.3) at r = 0, Eq. (5.4) gives With the same procedure, we find for t > t off which in the limit t t off becomes T (t, 0) = 1 2π (5.14) As a summary of scenario b), we can write for R/δv ⊥ < t < t off 1 2π for t > t off and r = 0 .

(5.15)
As a general remark, we point out that t a and t b are not only the saturated values of T in the respective scenarios but also the timescales over which the saturation is achieved. We can then define the saturation time To conclude this discussion, we notice that, as long asωR 2 R v −1 ⊥ , saturation occurs before the outgoing beam reaches the far-field zone. This is explained by the fact that the intensity saturates for t R/v p in the case of large average transverse velocity and for t R/δv ⊥ in the case of large transverse velocity dispersion.

C. Consistency of the paraxial approximation
To conclude this Section, we discuss the consistency of our results with the initial assumption of the paraxial limit. The paraxial approximation works very accurately when considering beams with divergences not larger than 1 rad [29]. For the echo wave, the divergence angle χ is related to the velocity v s , at which the echo spreads in the transverse direction, by χ = arctan(2v s ) .
On the other hand, the intensity saturates when the echo wave leaves the cross-sectional area of the outgoing beam, i.e., v s t sat ∼ R. It implies that v s is of the order of v p or δv ⊥ , depending on which parameter dominates the saturation. For any cold dark matter model, v p and δv ⊥ should not exceed 10 −3 . Hence the transversal velocity effects only contribute with small divergences, compatible with the paraxial approximation.

VI. INTENSITY MAP FOR THE CAUSTIC RING AND ISOTHERMAL MOD-ELS
In this Section, we compute the intensity of the echo on the z = 0 plane numerically. We consider two models for the local dark matter distribution: the isothermal sphere [24] and the caustic ring model [25].
We write the echo intensity as where t sat is defined in Eq. (5.16) and F (t, x ⊥ ) is a function whose maximum is of order 1 given by Here J(t, x , Ω) is given in Eq. (4.8) with the beam profile Eq. (5.1). Notice that, as follows from the discussion of Section V, in the limit δp z t 1, F has an asymptotic behavior. In particular, when the echo saturation is predominantly due to one of transverse average axion velocity or velocity dispersion, the maximum value of F tends to 1 for times larger than t sat .

A. The caustic ring model
The caustic ring model is a proposal for a complete description of the dark matter phasespace distribution in the halo of our galaxy [25]. The model is characterized by the presence of caustics in the galactic plane. These caustics have the shape of closed circular tubes, approximately centered at the galactic center, and whose cross-section is a section of the elliptic umbilic (D −4 ) catastrophe. Caustics of cold dark matter generically form after the growth of structure becomes non-linear. Cold, collisionless dark matter lies in sixdimensional phase space on a thin three-dimensional hypersurface, the "phase space sheet", the thickness of which is the primordial velocity dispersion of the dark matter particles. As dark matter particles fall in and out of a galactic gravitational potential well, their phase space sheet wraps up. Locations where the phase space sheet folds back onto itself are caustics. If the velocity field of the infalling particles is dominated by net overall rotation ( ∇ × v = 0), the caustics have the aforementioned tube shape [30].
Based on the location of two triangular features present in the IRAS [26], Planck [27] and Gaia [28] sky maps, it was determined that Earth is at a location close to a caustic ring, presumably inside the tube. In this case, the local dark matter velocity distribution is dominated by a single flow, called the big flow. From Table I of Ref. [28], we see that the big flow has a velocity of about 520 km/s in the galactic rest frame, implying a relative FIG. 1. Echo normalized intensity F (t, x, y = 0) in the case of saturation due to the transverse average axion velocity. All the parameters are the same as those of Fig. 4. We can clearly see the length of hot spot growing in the positive x-direction until mt off = 7 × 10 7 . After t off , the hot spot travels rigidly.
velocity with respect to us of order 290 km/s. The direction of the big flow is determined by the position of the triangles with an uncertainty of 0.01 radians. As a consequence, v p cannot be reduced to values smaller than 5 km/s. The velocity dispersion of the flow was determined in Ref. [31] to be at most of order 70 m/s. The local dark matter density in this model is 1 GeV/cm 3 or higher.
In this Section, we assume the axion phase space distribution Eqs. (2.33) and (2.37), with a velocity dispersion δv z = 70 m/s, negligible transverse velocity dispersion and a transverse average velocity v p = 5 km/s. We neglect daily and annual modulation effects, as they are negligible over the timescales considered in the plots presented here. 6 Although the results of this Section are presented for the caustic ring model, they are applicable to any situation in which a single cold flow dominates the local dark matter density.
With the chosen parameters, saturation happens due to transverse average velocity effects. Figure 4 shows the echo normalized intensity Eq. (6.2) at resonance, i.e. for ω = (m + p z )/2. With the chosen radius R = 4000/m, we have δp z t sat 1, so that the approximate formulas derived in Section V apply and the saturated value of F is approximately 1. In Fig. 4, we can clearly see the echo's "hot spot" spreading in the direction of the axion average velocity until the beam is turned off at mt off = 1.4 × 10 9 . Thereafter, the hot spot stops spreading and starts traveling rigidly. Fig 1 shows F along the line y = 0 at different times.

B. Isothermal halo model
For the isothermal halo model, we also use the axion momentum distribution Eqs. (2.33) and (2.37). In this case, the velocity dispersion is δv z = δv ⊥ / √ 2 = 270/ √ 3 km/s. The average axion velocity is non-zero as a consequence of Earth's motion with respect to the galactic rest frame and has typical value of 230 km/s, while the local dark matter density is approximately 0.45 GeV/cm 3 . Although the isothermal model does not accurately describe cold dark matter halos, we consider it here because the dark matter phase space distribution has a simple analytical form. Moreover, using a more realistic model, such as the NFW or Hernquist models, would not change our analysis, as locally, the axion momentum distribution can be approximated as a Gaussian, although the value of the parameters may change. Figure 5 shows F at z = 0 in a case of saturation due to velocity dispersion effects. The hot spot spreads out in all directions, while at the same time, the maximum intensity grows. This effect is also clearly visible in Fig. 2. Figure 3 shows F (t, 0 ) as a function of the beam frequency ω at four different times. F (t, 0 ) has been multiplied by the spectral intensity inherited from the beam dĨ 0 /dω| ω=Ω , whose maximum value has been normalized to 1. The beam spectral intensity Eq. (4.14) normalized to its maximum value is shown in black. The echo peaks atΩ = (m − p z )/2, while the beam is centered atω = (m+ p z )/2. The echo bandwidth is given by δωδp z /(2∆). It is then possible for the echo frequency band to be completely separated from the beam frequency band if Such a configuration may be advantageous to separate signal from noise if part of the latter comes from the outgoing beam itself or back-scattering in the atmosphere.  [32][33][34][35][36][37][38][39][40][41][42][43][44][45] and the grey region corresponds to bounds from the CAST experiment [46].

VII. ECHO POWER AND SENSITIVITY
In this Section, we estimate the sensitivity of the echo method assuming that the echo power is collected by a dish antenna with radius R c , located at the plane z = 0 and concentric with respect to the emitter. We leave an updated analysis taking into account the effect of the axion random phases to future work.
The signal to noise ratio is given by the Dick's radiometer equation where P c is the collected power, T n the noise temperature, B the bandwidth of the signal given roughly by B ≈ min(δω, δp z /2)/(2π) and t m the measurement time.
In a regime of constant signal power, the collected P c can be written in terms of the emitted power P 0 as Here t e is basically the minimum between t off and κ t ⊥ , where κ is a parameter that depends on the ratio R c /R and t ⊥ is defined as See Appendix D for a derivation of this result and an explicit expression for κ.
To determine t m , we have to take into account the fact that the echo spreads laterally at velocity v −1 ⊥ −1 . If t off < t ⊥ , the echo will continue to be received by the dish for times even larger than t off , and it will drop drastically for t > t ⊥ . In this case t m is given roughly by t ⊥ . On the other hand, when t off > t ⊥ , the echo spreads all over the dish before the outgoing beam is turned off. After turning off the beam, the echo intensity drops abruptly. In this case, t m is approximately t off . To summarize this reasoning, the measurement time is given roughly by the maximum between t off and t ⊥ .
To use Eqs. (7.1) and (7.2), we smooth t e , t m and B as We estimate the sensitivity by assuming a fixed amount of energy E is spent to search for the axion, covering an octave in axion mass range. For a source that shoots beam pulses with power P 0 and duration t off , the energy spent to cover a factor of two in axion mass is E ≈ mP 0 t off /(2δω). In terms of E, Eq. (7.1) is written as We see that the signal to noise ratio is maximized for t off → 0. However, as t off is limited by the uncertainty relation δωt off ≥ 1/2, Eq. (7.5) will be evaluated at t off = (2δω) −1/2 . The solid blue lines of Fig. 6 mark the parameter space where, in this approach, the echo method is sensitive. We used Eq. (7.5) assuming E = 10 MW year, s/n = 5, T n = 20 K, R = 50 m, R c = 100 m and δω = δp z /2. For the isothermal model (right panel), we assumed the outgoing beam pointing approximately in the direction of the axion average velocity, such that v p is negligible respect to δv ⊥ 220 km/s. For the caustic ring model, we took v p = 5 km/s, in agreement with the discussion in VI A. Fig. 6 shows a range of axion masses compatible with the working frequency of radio telescopes. The sensitivity shown by the solid blue lines of Fig. 6 assumes a fixed amount of consumed energy and does not specify how this energy is spent. Indeed, at this theoretical stage of our proposal, it is not yet clear how to optimize the efficiency of the energy delivery. For instance, one can restrict the time spent to cover an octave in axion mass range to a limited value t T . In such a case, the duration of the beam pulses t off is given by t off 2δωt T /m. It is not difficult to show that, with this constraint, the signal to noise is maximized for δω = m To have an idea of the order of magnitude, for R c = 100 m and t T = 1 year, δω is 844 Hz The dashed blue lines of Fig. 6 correspond to the sensitivities using a power source of P 0 = 10MW working for one year. The other parameters are the same as for the solid blue lines. Although the amount of energy spent in this time-restricted approach is the same, the sensitivity is worse. It is likely possible, though, to achieve a better sensitivity by setting up the experiment in a clever way. Thus, there is much room for improvement compared to the naive estimation used for the blue dashed lines of Fig. 6.
Finally, to give the reader an idea of the order of magnitude of t ⊥ , in Figure 7 we show the beam height at after a time t ⊥ (assuming t off > t ⊥ ) as a function of the radius of the receiver R c and v −1 −1 . The radius of the emitter is set to be R c /2. The range of axion masses is chosen as that for which the atmosphere is transparent. We also mark the height of the ionospheric layers for reference.

VIII. SUMMARY
In this manuscript, we have presented a detailed analysis of the signal of the echo method for axion dark matter detection proposed in Ref. [19]. We have found that the echo intensity grows in time until it saturates after a characteristic time that depends on the transverse axion velocity distribution. We have classified this velocity distribution effect into two cases: small and large dispersion, depending on whether the transverse velocity dispersion is smaller or larger than the transverse average velocity.
Using a Gaussian beam as a model and working in the paraxial approximation for the echo field, we have provided approximate analytical estimates of the saturation timescale and verified them by exact numerical integration for the caustic ring and the isothermal halo models. Moreover, we have shown that it is possible to achieve a separation in frequency between the beam and the echo bandwidths thanks to the average velocity of the axion flow along the beam's direction relative to the lab frame. This may help reduce the noise from atmospheric back-scattering and beam leakages.
Finally, we have derived an analytical expression for the power collected by a receiver dish concentric to the emitter. This expression allows us to optimize the experimental parameters in order to attain maximal signal to noise or minimal energy consumption. We have provided updated sensitivity estimates assuming a 100 m receiving dish. These experimental parameters are attainable with currently available technology. For instance, the receiver could be a radio telescope such as the Green Bank, Effelsberg or FAST, while the power source could be a high-power klystron, as used for radar transmitters or particle accelerators.
Our estimates agree with those of Ref. [19] when considering velocity distribution effects, and also with the fact that the outgoing beam divergence does not play any role in the signal. We leave for future work a detailed discussion of the role of the axion random phases and atmospheric noise, as well as a more practical description of possible experimental setups.
As we look for a solution that vanishes at t → −∞, ansatz (4.4) suggest the shift Ω → Ω−iη. We can see easily that (B3) has a pole at λ = −q − k 2 ⊥ /(2Ω) − iη. We perform the integrals in (B3) getting taking the lower part of Eq. (5.3), we have This result does not contain any effect from the axion velocity because the dish is big enough to collect all the power from the echo. In the second case, when R c /v ⊥ ≤ t ≤ t off , we can set the second erf function that appears in Eq. (5.3) equal to 1, since its argument is always large. For T c , we find where κ is defined as Notice that κ → 1 for R c /R → ∞.
With these results, the steady signal power can be written as where t e takes the value t off for t off < R c /v ⊥ and the value κ R c v −1 ⊥ for t off > R c /v ⊥ . Notice that if R c > R, κ is of order 1, so we can write the approximate result P c = 1 16