Axion dark matter-induced echo of supernova remnants

Axions are a theoretically promising dark matter (DM) candidate. In the presence of radiation from bright astrophysical sources at radio frequencies, nonrelativistic DM axions can undergo stimulated decay to two nearly back-to-back photons, meaning that bright sources of radio waves will have a counterimage ("gegenschein") in nearly the exact opposite spatial direction. The counterimage will be spectrally distinct from backgrounds, taking the form of a narrow radio line centered at $\nu = m_a/4\pi$ with a width determined by Doppler broadening in the DM halo, $\Delta \nu/\nu \sim 10^{-3}$. In this work, we show that the axion decay-induced echoes of supernova remnants may be bright enough to be detectable. Their non-detection may be able to set the strongest limits to date on axion DM in the $\sim 1-10 \, \mu$eV mass range where there are gaps in coverage from existing experiments.


I. INTRODUCTION
To date, dark matter (DM) has been detected only via its gravitational interactions, with additional nongravitational signatures possible in specific models of DM. In one such class of models, the DM is comprised of axions, which were first proposed as a solution to the strong-CP problem in QCD [1][2][3][4][5][6]. They also generically arise in string theory [7][8][9][10][11] and can appear as components of models which address other problems in the Standard Model (SM), such as the matter-antimatter asymmetry [12]. Axions generically couple to electromagnetic fields with a Lagrangian L ⊃ g aγγ aE · B [13] where g aγγ is the coupling strength, E and B are electric and magnetic fields, respectively, and a is the axion field. This axion-photon coupling can be leveraged in a wide array of axion detection strategies.
The most widely probed axion interaction involves the interconversion of axions and photons in the presence of a strong magnetic field. For instance, axion DM direct detection experiments (collectively referred to as "haloscopes") in the m a ∼ 1-100 µeV mass range search for DM axions converting to radio frequency photons inside of a resonant cavity that is tuned to the axion mass, enabling exquisite sensitivity to DM axions with very weak couplings to photons [14][15][16][17][18][19]. The same interaction can also be manifest in the magnetospheres of neutron stars, which are suffused with a plasma that alters the dispersion relation of photons to resonantly enhance the conversion of a massive DM axion to an in-medium radio photon, yielding a spectral line [20][21][22]. Finally, if axions exist (whether as DM or as an auxiliary field in the spectrum of some theory) they will be efficiently produced * yitians@mit.edu † Einstein Fellow; katelin.schutz@mcgill.ca ‡ nambrath@mit.edu § calvinl@mit.edu ¶ kmasui@mit.edu via thermal processes in the Sun with keV-scale energies; solar axion experiments ("helioscopes") like the CERN Axion Solar Telescope (CAST) search for the conversion of a solar axion to an X-ray photon in a strong B-field, placing constraints that are roughly independent of mass for m a 1 keV [23].
Axions can also decay to a pair of photons with lifetime τ = 64π/m 3 a g 2 aγγ , where m a is the axion mass. Given strong existing constraints, it is unlikely that spontaneous DM axion decays can be observed, even in nearby DM-rich dwarf galaxies [24]. On the other hand, the rate of stimulated axion decay can be significantly enhanced in the presence of radiation [25]. In the axion rest frame, the decay is stimulated by an incoming photon with a frequency ω = m a /2 and the two photons from the decay of the axion are emitted at the same frequency exactly back to back, traveling along the same axis as the incoming photon. Therefore, in the axion frame, incident light (for instance, from bright astrophysical radio sources) would appear to be both amplified in the forward direction [26] and reflected in the opposite direction [27], creating an "echo" or counterimage in the form of a narrow spectral line in the opposite direction to the source of photons. The latter effect, referred to as "axion gegenschein" (in analogy to solar gegenschein, which is created by dust reflecting light from the Sun), was recently studied in the context of individual bright radio point sources like Cygnus A inducing a spectral line coming from the antipodal part of the Galactic halo [28].
In this work, we expand upon the idea of gegenschein from axions in the ∼ 1 − 10 µeV mass range (corresponding to ∼GHz frequencies) using supernova remnants (SNRs) as the primary source of photons that stimulate axion decay. SNRs are promising sources because not only are they radio-bright as observed, but they were also substantially brighter in the past. This enhances the gegenschein signal due to the finite speed of light: we can simultaneously observe both i) the decay of relatively nearby axion DM that was stimulated by light that recently passed Earth and ii) the decay of distant axion DM that was stimulated by light that passed Earth a long time ago, when the SNR was orders of magnitude brighter. Since the light-crossing time of the Milky Way (MW) halo is ∼ 10 5 years, we can in principle observe gegenschein from SNRs integrated over these timescales in the form of an image that is stacked along the DM column that is antipodal to the SNR. In this work, we perform forecasts to determine the sensitivity of existing and planned radio telescopes to search for this effect. We consider the Five-hundred-meter Aperture Spherical Telescope (FAST) [29], the Canadian Hydrogen Observatory and Radio-transient Detector (CHORD) [30], and the Square Kilometre Array (SKA) [31], and find that future observations could be sensitive to DM axions with couplings that have yet to be experimentally constrained in terrestrial experiments. Our key findings are summarized in Figure 1. The rest of the paper is organized as follows. In Section II, we provide an overview of the previous work of Ref. [28] on the gegenschein induced by extragalactic radio point sources. We subsequently generalize this framework to account for time-dependent Galactic sources of finite size and distance. In Section III, we model SNR evolution and SNR radio spectra at early times when the luminosity was substantially higher. Given the specifications of various radio telescopes, which we discuss in Section IV, we perform forecasts in Section V and find that even under a variety of modeling assumptions, a null detection of SNR gegenschein can likely constrain new parameter space. Discussion and concluding remarks follow in Section VI. Note that throughout this work, we work in units where c = = 1 unless referring to quantities relevant for observations; in these cases, we provide the relevant units.

