Revisiting the analysis of axion-like particles with the Fermi-LAT gamma-ray observation of NGC1275

In this work, we re-analyze the Fermi-LAT observation of NGC 1275 to search for axion-like particle (ALP) effects and constrain ALP parameters. Instead of fitting the observed spectrum with ALP models, we adopt an alternative method for the analysis of this source which calculates the irregularity of the spectrum. With the newly used method, we find no spectral oscillation for the NGC 1275 and rule out couplings $g_{a\gamma}>3\times10^{-12}\,{\rm GeV^{-1}}$ around ALP mass of $m_a\sim$ 1 neV at 95\% confidence level, which is more stringent than the previous results. We also show that the constraints can be further improved by combining the observation of PKS 2155-304. We suggest that with more sources taken into account, we could obtain a much wider exclusion region.


Introduction
Historically, to solve the strong CP problem, Roberto Peccei and Helen Quinn proposed a new gauge symmetry [1]. The spontaneously broken of this symmetry results in a new particle, axion. Axion-like particles (ALPs) have properties similar to those of axion. One major difference between them is that both the mass m a and the coupling strength to photons g aγ of axion are proportional to the energy scale [1,2,3], while they are independent parameters for ALPs [4]. Although ALPs cannot account for the strong CP problem, both axions and ALPs could be good candidates to interpret the mysterious dark matter (DM) existed in our universe [5,6,7,8,9,10,11]. Thanks to the improved instrument sensitivity and abundant observation data accumulated in recent years, astronomical investigations on ALPs can be performed with gamma-ray observations. Considering the gamma-ray observations are generally sensitive only to the ALP parameter space, below we focus mainly on ALPs.
The coupling of ALPs to photons is described by a Lagrangian term of the form L aγ = g aγ E · Ba, where E, B and a represent electric, magnetic, and ALP fields, respectively. This term leads to a property of ALPs that they can convert into photons and vice versa in an external electromagnetic field (Primakoff effect [12]). Such a feature provides people a channel to detect ALPs indirectly by examining the behavior of photon beams after passing through magnetic fields. Via the Primakoff process, many laboratory experiments have been designed to search for signatures of the photon-ALPs conversion, such as ALPS, CAST and ADMX and so on [13,14,15]. Based on the same effect, astrophysical observations could be used to study ALPs as well. Gamma-ray emissions from distant astrophysical sources would pass through a series of magnetic field environments before they reaching the observer on the Earth.
Their spectra will be modulated if ALP fields exist and unique characteristics appear. More explicitly, the conversion between photons and ALPs will lead to oscillation or absorption structures in the observed spectrum, while the signals predicted by a standard astrophysical radiation mechanism usually have continuous spectrum. A wide range of the parameter space of ALP model has been probed according to observations from various target sources [16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39] The photon-ALP oscillation can also account for the problem of low observed opacity of the Universe for TeV photons. Very high energy (VHE) photons travelling in the space will interact with the extragalactic background light (EBL) and transform to electron/positron pairs, which prevent them from propagating a long distance. Therefore, the Universe is opaque to sufficiently distant TeV sources. However, previous researches showed that we can see VHE gamma-ray photons in the regime of optical depth τ > 2 even adopting a minimum EBL model [40,41,42,43]. This issue can be easily addressed if we take into account an ALP effect [44,45,46]. The VHE photons emitted from distant sources may convert into ALPs, which do not interact with EBL and propagate unimpededly in the universe, so that they can travel much a longer distance. These ALPs convert back into gamma-ray photons in the Milky Way's magnetic fields, making them detectable by the observer on the Earth.
Currently, studies on searching for ALPs or constraining the ALP properties using gamma-ray observations have been conducted over a wide variety of astrophysical environments. With the HESS observations of the distant BL Lac object PKS 2155-304 (z=0.116), HESS collaboration quantified the degree of irregularity of the spectrum and constrained the ALP parameters. As they didn't find any signal, more stringent upper limits on the g aγ can be provided for m a between 15 and 60 neV comparing to laboratory experiments [47]. Making use of 6 years of Fermi-LAT observations, Fermi-LAT collaboration searched for irregularities in the gamma-ray spectrum of the radio galaxy NGC 1275. Though no indication of ALP effect has been found, the parameter space of g aγ > 5 × 10 −11 GeV −1 for ALP masses 0.5 − 5 neV was ruled out according to the null result [48]. In addition, analyses of some Milky-Way sources (pulsars and supernova remnants) show intriguing indications of the spectral oscillation [49,50]. However, the best-fit parameters are in tension with the limit from the CAST helioscope and further analysis suggested that the spectrum irregularity in SNR may be due to systemic effects [51] (see however [52] and [53] for possible solution to the tension between CAST limit and the positive hint for ALPs). Meanwhile, the absence of counterpart gamma rays from a Milky-Way supernova explosion may probe the parameter space where ALPs could constitute the entire DM in the Universe [54].
In this paper, we re-analyze the Fermi-LAT observations of NGC 1275 to further investigate the ALP properties. Since previous work is based on six years of Fermi-LAT data [48], the exposure time of the source at present has been doubled. Moreover, in recent years, the gamma-ray emission of this source is in a flaring state with higher fluxes (see Fig. 2). We speculate the higher statistics of the data can provide more stringent constraints on the ALP parameters 1 . In addition, a different analysis method is used in our work. Rather than fit the source's spectrum with both null and ALP models and compare the goodness of fit between the two, we instead calculate the degree of irregularity of the spectrum. The method previously has been adopted in the ALP analyses of HESS data [47] and X-ray observations of galaxy clusters [55].

