Axion Echos from the Supernova Graveyard

Stimulated decays of axion dark matter, triggered by a source in the sky, could produce a photon flux along the continuation of the line of sight, pointing backward to the source. The strength of this so-called axion"echo"signal depends on the entire history of the source and could still be strong from sources that are dim today but had a large flux density in the past, such as supernova remnants (SNRs). This echo signal turns out to be most observable in the radio band. We study the sensitivity of radio telescopes such as the Square Kilometer Array (SKA) to echo signals generated by SNRs that have already been observed, and show that SKA could reach axion-photon couplings of order $g_{a\gamma\gamma} \sim \mathcal{O}(10^{-11}) \,\mathrm{GeV}^{-1}$ for axion masses $m_a \lesssim 10^{-5}\;\mathrm{eV}$. In addition, we show projections of the detection reach for signals coming from old SNRs and from newly born supernovae that could be detected in the future. Intriguingly, an observable echo signal could come from old"ghost"SNRs which were very bright in the past but are now so dim that they haven't been observed.


Introduction
An axionlike particle, a, is a periodic pseudoscalar field, a ∼ = a + 2πf a with f a the decay constant. The discrete shift symmetry acting on the field is its defining feature. It protects the small mass of the axion, m a , and imposes strong constraints on the forms of its couplings. Axions arise ubiquitously from top-down theory such as string models [1], and low-energy phenomenological models [2][3][4][5][6][7][8][9]. They serve as one of the most motivated scenarios of feebly coupled light particles beyond the standard model. In particular, an important observational target is the axion's coupling to photons in the standard model, which is usually parametrized as − g aγγ 4 aF µνF µν = g aγγ aE · B, (1.1) where the coupling g aγγ has energy dimension −1, inversely proportional to f a ; F µν is the field strength of electromagnetic U (1) em ;F µν = 1 2 µναβ F αβ is the dual field strength, and E, B are the corresponding electric and magnetic fields. Another attraction of an axion is that it could serve as a cold dark matter candidate [10][11][12]. The simplest canonical mechanism is misalignment: an axion is initially frozen at some random place in its field space due to the Hubble friction. When the Hubble expansion rate drops around its mass, the axion starts to oscillate around the minimum of its potential. The coherent oscillations redshift as cold matter and store energy in the axion field.
Given its theoretical importance, there has been a rapidly growing interest in searching for axions and their couplings over the years, cf.  to name a few. In this article, we will study a novel technique to probe the coupling of axion dark matter to photons, which has been underexplored, that is, searching for an "echo" 1 signal from stimulated axion dark matter decays. 2 The basic idea is depicted in Fig. 1: photons from a source traverse the axion dark matter halo. If the photon energy matches half of the axion mass, it induces stimulated decays of axion dark matter into two photons. Due to momentum conservation, these two photons travel in opposite directions with energies equal to half of the axion mass, m a /2. An observer could receive two fluxes of photons from opposite directions. One flux is along the line of sight from the source to the observer, which is a superposition of the continuum emission of the source (blue wavy line) and the line emission from axion decays with an angular frequency set by m a /2 (purple wavy line). Since the line emission from the stimulated decays is a much fainter signal, it is challenging to isolate it from the bright source continuum background. The other flux is the so-called "echo" signal also from stimulated axion decays, which is along the continuation of the line of sight to the direction opposite to the source. This one could potentially be a clean signal if there is no bright source in the opposite direction.
S O a a "echo" The echo technique was first applied and developed in Refs. [37,38], with the source being a powerful radio beam shooting from the Earth. Given the experimental challenges, it will be helpful to identify natural sources in the universe. So far the only astrophysical source that has been examined in the literature to trigger an echo signal is Cygnus A (Cyg A), an extragalactic active galactic nucleus [36]. In this work, we point out for the first time that the echo signal relies on the historical luminosity of the photon source, while the signal of photons from stimulated decays of axion dark matter that travel along with the source photon as studied in Ref. [39][40][41] require the brightest radio sources today. This opens up the possibility to probe axion DM with time-varying radio sources that may be dim today but were once very bright, such as supernova remnants (SNRs).
The analyses of SNRs will also be drastically different from Cyg A, and any other relatively stable sources over the timescale for light to travel across the dark matter halo. An SNR, as a time-dependent source, typically starts with a huge photon flux from the explosion of a supernova and undergoes several different phases before eventually merging with the interstellar medium. In this case, the echo signal could be induced by the flux from the SNR at time much earlier than that at which the observation is made. For a more intuitive understanding, one could consider an SNR emitting a series of pulses over time, as shown in the vertical direction of the cartoon in Fig. 2. The pulse emitted at earlier times could pass the observer and travel further before triggering the decay of axion dark matter and generating an echo photon, while the later-emitted pulse travels a shorter distance past the observer before triggering the axion stimulated decays, as demonstrated by the horizontal direction in Fig. 2. These echo photons could travel back and arrive at the observer at the same time t o . In other words, the observed signal is a collection of the echo photons generated by the pulses from the source at different times. More concretely, as we will show in a more rigorous way, the signal is an integration over the history of the SNR. It turns out that even a dim source today could still lead to an observable echo signal, since it was brighter earlier.
Throughout the paper, we work with natural units. Our paper is organized as follows. In Sec. 2, we derive our master formula to compute the axion echo signal from a time-dependent source, and comment on different effects of relative motions between SNR, the observer, and the dark matter halo. In Sec. 3, we review the basics and time evolution of SNRs that we will use in our analysis. In Sec. 4, we discuss the radio telescope, Square Kilometer Array, used for detecting the echo signal and the calculation of the signal to noise ratio. We present our results in Sec. 5 and conclude in Sec. 6. Finally, in Appendix A we check the uncertainties in the sensitivity of Square Kilometer Array to axion echoes trigged by the most promising SNR, G39.7-2.0. We have made publicly available the Python code we developed to aid our calculations; it can be found at github.com/ManuelBuenAbad/snr ghosts.

Echo signal
In this section, we will first derive the key formula to compute the flux density (also referred to as spectral irradiance in the literature), S ν , of the echo signal from a time-dependent source such as SNRs. We will then comment on the effects of relative motions between the source, the observer, and dark matter halo.