II. AXION GEGENSCHEIN FROM ASTROPHYSICAL SOURCES
In this Section, we review the expected gegenschein signal from astrophysical sources of radio waves following Ref. [28], which dealt with the case of a point source of constant brightness and effectively infinite distance. In that limit, one can take the flux from the source to be constant over the entire MW DM halo. These properties serve as a good approximation for the astrophysical sources considered in Ref. [28], namely bright radio galaxies. These sources are at cosmological distances that are much greater than the spatial size of our DM halo, and have variability on timescales that are longer than the light-crossing time of our Galaxy. The brightest radio galaxy in the Northern sky, with gegenschein that can be observed with a Southern hemisphere telescope, is Cygnus A. As we show below, for Galactic SNRs there are several key differences in the gegenschein signal due to the variability of the source luminosity and the fact that they are nearby (within the Galaxy).

A. Gegenschein from a distant point source
In general, in the limit where the axion has a high occupation number, the intensity of the gegenschein can be computed by solving the classical field equations for the axion coupled to electromagnetism. At leading order (ignoring the axion backreaction) the equation relating an ingoing to an outgoing wave is where A 0 and A 1 are the vector potentials of the ingoing and outgoing radiation. A nonrelativistic axion oscillates with a frequency ω ≈ m a and with an amplitude a 0 related to the axion energy density, ρ a = a 2 0 m 2 a /2. Fourier transforming and equating incoming and outgoing radiation, we can compute the flux density S g from axion gegenschein induced by a distant point source, where S ν,0 (ν a ) is the specific flux density of the source at frequency ν a = m a /4π and where the integral is along the DM column in the direction opposite the source assuming a dark matter density ρ(x d ). In this work, we take axion DM to be distributed as a Navarro-Frenk-White profile [32], as a function of galactocentric radius r, with scale radius r s of 16 kpc [33]. We take the local density at the solar position r = 8.22 kpc [34][35][36] to be 0.46 GeV cm −3 [33,37]. If axion DM is bound into com- Geometry of axion gegenschein for a general source at finite distance of finite (time-dependent) size. Finite-distance effects and the finite angular size of the axion-induced counterimage significantly affect our sensitivity estimates.
pact minihalos with M 1 M rather than a smooth distribution (see e.g. [38,39]), our signal is unlikely to be affected as the mass contained within our DM column is significantly greater than the mass of a minihalo and Poisson fluctuations are expected to be negligible in that limit. This is an advantage of this approach over terrestrial searches for local axion DM, since we are sensitive to the DM density over a larger spatial region and are therefore robust to local fluctuations in density.
For a source that is effectively at an infinite distance from Earth, the only geometric factor determining the total gegenschein signal strength is the DM column density. However, DM axions in a halo with nonzero velocity dispersion σ d ∼ 10 −3 are not at rest with respect to the observer, which complicates the form of the signal. The decay spectrum, which takes the form of a narrow line in the axion frame, is broadened by the Doppler effect, yielding a width ∆ν/ν ∼ σ d , corresponding to a ∼MHz width at GHz frequencies. Note that many radio observatories are optimized to be sensitive to spectral lines with similar widths because the 21 cm line is also Doppler broadened by a similar amount in other contexts [40]. An additional complication is that in the observer frame the two photons emitted from the decay are not emitted exactly back to back. Due to transverse motion of the axion, the angle between the photon emitted in the forward direction (which is in the same momentum state as the stimulating photon) and the one emitted in the backward direction is ∆θ ∼ 2σ d . The gegenschein signal will therefore be spatially smeared over an angular region of this characteristic size. In detail, these signatures can be computed using the axion DM phase space distribution. In this work, we adopt the phase space distribution that was empirically determined in Refs. [41,42] using accreted stars as kinematic tracers of substructure in the Solar neighborhood. For this distribution, we find that a typical velocity dispersion of σ d ∼ 116 km/s captures the effect of the full distribution on the gegenschein signal. Since σ d ∼ 10 −3 , the angular extent of the smearing will be a few arcminutes at the very minimum. Note that the total power of the observable gegenschein signal is not reduced by the smearing. Since we are embedded in the MW DM halo which covers the whole sky, signal loss from one direction of the sky can be compensated by a gain from neighboring regions up to negligible differences in DM density over a transverse distance ∼ σ d x d , where x d is the DM distance. Put another way, as depicted in Fig. 2, any "patch" of DM generates a cone of gegenschein, and the Earth sits at the intersection of many such cones for any given astrophysical source.

