Constraints on axion-like particles from the observation of Galactic sources by the LHAASO

Axion-like particles (ALPs) [1–4], a class of pseudo-scalar bosons, arise as a consequence of symmetry breaking in many extensions of the standard model. ALPs possess a broader parameter space and a rich phenomenology that is yet to be fully explored, compared to quantum chromodynamics axions addressing the strong CP problem [5–8]. The effective coupling between ALPs and photons can lead to ALP-photon oscillation in an external magnetic field. This phenomenon has attracted significant attention in astrophysics owing to the ubiquitous astrophysical magnetic fields [9–43].

The oscillation between ALPs and photons has the potential to induce irregularities in the gamma-ray spectrum. Detecting this phenomenon is typically more achievable at lower energies because detectors exhibit superior energy resolution for low-energy photons compared to high-energy photons. Nevertheless, the impact of ALP-photon oscillation on the astrophysical gamma-ray spectrum may also become pronounced at high energies. Interactions involving high-energy photons and low-energy background photons, such as those originating from the interstellar radiation field (ISRF) [44–46], cosmic microwave background (CMB) [47], and extragalactic background light [48–51], lead to the absorption of high-energy photons, thereby attenuating the observed gamma-ray spectra at high energies. The conversion of photons to ALPs may mitigate this absorption effect because of the interaction between ALPs and high-energy photons. Consequently, the presence of ALPs leads to a modification of the expected gamma-ray spectrum at high energies within the standard astrophysical framework. The observation of these effects would serve as evidence for the existence of ALPs. However, because there is no positive signal of such effects, we can only establish constraints on the ALP parameters based on the observation.

In recent years, significant progress in high-energy gamma-ray observation experiments has led to remarkable measurements of high-energy gamma-ray spectra [52–61]. Notably, the Large High Altitude Air Shower Observatory (LHAASO) [62] has contributed significantly to this field. In 2021, the LHAASO collaboration reported the detection of ultra-high energy gamma-rays from the Crab Nebula [53]. This measurement encompassed results from two detectors, namely, the Water Cherenkov Detector Array (WCDA) and Kilometer Square Array (KM2A), offering a precise gamma-ray spectrum of the Crab Nebula that spans more than three energy orders, from $ 500\; {\rm{GeV}} $ to $ 1.1\; {\rm{PeV}} $. In the same year, the LHAASO collaboration reported the detection of over 530 photons with energies above $ 100\; {\rm{TeV}} $ and up to $ 1.4\; {\rm{PeV}} $ from 12 ultra-high-energy gamma-ray sources [52]. The energy spectra of four Galactic sources, namely, the Crab Nebula, LHAASO J2226+6057, LHAASO J1908+0621, and LHAASO J1825-1326, were provided in the same report [52]. These measurements of high-energy gamma-ray spectra offer a promising opportunity to investigate the ALP-photon oscillation effect.

In this study, we utilize the LHAASO observations from the Crab Nebula, LHAASO J2226+6057, LHAASO J1908+0621, and LHAASO J1825-1326 to impose constraints on the ALP parameters. We consider the absorption effect induced by CMB and ISRF photons on ultra-high-energy photons within the Milky Way. Although several studies [30, 31, 39, 40] derived constraints on the ALP parameters from high energy gamma-ray observations of Galactic sources, we emphasize that the photons from the four Galactic sources considered in this study are more energetic. Furthermore, we employ the CL_s method [63, 64] to establish robust constraints on the ALP parameters, as in [37, 38]. To further enhance the constraints, we also conduct a combined analysis, incorporating Crab observations from many other experiments, including HAWC [65], $ {\rm{AS}}{\gamma} $ [66], HEGRA [67], MAGIC [68], HESS [69], and VERITAS [70].

Recently, the LHAASO collaboration reported updated selection criteria for photons and presented a new Crab spectrum within the energy range 10−1000 TeV using KM2A data from August 2021 to August 2022 [71]. We treat this measurement as independent and incorporate it into the combined analysis.

This paper is organized as follows. In Sec. II, we introduce the ALP-photon oscillation effect and the absorption of high-energy photons in the Milky Way. In Sec. III, we describe the process of fitting the gamma-ray spectra and the CL_s method. In Sec. IV, we present the constraints on the ALP parameters from the LHAASO observations of four Galactic sources and the combined constraint from the observations of the Crab Nebula from multiple experiments. Finally, we conclude in Sec. V.

