# The axion dark matter echo: a detailed analysis

Ariel Arza  Institute for Theoretical and Mathematical Physics, Lomonosov
Moscow State University (ITMP), 119991 Moscow, Russia  Elisa Todarello
Institut für Astroteilchenphysik, Karlsruhe Institute of Technology (KIT),
Hermann-von-Helmholtz-Platz 1, 76344 Eggenstein-Leopoldshafen, Germany

###### Abstract

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.

###### Contents

  1. I Introduction
  2. II Main features of the echo signal
    1. II.1 Energy and power of the echo
    2. II.2 Echo of a dish antenna beam
  3. III One-dimensional analysis
  4. IV Three-dimensional analysis
  5. V The role of the outgoing beam’s shape and the axion velocity distribution
    1. V.1 Small dispersion
    2. V.2 Large dispersion
    3. V.3 Consistency of the paraxial approximation
  6. VI Intensity map for the caustic ring and isothermal models
    1. VI.1 The caustic ring model
    2. VI.2 Isothermal halo model
  7. VII Echo power and sensitivity
  8. VIII Summary
  9. A Green’s function for the 1D computation
  10. B Green’s function for the paraxial equation
  11. C Approximated formula for $\delta p_{z}\,t\gg 1$
  12. D Details for the power collected by a receiving dish

##  I Introduction

The identity of dark matter is one of the biggest puzzles in modern science.
The QCD axion Peccei:1977hh ; Weinberg:1977ma ; Wilczek:1977pj , as well as
axion-like particles Svrcek:2006yi , are leading candidates for its
fundamental composition Preskill:1982cy ; Abbott:1982af ; Dine:1982ah ;
Arias:2012az . Under this assumption, several experimental efforts are being
carried out. See references Sikivie:2020zpn ; Irastorza:2018dyq  for recent
reviews. Most of these experiments are based on the axion to two photons
vertex Sikivie:1983ip  described by the interaction term
$gaF_{\mu\nu}\tilde{F}^{\mu\nu}$, where $a$ is the axion field, $F_{\mu\nu}$
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 Hertzberg:2018zte ;
Caputo:2018ljp ; Arza:2018dcy ; Caputo:2018vmy ; Wang:2020zur ; Arza:2020eik ;
Levkov:2020txo . 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. Arza:2019nta  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 Ghosh:2020hgd .

In Arza:2019nta , 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
Arza:2019nta , 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 Arza:2019nta
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
$\vec{x}\simeq 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 $\vec{x}\rightarrow 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.

###  II.1 Energy and power of the echo

The differential equation for the electromagnetic vector potential
$\vec{A}(x)$ with an axion background field $a(x)$ is given by

$\left(\partial_{t}^{2}-\nabla^{2}\right)\vec{A}=-g\partial_{t}a\vec{\nabla}\times\vec{A}\enspace,$
(1)