ALP propagation in magnetic fields
Our calculation of the ALP propagation refers mainly to [22] and the references therein.
The ALP-photon coupling is described by the Lagrangian For an initially polarized photon beam propagating through a single homogeneous magnetic field domain, the propagation equation is where A 1 (x 3 ) and A 2 (x 3 ) are the polarization states along x 1 and x 2 axis respectively, a(x 3 ) denotes the ALP state. M 0 is the photon-ALP mixing matrix The terms ∆ ⊥ , ∆ , ∆ a , ∆ aγ depend on m a , g aγ , photon energy E, strength of the transverse magnetic field B T and electron density n e . The absorption rate Γ is included in the mixing matrix to take into account the EBL absorption. We use the EBL model [41] to calculate the absorption rate.
For an initially unpolarized photon beam, it is convenient to work with the density matrix which obeys the von-Neumann-like commutator equation The final state of the photon-ALP beam is given by the solution to the above equation with T the full transfer matrix which can be determined by solving Eq. (2). The explicit solution can be found in [56]. The photon survival probability in the final state is thus given by where ρ 11 = diag(1, 0, 0), ρ 22 = diag(0, 1, 0) and the initial state ρ(0) = 1/2diag(1, 1, 0). From source to the observer, gamma rays pass through four main magnetic field environments: the fields in the jet (JET), the intracluster magnetic field (ICMF), the intergalactic magnetic field (IGMF), and the Galactic magnetic field of the Milky Way (GMF). For the ALP parameters and the Fermi-LAT energy range considered in this work, the JET and IGMF components are negligible.
The typical field strength of ICMF is in the range of 1 -10 µG. Faraday Rotation Measure (RM) observations and magnetohydrodynamic simulations show that the magnetic field is turbulent, which is usually modeled as a divergence-free homogeneous isotropic field with Gaussian turbulence with zero mean and a variance δ B . The power spectrum follows a power-law M(k) ∝ k q with the wave number k L < k < k H (k L = 2π/Λ max , k H = 2π/Λ min ). The radial profile of the magnetic field strength is B r = B 0 [n e (r)/n e (0)] η . For the profile of the electron density the Eq.(4) in [57] is used. To better compare with the previous work (i.e. Ajello et al. [48]), we use the following parameters to model the ICMF component, r max = 500 kpc, δ B = 10 µG, η = 0.5, Λ min = 0.7 kpc, Λ max = 35 kpc, q = −2.8. The uncertainty to the ALP effect induced by the choice of the values of these parameters have been discussed in [48,46,28].
The GMF consist of a large-scale regular component and a small-scale random component. Since the coherence length of the turbulent component is much smaller than the photon-ALP oscillation length, we do not include it in the calculation. For the regular magnetic field, we take into account the model developed by [58] that has been widely adopted in similar researches [48,54,49,31].