In this section, we provide a brief introduction to the ALP-photon oscillation effect in Galactic magnetic fields. The interaction between an ALP and photons can be described by the Lagrangian term

where $ g_{a \gamma} $ denotes the coupling coefficient between the ALP and photons, a denotes the ALP field, F denotes the electromagnetic field strength tensor, $ \tilde{F} $ denotes its dual tensor, $ \vec{E} $ denotes the photon electric field, and $ \vec{B} $ denotes the magnetic field. The propagation equation for a monochromatic ALP/photon beam can be written as [9]

where z denotes the distance along the propagation direction $ \hat{z} $, E denotes the energy of the ALP/photon, and ${\Psi}\equiv {(A_\perp,A_\parallel,a)}^{\rm T}$ , with $ A_\perp $ and $ A_\parallel $ representing the photon polarization amplitudes perpendicular and parallel to the transverse component of the external magnetic field $ B_t $, respectively. The matrix $ {\cal{M}}^{\prime} $ encompasses the ALP-photon oscillation effect and the absorption effects of high-energy photons and can be written as

The mixing matrix $ {\cal{M}} $ includes the interaction between photons and ALPs and environment effects and is expressed as

where $ \Delta_{\perp} = \Delta_{\rm{pl}} +2\Delta_{\rm{QED}} $, $ \Delta_{\parallel} = \Delta_{\rm{pl}} + 7/2\Delta_{\rm{QED}} $, $ \Delta_{a}=-m_a^2/ (2E) $, and $ \Delta_{a\gamma}=g_{a\gamma}B_t/2 $. Here, $ m_a $ is the mass of the ALP. The term $ \Delta_{\rm{pl}} = -\omega_{\rm{pl}}^{2}/(2E) $ describes the effective mass of photons in plasma with the typical frequency $ \omega_{\rm{pl}} $. $ \Delta_{\rm{QED}}=\alpha E/(45\pi)(B_{t}/B_{\rm{cr}})^2 $ is the QED vacuum polarization term, where α is the fine structure constant, $ m_e $ is the electron mass, and $ B_{\rm{cr}}\equiv m_e^2/|e| $ is the critical magnetic field. The off-diagonal element $ \Delta_{a\gamma}=g_{a\gamma}B_t/2 $ describes the ALP-photon mixing effect.

When ultra-high-energy photons propagate in the Galactic magnetic field, the absorption effect induced by CMB and ISRF photons through the pair production process $ \gamma + \gamma_{\rm{bkg}} \to e^+ + e^- $ [45, 72, 73] cannot be neglected. $ \Gamma_{\rm{BG}} $ representing these effects is calculated as

where $ E_{\rm{BG}} $ and $ n_{\rm{BG}} $ are the energy and number density of background radiation fields, respectively. In this analysis, we utilize the ISRF model presented in Ref. [74]. The term $ \hat{\sigma} $ is given by

where $ x\equiv 1-{\rm{cos}}{\theta_{\gamma\gamma}} $, and $ \theta_{\gamma\gamma} $ is the angle between incident photons. The cross section of the pair production $ \sigma_{\gamma\gamma} $ is given by

where $ \beta \equiv(1-4m_e^2/s)^{1/2} $, $ \sigma_{T} $ is the Thomson cross section [75, 76], and $ s=2xEE_{\rm{BG}} $ is the center-of-mass energy. For photons from the Crab Nebula with energies of $ 1\; {\rm{PeV}} $, this absorption would result in a loss of $ \sim $ 19% of the photon flux.

The generalized density matrix $ \rho \equiv \Psi \otimes \Psi^\dagger $ can be used to describe the polarized states of the ALP-photon system. This matrix ρ obeys the Von Neumann-like equation [10, 14]

The solution to Eq. (8) in a homogeneous magnetic field can be expressed as $ \rho(z) = {\cal{T}}(z) \rho(0) {\cal{T}}^\dagger(z) $, where the transfer function $ {\cal{T}}(z) $ is obtained from the solution of Eq. (2). High energy photons emitted from Galactic sources encounter the magnetic field of the Milky Way before reaching the Earth. The entire path can be divided into many pieces, with the magnetic field in each domain considered to be homogeneous. The total transfer matrix is given by

where $ {\cal{T}}_i (\Delta z_i) $ represents the transfer matrix in the i-th piece.

In this study, we focus on Galactic sources and investigate the ALP-photon oscillation effect in the Galactic field. The Galactic magnetic field consists of a regular component and turbulent component. As the latter component is small and can be safely ignored, we only consider the regular component. The Galactic magnetic model utilized in this study is the Jansson & Farrar model [77]. Additionally, we take the NE2001 model [78] for the Galactic electron distribution.