where the Coulomb gauge $\vec{\nabla}\cdot\vec{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 $\rho(\vec{p}\,)$. The axion field can
be written as

$a_{p}(t,\vec{x})={\sqrt{2\rho(\vec{p}\,)}\over
E_{p}}\sin(E_{p}t-\vec{p}\cdot\vec{x})\enspace,$ (2)

where $E_{p}=\sqrt{m^{2}+p^{2}}$. As $E_{p}=m+{\cal O}(p^{2})$, we will use
$E_{p}\approx 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
$\vec{A}^{(0)}$, the first order correction $\vec{A}^{(1)}$ is determined by

$\left(\partial_{t}^{2}-\nabla^{2}\right)\vec{A}^{(1)}=-g\partial_{t}a\vec{\nabla}\times\vec{A}^{(0)}\enspace.$
(3)

Let us expand the source field $\vec{A}^{(0)}$ as well as the correction
$\vec{A}^{(1)}$ in Fourier modes as

$\displaystyle\vec{A}^{(0)}(t,\vec{x})$ $\displaystyle=$
$\displaystyle\int{d^{3}k\over(2\pi)^{3}}\,\hat{e}\,{\cal
A}_{0}(\vec{k}\,)\,e^{-i(\omega(k)t-\vec{k}\cdot\vec{x})}$ (4)
$\displaystyle\vec{A}^{(1)}(t,\vec{x})$ $\displaystyle=$
$\displaystyle\int{d^{3}k\over(2\pi)^{3}}\,\hat{k}\times\hat{e}\,A^{(1)}(t,\vec{k}\,)\,e^{i\vec{k}\cdot\vec{x}}\enspace.$
(5)

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
$\hat{e}$. $A^{(1)}$ is polarized necessarily in the direction
$\hat{k}\times\hat{e}$, except for small corrections of the order of $p$. We
have also defined $\omega(k)\equiv|\vec{k}|$. Keeping only terms relevant to
the stimulated axion decay, i.e. the photon momentum modes $\vec{k}$ and
$\vec{q}=\vec{p}-\vec{k}$, Eq. (3) becomes

$\left(\partial_{t}^{2}+\omega(q)^{2}\right)A^{(1)}(t,\vec{q}\,)=-ig\,\sqrt{\rho(\vec{p}\,)\over
2}\,\omega(k)\,{\cal A}_{0}(\vec{k})^{*}e^{-i(m-\omega(k))t}\enspace.$ (6)

By looking at the source term of the above equation, we see that when
$m>\omega(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 $\omega(q)$ is equal or very close to $m-\omega(k)$, i.e. when
$\omega(k)\approx m/2$. Now we look for a resonant solution using the ansatz
$A^{(1)}(t,\vec{q}\,)={\cal A}(t,\vec{q}\,)e^{-i\omega(q)\,t}$ where ${\cal
A}(t,\vec{q}\,)$ varies slowly in time with respect to $A^{(1)}(t,\vec{q}\,)$.
Neglecting second derivatives in Eq. (6), we get

$\partial_{t}\,{\cal A}(t,\vec{q}\,)={g\over 2}\,\sqrt{\rho(\vec{p}\,)\over
2}\,{\cal
A}_{0}(\vec{k}\,)^{*}\,e^{i\,\epsilon(\vec{k},\vec{p}\,)\,t}\enspace,$ (7)

where

$\epsilon(\vec{k},\vec{p}\,)=\omega(k)+\omega(|\vec{p}-\vec{k}|)-m\enspace.$
(8)

Solving (7) with the initial condition ${\cal A}(0,\vec{q}\,)=0$, we obtain

${\cal A}(t,\vec{q}\,)={g\over 2}\,\sqrt{\rho(\vec{p}\,)\over 2}\,{\cal
A}_{0}(\vec{k}\,)^{*}\,e^{i\,\epsilon(\vec{k},\vec{p}\,)\,t}\left(\sin(\epsilon(\vec{k},\vec{p}\,)\,t/2)\over\epsilon(\vec{k},\vec{p}\,)/2\right)\enspace.$
(9)

In the limit $\epsilon t\gg 1$, we can substitute $\sin(\epsilon
t/2)/\epsilon\rightarrow\pi\delta(\epsilon)$ in Eq. (9). Integrating over all
axion momentum modes, we obtain the energy of the echo as

$U={\pi\over 4}g^{2}t\int d^{3}p\,{\partial^{3}\rho\over\partial p^{3}}\int
d^{3}k\,{\partial^{3}U_{0}\over\partial
k^{3}}\,\delta(\epsilon(\vec{k},\vec{p}\,))\enspace,$ (10)

where $U_{0}$ is the energy of the outgoing beam. In our non-relativistic
limit, $\epsilon\approx 2\omega(k)-m-p_{\parallel}$, where $p_{\parallel}$ is
the component of $\vec{p}$ in the direction of $\vec{k}$. Assuming that the
incident beam propagates mostly in a preferred direction, the direction of
$\vec{k}$ and $\vec{p}_{\parallel}$ are fixed, so

$U={\pi\over 4}\,g^{2}\,t\int dp_{\parallel}\,{\partial\rho\over\partial
p_{\parallel}}\int d\omega\,{\partial
U_{0}\over\partial\omega}\,\delta(\epsilon(\omega,p_{\parallel}))\enspace.$
(11)

Let’s perform the integral (11) in two scenarios: when the bandwidth of the
beam $\delta\omega$ is much bigger than the bandwidth of the axion $\delta p$
and in the opposite case. For $\delta\omega\gg\delta p$, the axion
distribution can be considered as if it were condensed at a single value
$\bar{p}_{\parallel}$. We then approximate $\partial\rho/\partial
p_{\parallel}=\rho\,\delta(p_{\parallel}-\bar{p}_{\parallel})$. We
straightforwardly find

$U={\pi\over 8}g^{2}\rho\,{\partial U_{0}\over\partial\omega_{*}}\,t\enspace,$
(12)

where

$\omega_{*}={m+\bar{p}_{\parallel}\over 2}\enspace.$ (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_{\rm off}$, the energy of the beam is $U_{0}=P_{0}\,t_{\rm off}$ and the
power of the echo is

$P={\pi\over 8}g^{2}\rho\,{\partial P_{0}\over\partial\omega_{*}}\,t_{\rm
off}\enspace.$ (14)

For $\delta\omega\ll\delta p$, we have $\partial
U_{0}/\partial\omega=U_{0}\,\delta(\omega-\bar{\omega})$. We find

$U={\pi\over
4}g^{2}\,U_{0}\,{\partial\rho\over\partial{p_{\parallel}}_{*}}\,t\enspace,$
(15)

as well as

$P={\pi\over
4}g^{2}\,P_{0}\,{\partial\rho\over\partial{p_{\parallel}}_{*}}\,t_{\rm
off}\enspace,$ (16)

where

${p_{\parallel}}_{*}=2\bar{\omega}-m\enspace.$ (17)

For the DFSZ model of the QCD axion, an input power of $1\,\text{kW}$ over a
bandwidth of $1\,\text{kHz}$ working for one hour111For instance, using a
klystron., provides a total echo power of $P\sim
10^{-18}\,\text{W}\left(m\over 10^{-4}~{}\text{eV}\right)^{2}$, which is in
principle not too difficult to detect with current radio astronomy technology.
Unfortunately, as discussed in Ref. Arza:2019nta , 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 $\Omega$, momentum
conservation implies $\Omega\sin(\chi)=mv_{\perp}$, where $\chi$ is the angle
formed by the trajectories of both photons and $v_{\perp}$ is the component of
the axion velocity perpendicular to the incident photon momentum. From energy
conservation, $\Omega$ is equal to $m/2$ plus small corrections, therefore
$\sin(\chi)\simeq 2v_{\perp}$. As the axion velocities are of the order of
$10^{-3}$ or smaller, we can write $\chi\simeq 2v_{\perp}$. 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}~{}\text{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}~{}\text{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.

###  II.2 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,\vec{x})$ and outgoing beam magnetic field
$\vec{B}^{(0)}(t,\vec{x})$, the echo vector potential field
$\vec{A}^{(1)}(t,\vec{x})$ is determined by

$\vec{A}^{(1)}(t,\vec{x})=-g\int d^{3}x^{\prime}\int
dt^{\prime}\,{\delta(t-t^{\prime}-|\vec{x}-\vec{x}\,^{\prime}|)\over
4\pi|\vec{x}-\vec{x}\,^{\prime}|}\,\partial_{t^{\prime}}a(t^{\prime},\vec{x}\,^{\prime})\vec{B}^{(0)}(t^{\prime},\vec{x}\,^{\prime})\enspace.$
(18)

We write the axion dark matter field, in a large volume $V$, as an expansion
in momentum modes as

$a(t,\vec{x})=\frac{a_{0}}{2}\,e^{-imt}\sqrt{V\over(2\pi)^{3}}\int
d^{3}p\,f_{a}(\vec{p}\,)\
e^{i(\vec{p}\cdot\vec{x}+\phi_{\vec{p}})}+c.c.\enspace.$ (19)

where $\vec{p}$ is the axion mometum in the experiment rest frame and
$\phi_{\vec{p}}$ random phases. These phases are not known, and are usually
modeled as uniformly distributed random numbers Foster:2017hbq ; Hui:2020hbq ;
Centers:2019dyn . In this article, we limit ourselves to considering the
ensemble average

$\left<e^{-i(\phi_{\vec{p}}-\phi_{\vec{p}\,^{\prime}})}\right>_{\rm
ens}={(2\pi)^{3}\over V}\delta^{3}(\vec{p}-\vec{p}\,^{\prime})\enspace,$ (20)

where $\left<\right>_{\rm 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}(\vec{p}\,)$ as

$\int d^{3}p\,|f_{a}(\vec{p}\,)|^{2}=1\enspace,$ (21)

such that the energy density in the volume $V$, averaged over time and random
phases is

$\rho=\frac{1}{2}m^{2}|a_{0}|^{2}\enspace.$ (22)

We consider a parabolic antenna of radius $R$ fed by a plane wave with
electric field amplitude $E_{0}(\omega)$, $\omega$ 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 $\omega R^{2}$ or larger from the
emission spot. Assuming that the center of the antenna is located at
$\vec{x}=0$ and that the central component of the beam propagates along
$\hat{z}$, the antenna’s electric field can be written in spherical
coordinates as

$\vec{E}^{(0)}(t,\vec{x},\omega)=\hat{x}\,i{E_{0}(\omega)\over 2}\omega
R^{2}\,{J_{1}(\omega R\sin(\theta))\over\omega
R\sin(\theta)}{e^{-i\omega(t-r)}\over r}\enspace,$ (23)

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 along $\hat{x}$. For
Eq. (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\gg\omega R^{2}$.

We perform integral (18) in the limit $\omega R\gg 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 $\omega R^{2}$, we can assume
$|\vec{x}\,^{\prime}|\gg R$. The field at the emission spot can then be found
using the limit $\vec{x}\rightarrow 0$. Second, the beamwidth of the outgoing
wave is of the order $(\omega R)^{-1}$, allowing us to approximate
$\theta\simeq 0$. If the outgoing beam is emitted continuously from $t=0$, the
echo vector potential at the emission spot reduces to

$\displaystyle\vec{A}^{(1)}(t,0,\Omega)$ $\displaystyle=$
$\displaystyle\hat{y}\,{gma_{0}^{*}E_{0}\over 4}\,e^{i\Omega t}{R\over
4\pi}\sqrt{V\over(2\pi)^{3}}\int
d^{3}p\,f_{a}(\vec{p}\,)^{*}e^{-i\phi_{\vec{p}}}\int_{0}^{\infty}dt^{\prime}\int_{0}^{t^{\prime}}dr^{\prime}\,e^{i\epsilon
r^{\prime}}$ (24)
$\displaystyle\delta(t-t^{\prime}-r^{\prime})\int_{0}^{\pi/2}d\theta^{\prime}\,J_{1}(\omega
R\,\theta^{\prime})\int_{0}^{2\pi}d\varphi^{\prime}\,e^{-i\vec{p}_{\perp}\cdot\hat{\rho}(\varphi^{\prime})r^{\prime}\theta^{\prime}}\enspace,$

where we have defined

$\displaystyle\Omega$ $\displaystyle=$ $\displaystyle m-\omega$ (25)
$\displaystyle\epsilon$ $\displaystyle=$ $\displaystyle 2\omega-m-p_{z}$ (26)
$\displaystyle\vec{p}_{\perp}$ $\displaystyle=$
$\displaystyle\hat{x}p_{x}+\hat{y}p_{y}\enspace,$ (27)

and
$\hat{\rho}(\varphi^{\prime})=\hat{x}\sin(\varphi^{\prime})+\hat{y}\cos(\varphi^{\prime})$.
After integration, we get

$\vec{A}^{(1)}(t,0,\Omega)=\hat{y}\,{gma_{0}^{*}E_{0}\over
8\omega}\,e^{i\Omega t}\sqrt{V\over(2\pi)^{3}}\int
d^{3}p\,f_{a}(\vec{p}\,)^{*}e^{-i\phi_{\vec{p}}}\int_{0}^{\xi/2}du\,e^{i\epsilon
u}\enspace,$ (28)

where $u=t-t^{\prime}$ and $\xi=\text{min}(t,mR/p_{\perp})$. To obtain the
result above we have extended the upper integration limit $d\theta^{\prime}$
to infinity. This is justified for $\omega R\gg 1$ and
$|\vec{p}_{\perp}|r^{\prime}\gg 1$, as then the largest contribution to the
integral comes for $\theta^{\prime}\ll\pi/2$.

To find the spectral intensity, we average the quantity
$\Omega^{2}A^{(1)}(t,0,\Omega)^{2}$ over fast time oscillations and over the
random phases $\phi_{\vec{p}}$. Assuming also that the axion momentum
distribution can be separated in a forward and transverse component
as222Although the axion phase space distribution cannot in general be
factorized, with this simplification, we can make analytical progress and gain
useful insights in the phenomenology of the signal.

$f_{a}(\vec{p}\,)=f_{a}(p_{z})f_{a}(\vec{p}_{\perp})\enspace,$ (29)

we have

${\partial I(t)\over\partial\Omega}={1\over 8}g^{2}\rho\,{dI_{0}(\omega)\over
d\omega}\int d^{2}p_{\perp}|f_{a}(\vec{p}_{\perp})|^{2}\int
dp_{z}|f_{a}(p_{z})|^{2}\int du\int
du^{\prime}e^{i\epsilon(u-u^{\prime})}\enspace.$ (30)

Consistently with Eq. (21), we normalize $f_{a}(p_{z})$ and $f_{a}(p_{\perp})$
as

$\int dp_{z}|f_{a}(p_{z})|^{2}=1$ (31)

and

$\int d^{2}p_{\perp}|f_{a}(\vec{p}_{\perp})|^{2}=1\enspace,$ (32)

respectively. From now on, we assume for $f_{a}(p_{z})$ the Gaussian shape

$|f_{a}(p_{z})|^{2}={1\over\sqrt{2\pi}\delta
p_{z}}e^{-{(p_{z}-\left<p_{z}\right>)^{2}\over 2\delta p_{z}^{2}}}\enspace.$
(33)

In our notation, $\left<\right>$ stands for averaging over the distribution
$|f_{a}(\vec{p}\,)|^{2}$, i.e. the expectation value for every function
$\hat{Q}(\vec{p}\,)$ is

$\left<\hat{Q}(\vec{p}\,)\right>=\int d^{3}p\
\hat{Q}(\vec{p}\,)\,|f_{a}(\vec{p}\,)|^{2}\enspace.$ (34)

Thus $\left<p_{z}\right>$ is the average value for $p_{z}$ and $\delta
p_{z}=\sqrt{\left<p_{z}^{2}\right>-\left<p_{z}\right>^{2}}$ the dispersion.
Notice that, with the normalization Eq. (21), if $\hat{Q}$ does not depend on
one of the momentum components, the corresponding integration trivially yields
1.

At resonance, i.e. when $\omega=(m+\left<p_{z}\right>)/2\equiv\omega_{*}$, we
have $\int dp_{z}|f_{a}(p_{z})|^{2}e^{i\epsilon(u-u^{\prime})}=e^{-{\delta
p_{z}^{2}\over 2}(u-u^{\prime})^{2}}$. Approximating

$e^{-{\delta p_{z}^{2}\over
2}(u-u^{\prime})^{2}}\rightarrow{\sqrt{2\pi}\over\delta
p_{z}}\delta(u-u^{\prime})\enspace,$ (35)

valid for $\delta p_{z}t\gg 1$, we get333It is true for almost all the axion
masses we are considering in this work

${\partial I(t)\over\partial\Omega_{*}}={1\over 8}\sqrt{\pi\over
2}{g^{2}\rho\over\delta
p_{z}}\,\text{min}\left(t,R\left<v_{\perp}^{-1}\right>\right)\,{dI_{0}\over
d\omega_{*}}\enspace,$ (36)

where $\vec{v}_{\perp}=\vec{p}_{\perp}/m$ is the axion transverse velocity and
$\Omega_{*}=(m-\left<p_{z}\right>)/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 ($\vec{v}=\vec{p}/m$), the isothermal
sphere Turner:1985si , and the caustic ring model Duffy:2008dk . In the
isothermal model, $|f_{a}(\vec{v}\,)|^{2}$ has a Maxwell-Boltzmann behavior
with dispersion $\delta v=270\,\text{km/s}$ and average velocity given by the
velocity of the Sun in the galactic rest frame, i.e.
$|\left<\vec{v}\,\right>|\simeq 230\,\text{km/s}$. In the caustic ring model,
the local dark matter is dominated by a single cold flow with velocity
dispersion $\delta v<70\,\text{m/s}$ and with an average velocity of
$|\left<\vec{v}\,\right>|\simeq 290\,\text{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

$|f_{a}(\vec{v}_{\perp})|^{2}={1\over\pi\delta
v_{\perp}^{2}}\,e^{-{(\vec{v}_{\perp}-\vec{v}_{p})^{2}\over\delta
v_{\perp}^{2}}}\enspace,$ (37)

where for the isothermal model we make $\delta v_{\perp}=\sqrt{2/3}\,\delta v$
in order to get the Maxwell-Boltzmann distribution for $|f_{a}(\vec{v})|^{2}$.
For the caustic ring model, as $\delta v_{\perp}$ is usually much smaller than
$v_{p}$, we take the limit $\delta v_{\perp}\rightarrow 0$, getting a delta
function centered at $\vec{v}_{\perp}=\vec{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
$\delta v_{\perp}\geq v_{p}$), however the direction of the big flow is
determined by the position of the IRAS Sikivie:2001fg , Planck Banik:2017ygz
and GAIA Chakrabarty:2020qgm  triangles with an uncertainty of $0.01$ radians.
It implies that $v_{p}$ cannot be reduced to values smaller than
$5\,\text{km/s}$, clearly much larger than $\delta v_{\perp}$.

To compute the quantity $\left<v_{\perp}^{-1}\right>$ that appears in Eq.
(36), we use Eq. (37). We explicitly get

$\left<v_{\perp}^{-1}\right>={\sqrt{\pi}\over\delta
v_{\perp}}e^{-{v_{p}^{2}\over 2\delta
v_{\perp}^{2}}}I_{0}\left({v_{p}^{2}\over 2\delta
v_{\perp}^{2}}\right)\enspace,$ (38)

where $I_{0}(x)$ is the modified Bessel function of the first kind of order 0.
For $v_{p}\gg\delta v_{\perp}$, compatible with the caustic ring model, we
have $\left<v_{\perp}^{-1}\right>\simeq v_{p}^{-1}$ while for $v_{p}\ll\delta
v_{\perp}$, which is compatible with the isothermal model when the outgoing
beam is pointed mainly in the direction of the axion wind, we have
$\left<v_{\perp}^{-1}\right>\simeq\sqrt{\pi}/\delta v_{\perp}$. It is also
useful to compute $\left<v_{\perp}\right>$, for the limit $v_{p}\gg\delta
v_{\perp}$ we have $\left<v_{\perp}\right>\simeq v_{p}$, while for
$v_{p}\ll\delta v_{\perp}$ we have $\left<v_{\perp}\right>={\sqrt{\pi}\over
2}\delta v_{\perp}$.

In this analysis, we have seen how the axion velocity distribution influences
the echo signal. Eq. (36) shows the effects explicitly. For
$t<R\left<v_{\perp}^{-1}\right>$ 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.1. After $t=R\left<v_{\perp}^{-1}\right>$, 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. (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 $\hat{z}$ 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 the $\hat{x}$-direction, which implies
that the echo is linearly polarized in the $\hat{y}$-direction. In this setup,
Eq. (3) can be written as

$(\partial_{t}^{2}-\partial_{z}^{2})A^{(1)}=-g\partial_{t}a\,B^{(0)}\enspace,$
(39)

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_{\rm off}$ and, for
simplicity, we assume that during this period it is emitted with a constant
amplitude 444We ignore transitions at $t=0$ and at $t=t_{\rm off}$.. This
means that the RHS of Eq. (39) is non-null in the regions

$\begin{array}[]{ccccccc}0&<&z&<&t,&\ \ \ \ \text{for}&t<t_{\rm off}\\\
t-t_{\rm off}&<&z&<&t,&\ \ \ \ \text{for}&t>t_{\rm off}\enspace.\end{array}$
(40)

We write the axion field in terms of its Fourier expansion as

$a(t,z)={a_{0}\over 2}\,e^{-imt}\sqrt{L\over 2\pi}\int
dp_{z}\,f_{a}(p_{z})\,e^{ip_{z}z+i\phi_{p_{z}}}\enspace,$ (41)

where $f_{a}(p_{z})$ is normalized as in Eq. (31). The axion energy density in
the length $L$, averaged over time and random phases555Analogously to the 3D
case, here we use the ensemble average
$\left<e^{-i(\phi_{p_{z}}-\phi_{p_{z}^{\prime}})}\right>_{\rm ens}={2\pi\over
L}\delta(p_{z}-p_{z}^{\prime})$ for the random phases $\phi_{p_{z}}$., is
related to $a_{0}$ as in Eq. (22).

The magnetic field of the incident beam is a plane wave given by

$B^{(0)}(t,z,\omega)={B_{0}(\omega)\over 2}\,e^{-i\omega(t-z)}\enspace,$ (42)

where the field amplitude $B_{0}$ is related to the beam intensity by

$I_{0}(\omega)={B_{0}(\omega)^{2}\over 2}\enspace.$ (43)

Plugging Eqs. (41) and (42) into Eq. (39) and keeping only terms relevant to
the stimulated decay, we have

$(\partial_{t}^{2}-\partial_{z}^{2})A^{(1)}=-i{g\over
4}ma_{0}^{*}B_{0}(\omega)\sqrt{L\over 2\pi}\int
dp_{z}\,f_{a}(p_{z})^{*}\,e^{i(\omega-p_{z})z}\,e^{i\Omega
t}\,e^{-i\phi_{p_{z}}}\enspace.$ (44)

For every $\omega$, the echo wave is sourced by a term that has frequency
$\Omega$, so we write the echo vector potential as

$A^{(1)}(t,z,\Omega)={1\over 2}{\cal
A}(t,z,\Omega)\,e^{i\Omega(t+z)}\enspace.$ (45)

As we look for resonant solutions, we assume $|\partial_{z}{\cal
A}|\ll\Omega{\cal A}$ and $|\partial_{t}{\cal A}|\ll\Omega{\cal A}$.
Neglecting second derivatives of ${\cal A}$, we get

$(\partial_{t}-\partial_{z}){\cal A}=-{g\over
4\Omega}ma_{0}^{*}B_{0}(\omega)\sqrt{L\over 2\pi}\int
dp_{z}\,f_{a}(p_{z})^{*}\,e^{i\epsilon z}\,e^{-i\phi_{p_{z}}}\enspace.$ (46)

The general solution of Eq. (46) is

${\cal A}(\Omega,t,z)=-{g\over 4\Omega}ma_{0}^{*}B_{0}(\omega)\sqrt{L\over
2\pi}\int dp_{z}\,f_{a}(p_{z})^{*}\,e^{-i\phi_{p_{z}}}\int dz^{\prime}\int
dt^{\prime}\,G(t-t^{\prime},z-z^{\prime})\,e^{i\epsilon z^{\prime}}\enspace,$
(47)

where $G(t,z)$ is the retarded Green’s function associated to the differential
operator $\partial_{t}-\partial_{z}$. The explicit calculation of $G(t,z)$ is
found in Appendix A. Using Eq. (97), we have

${\cal A}(\Omega,t,z)=-{g\over 4\Omega}ma_{0}^{*}B_{0}(\omega)\,\sqrt{L\over
2\pi}\int dp_{z}\,f_{a}(p_{z})^{*}\,e^{-i\phi_{p_{z}}}\int
dz^{\prime}\,e^{i\epsilon z^{\prime}}$ (48)

where the integration boundaries are

$\displaystyle\begin{cases}0<z^{\prime}<\frac{t+z}{2}&\quad\mathrm{for}\quad
t<t_{\rm off}-z\\\ \frac{t-t_{\rm
off}+z}{2}<z^{\prime}<\frac{t+z}{2}&\quad\mathrm{for}\quad t>t_{\rm
off}-z\enspace.\end{cases}$ (49)

The echo spectral intensity, at $z=0$, averaged over fast oscillations and
random phases is

${\partial I(t,\Omega)\over\partial\Omega}={1\over 8}g^{2}\rho\,{dI_{0}\over
d\omega}\int dp_{z}\,|f_{a}(p_{z})|^{2}\int dz^{\prime}\int
dz^{\prime\prime}\,e^{i\epsilon(z^{\prime}-z^{\prime\prime})}\enspace.$ (50)

Taking Eq. (33) for $|f_{a}(p_{z})|^{2}$ and the approximation Eq. (35), Eq.
(50), evaluated at resonance, becomes

${\partial I(t)\over\partial\Omega_{*}}={1\over 8}\sqrt{\pi\over
2}{g^{2}\rho\over\delta p_{z}}\,t_{<}\,{dI_{0}\over d\omega_{*}}\enspace,$
(51)

where $t_{<}$ is the minimum between $t$ and $t_{\rm 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

${\vec{B}}^{(0)}(t,\vec{x},\omega)=\hat{\epsilon}\,\frac{B_{0}(\omega)}{2}\,\eta(\vec{x},\omega)\
e^{-i\omega(t-z)}\enspace,$ (52)

where $\hat{\epsilon}$ is a constant polarization vector. The function
$\eta(\vec{x},\omega)$ defines the spatial shape of the beam. It is normalized
as $\eta(\vec{0},\omega)=1$, such that $B_{0}(\omega)$ is related to the beam
time-averaged intensity $I_{0}(\omega)$, at $\vec{x}=0$, by

$I_{0}(\omega)\equiv I_{0}(\vec{0},\omega)={B_{0}(\omega)^{2}\over
2}\enspace.$ (53)

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_{\rm off}$, it extends from
$z=t-t_{\rm off}$ to $z=t$. In other words, the integration domain in the $z$
direction is given by Eq. (49), while for the transverse directions, it is
determined by $\eta(\vec{x},\omega)$.

We describe the local dark matter axion field as a superposition of plane
waves, in the same way as in Section II.2 (see Eqs. (19)-(22)). In addition,
we also assume that $f_{a}(\vec{p}\,)$ can be factorized in its forward and
transverse parts (see Eqs. (29), (31) and (32)). Discarding non-resonant
terms, Eq. (3) reads

$\left(\partial_{t}^{2}-\nabla^{2}\right)\vec{A}^{(1)}=-i\hat{\epsilon}{g\over
4}ma_{0}^{*}B_{0}(\omega)e^{i\Omega(t+z)}\,\eta(\vec{x},\omega)\sqrt{V\over(2\pi)^{3}}\int
d^{3}p\ f_{a}(\vec{p}\,)^{*}\ e^{i\epsilon
z-i\vec{p}_{\perp}\cdot\vec{x}_{\perp}-i\phi_{\vec{p}}}\enspace,$ (54)

where $\vec{x}_{\perp}=(x,y,0)$, while $\Omega$ and $\epsilon$ are defined in
Eqs. (25) and (26), respectively. We now write the echo field as

$\vec{A}^{(1)}(t,\vec{x})=\frac{\hat{\epsilon}}{2}\,{\cal
A}(t,\vec{x},\Omega)\ e^{i\Omega(t+z)}\enspace.$ (55)

Notice that $\Omega$ can be written as $\Omega=\omega-p_{z}-\epsilon$. As at
resonance $\epsilon\rightarrow 0$, the frequencies of the echo are shifted
with respect to $\omega$ 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 $\omega$, 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. (55) into Eq. (54) and use the fact that the outgoing beam as
well as the echo field satisfy the paraxial approximation

$|\partial_{i}^{2}{\cal A}(t,\vec{x},\Omega)|\ll|\Omega\,\partial_{i}{\cal
A}(t,\vec{x},\Omega)|\enspace,$ (56)

with $i=t,z$. For the echo, we will check in Subsection V.3 that our results
are consistent with the paraxial assumption. In this limit, Eq. (54) reduces
to the paraxial equation with a source term

$\left(2i\Omega(\partial_{t}-\partial_{z})-\nabla_{\perp}^{2}\right){\cal
A}=-i{g\over
2}ma_{0}^{*}B_{0}(\omega)\,\eta(\vec{x},\omega)\sqrt{V\over(2\pi)^{3}}\int
d^{3}p\ f_{a}(\vec{p}\,)^{*}\ e^{i\epsilon
z-i\vec{p}_{\perp}\cdot\vec{x}_{\perp}-i\phi_{\vec{p}}}\enspace.$ (57)

The solution is

${\cal A}(t,\vec{x},\Omega)=-i{g\over
2}ma_{0}^{*}B_{0}(\omega)J(t,\vec{x},\Omega)\enspace,$ (58)

where

$J(t,\vec{x},\Omega)=\sqrt{V\over(2\pi)^{3}}\int d^{3}p\ f_{a}(\vec{p}\,)^{*}\
e^{-i\phi_{\vec{p}}}\int d^{4}x^{\prime}\
G(t-t^{\prime},\vec{x}-\vec{x}\,^{\prime})\,\eta(\vec{x}\,^{\prime},\omega)\
e^{i\epsilon
z^{\prime}-i\vec{p}_{\perp}\cdot\vec{x}\,^{\prime}_{\perp}}\enspace.$ (59)

As derived in Appendix B, the Green’s function for the paraxial equation is

$G(t,\vec{x}\,)=-\frac{\Theta(-z)}{4\pi z}\delta(t+z)\
e^{\frac{i\Omega}{2z}\vec{x}_{\perp}^{2}}\enspace.$ (60)

The echo spectral intensity averaged over fast time oscillations is given by

${\partial I(t,\vec{x})\over\partial\Omega}={1\over 2}g^{2}\rho{dI_{0}\over
d\omega}\,\Omega^{2}|J(t,\vec{x},\Omega)|^{2}\enspace.$ (61)

Although Eq. (61) is the most general formula of this work, its value at $z=0$
can be reduced notably in the limit $\delta p_{z}\,t\gg 1$. Since our goal is
the estimation of the echo amplitude, for which corrections ${\cal O}(p)$ are
irrelevant, we use the approximation $\Omega=\omega-
p_{z}-\epsilon\approx\omega$ for the following analytical computations. With
the axion momentum distribution Eq. (33), and after averaging over random
phases, the spectral echo intensity, at resonance, can be written as

${\partial I(t,\vec{x}_{\perp})\over\partial\Omega_{*}}={1\over
8}\sqrt{\pi\over 2}{g^{2}\rho\over\delta p_{z}}\,{\cal
T}(t,\vec{x}_{\perp})\,{dI_{0}\over d\omega_{*}}\enspace,$ (62)

where

${\cal T}(t,\vec{x}_{\perp})=8\omega_{*}^{2}\,\left<\int
dz^{\prime}\,|\xi(z^{\prime},\vec{x}_{\perp},\vec{p}_{\perp},\omega_{*})|^{2}\right>\enspace,$
(63)

and

$\xi(z^{\prime},\vec{x}_{\perp},\vec{p}_{\perp},\omega_{*})={1\over 4\pi
z^{\prime}}\int
d^{2}x_{\perp}^{\prime}\,\eta(\omega_{*},\vec{x}\,^{\prime})\,e^{-i\vec{p}_{\perp}\cdot\vec{x}\,^{\prime}_{\perp}}\,e^{-{i\omega_{*}\over
2z^{\prime}}(\vec{x}_{\perp}-\vec{x}\,^{\prime}_{\perp})^{2}}\enspace.$ (64)

The integration domain for $z^{\prime}$ is specified in Eq. (49), evaluated at
$z=0$. For a derivation of Eqs. (62) and (63) see Appendix C. Notice that Eq.
(62) is the three dimensional analog of Eq. (51). It means that all the 3D
effects enter through in ${\cal T}(t,\vec{x}_{\perp})$.

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

${dI_{0}\over
d\omega}={I_{0}\over\sqrt{2\pi}\delta\omega}e^{-{(\omega-\bar{\omega})^{2}\over
2\delta\omega^{2}}}$ (65)

where $I_{0}$ is the total intensity, $\bar{\omega}$ the central value of the
distribution, and $\delta\omega$ the dispersion.

Thus, the power collected over a surface $S$ (see also Appendix C) can be
computed as

$P_{c}={1\over 16}\sqrt{\pi\over 2}{g^{2}\rho\over\Delta}\,P_{0}\,{1\over
S_{0}}\int_{S}d^{2}x_{\perp}\,{\cal T}(t,\vec{x}_{\perp})\enspace,$ (66)

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

$\Delta=\sqrt{\delta\omega^{2}+{\delta p_{z}^{2}\over 4}}\enspace.$ (67)

##  V The role of the outgoing beam’s shape and the axion velocity
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,\phi,z)$ it is given by

$\eta(\omega,\vec{x}\,)=\frac{R}{w(z)}\,e^{-\frac{r^{2}}{2w(z)^{2}}}\,e^{i\omega\frac{r^{2}}{2{\cal
R}(z)}}\enspace.$ (68)

Here $w(z)$ is the radius at which the beam intensity is reduced by $1/e$
compared to its axial value and ${\cal R}(z)$ is the radius of curvature of
the beam wavefronts. Their mathematical expressions are

$\displaystyle w(z)$ $\displaystyle=R\sqrt{1+{z^{2}\over z_{R}^{2}}}\enspace,$
$\displaystyle{\cal R}(z)$ $\displaystyle={z^{2}+z_{R}^{2}\over z}\enspace,$

where $R=w(0)$ and $z_{R}=\omega 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 $\omega R\gg 1$ (always
satisfied for the parameters’ values relevant to this work), we recover the
behavior of the dish antenna beam Eq. (23) for $z\geq z_{R}$. It is worth
mentioning that, in Eq. (68), we have ignored the Guoy phase
$\psi(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}\pi R^{2}$, then the
surface parameter entering Eq. (66) must be given by $S_{0}=\pi 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 ${\cal
T}(t,\vec{x}_{\perp})$ defined in Eq. (63). According to Eq. (62), ${\cal 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 ${\cal T}(t,\vec{x}_{\perp})$ is equivalent to the behavior of
the echo intensity.

We start our analysis, by computing $\xi$ from Eq. (64) using the beam model
Eq. (68). We get

$|\xi(z^{\prime},\vec{x}_{\perp},\vec{p}_{\perp},\omega_{*})|^{2}={1\over
4\omega_{*}^{2}}\,e^{-{\left(\vec{x}_{\perp}-2\vec{v}_{\perp}z^{\prime}\right)^{2}\over
R^{2}}}\enspace.$ (69)

Using Eq. (69), it seems that Eq. (63) could be very difficult to integrate
analytically in both $dz^{\prime}$ and $d^{2}p_{\perp}$, 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 ${\cal
T}(t,\vec{x}_{\perp},\vec{v}_{\perp})$ such that ${\cal
T}(t,\vec{x}_{\perp})=\left<{\cal
T}(t,\vec{x}_{\perp},\vec{v}_{\perp})\right>$. From Eq. (63), it is found to
be

$\displaystyle{\cal T}(t,\vec{x}_{\perp},\vec{v}_{\perp})$ $\displaystyle=$
$\displaystyle{\sqrt{\pi}\over 2}{R\over v_{\perp}}e^{-{r^{2}\over
R^{2}}\sin(\varphi)^{2}}\times$ (70)
$\displaystyle\begin{cases}\left(\operatorname{erf}\left(r\cos(\varphi)\over
R\right)+\operatorname{erf}\left(v_{\perp}t-r\cos(\varphi)\over
R\right)\right)&\quad\mathrm{for}\quad t<t_{\rm off}\\\
\left(\operatorname{erf}\left(v_{\perp}t-r\cos(\varphi)\over
R\right)-\operatorname{erf}\left(v_{\perp}(t-t_{\rm off})-r\cos(\varphi)\over
R\right)\right)&\quad\mathrm{for}\quad t>t_{\rm off}\enspace,\end{cases}$

where $\varphi$ is the polar angle formed by $\vec{x}_{\perp}$ and
$\vec{v}_{\perp}$. Now what is left is to integrate Eq. (70) over the velocity
distribution, i.e., we have to find

${\cal T}(t,\vec{x}_{\perp})=\int
d^{2}p_{\perp}|f_{a}(\vec{p}_{\perp})|^{2}\,{\cal
T}(t,\vec{x}_{\perp},\vec{v}_{\perp})\enspace$ (71)

where $|f_{a}(\vec{v}_{\perp})|^{2}$ is given in Eq. (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.2, we find null effects coming
from the outgoing beam divergence. Our analysis will be split into two cases;
small dispersion ($\delta v_{\perp}\ll v_{p}$) and large dispersion ($\delta
v_{\perp}\gg v_{p}$).

Before discussing the 3D effects in detail, as a consistency check, notice
that when the transverse velocities are not important, i.e.
$\vec{v}_{\perp}\rightarrow 0$, we get from Eq. (70) that ${\cal T}(t,0)=t$
for $t<t_{\rm off}$ and ${\cal T}(t,0)=t_{\rm off}$ for $t>t_{\rm off}$. In
other words, we recover the result from the one dimensional analysis of
Section III.

###  V.1 Small dispersion

When $\delta v_{\perp}\ll v_{p}$, Eq. (37) can be approximated as
$|f_{a}(\vec{v}_{\perp})|^{2}=\delta^{2}(\vec{v}_{\perp}-\vec{v}_{p})$, then
the integration in Eq. (71) becomes trivial. We have

$\displaystyle{\cal T}(t,\vec{x}_{\perp})$ $\displaystyle=$
$\displaystyle{\sqrt{\pi}\over 2}{R\over v_{p}}e^{-{r^{2}\over
R^{2}}\sin(\phi)^{2}}\times$ (72)
$\displaystyle\begin{cases}\left(\operatorname{erf}\left(r\cos(\phi)\over
R\right)+\operatorname{erf}\left(v_{p}\,t-r\cos(\phi)\over
R\right)\right)&\quad\mathrm{for}\quad t<t_{\rm off}\\\
\left(\operatorname{erf}\left(v_{p}\,t-r\cos(\phi)\over
R\right)-\operatorname{erf}\left(v_{p}\,(t-t_{\rm off})-r\cos(\phi)\over
R\right)\right)&\quad\mathrm{for}\quad t>t_{\rm off}\enspace,\end{cases}$

where $\phi$ is the angle formed by $\vec{x}_{\perp}$ and $\vec{v}_{p}$. We
will analyze Eq. (72) in two scenarios; a) $t_{\rm off}<R/v_{p}$ and b)
$t_{\rm off}>R/v_{p}$, each characterized by two regimes; $t<t_{\rm off}$ and
$t>t_{\rm off}$.

In scenario a), we can Taylor expand ${\cal T}(t,\vec{x}_{\perp})$ in terms of
the small quantities $v_{p}\,t/R$ for $t<t_{\rm off}$, and $v_{p}\,t_{\rm
off}/R$ for $t>t_{\rm off}$. We find

${\cal T}(t,\vec{x}_{\perp})=\begin{cases}t\,e^{-{r^{2}\over
R^{2}}}&\quad\mathrm{for}\quad t<t_{\rm off}\\\ t_{\rm off}\,e^{-{r^{2}\over
R^{2}}\sin(\phi)^{2}}e^{-{\left(v_{p}\,t-r\cos(\phi)\right)^{2}\over
R^{2}}}&\quad\mathrm{for}\quad t>t_{\rm off}\enspace.\end{cases}$ (73)

We see that for $t<t_{\rm off}$, ${\cal T}(t,\vec{x}_{\perp})$, and therefore
the intensity, grows linearly in time, featuring a transverse Gaussian shape
with effective radius $R$. For $t>t_{\rm off}$, i.e., after the outgoing beam
is turned off, the Gaussian keeps its maximum value reached at $t=t_{\rm off}$
and moves rigidly with velocity $\vec{v}_{p}\,$.

In scenario b), ${\cal T}(t,\vec{x}_{\perp})$ behaves exactly as the upper
formula in Eq. (73) for $t<R/v_{p}$. Later, for $R/v_{p}<t<t_{\rm off}$, we
observe from Eq. (72) that the echo spreads along the direction of
$\vec{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 $\vec{x}_{\perp}$ when
$t>{r\cos(\phi)\over v_{p}}$. The maximum saturated value is found for
locations where $\phi=0$ and $R<r\cos(\phi)<v_{p}\,t$. The corresponding value
for ${\cal T}$ within this region is twice the value obtained at $r=0$. We
name it

$t_{a}=\sqrt{\pi}{R\over v_{p}}\enspace.$ (74)

Finally, for $t>t_{\rm off}$, the “sausage” shape moves rigidly at the speed
$\vec{v}_{p}$ having a length $v_{p}\,t_{\rm off}$ and a width $2R$.

###  V.2 Large dispersion

Let’s discuss now the case $\delta v_{\perp}\gg v_{p}$. As in the previous
discussion, we split the analysis in two scenarios; a) $t_{\rm off}<R/\delta
v_{\perp}$ and b) $t_{\rm off}>R/\delta v_{\perp}$, and each of them in two
time regimes; $t<t_{\rm off}$ and $t>t_{\rm off}$.

In scenario a), we Taylor expand ${\cal
T}(t,\vec{x}_{\perp},\vec{v}_{\perp})$, given in Eq. (70), in terms of the
small values $v_{\perp}t/R$ for $t<t_{\rm off}$ and $v_{\perp}t_{\rm off}/R$
for $t>t_{\rm off}$. The result is, of course, identical to Eq. (73) after
performing the substitutions $v_{p}\rightarrow v_{\perp}$ and
$\phi\rightarrow\varphi$. Now, integrating over the velocity distribution, we
straightforwardly obtain the upper formula in Eq. (73) for $t<t_{\rm off}$.
For $t>t_{\rm off}$, the integration yields

$\displaystyle{\cal T}(t,\vec{x}_{\perp})$ $\displaystyle=$ $\displaystyle
t_{\rm off}\,e^{-{r^{2}\over R^{2}}}{1\over\pi\delta
v_{\perp}^{2}}\int_{0}^{\infty}dv_{\perp}v_{\perp}\,e^{-\left({t^{2}\over
R^{2}}+{1\over\delta
v_{\perp}^{2}}\right)v_{\perp}^{2}}\int_{0}^{2\pi}d\varphi\,e^{{2v_{\perp}t\,r\over
R^{2}}\cos(\varphi)}$ (75) $\displaystyle=$ $\displaystyle{R^{2}\,t_{\rm
off}\over R^{2}+\delta v_{\perp}^{2}t^{2}}\,e^{-{r^{2}\over R^{2}+\delta
v_{\perp}^{2}t^{2}}}\enspace.$

This behavior can be separated into two sub-regimes; $t_{\rm off}<t<R/\delta
v_{\perp}$ and $t>R/\delta v_{\perp}$. In the first one, the intensity maximum
stays constant at the value obtained at $t=t_{\rm off}$, and also, its radial
extension $R$ does not change. In the second one, the radial extension becomes
$\delta v_{\perp}t$, i.e., the energy spreads radially at velocity $\delta
v_{\perp}$. On the other hand, the value of the intensity decreases as
$(\delta v_{\perp}t)^{-2}$.

Scenario a) can be summarized as follows

${\cal T}(t,x_{\perp})=\begin{cases}t\,e^{-{r^{2}\over
R^{2}}}&\quad\mathrm{for}\quad t<t_{\rm off}\\\ t_{\rm off}\,e^{-{r^{2}\over
R^{2}}}&\quad\mathrm{for}\quad t_{\rm off}<t<R/\delta v_{\perp}\\\ {t_{\rm
off}R^{2}\over\delta v_{\perp}^{2}t^{2}}\,e^{-{r^{2}\over\delta
v_{\perp}^{2}t^{2}}}&\quad\mathrm{for}\quad t>R/\delta
v_{\perp}\enspace.\end{cases}$ (76)

In scenario b), ${\cal T}(t,\vec{x}_{\perp})$ behaves exactly as the upper
formula of Eq. (73) for $t<R/\delta v_{\perp}$. After that, for $R/\delta
v_{\perp}<t<t_{\rm off}$, the energy spreads radially leading the intensity to
saturate for all the positions satisfying $r<\delta v_{\perp}t$. The
corresponding saturated value for ${\cal T}(t,\vec{x}_{\perp})$ can be found
taking the limit $\operatorname{erf}(x)\rightarrow 1$ for the second erf
function appearing in the upper formula of Eq. (70). Thus, Eq. (71) gives us