Fermi-LAT data analysis and calculation of the spectrum irregularity
The Fermi-LAT is a pair-conversion gamma-ray telescope sensitive to >30 MeV photons [59]. In this work, we analyze 12 years of Fermi-LAT Pass 8 data towards the direction of NGC 1275, which range from August 4, 2008 to August 4, 2020 (239557417 -616478947 in MET). Considering the relatively bad performance of LAT below 100 MeV and low statistics above 500 GeV, the data are extracted in the energy range of 100 MeV to 500 GeV with SOURCE event class (evclass=128, evtype=3). We adopt an region of interest (ROI) of 15 • centered on the radio position of the target source, NGC 1275. To eliminate the contamination from the Earth Limb, events with zenith angles larger than 90 • are ignored in the analysis (Z max = 90 • ). A data quality cut (DATA QUAL>0)&&(LAT CONFIC==1) is applied to guarantee that the data are suitable for scientific use. The instrument response functions (IRFs) of P8R3 SOURCE V2 is chosen to match the extracted data. To fit the observation data, both the target source and background sources are included in the model. For background sources, we consider all the 4FGL point sources 25 • around the NGC 1275 together with the Galactic and extragalactic diffuse components (described by gll iem v07.fits and iso P8R3 SOURCE V2 v1.txt respectively). The script make4FGLxml.py 2 is employed to generate xml model file.
The likelihood analysis are implemented by the command gtlike. We first perform a likelihood analysis using the whole data set covering all the energy range and time range (we call it global fit) to obtain a good estimation of the background parameters. During the global fit, the spectral parameters of the point sources within 5 • and normalization parameters of all the point sources within ROI and the two diffuse components are left free to vary, and energy dispersion is enabled to achieve better accuracy. Next, we divide the data set into multiple energy and time bins to derive the spectral energy distribution (SED) and light curve (LC) of NGC 1275. For SED we use 100 bins (0.1 -100 GeV) and for LC we use a weekly bin size. In each bin, a likelihood analysis is performed to derive the corresponding flux. All spectral indices are fixed in this procedure. If the fit is subject to a converge problem 3 , for background point sources we will only keep the prefactor of a 0.5 • nearby source free to vary. The resulting SED and LC are demonstrated in Fig. 1 and Fig. 2. As shown in Fig. 1, a LogParabola model can give a good overall estimation of the source spectrum.
2 https://fermi.gsfc.nasa.gov/ssc/data/analysis/user 3 In this case, the resulted error bar of the flux will be extremely small, so As mentioned above, in this work, we do not directly fit the observed spectrum with ALP model. Instead, we calculate the degree of irregularity (I) of the spectrum. The idea is, the more significant the ALP effect is, the higher degree of irregularity will appear in the spectrum. If a set of ALP parameters (m a , g aγ ) predicts a high value of I, but the one derived from actual observation (I obs ) is very low, then the parameters are disfavoured. In our analysis, considering the random nature of the turbulent magnetic fields and the statistical fluctuation, we will simulate 50,000 (see below) energy spectra for each pair of ALP parameters (m a , g aγ ) to derive the expected distribution of I. The distribution then is compared to the I obs to determine whether the ALP parameters are ruled out or not at the 95% confidence level (C.L.). Hence, how to quantify the spectral irregularity is crucial in our analysis.
In [47], the local power-law behavior is tested over three consecutive energy bins (triplet) of the spectrum, the residuals of the middle bins from the power law are quadratically summed to form the irregularity estimator (see [47] for details). However, this method cannot be applied directly to the analysis here. The small bin size of the NGC 1275 SED makes the spectral structure induced by ALP effect (after energy dispersion) is much wider than the bin size. Therefore each triplet behave always as a power-law and the irregularity cannot be evaluated correctly. Alternatively, we propose to define the level of irregularity as The spectrum is fitted in a series of energy windows of (0.5E mid , 1.5E mid ). In each window, a power law function is used to fit the spectrum. The narrow window width ensure the power law that the calculated I value is incorrect. approximation is reasonable. The E mid is the mid energy of each window with increments in steps of 0.5σ E starting from 200 MeV. Similar method called 'sliding window analysis' is often used in line-like signal searches [60,61]. The j and i denote the jth bin in the ith window. The φ pl is the best-fit PL spectrum, and φ, σ are the measured flux and corresponding uncertainty, respectively. The irregularity can also be estimated with the scatter of the fit residuals, as suggested in Ref. [19], where d is the degree of freedom and θ denotes all the spectral parameters of best-fit null model. One can notice that the I defined here is in fact a reduced-χ 2 of the best fit. The N = 100 is the bin number of the SED. A shortage of this method is that an intrinsic overall spectrum φ w/oALP need to be assumed first. However for NGC 1275, as we have shown in Figure 1, a Log-Parabola function can give a good overall estimate of the global spectrum (see also the I value derived in Sec. 4), thus would not introduce too much additional bias. For simplicity, we use Eq. (11) to quantify the level of irregularity in this work. For the upper limit bins (i.e. TS< 25), we ignore them in calculating the I (Eq.(11)), the corresponding bins will also be ignored in the below simulations. ALP parameter space of m a in (0.07 − 100) neV and g aγ in (0.1 − 7) × 10 −11 GeV −1 is divided into a (20 × 20) grid with logarithmic steps. For a given set of ALP parameters (m a , g aγ ), 500 random realizations of the magnetic fields are simulated considering the random nature of the turbulent ICMF. We calculate a photon survival probability P γγ for each of them. The model expected photon counts in each energy bin thus can be calculated by where D kk is LAT energy dispersion matrix and E(E) is the exposure. We use the built-in commands gtexpcube2 and gtdrm 4 of Fermitools software to obtain the E and D, respectively. The integration calculation in Eq.(12) is implemented using composite trapezoidal rule. Each interval of integration ∆E k is divided into 50 sample points (note that the interval ∆E k is already very small). We have tested that subdividing it into 100 points would not significantly improve the precision, since the modulated spectrum is rather smooth in a small energy interval after convolution by the instrument response function (see Fig. 9). To account for the observational effect, Poisson noises are added on µ k and a sample containing 100 pseudo spectra are produced for each B realization and ALP parameter pair. That is, for one set of (m a , g aγ ), 500 × 100 spectra are simulated. Each one of them is fitted with a LogParabola function and the I is derived. The I distribution for one (m a , g aγ ) pair is presented in Fig. 3. Our strategy of restricting ALP parameters is as follows. The null hypothesis is that the observed I is drawn from the probability density function (PDF) of a given set of ALP parameters. Thus, if the observed value is located lower than 5% of the PDF (integrate from the left side), the null can be rejected and the ALP parameters can be excluded with a confidence level of 95%, meaning that the spectral irregularity generated by the ALP effect is too large to produce the actual observation.

