Probing photon-ALP oscillations from the flat spectrum radio quasar 4C+21.35

Flat spectrum radio quasar (FSRQ) is the most luminous blazar at the GeV energies. In this paper, we probe the photon-axion-like particle (ALP) oscillation effect on the latest very-high-energy (VHE) $\gamma$-ray observations of the FSRQ 4C+21.35 (PKS 1222+216). The $\gamma$-ray spectra are measured by the collaborations Major Atmospheric Gamma Imaging Cherenkov Telescopes (MAGIC), Very Energetic Radiation Imaging Telescope Array System (VERITAS), and Fermi Large Area Telescope (Fermi-LAT), which cover two activity VHE flares of 4C+21.35 in 2010 and 2014. We show the spectral energy distributions (SEDs) of these two phases under the null and ALP hypotheses, and set the combined limit on the ALP parameter space. The 95% $\rm C.L.$ combined limit set by the FSRQ 4C+21.35 observations measured by MAGIC, VERITAS, and Fermi-LAT in the $m_a-g_{a\gamma}$ plane is roughly at the photon-ALP coupling $g_{a\gamma} \gtrsim 8\times 10^{-12} \rm \, GeV^{-1}$ for the ALP mass $[\,2\times 10^{-10}\, {\rm eV} \lesssim m_a \lesssim 2\times 10^{-8}\, \rm eV\,]$. Compared with the constraint of NGC 1275 set by Fermi-LAT, no stringent limit result is derived with the photon-ALP coupling $g_{a\gamma}$ from the FSRQ 4C+21.35, while this result could slightly broaden the ALP mass $m_a$ limit at the low-mass region.


I. INTRODUCTION
Axions are originally proposed in the Peccei-Quinn solution to solve the strong CP problem in quantum chromodynamics (QCD) [1][2][3][4]. Several extensions of the Standard Model (SM) at the high energies, such as string theories [5,6], predict the axion-like particles (ALPs). ALPs are ultralight pseudo-Nambu-Goldstone bosons (pNGBs) with the two-photon vertex g aγ . ALPs are also potential dark matter (DM) candidates [7][8][9]. The relevant interaction of photon-ALP in the external magnetic field can be described by the effective Lagrangian [10] L ALP = 1 2 ∂ µ a∂ µ a − 1 2 m 2 a a 2 − 1 4 g aγ aF µνF µν , where − 1 4 g aγ aF µνF µν = g aγ aE · B, with the ALP field a, the ALP mass m a , the coupling constant between the ALP and photons g aγ , the (dual) electromagnetic field tensor (F µν ) F µν , and the local electric and magnetic field vectors E, B. This photon-ALP coupling in the astrophysical magnetic field environments would lead to many potentially observable photon-ALP oscillation effects [11][12][13][14].
In this work, we focus our attention on the photon-ALP oscillation effect on the very-high-energy (VHE) γray photons emitted by the astrophysical sources which are far from the Milky Way, such as blazars [11,15]. Blazars are radio-loud active galactic nuclei (AGNs) with the relativistic jets pointing close to our line of sight in the extragalactic VHE sky [16,17]. According to the behaviors of their optical spectra, blazars can be divided into two classes, BL Lacertae (BL Lac) objects and the * lihaijun@bnu.edu.cn flat spectrum radio quasars (FSRQs) [18,19]. The bolometric luminosity of the FSRQ is significantly greater than that of the BL Lac blazar. In this case, we consider the photon-ALP oscillation in the blazar source region magnetic field and then further back-conversion in the Galactic magnetic field, which provides an effective way to reduce the extragalactic background light (EBL) absorption of the VHE γ-ray photon at the energy above 100 GeV [20,21]. Many studies  are performed to probe the effects of the photon-ALP oscillation on the VHE blazar spectra modifications or set the ALP limits in the m a − g aγ plane.
be ineffective to broaden the ALP parameter g aγ at the ALP mass [ 1 × 10 −10 eV < m a < 1 × 10 −7 eV ] region with the BL Lac blazar observations, which also depends on the astrophysical magnetic field models and the uncertainties of the VHE γ-ray observations. Additionally, only a few literatures discuss the photon-ALP oscillation effect from the FSRQ blazars [24,29,64,65]. Therefore, we turn our attention to the FSRQ object and expect to derive the more stringent constraint with the photon-ALP coupling g aγ or set the ALP limit at the different ALP mass region. Recently, the collaborations Very Energetic Radiation Imaging Telescope Array System (VERITAS) and Fermi -LAT reported the VHE γ-ray flare observations of the FSRQ blazar 4C+21.35 in 2014 [66]. The FSRQ 4C+21.35 (also known as PKS 1222+216), located at the redshift of z 0 = 0.432 1 , was discovered at the VHE region in a flaring state by MAGIC in 2010 [67]. Above two major γ-ray flares in June 2010 and February-March 2014 of 4C+21.35 are denoted as the phases Flare 2010 and Flare 2014, respectively. In this paper, we investigate the photon-ALP oscillation effect on the VHE γ-ray observations of the FSRQ 4C+21.35 measured by MAGIC, VERITAS, and Fermi -LAT with these two phases. We 1 http://tevcat2.uchicago.edu/sources/dDaGC0 also set the 95% C.L. combined limit on the ALP parameter (m a , g aγ ) space.
This paper is organized as follows. In Sec. II, we describe the astrophysical environments model for the VHE γ-ray photons propagating from the FSRQ blazar source region to the Milky Way. Our numerical results are presented in Sec. III. In Sec. IV, we comment on our results and conclude. The statistic method is given in Appendix A.