$\displaystyle{\cal T}(t,\vec{x}_{\perp})$ $\displaystyle=$
$\displaystyle{\sqrt{\pi}\over 2}R{1\over\pi\delta
v_{\perp}^{2}}\int_{0}^{\infty}\,dv_{\perp}\,e^{-{v_{\perp}^{2}\over\delta
v_{\perp}^{2}}}\int_{0}^{2\pi}d\varphi\,e^{-{r^{2}\over
R^{2}}\sin(\varphi)^{2}}\left(\operatorname{erf}\left(r\cos(\varphi)\over
R\right)+1\right)$ (77) $\displaystyle=$ $\displaystyle{\pi\over
2}{R\over\delta v_{\perp}}e^{-{r^{2}\over 2R^{2}}}I_{0}\left(r^{2}\over
2R^{2}\right)\enspace.$

The maximum value of $\cal T$ in this saturated regime is found at $r=0$ and
is

$t_{b}={\pi\over 2}{R\over\delta v_{\perp}}\enspace.$ (78)

We can also find analytically how the system transits to this saturated value.
Evaluating Eq. (70) at $r=0$, Eq. (71) gives

$\displaystyle{\cal T}(t,0)$ $\displaystyle=$
$\displaystyle\sqrt{\pi}{R\over\delta
v_{\perp}^{2}}\int_{0}^{\infty}dv_{\perp}\,e^{-{v_{\perp}^{2}\over\delta
v_{\perp}^{2}}}\operatorname{erf}\left(v_{\perp}t\over R\right)$ (79)
$\displaystyle=$ $\displaystyle{R\over\delta v_{\perp}}\arctan\left(\delta
v_{\perp}t\over R\right)\enspace.$