Results
Using the Eq. (11), we calculate the spectrum irregularity of NGC 1275 and obtain I obs = 1.07. Such a low irregularity does not indicate existing oscillation structure in the energy spectrum. This value is compared with the I distribution of 50,000 simulations predicted by the model without ALP effect. The resulting distribution is shown in Figure 3 as blue, which is peaked around I ∼ 1. Fitting the distribution with a χ 2 function, we can get the significance corresponding to the I obs = 1.07, being 0.8σ. In order to test this result, we also use Eq. (10) to calculate the irregularity, which also gives a low significance. We therefore could conclude that there is no significant irregularity behavior in the spectrum of NGC 1275. NGC 1275 is a source with dramatic variability (see Figure 2), and its emissions in different periods should have diverse energy spectra. One may imagine that the superposition of the different spectra will lead to a spectral irregularity. However the results here show that it is not the case. This is because the astrophysical radiation processes always produce emission with smooth and continuous spectrum, and their superposition should also be smooth in spectrum. Our result support that it is feasible to search for ALPs through calculating the irregularity even for the variable source. In the future work, we can further search other blazars to see if there are spectral oscillation in some sources. We can also analyze a sample of sources to explore the correlations between the I and other quantities, such as strength of magnetic field, redshift, which could not only gives hints to ALP properties but also the blazar physics.
Since there is no oscillation, we constrain the ALP parameters based on our results. As mentioned above, for each set of ALP parameters (m a , g aγ ), we examine whether its expected I distribution is consistent with the observed one. We exclude the ALP parameters with the distribution deviate from I obs by 95% (68%) probability. The exclusion region is shown in Fig-Figure 4: The excluded parameter space of ALPs based on the analysis of NGC 1275. The solid and dashed lines are the exclusion region at 68% and 95% confidence levels, respectively. The color represents the quantile of the distribution where the I obs is located (start from large I). For a comparison, the constraints from [48] is also plotted (cyan). ure 4. The solid line (dashed line) is for 95% (68%) C.L. As we can see from the figure, our approach can give relatively strong constraints, better than the previous works based on NGC 1275. We have ruled out g aγ > 3 × 10 −12 GeV −1 around ALP mass of m a ∼ 1 neV at 95% C.L. In particular, the 'hole-like' parameter region survived in the [48] can be excluded here 5 .
We need to clarify that during the simulation of giving the I distribution, Poisson fluctuation is only taken into account for the target sources (NGC 1275), but not included for neighboring sources and backgrounds. We test how would this affect our results in the following. To take into account all these fluctuations, we use the Fermipy python package to generate mock counts cubes of the ROI, and then perform likelihood analysis on the mock data. In our analysis framework, we simulate pseudo spectra under signal models, therefore for each (m a , g aγ ) pair we need to perform the whole simulation once, which is very time-expensive. Therefore, we only choose two sets of (m g , g aγ ) to demonstrate the influence. One is at the very center of the exclusion region and the other is close to the boundary.
Based on the best-fit model map of the ROI, we add Poisson fluctuations to it to generate pseudo counts cubes, then likelihood analysis is performed on these pseudo data. For each set of ALP parameters, we simulate 50,000 counts cubes (100 pseudo spectra per realization, totally 500 B realizations), and obtain the corresponding SEDs. The subsequent steps are the same as the standard method. The obtained I distributions are present in Fig. 5. As shown in the figure, the I values by average are smaller than those from the standard method, but the difference is not great. In particular, for the set of parameters close to the boundary of the previous exclusion region (left panel of Fig. 5), the new distribution is very close to the standard one. This may be due to that the NGC1275 is a very bright gamma-ray source, so that the influence by the background fluctuation is subdomi- nant. Therefore, we conclude that the differences between two methods would not significantly affect the results.