II. MODEL FOR THE FSRQ 4C+21.35
In this section, we introduce the magnetic field environments model for the FSRQ 4C+21.35 with the photon-ALP beam propagating from the source region to our Earth. Generally, the propagation process can be divided into the blazar source region, the extragalactic space, and the Milky Way region [11,15].
The main aim of this section is to derive the whole transport matrix T (s) of the photon-ALP beam in the above propagation process, which can be described by where T (s i ) region−i is the transport matrix of a single i-th region with the propagation distance s i . Then we can derive the final survival probability P γγ of the VHE γ-ray on the Earth [23,28] with where ρ(0) is the initial density matrix of the photon-ALP beam, and ρ(s) = T (s)ρ(0)T † (s) is the final density matrix.

A. Blazar source region
We begin with the photon-ALP propagation process in the FSRQ blazar source region, which can also be divided into three parts, the broad line region (BLR), the blazar jet region, and the host galaxy region. Compared with the BL Lac object, FSRQ is in a sense a more complicated version of the blazar source region.

Broad line region
For the photon-ALP beam propagates in the BLR of 4C+21.35, we consider the photon-ALP oscillation effect in the BLR magnetic field (BLRMF) and the BLR photon absorption effect due to the pair-production process [24] with the background photon γ BG . The optical depth of the BLR photon absorption in this region can be described by [68] τ BLR = 2π dµ dω dr(1 − µ)n(ω, µ, r) where µ = cos θ, θ is the scattering angle between the VHE photon E and the background photon ω, r is the distance from the BLR center, n(ω, µ, r) is the spectral number density of the BLR, and σ γγ (E, ω, µ) is the γγ pair-production cross-section [69] where σ T is the Thompson cross-section, β is the dimensionless parameter This cross-section implies that the absorption effect would become maximal for the photon energy In our analysis, the τ BLR distribution of the FSRQ 4C+21.35 is taken from Ref. [24]. Following Ref. [24], we adopt the homogeneous transverse magnetic field model for the BLRMF of 4C+21. 35 with the benchmark value B BLRMF = 0.14 G, which is approximately consistent with the results obtained in Refs. [70][71][72]. The radius of the 4C+21.35 BLR is adopted as R BLR = 0.23 pc.