With the same procedure, we find for $t>t_{\rm off}$

${\cal T}(t,0)={1\over 2\pi}{R\over\delta v_{\perp}}\left(\arctan\left(\delta
v_{\perp}t\over R\right)-\arctan\left(\delta v_{\perp}(t-t_{\rm off})\over
R\right)\right)\enspace,$ (80)

which in the limit $t\gg t_{\rm off}$ becomes

${\cal T}(t,0)={1\over 2\pi}{R^{2}\,t_{\rm off}\over\delta
v_{\perp}^{2}t^{2}}\enspace.$ (81)

As a summary of scenario b), we can write

${\cal T}(t,x_{\perp})=\begin{cases}t\,e^{-{r^{2}\over
R^{2}}}&\quad\mathrm{for}\quad t<R/\delta v_{\perp}\\\ {\pi\over
2}{R\over\delta v_{\perp}}e^{-{r^{2}\over 2R^{2}}}I_{0}\left(r^{2}\over
2R^{2}\right)&\quad\mathrm{for}\quad R/\delta v_{\perp}<t<t_{\rm off}\\\
{1\over 2\pi}{R^{2}\,t_{\rm off}\over\delta
v_{\perp}^{2}t^{2}}&\quad\mathrm{for}\quad t>t_{\rm off}\quad\mathrm{and}\quad
r=0\enspace.\end{cases}$ (82)

