The Axion-Photon Mixing and the Extragalactic Magnetic Background: Plateau Regimes, Resonances, and Non-Gaussian Boosts

The Axion-Photon Mixing and the Extragalactic Magnetic Bac kground: Plateau Regimes, Resonances, and Non-Gaussian Bo osts Andrea A ddazi, 1, 2 Yi-F u Cai, 3, 4 Salv atore Capozziello, 5, 6, 7 Qingyu Gan, 6, 7 and Gaetano Lam biase 8, 9 1 Center for The or etic al Physics, Col le ge of Physics, Sichuan University, China 2 INFN, L ab or atori Nazionali di F r asc ati, Via E. F ermi 54, I-00044 R oma, Italy 3 Dep artment of Astr onomy, Scho ol of Physic al Scienc es, University of Scienc e and T e chnolo gy of China, Hefei, Anhui 230026, China 4 CAS Key L ab or atory for R ese ar ches in Galaxies and Cosmolo gy, Scho ol of Astr onomy and Sp ac e Science, University of Scienc e and T e chnolo gy of China, Hefei, Anhui 230026, China 5 Dip artimento di Fisic a ”E. Pancini”, Università di Nap oli “F e derico II”, Compl. Univ. di Monte S. A ngelo, Edificio G, Via Cinthia, I-80126, Nap oli, Italy 6 Scuola Sup erior e Meridionale, Via Mezzo c annone 4, I-80134, Nap oli, Italy 7 Istituto Nazionale di Fisic a Nucle ar e, Sezione di Nap oli, Nap oli, Italy 8 Dip artimento di Fisic a “E.R. Caianiel lo”, Università di Salerno, I-84084 Fisciano (Sa), Italy 9 INFN- Gruppo Col le gato di Salerno, I-84084 Fisciano (Sa), Italy W e presen t an analytical treatmen t of Axion-Lik e-P article (ALP)—photon mixing with extra- galactic background light (EBL) attenuation for constant, Gaussian-sto chastic, and non-Gaussian magnetic field configurations—with direct implications for V ery High Energy (VHE) gamma-ra y observ ations suc h as LHAASO, HA W C, and CT A exp erimen ts. F or constant fields, w e deriv e exact probabilities and identify a perturbative plateau regime where photon surviv al scales as quartic order of magnetic field, isolating the four-point magnetic correlation as a sensitiv e prob e of non- Gaussianit y . F or Gaussian sto chastic fields, we obtain—for the first time—analytical formulas for non-exp onen tial-decay comp onents in the strong-attenuation regime. Con trary to the widely used domain-lik e mo del, photon surviv al is suppressed by 4-6 orders of magnitude, while b oth conv er- sion and surviv al probabilities exhibit distinct m ulti-peak structures from mass-equal resonance, sto c hastic resonance, and EBL attenuation. Extending to non-Gaussian fields, we show that non- Gaussianit y can enhance photon surviv al by several orders of magnitude relative to the Gaussian case, potentially explaining the unexp ectedly VHE photon even t observ ed b y LHAASO. Our results demonstrate that sto c hastic magnetic fields cannot b e reduced to domain-like coherence without losing essen tial physics, and that VHE gamma-ra y spectra enco de observ able information about b oth the p ow er sp ectrum and non-Gaussian structure of in tergalactic magnetic fields—critical as next-generation observ atories push to ward Pe V sensitivities. I. INTR ODUCTION Recen tly the Large High Altitude Air Show er Observ atory (LHAASO) rep orted the detection of more than 140 photons with energies abov e 3 T e V, and 9 photons ab o v e 10 T e V with high statistical significance from Gamma Ra y Burst 221009A at redshift 0 . 15 [ 1 , 2 ]. In addition, the Carp et exp eriment has rep orted a photon-like even t with an energy of appro ximately 300 T e V [ 3 ]. Within the framework of conv entional physics, very-high-energy photons abov e the T e V scale originating from distan t sources suc h as gamma-ray bursts are exp ected to b e strongly atten uated due to pair pro duction with extragalactic background ligh t (EBL), i.e., γ VHE γ EBL → e + e − , making their observ ation c hallenging. One prop osed new-physics explanation in v olv es axion-like particles (ALPs), whic h can oscillate in to photons in the presence of a magnetic field, thereby reducing the attenuation of gamma-ra y flux by the EBL [ 4 , 5 ]. Consequen tly , ALP-photon mixing in the presence of the EBL has been extensively in v estigated across a range of astroph ysical and theoretical con texts. These studies encompass mo del-building efforts [ 6 – 9 ], observ ational phenomena in A GNs and blazars [ 10 – 13 ], the interpretation of X-ray excesses [ 14 ], scenarios inv olving dark photons [ 15 ], photon p olarization effects [ 16 , 17 ], and com bined analyses with neutrino emissions [ 18 , 19 ]. In the studies of ALP-photon mixing—both with and without the influence of EBL—the treatmen t of magnetic fields constitutes a critical asp ect. Metho dologically , existing approaches can b e broadly divided into discretizations and in tegral methods. The discretization approach is typically implemented via domain-like mo dels, wherein the magnetic field is assumed to b e constant within eac h domain, allo wing for an exact solution to the ALP-photon mixing equations [ 6 , 8 – 10 , 13 , 14 ]. In con trast, the in tegral approac h is based on perturbative treatments that in tegrate ov er spatially v arying magnetic field configurations [ 20 – 23 ]. On the other hand, ALP-photon mixing can b e categorized based on the nature of the magnetic field configuration. The first type inv olves deterministic configurations, whic h are typically asso ciated with sp ecific astrophysical environmen ts such as the magnetosphere of neutron stars [ 24 ] or the regular large-scale magnetic fields in galaxy clusters and galaxies [ 25 , 26 ]. The second type concerns sto c hastic magnetic field configurations, whic h are more prev alen t in cosmic settings, including turbulent fields in 2 clusters, galaxies, and jets, as well as primordial magnetic fields and the magnetized cosmic w eb. T o characterize sto c hasticit y statistically , correlation functions are commonly employ ed. T he phenomenon of ALP-photon mixing in sto c hastic magnetic fields has b een studied extensively using sev eral forms of tw o-point correlation functions. These include mono chromatic and scale-in v arian t sp ectra [ 23 ], turbulen t sp ectra [ 10 , 13 , 21 , 27 , 28 ], as w ell as more general sp ectral shap es [ 22 ]. Concerning higher order correlation functions, effects of the magnetic field with Gaussian and non-Gaussian distributions on ALP-photon mixing hav e b een explored in Refs. [ 20 , 29 , 30 ]. In this work, we revisit ALP-photon mixing in the presence of EBL attenuation for b oth constant and sto c hastic magnetic field configurations. F or the constant field case, w e derive the exact solution which generalizes previous results in Refs. [ 6 , 7 , 17 ]. F or a sto chastic field with a Gaussian distribution, we adopt the integral approach and, for the first time, obtain the analytical form ula of the non-exp onen tial-decay term in the ALP-photon con v ersion and photon surviv al probabilities in the strong EBL region. W e clarify the p erturbative conditions for b oth constan t and sto c hastic magnetic fields and find that they are satisfied throughout most of the physically relev an t parameter space. In terestingly , the double oscillation pro cess γ → a → γ implies that, within the p erturbative regime, the photon surviv al probability is dominated b y terms at quartic order in the magnetic field strength, i.e., O ( B 4 ) . This provides a natural framework for inv estigating the impact of non-Gaussian magnetic fields on ALP-photon mixing. The pap er is organized as follows. In Section I I , we solve the ALP-photon mixing equations in a constant magnetic field. Section II I extends this analysis to Gaussian sto c hastic magnetic fields, where we apply the solution to realistic ph ysical scenarios and compare the predictions of the domain-lik e mo del with those of the integral approach. In Section IV , we inv estigate the photon surviv al probability in the presence of non-Gaussian magnetic fields. Finally , w e present a discussion of our findings in Section V and conclude in Section VI . I I. CONST ANT MA GNETIC FIELD An axion-like particle (ALP) field a , characterized b y mass m a and tw o-photon coupling g aγ , can oscillate into a photon (i.e., the electromagnetic p oten tial A µ ) in the presence of a transverse magnetic background field. F or relativistic ALPs with frequency ω satisfying m a ≪ ω , this oscillation giv es rise to an ALP-photon mixing system whose dynamics are gov erned b y the follo wing equations of motion [ 9 , 31 ]: ∂ l   A x A y a   = iK   A x A y a   = − i   ω + ∆ xx ∆ R 0 ∆ R ω + ∆ y y ∆ B 0 ∆ B ω + ∆ a     A x A y a   , (I I.1) where l the propagating direction. Giv en a homogeneous transv erse magnetic field, one ma y rotate the co ordinate frame so that the field aligns with the y -axis. In this frame, the relev ant mixing term is given by ∆ B = g aγ B / 2 . The ∆ a ≡ ∆ a ( ω ) = − m 2 a / 2 ω term corresp onds to the effective ALP mass con tribution. Regarding the photon terms, firstly , the F araday rotation term ∆ R couples the t wo photon p olarizations. While it is relev an t for the analysis of p olarized photon sources and irrelev an t to this study , w e therefore set ∆ R = 0 . Secondly , ∆ xx and ∆ y y include b oth real and imaginary parts, which w e assume to b e equal for simplicity and denote as ∆ xx = ∆ y y = ∆ γ = A − ib . The real part A con tains conv en tional contributions suc h as plasma effect, QED corrections and CMB effects (see Eq. I I I.22 ). F or v ery high energy photons propagating through the in tergalactic medium, it undergo es pair pro duction interacting with the optical and infrared EBL photons, i.e. γ VHE γ EBL → e + e − [ 32 ]. This leads to an EBL atten uation effect in the ALP-photon mixing, which can b e describ ed b y adding an imaginary part b = 1 / (2 λ γ ) with the photon mean free path λ γ ( ω ) . In this work, we adopt the semi-analytical expression used in [ 8 , 30 ], λ γ = 1 . 2 × 10 3  ω T e V  − 1 . 55 Mp c . (I I.2) The front factor 1 . 2 is c hosen suc h that λ γ = Mp c at ω = 100 T e V, serving as a practical b enc hmark p oint.F or later con venience, we define the oscillation length of ALP-photon mixing as l osc = 2 / p A 2 + ∆ 2 B in constant magnetic field and l osc = 2 / | A | in sto chastic magnetic field. The dynamics of ALP-photon mixing can b e solved in terms of the ev olution op erator U ( l ) . The evolution op erator is actually a 3 × 3 matrix defined as  A ⊥ ( l ) , A ∥ ( l ) , a ( l )  T = U ( l )  A ⊥ (0) , A ∥ (0) , a (0)  T , where we set l 0 = 0 as the initial point. The EoM for U ( l ) reads as ∂ l U ( l ) = iK U ( l ) , (I I.3) 3 with the initial condition U (0) = I 3 × 3 . Because K is l -indep endent, the explicit closed solution is simply U ( l ) = e iK l , and non-v anishing comp onents are giv en by U 11 = e − i ( ω +∆ γ ) l , U 22 = 1 2 s e − 1 2 i (2 ω +∆ γ +∆ a ) l  (∆ a − ∆ γ + s ) e 1 2 isl − (∆ a − ∆ γ − s ) e − 1 2 isl  , U 23 = U 32 = − ∆ B s e − 1 2 i (2 ω +∆ γ +∆ a ) l  e 1 2 isl − e − 1 2 isl  , U 33 = 1 2 s e − 1 2 i (2 ω +∆ γ +∆ a ) l  (∆ a − ∆ γ + s ) e − 1 2 isl − (∆ a − ∆ γ − s ) e 1 2 isl  , (I I.4) where s = q 4∆ 2 B + (∆ γ − ∆ a ) 2 . All other comp onents of U v anish. The relev ant probabilities, namely ALP-photon con version probabilit y , photon-ALP conv ersion probability , photon surviv al probability and ALP surviv al probability , are related to U by (see App endix A) P a → γ = U 13 U ∗ 13 + U 23 U ∗ 23 , P γ → a = 1 2 ( U 31 U ∗ 31 + U 32 U ∗ 32 ) , P γ → γ = 1 2 ( U 11 U ∗ 11 + U 12 U ∗ 12 + U 21 U ∗ 21 + U 22 U ∗ 22 ) , P a → a = U 33 U ∗ 33 . (I I.5) The square ro ot s can be ev aluated as s = ( α + iβ , I ≥ 0 , α − iβ , I < 0 , (I I.6) with α = p √ R 2 + I 2 + R √ 2 , β = p √ R 2 + I 2 − R √ 2 , R = 4∆ 2 B + A 2 − b 2 , I = − 2 bA. (I I.7) Com bining Eqs. ( I I.4 ), ( I I.5 ) and ( II.6 ), we obtain the relev an t probabilities after trav eling a distance d , P γ → a ( d ) = 1 2 P a → γ = 1 2 ∆ 2 B α 2 + β 2 e − bd { 2 cosh ( dβ ) − 2 cos ( dα ) } , P γ → γ ( d ) = 1 2 e − 2 bd + 1 4( α 2 + β 2 ) e − bd  A 2 + b 2 + α 2 + β 2  cosh ( dβ ) − 2 ( α | A | + β b ) sinh ( dβ ) −  A 2 + b 2 − α 2 − β 2  cos ( dα ) + 2 ( β | A | − α b ) sin ( dα )  , P a → a ( d ) = 1 2( α 2 + β 2 ) e − bd  A 2 + b 2 + α 2 + β 2  cosh ( dβ ) + 2 ( α | A | + β b ) sinh ( dβ ) −  A 2 + b 2 − α 2 − β 2  cos ( dα ) − 2 ( β | A | − α b ) sin ( dα )  . (I I.8) These expressions are v alid for b oth I ≥ 0 and I < 0 . In the limit b → 0 , they reduce to the con ven tional ALP-photon mixing probabilities in the absence of EBL effects. In the limit A → 0 , they reco v er the results of Refs. [ 7 , 17 ]. Using Eq. ( I I.8 ), w e compute the probabilities for sev eral representativ e parameter sets, with the results shown in Fig. 1 . F or certain parameter choices, b oth the ALP-photon conv ersion and photon surviv al probabilities remain appro ximately constant in the regime d ≫ λ γ , whic h corresp onds to the p erturbative regime of the mixing. T o clarify this b eha vior, we expand Eq. ( I I.8 ) in p ow ers of g aγ B λ γ and fo cus on the large-distance limit d ≫ λ γ , where exp onen tial-deca y terms are strongly suppressed. Then we obtain (see App endix B) P γ ↔ a ∼ ( g aγ B λ γ ) 2 1 1 + A 2 λ 2 γ + ( g aγ B λ γ ) 4 d λ γ 1 1 + A 4 λ 4 γ + ( g aγ B λ γ ) 6 d 2 λ 2 γ 1 1 + A 6 λ 6 γ + ..., P γ → γ ∼ ( g aγ B λ γ ) 4 1 1 + A 4 λ 4 γ + ( g aγ B λ γ ) 6 d λ γ 1 1 + A 6 λ 6 γ + ( g aγ B λ γ ) 8 d 2 λ 2 γ 1 1 + A 8 λ 8 γ + ..., (I I.9) 4 FIG. 1: The ALP-photon conv ersion probabilit y P a → γ (upp er panels) and the photon surviv al probabilit y P γ → γ (lo wer panels) in a constant magnetic field for v arious parameter sets. All dimensional quantities are normalized to λ γ . In the regime d ≪ λ γ , the probabilities reproduce the conv en tional ALP-photon mixing results in the absence of EBL attenuation. F or λ γ ≪ d ≪  1 + A 2 λ 2 γ  /  ∆ 2 B λ γ  , b oth probabilities lie within the p erturbative regime and remain approximately constant, exhibiting p ow er-la w scaling with the magnetic field as P a → γ ∼ B 2 and P γ → γ ∼ B 4 . A t larger distances, the probabilities decay exp onentially due to EBL absorption. Neglecting higher-order terms requires the p erturbativ e condition ( g aγ B λ γ ) 2 d λ γ 1 1 + A 2 λ 2 γ ≪ 1 . (I I.10) When this condition is satisfied, b oth the leading-order term of P γ ↔ a (at order B 2 ) and that of P γ → γ (at order B 4 ) b ecome indep endent of the distance d . Indeed, as shown in Fig. ( 1 ), for example with ∆ B λ γ = 0 . 1 , the probabilities P γ ↔ a and P γ → γ b egin to decay exponentially at d/λ γ > 100 for Aλ γ = 0 . 1 , and they remain constan t P γ ↔ a ≃ 10 − 4 and P γ → γ ≃ 10 − 8 at 1 ≪ d/λ γ < 1000 for Aλ γ = 10 . A similar distinct plateau profile of photon surviv al probability has b een observed in Refs. [ 7 , 10 , 17 ]. As clearly shown in Eq. ( I I.9 ) and Fig. ( 1 ), within the p erturbative region and at d ≫ λ γ , there exists a relation P γ → γ ∼ P 2 γ ↔ a , which can b e interpreted as a double conv ersion γ → a → γ . Ho wev er, no simply relation exists b etw een P γ ↔ a and P γ → γ b ey ond the p erturbativ e regime. Moreov er, we recall that in the absence (or for negligible) EBL effects, the perturbative regime of ALP-photon mixing corresponds to w eak mixing, c haracterized by the condition g aγ B l osc ≪ 1 [ 31 ]. In con trast, in the presence of strong EBL, the relev an t p erturbativ e parameter b ecomes g aγ B λ γ ≪ 1 , as justified by the expansion in Eq. ( I I.9 ) and the plateau b ehavior observ ed in Fig. 1 . In general, the magnetic field is not coheren t across the entire propagation path. T o account for this, a domain-like approac h has b een developed to mo del the decoherence of the magnetic field. In this framework, the total space is partitioned into N domains, each with a size equal to the t ypical correlation length λ B of the v arying magnetic field. Within each domain, the magnetic field v ector is assumed to ha ve a fixed but randomly orien ted direction. F or 5 simplicit y , w e adopt a reduced v ersion of the domain-like mo del in whic h only the y -comp onen t of the magnetic field is retained, implemented via vector pro jection. Specifically , we replace Eq. ( II.1 ) with ∆ B → ∆ B cos ϑ , where ϑ is a random v ariable uniformly distributed within [0 , π ] . Let F (0) γ and F (0) a denote the initial photon and ALP fluxes in the first domain. The iterative relation of photon and ALP fluxes b etw een patch n and patch n + 1 is F ( n +1) γ = F ( n ) γ P ( n ) γ → γ + F ( n ) a P ( n ) a → γ , F ( n +1) a = F ( n ) γ P ( n ) γ → a + F ( n ) a P ( n ) a → a , (I I.11) where P ( n ) denote the probability in n -th patch. F or an initial state consisting purely of ALP (photon), the ov erall ALP-photon conv ersion probability (photon surviv al probabilit y) after tra v ersing N patc hes is giv en by P cell a → γ = F ( N ) γ /F (0) a ( P cell γ → γ = F ( N ) γ /F (0) γ ). Later w e will compare the ALP-photon mixing probabilities in s tochastic magnetic field using domain-like approac h with those obtained using the in tegral approach in next section. I II. STOCHASTIC MAGNETIC FIELD The p erturbative expansion of the probabilities in the constant magnetic field case, given in Eq. ( I I.9 ), motiv ates a similar treatment for stochastic magnetic field configurations. In this scenario, one m ust expand the solution up to O ( B 2 ) for P a → γ and up to O ( B 4 ) for P γ → γ (see App endix B). T o implemen t this, w e first replace the constan t mixing term ∆ B in Eq. ( I I.1 ) with a spatially v arying coun terpart and express the kernel matrix as K =   ω + ∆ γ 0 ∆ B x ( l ) 0 ω + ∆ γ ∆ B y ( l ) ∆ B x ( l ) ∆ B y ( l ) ω + ∆ a   , (I II.12) where ∆ B x ( l ) = 1 2 g aγ B x ( l ) and ∆ B y ( l ) = 1 2 g aγ B y ( l ) . The EoM and initial condition of the ev olution op erator remain Eq. I I.3 , i.e. ∂ l U ( l ) = iK ( l ) U ( l ) and U ( l 0 = 0) = I 3 × 3 . Due to the l -dep endent K ( l ) , the solution can no longer b e expressed in closed analytic form as in the constant magnetic field case. Consequently , we employ the p erturbative framew ork introduced in Refs. [ 21 , 23 , 31 , 33 ]. The k ernel is split as K ( l ) = K 0 + δ K ( l ) , where K 0 and δ K ( l ) contain the diagonal and off-diagonal comp onen ts, resp ectiv ely . Solving the EoM iteratively yields a Dyson p erturbation series U ( l ) = e ilK 0 + e ilK 0 i Z l 0 dl 1 e − il 1 K 0 δ K ( l 1 ) e il 1 K 0 + e ilK 0 i Z l 0 dl 1 e − il 1 K 0 δ K ( l 1 ) e il 1 K 0 i Z l 1 0 dl 2 e − il 2 K 0 δ K ( l 2 ) e il 2 K 0 + e ilK 0 i Z l 0 dl 1 e − il 1 K 0 δ K ( l 1 ) e il 1 K 0 i Z l 1 0 dl 2 e − il 2 K 0 δ K ( l 2 ) e il 2 K 0 i Z l 2 0 dl 3 e − il 3 K 0 δ K ( l 3 ) e il 3 K 0 + e ilK 0 i Z l 0 dl 1 e − il 1 K 0 δ K ( l 1 ) e il 1 K 0 i Z l 1 0 dl 2 e − il 2 K 0 δ K ( l 2 ) e il 2 K 0 i Z l 2 0 dl 3 e − il 3 K 0 δ K ( l 3 ) e il 3 K 0 i Z l 3 0 dl 4 e − il 4 K 0 δ K ( l 4 ) e il 4 K 0 + O  δ K 5  . (I II.13) In addition, the spatially v arying magnetic field is assumed to b e sto c hastic and non-helical, characterized by the t wo-point correlation function ⟨ B i ( x ) B j ( x ′ ) ⟩ = 1 (2 π ) 3 Z d 3 k e i k · ( x − x ′ )  δ ij − k i k k j k  P B ( k ) , (I II.14) where k = | k | and P B ( k ) is the p ow er sp ectrum. Note that all l -in tegrals in Eq. ( II I.13 ) are independent from the azim uthal angle φ ; th us w e can immediately integrate on φ and obtain ⟨ B x ( l ) B x ( l ′ ) ⟩ = ⟨ B y ( l ) B y ( l ′ ) ⟩ = π (2 π ) 3 Z k 2 dk Z 1 − 1 dte ikt ( l − l ′ )  1 + t 2  P B ( k ) ⟨ B x ( l ) B y ( l ′ ) ⟩ = ⟨ B y ( l ) B x ( l ′ ) ⟩ = 0 (I II.15) where t = cos θ with the θ the polar angle with resp ect to the l -direction. 6 FIG. 2: The ALP-photon mixing probabilities P a → γ and P γ → γ in the stochastic magnetic field with Gaussian dis- tribution and the mono c hromatic spectrum. All dimensional quantities are properly normalized to λ γ . Eq. I II.21 pro vides an appro ximation (colored dash lines) whic h that holds in the regime d/λ γ ≫ 1 . In the first plot illustrating ALP-photon con version, the initial linear gro wth is a signature of stochastic resonance (see Section II I for further details). The gray dashed lines corresp ond to e − d/λ γ and e − d/ (2 λ γ ) , ab ov e whic h the exp onen tial-decay component of the P γ → γ b ecomes sub dominan t and can be safely neglected compared to the non exp onen tial-decay contribution. F rom the perturbative expansion of U ( l ) in Eq. ( I II.13 ), w e extract the comp onen ts U 13 and U 23 and obtain the ALP-photon con version probability to O ( B 2 ) : P a → γ ( d ) = 1 4 g 2 aγ e − 2 bd Z d 0 dl 1 Z d 0 dl 2 e i ( l 1 − l 2 ) A e ( l 1 + l 2 ) b ( ⟨ B x ( l 1 ) B x ( l 2 ) ⟩ + ⟨ B y ( l 1 ) B y ( l 2 ) ⟩ ) = 1 4 g 2 aγ 1 (2 π ) 3 Z dk 2 π k 2 P B Z 1 − 1 dt  1 + t 2  1 + e − 2 bd − e − bd 2 cos (( A + kt ) d ) ( A + kt ) 2 + b 2 . (I I I.16) The complete expression after integrating ov er t is presen ted in Eq. ( C.1 ) in App endix C. The first term in the nominator corresponds to the non exp onential-deca y comp onent of P a → γ . The con v ersion probabilit y γ − ALP is given by P γ → a = 1 2 P a → γ . In the limit b → 0 , the result consistently reduces to the con ven tional ALP-photon con version in the sto c hastic magnetic field without EBL attenuation [ 33 ]. T urning to the photon surviv al probability , the p erturbativ e expansion must b e extended to quartic order in B , i.e., O ( B 4 ) : P γ → γ = P (0) γ → γ + P (1) γ → γ + P (2) γ → γ , (I II.17) 7 where P (0) γ → γ = e − 2 bd , P (1) γ → γ = − 1 2 1 4 g 2 aγ e − 2 bd Z d 0 dl 1 Z l 1 0 dl 2 e ( l 1 − l 2 ) b  e i ( l 1 − l 2 ) A + e − i ( l 1 − l 2 ) A  ( ⟨ B x ( l 1 ) B x ( l 2 ) ⟩ + ⟨ B y ( l 1 ) B y ( l 2 ) ⟩ ) , P (2) γ → γ = 1 2 1 16 g 4 aγ e − 2 bd  Z d 0 dl 1 Z l 1 0 dl 2 Z d 0 dl 3 Z l 3 0 dl 4 e i ( l 1 − l 2 ) A +( l 1 − l 2 ) b e − i ( l 3 − l 4 ) A +( l 3 − l 4 ) b ( ⟨ B x ( l 1 ) B x ( l 2 ) B x ( l 3 ) B x ( l 4 ) ⟩ + ⟨ B y ( l 1 ) B y ( l 2 ) B y ( l 3 ) B y ( l 4 ) ⟩ + ⟨ B x ( l 1 ) B y ( l 2 ) B x ( l 3 ) B y ( l 4 ) ⟩ + ⟨ B y ( l 1 ) B x ( l 2 ) B y ( l 3 ) B x ( l 4 ) ⟩ ) + Z d 0 dl 1 Z l 1 0 dl 2 Z l 2 0 dl 3 Z l 3 0 dl 4 e ( l 1 − l 2 + l 3 − l 4 ) b  e i ( l 1 − l 2 + l 3 − l 4 ) A + e − i ( l 1 − l 2 + l 3 − l 4 ) A  ( ⟨ B x ( l 1 ) B x ( l 2 ) B x ( l 3 ) B x ( l 4 ) ⟩ + ⟨ B x ( l 1 ) B y ( l 2 ) B y ( l 3 ) B x ( l 4 ) ⟩ + ⟨ B y ( l 1 ) B y ( l 2 ) B y ( l 3 ) B y ( l 4 ) ⟩ + ⟨ B y ( l 1 ) B x ( l 2 ) B x ( l 3 ) B y ( l 4 ) ⟩ ) } . (I II.18) Similarly to constant magnetic field case, we find that the non exp onential-deca y term only app ears at forth order B 4 of probability expansion, namely P (2) γ → γ (see App endix B). T o ev aluate the four p oint correlation function, we further assume that the sto c hastic magnetic field is Gaussian. Then we are allow ed to apply the Wick theorem to decomp ose four point correlation functions into pairings of t wo p oint correlation functions: ⟨ B x ( l 1 ) B x ( l 2 ) B x ( l 3 ) B x ( l 4 ) ⟩ = ⟨ B y ( l 1 ) B y ( l 2 ) B y ( l 3 ) B y ( l 4 ) ⟩ = ⟨ B x ( l 1 ) B x ( l 2 ) ⟩ ⟨ B x ( l 3 ) B x ( l 4 ) ⟩ + ⟨ B x ( l 1 ) B x ( l 3 ) ⟩ ⟨ B x ( l 2 ) B x ( l 4 ) ⟩ + ⟨ B x ( l 1 ) B x ( l 4 ) ⟩ ⟨ B x ( l 2 ) B x ( l 3 ) ⟩ , ⟨ B x ( l 1 ) B y ( l 2 ) B x ( l 3 ) B y ( l 4 ) ⟩ = ⟨ B x ( l 1 ) B x ( l 3 ) ⟩ ⟨ B x ( l 2 ) B x ( l 4 ) ⟩ , (I II.19) where w e ha ve used Eq. ( I II.15 ). In principle, one could substitute Eqs. ( I I I.19 ) and ( II I.15 ) in to Eq. ( I I I.18 ) and ev aluate the full integral. Ho w- ev er, the resulting expression is lengthy and tedious. Since our primary interest lies in the non-exponential-deca y con tribution, w e observ e that it emerges only from the following integral structure: R d 0 dl 1 R l 1 0 dl 2 R d 0 dl 3 R l 3 0 dl 4 , and w e obtain P (4) γ → γ → 1 2 1 16 g 4 aγ e − 2 bd Z d 0 dl 1 Z l 1 0 dl 2 Z d 0 dl 3 Z l 3 0 dl 4 e i ( l 1 − l 2 ) A +( l 1 − l 2 ) b e − i ( l 3 − l 4 ) A +( l 3 − l 4 ) b × (2 ⟨ B x ( l 1 ) B x ( l 2 ) ⟩ ⟨ B x ( l 3 ) B x ( l 4 ) ⟩ + 4 ⟨ B x ( l 1 ) B x ( l 3 ) ⟩ ⟨ B x ( l 2 ) B x ( l 4 ) ⟩ + 2 ⟨ B x ( l 1 ) B x ( l 4 ) ⟩ ⟨ B x ( l 2 ) B x ( l 3 ) ⟩ ) → 1 2 1 16 g 4 aγ e − 2 bd 2 Z d 0 dl 1 Z l 1 0 dl 2 Z d 0 dl 3 Z l 3 0 dl 4 e i ( l 1 − l 2 ) A +( l 1 − l 2 ) b e − i ( l 3 − l 4 ) A +( l 3 − l 4 ) b ⟨ B x ( l 1 ) B x ( l 2 ) ⟩ ⟨ B x ( l 3 ) B x ( l 4 ) ⟩ → 1 16 g 4 aγ π 2 (2 π ) 6 Z dk k 2 P B ( k ) Z dk ′ k ′ 2 P B ( k ′ ) Z 1 − 1 dt  1 + t 2  Z 1 − 1 dt ′  1 + t ′ 2  e id ( kt + k ′ t ′ ) ( A − ib + kt ) 2 ( A + ib − k ′ t ′ ) 2 , (I I I.20) where → signifies that we hav e discarded all terms exhibiting exp onential decay . The analytical expression after angular integrations ov er t and t ′ in the last line is given in Eq. ( C.2 ) in App endix C. F or the sake of completeness, we remark that the ALP surviv al probability is dominated by the unity term, P a → a ∼ 1 , within the p erturbativ e region. In this w ork, w e consider the mono chromatic sp ectrum of the magnetic field P B ( k ) = π 2 B 2 rms δ ( k − k B ) /k 2 B , where B rms = p ⟨ B 2 ⟩ and λ B = 2 π/k B are the a v erage strength and correlation length of the stochastic magnetic field, resp ectiv ely . W e use Eqs. ( I II.16 ) and ( I II .20 ) to ev aluate the probabilities for several representativ e parameter sets in Fig. 2 . F or ALP-photon conv ersion probability P a → γ , when d ≪ λ γ suc h that EBL attenuation effect is negligible, there is a linear growth with propagation distance as P a → γ ∼ d for λ B ≲ l osc ≃ | A | − 1 and the conv en tional quadratic gro wth P a ↔ γ ∼ d 2 for λ B ≳ l osc . This linear growth reflects the phenomenon of sto c hastic resonance in particle mixing within a sto c hastic background [ 23 , 33 ]. In the regime d ≫ λ γ , the con version probability P a → γ approac hes a constan t v alue that generally decreases with λ B , scaling approximately as P a → γ ∼ λ B . F or P γ → γ , we hav e extracted only the non-exp onential-deca y comp onent (see Eq. ( I II.20 )). Consequen tly , the results are reliable only at sufficiently large distances d ≫ λ γ , where the exponential-deca y con tribution is strongly suppressed. This corresponds to the region ab ov e the dashed lines in Fig. 2 . In the limit of large correlation length λ B , the photon surviv al probability in a sto c hastic magnetic field approaches that of the constan t magnetic field case, as exp ected. As λ B decreases, P γ → γ is suppressed, scaling appro ximately as P γ → γ ∼ λ 6 B . Moreov er, in con trast to P a → γ , which saturates at large distances, the photon surviv al probability decays with distance as P γ → γ ∼ 1 /d 2 . The different scaling behaviors of P a → γ and P γ → γ with resp ect to λ B and d lead to a mo dified relation P γ → γ ≲ P 2 a ↔ γ in the sto c hastic case, in con trast to the constan t magnetic field scenario where P γ → γ ∼ P 2 γ ↔ a (see Figs. 1 and 2 ). Motiv ated by the key features observed in Fig. 2 and the analytical estimates of the non-exp onential-deca y mo de deriv ed in Appendix B, we construct simple parametric form ulas that appro ximate the probabilities at sufficiently 8 large distances within the p erturbativ e regime: P param a → γ ≃ O (1) × g 2 aγ B 2 λ 2 γ 1 1 + λ 2 γ A 2 + λ γ k B , P param γ → γ ≃ O (1) × g 4 aγ B 4 λ 4 γ 1 (1 + k 2 B d 2 )  1 + λ 4 γ ( | A | − k B ) 4  . (I II.21) As as sho wn in Fig. 2 , b oth P param a → γ and P param γ → γ matc h the numerical results with fine accuracy . Although P param a → γ deviates in the resonance region | A | = k B , the discrepancy remains within an order of magnitude ( O (10) ), which is sufficien t for leading-order estimation purposes. In fact, the exact form of the non-exp onential decay derived from Eq. ( II I.16 ) inv olv es an arctangent function (see Eq. ( C.1 ) in Appendix C). In a realistic cosmic environmen t, the mixing elements in Eq. ( I I.1 ) are given in prop er units b y [ 9 , 30 ] ∆ pl ≃ − 10 − 11  T e V ω   n e 10 − 7 cm − 3  Mp c − 1 , ∆ QED ≃ 4 × 10 − 8  ω T e V   B 10 − 9 Gs  2 Mp c − 1 , ∆ CMB ≃ 8 × 10 − 2  ω T e V  Mp c − 1 , ∆ a ≃ − 8  m a 10 − 8 e V  2  T e V ω  Mp c − 1 , ∆ B ≃ 5 × 10 − 3  g aγ 10 − 12 Ge V − 1   B 10 − 9 Gs  Mp c − 1 . (I II.22) Note that A = ∆ pl + ∆ QED + ∆ CMB − ∆ a . Here ∆ pl accoun ts for plasma effects in the in tergalactic medium, and ∆ QED denotes the QED v acuum p olarization correction for a photon in the magnetic background. In the intergalactic region, b oth terms are t ypically muc h smaller than the ALP mass term ∆ a . The only contribution that can induce a sizable effect is ∆ CMB , whic h accounts for the electromagnetic energy densit y provided by CMB [ 34 ]. In particular, A approac hes to zero at energy ω ≃ m a /  10 − 9 e V  T e V , leading to a large oscillation length l osc as shown in left panel in Fig. 3 . Let us comment on the p erturbativ e condition in realistic settings. F rom the analysis in App endix B, the p erturba- tiv e regime for a stochastic magnetic field requires ∆ 2 B λ γ d ≪ 1 . This condition can b e satisfied in most of the parameter space of our in terests for ALP-photon mixing ov er the intergalactic scale. F or instance, taking g aγ = 10 − 12 Ge V − 1 , B rms = 10 − 9 Gs and d = Gp c, the perturbative condition implies λ γ ≪ 40 Mp c, corresp onding to photon energy ω ≳ 10 T e V (ligh t green line in left panel in Fig. 3 ). A smaller pro duct g aγ B rms further enlarges the p erturbative region. In Fig. 3 , w e compare the ALP-photon con version and photon surviv al probabilities in the sto c hastic magnetic field treated b y the simplified domain-like mo del (Eq. ( II.11 )) and the integral approac h developed in this w ork (Eqs ( II I.16 ) and ( II I.20 ) ). F or comparison, the unphysical scenario of a constan t magnetic field extending across the entire in tergalactic region is shown as a comparative reference. F or the constant magnetic field, b oth P a → γ and P γ → γ exhibit p eaks caused by the mass-equal resonance at | A | = 0 . In con trast to the simple p eak structure with its well-defined resonance pattern, the scenario b ecomes considerably more complex when a sto chastic magnetic field is in tro duced. Within the domain-lik e approach, b oth P a → γ and P γ → γ exhibit a single-peak structure similar to the constan t-field case for m a ≳ 10 − 7 e V. How ever, for m a ≲ 10 − 8 e V, the probabilities display a distinct feature: they decrease with increasing energy , follow ed by a small p eak around ω ≃ 100 T e V. This feature may b e attributed to the resonance condition | A | ≃ k B , i.e., l osc ≃ λ B . Such a p eak in the photon surviv al probabilit y , as seen in the simplified domain-lik e mo del (right panel of Fig. 3 ), has also b een observed in more sophisticated domain-like treatments [ 9 , 30 ]. F or the in tegral approac h, we again assume a mono c hromatic s pectrum for the Gaussian magnetic field and present the results (thick curves) within the p erturbativ e regime, as sho wn in Fig. 3 . The amplitude of P a → γ is comparable to that obtained in domain-lik e mo del across the considered mass cases. Combining the left panel in Fig. 3 with upp er ro w in Fig. 2 , the p eak or knee structure in P a → γ originates from the sto chastic resonance when | A | − 1 ≳ k − 1 B , i.e. λ B ≲ l osc . F or photon surviv al probability , tw o notable features emerge: i) P γ → γ is suppressed by several orders of magnitude relative to the domain-like mo del; ii) the p eak structures of all considered masses cases could b e explained b y the resonance at | A | ≃ k B . In particular, a clear double-p eak structure emerges for m a = 10 − 7 e V. As shown in Fig. 3 , within the p erturbative regime the domain-lik e treatment yields the same quadratic amplitude relation as in the constant-field case, namely P γ → γ ≃ P 2 a → γ . This b eha vior can b e interpreted as arising from a double oscillation pro cess γ → a → γ . In con trast, the in tegral approach rev eals that this quadratic relation no longer holds, 9 FIG. 3: Left: The c haracteristic length scale 1 / | A | , in the oscillation Hamiltonian, for different ALP mass cases and photon mean free path λ γ as a function of energy . Middle and Righ t: the ALP-photon mixing probabilities P a → γ and P γ → γ in three magnetic field configurations: a constant magnetic field filling the en tire space, a simplified domain-lik e mo del, and a sto chastic magnetic field. The parameters are set to B rms = 10 − 9 Gs, g aγ = 10 − 12 Ge V − 1 and d = Gp c. F or the stochastic magnetic field, its distribution is Gaussian and the spectrum is monochromatic with correlation length λ B = 1 Mp c. The blac k dotted line in right panel denotes the photon surviv al probabilit y in the absence of ALP mixing. The thick curv es are constrained to ω ≳ 20 T e V where the p erturbative condition is w ell satisfied in the sto chastic case. and the p eak structures of P a → γ and P γ → γ differ significantly . By consisten tly preserving the sto chastic nature of the magnetic field from the outset, the in tegral approach captures the nontrivial interpla y b et w een EBL attenuation and ALP-photon mixing in a stochastic magnetic environmen t. These results highlight the imp ortance of moving b ey ond simplified domain-like approximations when mo deling magnetic turbulence and motiv ate further in vestigation in future work. IV. NON-GA USSIAN MAGNETIC FIELD In the previous section, w e considered a Gaussian random magnetic field. F rom now on, w e consider the case of a non-Gaussian distribution. In this case, the four-p oin t correlation function of the magnetic field can no longer b e fully decomposed into products of tw o-point functions via Wic k’s theorem (Eq. ( I I I.19 )). Instead, it includes an additional nonv anishing connected comp onen t, ⟨ B ( l 1 ) B ( l 2 ) B ( l 3 ) B ( l 4 ) ⟩ c , which corresponds to the trisp ectrum in momen tum space. A detailed treatmen t of the trispectrum is challenging, as its sp ecific form depends on the underlying magnetogenesis mechanism—whether of primordial origin or generated by astrophysical pro cesses. T o qualitativ ely assess the impact of magnetic non-Gaussianity on the photon surviv al probabilit y in ALP-photon mixing, we therefore adopt three heuristic approaches. F or a qualitative analysis, we treat the magnetic vector field sim ply as a scalar in three-dimensional space, denoted B ( x ) . F ollowing the deriv ation leading to Eqs. ( I I I.20 ), the non-exp onential-deca y contribution arises from an integral of the form R d 0 dl 1 R l 1 0 dl 2 R d 0 dl 3 R l 3 0 dl 4 . F rom this, the relev an t part of P γ → γ can be extracted as P NG γ → γ = 4 1 2 e − 2 bd 1 16 g 4 aγ Z d 0 dl 1 Z l 1 0 dl 2 Z d 0 dl 3 Z l 3 0 dl 4 e i ( l 1 − l 2 ) A +( l 1 − l 2 ) b e − i ( l 3 − l 4 ) A +( l 3 − l 4 ) b ⟨ B ( l 1 ) B ( l 2 ) B ( l 3 ) B ( l 4 ) ⟩ c , (IV.23) where the ov erall prefactor 4 accounts for scalar approximation B x ( x ) = B y ( x ) = B ( x ) . Non-Gaussian c ase 1 . The simplest case is to assume a constant 4-p oin t function ⟨ B ( l 1 ) B ( l 2 ) B ( l 3 ) B ( l 4 ) ⟩ c = κB 4 rms . (IV.24) Under this assumption, we obtain P NG γ → γ = 4 1 2 e − 2 db 1 16 g 4 aγ κB 4 rms Z d 0 dl 1 Z l 1 0 dl 2 Z d 0 dl 3 Z l 3 0 dl 4 e i ( l 1 − l 2 ) A +( l 1 − l 2 ) b e − i ( l 3 − l 4 ) A +( l 3 − l 4 ) b → 1 8 κg 4 aγ B 4 rms 1 ( A 2 + b 2 ) 2 , (IV.25) 10 FIG. 4: The photon surviv al probabilit y P γ → γ in the stochastic magnetic field with a mono chromatic spectrum with λ B = 1 Mp c, for Gaussian (thick line) and non-Gaussian distributions (case-1 for dashed line, case-2 for dotted line and case-3 for thin line). The left panel is at fixed ω = 100 T e V, and the middle and right panels is fixed at d = Gp c. The middle panel corresp onds to the weak non-Gaussianit y , where parameters are set to κ = 0 . 01 for case-1, κ = 0 . 01 , λ ρ = 0 . 1 Mpc for case-2, and g N L = 0 . 01 for case-3. The righ t panel corresp onds to strong non- Gaussianit y , where parameters are set to κ = 100 for case-1, κ = 1 , λ ρ = 50 Mp c for case-2, and g N L = 100 for case-3. Note that non-gaussianities highly b o ost the photon surviv al probability with mo del dep enden t profiles. where → denotes the non exp onen tial-decay part ( db < 1) . The parameter κ is the kurtosis of magnetic field, κ =  B 4  /  B 2  2 − 3 , which is commonly used to quantify intermittency in MHD simulations [ 35 ]. Non-Gaussian c ase 2 . In the literature, magnetic non-Gaussianit y is often c haracterized via the magnetic energy densit y ρ ( x ) = B 2 ( x ) or the stress-energy tensor [ 36 , 37 ]. A dopting a similar strategy , w e introduce an ansatz that incorp orates a Dirac delta function b y hand: ⟨ B ( l 1 ) B ( l 2 ) B ( l 3 ) B ( l 4 ) ⟩ c = ⟨ ρ ( l 1 ) ρ ( l 2 ) ⟩ c δ ( l 1 /l 3 − 1) δ ( l 2 /l 4 − 1) + ⟨ ρ ( l 1 ) ρ ( l 3 ) ⟩ c δ ( l 1 /l 2 − 1) δ ( l 3 /l 4 − 1) + ⟨ ρ ( l 1 ) ρ ( l 2 ) ⟩ c δ ( l 1 /l 4 − 1) δ ( l 2 /l 3 − 1) , (IV.26) where ⟨ ρ ( l 1 ) ρ ( l 2 ) ⟩ c = ⟨ ρ ( l 1 ) ρ ( l 2 ) ⟩ − ⟨ ρ ( l 1 ) ⟩ ⟨ ρ ( l 2 ) ⟩ is the connected t wo p oin t correlation function of the magnetic energy density . This quan tit y characterizes the spatial distribution of magnetic energy . A long-range connected correlation ⟨ ρρ ⟩ c indicates a relatively smo oth and homogeneous field configuration, whereas a short-range correlation corresp onds to a tangled and intermitten t field. In general, fluctuations in the magnetic energy deca y rapidly b eyond the intrinsic correlation length of the magnetic field. In momentum space, the p o wer sp ectrum of the magnetic energy is defined as ⟨ ρ ( x 1 ) ρ ( x 2 ) ⟩ c = Z d 3 k e i ( x 1 − x 2 ) · k P ρ ( k ) . (IV.27) Substituting Eqs. ( IV.26 ) and ( IV.27 ) into Eq. ( IV.23 ), w e find that only the com bination δ ( l 1 /l 3 − 1) δ ( l 2 /l 4 − 1) con tributes to the non-exponential-deca y part P NG γ → γ → 4 1 2 e − 2 db 1 16 g 4 aγ Z d 0 dl 1 Z l 1 0 dl 2 Z d 0 dl 3 Z l 3 0 dl 4 e 2( l 1 − l 2 ) b ⟨ ρ ( l 1 ) ρ ( l 2 ) ⟩ c δ ( l 1 /l 3 − 1) δ ( l 2 /l 4 − 1) → π 4 g 4 aγ Z dk k 2 P ρ ( k ) Z 1 − 1 dt − e idkt (1 − 2 bd − idk t ) (2 b + ikt ) 4 . (IV.28) Its expression after integration o ver t is given in Eq. ( C.3 ) in App endix C. In this work, we consider a mono chromatic sp ectrum of the magnetic energy density , P ρ ( k ) = κB 4 rms δ ( k − k ρ ) /k 2 with characteristic wa v elength λ ρ = 2 π /k ρ . W e set λ ρ = 0 . 1 Mp c for intergalactic magnetic field with t ypical λ B = 1 Mp c. Non-Gaussian c ase 3 . W e follo w the standard p erturbative approac h emplo y ed in the analysis of non-Gaussianit y in primordial curv ature p erturbations and the CMB trisp ectrum. T o this end, we first rescale B to a dimensionless quan tity B ( x ) = B ( x ) /B rms . A common parametrization of non-Gaussianity is given by a lo cal expansion in terms of underlying Gaussian random fields [ 38 , 39 ]: B ( x 1 ) = B G + 1 2 √ τ N L  B 2 G −  B 2 G  + 9 25 g N L B 3 G + ..., (IV.29) 11 where τ N L and g N L are dimensionless constants. The connected 4-p oint correlation function is ⟨B ( x 1 ) B ( x 2 ) B ( x 3 ) B ( x 4 ) ⟩ c = 1 (2 π ) 12 Z d 3 k 1 Z d 3 k 2 Z d 3 k 3 Z d 3 k 4 (2 π ) 3 δ (3) ( k 1 + k 2 + k 3 + k 4 ) × e i ( k 1 · x 1 + k 2 · x 2 + k 3 · x 3 + k 4 · x 4 ) T ( k 1 , k 2 , k 3 , k 4 ) (IV.30) with T ( k 1 , k 2 , k 3 , k 4 ) = τ N L X a  = b,c & b 10 T e V photons from GRB 221009A observ ed b y LHAASO in v oking ALPs. Our results establish that sto c hastic magnetic fields cannot b e reduced to domain-like coherence without losing essen tial physics, and that high-energy γ -ray sp ectra enco de statistically accessible information ab out b oth the p o wer sp ectrum and the non-Gaussian structure of in tergalactic magnetic fields. The quartic scaling P γ → γ ∼ B 4 within the perturbative regime naturally connects photon surviv al to four-p oint magnetic correlations, offering a direct observ ational window into magnetic non-Gaussianity . As next-generation observ atories push tow ard P e V sensitivities, the effects identified here will b ecome essen tial for distinguishing b et ween astrophysical atten uation and new physics in the extragalactic very-high-energy sky . The distinctiv e sp ectral signatures predicted b y our framew ork—peak structures from m ultiple resonances, non-exponential plateaus, and non-Gaussian enhancements—pro vide concrete targets for LHAASO, CT A, and future exp eriments to prob e b oth ALP ph ysics and the statistical prop erties of in tergalactic magnetic fields. A c knowledgemen ts. AA w ork is supp orted b y the National Science F oundation of China (NSFC) through the gran t No. 12350410358; the T alent Scientific Researc h Program of College of Ph ysics, Sich uan Universit y , Grant No. 1082204112427; the F ostering Program in Disciplines P ossessing No v el F eatures for Natural Science of Sich uan Univ ersity , Gran t No.2020SCUNL209 and 1000 T alent program of Sich uan pro vince 2021. YFC is supported in part by the National Key R&D Program of China (2021YF C2203100), by NSFC (12433002), b y CAS young interdisciplinary inno v ation team (JCTD-2022-20), by 111 Pro ject (B23042), and b y USTC F ellowship for International Co op eration. SC, QG, and GL ackno wledge the Istituto Nazionale di Fisica Nucleare (INFN), Sezione di Nap oli, iniziative sp e cifiche QGSKY and MOONLIGHT2. This paper is based up on w ork from COST Action CA21136 A ddr essing observational tensions in c osmolo gy with systematics and fundamental physics (CosmoV erse) supp orted by Europ ean Co op eration in Science and T echnology . App endix A In this app endix section, we deriv e Eq. ( II.5 ), which gives the relation b et ween the relev ant probabilities in ALP- photon mixing and the ev olution op erator. W e b egin with a brief review of the w av efunction, evolution op erator and densit y matrix [ 9 , 41 ]. Given a quantum state | Ψ( t ) ⟩ , its dynamics satisfy the Schrodinger equation i d dt | Ψ( t ) ⟩ = H ( t ) | Ψ( t ) ⟩ . (A.1) In terms of the evolution op erator U ( t, t 0 ) giv en b y | Ψ( t ) ⟩ = U ( t, t 0 ) | Ψ( t 0 ) ⟩ , its equation reads i d dt U ( t, t 0 ) = H ( t ) U ( t, t 0 ) . (A.2) 14 Define the densit y matrix ρ ( t ) = | Ψ( t ) ⟩ ⟨ Ψ( t ) | , its relation to ev olution operator is ρ ( t ) = U ( t, t 0 ) ρ ( t 0 ) U † ( t, t 0 ) . It satisfies i d dt ρ ( t ) = H ( t ) ρ ( t ) − ρ ( t ) H † ( t ) . (A.3) Note that EoM for the wa v efunction (Eq. ( A.1 )), the evolution op erator (Eq. ( A.2 )) and density matrix (Eq. A.3 ) are equiv alen t. T o illustrate the physical meaning of density matrix, we consider the initial pure state | Ψ( t 0 ) ⟩ = | φ ⟩ with density matrix ρ ( t 0 ) = | φ ⟩ ⟨ φ | . Using Eq. ( A.3 ), one obtains the density matrix ρ ( t ) at arbitrary time t . The probabilit y of detecting the sate | ψ ⟩ at t is P ( t ) = ⟨ ψ | ρ ( t ) | ψ ⟩ = T r ( | ψ ⟩ ⟨ ψ | ρ ( t )) , (A.4) where T r is the trace op erator. In other words, P ( t ) can b e interpreted as the con version probability from the initial state | φ ⟩ to the final state | ψ ⟩ . Applying this formalism to the APL-photon mixing, one can see that the Sc hro dinger equation Eq. A.1 is exactly same as Eq. ( I I.1 ) b y replacing time t by propagation distance l , and the Hamiltonian H with kernel matrix − K . A general quan tum states can b e written as | Ψ( l ) ⟩ = c 1 ( l ) | A x ⟩ + c 2 ( l ) | A y ⟩ + c 3 ( l ) | a ⟩ , (A.5) where | A x ⟩ =   1 0 0   , | A y ⟩ =   0 1 0   , | a ⟩ =   0 0 1   (A.6) form the basis representing the eigenstates of t wo photon p olarization states and the ALP state. Assuming the initial condition is a pure ALP state, | Ψ( l 0 ) ⟩ = | a ⟩ , the corresp onding density matrix is ρ ( l 0 ) = | Ψ( l 0 ) ⟩ ⟨ Ψ( l 0 ) | =   0 0 0 0 0 0 0 0 1   . (A.7) One can either solve Eq. ( A.3 ) to directly obtain ρ ( l ) , or equiv alen tly , first solve U ( t, t 0 ) from Eq. ( A.2 ) and then compute ρ ( l ) = U ( l , l 0 ) ρ ( l 0 ) U † ( l, l 0 ) . A t distance l , the probabilit y to obtain the states | A x ⟩ and | A y ⟩ is P a → γ ( l ) = ⟨ A x | ρ ( l ) | A x ⟩ + ⟨ A y | ρ ( l ) | A y ⟩ = T r ( | A x ⟩ ⟨ A x | ρ ( l )) + T r ( | A y ⟩ ⟨ A y | ρ ( l )) = T r     1 0 0 0 0 0 0 0 0   U ( l , l 0 ) ρ ( l 0 ) U † ( l, l 0 )   + T r     0 0 0 0 1 0 0 0 0   U ( l , l 0 ) ρ ( l 0 ) U † ( l, l 0 )   = U 13 U ∗ 13 + U 23 U ∗ 23 , (A.8) where in last step we used the matrix form of U and ρ ( l 0 ) to take the trace op eration. The probabilit y P a → γ ( l ) can b e interpreted as the conv ersion probability from ALP to photon with t wo p olarization modes after mixing system tra vels ov er a distance d = l − l 0 . Similarly , the ALP surviv al probability is P a → a = T r  | a ⟩ ⟨ a | U ( l , l 0 ) ρ ( l 0 ) U † ( l, l 0 )  = U 33 U ∗ 33 . (A.9) No w consider an initially unp olarized photon state with density matrix given by ρ ( l 0 ) = 1 2 | A x ⟩ ⟨ A x | + 1 2 | A y ⟩ ⟨ A y | = 1 2   1 0 0 0 1 0 0 0 0   , (A.10) using the same pro cedure, w e obtain the photon surviv al probabilit y and photon-ALP con v ersion probability as P γ → γ = T r  | A x ⟩ ⟨ A x | U ( l , l 0 ) ρ ( l 0 ) U † ( l, l 0 )  + T r  | A y ⟩ ⟨ A y | U ( l , l 0 ) ρ ( l 0 ) U † ( l, l 0 )  = 1 2 ( U 11 U ∗ 11 + U 12 U ∗ 12 + U 21 U ∗ 21 + U 22 U ∗ 22 ) , P γ → a = T r  | a ⟩ ⟨ a | U ( l , l 0 ) ρ ( l 0 ) U † ( l, l 0 )  = 1 2 ( U 31 U ∗ 31 + U 32 U ∗ 32 ) . (A.11) 15 App endix B In this app endix section, w e derive the p erturbation expansion of the relev an t probabilities with resp ect to the dimensionless parameter g aγ B λ γ . An order-of-magnitude analysis is p erformed to determine the p erturbativ e regime. F or ALP-photon mixing in the constant magnetic field, expanding Eq. ( I I.8 ) yields P γ ↔ a ∼ ( g aγ B λ γ ) 2 1 1 + A 2 λ 2 γ  1 + e − d λ γ + e − d 2 λ γ  + ( g aγ B λ γ ) 4 1  1 + A 2 λ 2 γ  3  1 + d λ γ   1 + A 2 λ 2 γ  + e − d λ γ + e − d 2 λ γ  + ( g aγ B λ γ ) 6 1  1 + A 2 λ 2 γ  5  1 + d λ γ + d 2 λ 2 γ   1 + A 2 λ 2 γ + A 4 λ 4 γ  + e − d λ γ + e − d 2 λ γ  + o  ( g aγ B λ γ ) 8  , (B.1) and P γ → γ ∼ e − d λ γ + ( g aγ B λ γ ) 2 1  1 + A 2 λ 2 γ  2 n e − d λ γ + e − d 2 λ γ o + ( g aγ B λ γ ) 4 1  1 + A 2 λ 2 γ  4 n  1 + A 2 λ 2 γ + A 4 λ 4 γ  + e − d λ γ + e − d 2 λ γ o + ( g aγ B λ γ ) 6 1  1 + A 2 λ 2 γ  6  1 + d λ γ   1 + A 2 λ 2 γ + A 4 λ 4 γ + A 6 λ 6 γ  + e − d λ γ + e − d 2 λ γ  + ( g aγ B λ γ ) 8 1  1 + A 2 λ 2 γ  8  1 + d λ γ + d 2 λ 2 γ   1 + A 2 λ 2 γ + A 4 λ 4 γ + A 6 λ 6 γ + A 8 λ 8 γ  + e − d λ γ + e − d 2 λ γ  + o  ( g aγ B λ γ ) 10  , (B.2) and P a → a ∼ 1 + ( g aγ B λ γ ) 2 1  1 + A 2 λ 2 γ  2  1 + d λ γ   1 + A 2 λ 2 γ  + e − d 2 λ γ  + o  ( g aγ B λ γ ) 4  . (B.3) In the strong EBL attenuation regime, d ≫ λ γ , the exp onential-deca y term is muc h suppressed, then we obtain P γ ↔ a ∼ ( ( g aγ B λ γ ) 2 + ( g aγ B λ γ ) 4 d λ γ + ( g aγ B λ γ ) 6 d 2 λ 2 γ + ... Aλ γ ≲ O (1) , ( g aγ B λ γ ) 2 1 ( Aλ γ ) 2 + ( g aγ B λ γ ) 4 d λ γ 1 ( Aλ γ ) 4 + ( g aγ B λ γ ) 6 d 2 λ 2 γ 1 ( Aλ γ ) 6 + ... Aλ γ ≫ 1 , P γ → γ ∼ ( ( g aγ B λ γ ) 4 + ( g aγ B λ γ ) 6 d λ γ + ( g aγ B λ γ ) 8 d 2 λ 2 γ + ... Aλ γ ≲ O (1) , ( g aγ B λ γ ) 4 1 ( Aλ γ ) 4 + ( g aγ B λ γ ) 6 d λ γ 1 ( Aλ γ ) 6 + ( g aγ B λ γ ) 8 d 2 λ 2 γ 1 ( Aλ γ ) 8 + ... Aλ γ ≫ 1 , P a → a ∼ ( 1 + ( g aγ B λ γ ) 2 d λ γ + ... Aλ γ ≲ O (1) , 1 + ( g aγ B λ γ ) 2 d λ γ 1 ( Aλ γ ) 2 + ... Aλ γ ≫ 1 . (B.4) It shows that the p erturbation condition is ( g aγ B λ γ ) 2 d λ γ ≪ 1 for Aλ γ ≲ O (1) , and ( g aγ B λ γ ) 2 d λ γ 1 ( Aλ γ ) 2 ≪ 1 for Aλ γ ≫ 1 . F or order-of-magnitude estimates, these t wo region Aλ γ ≲ O (1) and Aλ γ ≫ 1 shown ab ov e can b e com bined into a compact expression given in Eq. ( I I.9 ). In the short distance d < λ γ , the EBL absorption can b e neglected, and Eq. ( B.1 ) reduces to P γ ↔ a ∼ ( ( g aγ B λ γ ) 2 + ( g aγ B λ γ ) 4 + ( g aγ B λ γ ) 6 + ... Aλ γ ≲ O (1) , ( g aγ B λ γ ) 2 1 ( Aλ γ ) 2 + ( g aγ B λ γ ) 4 1 ( Aλ γ ) 4 + ( g aγ B λ γ ) 6 1 ( Aλ γ ) 6 + ... Aλ γ ≫ 1 . (B.5) 16 The p erturbation condition b ecomes ( g aγ B λ γ ) 2 ≪ 1 for Aλ γ ≲ O (1) , and ( g aγ B λ γ ) 2 1 ( Aλ γ ) 2 ≪ 1 for Aλ γ ≫ 1 . F or sufficien t weak EBL attenuation where λ γ is large enough at Aλ γ ≫ 1 , the p erturbation condition g aγ B l osc ≪ 1 ( l osc ∼ | A | − 1 ) corresp onds to w eak mixing angle of ALP-photon oscillation [ 31 ]. Regarding the photon surviv al probabilit y , the unitary conserv ation in weak EBL regime requires P γ → a + P γ → γ ≃ 1 . Therefore, the leading order of photon surviv al probability is around unity for d < λ γ . A similar analysis applies to sto c hastic magnetic field. Although the general treatment in volv es complicate integrals o ver R d 3 k and R dl , w e adopt tw o simplifications that would not alter the main conclusion. First, w e study the ALP-photon mixing in t wo spatial dimensions instead of three. Second, we fo cus on the mono c hromatic spectrum of magnetic field, with correlation function ⟨ B ( l ) B ( l ′ ) ⟩ ∼ e ik B ( l − l ′ ) , whic h allows the magnetic field to b e treated effectiv ely as an oscillatory field B ( l ) ∼ B e ik B l .Note that here B ( l ) is treated as a complex field; consequently , the complex conjugate B ∗ ( l ) must b e included in the evolution op erator where appropriate. Within these simplifications, from Eq. ( I I I.13 ) w e obtain U 12 ∼ g aγ e − bd Z d 0 dl 1 e l 1 b e − il 1 A B ( l 1 ) + g 3 aγ e − bd Z d 0 dl 1 Z l 1 0 dl 2 Z l 2 0 dl 3 e ( l 1 − l 2 + l 3 ) b e − i ( l 1 − l 2 + l 3 ) A B ( l 1 ) B ∗ ( l 2 ) B ( l 3 ) + g 5 aγ e − bd Z d 0 dl 1 Z l 1 0 dl 2 Z l 2 0 dl 3 Z l 3 0 dl 4 Z l 4 0 dl 5 e ( l 1 − l 2 + l 3 − l 4 + l 5 ) b e − i ( l 1 − l 2 + l 3 − l 4 + l 5 ) A B ( l 1 ) B ∗ ( l 2 ) B ( l 3 ) B ∗ ( l 4 ) B ( l 5 ) + ... ∼ g aγ B e − bd Z d 0 dl 1 e l 1 b e − il 1 ( A − k B ) + g 3 aγ B 3 e − bd Z d 0 dl 1 Z l 1 0 dl 2 Z l 2 0 dl 3 e ( l 1 − l 2 + l 3 ) b e − i ( l 1 − l 2 + l 3 )( A − k B ) + g 5 aγ B 5 e − bd Z d 0 dl 1 Z l 1 0 dl 2 Z l 2 0 dl 3 Z l 3 0 dl 4 Z l 4 0 dl 5 e ( l 1 − l 2 + l 3 − l 4 + l 5 ) b e − i ( l 1 − l 2 + l 3 − l 4 + l 5 )( A − k B ) + ..., (B.6) and U 11 ∼ e − bd + g 2 aγ B 2 e − bd Z d 0 dl 1 Z l 1 0 dl 2 e ( l 1 − l 2 ) b e − i ( l 1 − l 2 )( A − k B ) + g 4 aγ B 4 e − bd Z d 0 dl 1 Z l 1 0 dl 2 Z l 2 0 dl 3 Z l 3 0 dl 4 e ( l 1 − l 2 + l 3 − l 4 ) b e − i ( l 1 − l 2 + l 3 − l 4 )( A − k B ) + g 6 aγ B 6 e − bd Z d 0 dl 1 Z l 1 0 dl 2 Z l 2 0 dl 3 Z l 3 0 dl 4 Z l 4 0 dl 5 Z l 5 0 dl 6 e ( l 1 − l 2 + l 3 − l 4 + l 5 − l 6 ) b e − i ( l 1 − l 2 + l 3 − l 4 + l 5 − l 6 )( A − k B ) + g 8 aγ B 8 e − bd Z d 0 dl 1 Z l 1 0 dl 2 Z l 2 0 dl 3 Z l 3 0 dl 4 Z l 4 0 dl 5 Z l 5 0 dl 6 Z l 6 0 dl 7 Z l 7 0 dl 8 × e ( l 1 − l 2 + l 3 − l 4 + l 5 − l 6 + l 7 − l 8 ) b e − i ( l 1 − l 2 + l 3 − l 4 + l 5 − l 6 + l 7 − l 8 )( A − k B ) + .... (B.7) In the strong EBL region d ≫ λ γ , w e drop the exp onential-deca y terms and obtain P a → γ ∼ | U 12 | 2 ∼    ( g aγ B λ γ ) 2 + ( g aγ B λ γ ) 4 d λ γ + ( g aγ B λ γ ) 6 d 2 λ 2 γ + ... e Aλ γ ≲ O (1) , ( g aγ B λ γ ) 2 1 ( e Aλ γ ) 2 + ( g aγ B λ γ ) 4 d λ γ 1 ( e Aλ γ ) 4 + ( g aγ B λ γ ) 6 d 2 λ 2 γ 1 ( e Aλ γ ) 6 + ... e Aλ γ ≫ 1 , (B.8) and P γ → γ ∼ | U 11 | 2 ∼    ( g aγ B λ γ ) 4 + ( g aγ B λ γ ) 6 d λ γ + ( g aγ B λ γ ) 8 d 2 λ 2 γ + ... e Aλ γ ≲ O (1) , ( g aγ B λ γ ) 4 1 ( e Aλ γ ) 4 + ( g aγ B λ γ ) 6 d λ γ 1 ( e Aλ γ ) 6 + ( g aγ B λ γ ) 8 d 2 λ 2 γ 1 ( e Aλ γ ) 8 + ... e Aλ γ ≫ 1 , (B.9) where e A ∼ A − k B . Similarly to the constant magnetic field case (Eq. ( B.4 )), the p erturbativ e condition is ( g aγ B λ γ ) 2 d λ γ 1 1+ ( e Aλ γ ) 2 ≪ 1 in the strong EBL attenuation regime d ≫ λ γ . F or typical in tergalactic parameters, e.g. g aγ = 10 − 12 Ge V − 1 , B = 10 − 9 Gs, d = Gp c and ω ≳ 20 T e V, this conserv ativ ely satisfies ( g aγ B λ γ ) 2 d λ γ ≲ 1 . Reducing 17 g aγ or B b y an order of magnitude ensures that the p erturbative regime is well satisfied, extending the v alid energy range down to ω ∼ 0 . 1 T e V. Moreo ver, the leading contribution to P γ → γ in Eq. ( B.9 ), ( g aγ B λ γ ) 4 /  1 + λ 4 γ ( A − k B ) 4  (v alid for all e Aλ γ ), partially justifies the parameterized expression of P γ → γ giv en in Eq. ( II I.21 ). App endix C In this app endix section, w e p resen t the complete expression of several formula quoted in the main text. The full expression of P a → γ in the Gaussian magnetic field case is (Eq. ( I I I.16 ) ) P a → γ = 1 4 g 2 aγ 1 (2 π ) 3 Z dk 2 πk 2 P B Z 1 − 1 dt  1 + t 2  1 + e − 2 bd − e − bd 2 cos (( A + k t ) d ) ( A + kt ) 2 + b 2 = 1 4 g 2 aγ 1 (2 π ) 3 Z dk 2 πk 2 P B 1 bk 3 ×   A 2 − b 2 + k 2   arctan  A + k b  − arctan  A − k b  + b  2 k + A ln  ( A − k ) 2 + b 2 ( A + k ) 2 + b 2  + 1 4 g 2 aγ e − 2 bd 1 (2 π ) 3 Z dk 2 πk 2 P B 1 bk 3 ×   A 2 − b 2 + k 2   arctan  A + k b  − arctan  A − k b  + b  2 k + A ln  ( A − k ) 2 + b 2 ( A + k ) 2 + b 2  + 1 4 g 2 aγ e − db 1 (2 π ) 3 Z 2 π k 2 P B dk 2 bk 3  b d sin ( d ( A − k )) − b d sin ( d ( A + k )) +  A 2 − b 2 + k 2  cosh ( db ) Im ( Ci ( d ( A − ib − k ))) − 2 Ab cosh ( db ) Re ( Ci ( d ( A − ib − k ))) −  A 2 − b 2 + k 2  cosh ( db ) Im ( Ci ( d ( A − ib + k ))) + 2 Ab cosh ( db ) Re ( Ci ( d ( A − ib + k ))) −  A 2 − b 2 + k 2  sinh ( db ) Re ( Si ( d ( A − ib − k ))) − 2 Ab sinh ( db ) Im ( Si ( d ( A − ib − k ))) +  A 2 − b 2 + k 2  sinh ( db ) Re ( Si ( d ( A − ib + k ))) + 2 Ab sinh ( db ) Im ( Si ( d ( A − ib + k )))  , (C.1) where Si ( z ) = R z 0 sin ( t ) /tdt , Ci ( z ) = − R ∞ z cos ( t ) /tdt , arctan ( z ) is the in verse tangent function, Re is the real part and Im is the imaginary part. Here, the first in tegral term corresp onds to the non exp onen tial-decay part, whic h dominates in the strong EBL attenuation regime. The expression of the non exp onen tial-decay part of P γ → γ in the Gaussian magnetic field is (Eq. II I.20 ) P γ → γ = 1 16 g 4 aγ π 2 (2 π ) 6 Z dk k 2 P B ( k ) Z dk ′ k ′ 2 P B ( k ′ ) Z 1 − 1 dt  1 + t 2  Z 1 − 1 dt ′  1 + t ′ 2  e id ( kt + k ′ t ′ ) ( A − ib + kt ) 2 ( A + ib − k ′ t ′ ) 2 = 1 16 g 4 aγ π 2 (2 π ) 6 Z dk k 2 P B ( k ) Z dk ′ k ′ 2 P B ( k ′ ) 1 d 2 k 3 k ′ 3 × ( e − ikd  − A (1 − 2 bd ) + k + i  b +  A 2 − b 2 + k 2  d  A − k − ib + e ikd  A (1 − 2 bd ) + k − i  b +  A 2 − b 2 + k 2  d  A + k − ib + de − iAd  2 b +  A 2 − b 2 + k 2  d + i 2 A (1 − bd )  e − bd ( Ei [( i ( A − k ) + b ) d ] − Ei [( i ( A + k ) + b ) d ]) o × ( e − ik ′ d  A (1 − 2 bd ) + k ′ + i  b +  A 2 − b 2 + k ′ 2  d  A + k ′ + ib + e ik ′ d  A (1 − 2 bd ) − k ′ + i  b +  A 2 − b 2 + k ′ 2  d  − A + k ′ − ib + de iAd  2 b +  A 2 − b 2 + k ′ 2  d − i 2 A (1 − bd )  e − bd  Ei  −  i  A − k ′  − b  d  − Ei  −  i  A + k ′  − b  d  o , (C.2) where Ei ( z ) = − R ∞ − z e − t /tdt . The expression of the non exp onential-deca y part of P γ → γ in the sto c hastic magnetic field with non-Gaussian case 18 2 is (Eq. ( IV.28 )) P NG γ → γ = 4 1 2 1 16 g 4 aγ 2 π Z k 2 dk P ρ ( k ) Z 1 − 1 dt − e idkt (1 − 2 bd − idkt ) (2 b + ikt ) 4 = 4 1 2 1 16 g 4 aγ 2 π Z k 2 dk P ρ ( k )  e − 2 bd d 3 3 k i ( Ei (2 bd − idk ) − Ei (2 bd + idk )) +2 ( − 12 b 2 + 16 b 3 d + 16 b 4 d 2 + k 2 + 4 bdk 2 + 8 b 2 d 2 k 2 + d 2 k 4 ) 3(4 b 2 + k 2 ) 3 cos ( dk ) +2 (8 b 3 − 16 b 4 d − 32 b 5 d 2 − 6 bk 2 − 16 b 3 d 2 k 2 + dk 4 − 2 bd 2 k 4 ) 3 k (4 b 2 + k 2 ) 3 sin ( dk )  (C.3) The expression of the non exp onential-deca y part of P γ → γ in the sto c hastic magnetic field with non-Gaussian case 3 is (Eq. ( IV.34 )) P NG γ → γ = 1 2 9 π 2 g 4 aγ B 4 rms 54 25 g N L Z +1 − 1 dt  1 A 4 + 2 A 2 ( b 2 − k 2 B t 2 ) + ( b 2 + k 2 B t 2 ) 2 + e 2 ik B td ( A 2 + ( b + ik B t ) 2 ) 2 + 1 ( b 2 + ( A + k B t ) 2 ) 2  = 1 2 9 π 2 g 4 aγ B 4 rms 54 25 g N L  1 2 Ab ( A 2 + b 2 ) k B  ( iA + b ) arctanh  k B A + ib  − ( iA − b ) arctanh  k B A − ib  − k B ( − A 2 + b 2 + k 2 B ) cos (2 dk B ) − b ( A 2 + b 2 + k 2 B ) sin (2 dk B ) A 2 k B ( A 2 + b 2 − 2 Ak B + k 2 B )( A 2 + b 2 + 2 Ak B + k 2 B ) + ( − 1 + 2 iAd ) e − 2 bd 4 A 3 k B e 2 iAd Ei ( − 2 id ( A + ib − k B )) + ( − 1 − 2 iAd ) e − 2 bd 4 A 3 k B e − 2 iAd Ei (2 id ( A − ib − k B )) + (1 − 2 iAd ) e − 2 bd 4 A 3 k B e 2 iAd Ei ( − 2 id ( A + ib + k B )) + (1 + 2 iAd ) e − 2 bd 4 A 3 k B e − 2 iAd Ei (2 id ( A − ib + k B )) − 1 2 b 3 k B  b ( A − k B ) b 2 + ( A − k B ) 2 − b ( A + k B ) b 2 + ( A + k B ) 2 + arctan  A − k B b  − arctan  A + k B b  , (C.4) where arctanh ( z ) is the in verse function of tanh ( z ) . [1] Zhen Cao et al. A tera–electron volt afterglo w from a narrow jet in an extremely bright gamma-ra y burst. Scienc e , 380(6652):adg9328, 2023. , doi:10.1126/science.adg9328 . [2] Zhen Cao et al. V ery high-energy gamma-ray emission b ey ond 10 T e V from GRB 221009A. Sci. A dv. , 9(46):adj2778, 2023. arXiv:2310.08845 , doi:10.1126/sciadv.adj2778 . [3] D. D. Dzhappuev et al. Carp et-3 detection of a photonlik e air sho w er with estimated primary energy ab o ve 100 T e V in a spatial and temporal coincidence with GRB 221009A. Phys. R ev. D , 111(10):102005, 2025. , doi:10.1103/PhysRevD.111.102005 . [4] Giorgio Galan ti, Lara Na v a, Marco Roncadelli, F abrizio T av ecchio, and Giacomo Bonnoli. Observ ability of the V ery- High-Energy Emission from GRB 221009A. Phys. R ev. L ett. , 131(25):251001, 2023. , doi:10.1103/ PhysRevLett.131.251001 . [5] Giorgio Galanti and Marco Roncadelli. Is Lorentz inv ariance violation found? 4 2025. . [6] Alessandro De Angelis, Giorgio Galanti, and Marco Roncadelli. Relev ance of axion-like particles for v ery-high-energy astroph ysics. Phys. R ev. D , 84:105030, 2011. [Erratum: Ph ys.Rev.D 87, 109903 (2013)]. , doi:10.1103/ PhysRevD.84.105030 . [7] Csaba Csaki, Nemanja Kalop er, Marco Peloso, and John T erning. Sup er GZK photons from photon axion mixing. JCAP , 05:005, 2003. arXiv:hep- ph/0302030 , doi:10.1088/1475- 7516/2003/05/005 . [8] Alessandro Mirizzi and Daniele Montanino. Sto chastic con versions of T e V photons into axion-like particles in extragalactic magnetic fields. JCAP , 12:004, 2009. , doi:10.1088/1475- 7516/2009/12/004 . [9] Giorgio Galanti and Marco Roncadelli. Extragalactic photon–axion-lik e particle oscillations up to 1000 T e V. JHEAp , 20:1–17, 2018. , doi:10.1016/j.jheap.2018.07.002 . [10] Dieter Horns, Luca Maccione, Man uel Mey er, Alessandro Mirizzi, Daniele Montanino, and Marco Roncadelli. Hardening of T e V gamma sp ectrum of AGNs in galaxy clusters by conv ersions of photons into axion-like particles. Phys. R ev. D , 86:075024, 2012. , doi:10.1103/PhysRevD.86.075024 . [11] Alessandro De Angelis, Marco Roncadelli, and Oriana Mansutti. Evidence for a new light spin-zero b oson from cosmological gamma-ra y propagation? Phys. R ev. D , 76:121301, 2007. , doi:10.1103/PhysRevD.76.121301 . [12] Hai-Jun Li, W ei Chao, Xiu-Hui T an, and Y u-F eng Zhou. Axion effects on gamma-ray sp ectral irregularities. I I: Implications of EBL absorption. Phys. L ett. B , 862:139334, 2025. , doi:10.1016/j.physletb.2025.139334 . [13] Man uel Meyer, Daniele Mon tanino, and Jan Conrad. On detecting oscillations of gamma rays into axion-like particles in turbulen t and coherent magnetic fields. JCAP , 09:003, 2014. , doi:10.1088/1475- 7516/2014/09/003 . 19 [14] Carmelo Ev oli, Matteo Leo, Alessandro Mirizzi, and Daniele Montanino. Reionization during the dark ages from a cosmic axion background. JCAP , 05:006, 2016. , doi:10.1088/1475- 7516/2016/05/006 . [15] Xiao-Jun Bi, Y u Gao, Junguang Guo, Nic k Houston, Tianjun Li, F angzhou Xu, and Xin Zhang. Axion and dark photon limits from Crab Nebula high energy gamma-ra ys. Phys. R ev. D , 103(4):043018, 2021. , doi:10.1103/ PhysRevD.103.043018 . [16] Run-Min Y ao, Xiao-Jun Bi, Jin-W ei W ang, and Peng-F ei Yin. Optical circular p olarization induced by axionlik e particles in blazars. Phys. R ev. D , 107(4):043031, 2023. , doi:10.1103/PhysRevD.107.043031 . [17] Qing-Hong Cao, Zuow ei Liu, and Jun-Chen W ang. Distinct photon-ALP propagation mo des. JCAP , 01:099, 2025. arXiv: 2307.15602 , doi:10.1088/1475- 7516/2025/01/099 . [18] Ariane Dekker, Gonzalo Herrera, and Dimitrios Kantzas. Axion-like particle limits from multi-messenger sources. 6 2025. arXiv:2506.14659 . [19] Christopher Eckner and F rancesca Calore. First constraints on axionlike particles from Galactic sub-P e V gamma rays. Phys. R ev. D , 106(8):083020, 2022. , doi:10.1103/PhysRevD.106.083020 . [20] Pierluca Carenza, Carmelo Ev oli, Maurizio Giannotti, Alessandro Mirizzi, and Daniele Mon tanino. T urbulent axion-photon con versions in the milky-w ay , 2021. . [21] Alessandro Mirizzi, Georg G. Raffelt, and Pasquale D. Serpico. Signatures of Axion-Lik e P articles in the Sp ectra of T e V Gamma-Ra y Sources. Phys. R ev. D , 76:023001, 2007. , doi:10.1103/PhysRevD.76.023001 . [22] M. C. David Marsh, James H. Matthews, Christopher Reynolds, and Pierluca Carenza. F ourier formalism for relativistic axion-photon con v ersion with astrophysical applications. Phys. R ev. D , 105(1):016013, 2022. , doi: 10.1103/PhysRevD.105.016013 . [23] Andrea Addazi, Salv atore Cap ozziello, Qingyu Gan, Gaetano Lam biase, and Rome Saman ta. Excess radiation from axion- photon conv ersion. Phys. R ev. D , 112(6):063025, 2025. , doi:10.1103/ww3t- xtcz . [24] Jamie I. McDonald and Sebastian A. R. Ellis. Resonant conv ersion of gravitational wa v es in neutron star magnetospheres. Phys. R ev. D , 110(10):103003, 2024. , doi:10.1103/PhysRevD.110.103003 . [25] Giorgio Galanti and Marco Roncadelli. Beha vior of axionlike particles in smo othed out domainlike magnetic fields. Phys. R ev. D , 98(4):043018, 2018. , doi:10.1103/PhysRevD.98.043018 . [26] Sergey T roitsky . T ow ards a model of photon-axion conv ersion in the host galaxy of GRB 221009A. JCAP , 01:016, 2024. arXiv:2307.08313 , doi:10.1088/1475- 7516/2024/01/016 . [27] M. Kac helriess and J. Tjemsland. On the origin and the detection of characteristic axion wiggles in photon spectra. JCAP , 01(01):025, 2022. , doi:10.1088/1475- 7516/2022/01/025 . [28] Stephen Angus, Joseph P . Conlon, M. C. David Marsh, Andrew J. P ow ell, and Luk as T. Witko wski. Soft X-ray Excess in the Coma Cluster from a Cosmic Axion Background. JCAP , 09:026, 2014. , doi:10.1088/1475- 7516/ 2014/09/026 . [29] Pierluca Carenza, Ramkishor Sharma, M. C. Da vid Marsh, Axel Branden burg, and Eik e Rav ensburg. Magnetoh ydro dynam- ics predicts heavy-tailed distributions of axion-photon conv ersion. Phys. R ev. D , 108(10):103029, 2023. , doi:10.1103/PhysRevD.108.103029 . [30] Daniele Montanino, F ranco V azza, Alessandro Mirizzi, and Matteo Viel. Enhancing the Sp ectral Hardening of Cosmic T e V Photons b y Mixing with Axionlik e Particles in the Magnetized Cosmic W eb. Phys. R ev. L ett. , 119(10):101101, 2017. arXiv:1703.07314 , doi:10.1103/PhysRevLett.119.101101 . [31] Georg Raffelt and Leo Sto dolsky . Mixing of the photon with lo w-mass particles. Physic al R eview D , 37(5):1237, 1988. [32] Alberto F ranceschini and Giulia Ro dighiero. The extragalactic background light revisited and the cosmic photon-photon opacit y. Astr on. Astr ophys. , 603:A34, 2017. , doi:10.1051/0004- 6361/201629684 . [33] Andrea Addazi, Salv atore Cap ozziello, and Qingyu Gan. Resonant graviton–photon conv ersion with stochastic magnetic field in the expanding universe. Phys. Dark Univ. , 48:101844, 2025. , doi:10.1016/j.dark.2025.101844 . [34] Alexandra Dobrynina, Alexander Karta vtsev, and Georg Raffelt. Photon-photon dispersion of T e V gamma rays and its role for photon-ALP con version. Phys. R ev. D , 91:083003, 2015. [Erratum: Phys.Rev.D 95, 109905 (2017)]. , doi:10.1103/PhysRevD.91.083003 . [35] Grzegorz K ow al, A. Lazarian, and A. Beresny ak. Density Fluctuations in MHD turbulence: Sp ectra, intermittency and top ology. Astr ophys. J. , 658:423–445, 2007. arXiv:astro- ph/0608051 , doi:10.1086/511515 . [36] Iain Brown and Rob ert Crittenden. Non-Gaussianit y from cosmic magnetic fields. Phys. R ev. D , 72:063002, 2005. arXiv: astro- ph/0506570 , doi:10.1103/PhysRevD.72.063002 . [37] P . Jishn u Sai, S. R. Haridev, and Ra jeev Kumar Jain. Inflationary trisp ectrum of gauge fields from scalar and tensor exc hanges. 9 2025. . [38] Da vid Seery , Martin S. Sloth, and Filipp o V ernizzi. Inflationary trisp ectrum from gra viton exchange. JCAP , 03:018, 2009. arXiv:0811.3934 , doi:10.1088/1475- 7516/2009/03/018 . [39] N Aghanim, C Armitage-Caplan, M Arnaud, M Ashdown, F A trio-Barandela, J Aumon t, C Baccigalupi, RB Barreiro, et al. Planc k 2013 results. xxiv. constraints on primordial non-gaussianity . Astr onomy & Astrophysics , 571:A24, 2014. [40] Pranjal T riv edi, Kandaswam y Subramanian, and T. R. Seshadri. Primordial magnetic field limits from the CMB trisp ec- trum: Scalar mo des and Planc k constraints. Phys. R ev. D , 89(4):043523, 2014. , doi:10.1103/PhysRevD. 89.043523 . [41] Alexander D. Dolgo v and Damian Ejlli. Conv ersion of relic gra vitational wa ves into ph otons in cosmological magnetic fields. JCAP , 12:003, 2012. , doi:10.1088/1475- 7516/2012/12/003 . [42] A. Addazi, S. Capozziello and Q. Gan, Exotic P e V atrons as sources of ultra-high-energy gamma ra ys, JHEAp 52 (2026), 100561 doi:10.1016/j.jheap.2026.100561, arXiv:2510.00254 [astro-ph.HE]].