Exploring axion-like particle from observation of FSRQ Ton 599 by Fermi-LAT

Axion-like particles (ALPs), a class of extremely light pseudoscalar bosons, are anticipated in various extensions of the Standard Model [1−3]. Unlike the axions associated with the solution to the strong CP problem within the standard model [4−7], ALPs offer a broader parameter space that is yet to be fully investigated. The effective coupling between ALPs and photons can induce ALP-photon oscillations in the presence of external magnetic fields. Given the prevalence of astrophysical magnetic fields, the phenomenon of ALP-photon oscillation has attracted significant attention in the field of astrophysics [8−43]. The oscillation between ALPs and photons may lead to irregularities in high energy γ-ray spectra. Various astrophysical sources, including blazars [26, 28, 33, 35, 44], gamma-ray bursts [36, 45−48], galaxy clusters [49−51], supernova remnants [31, 52], and pulsars [53−55], have been utilized to investigate the properties of ALPs. Notably, extragalactic sources, from which high energy photons traverse magnetic field environments on larger spatial scales compared to Galactic sources, hold great promise for investigating the effects of ALPs.

Among extragalactic sources, blazars have garnered significant attention for their prominence in the extragalactic γ-ray sky. Blazars are a subclass of active Galactic nuclei (AGNs) characterized by the presence of relativistic jets typically oriented toward Earth. They are further classified into two categories: BL Lacertae (BL Lac) objects and flat-spectrum radio quasars (FSRQs), based on the rest-frame equivalent width of the emission lines observed in their optical spectra. FSRQs show strong broad emission lines with equivalent widths exceeding 5 $ \mathring{A} $, whereas BL Lacs exhibit either absent or weak emission lines. The presence of these broad emission lines suggests that FSRQs contain rapidly moving gas clouds near the central black hole, known as the broad-line region (BLR). Furthermore, infrared observations of FSRQs indicate the presence of a dusty torus (DT) located beyond the BLR. In FSRQs, both the BLR and DT play critical roles in reprocessing photons from the accretion disk, resulting in the emission of low-energy radiation. The BLR primarily emits a high density of ultraviolet photons, whereas the DT generates infrared photons. High-energy photons emitted from the emission region can interact with these ultraviolet and infrared photons, potentially being absorbed through the $ e^+e^- $ pair production process. This intricate astrophysical environment of FSRQs surpasses that of BL Lac objects, prompting previous investigations to primarily focus on BL Lac sources.

Despite the increased complexity of FSRQs compared to BL Lacs, they can serve as important targets for investigating ALP-photon oscillations. The detection of very high energy photons from FSRQs may provide evidence of ALP-photon oscillations [56, 57], offering a potential solution to avoiding the substantial absorption effects induced by the BLR and DT. The analysis of the γ-ray spectra of FSRQs can also be used to set constraints on the ALP parameters [20, 58−61]. The rate of ALP-photon oscillations is significantly influenced by the intensity of external magnetic fields. The magnetic field strength within the emission region of FSRQs is typically in the range of $ {\cal{O}}(1) $ to $ {\cal{O}}(10) $ G, whereas that of many BL Lacs typically ranges from $ {\cal{O}}(0.1) $ to $ {\cal{O}}(1) $ G [20, 62, 63]. In the Third LAT AGN Catalogs, the median magnetic field strength values of FSRQs are nearly an order of magnitude larger than those of BL Lacs, as shown in Ref. [64]. Therefore, FSRQs have the potential to exhibit a significant ALP-photon oscillation rate.

In this study, we consider the observations of FSRQ Ton 599. This source is positioned at Right Ascension (RA) = 179.88$ ^{\circ} $ and Declination (Dec) = 29.25$ ^{\circ} $, with a redshift of z = 0.725. Ton 599 has been observed by various instruments including EGRET [65], Fermi-LAT [66], and VERITAS [67]. Ton 599 presents significant variability in both optical and γ-ray energy bands, with observations of this source typically classified into quiescent and bursting phases. Recently, in Ref. [68], a detailed analysis of data from five distinct periods, encompassing both quiescent and bursting phases, was conducted using observations from Fermi-LAT. The multi-wavelength analysis indicates a magnetic field strength in the emission region on the order of Gauss. In this work, we set constraints on the ALP parameters based on these results of Ton 599.