Joint constraints with PKS 2155-304
We also analyze the source PKS 2155-304, which has been widely used to study the ALP properties. We analyze 12-years Fermi-LAT data of PKS 2155-304 (4FGL J2158.8-3013) and constrain the ALP parameters with the same procedure as described in Sec.3. For the ICMF parameters, we adopt the values in [28] (σ B = 10 µG, Λ min = 0.5 kpc, Λ max = 20 kpc, q = −11/3, η = 0.5). By calculating the irregularity of the spectrum, we obtain I = 1.2, corresponding to a significance of 1.5 σ. The constraints of the ALP parameter based on this source are shown in the left panel of Fig. 6. The results are consistent with that presented in [28].
One of the advantages of the approach used in this article is that we can easily combine multiple sources together to improve the constraints. Here we show the example of combining NGC 1275 and PKS 2155-304. The combined irregularity is defined as I obs = χ 2 i / d i . For any given ALP parameters (m a , g aγ ) or a w/oALP model, the distribution of I com can be easily derived by the distribution of I NGC1275 and I PKS2155−304 generated before. As a result, by comparing the observed I com with the model expected distribution of I com , we derive the combined constraints (see right panel of Fig. 6). We can see that, by combining these two sources together, the constraints are further improved, better than the results reported before [48,53]. We suggest that with more sources taken into account it will hopefully to further improve the results.