The survival probability of the photon can be expressed as [9, 13]

where $ \rho(0) = {\rm{diag}}(1/2, 1/2, 0) $, $ \rho_{11}={\rm{diag}} (1, 0, 0) $, and $\rho_{22}= {\rm{diag}} (0, 1, 0)$ for unpolarized photons. Considering the ALP-photon oscillation and the absorption effects described above, the observed photon energy spectrum is given by

where $ \left.\frac{{\rm{d}}N}{{\rm{d}}E}\right\vert_{\rm{int}} $ is the intrinsic spectrum of the source. For the sources considered in this study, we set their intrinsic spectra to be a log-parabolic function given by $ F_0(E/E_0)^{-\Gamma - b{\rm{log}} (E/E_0)} $, where $ F_0 $, Γ, and b are free parameters, and $ E_0 $ is set as $ 10 \; {\rm{TeV}} $. This choice is based on the results in Refs. [52] and [53], where the log-parabolic spectrum provides better fitting results.

In this section, we briefly introduce the analysis method used to set constraints on the ALP parameters. The best-fit spectrum can be obtained by minimizing the $ \chi^2 $ function

where $ \chi^2_{j} $ denotes the $ \chi^2 $ function of the j-th source. $ \chi^2_{j} $ is given by

where $ \widetilde{\Phi}_i $, $ \Phi_i $, and $ \delta \Phi_i $ represent the predicted value, observed value, and experimental uncertainty of the photon flux in the i-th energy bin, respectively.

For given $ m_a $ and $ g_{a\gamma} $, we define the test statistic (TS) as

where $ \chi^2_{\rm{Null}} $ represents the best-fit $ \chi^2 $ value under the null hypothesis without the ALP-photon oscillation effect, $ \chi^2_{\rm{ALP}} $ represents the best-fit $ \chi^2 $ value under the alternative hypothesis including the ALP-photon oscillation effect with the given two parameters $ m_a $ and $ g_{a\gamma} $, and $ (\hat{F_0}, \hat{\Gamma}, \hat{b}) $ and $ (\hat{F_0}^{\prime}, \hat{{\Gamma}^{\prime}}, \hat{{b}^{\prime}}) $ denote the best-fit values of the parameters of the intrinsic spectrum under the null and alternative hypotheses, respectively.

Owing to the non-linear impact of ALPs on the photon spectrum, the application of Wilks' theorem [79] is unsuitable here, as discussed in [23]. Therefore, the TS distribution cannot be adequately described by a $ \chi^2 $ distribution, necessitating Monte Carlo simulations to obtain a realistic TS distribution. In this study, we employ the $ {\rm{CL}}_{\rm{s}} $ method [63, 64, 80] to establish constraints on the ALP parameters. The constraints are derived following the same procedure as that described in Refs. [37, 38]. Here, we only offer a concise introduction to this method.

We explore the $ (m_a, g_{a\gamma}) $ parameter space and assess the exclusion of each parameter point using the $ {\rm{CL}}_{\rm{s}} $ method in accordance with observations. For each parameter point, we generate two mock datasets, denoted as $ \{{\rm{d}}\}_{\rm{s}}+{\rm{b}} $ and $ \{{\rm{d}}\}_{\rm{b}} $, based on the expected spectra with and without ALPs, respectively. Both mock datasets $ \{{\rm{d}}\}_{\rm{s}}+{\rm{b}} $ and $ \{{\rm{d}}\}_{\rm{b}} $ consist of 1000 samples. For a specific mock data sample, in each energy bin, the photon flux is randomly generated from a Gaussian distribution, where the mean value and deviation are set as the expected flux and experimental uncertainty, respectively. Utilizing $ \{{\rm{d}}\}_{\rm{s}}+{\rm{b}} $ and $ \{{\rm{d}}\}_{\rm{b}} $, we obtain two TS distributions $ \{{\rm{TS}}\}_{\rm{b}} $ and $ \{{\rm{TS}}\}_{\rm{s+b}} $ using Eq. (14). Given the TS value $ {\rm TS_{obs}} $ obtained from the actual observed data, the CL_s value is defined as