This paper is organized as follows. In Section II, we introduce the ALP-photon oscillation effect. In Section III, we introduce the astronomical environment of Ton 599 and calculate the survival probability of photons from the FSRQ. In Section IV, we describe the process of fitting the gamma-ray spectra and the statistical method. In Section V, we present the constraints on the ALP parameters from the Fermi-LAT observations of FSRQ Ton 599 and combined constraint from the observations of all epochs. Finally, we present the conlusions in Section VI.

The two key parameters in the effective theory influencing the ALP-photon oscillation are the ALP mass $ m_a $ and ALP-photon coupling $ g_{a\gamma} $. The state of the ALP-photon system can be characterized by the density matrix $ \rho \equiv \Psi \otimes \Psi^\dagger $, where $ {\Psi}\equiv {(A_\perp,\; A_\parallel,\; a)}^T $. Here, a represents the ALP field, and $ A_\perp $ and $ A_\parallel $ denote the photon polarization amplitudes perpendicular and parallel to the transverse component of the external magnetic field $ B_t $, respectively. When the system with energy $ E \gg m_a $ traverses a homogeneous magnetic field, the density matrix ρ obeys a von Neumann-like commutator equation [10, 14], expressed as follows:

where z denotes the distance along the propagation direction, and $ {\cal{M}}_0 $ is the mixing matrix encompassing various ALP and electromagnetic effects. $ {\cal{M}}_0 $ is given by

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

High energy photons, emitted from extragalactic sources, pass through various astrophysical magnetic fields on their way to Earth. The entire path can be separated into multiple segments, with the magnetic field in each segment assumed to be constant. By solving Eq. (1), the survival probability of the photon can be expressed as [8, 13]:

where $ {\cal{T}}(\boldsymbol{z}) \equiv \prod \limits_i^n {\cal{T}}_i (\boldsymbol{z}) $, and $ {\cal{T}}_i(\boldsymbol{z}) $ is the transferring matrix obtained from the i-th segment. The polarization of very high energy γ-rays is typically unmeasurable; therefore, the γ-ray photons emitted from the source are assumed to be unpolarized, and $ \rho(0) $ is taken to be $ \text{diag}(1/2,1/2,0) $ in this case.

In this section, we discuss the various astrophysical environments that influence photons originating from the source Ton 599, considering absorption effects and ALP-photon oscillations. These environments include the BLR, DT, blazar jet, extragalactic space, and Galactic region. Specifically, we consider ALP-photon oscillations in the blazar jet and Galactic regions, absorption effects arising from background photons in the extragalactic region, and background photons emitted from the BLR and DT.

The BLR is a high-velocity gas structure surrounding the central black hole in the AGN. The gas within the BLR emits broad spectral lines and typically exhibits rapid rotational motion around the black hole. The DT is a large-scale toroidal structure composed of dust located in the equatorial plane of the AGN. Photons originating from the accretion disk interact with both the BLR and DT. The BLR and DT reprocess these photons, leading to the emissions of ultraviolet and infrared light, respectively. These ultraviolet and infrared photons subsequently interact with high-energy photons emitted from the central emission zone, absorbing them and thereby altering the observed spectrum of γ-ray photons.

In this study, we characterize the BLR as an infinitesimally thin spherical shell located at a distance $ R_{\rm BLR} $ from the central black hole. The BLR reprocesses the disk radiation $ L_{\rm disk} $ with a fraction denoted by $ \xi_{\rm BLR} $, which is assumed to be 0.1. The emission from the BLR occurs at a single energy corresponding to the Mg II emission line. Moreover, the DT is modeled as a ring-shaped structure located at a distance $ R_{\rm DT} $ from the central black hole. The DT reprocesses the disk radiation with a fraction of $ \xi_{\rm DT} = 0.5 $, and emits radiation at a single energy determined by its temperature $ T_{\rm DT} $. The parameter values for the BLR and DT are adopted from Ref. [68], as outlined in Table 1. The radius of the emission region and distance between the emission region and central black hole, denoted as R and $ r_{\rm VHE} $, are also listed. Notably, the value of $ R_{\text{BLR}} $ is not provided in Ref. [68]. Therefore, we have estimated an average distance based on the inner and outer radii of the BLR $ R_{\rm BLR_{\text{in}}} $ and $ R_{\rm BLR_{\text{out}}} $ given in Ref. [68]. In our analysis, we utilize the open-source package agnpy [69] to calculate the absorption of γ-rays by ultraviolet and infrared photons emitted from the BLR and the DT.