As a general remark, we point out that $t_{a}$ and $t_{b}$ are not only the
saturated values of $\cal T$ in the respective scenarios but also the
timescales over which the saturation is achieved. We can then define the
saturation time

$t_{\text{sat}}=\left({1\over t_{a}^{2}}+{1\over
t_{b}^{2}}\right)^{-1/2}\enspace.$ (83)

To conclude this discussion, we notice that, as long as $\bar{\omega}R^{2}\gg
R\left<v_{\perp}^{-1}\right>$, saturation occurs before the outgoing beam
reaches the far-field zone. This is explained by the fact that the intensity
saturates for $t\gtrsim R/v_{p}$ in the case of large average transverse
velocity and for $t\gtrsim R/\delta v_{\perp}$ in the case of large transverse
velocity dispersion.

###  V.3 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\,\text{rad}$ hecht . For the echo wave, the divergence angle $\chi$ is
related to the velocity $v_{s}$, at which the echo spreads in the transverse
direction, by

$\chi=\arctan(2v_{s})\enspace.$ (84)

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_{\text{sat}}\sim
R$. It implies that $v_{s}$ is of the order of $v_{p}$ or $\delta v_{\perp}$,
depending on which parameter dominates the saturation. For any cold dark
matter model, $v_{p}$ and $\delta v_{\perp}$ 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 models

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 Turner:1985si  and the caustic ring model Duffy:2008dk .

We write the echo intensity as

$\frac{dI(t,\vec{x}_{\perp})}{d\Omega}=\frac{1}{8}\sqrt{\frac{\pi}{2}}\frac{g^{2}\rho}{\delta
p_{z}}\,t_{\text{sat}}F(t,\vec{x}_{\perp})\frac{dI_{0}}{d\omega}\enspace,$
(85)

where $t_{\text{sat}}$ is defined in Eq. (83) and $F(t,\vec{x}_{\perp})$ is a
function whose maximum is of order 1 given by

$F(t,\vec{x}_{\perp})=\frac{1}{t_{\text{sat}}}\,\frac{\delta
p_{z}}{\sqrt{2\pi}}\,8\Omega^{2}\,|J(t,\vec{x}_{\perp},\Omega)|^{2}\enspace.$
(86)

Here $J(t,\vec{x}\,,\Omega)$ is given in Eq. (59) with the beam profile Eq.
(68). Notice that, as follows from the discussion of Section V, in the limit
$\delta p_{z}\,t\gg 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_{\text{sat}}$.

Figure 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_{\rm off}=7\times 10^{7}$. After $t_{\rm
off}$, the hot spot travels rigidly.

###  VI.1 The caustic ring model

The caustic ring model is a proposal for a complete description of the dark
matter phase-space distribution in the halo of our galaxy Duffy:2008dk . 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 six-dimensional 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 $(\vec{\nabla}\times\vec{v}\neq 0)$, the
caustics have the aforementioned tube shape Natarajan:2005ut .

Based on the location of two triangular features present in the IRAS
Sikivie:2001fg , Planck Banik:2017ygz  and Gaia Chakrabarty:2020qgm  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. Chakrabarty:2020qgm , we see that the big flow has a velocity of about
520 km/s in the galactic rest frame, implying a relative 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\,\text{km/s}$. The velocity dispersion of the flow was determined in Ref.
Banik:2015vts  to be at most of order 70 m/s. The local dark matter density in
this model is 1 GeV/cm3 or higher.

In this Section, we assume the axion phase space distribution Eqs. (33) and
(37), with a velocity dispersion $\delta 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.666As will be explained in
Section VII, the measurement time $t_{m}$ can be at most of order 1 s. Earth’s
revolution velocity is $v_{rev}\sim 30$ km/s, while its rate of change is
$\dot{v}_{rev}\sim\omega_{rev}v_{rev}$, where $\omega_{rev}=2\pi/\mathrm{yr}$.
During a time $t_{m}$, $v_{rev}$ changes by $\Delta v_{rev}\sim 2\times
10^{-7}\,v_{rev}$. Assuming that all $\Delta v_{rev}$ contributes to $v_{p}$,
during a time $t_{m}$ it makes the echo hot spot move by about 6 mm, which is
negligible compared to the receiver’s size. The same argument applies to
Earth’s rotation velocity. In conclusion, we can consider the axion’s velocity
relative to the beam as a constant during any given data collection event.
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. (86) at
resonance, i.e. for $\bar{\omega}=(m+\langle p_{z}\rangle)/2$. With the chosen
radius $R=4000/m$, we have $\delta p_{z}t_{\text{sat}}\gg 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_{\rm off}=1.4\times 10^{9}$. Thereafter, the hot spot stops
spreading and starts traveling rigidly. Fig 1 shows $F$ along the line $y=0$
at different times.

###  VI.2 Isothermal halo model

For the isothermal halo model, we also use the axion momentum distribution
Eqs. (33) and (37). In this case, the velocity dispersion is $\delta
v_{z}=\delta v_{\perp}/\sqrt{2}=270/\sqrt{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/cm3. 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,\vec{0}\,)$ as a function of the beam frequency $\omega$
at four different times. $F(t,\vec{0}\,)$ has been multiplied by the spectral
intensity inherited from the beam $d\tilde{I}_{0}/d\omega|_{\omega=\Omega}$,
whose maximum value has been normalized to 1. The beam spectral intensity Eq.
(65) normalized to its maximum value is shown in black. The echo peaks at
$\bar{\Omega}=(m-\langle p_{z}\rangle)/2$, while the beam is centered at
$\bar{\omega}=(m+\langle p_{z}\rangle)/2$. The echo bandwidth is given by
$\delta\omega\delta p_{z}/(2\Delta)$. It is then possible for the echo
frequency band to be completely separated from the beam frequency band if