where $ \rm CL_{s+b} $ and $ \rm CL_b $ represent the probabilities of finding a TS value larger than $ {\rm TS_{obs}} $ according to the distributions $ \{{\rm{TS}}\}_{\rm{s+b}} $ and $ \{{\rm{TS}}\}_{\rm{b}} $, respectively. If CL_s is less than 0.05, this parameter point is considered to be excluded at the 95% confidence level (C.L.).

In Fig. 1, we show the TS distributions $ \{{\rm{TS}}\}_{\rm{s+b}} $ and $ \{{\rm{TS}}\}_{\rm{b}} $ for two specific parameter points $(m_a, g_{a\gamma})= (8\times 10^{-7}\, {\rm{eV}}, 3\times 10^{-10}\, {\rm{GeV}}^{-1})$ and $(3\times 10^{-7}\, {\rm{eV}}, 10^{-10}\, {\rm{GeV}}^{-1})$ as examples. The corresponding CL_s values for these two parameter points are 0.0 and 0.16, respectively. These results indicate that the first parameter point can be excluded at the 95% C.L., whereas the second parameter point is still allowed.

Figure 1.  (color online) TS distributions $ \{{\rm{TS}}\}_{\rm{b}} $ and $ \{{\rm{TS}}\}_{\rm{s+b}} $ for two parameter points $(m_a, g_{a\gamma})= (8\times 10^{-7}\; {\rm{eV}}, 3\times 10^{-10}\; {\rm{GeV}}^{-1})$ and $(3\times 10^{-7}\; {\rm{eV}}, 10^{-10}\; {\rm{GeV}}^{-1})$ are shown in the left and right panels, respectively. The vertical black solid lines represent $ \rm TS_{obs} $.

In this section, we present the constraints on the ALP parameters derived from the LHAASO observations of four Galactic sources, namely, the Crab Nebula, LHAASO J2226+6057, LHAASO J1908+0621, and LHAASO J1825-1326. Note that LHAASO J2226+6057, LHAASO J1908+0621, and LHAASO J1825-1326 have multiple potential candidates, as detailed in Extended Data Table 2 in Ref. [52]. The distances to these sources are important for determining the propagation length of photons within the Galactic magnetic field, thereby influencing the constraints on the ALP parameters. For consistency, the distances used in Ref. [52] for spectral fitting are adopted here, which are 0.8 kpc for LHAASO J2226+6057, 3.4 kpc for LHAASO J1908+0621, and 3.1 kpc for LHAASO J1825-1326.

In our analysis, we calculate the best-fit $ \chi^2 $ values for the four sources under the null hypothesis, yielding $ \chi^2/d.o.f. $ = 1.98, where $ d.o.f. $ denotes the number of degrees of freedom. These results indicate that the observed data align well with the null hypothesis without the presence of ALPs. The corresponding best-fit spectra are depicted by the blue lines in Fig. 2. Furthermore, to illustrate the influence of ALPs on the spectrum, we include the best-fit spectra for three ALP parameter points in Fig. 2.

Figure 2.  (color online) Best-fit spectra for the Crab Nebula, LHAASO J1825-1326, LHAASO J1908+0621, and LHAASO J2226+6057. The solid and dashed lines represent the spectra under the null and alternative hypotheses, respectively. The black points denote the photon spectra measured by the LHAASO [52, 53].

Using the CL_s method, we scan the parameter space with $ m_a \in [10^{-8},10^{-5}]\; {\rm{eV}} $ and $ g_{a\gamma} \in [10^{-11},10^{-9}]\; {\rm{GeV}}^{-1} $ and establish constraints at the 95% C.L. for each source, as illustrated in Fig 1. The solid lines in black, purple, red, and blue represent the constraints derived from the LHAASO observations of the Crab Nebula, LHAASO J2226+6057, LHAASO J1908+0621, and LHAASO J1825-1326, respectively. The most stringent constraint of $ g_{a\gamma} $ is approximately $ 1.4 \times 10^{-10} $ $ {\rm{GeV}}^{-1} $ with an ALP mass of $ m_a \sim 4\times 10^{-7} \; {\rm{eV}} $.

Notably, the constraints from the Crab Nebula are considerably more stringent than those from the other sources. This can be attributed to two advantages of the LHAASO observations of the Crab Nebula. First, the spectrum of the Crab Nebula encompasses both WCDA and KM2A results and is precise in the energy regions of $ {\cal{O}}(1) $ TeV. In contrast, the spectra of the other sources only include the KM2A observations at energies above $ {\cal{O}}(10) $ TeV. This indicates that it is easier to precisely determine the intrinsic spectrum of the Crab Nebula and investigate the effects of ALP-photon oscillation in the Crab spectrum at lower energies. Second, the highest energy bin of the Crab spectrum reaches $ 1\; {\rm{PeV}} $, surpassing that of the other sources. As previously mentioned, the compensation of ALP-photon oscillation to the absorption effect may be considerably more significant for higher energy photons. These characteristics enable the Crab Nebula to provide more stringent constraints.