Name of parameters	Values
$ {L}_{\rm disk} $	$ 4.5 \times 10^{45} \; {\rm erg/s} $
$ \xi_{\rm BLR} $	$ 0.1 $
$ R_{\rm BLR} $	$ 3.89 \times 10^{-5} \; {\rm kpc} $
$ \xi_{\rm DT} $	$ 0.5 $
$ R_{\rm DT} $	$ 1.62 \times 10^{-3} \; {\rm kpc} $
$ T_{\rm DT} $	$ 1000 \; {\rm K} $
$ {R} $	$ 2.59 \times 10^{-6} \; {\rm kpc} $
$ r_{\rm VHE} $	$ 2.59 \times 10^{-4} \; {\rm kpc} $

Table 1.  Parameters of Ton 599 sourced from [68].

For FSRQs such as Ton 599, the photons emitted from the vicinity of the central black hole traverse several astrophysical magnetic fields before reaching Earth. One of the initial regions encountered is the blazar jet, where ALP-photon oscillations may occur within the blazar jet magnetic field (BJMF). The BJMF can be typically described as having both toroidal and poloidal components. At larger distances from the central black hole, the toroidal component dominates [70, 71], while the poloidal component becomes negligible. Therefore, in this study, we focus exclusively on the toroidal component of the magnetic field. The strength of the magnetic field in the BJMF and electron density can be characterized as follows [70−72]:

where $ r_{\rm VHE} $ is the distance between the emission region and central black hole, and $ B_{0} $ and $ n_{0} $ represent the magnetic field strength and electron density at the emission region, respectively. Furthermore, the magnetic field is assumed to diminish beyond the maximum scale of the jet, which is taken to be 1 kpc.

In this study, we consider the observations of Ton 599 spanning five epochs within the period from MJD 54686 to MJD 60008 [68]. Each epoch encompasses a 100-day duration. Among these five epochs, one is categorized as a quiescent epoch, denoted as epoch A, whereas the remaining four are flaring epochs, denoted as epochs B, C, D, and E. During epochs B, C, and E, gamma-ray flares were concurrently observed with optical flares, whereas in epoch D, the corresponding gamma-ray flares were either weak or absent. The parameter values of $ B_0 $ and $ n_0 $ for five epochs of Ton 599 obtained from multi-wavelength analyses in Ref. [68] are listed in Table 2.

Parameters	EpochA	EpochB	EpochC	EpochD	EpochE
$ n_{0}/(\times 10^3 {\rm cm^{-3}}) $	0.98	1.19	1.31	0.89	0.84
$ B_{0} /({\rm Gauss}) $	1.63	1.51	1.51	2.05	1.73

Table 2.  Values of parameters $ B_{0} $ and $ n_{0} $ during five epochs. $ B_{0} $ and $ n_{0} $ are the magnetic field strength and electron density at the emission site, respectively.

High energy photons are expected to be absorbed by the extragalactic background light (EBL) through the $ e^+e^- $ pair production process. We adopt the results presented in [73] to model the EBL.

The extragalactic magnetic field remains undetectable, with only an upper limit of approximately $ \; {\cal{O}}(1) $ nG. If the strength of the extragalactic magnetic field is significantly below this value, the ALP effect can indeed be neglected. Previous studies, such as Refs. [74] and [75], have indicated that for a magnetic field of $ {\cal{O}}(1) $ nG, the ALP effect may significantly affect the EBL attenuation at high energies. For Ton 599, the photon survival probability in the presence of an extragalactic magnetic field of $ 1\; \text{nG} $ would exceed that without such a field by a maximum relative enhancement of $ \sim30 $% at energies above $ \sim 10 $ GeV, for some parameter points situated near the boundary of the excluded region. However, this increase in the photon flux is relatively modest in comparison to the uncertainties associated with the data at high energies for Ton 599. Consequently, the potential influence of an extragalactic magnetic field is not considered in our analysis.

The ALP-photon oscillation may occur within the Galactic magnetic field. The Galactic magnetic field model adopted for this investigation is the Jansson and Farrar model, as detailed in [76]. Additionally, we employ the NE2001 model [77] to characterize the distribution of the Galactic electron density. Notably, the Galactic magnetic field comprises both regular and turbulent components. Given the negligible impact of the turbulent component, we focus exclusively on the regular component.