B. Gegenschein geometry of a general source
In contrast to the example in the previous Subsection, the gegenschein image of a generic source will be complicated by a number of factors. The source may be (1) at a finite distance, (2) spatially extended over a large solid angle, and (3) varying substantially in spatial size and radio brightness on timescales corresponding to the light crossing time of our Galaxy. In this Subsection, we generalize Eq. (2) to account for all of these effects.
In the absence of any DM velocity dispersion, the gegenschein image of a source of finite angular size will be the same angular size as the original, albeit flipped and in the antipodal location. The DM velocity dispersion introduces a nontrivial blurring effect for sources at finite distance, beyond the simple ∆θ ∼ 2σ d effect described in the previous Subsection. In particular, the blurring effect depends on the ratio between the distance to the decaying DM and to the source, with the geometry depicted in Fig. 2. Intuitively, in the limit where the source is closer to us than the decaying DM is, the decaying DM emitting photons at an angle θ d ∼ 2σ d relative to the incoming radiation can deviate substantially from the antipodal axis, with large values of where x ds denotes the distance between the decaying DM and the source, and x s is the distance from the Earth to the source. This effect will cause the image of very nearby radio emitters to be blurred over a very large solid angle. Just like in the case of a distant point source, the total gegenschein signal power is conserved even after accounting for velocity dispersion effects, however this only holds up until a point of saturation where the gegenschein image fills the entire sky. To see this, consider the unphysical case of static DM in an infinitely large halo with the Earth and source separated by a finite distance x s . In this case, all of the initial source radiation will stimulate axion decay at some point and all the corresponding gegenschein will be focused directly back to the source position, as if the source were surrounded by a perfect spherical mirror. In this case, all the gegenschein passes through the spherical shell centered at the source with a radius x s , meaning that the gegenschein has a chance of being observed. However, once velocity dispersions are introduced, the gegenschein is not perfectly focused back to the source. Instead, the gegenschein emission from each patch of axion DM is beamed back towards the source with a conic geometry where the opening angle is ∼ 2σ d and with the axis pointing towards the source. The question then becomes whether most of the emission passes through the spherical shell of radius x s centered at the source (in which case the emission is in principle observable) or whether most of the emission misses that region entirely. In the first case, then we see some finite blurring with very little loss of total power. In the second case, then the power drops off like 1/x 2 d , corresponding to a geometric loss of flux. If we require the image of a point source to be no larger than some θ i0 , for instance the width of a beam or field of view, then this translates to a lower bound on the Earth-source distance given the echoing DM distance x d , where the approximation holds for small θ i0 . Note that θ i0 must be greater than 2σ d and that in this formula we have neglected the intrinsic size of the source and have only focused on finite Earth-source distance effects. For a typical DM column depth of x d ∼ few kpc (corresponding to the distance light can travel over the typical age of a SNR), most Galactic sources should satisfy this bound. For instance, if θ i0 = 30 arcminutes, then any source farther than ∼ 100 pc will satisfy this criterion.
Having established the angular extent of the gegenschein due to finite distance effects, here we construct a general formula for the observed gegenschein intensity I g (n, t ) generated by an extended, time-varying source with specific luminosity per unit volume p ν (x s , t), inside an axion halo with density profile ρ(x d ). The location vectors x s and x d point to the source and the DM from the observer, respectively, as shown in Fig. 2. We assume for simplicity that p ν (x s , t) uniformly illuminates the region surrounding the source and that there is no opacity in any direction. Since the gegenschein image subtends a small angular region on the sky (based on the criterion of Eq. (5)), any anisotropy in source radiation can be absorbed into the local normalization of the source intensity in the gegenschein direction. Consider a source volume d 3 x s . The specific flux density seen by a DM volume where t = t + x ds due to the finite speed of light, and Let dA ⊥ be the cross sectional area of the d 3 x d volume facing the source with d 3 x d = dA ⊥ dx ds , then the total gegenschein power emitted by dark matter at time t is dp g,tot (t )= g 2 The gegenschein intensity observed by a telescope at time t is then where t = t + x d again due to the finite speed of light and f (θ d ) is the gegenschein angular distribution as determined by the DM phase space projected onto the celestial sphere, with θ d being the angle between the stimulating ray and the gegenschein ray, as shown in Fig. 2. Dividing by the solid angle of the volume d 3 x d as seen on Earth, dΩ i , and then integrating, we obtain the gegenschein intensity This is the most general expression we consider for the gegenschein intensity of an extend, time-varying source. Assuming the source is not extended in depth, one can further reexpress the power density of the source p ν as the observed specific intensity I ν as wheren s is the direction of x s . In the limit where the image size θ i 1, we have x ds ≈ x d + x s , and Eq. (4) simplifies to θ i = θ d · x ds /x s , which allow us to re-express f (θ d ) into a widened gaussian function of θ i with the same normalization as f : Eq. (9) then simplifies to where dΩ s is a unit solid angle in then s direction. Note that θ i is the angle between the viewing directionn and the countersource direction −n s . This expression has the following interpretation: For each layer of DM at the distance x d , the gegenschein image is like the source image I ν smeared by an gaussian kernel h(θ i ) whose width depends on x d and the source distance x s . These images are then stacked together weighted by ρ(x d ) to produce the final gegenschein image.

III. SUPERNOVA REMNANTS
Supernova remnants are promising sources for axion gegenschein detection not only because they are radio bright as observed, but they were also much brighter in the past. 1 The travelling time of the gegenschein echo in the Milky Way halo allows us to observe images stimulated by radio waves that passed the Earth up to 10 5 years ago, which means that for many SNRs all of the brightness history will be manifest in the expected gegenschein signal. As we discuss below, the radio luminosity decreases steeply with time in the Sedov-Taylor phase of the SNR's evolution, which typically lasts for ∼ 3 × 10 4 years [44], and the majority of the gegenschein power from SNRs originates from when the source was young. The total integrated gegenschein luminosity along the line of sight is thus expected to be much greater than one would expect if one were to assume the radio luminosity of the source to be constant at its present value.
In order to determine how the luminosity of SNRs evolves with time, the key quantities of interest are the brightness-diameter (Σ − D) relation and the expansion dynamics of SNRs, which we review below. Each SNR has its own properties that determine the evolution, which we account for in our analysis. From the perspective of inducing a strong gegenschein signal, the ideal SNR source is one in the Sedov-Taylor phase of its evolution that is both bright and old, with the latter criterion providing a long lever arm in time to scale back the synchrotron brightness to when the SNR was much brighter. In the remainder of this Section, we first consider general properties of SNRs and then present the detailed considerations that enter the analyses of particularly strong individual candidate sources.

A. Evolution
The evolution of a SNR is roughly separated into four phases: the free expansion phase, the Sedov-Taylor phase (sometimes referred to as the adiabatic phase), the snowplough (or radiative) phase, and the terminal phase. There are transition periods in between phases, however, for a simple model of SNR evolution, we do not explicitly consider them.
The free expansion phase is relatively short-lived, lasting around 100 years. In this period, the evolution dynamics are driven by the supernova ejecta, which travel unencumbered by the surrounding interstellar medium (ISM). This results in a near constant expansion of the shock front, R = v sh t, where v sh is the shock velocity. This phase continues until the ISM mass swept up by the shock front is comparable to the ejecta mass, which slows down the shock front and transitions the evolution into the Sedov-Taylor phase that begins around 1000 years after the supernova explosion.
After the transition, a SNR in a roughly uniform ISM environment can be described by the Sedov-Taylor solution [45], during which a fixed fraction of the total shock energy goes into a conserved bulk kinetic energy E (i.e. we assume no heat loss from the system). We can therefore construct a dimensionless quantity, ξ = R(ρ ISM /Et 2 ) 1/5 , that can be recast as a relationship be-tween the shock radius R and time t, These relations allow us to scale known SNR parameters in the Sedov-Taylor phase back to its beginning, when the SNR is significantly brighter. For some SNRs, the evolution in this phase is complicated by the presence of clouds scattered throughout the ISM that are denser and colder than the bulk of the ISM. Collision and the subsequent evaporation of the clouds will slow down the shock front, however, after taking into account the energy loss, it has been shown that one can construct a new similarity solution whose expansion dynamics are governed by the same relation, R ∝ t 2/5 [46]. Once the age of the remnant becomes comparable to the radiative cooling time of the gas within, the assumption of adiabaticity (no heat loss) becomes a poor one and the SNR transitions into the radiative phase where only momentum is conserved.