$|\bar{\Omega}-\bar{\omega}|>\delta\omega+\frac{\delta\omega\delta
p_{z}}{2\Delta}\enspace.$ (87)

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.

Figure 2: Echo normalized intensity $F(t,x,y=0)$ in the case of saturation due
to velocity dispersion effects. All the parameters are the same as those of
Fig. 5. The spreading of the hot spot is clearly visible. Figure 3: Frequency
profile of $F(t,\vec{0}\,)$ times the spectral intensity
$d\tilde{I}_{0}/d\omega|_{\omega=\Omega}$ (colored lines) and rescaled beam
frequency profile $d\tilde{I}_{0}/d\omega$ (black line) at resonance
$\omega=\omega_{*}$ in the isothermal halo model. All the parameters are the
same as those of Fig. 5, with the addition of $\delta\omega=2\times 10^{-5}m$.

Figure 4: Echo normalized intensity $F(t,\vec{x}_{\perp})$ at resonance
$\omega=\omega_{*}$ in the case of saturation due to the transverse average
axion velocity. The beam is turned on at $t=0$ and turned off at $mt_{\rm
off}=1.4\times 10^{9}$. The radius of the emitter is set to $mR=4000$. The
axion velocity components are those of the caustic ring model, with minimum
transverse velocity $v_{x}=5$ km/s and $v_{z}=290$ km/s. The velocity
dispersion is $70$ km/s in the $z$ direction. The saturation time is
$mt_{\text{sat}}\approx mt_{a}=4.3\times 10^{8}$, so that $\delta
p_{z}t_{\text{sat}}=99$ and the approximation discussed in Section V holds.
Figure 5: Echo normalized intensity $F(t,\vec{x}_{\perp})$ at resonance
$\omega=\omega_{*}$ in the case of saturation due to velocity dispersion
effects. The beam is turned on at $t=0$ and turned off at $mt_{\rm
off}=7\times 10^{5}$. The radius of the emitter is set to $mR=100$. The axion
velocity components are those of the isothermal model, with $v_{x}=10$ km/s
and $v_{z}=230$ km/s, while the velocity dispersion is $270/\sqrt{3}$ km/s in
each direction. The saturation time is $mt_{\text{sat}}\approx
mt_{b}=2.1\times 10^{5}$, so that $\delta p_{z}t_{\text{sat}}=111$ and the
approximation discussed in Section V holds.

##  VII Echo power and sensitivity

Figure 6: Expected sensitivity of the echo method in the space of parameters
$(m,g)$ for the caustic ring (left) and isothermal model (right) of the
galactic halo. These plots assume that the method consumes an energy of
$10\,\text{MW\,year}$ per factor of two in axion mass range. The solid blue
lines correspond to the scenario where this energy is spent in an efficient
way while the dashed blue lines show the sensitivity for a fixed power source
of $10\,\text{MW}$ working for one year (see the text for more details). The
green regions are current bounds from axion dark matter haloscopes
DePanfilis:1987dk ; Hagmann:1990tj ; ADMX:2018gho ; ADMX:2019uok ;
ADMX:2018ogs ; Lee:2020cfj ; Jeong:2020cwz ; CAPP:2020utb ; HAYSTAC:2018rwy ;
HAYSTAC:2020kwv ; Alesini:2019ajt ; Alesini:2020vny ; McAllister:2017lkb ;
CAST:2021add  and the grey region corresponds to bounds from the CAST
experiment CAST:2017uph .

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

$s/n={P_{c}\over T_{n}}\sqrt{t_{m}\over B}\enspace,$ (88)

where $P_{c}$ is the collected power, $T_{n}$ the noise temperature, $B$ the
bandwidth of the signal given roughly by
$B\approx\text{min}(\delta\omega,\delta p_{z}/2)/(2\pi)$ 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

$P_{c}={1\over 16}\sqrt{\pi\over
2}{g^{2}\rho\over\Delta}\,t_{e}\,P_{0}\enspace.$ (89)

Here $t_{e}$ is basically the minimum between $t_{\rm off}$ and
$\kappa\,t_{\perp}$, where $\kappa$ is a parameter that depends on the ratio
$R_{c}/R$ and $t_{\perp}$ is defined as

$t_{\perp}=R_{c}\left<v_{\perp}^{-1}\right>\enspace.$ (90)

See Appendix D for a derivation of this result and an explicit expression for
$\kappa$.

To determine $t_{m}$, we have to take into account the fact that the echo
spreads laterally at velocity $\left<v_{\perp}^{-1}\right>^{-1}$. If $t_{\rm
off}<t_{\perp}$, the echo will continue to be received by the dish for times
even larger than $t_{\rm off}$, and it will drop drastically for
$t>t_{\perp}$. In this case $t_{m}$ is given roughly by $t_{\perp}$. On the
other hand, when $t_{\rm off}>t_{\perp}$, 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_{\rm
off}$. To summarize this reasoning, the measurement time is given roughly by
the maximum between $t_{\rm off}$ and $t_{\perp}$.

To use Eqs. (88) and (89), we smooth $t_{e}$, $t_{m}$ and $B$ as

$\displaystyle t_{e}$ $\displaystyle={\kappa\,t_{\perp}t_{\rm
off}\over\sqrt{\kappa^{2}\,t_{\perp}^{2}+t_{\rm off}^{2}}}$ $\displaystyle
t_{m}$ $\displaystyle=\sqrt{t_{\rm off}^{2}+t_{\perp}^{2}}$ $\displaystyle B$
$\displaystyle={1\over 4\pi}{\delta\omega\delta p_{z}\over\Delta}\enspace.$
(91)

Figure 7: Beam height at $t_{\perp}$ for $R=R_{c}/2$ as a function of $R_{c}$
and $\langle v^{-1}\rangle^{-1}$. The horizontal blue lines mark the value of
$\langle v^{-1}\rangle^{-1}$ for minimal axion transverse velocity,
specifically for $v_{p}=0$ and $v_{p}=5$ km/s for the isothermal and caustic
ring models, respectively. For the caustic ring model we set $\delta
v_{\perp}\simeq 0$, while for the isothermal sphere model with $\delta
v_{\perp}=270\sqrt{2/3}$ km/s. The black lines indicate the height of the
lower end of the ionosphere’s layers: the D-layer starting at a height of 60
km, the E-layer starting at 90 km and the F-layer starting at 150 km.

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_{\rm off}$, the
energy spent to cover a factor of two in axion mass is $E\approx mP_{0}t_{\rm
off}/(2\delta\omega)$. In terms of $E$, Eq. (88) is written as

$s/n={\pi\over 4\sqrt{2}}{g^{2}\rho\,t_{\perp}E\over mT_{n}\sqrt{\delta
p_{z}}}\sqrt{\delta\omega\over\Delta\,t_{m}}\enspace.$ (92)