After incorporating the aforementioned effects, the photon spectrum at Earth can be expressed as follows:

where $ {\rm d}\Phi_{\rm int}/{\rm d}E $ represents the intrinsic spectrum, and $ P_{\gamma\gamma} $ denotes the survival probability of the photon. The survival probability $ P_{\gamma\gamma} $ is calculated numerically using Eqs. (1) and (3). To perform these calculations, we utilize the open-source package gammaALP [78].

In Fig. 1, we illustrate the survival probability of photons in epoch E across various scenarios. The green and purple dash-dotted lines correspond to the cases without and with the absorption in Ton 599 under the null hypothesis, respectively, in the absence of ALP effects. The absorption effect resulting from BLR and DT can significantly suppress the $ \gamma $-ray flux above $ \sim 40 $ GeV.

Figure 1.  (color online) Survival probability of photons from Ton 599 during epoch E. The solid lines represent the results of four specific ALP parameter points including the BLR and DT contributions. The corresponding parameter values of these parameter points are detailed in Table 3. The dashed-dotted lines represent the results without ALP effects: one line includes contributions from the BLR and DT, whereas the other line excludes these contributions.

We display the survival probability of photons for four ALP parameter points, labeled as I, II, III and IV, which are excluded in the subsequent analysis. The corresponding parameter values for these points are listed in Table 3. The conversion between the photon and ALP induces oscillatory behavior in the γ-ray spectrum below $ \sim 10 $ GeV for parameter points I and IV with $ m_a\sim 10^{-9} $ eV. For parameter point II with a smaller ALP mass, the survival probability is almost energy-independent. Conversely, for parameter point III with a large ALP mass, the survival probability oscillates rapidly while almost maintaining an averaged value across large energy bins. These patterns will be helpful in our understanding of the constraint results in Section V.

Parameters	I	II	III	IV
$ m_a /\text{eV} $	$ 2\times10^{-9} $	$ 2\times10^{-10} $	$ 2\times10^{-8} $	$ 2\times10^{-9} $
$ g_{a\gamma}/\text{GeV}^{-1} $	$ 5\times10^{-11} $	$ 5\times10^{-11} $	$ 5\times10^{-11} $	$ 5\times10^{-12} $

Table 3.  Parameter values of four specific ALP parameter points considered in Fig. 1. These parameter points are excluded in the subsequent analysis.

Our analysis indicates that the photons from Ton599 with energies above $ \sim40 $ GeV undergo substantial absorption by the BLR and DT. However, the MAGIC collaboration has reported the detection of photons with energies surpassing 100 GeV from Ton 599 [79]. This observed high energy γ-ray flux amounts to approximately $ 1.5 \times 10^{-10} \;\rm{ph/cm^2/s} $, equivalent to 0.3 Crab units. Note that this result is derived from the one-hour data collected by MAGIC on December 15, 2017, which may correspond to a period of high state for Ton 599. In contrast, our calculations in this analysis are based on observations from different epochs [68], each spanning a 100-day duration. Given that the uncertainties associated with the spectral data points provided in Ref. [68] are substantial above 10 GeV, the precise cutoff of the γ-ray spectrum induced by absorption effects would not significantly affect the constraints on the ALP parameters, resulting from the oscillation effects at lower energies.

In this section, we introduce the analysis method employed to set constraints on the ALP parameters. In this study, we consider the intrinsic spectrum of Ton 599 to follow a log parabola (LP) model given by

where $ N_0 $, α, and β are free parameters, with $ E_0 $ set to100 $ \text{MeV} $.

The determination of the best-fit spectrum involves minimizing the $ \chi^2 $ function, defined as

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

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

For a given set of ALP parameters $ m_a $ and $ g_{a\gamma} $, we define the test statistic (TS) as follows:

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

To set constraints on the ALP parameters, the distribution of the TS must be understood. In cases where ALP parameters have a non-linear impact on the photon spectrum, Wilks' theorem is not applicable, and the TS distribution cannot be described by a $ \chi^2 $ distribution. Therefore, Monte Carlo simulations are required to derive constraints on ALP parameters. In this study, we employ the $ {{\rm CL}_s} $ method [80−82] to establish constraints on these parameters. This method has been used in our previous works [35, 36, 52, 83], and a detailed description of the method can be found in Refs. [35, 36].