Projected constraints on ALP parameters with improved energy resolution
We would like to discuss how the energy resolution of the instrument would alter the results here. The irregular structure caused by the ALP effect would vary rapidly with energy, especially when the turbulent magnetic fields are taken into account [47]. Poor energy resolution will eliminate these structures, leading to a loss of the sensitivity. A high energy resolution may thus be crucial in the ALP research. At present, in the range of GeV-TeV energy band, the DArk Matter Particle Explorer (DAMPE) [62] has the best energy resolution, which can achieve <1.5% above ∼ 1 GeV. However, the effective area (a eff ) of DAMPE is small. We try to use the a eff and energy dispersion of DAMPE to perform a simulation, we find that no ALP parameters in the parameter region considered in this paper can be constrained. Therefore, we further consider the following case. We assume that in the future, an instrument can observe the same amount of photons of NGC 1275 as Fermi-LAT in the energy range of 1 GeV-1 TeV, but has an energy resolution of 1.5%. We examine what ALP parameters can be probed by such a instrument. The simulation is similar to that in Sec. 3, except we adopting an energy dispersion of a Gaussian function with σ E = 0.015E true , where E true is the true energy of the incident photon. We assume that no spectral irregularity is observed, namely I obs ∼ 1. The results are presented in Figure 7. In the region of m a < 1 neV, the constraints are weaker than above because we only consider > 1 GeV data. It can be seen that the exclusion region is not enlarged too much although the energy resolution has been improved significantly. Similar results has also been demonstrated in [31]. This results indicate that the effect of energy resolution on ALP search is not as strong as we imagined.