Flux density of echo signal
Since we will consider SNRs in our galaxy, we ignore the expansion of the universe. The particle number density n is related to its phase space density f (p) as dn = DPf , with DP ≡ gd 3 p/(2π) 3 and g the degrees of freedom of the particle. In our scenario, the evolution of the photon phase space density f γ is governed by the Boltzmann equation: where C[f γ ] is the collision term for photons. Taking into account both the stimulated decays of axion dark matter and the inverse process, a(p a ) → γ(p 1 ) + γ(p 2 ) and γ(p 1 ) + γ(p 2 ) → a(p a ), the collision term is given by: where E's (p's) are the energies (momenta) of the particles involved (subscripts a, 1 and 2 denote axion, photon 1 and 2 respectively); f 's are the phase space densities; and |M| 2 is the averaged matrix element squared of the processes. It can be shown, for the axion and photon occupation numbers in our study (which are functions of the dark matter density, the luminosity of the source, and the axion mass and photon energy), that the backward process ∝ f 1 f 2 is negligible compared to that of the forward process ∝ f a (1 + f 1 + f 2 ). Furthermore, in a leading-order approximation, we take the axion dark matter to be perfectly cold and thus we can ignore its motion, i.e. f a = n a (2π) 3 δ 3 (p a ). For photons coming from directional sources, the phase space density peaks in a given direction, f 1 = f γ (|p 1 |)h(p 1 ), with the direction function h(p 1 ) in general anisotropic andp 1 being the direction of the original photon propagation. For photons from point sources such as SNRs, h(p 1 ) ≈ δ 2 (Ω 1p1 + Ω * n * ), where Ω 1 is the corresponding phase space solid angle ofp 1 , and Ω * is the solid angle of the light source as seen from the Earth, in the directionn * . Said differently, for point sources the phase space solid angle Ω 1 is that of the light source Ω * , and the direction of motion of the photonsp 1 is (obviously) antiparallel to the vectorn * going from the observer to the source.
Finally, due to momentum conservation, integrating Eq. (2.1) and Eq. (2.2) over DP a and DP 2 gives: where the axion (spontaneous) decay width is given by Γ a ≡ g 2 aγγ m 3 a /(64π), ρ a is the dark matter energy density ρ a = n a m a , and E a = m a /2. The three terms in the parentheses correspond to the isotropic spontaneous decay, the stimulated decay photon along the original photon direction (as mentioned above, antiparallel ton * ), and the echo photon going in the opposite direction (i.e. , parallel ton * , pointing toward the source), respectively. In what follows, we focus on the latter term, the echo signal.
For pointlike sources, as opposed to extended ones, the relevant observable is the flux density S ν , which is related to the phase space density as S ν = E 3 2π 2 dΩ f (E, Ω), with ν = E/2π the photon frequency. It can further be shown that the flux of the echo signal is much smaller than that of the source, i.e. we stay in the perturbative response regime. In other words, the echo stimulated by an echo (as opposed to by the original source itself) is negligible. 3 We can therefore integrate Eq. (2.3) over the solid angle and then over time, in order to obtain the flux density of the echo signal S νa,e observable today, and relate it to the source's flux density S νa,s (t), which is a function of time: Here ν a = E a /(2π) = m a /(4π). Note the crucial retarded time argument in the source's flux density S νa,s (t). This reflects the fact that the portion of the echo emitted by axion dark matter at a distance x along the line of sight (and being observed today, when the source has an age t age ) was produced by light that first passed the observer's location a time 2x ago. Note also that the axion dark matter density is evaluated along the line of sight in the direction opposite to the source, hence the name of "echo". Throughout the rest of this paper, we will assume that the dark matter is distributed according to the Navarro-Frenk-White profile [43], with a scale radius of 20 kpc [44,45] and a local density of 0.4 GeV cm −3 in the solar neighborhood [46,47]. We can immediately see from Eq. (2.4) that, under the assumption of perfectly cold dark matter, the echo signal looks like a line at photon energy equal to E a = m a /2. Let us now relax that assumption. Indeed, the dark matter has a velocity dispersion σ v ≈ 5×10 −4 c (c being the speed of light) [48]. As a result, the echo signal is not exactly a delta function, but a distribution with a finite width of order ∼ σ v E a . Finally, it can be explicitly shown that for dark matter velocity distributions following the standard halo model, the echo signal is approximately distributed as a Gaussian with standard deviation σ v E a .
For narrow signals, the observationally relevant quantity is the flux density averaged over the bandwidth of interest ∆ν. This can be obtained from Eq. (2.4) simply by integrating the finite distribution, and dividing by ∆ν. The result is Here f ∆ is the fraction of the total (integrated) signal flux falling within the bandwidth ∆ν. Anticipating the results from Sec. 4, we point out here that there is an optimal ∆ν (and therefore f ∆ ) that will maximize the ratio of the signal to the noise, with the latter scaling as ∝ √ ∆ν. Therefore, approximating the signal distribution as a Gaussian function, we find that the optimal numbers are [36]: In Eq. (2.5) we have included the optical depth τ . The optical depth due to free-free absorption can be computed as follows [49]: where EM = n 2 e d is the emission measure of electrons along the photons' propagation, and T e is the temperature of electrons. There have been many attempts to model the electron number density in the galaxy, e.g. Ref. [50]. When both the echo signal and the source photon come from a direction away from the galactic center, one could estimate the electron number density with a single scale height h using a simple exponential function n e = n e,0 e −|z|/h . A fit from [50] with pulsar measurements (model M1 therein) gives an electron density in the galactic midplane as n e,0 ≈ 0.015 cm −3 . This translates to an EM: which makes the optical depth negligible for the entire range of frequency we investigate. However, when the echo signal comes from the vicinity of the galactic center, this expression greatly underestimates the EM. From Ref. [51], the emission measure is estimated to be EM ∼ 10 5 cm −6 pc, based on the radio flux from Sgr A * . Using this, we can find the optical depth at the galactic center to be ∼ O(1) at ν ≈ 300 MHz. Note that the emission measure in [51] is based on the radio observation of Sgr A * and the 7 arcmin halo surrounding it. When echo signals come straight from the galactic center, the background is very strong leading to poor sensitivity. In what follows, we do not discuss signals coming directly from the galactic center, therefore we neglect the optical depth.

Motions of source, observer, and dark matter
In the derivation above, we implicitly assume that the source, the observer on Earth, and the dark matter halo are at rest. However, in reality, they all move with respect to each other, and their relative motions could cause a reduction of the echo signal, which we will estimate in this section. For simplicity, we will study one motion at a time. The general motions of the system could be more complicated but the overall effect is a combination of the three effects discussed below. First, we consider the coherent motion of the source. An echo photon points back toward the position the source had at the time when it emitted the photon that triggered the stimulated decays. Since the signal is a stack of the echo photons, the motion of the source causes an aberration of the echo image. During the time relevant for generating the echo signal, the source travels a distance d s . If d s is much below the distance D between the source and the observer, the aberration angle is estimated to be where v s is the speed of the source, for which we picked a characteristic value as a benchmark, considering SNRs in our galaxy. For a SNR with size ∼ θ ab , the image will be blurred and enlarged by an order one factor. On the other hand, the motion of the Earth away from the source could lead to a reduction of the flux. The change in the total echo photon flux received could be estimated as where d o is the distance the observer travels and v o is the observer's velocity. One can see that this is a negligible effect. Finally, we consider the peculiar motion of the dark matter. The peculiar motion could cause the echo photon to move out of the telescope's receiver dish, as shown in panel (a) of Fig. 3. More specifically, due to the peculiar motion of dark matter particles, an echo photon could depart from its original route and move over a distance of order σ v t age ∼ 10 light years, for a SNR of 10 4 years old, in the direction transverse to the line of sight. This is significantly larger than a telescope dish size and the photon could not be received. To overcome this difficulty, we need to consider a wider collecting solid angle, as illustrated in panel (b) of Fig. 3. Let us consider an echo photon that could not hit the detector without peculiar motion (the original route is denoted by the dashed red line). With peculiar motion (represented by the orange dashed line; for simplicity, we assume that the peculiar motion is in the transverse direction), the echo photon will hit the detector along the red solid route. Thus increasing the collecting angle could help collect echo photons nearby to overcome the loss due to the peculiar motion. A straightforward computation based on the geometry in panel (b) tells us that as long as σ v D/x, the new collecting angle we need is given by: where, in the second step, we take the maximum distance, x, which an echo photon being collected today could travel, to be t age /2. Note that the aberration angle θ ab in Eq. (2.9) and δ are of similar size since the rotational speed of a source not far away from the solar system and the dark matter peculiar velocity are of the same order.  For illustration, we show the peculiar motion in the direction transverse to the line of sight. The final route of the echo photon is depicted by the red solid arrow. Right: enlarging the collecting angle to δ to collect echo photons from neighboring points in order to compensate for the loss due to peculiar motions. The blue dot represents the source; the observer is shown as a black dot. More details can be found in the main text.