B. Radio brightness
For the majority of its lifetime after the initial supernova explosion, SNRs are radio bright due to synchrotron radiation produced by energetic electrons and magnetic fields that are generated as the shock front passes through the surrounding medium. The evolution of the synchrotron radiation flux can be well-modeled to be consistent with observations in the Sedov-Taylor (adiabatic) phase. However in the early stage of evolution where the majority of the total integrated flux originates, it is difficult to determine the precise light curve without incurring modeling uncertainties. In fact, observations of core-collapse supernova light curves suggest that the peak radio luminosity can span 4 orders of magnitude [47]. Rather than adopting purely empirical values for the peak luminosity of a typical young SNR in the context of a large population of young SNRs, which would incur large systematic uncertainty in our estimates, we instead use a combination of present-day measurements of older SNRs in combination with SNR evolution theory to scale back SNR synchrotron emission to its plausible luminosity at early times. In the following discussion, we use a set of simple assumptions to describe a fiducial model, while also presenting alternatives to reflect the effects of modelling uncertainty on the gegenschein constraining power of SNRs. We will discuss the uncertainty associated with model parameters in detail in Sec.V.
We assume SNR radio emission originates from synchrotron radiation by an ensemble of relativistic electrons from the background ISM that have entered the shock front. Their differential energy spectrum is where γ = E e /m e is the electron Lorentz factor, and K e is a normalization factor. Then the specific synchrotron radiation flux at frequency ν observed from distance d away (up to numerical factors depending on the electron power law index p) is where B is the magnetic field strength, 2 and V is the volume in which both the relativistic electrons and the magnetic field are present. Per the conventions of the literature, we use α ≡ (p − 1)/2 to denote the frequency power law index. Historically, equipartition of energy between the ionized particles and fields along with the luminosity has been used to estimate the magnetic field amplitude. However, one must still determine the electron energy density independently from the luminosity to obtain a Σ − D relation. For our fiducial model, we follow the classical analysis of Ref. [48], where it is assumed that electrons enter the shock interior constantly so that V ∝ R 3 , but are decelerated as E e ∝ R −1 due to collision with the expanding magnetic field perturbations [49]. This effect does not change the electron spectral shape, but only modifies the normalization such that which implies V K e (R) ∝ R 1−p . Alternatively, we simply assume that the relativistic electrons evolve adiabatically, and that they take up a fixed small fraction of the total supernova energy throughout the adiabatic phase (this assumption is supported by e.g. simulations in [50]) so that which will result in a slightly shallower projection than our fiducial model. We project our constraints under this alternative electron model separately from the fiducial model. To determine the flux density evolution, we further need to model the evolution of the magnetic field strength in the SNR, which can be quite different during the various evolution phases. Below, we discuss the magnetic field and flux density evolution of SNRs in different phases.

Sedov-Taylor phase
Historically, the magnetic field was assumed to evolve such that the total flux through the shock front is preserved (see e.g. Ref. [48]), yielding the scaling relation B ∝ R −2 . However, more recent radio observations [51], X-ray observations [52][53][54][55], and numerical simulations [56][57][58] suggest that the magnetic field amplitude near the shock front is much greater than one would expect from compression of interstellar magnetic fields, and alternate descriptions of amplification mechanisms are needed. The precise mechanism for magnetic field amplification (MFA) is still an open question. Notably, the diffusive shock acceleration of cosmic rays allows for many mechanisms to produce magnetic fields with a range of amplitude dependence on gas density and shock velocity.
Amplification resulting from resonant streaming instabilities allows the magnetic field to saturate to B 2 ∝ v 2 sh [59], while non-resonant mode amplification saturates to sh [60]. In the Sedov-Taylor phase, these dependences translate to B ∝ R −1.5 and B ∝ R −2.25 respectively. Moreover, it has been shown that non-resonant modes dominate the free expansion and early adiabatic phase, while resonant modes are relevant in later adiabatic phase [61]. For simplicity, we will use a fiducial relation of B ∝ R −2 that falls centrally within the two limits for the entire adiabatic phase.
All together, we thus find that the synchrotron specific flux depends on the SNR radius R as assuming the electron spectrum evolution in Eq. (16). For the alternative model using Eq. (17), we have For a typical SNR with spectral index α = 0.5, the corresponding Σ − D relations are Σ ∼ D 6 and Σ ∼ D 5 , respectively, which are in good agreement with recent observations of galactic SNRs [62]. In our forecast, we will present Eq. (18) as our fiducial model, while showing estimates with the alternative electron evolution model of Eq. (19) as well.