The $ {\rm CL}_s $ method can be briefly described as follows. To test whether a specific ALP parameter point $ (m_a, g_{a\gamma}) $ can be excluded, we first generate a set of mock data based on the expected spectrum without ALPs, denoted as $ \{d\}_{b} $. Next, we generate a second dataset using a similar method, but based on the energy spectrum with ALPs. The photon flux for each energy bin is drawn from a Gaussian distribution, with the expected flux as the mean and experimental uncertainty as the standard deviation. We then calculate two TS distributions for the specific $ (m_a, g_{a\gamma}) $ point: one from the dataset $ \{d\}_{s+b} $, which includes the signal and background, and one from $ \{d\}_b $, which contains only the background. These distributions are labeled as $ \{ {\rm TS}\}_b $ and $ \{{\rm TS}\}_{s+b} $, respectively. The observed TS value, denoted as $ {{\rm TS}_{obs}} $, is calculated from the actual observed data. The $ {{\rm CL}_s} $ value is defined as follows:

where $ {\rm CL}_{s+b} $ and $ {\rm CL}_b $ represent the probabilities of obtaining a TS value greater than $ {\rm TS}_{obs} $ according to the distributions $ \{{\rm TS}\}_{s+b} $ and $ \{{\rm TS}\}_b $, respectively. If $ {\rm CL}_s $ is less than 0.05, the corresponding parameter point is considered excluded at a 95% confidence level (C.L.).

In this section, we present the constraints on the ALP parameters derived from the observational results of the Fermi-LAT for the FSRQ Ton 599. Our analysis indicates that the observational data are consistent with the null hypothesis, assuming the absence of ALPs. The best-fit spectra under the null hypothesis are shown as blue dashed lines in Fig. 2, with five subfigures corresponding to epochs A, B, C, D, and E. For comparison, we also display the best-fit spectra under the ALP hypothesis for the four parameter points, for which the survival probability of photons have been shown in Fig. 1.

Figure 2.  (color online) Best-fit spectra for the five observations of Ton 599. The blue dashed lines and red solid lines represent the best-fit spectra under the null hypothesis, with and without considering the absorption effects of the BLR and DT, respectively. The yellow and green dashed lines correspond to the best-fit spectra under the ALP hypothesis, for four selected parameter points, which are the same as those in Fig. 1.

We perform a scan of the parameter space with $ m_a \in [10^{-10}, 10^{-6}]\; \text{eV} $ and $ g_{a\gamma} \in [10^{-13}, 10^{-10}]\; \text{GeV}^{-1} $ and establish constraints at the 95% C.L. using the $ \rm CL_s $ method. The constraints derived from observations of the five epochs are shown in Fig. 3. Owing to the large flux uncertainties in the quiescent epoch A, it cannot provide effective constraints on the ALPs. The most stringent constraints are obtained from epochs B and E. For $ m_a = 8\times 10^{-10}\; \text{eV} $, the constraints on $ g_{a\gamma} $ approach $ 1.3 \times 10^{-12} \text{GeV}^{-1} $ in epoch B. For $ m_a = 9.2\times10^{-10}\; \text{eV} $, the constraints on $ g_{a\gamma} $ approach $ 1.3 \times 10^{-12}\; \text{GeV}^{-1} $ in epoch E. Furthermore, we present the combined constraints obtained from all observations in the five epochs in Fig. 3. These combined constraints are more stringent than those derived from any single epoch, with the best constraint on $ g_{a\gamma} $ reaching approximately $ 1.1 \times 10^{-12}\; \text{GeV}^{-1} $ for $ m_a \sim 9 \times 10^{-10}\; \text{eV} $.

Figure 3.  (color online) TS map in the $ m_a-g_{a\gamma} $ plane based on the individual and combined analyses of five epochs. The solid black lines represent the 95% C.L. constraints established using the $ \rm CL_s $ method in this study. The dash black line denotes the constraints from the CAST experiment [84].