Given that the constraints from the individual sources complement each other in the parameter space, we present the combined analysis result in Fig. 3. The green region represents the combined constraint of these sources. This improves the constraints from a single source for ALP masses above $ 10^{-6}\; {\rm{eV}} $.

Figure 3.  (color online) Constraints at the 95% C.L. in the $ m_a-g_{a\gamma} $ plane resulting from the LHAASO observations of the Crab Nebula, LHAASO J2226+6057, LHAASO J1908+ 0621, and LHAASO J1825-1326. The black dashed line represents the constraints obtained from CAST observations, indicating $g_{a\gamma} < 6.6\times 10^{-11}\; {\rm{GeV}}^{-1}$ [81].

Because the observation of the Crab Nebula provides the most stringent constraints among the four sources, we conduct an analysis combining the observations from the LHAASO and other experiments, including HAWC [65], $ {\rm{AS}}\gamma $ [66], HEGRA [67], MAGIC [68], HESS [69], and VERITAS [70] for the Crab Nebula. $ \chi^2 $ for this analysis is defined as [31]

where the subscript k denotes the k-th experimental data, i denotes the i-th energy bin, and $ \widetilde{\Phi} $, Φ, and $ \delta \Phi $ are the expected value, observed value, and uncertainty of the photon flux, respectively. Because the high-energy gamma-ray experiments under consideration typically have an energy resolution of approximately $ {\cal{O}}(10) $ percent, the spectra measured by different experiments may not precisely match. To derive consistent results for all experiments while accounting for this effect, we introduce additional free parameters f to scale the energies of all the experiments except for the LHAASO and add their Gaussian contributions to the $ \chi^2 $ function. As experimental data is often presented in the form of $E^n \dfrac{{\rm d} N}{{\rm d} E}$, Φ and $ \delta \Phi $ are also scaled by a factor of $ f^{n-1} $, with n=2 in this study. We take the deviations of scale factors $ \delta f $ according to experimental energy resolutions, with values of 0.15 for HEGRA, MAGIC, HESS, and VERITAS, 0.14 for HAWC, and 0.12 for $ \rm{AS \gamma} $. This approach enables us to accommodate the uncertainties arising from energy reconstruction and conduct a more comprehensive analysis incorporating data from multiple experiments.

We achieve a best-fit $ \chi^2/d.o.f. $ of 1.39 under the null hypothesis, with the corresponding spectrum illustrated in Fig. 4. Furthermore, the best-fit spectra for three ALP parameter points and the observational data of the Crab Nebula from various experiments are depicted in Fig. 4. Utilizing the CL_s method, we establish 95% C.L. constraints on the ALP parameters, as depicted by the red solid line in Fig. 5. Notably, the constraints exceed the results obtained solely from the LHAASO observations. For comparison, constraints derived from observations by ASγ, HAWC, HEGRA, and MAGIC, as reported in [31], are also included in Fig. 5. It is evident that our constraints are more stringent than those in Ref. [31], and the most stringent constraint at $ m_a \sim 4 \times 10^{-7}\; {\rm{eV}} $ is similar to the CAST constraint in [81].

Figure 4.  (color online) Best-fit photon spectra for the Crab Nebula. The solid and dashed lines represent the spectra under the null and alternative hypotheses, respectively. The ALP parameter points I, II, and III represent $(m_a, g_{a\gamma})= ({1\times 10^{-6}} $$ {{\rm{eV}}, \; 1\times 10^{-9}\; {\rm{GeV}}^{-1}})$, $(8\times 10^{-6} \;{\rm{eV}},$$\; 1\times 10^{-9} \;{\rm{GeV}}^{-1})$, and $({1\times 10^{-6}} $$ {{\rm{eV}}, \; 5\times 10^{-11}\; {\rm{GeV}}^{-1}})$, respectively. The data points represent observations from the LHAASO [53], HAWC [65], $ {\rm{AS}}{\gamma} $ [66], HEGRA [67], MAGIC [68], HESS [69], and VERITAS [70].