Pre-Sedov-Taylor phases
As mentioned above, due to the steep decrease of radio luminosity with time, the majority of the integrated flux comes from the earliest stages of emission. For SNRs that lack a dense circumstellar medium (CSM) interacting with the initial shock front (e.g. type Ia supernovae), there is an extended brightening phase that continues into the Sedov-Taylor phase [63]. We therefore do not consider such SNRs for their relatively small integrated flux, and focus on core-collapse (CC) supernovae (SN), which have dense CSM environments that immediately interact with the high velocity shock front and amplify magnetic fields on a time scale t MFA < 100 years.
This timescale is often shorter than the onset time of Sedov-Taylor phase t ST ∼ 1000 years. Between t MFA and t ST , the SNR is typically in a transition phase where the power law index of R(t) is greater in magnitude than 2/5, the Sedov-Taylor phase value [64]. This implies the relativistic electron energy E e ∼ V K e will have a steeper t dependence. However, the magnetic field amplitude B(t) is expected to have a weaker t dependence, since B ∼ v sh in the early Sedov-Taylor phase, and v sh has a weaker dependence on t, and tends toward a constant as we trace back to the free expansion phase. Simulations in Ref. [63] suggest that for SN that are dense enough environments (n H > 0.5 cm −3 ), the power law index of S ν (t) stays nearly the same when the SNR is in the pre-Sedov-Taylor phases (up to ∼ 100 years after the SN). For our simple model, it would therefore be reasonable to assume that the S ν (t) power law stays the same in this transition period.
Prior to reaching the peak magnetic field strength at time ∼ t MFA , our knowledge about the SNR's light curve is more uncertain, and may depend on the details of the explosion. According to our previous theoretical assumptions, the relativistic electron energy E e ∼ V K e will decrease steeply with t as R 1−p ∼ t 1−p . Theoretical analyses of post-shock turbulence show that magnetic field energy grows linearly in time before reaching the peak value, i.e. B ∼ B 2 ∼ t. The combination of these effects means that the flux density for a typical SNR depends weakly on time. Observations of young SNR light curves [47] suggest that radio luminosities decrease slightly in the free expansion phase. To avoid introducing new parameters, we conservatively assume that the radio luminosity (and thus flux density) remains constant before t MFA .
To determine the magnetic field amplification timescale t MFA , we refer to the classical magnetohydrodynamics simulations in Ref. [56], which found a log-linear t MFA -v sh scaling. We apply the scalings in that reference to v sh = 5000 km/s, the typical initial shock velocity for a core-collapse SN [65], which corresponds to an amplification timescale of around 100 years after the SN. This timescale is further supported by the more recent simulation in Ref. [63] where the SNR luminosity of a CC SN in a dense CSM was shown to be decreasing since D = 0.5 pc, which corresponds to ∼ 70 years after the SN. Because we can only estimate the timescale for the magnetic field to be amplified up to order unity factors, we will show constraints on axion DM projected for t MFA equal to 30, 100, and 300 years in Fig. 4. As a cross check, we calculated the luminosity distribution at t MFA for catalogued SNR ensemble that we describe in the following Section. The mean luminosity is 10 25.7±1.1 erg s −1 Hz −1 at 6.3 GHz, which is reasonable given that the observed peak luminosity in young core-collapse SN is 10 25.5±1.6 erg s −1 Hz −1 with measurements from 4 to 10 GHz [47]. However, we expect this comparison to be rough, since the ensemble of older SNRs in the catalog described below are mostly Galactic as opposed to the young extragalactic SNRs of Ref. [47]. These two populations are subject to different selection biases regarding how well their various properties (including age and distance) can be measured and therefore the two cannot be compared directly.
In summary, our model fixes the evolution history of the SNR's flux density given the observed flux density S ν,0 (ν), spectral index α, magnetic field amplification time t MFA , and the age t 0 . The expected gegenschein flux scales as where ν a is the axion decay resonant frequency, ρ(x d ) the DM density, and x d parametrizes the DM column depth.

C. Sources
The Galactic SNR candidates are gathered from Green's SNR catalog [66,67], and SNRcat [68,69], combining information such as size, age, distance, radio spectral index and flux density when necessary. Candidates that do not have either upper or lower bound information on either distance or age are left out. They are then run through the constraint projection pipeline for each telescope we consider, and the ones with the strongest projected constraints are individually checked for their properties (some of which are not found directly in the catalogs) such as SN type, initial shock velocity, and spectral index for some SNRs. We then make projections based on the updated information after a detailed review of the literature. In the following, we discuss each of the three candidate sources (Vela SNR, W28, and W50) that provide the strongest projected constraints on axion DM due to the combined factors of their gegenschein power and image sizes. There are many other candidate sources of similar strength in the catalogs, but we focus only on these sources for the purposes of this exploratory study. We also discuss the assumptions and uncertainties involved when using them to project a constraint. We have estimated the effects of proper motion for all of these sources and find that the angular smearing due to the DM velocity dispersion dominates over proper motion effects, which we therefore neglect. We leave discussion about how other uncertainties affect our forecast to Section V. Note that we focus primarily on Galactic SNRs rather than young SNRs with empirically determined light curves, because in spite of their high luminosity young SNRs are generally extragalactic and therefore are too far away to impart a considerable flux on the DM halo.
For all the sources we consider below, the size of the gegenschein image of a SNR is on the order of 10 arcmin., and varies slightly from source to source. The angular size of the image is predominantly determined by the distance and age of the SNR instead of the current observed size, since the total gegenschein flux is dominated by contributions that were stimulated by light coming from the SNR when it was at the very early stages of its evolution, when the source was small. The extent of the gegenschein image thus comes primarily from the blurring effect due to the DM velocity dispersion, amplified by the ratio x ds /x s as shown in Eq. (4). The value of x ds in question is determined by the gegenschein echoing time of the earliest SNR flux, which is in turn determined by the age of the SNR. For a 20,000-year-old SNR that is 2 kpc away, the ratio is x ds /x s ∼ 2.5. Multiplying this with the innate DM blurring size of 2∆θ = 4σ d , we obtain a gegenschein image of 12.5 arcminutes.
Vela SNR (G263.9-03.3) is the closest SNR to earth, with a distance of 287 +19 −17 pc [70]. It is the remnant of a CC supernova with explosion energy around 1.4 × 10 50 erg, fairly small compared to the typical energy of 10 51 erg. The age of the remnant is estimated to be around 1.2 × 10 4 years. Combining estimates of hydrodynamical age [71,72], and characteristic age of the Vela pulsar [73], the uncertainty is on the order of 2 × 10 3 years. Its radio spectrum is well fit by a power law relation, with spectral index α = 0.74 ± 0.04, and flux density 610 Jy at 1 GHz [74,75]. This determination of flux has already excluded Vela X, which is believed to be the pulsar wind nebula of the Vela pulsar and should not be included in the flux of Vela SNR.
While nominally our best source, W50 (SNR G039.7-02.0) as a candidate source for our constraint projection carries more uncertainty than Vela SNR. It is situated near the equator (at declination δ = 5.1°), so its gegenschein image can be observed by telescopes of both the northern and southern hemisphere. Its distance to Earth is around 5 kpc, with recent estimates ranging from 4.5 kpc to 5.5 kpc [76][77][78], while its age estimate varies from 3 × 10 4 years to 10 5 years [68]. SNRs of this age should have exited the Sedov-Taylor phase, and entered the pressure-driven snowplough or snowplough phase, and have begun to lose significant energy radiatively. Nevertheless, assuming a typical transition time (out of the Sedov-Taylor phase) of 2.9 × 10 4 years [44], whether or not one includes the transition to snowplough phase is a very small effect in comparison to the alreadylarge uncertainty in age, and a small effect overall on the constraint projection (as later shown in Fig. 4). The radio spectral index of W50 is also subject to large uncertainties, with measured values ranging from 0.5 to 0.8 in different part of the remnant [79]. We assume the value listed in Green's catalog of 0.7. Effects of these variations on the constraint projection are discussed in Section V. Finally, we assume a log-linear radio spectrum with flux density of 85 Jy at 1 GHz [66,67]. Possible absorption features along the line of sight are observed for a portion  (2) the choice of electron model, as well as upper and lower limits of measured quantities, including (3) the electron spectral index α, (4) the SNR age, and (5) the distance to the SNR. For each thin line, only one of the above parameters are changed from the fiducial model. Note that the late time end points typically correspond to the fiducial age of the SNR and observed luminosity today, but vary for models using the upper and lower limits of the age and distance (the latter of which affects the observed luminosity). of the remnant [79], however, they should not have an effect on a younger SNR when the majority of the radio flux is emitted.
W28 (SNR G006.4-00.1) is estimated to be 1.6-2.2 kpc from the Earth [80] and 3.3-3.6 × 10 4 years in age [81]. It may also have just exited the Sedov-Taylor phase, however we do not consider this change in evolution due to the proximity to the typical transition time and the fact that the age has a very small effect on the final constraint. The radio spectral index above 70 MHz can be fitted from the summary of studies in Ref. [82] to be α = 0.42 ± 0.02 with a flux density of 310 Jy at 1 GHz.
In Fig. 3, we show the assumed luminosity evolution corresponding to the fiducial models we use for our SNR candidates, as well as the variation due to different modeling assumptions and uncertainties in measured quantities described above. We emphasize that our simplified modeling of SNR brightness evolution is likely to be too conservative in predicting the size of the axion gegenschein signal.