In our analysis, the constraints are effectively imposed on the ALP mass in the order of $ 10^{-9} $ eV. This behavior can be understood as follows. Because the intrinsic spectra parameters are unknown and are treated as free parameters in the analysis, the overall attenuation of the spectra caused by the ALPs does not manifest distinguishable effects in the observations; hence, it does not yield constraints. The constraints arise from the oscillatory behaviors across various energies in the observed spectrum below approximately 10 GeV, which are induced by the significant energy dependent ALP-photon oscillations. As shown in Fig. 1, the survival probability of photons exhibits significant oscillatory patterns for the two ALP parameter points with $ m_a\sim10^{-9} $ eV.

The survival probability of photons is determined by the mixing matrix, as given in Eq. (2), and the propagation distance L. The conversion probability between the photon and ALP within a constant magnetic field can be expressed as follows:

where θ represents the mixing angle given by

and $ \Delta_{\rm{osc}} $ represents the oscillation wave number given by

In the scenario considered for photons with energies below $ \sim {\cal{O}}(1) $ GeV, the oscillatory behaviors in the spectrum across different energies can be triggered by the mass term $ \Delta_a= -m_a^2/2E $ within $ \Delta_{\rm{osc}} $.

When $ m_a \sim {\cal{O}}(10^{ -9}) $ eV in the considered scenario, the contribution to the oscillation rate from the mass term is comparable to the mixing term $ \Delta_{a\gamma} = g_{a\gamma} B/2 $ in the jet, leading to pronounced energy-dependent oscillations. As the ALP masses decrease, the energy-dependent impact of $ \Delta_a $ diminishes, resulting in a primarily global attenuation in the spectrum, which does not yield significant constraints. However, for larger ALP masses, the phase term $ \Delta_{\rm{osc}} L $ would be very large, leading to extremely rapid oscillatory patterns that may not be discerned in observations due to limited resolutions. Moreover, the substantial mass term can suppress the contribution from the mixing term $ \Delta_{a\gamma} $, as shown in Eq. (12), resulting in no constraints in the very high ALP mass region. These qualitative discussions are consistent with the results shown in Fig. 1, providing an understanding of the distribution of excluded regions in the parameter space.

In Fig. 4, we provide a comparison of the results obtained in this study with those from other experimental studies. The shaded purple region represents the comprehensive constraint at the 95% C.L. derived from the results of this study. The shaded blue region represents the parameter space excluded by the CAST experiment [84]. Additionally, constraints from various sources are included: the Fermi-LAT observation of NGC 1275 [23] (brown dashed line), ARGO-YBJ and Fermi-LAT observations of Mrk 421 [33] (purple dashed line), H.E.S.S. observation of PKS 2155-304 [18] (blue dashed line), observation of SN1987A [43] based on the Solar Maximum Mission $ \gamma $-ray data (yellow dashed line), and polarization measurements of thermal radiation from magnetic white dwarf stars [41] (green dashed line). The constraint of FSRQ 4c+21.35 from the observations of MAGIC, VERITAS, and Fermi-LAT, as reported in Ref. [60], is depicted by the orange dashed line. Our study imposes stricter constraints compared to other works within the range $ m_a \in [4\times10^{-10} ,2\times10^{-9}]\; \text{eV} $.

Figure 4.  (color online) Comparison of the constraints derived in this study with those from other studies. More constraints can be found in Ref. [85].

In this study, we investigate the impact of ALP-photon oscillations on the γ-ray observations of the FSRQ Ton 599. By analyzing the Fermi-LAT observations of Ton 599 in 2023, we establish constraints on the ALP parameter space at a 95% confidence level, corresponding to a photon-ALP coupling of approximately $ 1.1 \times 10^{-12} $ $ {\rm GeV}^{-1} $, applicable to ALP masses within the range of $ [4\times10^{-10}-2\times10^{-9}]\; \text{eV} $.

Our study comprehensively considers the internal and external astrophysical environments of the FSRQ Ton 599, including the BLR, DT, blazar jet, extragalactic region, and Galactic region. We determined that the BLR and DT have minor impacts on the results of spectral energy distribution fitting, leading to the conclusion that they do not influence the final outcomes significantly.

The propagation of high energy photons from FSRQ blazars through astrophysical magnetic fields on larger spatial scales provides a valuable opportunity to study ALP-photon oscillations. With the advancement of observational techniques and development of high-precision scientific instruments, research on FSRQ blazars and ALPs will continue to progress. Future large-scale facilities such as LHAASO, MAGIC, HESS, CTA, and DAMPE will enhance the precision of measurements and contribute to further advancements in our research efforts.