Blazar jet region
Similar to the BL Lac object, we consider the photon-ALP oscillation effect in the blazar jet magnetic field (BJMF) for the FSRQ 4C+21.35 blazar jet region. Generally, this magnetic field can be modeled as the poloidal (B ∝ r −2 ) and toroidal (B ∝ r −1 ) components [73]. Following Refs. [33,41], we adopt the jet transverse magnetic field model B BJMF (r) as [74,75] and the electron density distribution model n BJMF el (r) as [76] n BJMF where B BJMF 0 and n BJMF 0 are the core magnetic field strength and the electron density at R BLR , respectively. Following Ref. [24], here we take B BJMF 0 = B BLRMF = 0.14 G, n BJMF 0 = 10 2 cm −3 , η BJMF = −1, and ξ BJMF = −2. For the distance r > 6.7 kpc, we take the magnetic field B BJMF = 0. Moreover, we introduce the photon energy transformation between the laboratory frame E L and the co-moving frame E j where δ D is the Doppler factor with the benchmark value 20 for the blazar jet region of 4C+21.35 [71,72,77].

Host galaxy region
FSRQ is hosted by the elliptical galaxy, the magnetic field of this region is rarely known. As discussed in Ref. [24], the photon-ALP oscillation effect in this part plays a very minor role and can be totally neglected. Therefore, for the photon-ALP beam propagates in the host galaxy region of 4C+21.35, we do not consider the photon-ALP oscillation effect in the magnetic field.
Additionally, for the blazars located in the rich environment cluster, one should consider the photon-ALP oscillation effect in the turbulent inter-cluster magnetic field (ICMF) [28,78]. This magnetic field model B ICMF (r) can be described by with the electron density distribution model n ICMF el (r) where B However, there is no definite evidence that the FSRQ 4C+21.35 considered in this work is located in the rich environment cluster. Thus, we do not take into account the photon-ALP oscillation effect in this ICMF.

B. Extragalactic space
For the photon-ALP beam propagates in the extragalactic space, we also consider the EBL photon absorption effect due to the pair-production process. The optical depth of the EBL photon absorption in this region can be described by [79,80] with where z 0 is the redshift of the source, E th is the threshold energy,σ(E, ω, z) is the integral pair-production crosssection, dn(z)/dω is the proper number density of the EBL, H 0 67.4 km s −1 Mpc −1 is the Hubble constant, Ω m 0.315 is the matter density, and Ω Λ 0.685 is the dark energy density [81]. In this work, we take the EBL model from Ref. [80].
Moreover, we do not consider the photon-ALP oscillation effect in the extragalactic magnetic field. The strength of this magnetic field is very weak (current limit: 10 −7 nG B ext 1.7 nG [82,83]) and can not be explicitly measured.

C. Milky Way
For the photon-ALP beam propagates in the Milky Way, we also consider the photon-ALP oscillation effect in the Galactic magnetic field. Following Refs. [84,85], the Galactic magnetic field can be modeled as the disk and halo components (parallel to the Galactic plane), and the "X-field" component (out-of-plane at the Galactic center). One can find the latest version about this magnetic field model in Ref. [86].