IV. OBSERVATORIES
Radio telescopes occupy a wide range of designs optimized for different types of observations. We will broadly group them into the following categories: imaging inter-ferometers, compact mapping interferometers, and single dishes. Here we provide a basic overview of these classes of radio telescopes. We elaborate more on specific interferometers and single-dish telescopes in Subsections IV B and IV A, respectively.
Imaging interferometers include telescopes such as the Very Large Array (VLA), Giant Meter wave Radio Telescope, MeerKAT, the Murchinson Widefield Array (MWA), the Australian Square Kilometer Array Pathfinder (ASKAP), the Low Frequency Array (LO-FAR), and the upcoming SKA. These are composed of an array of modestly sized pointable antennas, with interferometry performed on all pairs of antennas (with each pair acting as a baseline). Every baseline measures a single Fourier mode on of the sky determined by the physical separation of the two antennas in units of the wavelength. Antennas are typically arranged to have relatively long baselines to achieve higher imaging resolution. Critically, extended sources are invisible to long baselines since they lack the relevant Fourier components. For this reason, the sensitivity of imaging interferometers to axion gegenschein from SNRs is highly suppressed by their relatively large angular extent.
Compact mapping interferometers include the Hydrogen Epoch of Reionization Array (HERA), the Canadian Hydrogen Mapping Experiment (CHIME) and its successor CHORD, and Hydrogen Intensity and Real-time Analysis eXperiment (HIRAX). Unlike imaging interferometers, compact mapping interferometers are designed to be sensitive primarily to large structures spanning a large angle on the sky; their ability to study individual sources in detail is limited by their poor angular resolution, but this is not a problem for scientific goals where small-scale characterization of sources is not necessary, e.g. intensity mapping. Many of these telescopes are therefore optimized for wide-field surveys in that they have large fields-of-view, modest point-source sensitivity, and are unable to continuously point to a particular point on the sky for prolonged periods. This is less ideal for studying individual sources with known positions, where the very large field-of-view is not necessary.
Monolithic, single-dish telescopes are optimized for sensitivity to point sources. These include the largest collecting area telescopes such as the Green Bank Telescope, the now decommissioned Arecibo Telescope, and FAST. Single dishes do not resolve out extended sources, and their simplicity and well-controlled systematics make them complementary to modern interferometers. In particular, the unprecedented sensitivity of FAST makes it one of the most promising facilities to observe the axion gegenschein from carefully-selected, bright point sources.
These various radio telescopes are equipped with many digital backends optimized for a wide variety of scientific goals. One key science goal of many such telescopes is the observation of the 21-cm line. Since the Doppler shifting that affects 21-cm observations is similar to the Doppler shift in the axion gegenschein line from the DM velocity dispersion, backends for 21-cm observations are well-suited to the linewidths expected from axion gegenschein. The bandwidth that these telescope backends can process has steadily increased in recent years: modern instrumentation allows O(1) fractional bandwidths to be observed simultaneously, vastly extending the mass reach of an astrophysical axion search.
As discussed above in Sec. III C the SNRs of interest generally create gegenschein images with an angular extent of ∼ 10 arcmin. due to the DM velocity dispersion. For comparison, the scale of the half-power beam width (HPBW) for a telescope with baseline or aperture size D = 300 m, at a frequency ν = 300 MHz is about For single dish telescopes, this means that a larger collecting area cannot increase the received power for this particular frequency, since the beam area decreases like D −2 , although it does increase the received power for lower frequencies. Note that if single-dish telescopes are equipped with multipixel receivers that can collect more total light, the signal-to-noise will increase like the square root of the number of pixels. Meanwhile, for interferometers, baselines longer than 300 m would receive only a fraction of the signal power at ν = 300 MHz because these baselines "resolve out" the image, and hence do not contribute to the measured signal. At frequencies where the beam is saturated, single dish telescopes like FAST win over interferometers of comparable point source effective area such as SKA, due to the fact that the extended image is still visible to a single-dish telescope. Arrays like SKA may still be treated as an incoherent collection of individual pointing telescopes, but in that case the sensitivity of an incoherent collection is reduced.
In addition to specific instrumental considerations, all radio telescopes must contend with contamination which increases the noise floor, which consists primarily of synchrotron radiation from our Galaxy. To estimate this, we use the 408 MHz all-sky continuum survey by Haslam et al. as an estimate for synchrotron emission [83]. The Haslam maps need to be extrapolated from the observation frequency of 408 MHz to a range of synchrotron frequencies as T (ν) ∝ ν −βs .
At low frequency from 45 MHz to 408 MHz, the spectral index β s is measured to have an average of ∼ 2.5 [84], while at higher frequencies between 408 MHz and 23 GHz direct measurements [85] and measurements of the cosmic ray spectral index [86] suggest β s is closer to 3. Since our strongest bounds are slightly below 408 MHz, we will assume β s = 2.5 through our entire frequency range from 70 MHz to 15 GHz. At higher frequencies, the total noise will be dominated by other sources such as the receiver temperature, and our constraint projection should only be slightly conservative.