Figure 5.  (color online) Constraint at the 95% C.L. in the $ m_a-g_{a\gamma} $ plane derived from observations of the Crab Nebula, combining the data from seven experiments, that is, the LHAASO [53], HAWC [65], $ {\rm{AS}}{\gamma} $ [66], HEGRA [67], MAGIC [68], HESS [69], and VERITAS [70], is depicted by the red solid line. The blue solid line represents the combined constraint at the 95% C.L. with the inclusion of the new LHAASO KM2A measurement [71] derived in this study. For comparison, the constraint from the observations of ASγ, HAWC, HEGRA, and MAGIC, as reported in Ref. [31], is depicted as the orange dashed line. The shaded region represents the parameter space excluded by CAST [81]. The constraints from the NGC 1275 observation of Fermi-LAT [23] (brown dashed line), the Mrk 421 observation of MAGIC [33] (purple dashed line), the PKS 2155-304 observation of H.E.S.S. [18] (gray dashed line), and the polarization measurements of thermal radiation from magnetic white dwarf stars [43] (green dashed line) are also presented.

The full LHAASO array has been operational since July 2021. Recently, the LHAASO collaboration reported the refined selection criteria of the detector KM2A for very high and ultra-high-energy photons, using data from August 2021 to August 2022 [71]. With these updated data selections, the LHAASO collaboration presented a new spectrum of the Crab Nebula. In comparison to the 2021 spectrum, the highest energy bin of this new spectrum was below 1 PeV. Because only KM2A data is utilized in the analysis, the energy of this new spectrum is presented above 10 TeV. Given the differing datasets and selections of these two analyses, we treat the new LHAASO measurement as independent and integrate it into the combined analysis for the Crab Nebula. The constraint is illustrated in Fig. 5. As shown, this new constraint exhibits a slight improvement over the constraint derived in the previous analysis.

Fig. 5 shows constraints that are more stringent than the CAST constraint, including those derived from the NGC 1275 observation of Fermi-LAT [23], the Mrk 421 observation of MAGIC [33], the PKS 2155-304 observation of H.E.S.S. [18], and the polarization measurements of thermal radiation from magnetic white dwarf stars [43], for comparative purposes. It is expected that future precise measurements of the high-energy gamma-ray spectra of the LHAASO will further improve the constraints from extragalactic sources. It is important to note that although these constraints are more stringent than those derived from the high-energy Galactic sources in our study, they are contingent upon specific magnetic field models and the astrophysical environments of different sources. The constraints derived in our study are independent and hold significance.

The oscillation effect between ALPs and photons is heavily reliant on the strength and structure of astrophysical magnetic fields. As a result, the main theoretical uncertainties arise from the uncertainties in the magnetic field model. In the previous analysis, we employ the Jansson & Farrar model as a benchmark Galactic magnetic field model. To investigate the impact of the Galactic magnetic field model, we conduct analyses using two alternative models, namely, the Pshirkov bisymmetric and Pshirkov axisymmetric models [82]. Our findings depicted in Fig. 6 indicate that these magnetic field models yield consistent results with our previous findings.

Figure 6.  (color online) Constraints derived from observations of the Crab Nebula with various Galactic magnetic field models, including the Jansson & Farrar [77], Pshirkov bisymmetric [82], and Pshirkov axisymmetric models [82].

In this study, we investigate the impact of the ALP-photon oscillation effect on the gamma-ray spectra of four Galactic sources, namely, the Crab Nebula, LHAASO J2226+6057, LHAASO J1908+0621, and LHAASO J1825-1326, measured by the LHAASO. We consider the compensation of ALP-photon oscillation to the absorption effect for high-energy photons and utilize the CL_s method to set constraints on the ALP parameters.

Among the four sources, the Crab Nebula provides considerably more stringent constraints than the other sources owing to its energy spectra spanning a wide range over three orders. By combining the data from the four sources, we find that ALP-photon couplings larger than $ 1.4\times 10^{-10} $ $ {\rm{GeV}}^{-1} $ can be excluded for an ALP mass of $ \sim 4\times 10^{-7} \; {\rm{eV}} $ at the 95% C.L. Furthermore, we perform a combined analysis for the observations of the Crab Nebula from the LHAASO and other experiments. Our analysis sets a limit on $ g_{a\gamma} $ of approximately $ 7.5\times 10^{-11} $ $ {\rm{GeV}}^{-1} $ for an ALP mass of $ \sim 4 \times 10^{-7}\; {\rm{eV}} $, which is similar to the CAST constraint.