Constrains with regular magnetic field
The ALP-photon mixing depends largely on the morphology and strength of the magnetic field in which the γ-ray emitter resides. Following Ref. [48], a purely turbulent ICM is considered when deriving the above results. Ref. [33] pointed out that the use of a purely turbulent magnetic field in the Perseus cluster is not supported by any Faraday Rotation Measure observations. Therefore, they adopted a regular magnitude field to derive the  constraints on ALP-photon coupling. Such a regular field is associated with the X-ray cavity near NGC 1275 and is capable of reproducing its RM observation of 6500 ∼ 7500 rad m −2 [63]. It is found that the constraints are relaxed by orders of magnitude due to the reduction of spectral irregularity induced by ALP (please see [33] for details).
Here we test how would the limits change when using the regular magnetic field. Following the magnetic field configuration described in [64] (namely the one adopted in [33]), we derive the result presented in Fig. 8. As is shown, while the exclusion region is not the same as the turbulent one, a large amount of ALP parameters can be constrained. Furthermore, around m a ∼ 1.0 neV, the 90% upper limits of g aγ is down to ∼ 10 −12 GeV −1 , close to the boundary that ALP could account for all dark matter [48].
The result of Fig. 8 is obviously inconsistent with those re- Figure 8: The excluded parameter space of ALPs based on the regular magnetic field of Perseus [33]. The solid and dashed lines are the exclusion region at 68% and 95% confidence levels, respectively. The color represents the quantile of the distribution where the I obs is located (start from large I). For a comparison, the constraints from [48] is also plotted (cyan). The red point corresponds to the parameters we use to examine the results (see text for details).
ported in [33]. To examine whether the excluded parameters are really disfavored by the observation, for a representative pair of ALP parameters (magenta point in Fig. 8) we compare the model expected ALP spectrum with observation in Fig. 9. We find that although the regular magnetic field would not introduce dramatic oscillations in the spectrum, there is a rapid decrease of P γγ around the critical energy (∼GeV) followed by a constant P γγ ∼ 0.6, making the modulated spectrum deviate from the observed LogParabola-like one. For many ALP parameters, the flux varies at a level of 30%, greater than the error bars of the Fermi-LAT measurements. As a result, many ALP parameters can still be constrained.
Since an alternative method has been used in our analysis, one may wonder if this is the reason for the inconsistency. For the selected pair of parameters, we also use the χ 2 fit to compare models. We fit the observed spectrum with both ALP and w/o-ALP model and derive their minimum χ 2 . The best fitted models are shown in the right panel of Fig. 9. We obtain a ∆χ 2 ∼ 327. Such a large ∆χ 2 we believe can exclude the ALP model, although a null distribution based on Monte Carlo simulation is required to give precise judgement.
We note that recently Meyer et al. released their code gam-maALPs for ALP calculation online 6 . To check the correctness of our program and to ensure the discrepancy is not due to errors in our calculation, we also use their code to derive corresponding results. Not strangely, since our calculation procedure mainly refers to Meyer et al. [22] (and the references therein), we find that our results consistent with theirs very well. Therefore, we tend to think that the inconsistency is caused by other reasons (e.g. choice of model parameters), which need to be further investigated. The code for our calculation of the ALPphoton conversion is also put online now 7 .
We have confirmed that the inconsistency between the two works is mainly due to the longer data set used in our analysis. In Ref. [33], they used 6 years of Fermi-LAT data with EDISP3 event type, while we are using 12-year data without the event-type selecion (evtype=3). Together with the fact that NGC 1275 is in a flaring state in recent years (see fig), therefore the data amount in Ref. [33] is ∼ 1/20 of ours. The non-linear dependence of the ALP effect on the ALP parameters further makes the constraints different for ∼2 orders of magnitude.

Summary
ALP-photon oscillation in external magnetic fields will induce spectral irregularity of high-energy gamma-ray sources.
The degree of irregularity, I, can be quantified with Eq. (11). A low observed I obs in the spectrum of astrophysical source disfavor the ALP parameters predicting high values of irregularities. In this work, we revisit 12 years of Fermi-LAT observation towards NGC 1275, a source which has been widely considered in ALP studies. We calculate the spectral irregularity and then compare it to Monte Carlo simulations to test whether the corresponding ALP parameters are excluded. We find NGC 1275 has a low degree of irregularity of I obs = 1.07, indicating no significant evidence for ALPs. We thus can place constraints on ALP parameters more stringent than previous works and rule out g aγ > 3 × 10 −12 GeV −1 around ALP mass of m a ∼ 1 neV at 95% C.L. The "holelike" feature in the parameter region not probed in [48] has also been excluded. We show that the constraints can be further improved by combining the observation of PKS 2155-304, reaching g aγ ∼ 2 × 10 −12 GeV −1 (assuming δ B = 10 µG for PKS 2155-304). It is thus worthwhile to combine more other sources in future work to enlarge the parameter region. We also find that a better instrument energy resolution would not enhance the constraints too much.
It has been pointed out that a purely turbulent magnetic field is difficult to reconcile with the rotation measure of Perseus cluster, a regular magnetic field component should be included in the ALP analysis [33]. We therefore also derive constraints on g aγ with the purely regular field. We however obtain different results comparing to the Ref. [33]. We find that although the regular magnetic field would not introduce dramatic irregularity in the spectrum, a prominent step around E c , which is induced by the transition to the strong mixing regime, makes the modulated spectrum not consistent with the observed LogParabola-like one. As a result, many ALP parameters can still be constrained, even down to g aγ ∼ 2×10 −11 GeV −1 around m a ∼ 1 neV.