Constraints on axion-like particle properties with very high energy gamma-ray observations of Galactic sources

Axion-like particles (ALPs) can oscillate to photons and vice versa in electromagnetic fields. The photon-ALP oscillation provides an attractive solution to the apparent transparency of the Universe to TeV photons. The allowed parameter regions for the ALP mass $m_{\rm a}\leq 10^{-7}$ eV have been tightly constrained by the Fermi-LAT and H.E.S.S observations of some extragalactic sources. In this work we show for the first time that the H.E.S.S observations of some TeV sources in the Galactic plane exclude the highest ALP mass region (i.e., $m_{\rm a}\sim {\rm a ~few\times 10^{-7}}$ eV) that accounts for the TeV transparency of the Universe. The upcoming CTA observations of the Galactic TeV sources are shown to be able to improve the constraints significantly.


Introduction
The axion, a type of pseudoscalar particle beyond the standard model (SM), is postulated to solve a puzzle in quantum chromodynamics (QCD) known as "the strong CP problem" [1][2][3]. Axions, if exist and have low mass within a specific range, can be a possible component of cold dark matter (DM) [4,5]. More generalized axionlike-particles (ALPs), which have similar properties to axions, are predicted in several string-theory-motivated extensions of the SM as an alternative to axions to solve the DM problem [6][7][8][9]. Axions and ALPs have an interesting property that they can oscillate into photons and vice versa in external magnetic fields via the Primakoff process [10]. The Lagrangian of the interaction can be written as L = g aγ E · Ba, where E, B and a are electric, magnetic and axion (or axion-like) fields, and g aγ is the coupling constant. The axion mass m a is proportional to the coupling constant g aγ , while for ALPs these two parameters are not correlated any longer.
In the literature, ALPs have been adopted to explain the observed transparency of the Universe for TeV photons [29][30][31][32][33]. TeV photons may interact with extragalactic background light (EBL) in their routes to the Earth, thus the Universe is opaque to sufficiently distant TeV sources. The TeV opacity is governed by the source distance, EBL photon density and the photon energy. Interestingly, some studies showed that the opacity of the Universe is lower than the expectation [30,[34][35][36][37] (see however [38][39][40] for different opinion). The TeV -1 -

JCAP06(2019)042
transparency can be easily explained if ALPs exist [29][30][31][32][33]. TeV photons can oscillate into ALPs and vice versa in the magnetic fields near the TeV source. Once produced, ALPs will unimpededly pass through long distances, then convert back to TeV photons before reaching the earth, thus lower the opacity of the universe. The ALP parameter space accounting for the transparency has been calculated in [41]. Both the aforementioned laboratory experiments and the γ−ray observations of some extragalactic sources have effectively narrowed down the parameter space. In particular, with the γ-ray spectrum of NGC 1275, the Fermi-LAT collaboration has excluded the coupling g aγ above 5 × 10 −11 GeV −1 for ALP masses 0.5-5 neV [23] (see also [28] for the additional constraint set by the Fermi-LAT observation of PKS 2155-304). Moreover, the H.E.S.S. observation of PKS 2155-304 has yielded an upper limit of the coupling g γa < 2.1 × 10 −11 GeV −1 for m a ∼ 15-60 neV [19] with the conservative galaxy cluster magnetic field model. Nevertheless, the high end part of the allowed parameter space (i.e., m a > 10 −7 eV) for the TeV transparency has not been independently probed by the gammaray data, yet. To surpass the CAST limits (i.e., g aγ < 6 × 10 −11 GeV −1 ), we need to observe the photons at the energies of ∼ 2 TeV (m a /10 −7 eV) 2 (g aγ /5 × 10 −11 where B T is the strength of the large scale magnetic filed. Therefore the TeV observations of the sources in the Galactic plane are ideal probe of the photon-ALP oscillation effect for m a ∼ 10 −7 eV. In this work we focus on the bright H.E.S.S sources in the Galactic plane and search for the potential irregularities in the gamma-ray spectra caused by the photon-ALP oscillation. 2 Photon-ALP oscillation in the Milky Way magnetic field ALPs can convert into photons (and vice versa) in the external magnetic fields. For a coherent magnetic field with size l and a transversal strength B T , the survival probability for an initially polarized incoming photon with energy E γ is approximated as [10,42] P ALP = 1 − P γ→a (2.1) The oscillation behavior becomes most significant near the characteristic energy E c which is defined as where w 2 pl = 4παn e /m e is the plasma frequency. The Milky Way magnetic fields are of the order of B T ∼ 1 µG, and the thermal electron density n e ∼ 0.1 cm −3 . With ALP mass m a = 10 −7 eV and the coupling g aγ ∼ 10 −10 GeV −1 (the current upper limit), we have E c ≈ 1 TeV, at which the H.E.S.S is sensitive. In our analysis, rather than directly adopting eq. (2.1), we solve the evolution equation for photon-ALP beam as in refs. [10,42] to calculate the survival probability. For simplicity, the source photons are assumed to be unpolarized and conservative results are yielded.
The Galactic magnetic field (GMF) consists of a large-scale regular component and a small-scale random component. In our work, we neglect the later since its coherence length is much smaller than the photon-ALP oscillation length. For the regular magnetic -2 -JCAP06(2019)042 field component, we take the model developed by Jansson & Farrar [43] (Bfield1) which has been widely adopted in ALP researches [23][24][25][26][27][28]. In this paper, the results are based on Bfield1 unless particularly specified. Two other GMF models developed by Sun et al. [44] (Bfield2) and Pshirkov et al. [45] (Bfield3) have been employed to check the uncertainties of the constraints. We use the modified parameters as updated with the polarized synchrotron and dust emission data measured by Planck satellite [46]. The updated parameters are b disc 6 = −3.5, B X = 1.8 for Bfield1 and B c = 0.5 µG for Bfield2. The magnetic fields inside the sources themselves may cause the spectral oscillation which is different from that induced by GMF. For the representative PWN among our sample, HESS J1825-137, the magnetic field and electron density in the source region are B = 100 µG and n e = 10 −2 cm −3 , with a scale of R = 35 pc [47]. For the SNRs, the parameters of HESS J1640-465 are B = 35 µG [48], n e = 130 cm −3 and R = 55 pc [49]. In both cases, we find the oscillation effect induced by the magnetic fields inside sources is negligible in the parameter space of our interest (g aγ < 10 −10 GeV −1 , 1 neV < m a < 10 4 neV). So we don't take the magnetic fields in the source region into account.

The source sample and methods
High Energy Stereoscopic System (H.E.S.S.) is a system of Imaging Atmospheric Cherenkov Telescopes that observe gamma rays in the energy range from tens of GeV to tens of TeV [58]. In this work we focus on the H.E.S.S. observation of bright Galactic sources, though some other devices, such as VERITAS [59], MAGIC [60] and HAWC [61], are also sensitive to the TeV γ-rays. We collect the H.E.S.S. observation spectra from the literature. 1 The oscillation effect due to the ALP-photon conversion is significant only when the γ-rays have traversed considerable magnetic fields. To increase the amplitude of the oscillation effect, we limit our sample to those with distances greater than 1 kpc. In addition, to guarantee sufficient statistics, we select only the bright sources. Totally 10 sources are considered in our work. Most of them are supernova remnants (SNRs) or pulsar wind nebulae (PWNe). The basic information of these sources are listed in table 1.
To examine whether the oscillation signals exist, we fit each H.E.S.S. spectrum with two types of models, the background model (the null model) and the ALP model (the alternative model). The background model is for the case that intrinsically no ALP exists and these TeV sources can be well fitted with smooth astrophysical spectra. The signal model is the superposition of the ALP effect on the background spectrum (i.e., the modulation on the spectrum).
In this work, the intrinsic spectrum is described by two types of functions, a power law (PL) or an exponential cutoff power law (ECPL), where N 0 , α, β and E cut are independently set free in both fits for the background and ALP models. For each source, the function type of the intrinsic spectrum is chosen as that in the literature [48,[50][51][52][53][55][56][57]. The spectral types and the corresponding parameter values JCAP06(2019)042  reported therein are also presented in table 1. It is worth noting that the spectra adopted here are only the phenomenological descriptions to the observations. The true spectra of the sources may be not exactly described by them. For our purpose, the PL or ECPL assumption, widely adopted in literature [e.g., 23, 25-28], is reasonable/sufficient. For the ALP model, the spectrum is expressed as where the P ALP (g aγ , m a , E) is the survival probability.
The energy dispersion of ground based Cherenkov telescope like H.E.S.S. is usually large and thus should be considered in the ALP analysis. We approximate the energy dispersion function to be a Gaussian distribution with its σ being the energy resolution of the instrument. The H.E.S.S. energy resolution reported in the literature ranges from 10% to 20% [19,[62][63][64][65]. In our analysis, we adopt the energy-dependent energy resolution presented in ref. [63]. Two caveats should be made about the energy dispersion we use in the analysis. First, the energy dispersion will depend on the zenith angle of the observation which is not considered in our analysis. Second, the dispersion in ref. [63] only applies to one of the two independent analysis chains adopted within the H.E.S.S. collaboration which might not be the chain with which the spectra in the respective publications are derived [66]. However, we will show in section 5.2 that the inaccurate use of the energy resolution within 10% to 20% range has only minor effects on the results.
After convolving the energy dispersion function D(E , E), we get the final model ex- The χ 2 fit is utilized to analyze the spectrum of each source. The χ 2 value is defined as where E i is the geometrical central energy,φ i and δ i are the observed flux and its uncertainty in the i-th bin, respectively. To validate the χ 2 fit, we only consider the bins withφ i > 3δ i . 2 Smaller χ 2 indicates the selected model can fit the observation better, thus whether the ALP model is favored or not can be determined through the quantity ∆χ 2 = χ 2 wALP − χ 2 w/oALP . For deriving the exclusion region of the ALP parameters, we define λ(m a , g aγ ) as the difference between the best fit χ 2 value for each set of ALP parameters (m a , g aγ ) and the best fit over the whole parameter space. λ thr is the threshold λ value, above which the ALP parameters are excluded at a 95% confidence level (C.L.). The threshold value is set to be λ thr = 13.9 in our analysis according to Monte-Carlo simulations (see appendix A for details).

Results
With the procedure described above, we fit the spectrum of each source with two sets of spectral models. For the ALP model, we scan a grid of m χ and g aγ and derive the ∆χ 2 for each set of parameters. In the left panel of figure 1, the spectrum of a representative source, HESS J1640-465, is shown, together with the best-fit background model (red line) and ALP model (green line). The map of ∆χ 2 value as a function of ALP mass m a and photon-ALP coupling constant g aγ is presented in the right panel of figure 1. A negative ∆χ 2 means that the ALP model results in an improved fit to the observed spectrum than the background one. Meanwhile, a positive ∆χ 2 suggests that the fit with ALP effect is worsen, indicating such a hypothesis is disfavored. As is shown, for part of the ALP parameters (blue region), the inclusion of the modulation effect does improve the goodness of fit. The highest TS value of the ALP signal is TS∼7.8 for HESS J1731-347 (the test statistics, i.e. TS, is defined as TS = −∆χ 2 ). According to the Monte-Carlo simulation in appendix A, which shows that the null distribution can be described with a χ 2 distribution with ∼ 6.7 degree of freedom (d.o.f.), this TS value corresponds to a significance of only ∼ 1.0σ. Combining 10 sources together, we get TS = 10.3 which corresponds to a low significance level of 1.4σ. Moreover the most favored parameters (m a ∼ 6 × 10 −7 eV, g aγ ∼ 5 × 10 −10 GeV −1 ) are excluded by 2 According to Poisson statistics, this condition requires the photon number Ni 10.  [19]. The light blue region is the parameter space where the low gamma-ray opacity of the universe can be explained by the ALPs [41]. The benchmark result with Bfield1 model derived in this work is shown as yellow.
In the plot, we also present the projected exclusion region assuming 50-hour observations of the 10 Galactic sources by future CTA mission. The two dotted lines are for different Milky Way magnetic field models (red: Bfield2, blue: Bfield3). Also plotted is the upper limits of the photon-ALP coupling set by CAST [15] experiment (horizontal line). the CAST limit [15] (horizontal line in figure 1). Thus, in the following, we focus on the exclusion of parameter region (roughly the red region in the ∆χ 2 map of figure 1). Below the CAST limits, each single source can not impose stringent constraints on ALP properties (see appendix B for ∆χ 2 maps of other 9 sources).
By summing the ∆χ 2 maps of the 10 sources, we set the combined constraints on the ALP parameters. The results are shown in figure 2. The yellow region is excluded by the combined analysis at 95% confidence level. In comparison with the previous constraints by other targets and instruments [19,23,28], our results narrow down the ALP properties at higher m a .
The above results are based on the Galactic magnetic field model of Jansson & Farrar [43]. The constraints are more stringent for the model of Sun et al. [44] and the model of Pshirkov et al. [45]. The stronger constraints of these two GMF models are mainly due to the higher magnetic strength in the disk component, though for some sources, e.g. HESS J1640-465, the Bfield1 model has a better field orientation for the ALP-photon mixing. The regions, excluded by these two models, cover the high ALP mass space that interprets the TeV transparency (the light blue region in figure 2).

Validating the fitting method
The raw event data of H.E.S.S. are not publicly available, we use the published spectral points in the literature to carry out our analysis. However, the actual H.E.S.S. analysis is    based on a forward folding technique, which convolves the source spectrum with the energy dispersion of H.E.S.S. to fit the raw data in the entire energy range. The different fitting method used here may introduce additional bias. To test this, we compare our best-fit parameters without axion-like particles to the best-fit spectra in the H.E.S.S. publications. The comparisons are shown in figure 3. The left and right panels are for the normalization N 0 and the spectral index α respectively. We find that the parameters derived with our method are well consistent with those reported in the literature, demonstrating that our method is valid to derive the spectral parameters. We would like to point out that, for many sources in the sample, the best fits without ALPs lead to χ 2 /d.o.f. significantly smaller than 1 (i.e.,

Related uncertainties
Our analysis suffers from some systematic uncertainties. Among them, the magnetic models and the flux measurement uncertainty likely play the most important roles.
Galactic magnetic field. To test how our results depend on the choice of magnetic model, we have taken into account three GMF models developed by Jansson & Farrar [43] (Bfield1), Sun et al. [44] (Bfield2) and Pshirkov et al. [45] (Bfield3). As shown in figure 2, different GMF models would lead to very different exclusion regions. The improvement comparing to the CAST limits is marginal for the result of Bfield1, while the other two models yields much more stringent constraints on the ALP parameters. However, Bfield1 is preferable over the other two models in many respects. First, Jansson & Farrar [43] made use of vastly more Faraday rotation measures (RM) data, which are also binned into much finer RM pixels, in their work of deriving Bfield1. They included an additional out-of-plane component, which has been observed in external edge-on galaxies and can significantly improve the model fit. They also enforced that the field is divergenceless. Second, the Bfield1 is obtained by fitting both RM and polarized synchrotron radiation (PI) observations, while the later is not used when deriving Bfield3. The PI observations are crucial to determine the transversal magnetic field, that is the component contributing to ALP-photon mixing. In addition, in the procedure of  Figure 4. Left: the influence on the constraints due to the adopted energy resolution in the analysis. Here we show the constraints for three constant energy resolution of H.E.S.S. The benchmark result using the energy-dependent energy resolution in [63] is shown as red line. Right: the results taking into account a 20% uncertainty on flux measurements on basis of the energy resolution of [63]. The shaded area is for our benchmark GMF model (i.e. Bfield1).
deriving Bfield3 model, a large part of the Galactic plane region was masked [45], thus may bias the determination of the magnetic field model in the region where our target sources locate. Finally, quantitative comparisons show that Bfield1 model leads to vastly better figure of merit than others for the same set of RM+PI observations (χ 2 red = 1.096, 1.672 and 2.633 for Bfield1, Bfield2 and Bfield3, respectively, for nearly 6600 degree of freedom) [43,67]. Thus the fiducial Galactic magnetic field model adopted in our main analysis (i.e., Bfield1) is preferable. Nevertheless, it should be noticed that the magnetic field distribution in the Galaxy is still to be better understood, and the constraints on the ALP properties may be subject to change in the future.
Energy resolution. The energy resolution of H.E.S.S. reported in the literatures ranges from ∼10% to ∼20% [19,[62][63][64][65]. In this work, we use the energy-dependent energy resolution of H.E.S.S. presented in ref. [63] to derive the main results. It is straightforward to speculate that, for searching for the fine structure in spectra deduced by ALP effect, inaccurate use of the energy resolution may cause systematics to the results. We examined such an effect via adopting three constant values of energy resolution, 10%, 15% and 20%. Figure 4 shows the results. The modification of the energy resolution between 10% and 20% has only marginal effects on the constraints and our fiducial result shown in red is quite consistent with these constraints.
Energy binning. The spectra used by us are collected from the literature, we could not examine the systematics associated with the energy binning. As is shown in ref. [19], which also focused on studying ALP signals with H.E.S.S. data, varying the binning does lead to a certain level of fluctuations on the irregularity caused by the ALP modulation. However, the method used in their analysis, which measures the level of irregularity in the spectrum of extragalactic source, is not the same as that in this work. In our previous work of adopting also χ 2 fit to analyze the Fermi-LAT data of the Galactic sources, we find that the energy binning has only a minor influence on the results [26].   [63]. Right panel: the survival probability for the source HESS J1640-465 for ALP mass of 150 neV and assuming the distance to be 8.6 kpc or 13 kpc. The oscillation behaviors differ substantially for two assumed distances, but the magnitudes of the prominent structures around 1 TeV do not differ too much, which may interpret why the constraints do not depend on the choice of distance significantly. All these results are for Bfield2.
Flux measurements. A number of instrument effects (e.g. effective area) may generate systematic uncertainties in the flux measurements. In this work, we use a simple way to estimate how such uncertainties would influence the final constraints. We simply assume the total systematic uncertainty on the flux measurements is 20% [62] and add it in quadrature to the statistical errors. With this additional flux uncertainty included, we get weaker constraints on the ALP parameters as expected (see right panel of figure 4). The parameter space below the CAST limit can not be probed any more if taking into account the 20% systematic uncertainties for Bfield1 model.
Distance of the sources. The distance of the object is another parameter to affect the final constraints. The measured distances for the Galactic TeV sources suffer from large uncertainty. Therefore, we have also tested the systematic uncertainty due to choice of distance. To do that, we search for the distance uncertainty information in TeVCat, 3 SNRcat 4 and the references listed in table 1. We find the uncertainty information for six of the ten sources, which are listed in the fifth column of table 1. For other sources, which are lack of such information in these catalogs and references, the uncertainty is assumed to be 20%. Then, we repeat the same analysis procedure as that in the main text 100 times. For each, a set of distances for the 10 sources are uniform-randomly assigned within the uncertainty ranges. In figure 5, we demonstrate the 2σ coverage due to varying the distances. We find that the effects are not sizeable. To understand why, we plot the conversion probability for HESS J1640-465 in the right panel of figure 5 assuming distances of 8.6 kpc and 13 kpc, respectively. As is shown, though the resulting conversion probabilities vary substantially for two assumed distances, the magnitudes of the most prominent oscillation structures around 1 TeV do not differ too much. For the corresponding ALP parameters, such a structure we believe plays the leading role in constraining the ALP model.

Testing the effect induced by spectral variation among sources
In this work, we employ two phenomenological function (PL or ECPL) to model the intrinsic spectra of the sources. As stated before, the intrinsic spectra may be not exactly described by them. Especially, due to the spatial extension of some sources, the gamma-ray photons may come from several different sub-zones. Also, some sources (e.g. X-ray binary) are time variable. The sum of radiation from different regions/phases may cause the irregularity in the intrinsic spectra.
Here, we take HESS J1825-137 as an example to discuss how the temperally or spatially spectral variation would influence the final results. The HESS J1825-137 has an extension of > 1 • and is potentially the largest gamma-ray PWN currently known [54]. It is also one of the brightest sources in the Milky Way at the VHE range. As a result of the cooling of the electrons, the spatially resolved spectra are softened with increasing distance from the pulsar [54]. The spectrum of HESS J1825-137 used in this work is in fact the one averaged over the full PWN.
In figure 6, we use spectra of HESS J1825-137 from two different regions to perform the ALP analysis, one is for the radiation within the core 0.4 • region and another is for the full HESS J1825-137 emission. The later can be treated as an overlap of the spectra of the core and outer (> 0.4 • , which is not presented in the literature) regions. By comparing the two results, we find that the average spectrum gives wider exclusion region and does not introduce any artificial signal (no significant blue region in the right panel of figure 6).
We would like to suggest that although in this case we don't find any evidence of the artificial signal caused by averaging the spectra over the entire sources, it is still possible for the superposition of the spectra to account for the observed oscillation structure, as suggested in [26].

CTA perspective
The upcoming next generation Cherenkov telescope, CTA [68], with significantly enhanced sensitivity from 20 GeV to 300 TeV, will make a great improvement on the ALP searches. Some previous predictions on CTA sensitivity for ALP detection have been studied in refs. [69,70]. Here we perform Monte-Carlo simulations to derive the projected exclusion regions that CTA can provides through observing Galactic sources.
We use ctools [71] to simulate event data and generate the Monte-Carlo SEDs of the 10 sources. The ctools is a software package developed for the scientific analysis of Cherenkov Telescope Array (CTA) data. We simulate events from 50-hour CTA observation for each source. The South 50h CTA instrument response function (IRF) [72] taken from the latest prod3b calibration database 5 is used in order to match the fact that the 10 sources are in the southern sky. The input spectral types and parameters of the 10 sources are based on the best-fit null model obtained in the above analyses of H.E.S.S. spectra. We add a IRFbackground component into the model file to consider the cosmic-ray background. After the pseudo event data are generated, we perform a binned analysis using also the ctools package to derive the SEDs. The data are binned into 200 × 200 spatial bins with a bin size of 0.02 • and 30 logarithmically spaced energy bins from 0.2 TeV to 50 TeV. We do not take into account the data outside this energy range where CTA is also sensitive to avoid extrapolation of the spectra. We take the energy dispersion into account for both simulating events and deriving SEDs by set edisp=yes. The details related to the usage of the ctools 5 https://www.cta-observatory.org/science/cta-performance/. can be found in the ctools website. 6 One realization of the Monte-Carlo spectrum of HESS 1640-465 from such a procedure is shown in figure 7.
Once the Monte-Carlo SEDs are prepared, the same χ 2 analysis procedure is used to fit the simulated spectra to derive the combined ∆χ 2 map for each simulation. In the ALP analysis of fitting the SEDs, the energy dispersion function is approximated by a Gaussian distribution with its σ being the energy-dependent energy resolution of the CTA. Totally, 100 simulations are performed. The mean exclusion regions derived from Monte-Carlo simulations for 50-hour CTA observations are plotted in figure 2 (gray region). Please note that our CTA projections do not take into account any systemic errors that would dominate the uncertainty at the low energy. Therefore our results likely have overestimated the prospect.
Our simulation shows that if no ALP signal is found in the future CTA observations of these bright Galactic gamma-ray sources, the allowed region for the ALP model, especially those can account for the low opacity of the Universe, will be further constrained. It is also reasonable to speculate that the CTA observation of some extragalactic sources, such as PKS 2155-304, would notably improve the current limits around the ALP masses of m a ∼ 10 −8 -10 −7 eV.

Summary
The ALP is one promising type of cold dark matter candidate, which can also solve the puzzle of the TeV transparency of the Universe. In this work, for the first time, we constrain the ALP parameters with H.E.S.S. observations of Galactic sources. We suggest that the Galactic TeV sources, which are usually observed by ground based Cherenkov telescope, can set constraints on the ALP parameters stronger than the CAST limits, though the constraints are sensitively dependent on GMF models. In this work, we also perform Monte-Carlo simulations to show that the next generation Cherenkov telescope, CTA, can probe a wider region of the parameter space. We expect that the ALP interpretation of the low opacity of the universe will be unambiguously tested in the near future.
Note that some on-orbit/future space-borne instruments, such as the Dark Matter Particle Explorer (DAMPE) [73,74] and High Energy cosmic-Radiation Detection Facility (HERD) [75], are also sensitive to the photons in TeV energy range. With the significantly higher energy resolution, they would help us better understand the ALP properties. Together with the upcoming ALP-II [76] and IAXO [77] experiments, significant progresses on the probe of ALP are expected in the next decade.

A Null distribution and 95% confidence level
We define the quantity λ(m a , g aγ ) as the ∆χ 2 between the best fit with given ALP mass and coupling (i.e., only fit the nuisance parameters) and that with all parameters free to vary, λ(m a , g aγ ) =χ 2 (ma,gaγ ) −χ 2 global . Usually, the 95% confidence level (C.L.) exclusion region corresponds to the parameters leading to an increase of theχ 2 by 5.99 compared with the best fit (i.e., λ > λ thr = 5.99 in this case for two free parameters). However, due to the non-linear dependence of the spectral irregularities on the ALP parameters and possible systematics in the observations, this threshold value λ thr may be biased. In this work we adopt the λ thr derived from Monte-Carlo simulations to set the 95% C.L. limits.
The λ thr values are in fact different for different sets of ALP parameters (m a , g aγ ) and we need to simulate the distribution of λ(m a , g aγ ) (alternative distribution) for the complete parameter space, which is however computational expensive. Following ref. [23], in this work we assume that the probability distribution of the alternative hypothesis can be approximated with the null distribution. It is found that such an assumption would yield conservative limits [23].
Accurate simulations of the observed spectra of the 10 sources require the H.E.S.S. exposures to these sources, which are however unknown. We therefore generate the simulated spectra in the following way. For each source, the simulated spectrum is set to have the same energy bins as the observed one. In each energy bin, the nominal value of the flux is randomly generated based on the Gaussian distribution with its mean being the model expected flux calculated by the best fit null model and the sigma being the uncertainty of the observation. The error bar of the flux in the simulated spectrum is set to the same value as the observed one. With this method, we generate 500 sets of the Monte-Carlo null spectra of these 10 sources. For each set of the simulated spectra, we perform the same analysis as used in the main article to derive the best fitχ 2 for both null model and ALP model. The distribution of the T S =χ 2 null −χ 2 wALP for 500 simulations is shown in figure 8. A non-central χ 2 fit to the null distribution results in degree of freedom (d.o.f.) d = 6.7 ± 0.2 and non-centrality s ≈ 0 (indicating that a standard χ 2 function can fit the distribution well). From the best-fit JCAP06(2019)042 distribution, we derive the threshold λ thr = 13.68 ± 0.22 above which the ALP parameters could be constrained at 95% confidence level. The 1σ error bar is due to the limited number of our simulations and is derived from bootstrapping the null distribution 10 4 times. For the sake of conservativeness, the upper bound λ thr = 13.90 is used to set our limits.

B Individual search results
The maps of ∆χ 2 as a function of ALP mass m a and photon-ALP coupling constant g aγ for other 9 sources of our sample (except HESS J1640-465, which has been presented in the main text) are plotted in figure 9. w/oALP Figure 9. The maps of ∆χ 2 as a function of ALP mass m a and photon-ALP coupling constant g aγ for 9 sources among our sample. The upper limit of g aγ set by CAST [15] experiment is also plotted as a reference (dashed horizontal line).