III. RESULTS
With the magnetic field environments model for the FSRQ 4C+21.35 considered in Sec. II, we can derive the final survival probability of the VHE γ-ray on the Earth, which is shown in Fig. 2. The solid and dot-dashed lines represent the photon survival probability with/without the photon-ALP oscillation effect, respectively. We show the final photon survival probability distributions for five typical ALP parameter sets (with m neV =1, 10, 100 and g 11 = 1, 10, 100) in the m a − g aγ plane. The behavior of the photon-ALP oscillation varies dramatically with the ALP mass m a and the coupling constant g aγ . Using the final photon survival probability P γγ in Eq. (3), we can derive the expected VHE γ-ray spectrum where Φ int (E) is the γ-ray intrinsic spectrum. Following Ref. [50], the VHE blazar intrinsic spectrum model can be described by a simple and smooth concave function with several free parameters. In Ref. [43], we explored the effects between the γ-ray intrinsic spectrum model selections and the ALP limits, which showed an insignificant relationship between them. For the phase Flare 2010 of 4C+21.35, we use the power law with exponential cut-off (EPWL; three free parameters) model while for the phase Flare 2014, we use the power law with super-exponential cut-off (SEPWL; four free parameters) model with the normalization constant N 0 , the normalization energy E 0 = 1 GeV, the spectral index Γ , and the free parameters E c and d.
Then the χ 2 value can be described by with the spectral point number N , the expected γ-ray spectrum Φ exp (E i ), the detected γ-ray spectrumφ i , and the corresponding uncertainty of the observation δ i . The best-fit γ-ray spectral energy distributions (SEDs) of the two phases Flares 2010 and 2014 of 4C+21.35 under the null hypothesis are shown in Fig. 3.
For each ALP parameter set in the m a − g aγ plane, we can derive the value of the best-fit χ 2 ALP . We also give the minimum best-fit γ-ray SEDs for the two phases of 4C+21.35 under the ALP hypothesis in Fig. 3 for comparisons. Following Refs. [41,42], we show the result for these two phases combined together, which is essential for the multi-phase analysis. The χ 2 ALP distribution in the m a − g aγ plane for the two phases Flares 2010 and 2014 of 4C+21.35 combined is shown in Fig. 4.
Following Refs. [30,41], in order to set the 95% C.L. limit on the ALP parameter space, 400 sets of the γ- We show the 95% C.L. photon-ALP combined limit set by the FSRQ 4C+21.35 observations measured by MAGIC, VERITAS, and Fermi -LAT in Fig. 1 (one can also find in Fig. 5) with the photon-ALP coupling at g aγ 8×10 −12 GeV −1 for the ALP mass [ 2×10 −10 eV m a 2 × 10 −8 eV ]. The other latest photon-ALP limits at the ALP mass [ 1 × 10 −13 eV < m a < 1 × 10 −6 eV ] region are also shown for comparisons. Compared with the constraint of NGC 1275 set by Fermi -LAT [30], we could not obtain a more stringent limit result with the photon-ALP coupling g aγ from the FSRQ 4C+21.35, while this result could slightly broaden the ALP mass m a limit at the mass m a 2 × 10 −10 eV region.
As discussed in Ref. [41], the magnetic field environments model parameters B BLRMF (B BJMF 0 ) and R BLR would significantly affect the final limit results and can not be completely determined in this work. For comparison, we give the ALP limit results with the large values of these two model parameters. In Fig. 5, the 95% C.L. photon-ALP combined limits set by the FSRQ 4C+21.35 observations with the parameters B BLRMF (B BJMF 0 ) = 0.20 G (χ 2 ALP distribution is shown in Fig. 4) and R BLR = 0.30 pc are given, which show the slightly more stringent limits.

IV. CONCLUSION
In this paper, we have presented the photon-ALP oscillation effect on the latest VHE γ-ray observations of the FSRQ 4C+21.35 measured by VERITAS and Fermi - LAT in 2014. Combined with the 4C+21.35 observations measured by MAGIC and Fermi -LAT in 2010, the 95% C.L. limit on the ALP parameter space is roughly at the photon-ALP coupling g aγ 8 × 10 −12 GeV −1 for the ALP mass [ 2×10 −10 eV m a 2×10 −8 eV ]. Compared with the constraint of NGC 1275 set by Fermi -LAT, no stringent limit result is obtained with the photon-ALP coupling g aγ from the FSRQ 4C+21. 35. However, this result shows a new photon-ALP excluded region at the low-mass (m a 2 × 10 −10 eV) region. Additionally, we also test the effects of the magnetic field environments model parameters B BLRMF (B BJMF 0 ) and R BLR on the final ALP limits.
The results in this paper are of great importance to the next generation VHE γ-ray measurements, like the Large High Altitude Air Shower Observatory (LHAASO) [87], Tunka Advanced Instrument for Gamma-ray and Cosmic ray Astrophysics-Hundred Square km Cosmic ORigin Explorer (TAIGA-HiSCORE) [88], High Energy cosmic-Radiation Detection (HERD) [89], Cherenkov Telescope Array (CTA) [90], and Gamma-Astronomy Multifunction Modules Apparatus (GAMMA 400) [91], which would also provide more accurate data for the FSRQ blazars to set the ALP limits.

ACKNOWLEDGMENTS
The author would like to thank Ciaran A.J. O'Hare for his public webpage https://github.com/cajohare/ AxionLimits [52] from which a number of photon-ALP limits are given. This work is supported by the National Natural Science Foundation (NNSF) of China (Grants No. 11775025 and No. 12175027).