A. Single dish telescopes
In this Subsection, we describe the sensitivity of singledish telescopes to SNR-induced axion gegenschein. We specifically consider FAST, a single-dish telescope with an illuminated area of A illu = 70700m 2 , corresponding to a diameter of 300 m. The design frequency coverage ranges from 70 MHz to 3 GHz, and up to 8 GHz with future upgrades. The receiver can be moved around in the focal plane, and maintain an aperture efficiency of η A = 0.7 with the zenith angle θ ZA < 26.4°. Beyond this angle, η A decrease linearly to ∼ 0.5 when the zenith angle reach a maximum of θ ZA = 40° [29].
FAST is currently equipped with a 19-beam L-band receiver with a frequency range of 1050-1450 MHz. The receiver temperature is around 7-9 K, and the measured system temperature T sys is around 20 K at 1400 MHz for all 19 beams. This includes the contribution from the galactic synchrotron radiation background off the galactic plane. Since this background is frequency dependent, as described above, we adjust T sys accordingly. Meanwhile, the system temperature will also increase linearly due to the emission and scattering of surrounding terrain if the zenith angle θ ZA > 20°, reaching 26 K at θ ZA = 40° [87]. The beam widths of the 19 beams are fit from values in Ref. [87].
For the other available frequencies that can be accessed with receivers on FAST, we will assume a single beam with receiver temperature T sys = 20 K (at 1400 MHz) and adjust the galactic synchrotron radiation for frequency. The aperture efficiency is assumed to be η A = 0.7, and the beam HPBW is taken as d = 1.22λ/D where D is the illuminated diameter 300 m.
The spectral feature is determined by the DM velocity dispersion, which we approximate here as a Gaussian with width σ d = 116 km/s. To achieve the optimal signal-to-noise ratio for a signal of this form, we adopt the signal-to-noise maximizing spectral line windowing as in Ref. [28], taking ∆ν = 2.17ν d σ d . For this choice of ∆ν, the fraction of the signal power within the window is f ∆ = 0.721. The signal power received by a single beam is thus wheren is the sky direction, I g (n) is the spatial intensity distribution of the (counter)source, and b(n) is the beam envelope, assumed to be Gaussians of width σ = θ HPBW and normalized to have a maximum value of unity. Over an integration of duration t obs , we wish to measure a statistically significant increase in the total power due to P S , relative to the thermal noise fluctuations in the telescope averaged over t obs . We denote this quantity σ N (though it is called P N in Ref. [28]): where k B is Boltzmann's constant, ∆ν is the bandwidth of the integration, and T sys is the system temperature which quantifies the search sensitivity. For the L-band receiver, the total signal-to-noise ratio of the 19 beams combined is the root mean square sum of that of each individual beams, ignoring correlated noise fluctuations between beams: where the index i runs over the beams.

B. Interferometers
In this Subsection, we describe the sensitivity of shortand long-baseline interferometers to SNR-induced axion gegenschein. As a representative of the compact mapping interferometers, we focus on CHORD, a proposed interferometer that consists a 24×22 rectangular array of ultra-wideband dishes that operate from 300 MHz to 1500 MHz. The distances between adjacent telescopes are approximately 9 m and 7 m, along the 24-site and 22-site direction, respectively. Each dish is 6m in diameter, and the total receiving area is 14400 m 2 . The aperture efficiency for each dish is taken to be η A = 0.5 and the system temperature is taken to be 30 K. Meanwhile, for the long-baseline interferometers, we consider the case of SKA1 (SKA Phase 1). SKA1 consists of SKA1low, which is an array of about 131000 antennas spread among 512 stations, and SKA1-mid, which is an array of 197 dishes. The operating frequency is 50-350 MHz for SKA1-low, and 350 MHz-15.3 GHz for SKA1-mid. We will project a signal-to-noise ratio using the single dish/station sensitivity information in Ref. [88], while taking into account the reduced power due to an extended source as follows. 3 Not all baselines in an interferometer will receive the full signal power from an extended source. The exact fraction of the total flux that remains unresolved on each baseline depends on the spatial structure of the gegenschein signal. Here, we make the simplifying assumption that on each baseline, the gegenschein signal is either completely unresolved (detectable) or resolved (undetectable). This is similar to assuming that the high spatial frequency modes of the gegenschein are washed out completely by the DM velocity dispersion, leaving no flux detectable on long baselines, and leaving all the flux detectable on short baselines. We take a baseline of length D to be usable when For SKA1, we assume that the baselines between stations in SKA1-low and baselines between dishes in SKA1-mid are not usable. More formally, the total signal power received by an interferometer can be expressed as by generalizing Eq. (22) to Fourier modes on the sky besides the monopole moment. However, we approximate this calculation by modifying Eq. (22): where N usable and N total denotes the usable and total number of baselines, respectively. However, the noise power received is the same as in Eq. (23). Therefore, the signal-to-noise ratio is where the sensitivity S is defined as For CHORD, we calculate the sensitivity with the aperture efficiency and the geometric area provided above. For SKA, we used the single station/dish sensitivity of Ref. [88], which gives a more detailed analysis of the frequency dependence of the effective area and noise temperature of the instrument.