Shadow of the Earth
While the relative motions of dark matter, source, and the observer can cause the reduction or blurring of the signal, as discussed above, they can also help overcome the effects of the shadow of the Earth. When both the source and the observer are at rest in the dark matter rest frame, the shadow of the observer, the Earth in our case, will form. The shadow could be detrimental to the signal in two ways: (a) stimulated decays of axions located within the Earth's shadow are prevented; (b) the echo photons from axion decays outside the shadow will not pass the Earth. However, neither of the two effects poses a problem due to the motion of the Earth in the dark matter rest frame.
The shadow at any given moment starts from the Earth, extending to the direction opposite to the source. Given that the angular diameter size of the sources are much larger than that of the Earth, the length of the shadow can be estimated as where R ⊕ is the radius of the Earth and θ s is the angular diameter size of the source. Let us take Cassiopeia A (G111.7-2.1 from Green's SNR catalog [52][53][54][55]) as an example. Its associated shadow size is sd ≈ 0.06 au. However, this estimate is obtained by using the current size of the SNR. In the past, the size was smaller and therefore the shadow it cast was larger. Nevertheless, even taking the SNR expansion into account, the length of the shadow barely reached ∼ 1 au when the SNR was around 10 years old. This leads to a negligible reduction of the integral in Eq. (2.5), which has t age ∼ O(10 4 ) year ∼ O(a few) kpc for most cases of interest. Now let us look at the seemingly more fatal challenge from effect (b). If everything is still in the dark matter rest frame, stimulated photons trace back to the source along the exact route of the original source photons. An obvious condition for echo photons to be generated is that the source photons should not be blocked by the Earth. This leads to the simple observation that no echo photons will be captured by a telescope on the Earth regardless of either the size of the shadow or the finite size of the source if the Earth does not move! Fortunately, the Earth, together with the entire solar system, moves with respect to the dark matter halo. It only takes a negligible amount of time for the Earth to move out of its own shadow, t ∼ R ⊕ /v ∼ 21 s, as the solar system revolves around the galactic center with a circular speed v ∼ 10 −3 . Thus effect (b) of the shadow is also not a problem for observing the echo signal.

Supernova remnants
In this section, we will discuss the properties of SNRs, as well as their time evolution, both of which will be used in our analysis. We do not intend to provide a thorough review on this vast subject and will instead refer interested readers to Refs. [52][53][54][55][56][57].
For the frequency range we are interested in, the photon spectrum of a SNR follows a simple power law: where α is the spectral index and E γ is the photon energy. We further choose the reference (or "pivot") frequency to be 1 GHz, following the convention in the literature. The spectral luminosity is then given by: where D is the distance of the source, and spherical symmetry of the source is assumed. From the Green catalog of SNRs [52][53][54][55], one could see that spectral indices are centered around α ≈ 0.5 with considerable scattering, as shown in Fig. 4. We note that among the SNRs in the catalog, most are shell type, with less than 10 filled type SNR and about 20 of the composite type. The latter two types typically have a spectral index smaller than 0.5 but only constitute a small fraction of all the observed SNRs.

Time evolution
An important input for our analysis is the time evolution of the spectral luminosity L ν of SNRs in the radio band, also called radio light curves. While observed SNR light curves exhibit a variety of forms, they typically show an initial exponential rise and subsequent power-law decays. SNR light curves have been the subject of intense study for decades [57][58][59]. Four distinct phases, first proposed in Ref. [58], have been identified: 1. Free expansion phase: in this phase, the swept-up gas mass is much smaller than the supernova mass. Sometimes it is also called the ejecta-dominated phase [60]. During this phase, the shock wave expands with a roughly constant velocity ∼ O(0.01) (in units of the speed of light), leading to r ∝ t, with r the SNR radius, and t the time elapsed after the explosion [56]. The optical depth due to external free-free absorption and synchrotron radiation drops over time, which results in an exponential rise in the luminosity shortly after the explosion [56,59,61]. This phase is estimated to last O(few × 100) years by comparing the initial supernova mass and the mass swept up [56,[60][61][62].

2.
Adiabatic expansion or Sedov-Taylor phase: the energy loss rate is much smaller compared to the expansion rate, so energy conservation approximately holds. The driving force is the pressure behind the shock wave, with surrounding material being pushed away. Hydrodynamic estimates show an expansion as r ∝ t 2/5 in this phase [56]. This phase can last up to ∼ O(few × 10 4 ) years after the explosion [62].
3. Snow plough or radiative phase: after the remnant has dissipated half of its initial energy, the pressure driving the shock wave drops. The shell sweeps with a constant momentum and the expanding radius scales as r ∼ t 1/4 . This phase typically ends around 10 6 years after the explosion [62,63].
4. Dispersion phase: the shock velocity falls below the speed of sound and the density drops to the same as that of the interstellar medium (ISM). It marks the merger of the SNR with the ISM.
Parameter mean (µ) standard deviation (σ) log 10 L ν,pk erg s −1 Hz −1 25. 5 1.5 log 10 t pk days 1.7 0.9 Table 1. Mean and standard deviations of the SNRs peak light curve parameters L ν,pk and t pk in the radio band, describing the free expansion phase; see Eq. (3.3). The parameters are log-normal distributed. While these parameters were derived from observations in a range of 2-10 GHz, they are clustered around 8 GHz. Therefore, we take the L ν,pk values in this table to denote the spectral luminosity at ν = 8 GHz; the values at other frequencies can be deduced following the power law ν −α . For example, for a SNR with spectral index of α = 0.7, its expected peak luminosity at 1 GHz (corresponding to the mean of the distribution) is L 1GHz,pk ≈ 10 26 erg s −1 Hz −1 .
Strictly speaking, one needs to implement full hydrodynamical simulations of the SNR time evolution for a precise computation of its properties. However, in our work, we only intend to provide a proof-of-principle computation of axion dark matter echoes triggered by SNRs. We will thus restrict ourselves to consider the contributions to the echo signal of the first two phases. In addition, we will apply some simple parametrizations of these two phases, which have been used in the literature.
In the first free expansion phase, we follow Ref. [59] and take the following phenomenological model for the light curve: where L ν,pk is the peak luminosity of the SNR, and t pk is the corresponding peak time when L ν,pk is reached. Initially, the exponential rise dominates till the luminosity reaches its peak at t pk . Afterward, the light curve falls off following a power law. While there are large uncertainties associated with SNRs in general, typical value ranges for the peak parameters L ν,pk and t pk can be derived from observations, as demonstrated in Ref. [59]. After analyzing about 300 recent supernovae in the radio frequency 2-10 GHz, the authors of Ref. [59] found that the peak parameters are well described by a log-normal distribution, with means and standard deviations shown in Table 1. However, these measurements cluster around ∼ 8 GHz, with a relatively narrow width (standard deviation of ∼ 2 GHz). We therefore take the L ν,pk values quoted in Table 1 to be those for L 8GHz,pk . Following Eq. (3.2) we can deduce the peak luminosity at other frequencies, particularly 1 GHz, by applying the ν −α scaling. As the free expansion phase comes to an end, the SNR transitions into the adiabatic phase at a time t tran . The spectral luminosity follows a single power law during the adiabatic expansion as: with the adiabatic index [56] Alternatively, for SNRs still in the second phase, we can write the adiabatic evolution in terms of the SNR's spectral luminosity today, L ν,0 , as well as its age: Thus, the light curve of an SNR during its first two phases can be described entirely by anchoring the light curve with early luminosity {L ν,pk ; t pk ; t tran ; t age ; γ(α)}, anchoring with late luminosity, {L ν,0 ; t pk ; t tran ; t age ; γ(α)}, or any other combination such as luminosity at early and late times with the age being inferred, {L ν,pk ; t pk ; t tran ; L ν,0 ; γ(α)}.
In summary, the model for the SNR light curve we will use throughout the rest of this paper is: This is shown in Fig. 5. For the model to make sense, the light curve should reach the peak luminosity before its transition to the adiabatic phase, i.e. t pk > t tran .

Supernovae rate
In our analysis, we will not only consider the SNRs that have already been observed and collected in the Green catalog but also make projections for possible SNRs that could be observed in the future. Thus, in this section, we will estimate the total number of SNRs in our galaxy. It is estimated that the supernovae formation rate is about 0.02 to 0.03 per year [64]. Considering a time ∼ 10 5 years, the time scale for a SNR going through the first three phases before merging with the ISM, there should be 2000 ∼ 3000 SNRs in the Milky Way. Yet there are only 294 SNRs that have been successfully identified so far [55]. Explanations of the deficit in observation and completeness of the SNR catalogs are discussed in [52,54,63,65].
We will further estimate the number of SNRs in our neighborhood. In Ref. [54], the distribution of SNRs inside the Milky Way is estimated using an empirical two-dimensional density distribution: where Σ(R) is the surface density of SNRs at radius R from the galactic center and R = 8.5 kpc is the distance from the Sun to the galactic center. The best fit in Ref. [54] is a = 1.09, b = 3.87. In Ref. [66], however, the best fit is a = 2.00 ± 0.67, b = 3.53 ± 0.77 with the same fitting function. In either case, for a total of 2000 SNRs with ages within 10 5 years, we expect O(10) SNRs within 1 kpc of the Sun. As we will show that old SNRs nearby could potentially lead to strong signals, this estimate implies that we could look for echoes resulting from old SNRs as close as within 1 kpc.

Extra-old SNRs
In our analysis, we will consider SNRs as old as or even older than a few ×10 4 years. We refer to these SNRs as extra-old SNRs. An extra-old SNR's evolution may no longer be in the free expansion or adiabatic phase. If an old SNR exits the adiabatic phase into the third phase, the light curve model in Eq. (3.7) could no longer apply to the subsequent evolution of its luminosity. However, these SNRs can still contribute significantly to the echo signals. This is because the signal is an integral of SNR's flux density when it was much younger. There are two potential benefits from looking at extra-old SNRs. First, if an old SNR is further away from the galactic center than the Earth, their wavefronts pass the Earth and propagate toward a dark matter dense region, i.e. close to the galactic center, which could enhance the integral in Eq. (2.5). 4 Second, the extra-old SNRs have higher statistics than those still in the first two phases. Assuming a uniform formation rate of supernovae in a fixed volume around the Earth, we can expect 10 times more extra-old SNRs with age 10 5 years than younger SNRs < 10 4 years.
Observationally, since an extra-old SNR is long past its peak time, there could be no historical record of its explosion. Direct observation of the radio emission from such a source today might be challenging. We may not have a complete catalog of all the old SNRs today. Therefore, we should look for echoes even if there is no observed SNR in the opposite direction. A positive identification of the echo signal from extra-old SNRs can in turn motivate searches for these missing SNRs.

Detecting the axion echoes
In this section, we will discuss the procedure to estimate the signal-to-noise ratio for the echo signal at a telescope on Earth. Detailed computations, which we will show later, demonstrate that the best sensitivity could be achieved in the radio band, due to a combination of large Bose enhancement to boost the signal and relatively low background. Thus we will focus on the radio telescopes. In the near future, one of the most powerful radio telescopes will be Square Kilometer Array (SKA). We choose SKA as a benchmark experiment to evaluate the prospect of detecting an echo signal triggered by SNRs.
SKA, with two sites in Australia and Africa, could provide over a million square meters of collecting area through thousands of connected telescopes when both phases (SKA1 and SKA2) are completed [67]. Since SKA1 is being designed now, we will only focus on it for the rest of our work. SKA1 is further divided into two frequency ranges: SKA1-low and SKA1-mid. Currently, the setup of SKA1-low is designed to have [67,68]: • frequency range: 50 MHz -350 MHz; • number of stations: 512; • effective diameter of each station: 38 m; • longest baseline (separation between two stations): 80 km, while SKA1-mid has [68] • frequency range: 350 MHz -15.4 GHz; • 133 dishes with a diameter of 15 m (SKA dishes), 64 dishes with a diameter of 13.5 m (MeerKAT dishes); • longest baseline: 150 km.
For SKA1-low, each station will include 256 antennas, and work as a single telescope. For SKA1-mid, both types of dishes are offset Gregorian reflectors with four frequency bands. SKA could operate in both the single dish and the interferometer modes, each with its own advantages and disadvantages. We will consider both modes in our estimates and review their basic formulas below.

Single dish mode
For the single dish mode, the angle corresponding to the Rayleigh criterion, or equivalently the angular resolution of the telescope, is given by [56] θ res = 1.22 where λ is the wavelength and d is the diameter of the dish. The angular size of the signal, θ echo is: where θ s is the angular size of the SNR source, θ ab is the aberration angle in Eq. (2.9), and 2δ, given by Eq. (2.11), is the extended angular size to overcome the loss due to the peculiar motion of dark matter.
In the case when θ echo < θ res (the inner structure of the signal will not be resolved in this case), the noise power of a single telescope is given by: Here the optimal search bandwidth ∆ν is given in Eq. (2.6); and t obs is the observation time, for which we take 100 hours as a benchmark for both SKA1-low and SKA1-mid. The factor of 2 in the numerator accounts for the polarization degrees of freedom of the photon. The denominator accounts for the noise reduction due to independent measurements: ∆ν t obs measurements in time, and those same 2 photon polarizations [39,56]. Note that this expression is independent of θ echo . The total system noise temperature, T sys , has multiple contributions: where T cmb is the CMB temperature; T gal is the synchrotron background from the Milky Way or extragalactic sources, evaluated at the position of the echo in the sky; T atm is the atmospheric temperature; T rcv is the receiver temperature; and T spl is the spill-over temperature. 5 These temperatures could all depend on the frequencies being probed. We take T cmb = 2.725 K and ignore its small frequency dependence. T gal is computed from the 408-MHz Haslam all-sky map [69] by converting from 408 MHz to the frequency being observed. 6 This is done by using the power law fit, T gal (ν) = T 408 (ν/0.408) p , with exponent p = −2.55 [69] and frequency ν in unit of GHz. Note that T gal is the dominant contribution to T sys for SKA1-low since it probes a lower frequency range. We infer T atm from subtracting T cmb andT gal (T gal averaged over the sky probed by the telescope) from the total sky temperature (a sum of T cmb ,T gal and T atm ) provided in Fig. 19 of Ref. [68] for SKA1-low and in Fig. 6 of Ref. [70] for SKA1-mid. For SKA1-low, T rcv = 40 K (Fig. 9 of [71]). For SKA1-mid, the function of T rcv for SKA dishes and MeerKAT dishes is provided in [70]. For the overlapping frequency range of both types of dishes, we take a weighted mean of their T rcv using their number of dishes as weights; we use instead T rcv of the SKA dish for the rest of the frequency range. T spl is set to be zero for SKA1-low and 3 K for SKA1-mid.
In the case when θ echo > θ res , the noise power of a single telescope is modified to be: 7 whereT sys is the average system temperature over the solid angle of the signal ∝ θ 2 echo . The enhanced signal collecting area expanded by the angle θ echo instead of θ res has two competing effects: it leads to an increase in both the noise power and the number of independent measurements by a factor of (θ echo /θ res ) 2 . This results in the additional factor of (θ echo /θ res ) in the noise power in Eq. (4.5), compared to Eq. (4.3). The total noise for N telescopes operating in single dish (SD) mode is On the other hand, the signal power is given by: where S e is the total axion echo flux (i.e. irradiance), ∆ν is the frequency bandwidth, f ∆ ≈ 0.84 the fraction of the flux that falls within that band (see Eq. (2.6), and Ref. [36]), S νa,e the average flux density (i.e. spectral irradiance) over the bandwidth (see Eq. (2.5)), A tot the total geometric collecting area of the telescopes, and η is the ν-dependent detector efficiency. Note that the number of photon polarizations is already included in S e and S νa,e , because of the phase space integral over the momentum measure DP (see Sec. 2). We extract η from the sensitivities of SKA1 provided in Ref. [70]. 8 We show the efficiency in the left panel of Fig. 6. Finally, the signal-to-noise ratio in the single dish mode is simply given by: S νa,e S ν, noi; SD , (4.8) with S ν, noi; SD as the noise flux density for an array in single dish mode [72]: with A unit ≡ A tot /N as the average area of each individual telescope in the array.

Interferometer mode
A single pair of telescopes (a two-element interferometer) has an angular resolution [56] where B is the length of the interferometer baseline (the separation between the two telescopes). Since the baseline is usually much longer than the diameter size of a single dish, the interferometer usually has a much better angular resolution than the single dish. In addition, for N telescopes, there are N (N − 1)/2 pairs of two-element interferometers, which could improve the signal-to-noise ratio by an additional factor of √ N , compared to the single dish mode. On the other hand, for an image with angular size θ echo > θ b , the visibility function R ∝ [sin(πθ echo /θ b )] /(πθ echo /θ b ) falls off fast to almost zero (Sec. 9.2.3.2 in Ref. [56]). Thus not all telescope pairs could contribute to the detection. In particular, those pairs of telescopes far separated from each other will not be sensitive to objects with large angular sizes. Therefore, for an extended echo signal, the single dish mode could still potentially perform better than the interferometer mode.
From the SKA1 design summary [67], we extract the configuration of the telescopes. We then count the baseline length of each pair to obtain a distribution of the baselines. We parametrize it with g(B), the number of baselines shorter than B. This function is presented in the right panel of Fig. 6. As discussed above, a pair of telescopes contribute to the detection only if B < λ/θ echo .
The RMS noise level for the interferometry (IN) mode is similar to Eq. (4.3), after correcting it by the increased noise from the multiple telescopes and dividing by the statistics of multiple independent measurements P noi; IN = P noi; unit N N pairs , (4.11) where N is the number of telescopes active in the interferometry mode and making up N pairs , the number of pairs of telescopes that contribute to the detection: We denote the minimum baseline length in SKA1 B min (if B min > λ/θ echo , the interferom-eter mode will fail to work). The signal-to-noise ratio is: where S ν, noi; IN is the noise flux density [72]: (4.14) For a pointlike source, this implies a sensitivity equivalent to an antenna of aperture N pairs A ∼ N (N − 1)A ∼ N A when the number of telescopes N is large.

Results
Having already developed the formalism for the computation of the axion echo signal strength, as well as having established our treatment of the SKA telescope's properties and detection abilities, we devote this section to their application to a variety of SNRs.
We will endeavor to demonstrate the effectiveness of this technique in axion dark matter detection efforts. To do this, we calculate the reach of the axion-photon coupling for various axion masses, focusing on SKA1 as an example radio telescope. We look first at the discovery potential of SNRs that are already observed, and afterward, we consider the reach projections of potential future SNRs.

Constraints from the Green Catalog
We analyzed the SNRs in a catalog 9 of 294 supernova remnants [53,55] (henceforth the "Green catalog" or GC). We further supplemented the GC catalog with another SNR catalog 10 [73], which includes age estimates for many SNRs also found in the GC. We then selected those SNRs from the GC with a known distance, flux density today, spectral index, and age, after which we are left with a total of 60 SNRs that can be used for our computation. For an SNR whose age is uncertain and known only within a range, we took the geometric mean of the upper and lower bounds as a benchmark value. Subsequently, we calculated the SKA1 telescope signal-to-noise ratio of the SNRs' corresponding axion echoes in two separate cases: 1. "Adiabatic-only": We considered contributions to the echo coming only from the portion of the line of sight corresponding to the adiabatic phase of the SNR's history, thereby excluding the free expansion phase. In this conservative case, the only free parameter is then the duration of the adiabatic phase or, alternatively, the age t tran at which the SNR transitioned from the free to the adiabatic expansion, typically ∼ O(100 years) [56].
2. "Free + Adiabatic": We included contributions to the echo coming from the free expansion phase. The signal is then enhanced with respect to the previous case, at the cost of one more unknown parameter. Indeed, the echo signal depends on a choice of two out of three parameters {t tran , t pk , L ν,pk }. The third one can be calculated in terms of the other two by using Eq. (3.7), for a SNR with known age and luminosity today.  1GHz,s is the flux density of the source today, for which there is roughly ∼ 10% uncertainty [53,55,74]. In our computations we take the geometric mean of the age range given in [73], namely ≈ 5.5 × 10 4 years. For the distance D and the spectral index α, we use the values found in Ref. [53,55], while for the adiabatic expansion index γ, we assume Eq. (3.5), the Sedov-Taylor formula. In parenthesis we quote the uncertainty ranges of these parameters, according to various studies [75][76][77][78]. For more details about the impact of these uncertainties on sensitivities, see Appendix A.
In Fig. 7, we show the expected SKA1 reach for the most promising SNR in the GC, G39.7-2.0 (also known as W50), whose properties are listed in Table 2. In this case, we anchor the light curve with its flux density today and extrapolate backward. We assume that the SNR is still in its adiabatic phase, therefore the light curve evolves according to Eq. (3.6). Since we are interested in demonstrating the sensitivity of the axion echo technique for axion dark matter detection by SKA, we show the reach setting the signalto-noise ratio s/n = 1 for illustrative purposes. In our figure, we have combined SKA1-low and SKA1-mid in a single curve, and have used either the single-dish or the interferometer mode, depending on which mode is the most sensitive at the particular frequency under consideration. We have assumed a typical value of t tran ≈ 200 years. The conservative case of adiabatic-only contribution to the axion echo is shown in blue, while the reach whose computation also included contributions from the free expansion is shown in red. For the latter, we have further assumed a peak time of t pk = 3 × 10 3 days, for which Eq. (3.7) yields a peak luminosity of L 1GHz,pk = 3 × 10 30 erg s −1 Hz −1 .
From this figure we can see that SKA1 can be sensitive to echoes produced by axion dark matter decays stimulated by SNR G39.7-2.0 for a range of axion masses between 4 × 10 −7 and 2 × 10 −6 eV, down to axion couplings of g aγγ = 2 × 10 −11 GeV −1 , a factor of ∼ 3 below the current CAST limit. While the precise value of g aγγ that a hypothetical . The sensitivity reach for axion dark matter coupling g aγγ of SKA1, for the echo produced by SNR G39.7-2.0. We take a signal-to-noise ratio of s/n = 1. SKA1-low and SKA1-mid, in both their single dish and interferometer modalities, are combined into a single curve. We assume a typical value of t tran = 200 years. The reach for the conservative "adiabatic-only" case is shown in blue. The reach for the "free+adiabatic" case is plotted in red. For the latter case, we further assume t pk = 3, 000 days, which yields L 1GHz,pk = 3 × 10 30 erg s −1 Hz −1 . The grey regions are existing bounds from previous literature: CAST [14,15], ADMX [16][17][18][19][20][21][22], RBF+UF [23,24], CAPP [25][26][27], HAYSTACK [28], QUAX [29,30], and neutron stars [31][32][33]. These bounds are taken from github.com/cajohare/AxionLimits. discovery could give would depend on the SNR properties (in particular, during the early free expansion phase), the computation demonstrates the potential of the echo search strategy to probe axion dark matter. In Fig. 8, we explore the effect that a different choice of values for t tran and t pk has on the SKA sensitivity to axion echoes, for an axion mass of m a = 10 −6 eV or, equivalently, a photon frequency of ν a ≈ 0.1 GHz. The green contours represent the g aγγ sensitivity reach for s/n = 1, and the golden star represents the choice of t tran and t pk used in Fig. 7. In blue we show the corresponding contours for the peak luminosity at 1 GHz, L 1GHz,pk . Their corresponding deviation from the typical SNRs values described in Table 1 are shown in purple, at +2σ and +3σ levels: it can be seen that G39.7-2.0 is an atypically bright SNR (around the peak time t pk shortly after the explosion), which is why it is such a good candidate for axion echo discovery. 11 Note that G39.7-2.0 is not the brightest SNR today.  Figure 8. The sensitivity reach for axion dark matter coupling g aγγ of SKA1 at m a = 10 −6 eV (ν a ≈ 0.1 GHz), as a function of light curve parameters t tran and t pk , for the echo produced by SNR G39.7-2.0. We take a signal-to-noise ratio of s/n = 1. The blue curves represent the values of L 1GHz,pk (in units of erg s −1 Hz −1 ) that satisfy the light curve equation Eq. (3.7). How typical these values are is represented by the purple curves, which denote the number of standard deviations from the mean of L 1GHz,pk [59]. The horizontal red lines in turn represent the mean (solid), ±1σ (dashed), and ±2σ (dotted) values for t pk (idem). The golden star represents the light curve parameters assumed in Fig. 7.
In fact, its current flux density is only about 1/30 of that of another SNR, Cassiopeia A, the brightest extrasolar radio source. It is the luminosity in the early times of G39.7-2.0 that makes it a more promising source than Cassiopeia A. The horizontal red lines in turn represent the mean (solid), ±1σ (dashed), and ±2σ (dotted) values for t pk (see Table 1). It is evident from this figure that satisfying Eq. (3.7) for typical values of t tran and t pk yields large L 1GHz,pk and even stronger sensitivities to g aγγ .

Projection of unobserved SNR
In this section, we want to first explore echo signals triggered by extra-old SNRs with ages above 10 4 years, which are briefly discussed in Sec. 3.3. They could be long past their peak times and have not been observed so far. But they could still potentially lead to observable axion echo signals, as we will show. In addition, we will also consider an opposite scenario: young supernova explosion that could happen close to us in the near future. In this case, if the peak flux density is very high (e.g., comparable or higher than that of a solar burst), it could also lead to an observable echo signal. multiply our L 1GHz,pk by the conversion factor 8 GHz 1 GHz −α in order to compute the corresponding L 8GHz,pk and obtain the standard deviations (see discussion in Sec. 3).

Benchmark
A (adiabatic only) B (free+adiabatic) C (old) D (new)  Table 3. Astrophysical properties of three hypothetical SNRs/supernovae in the four benchmarks. Note that the SNR in benchmarks A and B are the same. For benchmark A, we only consider its adiabatic phase while for benchmark B, we consider both the free expansion and adiabatic phase. The quantities with an asterisk symbol means that they are derived based on other inputs. Note that θ s is the size of the SNR based on an estimate of its radius. Luminosity is given in the cgs unit of erg · s −1 · Hz −1 .
We consider the following four benchmark SNRs/supernovae shown in Table 3. In benchmark A and B, we take a SNR with similar properties to G6.4-0.1, an observed SNR in GC which could lead to even stronger echo signals, compared to G39.7-2.0 discussed in the previous section. But G6.4-0.1's signal is in the northern hemisphere and could not be observed by SKA1. The hypothetical SNR in the first two benchmarks are, however, taken to be located in the northern hemisphere so that the echo in the southern hemisphere can be observed by SKA1. We fix its current flux density to be the same as that of G6.4-0.1 and extrapolate backward, assuming the light curve model in Eq. (3.7). In benchmark A, we only consider the signal from the adiabatic phase. In benchmark B, we include the free expansion phase, assuming t tran = 30 t pk . In benchmarks C and D, instead of fixing the flux density today and extrapolating backward, we model the first two phases of the expansion by anchoring the early stages of the light curves. We choose L 1GHz, pk and t pk to be within the total 2σ range shown in Ref. [59] and quoted in Table 1. Note that in [59], the analysis is implemented in the frequency range 4 − 10 GHz except for SN1987A. As described in Sec. 3 and the previous subsection, we convert this to the fiducial frequency 1 GHz used here using the spectral index. The only differences in the inputs of benchmarks C and D are their ages and galactic coordinates. In benchmark C, even though the echo comes from a region around the galactic center, it is still more than 5 • away from it. Therefore, we neglect the optical depth as discussed in Sec. 3.3. Note that the SNRs in benchmark A, B, and C are extra-old ones with ages above 10 4 years while the SNR in benchmark D is a baby one, which is still in the very early stage of the free expansion phase shortly after the explosion.  Figure 9. We show SKA1's sensitivities to axion DM for four benchmarks in Table 3. We set the signal-to-noise ratio, s/n= 1. The blue curves only include the second phase (adiabatic expansion) of SNR, while the red curves also include the first phase (free expansion) as well. The solid (dashed) curves correspond to the interferometry mode (single-dish mode) of SKA1, respectively. In benchmarks A and B, the luminosity today is taken to be the same as SNR G6.4-0.1 in GC [53,55]; in benchmarks C and D, we choose L pk and t pk to be within the total 2σ range based on recent supernova explosions [59]. See the main text for more details. The grey regions are the existing bounds: CAST [14,15], ADMX [16][17][18][19][20][21][22], RBF+UF [23,24], CAPP [25][26][27], HAYSTACK [28], QUAX [29,30], and neutron stars [31][32][33]. These bounds are taken from github.com/cajohare/AxionLimits.
For benchmarks C and D, we estimate the SNR/supernova radius as a function of their ages. During the free expansion phase, the speed of homologous expansion is estimated to be v hom ≈ 0.04c [60]. The radius is given by The sensitivities of these four benchmarks are presented in Fig. 9. For all these benchmarks, we see that for the low frequency range (corresponding to m a around 10 −6 eV), searching for echo signals at SKA1 could provide a better reach compared to CAST. A couple of comments are in order: • At low frequencies, the interferometry mode mostly dominates because of the lower array noise suppressed by N pairs . At high frequencies, the resolution is enhanced, which leads to a suppressed visibility function as we discussed in Sec. 4. As a result, the number of baselines that can observe the echo signal drops as the frequency increases. On the other hand, the number of dishes/stations that can observe the signal is constant in the entire frequency range of SKA1-low (SKA1-mid). One could verify that when the number of active baselines, N pairs , drops to be equal to that of dishes/stations, the two modes have equal sensitivities. At even higher frequencies, the single dish mode dominates. This explains the crossover of the two curves in all four plots in Fig. 9. The single dish mode is almost always better at frequencies higher than 1 GHz except for baby supernovae as in benchmark D. For a very young supernova close by, which could be in its free expansion phase, it could have a small source size as well as small aberration angle and δ in Eq. (2.11) due to dark matter peculiar velocity. Thus the interferometry mode may not suffer from the degrading of the visibility function, which is present for more extended sources.
• In benchmarks A and B, the peak luminosity is just above the 2σ range in [59], just as SNR G39.7-2.0 discussed in the previous section. This may not be an issue since Ref. [59] is a statistical study based on very recent supernovae explosions, while we consider a very old SNR in benchmarks A and B. As a result of selection bias, the samples of [59] may not be representative of the supernovae leading to old SNRs. On the other hand, in benchmark C, we choose the peak luminosity to be within the 2σ range in Ref. [59], and show that it could still lead to a strong echo signal, even though its luminosity today is at the tail of all well-measured SNRs [55].
Lastly, we show the dependence of the reach for g aγγ on SNRs' intrinsic properties. We vary two parameters at a time and fix the other parameters to be the same as in benchmark C. We compute SKA1's reach of g aγγ at m a = 10 −6 eV. We show the g aγγ contours in each two-dimensional slice of the parameter space. The results are collected in Fig. 10. From the figure, we make the following observations: • The signal is most sensitive to the peak luminosity, L pk , and the peak time, t pk .
• An SNR nearby leads to a better reach, if the luminosity is fixed. This is expected as it is the flux density of the source (inversely proportional to the distance squared) that determines the signal flux density as shown in Eq. (2.5).
• When the echo comes from the region around the galactic center, old SNRs lead to stronger signals. This is due to the wavefront of old SNR propagating into the dark matter-dense region of the galaxy. One should however ignore the reach at echo = b echo = 0, since the optical depth could be large from the galactic center.
• The reach has little dependence on the specific transition time, t tran , from the freeexpansion phase to the adiabatic phase, when both phases are considered fixing L pk and t pk .

Conclusions
In this paper we discuss the axion echo radio signal, namely the photons from axion stimulated decays that travel back to the radio source. We stress that there are two novel features of studying the echo signals compared to the forward-going radio signals also from axion stimulated decays: the low background, and its dependence on the history of the radio sources. The latter enables us to look for time-varying radio sources. We identify SNRs as good candidates and demonstrate that they could probe interesting parameter regions of axion dark matter. Even for the SNRs that are not the brightest radio sources we observe today, they can generate axion echoes stronger than the bright constant radio sources such as the radio galaxy Cygnus A.
In discussing the candidate SNRs, we model their light curves with inputs from SNR catalogs as well as existing statistical studies based on newly exploding supernovae. We show the parameter region of (m a , g aγγ ) probed by SKA1 with known SNRs and with hypothetical but otherwise statistically normal ones that could be discovered in the future. We observe that SKA1 can reach g aγγ as low as ∼ 10 −11 GeV −1 around m a ∼ 10 −6 eV. We also discuss the theoretical uncertainties in the modeling of the SNR light curves, the dependence of the SKA1 reach on SNR parameters, as well as the interplay with the two operation modes at SKA1.
There are several important directions to expand our work: • As we have emphasized, our work intends to provide a proof of concept estimate to demonstrate the power of searching for axion echo signals triggered by SNRs. More precise modeling of the SNR light curves could improve the estimate and put it on more solid ground.
• There is the very intriguing possibility that even for an SNR that has not been observed yet (i.e. , it fails to pass the selection criteria of SNR searches), it could still lead to an observable echo signal. This could result in the interesting scenario of observing a radio line signal with no bright source in either direction along the line of sight as if the source is a "ghost". In other words, searching for the echo signal could not only benefit the axion dark matter search but also the SNR search.
• We only focus on SKA1 as an example of a powerful radio telescope. One could extend the work to consider other radio telescopes. Note that with SKA1 to be constructed,  Figure 10. We show SKA1's reach of g aγγ and its dependence on the SNR properties. We set the signal-to-noise ratio, s/n= 1, and compute the sensitivity with fixed axion mass at m a = 10 −6 eV (ν ≈ 0.1 GHz). The contours are computed with both interferometry mode and single dish mode, with the former dominating at this frequency most of the time. The contours are labeled by g aγγ /(10 −10 GeV −1 ). The dashed circle in the top left corner block corresponds to the 2σ contour of (t pk , L pk ) from [59]. Note that we vary two parameters in each panel, fixing all the other parameters to be the same as in benchmark C in Table 3. We consider the first two phases of the SNR to compute the signal strength. The galactic coordinates are for the echo signal instead of the SNR. The small black dot at the center in the − b slice marks the galactic center, where the optical depth could modify the SKA1 reach.
studies of its potential in addressing fundamental questions and unveiling new physics will be useful and well timed. With SKA2's design currently under discussion, we leave it for future work to study its capability to detect axion dark matter in the mass range of 10 −7 eV to 10 −4 eV.
• The 408-MHz Haslam all-sky map still contains faint point sources that contributes to the background noise. With better measurements, the point sources could be removed in the future, which could potentially further suppress the noise and deepen the reach of SKA in g aγγ .
Note added -During the late stages of the completion of this work, we became aware of other work [79] by Leung, Masui, Nambrath, Schutz, and Sun on the same subject. We note that the approaches and analyses of both works, while sharing common components, differ and complement each other in several aspects. phase of the light curve history in the computation of the echo signal (the conservative case we dubbed "adiabatic-only" in the main body of this paper). As described in Sec. 5, the sensitivity could be improved with the inclusion of the echo flux from the free expansion phase of the light curve.
• Flux density (S ν we used in Sec. 5.1 is taken from the Green catalog [53,55], and is 85 Jy at ν = 1 GHz. Reference [74] quantifies the uncertainty in the flux density of this SNR and it is roughly ∼ 10% for measurements at wavelengths of 6, 11, and 21 cm (frequencies of 5, 2.7, and 1.4 GHz respectively). Since the echo signal scales linearly with the flux density (Eq. (2.5)) and, as we shall see below, the impact from uncertainties in other quantities far outweighs that from a 10% uncertainty of S (0) ν , we will not consider this effect any further.
• Distance (D): The Green catalog quotes D = 4.9 kpc [53,55]. Recent studies place the exact value between 4.5 kpc and 5.5 kpc [75][76][77]. Since, however, it is the flux density S ν that is relevant for the computation of the axion echo (Eq. (2.5)) and of the signal strength (Eqs. (4.8) and (4.13)), an uncertainty in D merely translates into an uncertainty in the spectral luminosity L ν . We therefore do not consider the effect of this uncertainty any further.
• Age (t age ): Ref. [73] and the public catalog SNRcat 12 cite an age between 3 × 10 4 and 10 5 years. This age roughly coincides with the typical end of the adiabatic expansion and the beginning of the snow plough phase (Ref. [62] cites a benchmark of ∼ 5 × 10 4 years). During the snow plough phase, the SNR rapidly loses energy, which translates into a weaker stimulation of axion decays and consequently a weaker signal. Nevertheless, since the estimated age is right around the transition between both phases, and in light of the larger uncertainty in the spectral index (see next item), we refrain from modeling the snow plough phase (which will complicate the SNR light curve model by introducing more parameters). In the top left panel of Fig. 11, we show in red the impact of the uncertainty in the age of SNR G39.7-2.0, considering echoes from the adiabatic phase only. The dashed black line corresponds to our fiducial age of t age ≈ 5.5 × 10 4 years, the geometric mean of the upper and lower ends of the age.
• Spectral index (α): The Green catalog quotes α = 0.7 [53,55], but its precise value could lie between 0.5 and 0.8 [78]. Since α is an exponent, the uncertainty in its value has the largest impact on the sensitivity reaches. It not only affects the dependence of the reach on the axion mass (since S ν ∝ ν −α as shown in Eq. (3.1), and ν = m a /4π), but it also directly impacts the adiabatic expansion index γ, which according to the Sedov-Taylor formula (Eq. (3.5)) is given by γ = 4 5 (2α + 1). The uncertainty in α then translates into an uncertainty in γ, which lies between 1.6 and 2.08. Since the axion echo is an integral over the SNR history (Eq. (2.5)), the value of the adiabatic index is crucial in determining whether the echo signal is observable.  Figure 11. The uncertainties in the sensitivity reach for g aγγ of SKA1, for the echo produced by SNR G39.7-2.0, as a result of the uncertainties of the SNR parameters. We take a signal-tonoise ratio of s/n = 1. SKA1-low and SKA1-mid, in both their single dish and interferometer modalities, are combined into a single curve. We take the conservative "adiabatic-only" case (i.e. echoes produced by the adiabatic phase of the light curve only). In each panel, we vary a single parameter while keeping the other parameters fixed at the fiducial values used in the main text.
Top left: in red, varying t age between 3×10 4 and 10 5 years. Top right: in yellow, varying α between 0.5 and 0.8. Bottom left: in green, varying t tran between 100 and 300 years. Bottom right: in blue, the combined effect of the previous three cases; it is clearly dominated by the uncertainty in α.
The effect that the uncertainty in α has on the sensitivity reach is shown as the yellow band in the top right panel of Fig. 11; the dashed black line takes the Green catalog value of α = 0.7 and corresponding γ = 1.92.
• Transition time (t tran ): The time at which the SNR transitioned from the free expanding phase to the adiabatic phase is typically of the order of O(few × 100) years (Ref. [62] quotes a benchmark of ∼ 190 years). Without a dedicated study of the astrophysical properties of SNR G39.7-2.0 (possibly including astronomical observations and numerical modeling of the supernova remnant evolution), it is hard to know with precision at what time this transition roughly took place. In the simplest analytical estimates, t tran is a function of the interstellar medium number density, the energy of the supernova explosion, and the mass of the ejecta [60,62]. We can however quantify the impact of this uncertainty on the sensitivity reach by simply varying t tran as a free parameter. In the bottom left panel of Fig. 11, we show in green the resulting reach band, allowing t tran to vary between 100 and 300 years; the dashed black line assumes t tran = 200 years, as was used in Fig. 7.
Finally, in the bottom right panel of Fig. 11 we show in blue the sensitivity band resulting from the combination of all these variations. This band serves as a quantitative proxy of the uncertainty in the sensitivity reach for SNR G39.7-2.0. Clearly, the uncertainty in the spectral index α (and consequently on the adiabatic index γ) has the dominant effect. A more precise determination of the spectral index of G39.7-2.0 and of its expansion history could improve estimating the sensitivity of SKA1 to axion echoes.