We see that the signal to noise ratio is maximized for $t_{\rm off}\rightarrow
0$. However, as $t_{\rm off}$ is limited by the uncertainty relation
$\delta\omega t_{\rm off}\geq 1/2$, Eq. (92) will be evaluated at $t_{\rm
off}=(2\delta\omega)^{-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. (92) assuming
$E=10\,\text{MW\,year}$, $s/n=5$, $T_{n}=20\,\text{K}$, $R=50\,\text{m}$,
$R_{c}=100\,\text{m}$ and $\delta\omega=\delta 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 $\delta v_{\perp}\simeq 220\,\text{km}/\text{s}$. For the caustic
ring model, we took $v_{p}=5\,\text{km/s}$, in agreement with the discussion
in VI.1. 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_{\rm off}$ is given by $t_{\rm
off}\simeq 2\delta\omega t_{T}/m$. It is not difficult to show that, with this
constraint, the signal to noise is maximized for $\delta\omega={m\over
2}\sqrt{\delta v_{z}t_{\perp}\over t_{T}}$ and $t_{\rm off}=\sqrt{\delta
v_{z}t_{\perp}t_{T}}$. To have an idea of the order of magnitude, for
$R_{c}=100\,\text{m}$ and $t_{T}=1\,\text{year}$, $\delta\omega$ is
$844\,\text{Hz}\,\left(m\over 10^{-5}\text{eV}\right)$ for the isothermal
model and $89.2\,\text{Hz}\,\left(m\over 10^{-5}\text{eV}\right)$ for the
caustic ring model. On the other hand, $t_{\rm off}$ is about $1.12\,\text{s}$
for the isothermal model and about $370\,\text{ms}$ for the caustic ring
model.

The dashed blue lines of Fig. 6 correspond to the sensitivities using a power
source of $P_{0}=10\text{MW}$ 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_{\perp}$,
in Figure 7 we show the beam height at after a time $t_{\perp}$ (assuming
$t_{\rm off}>t_{\perp}$) as a function of the radius of the receiver $R_{c}$
and $\langle v^{-1}\rangle^{-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. Arza:2019nta . 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. Arza:2019nta  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.

Acknowledgements: We would like to thank Gianni Bernardi, Richard Bradley,
Fritz Caspers, Jordi Miralda, Alessio Notari, Marco Regis, Thomas Schwetz-
Mangold, Pierre Sikivie, Lei Wu, Qiaoli Yang and Bin Zhu for useful
discussions. This project has received support from the European Union’s
Horizon 2020 research and innovation programme under the Marie Sklodowska-
Curie grant agreement No 860881-HIDDeN.

##  Appendix A Green’s function for the 1D computation

The retarded Green’s function $G(t,z)$ associated with the operator
$\partial_{t}-\partial_{z}$ can be found by solving the equation

$(\partial_{t}-\partial_{z})G(t,z)=\delta(t)\delta(z)\enspace.$ (93)

To solve Eq. (93), we expand $G(t,z)$ in Fourier space as

$G(t,z)={1\over(2\pi)^{2}}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\,dk\,d\omega\,e^{-i(\omega
t-kz)}\hat{G}(\omega,k)\enspace.$ (94)

Plugging (94) into (93), we get

$\hat{G}(\omega,k)={i\over\omega+k}\enspace,$ (95)

and therefore

$G(t,z)={i\over(2\pi)^{2}}\int_{-\infty}^{\infty}dk\,e^{ikz}\int_{-\infty}^{\infty}\,d\omega{e^{-i\omega
t}\over\omega+k}\enspace.$ (96)

To find the retarded solution, we want $G(t,z)$ to vanish for $t<0$. To do
that, we make the shift $\omega\rightarrow\omega+i\epsilon$ in the denominator
of the integrand in (96). Using the Residue theorem, we find $G(t,z)=0$ for
$t<0$ and

$\displaystyle G(t,z)$ $\displaystyle=$
$\displaystyle{i\over(2\pi)^{2}}\int_{-\infty}^{\infty}dk\,e^{ikz}(-2\pi
i)e^{ikt}$ $\displaystyle=$ $\displaystyle{1\over
2\pi}\int_{-\infty}^{\infty}dk\,e^{ik(z+t)}$ $\displaystyle=$
$\displaystyle\delta(t+z)\enspace,$

for $t>0$. So we have

$G(t,z)=\Theta(t)\delta(t+z)\enspace.$ (97)

##  Appendix B Green’s function for the paraxial equation

The Green’s function $G(t,\vec{x}\,)$ satisfies the equation

$\left(2i\Omega(\partial_{t}-\partial_{z})-\nabla_{\perp}^{2}\right)G(t,\vec{x})=\delta^{3}(\vec{x})\delta(t)\enspace.$
(98)

To solve it, we write $G(t,\vec{x})$ as the Fourier expansion

$G(t,\vec{x})={1\over(2\pi)^{4}}\int
d^{2}k_{\perp}e^{i\vec{k}_{\perp}\cdot\vec{x}_{\perp}}\int dq\ e^{iqz}\int
d\lambda\ e^{-i\lambda t}\tilde{G}(k_{\perp},q,\lambda)\enspace.$ (99)

Plugging (99) into (98), we find

$\tilde{G}(k_{\perp},q,\lambda)={1\over
2\Omega}\left({1\over\lambda+q+{k_{\perp}^{2}\over 2\Omega}}\right)\enspace.$
(100)

As we look for a solution that vanishes at $t\rightarrow-\infty$, ansatz (55)
suggest the shift $\Omega\rightarrow\Omega-i\eta$. We can see easily that
(100) has a pole at $\lambda=-q-k_{\perp}^{2}/(2\Omega)-i\tilde{\eta}$. We
perform the integrals in (100) getting

$\displaystyle G(t,\vec{x})$ $\displaystyle=$
$\displaystyle{1\over(2\pi)^{4}}{1\over 2\Omega}\int
d^{2}k_{\perp}e^{i\vec{k}_{\perp}\cdot\vec{x}_{\perp}}\int dq\ e^{iqz}\int
d\lambda\,{e^{-i\lambda t}\over\lambda+q+{k_{\perp}^{2}\over
2\Omega}+i\tilde{\eta}}$ (101) $\displaystyle=$
$\displaystyle{1\over(2\pi)^{4}}{1\over 2\Omega}(-2\pi i)\Theta(t)\int
d^{2}k_{\perp}e^{i\vec{k}_{\perp}\cdot\vec{x}_{\perp}}e^{i{k_{\perp}^{2}\over
2\Omega}t}\int dq\ e^{iq(z+t)}$ $\displaystyle=$
$\displaystyle{1\over(2\pi)^{3}}{1\over 2\Omega}(-2\pi
i)\Theta(t)\delta(t+z)\int
d^{2}k_{\perp}e^{i\vec{k}_{\perp}\cdot\vec{x}_{\perp}}e^{i{k_{\perp}^{2}\over
2\Omega}t}$ $\displaystyle=$ $\displaystyle{1\over(2\pi)^{3}}{1\over
2\Omega}(-2\pi
i)\Theta(t)\delta(t+z)\int_{0}^{\infty}dk_{\perp}k_{\perp}e^{i{k_{\perp}^{2}\over
2\Omega}t}\int_{0}^{2\pi}d\phi\,e^{ik_{\perp}x_{\perp}\cos(\phi)}$
$\displaystyle=$ $\displaystyle{1\over(2\pi)^{2}}{1\over 2\Omega}(-2\pi
i)\Theta(t)\delta(t+z)\int_{0}^{\infty}dk_{\perp}k_{\perp}e^{i{k_{\perp}^{2}\over
2\Omega}t}J_{0}(k_{\perp}x_{\perp})$ $\displaystyle=$
$\displaystyle-{\Theta(-z)\over 4\pi z}\delta(t+z)e^{{i\Omega\over
2z}x_{\perp}^{2}}\enspace.$

##  Appendix C Approximated formula for $\delta p_{z}\,t\gg 1$

Defining

$\xi(z^{\prime},\vec{x}_{\perp},\vec{p}_{\perp},\Omega)={1\over 4\pi
z^{\prime}}\int
d^{2}x_{\perp}^{\prime}\eta(\omega,\vec{x}\,^{\prime})e^{-i\vec{p}_{\perp}\cdot\vec{x}\,^{\prime}_{\perp}}\,e^{-{i\Omega\over
2z^{\prime}}|\vec{x}_{\perp}-\vec{x}\,^{\prime}_{\perp}|^{2}}\enspace,$ (102)

the function $J(t,\vec{x},\Omega)$ defined in Eq. (59), at $z=0$, can be
written as

$J(t,\vec{x}_{\perp},\Omega)=\sqrt{V\over(2\pi)^{3}}\int
d^{3}pf_{a}(\vec{p})^{*}e^{-i\phi_{\vec{p}}}\int dz^{\prime}e^{i\epsilon
z^{\prime}}\xi(z^{\prime},\vec{x}_{\perp},\vec{p}_{\perp},\Omega)\enspace.$
(103)

The general formula for the spectral intensity Eq. (61) averaged over random
phases gives us

${\partial I(t,\vec{x}_{\perp})\over\partial\Omega}={1\over
2}g^{2}\rho{dI_{0}\over d\omega}\Omega^{2}\int dz^{\prime}\int
dz^{\prime\prime}\left<\xi(z^{\prime})\xi(z^{\prime\prime})^{*}\right>\int
dp_{z}|f_{a}(p_{z})|^{2}e^{i\epsilon(z^{\prime}-z^{\prime\prime})}\enspace.$
(104)

Using Eq. (33), the integral in $p_{z}$ is easily evaluated, it gives

$\int
dp_{z}|f_{a}(p_{z})|^{2}\,e^{i\epsilon(z^{\prime}-z^{\prime\prime})}=e^{i(2\omega-m-\left<p_{z}\right>)(z^{\prime}-z^{\prime\prime})-\delta
p_{z}^{2}(z^{\prime}-z^{\prime\prime})^{2}/2}\enspace.$ (105)

At resonance, i.e., when $\omega=\omega_{*}$, we can use approximation Eq.
(35) in the limit $\delta p_{z}t\gg 1$, getting

${\partial I(t,\vec{x}_{\perp})\over\partial\Omega_{*}}={1\over
8}\sqrt{\pi\over 2}{g^{2}\rho\over\delta p_{z}}\,{\cal
T}(t,\vec{x}_{\perp})\,{dI_{0}\over d\omega_{*}}\enspace,$ (106)

where ${\cal T}(t,\vec{x}_{\perp})$ is defined as

$\displaystyle{\cal T}(t,\vec{x}_{\perp})$ $\displaystyle=$ $\displaystyle
8\omega_{*}^{2}\int d^{2}p_{\perp}|f_{a}(\vec{p}_{\perp})|^{2}\int
dz^{\prime}\,|\xi(z^{\prime},\vec{x}_{\perp},\vec{p}_{\perp},\omega_{*})|^{2}$
(107) $\displaystyle=$ $\displaystyle 8\omega_{*}^{2}\,\left<\int
dz^{\prime}\,|\xi(z^{\prime},\vec{x}_{\perp},\vec{p}_{\perp},\omega_{*})|^{2}\right>\enspace.$

To calculate the total intensity, we must integrate Eq. (104) also over
frequency, so the knowledge of $dI_{0}/d\omega$ is required. We will use the
simple Gaussian distribution defined in Eq. (65). Now we have to integrate Eq.
(105) as

$\int d\omega{dI_{0}\over d\omega}\int
dp_{z}|f_{a}(p_{z})|^{2}\,e^{i\epsilon(z^{\prime}-z^{\prime\prime})}=e^{i(2\bar{\omega}-m-\left<p_{z}\right>)(z^{\prime}-z^{\prime\prime})-2\Delta^{2}(z^{\prime}-z^{\prime\prime})^{2}/2}\enspace,$
(108)

where $\Delta=\sqrt{\delta\omega^{2}+\delta p_{z}^{2}/4}$. The function
$\Omega^{2}\left<\xi(\Omega)\xi(\Omega)^{*}\right>$ was simply evaluated at
$\bar{\omega}$ since corrections in $p_{z}$ and $\delta\omega$ are small.

If the peaks of the momentum and frequency distributions satisfy the resonance
condition for stimulated axion decay, i.e., $\bar{\omega}=\omega_{*}$, then
Eq. (108) can be approximated as
$\sqrt{\pi/2}\delta(z^{\prime}-z^{\prime\prime})/\Delta$ in the limit $\Delta
t\gg 1$. The total intensity is found to be

$I(t,\vec{x}_{\perp})={1\over 16}\sqrt{\pi\over
2}{g^{2}\rho\over\Delta}\,{\cal T}(t,\vec{x}_{\perp})\,I_{0}\enspace.$ (109)

Finally, the power $P_{c}$ collected over a arbitrary surface $S$ can be
calculated integrating Eq. (109) over $S$. Defining the outgoing beam power
$P_{0}=I_{0}S_{0}$, where $S_{0}$ is the effective cross-section of the
outgoing beam at $t=0$, we have

$P_{c}(t)={1\over 16}\sqrt{\pi\over 2}{g^{2}\rho\over\Delta}\,P_{0}\,{1\over
S_{0}}\int_{S}d^{2}x_{\perp}{\cal T}(t,\vec{x}_{\perp})\enspace.$ (110)

##  Appendix D Details for the power collected by a receiving dish

In this Appendix, we present the calculation of the power collected by a dish
located on the plane $z=0$, concentric to the emission surface. We start
writing Eq. (66) as

$P_{c}={1\over 16}\sqrt{\pi\over 2}{g^{2}\rho\over\Delta}\,\left<{\cal
T}_{c}(t,\vec{v}_{\perp})\right>\,P_{0}\enspace,$ (111)

where

${\cal T}_{c}(t,\vec{v}_{\perp})={1\over S_{0}}\int_{S}d^{2}x_{\perp}\,{\cal
T}(t,\vec{x}_{\perp},\vec{v}_{\perp})\enspace.$ (112)

Here $S_{0}=\pi R^{2}$ and $S$ is the integration domain that corresponds to
the receiving dish surface. From now on, we will assume a circular cross-
sectional area with radius $R_{c}$ for the receiving dish.

We want to analyze the cases in which the dish collects a steady-state signal.
To find them, first notice that as the echo spreads at velocity $\sim
v_{\perp}$, it takes a time $R_{c}/v_{\perp}$ to spread out of the dish. If
$t_{\rm off}<R_{c}/v_{\perp}$, i.e., when the outgoing beam is turned off
before the echo spreading reaches the dish boundaries, a steady-state signal
is received for a time defined by $t_{\rm off}\leq t\leq R_{c}/v_{\perp}$. On
the other hand, if $t_{\rm off}>R_{c}/v_{\perp}$, i.e. after the echo passes
the dish boundary, a steady-state signal is found for $R_{c}/v_{\perp}\leq
t\leq t_{\rm off}$.

In the first case, when $t_{\rm off}\leq t\leq R_{c}/v_{\perp}$, as all the
echo lies inside the dish, the integral Eq. (112) can be performed over a
infinite space. Defining $x=r\cos(\varphi)$ and $y=r\sin(\varphi)$ and taking
the lower part of Eq. (70), we have

$\displaystyle{\cal T}_{c}(\vec{v}_{\perp})$ $\displaystyle=$
$\displaystyle{1\over S}{\sqrt{\pi}\over 2}{R\over
v_{\perp}}\int_{-\infty}^{\infty}dy\,e^{-{y^{2}\over
R^{2}}}\int_{-\infty}^{\infty}dx\left(\operatorname{erf}\left(v_{\perp}t-x\over
R\right)-\operatorname{erf}\left(v_{\perp}(t-t_{\rm off})-x\over
R\right)\right)$ (113) $\displaystyle=$ $\displaystyle{1\over
S}{\sqrt{\pi}\over 2}{R\over v_{\perp}}\sqrt{\pi}R\,2v_{\perp}t_{\rm off}$
$\displaystyle=$ $\displaystyle t_{\rm off}\enspace.$

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_{\perp}\leq t\leq t_{\rm off}$, we can set
the second erf function that appears in Eq. (70) equal to $1$, since its
argument is always large. For ${\cal T}_{c}$, we find

$\displaystyle{\cal T}_{c}(\vec{v}_{\perp})$ $\displaystyle=$
$\displaystyle{1\over S}{\sqrt{\pi}\over 2}{R\over
v_{\perp}}\int_{0}^{R_{c}}dr\,r\int_{0}^{2\pi}d\varphi\,e^{-{r^{2}\over
R^{2}}\sin(\varphi)^{2}}\left(\operatorname{erf}\left(r\cos(\varphi)\over
R\right)+1\right)\enspace,$ (114) $\displaystyle=$ $\displaystyle{1\over
S}{\sqrt{\pi}\over 2}{R\over
v_{\perp}}\int_{0}^{R_{c}}dr\,r\,2\pi\,e^{-{r^{2}\over
2R^{2}}}I_{0}\left(r^{2}\over 2R^{2}\right)$ $\displaystyle=$
$\displaystyle{1\over S}{\sqrt{\pi}\over 2}{R\over
v_{\perp}}2\pi{R_{c}^{2}\over 2}\,e^{-{R_{c}^{2}\over
2R^{2}}}\left(I_{0}\left(R_{c}^{2}\over
2R^{2}\right)+I_{1}\left(R_{c}^{2}\over 2R^{2}\right)\right)$ $\displaystyle=$
$\displaystyle\kappa{R_{c}\over v_{\perp}}$ (115)

where $\kappa$ is defined as

$\kappa={\sqrt{\pi}\over 2}{R_{c}\over R}\,e^{-{R_{c}^{2}\over
2R^{2}}}\left(I_{0}\left(R_{c}^{2}\over
2R^{2}\right)+I_{1}\left(R_{c}^{2}\over 2R^{2}\right)\right)\enspace.$ (116)

Notice that $\kappa\rightarrow 1$ for $R_{c}/R\rightarrow\infty$.

With these results, the steady signal power can be written as

$P_{c}={1\over 16}\sqrt{\pi\over
2}{g^{2}\rho\over\Delta}t_{e}\,P_{0}\enspace,$ (117)

where $t_{e}$ takes the value $t_{\rm off}$ for $t_{\rm off}<R_{c}/v_{\perp}$
and the value $\kappa\,R_{c}\left<v_{\perp}^{-1}\right>$ for $t_{\rm
off}>R_{c}/v_{\perp}$. Notice that if $R_{c}>R$, $\kappa$ is of order $1$, so
we can write the approximate result

$P_{c}={1\over 16}\sqrt{\pi\over 2}{g^{2}\rho\over\Delta}\min\left(t_{\rm
off},\kappa R_{c}\left<v_{\perp}^{-1}\right>\right)P_{0}\enspace.$ (118)

## References

  * (1) R. Peccei and H. R. Quinn, CP Conservation in the Presence of Instantons,  Phys. Rev. Lett. 38 (1977) 1440–1443.
  * (2) S. Weinberg, A New Light Boson?,  Phys. Rev. Lett. 40 (1978) 223–226.
  * (3) F. Wilczek, Problem of Strong $P$ and $T$ Invariance in the Presence of Instantons,  Phys. Rev. Lett. 40 (1978) 279–282.
  * (4) P. Svrcek and E. Witten, Axions In String Theory,  JHEP 06 (2006) 051, [hep-th/0605206].
  * (5) J. Preskill, M. B. Wise, and F. Wilczek, Cosmology of the Invisible Axion,  Phys. Lett. B 120 (1983) 127–132.
  * (6) L. Abbott and P. Sikivie, A Cosmological Bound on the Invisible Axion,  Phys. Lett. B 120 (1983) 133–136.
  * (7) M. Dine and W. Fischler, The Not So Harmless Axion,  Phys. Lett. B 120 (1983) 137–141.
  * (8) P. Arias, D. Cadamuro, M. Goodsell, J. Jaeckel, J. Redondo, et al., WISPy Cold Dark Matter,  JCAP 06 (2012) 013, [1201.5902].
  * (9) P. Sikivie, Invisible Axion Search Methods,  Rev. Mod. Phys. 93 (2021), no. 1 015004, [2003.02206].
  * (10) I. G. Irastorza and J. Redondo, New experimental approaches in the search for axion-like particles,  Prog. Part. Nucl. Phys. 102 (2018) 89–159, [1801.08127].
  * (11) P. Sikivie, Experimental Tests of the Invisible Axion,  Phys. Rev. Lett. 51 (1983) 1415–1417. [Erratum: Phys.Rev.Lett. 52, 695 (1984)].
  * (12) M. P. Hertzberg and E. D. Schiappacasse, Dark Matter Axion Clump Resonance of Photons,  JCAP 11 (2018) 004, [1805.00430].
  * (13) A. Caputo, C. P. n. Garay, and S. J. Witte, Looking for Axion Dark Matter in Dwarf Spheroidals,  Phys. Rev. D 98 (2018), no. 8 083024, [1805.08780]. [Erratum: Phys.Rev.D 99, 089901 (2019)].
  * (14) A. Arza, Photon enhancement in a homogeneous axion dark matter background,  Eur. Phys. J. C 79 (2019), no. 3 250, [1810.03722].
  * (15) A. Caputo, M. Regis, M. Taoso, and S. J. Witte, Detecting the Stimulated Decay of Axions at RadioFrequencies,  JCAP 03 (2019) 027, [1811.08436].
  * (16) Z. Wang, L. Shao, and L.-X. Li, Resonant instability of axionic dark matter clumps,  JCAP 07 (2020) 038, [2002.09144].
  * (17) A. Arza, T. Schwetz, and E. Todarello, How to suppress exponential growth—on the parametric resonance of photons in an axion background,  JCAP 10 (2020) 013, [2004.01669].
  * (18) D. Levkov, A. Panin, and I. Tkachev, Radio-emission of axion stars,  Phys. Rev. D 102 (2020), no. 2 023501, [2004.05179].
  * (19) A. Arza and P. Sikivie, Production and detection of an axion dark matter echo,  Phys. Rev. Lett. 123 (2019), no. 13 131804, [1902.00114].
  * (20) O. Ghosh, J. Salvado, and J. Miralda-Escudé, Axion Gegenschein: Probing Back-scattering of Astrophysical Radio Sources Induced by Dark Matter, 2008.02729.
  * (21) J. W. Foster, N. L. Rodd, and B. R. Safdi, Revealing the Dark Matter Halo with Axion Direct Detection,  Phys. Rev. D 97 (2018), no. 12 123006, [1711.10489].
  * (22) L. Hui, A. Joyce, M. J. Landry, and X. Li, Vortices and waves in light dark matter,  JCAP 01 (2021) 011, [2004.01188].
  * (23) G. P. Centers et al., Stochastic fluctuations of bosonic dark matter, 1905.13650.
  * (24) M. S. Turner, Cosmic and Local Mass Density of Invisible Axions,  Phys. Rev. D 33 (1986) 889–896.
  * (25) L. Duffy and P. Sikivie, The Caustic Ring Model of the Milky Way Halo,  Phys. Rev. D 78 (2008) 063508, [0805.4556].
  * (26) P. Sikivie, Evidence for ring caustics in the Milky Way,  Phys. Lett. B 567 (2003) 1–8, [astro-ph/0109296].
  * (27) N. Banik, A. J. Christopherson, P. Sikivie, and E. M. Todarello, New astrophysical bounds on ultralight axionlike particles,  Phys. Rev. D 95 (2017), no. 4 043542, [1701.04573].
  * (28) S. S. Chakrabarty, Y. Han, A. Gonzalez, and P. Sikivie, Implications of triangular features in the Gaia skymap for the Caustic Ring Model of the Milky Way halo, 2007.10509.
  * (29) E. Hecht, Optics, 3rd ed. Addison-Wesley, 1998.
  * (30) A. Natarajan and P. Sikivie, The inner caustics of cold dark matter halos,  Phys. Rev. D 73 (2006) 023510, [astro-ph/0510743].
  * (31) N. Banik and P. Sikivie, Evolution of Velocity Dispersion along Cold Collisionless Flows,  Phys. Rev. D 93 (2016), no. 10 103509, [1511.05947].
  * (32) S. De Panfilis, A. C. Melissinos, B. E. Moskowitz, J. T. Rogers, Y. K. Semertzidis, et al., Limits on the Abundance and Coupling of Cosmic Axions at 4.5-Microev $<$ m(a) $<$ 5.0-Microev,  Phys. Rev. Lett. 59 (1987) 839.
  * (33) C. Hagmann, P. Sikivie, N. S. Sullivan, and D. B. Tanner, Results from a search for cosmic axions,  Phys. Rev. D 42 (1990) 1297–1300.
  * (34) ADMX, N. Du et al., A Search for Invisible Axion Dark Matter with the Axion Dark Matter Experiment,  Phys. Rev. Lett. 120 (2018), no. 15 151301, [1804.05750].
  * (35) ADMX, T. Braine et al., Extended Search for the Invisible Axion with the Axion Dark Matter Experiment,  Phys. Rev. Lett. 124 (2020), no. 10 101303, [1910.08638].
  * (36) ADMX, C. Boutan et al., Piezoelectrically Tuned Multimode Cavity Search for Axion Dark Matter,  Phys. Rev. Lett. 121 (2018), no. 26 261302, [1901.00920].
  * (37) S. Lee, S. Ahn, J. Choi, B. R. Ko, and Y. K. Semertzidis, Axion Dark Matter Search around 6.7 $\mu$eV,  Phys. Rev. Lett. 124 (2020), no. 10 101802, [2001.05102].
  * (38) J. Jeong, S. Youn, S. Bae, J. Kim, T. Seong, et al., Search for Invisible Axion Dark Matter with a Multiple-Cell Haloscope,  Phys. Rev. Lett. 125 (2020), no. 22 221302, [2008.10141].
  * (39) CAPP, O. Kwon et al., First Results from an Axion Haloscope at CAPP around 10.7 $\mu$eV,  Phys. Rev. Lett. 126 (2021), no. 19 191802, [2012.10764].
  * (40) HAYSTAC, L. Zhong et al., Results from phase 1 of the HAYSTAC microwave cavity axion experiment,  Phys. Rev. D 97 (2018), no. 9 092001, [1803.03690].
  * (41) HAYSTAC, K. M. Backes et al., A quantum-enhanced search for dark matter axions,  Nature 590 (2021), no. 7845 238–242, [2008.01853].
  * (42) D. Alesini et al., Galactic axions search with a superconducting resonant cavity,  Phys. Rev. D 99 (2019), no. 10 101101, [1903.06547].
  * (43) D. Alesini et al., Search for invisible axion dark matter of mass m${}_{a}=43~{}\mu$eV with the QUAX–$a\gamma$ experiment,  Phys. Rev. D 103 (2021), no. 10 102004, [2012.09498].
  * (44) B. T. McAllister, G. Flower, E. N. Ivanov, M. Goryachev, J. Bourhill, et al., The ORGAN Experiment: An axion haloscope above 15 GHz,  Phys. Dark Univ. 18 (2017) 67–72, [1706.00209].
  * (45) CAST, A. A. Melcón et al., First results of the CAST-RADES haloscope search for axions at 34.67 $\mu$eV, 2104.13798.
  * (46) CAST, V. Anastassopoulos et al., New CAST Limit on the Axion-Photon Interaction,  Nature Phys. 13 (2017) 584–590, [1705.02290].