V. FORECASTS
In Fig. 1 we show projected constraints assuming a null detection for the axion-photon coupling g aγγ using our fiducial model for SNR radio brightness history. We assume a detection threshold of a signal-to-noise ratio of unity. Bands correspond to the uncertainty in the projection coming from modeling choices (for instance, the treatment of the magnetic field amplification and relativistic electrons) and measurement uncertainty in quantities entering into the SNR flux evolution (the spectral index, age, and distance). To obtain the combined uncertainty estimates, at each frequency, we compute the smallest and largest value of constrained g aγγ if we vary one parameter from the fiducial model. Since the dominant uncertainties are systematic theory uncertainties, their effects on the projected constraints are degenerate with each other. The parameters and assumptions for the sources, Vela SNR, W28 and W50 are described in Sec. III C. One can see that despite astrophysical uncertainties, the null detection of SNR gegenschein images can place a competitive bound on the axion-photon coupling with existing telescopes like FAST.
In Fig. 4, we break down our uncertainty estimates based on individual parameters that enter in our radio brightness modeling and projection. While there are two parameters describing the SNR's brightness evolution, magnetic field amplification time t MFA and spectral index α, parameters such as the SNR's age and distance are . The "electron model" panel shows two lines for each candidate source, which corresponds to the different choice of electron energy evolution model employed in our brightness modelling. The more constraining bound corresponds to Eq. (16), our fiducial model, while the more conservative bound corresponds to Eq. (17). Despite uncertainties in the age and distance of the SNR, the uncertainty in the reach of our axion search is not dominated by these two parameters; instead it is dominated by uncertainties on the parameters relevant to the SNR radio brightness modeling: the magnetic field amplification time tMFA, source spectral index α, and the choice of electron evolution model. In all of the panels, we assume 100 hours of observing time on FAST.
also important in determining the final projection. One might reasonably expect that large uncertainties from astrophysical measurements would significantly affect the constraints; however, this turns out to not be the case. For example, the relatively large uncertainty in the age of W50 (3-10×10 4 years) affects the projected gegenschein power in two opposite ways: with a larger age we expect the initial brightness to be greater given our model, but also expect the distance of the echoing DM to be further, which means increased angular imaged size and decreased intensity. For W50, these two factors affect the total received gegenschein power in ways that nearly compensate, meaning that the final projected sensitivity is not very sensitive to uncertainties in SNR age. Similarly, SNR distance (given a fixed observed intensity) affects the projected constraint primarily through its effect on the expected gegenschein image size, and there are compensating effects between the increased SNR flux and increased gegenschein image smearing. The uncertainty is present mostly at high frequencies, where the beam of the telescope does not fully capture the gegenschein image and the expected image size affects the projection significantly. However, we still see that uncertainties in the measured distance do not significantly affect our projected sensitivity. As shown in Fig. 4, the biggest uncertainties stem from the modeling the radio brightness evolution of SNR. The relevant parameters are the magnetic field amplification time t MFA , the spectral index α which determines the brightness decay power law, and the choice of the electron energy model. The uncertainty on the spectral index mostly stems from variation in different part of the SNR, and not from the measurement uncertainty of a single spot, hence it is reflective of our model's simplification rather than measurement uncertainty. Better understanding of the radio brightness evolution of SNRs thus can greatly reduce the uncertainty on the projected constraint. Comparison of the constraining power of different telescopes using W50 (which is close to the equator and can potentially be seen by telescopes in both hemispheres) as a source, assuming 100 hours of observing time. Note that the increase in FAST's sensitivity around 1 GHz is due to its Lband 19-beam receiver. Despite having smaller illuminated area, single-dish telescopes like FAST can be more sensitive to extended sources like the gegenschein image of W50. Long baselines in interferometers such as CHORD or SKA1 only take in a fraction of the total gegenschein power, and thus these telescopes have much less effective area than they would in a point-source observation. This suggests that single-dish telescopes like FAST are best-suited for a search for axion gegenschein from individual SNRs. Note that W50 is not in CHORD's field of view, but is included for a sensitivity comparison.

VI. DISCUSSION
We have shown that the echoes of SNRs from stimulated axion decay may be detectable in the form of a spatially extended radio emission line coming from the antipodal direction of the SNR. Compared to other sources, nearby Galactic SNRs generate particularly bright axion echoes because their large temporal variation in brightness translates to a large spatial variation in the brightness of axion gegenschein along the DM column.
Nearby SNRs have a higher flux, which enhances stimulated decay. However, due to the DM velocity dispersion in the halo, nearby sources have their gegenschein images substantially blurred over a large angle that is parametrically enhanced for small source-observer distances.
Primarily for this reason, we find that single-dish telescopes like FAST and short-baseline interferometers like CHORD are best suited for constraining SNR-induced axion gegenschein as compared with long-baseline interferometers like SKA, as shown in Figure 5.
Making projections for the axion-induced signal required some astrophysical modeling of SNRs that are observed at present day in relatively late stages of their evolution. Guided by a mixture of theory and simulation, which were bolstered by population-level studies of SNRs (i.e. capturing SNRs with a wide range of ages), we back-evolved the observed synchrotron radiation of SNRs to a time roughly one hundred years after the supernova explosion making a series of conservative assumptions. We additionally estimated the size of uncertainties in our back-evolution, with the most important sources of uncertainty being the treatment of relativistic electrons, the timescale for magnetic field amplification, and measurement uncertainties of the spectral indices of SNRs at the present day. In spite of the conservative modeling choices we made and in spite of the uncertainties, we still find that it may be possible in the near future to constrain new axion parameter space with broadband radio observations. Some of the parameter space we can access has already been explored using terrestrial experiments. However, because our signal depends on a deep DM column, our projection is fairly robust to local deviations in the DM density, which may pose an issue for terrestrial experiments in a scenario in which the axions are largely bound into minihalos. We are also able to explore parameter space that may be interesting for stellar cooling anomalies, where axions with g aγγ ∼ few×10 −11 GeV −1 may provide an additional energy-loss channel beyond usual SM channels [89].
Note added: During the late stages of completing this manuscript, we became aware of a similar study [90]. We note that the approach and analyses of both works, while sharing common components, also differ and complement each other in several ways. early stages of this work, KS was supported by a Pappalardo Fellowship in the MIT Department of Physics and by NASA through the NASA Hubble Fellowship grant HST-HF2-51470.001-A awarded by the Space Telescope Science Institute, which is operated by the Associ-ation of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555. CL was supported by the U.S. Department of Defense (DoD) through the National Defense Science & Engineering Graduate Fellowship (NDSEG). KM was supported by NSF grants AST-2018490 and AST-2006911.