Constraints on axionlike particles from a combined analysis of three flaring $\textit{Fermi}$ flat-spectrum radio quasars

Many theories beyond the Standard Model of particle physics predict the existence of axionlike particles (ALPs) that mix with photons in the presence of a magnetic field. Searching for the effects of ALP-photon mixing in gamma-ray observations of blazars has provided some of the strongest constraints on ALP parameter space so far. Previously, only individual sources have been analysed. We perform a combined analysis on $\textit{Fermi}$ Large Area Telescope data of three bright, flaring flat-spectrum radio quasars, with the blazar jets themselves as the dominant mixing region. For the first time, we include a full treatment of photon-photon dispersion within the jet, and account for the uncertainty in our B-field model by leaving the field strength free in the fitting. Overall, we find no evidence for ALPs, but are able to exclude the ALP parameters $m_a\lesssim200$ neV and $g_{a\gamma}\gtrsim 5 \times 10^{-12}$ GeV$^{-1}$ with 95\% confidence.


I. INTRODUCTION
Axions are very light pseudoscalar particles beyond the Standard Model (SM) of particle physics [1][2][3], which provide a theoretical solution to the strong CP problem [4].Importantly, axions would couple to photons in the presence of an external magnetic field, with a coupling g aγ , proportional to its mass m a [5,6].This coupling would lead to oscillations between photons and axions, comparable to those between neutrino statesan effect that has been the basis for many experimental axion searches (e.g., [7]).So far, none have been found.
Axionlike particles (ALPs) are similar particles in which the m a /g aγ relation is relaxed.Such particles commonly arise in string theories, or as pseudo-Nambu-Goldstone bosons in other SM extensions [8][9][10][11].ALPs would no longer necessarily solve the strong CP problem, but they are good candidates to make up all or some of the dark matter content of the Universe [12][13][14][15].This makes them interesting targets for direct and indirect searches too (e.g, [7,16,17]).In particular, ALP-photon mixing in the various magnetic fields found in space could affect observations of astrophysical sources (e.g., [18,19]).X-and gamma-ray observations of blazars have been used to set some of the strongest constraints on ALP parameter space so far for masses m a 100 neV [20][21][22][23][24][25].
Blazars are active galactic nuclei (AGN) producing jets of relativistic plasma, which are pointed towards us (within a few degrees).This means their emission is strongly enhanced by relativistic effects; blazars make up some of the brightest gamma-ray sources in the sky [26], though they emit across the entire electromagnetic spectrum, from radio to gamma rays.While the detailed emission mechanisms of blazars are still unclear, the low energy emission is usually considered to be synchrotron photons emitted by electrons in the plasma.The high energy emission is then thought to be inverse-Compton (IC) emission from these same electrons up-scattering either their own synchrotron photons (synchrotron self-Compton), or other background photon fields (external Compton) [27].Hadronic models are also possible, for the high energy peak in particular, (e.g., [28,29]), though these models may require super-Eddington jets [30].Significantly for ALP searches, a smooth nonthermal distribution of electrons (as produced by, e.g., shock acceleration [31,32]) would produce intrinsically smooth gamma-ray spectra.The presence of ALPs, however, could produce oscillatory spectral features, as the ALP-photon oscillation length could be energy-dependent for some astrophysical environments along the line of sight to the source [33].Looking for these irregularities in individual blazar spectra (NGC 1275 and PKS 2155-304), using their magnetized cluster environments as the mixing region, has been the basis for constraints with Fermi Large Area Telescope (LAT), High Energy Stereoscopic System, and Chandra observations [20][21][22].These searches require good statistics in the gamma-ray data, which is why bright blazars make good targets-especially when in a flaring state.
Here, we perform a similar search with a combined analysis of Fermi -LAT data for three bright, flaring flatspectrum radio quasars (FSRQs; 3C454.3,CTA 102, and 3C279), using the blazar jets themselves as the main mixing regions.It has been suggested that the strong field in the jet could lead to ALP-photon mixing at higher masses than previously probed by gammaray searches [34][35][36][37][38][39][40], though so far no search has been performed using it as a mixing region.By combining observations from multiple sources, it should also be possible to strengthen the constraints-within the parameter space probed by all the sources, if an ALP signature is seen in one source, it should be seen in the others too.
In Sec.II we outline the data selection and spectral analysis performed on the three sources.Then, in Sec.III we describe the jet and photon-field modeling required to calculate the ALP-photon oscillations produced within the sources.In Sec.IV, we then discuss the fitting and statistical analysis used to compare the ALP and no-ALP hypotheses and place limits on ALP parameter space, before presenting the results of our analysis in Sec.V. Details of the field structure parameters and the spectral energy distribution (SED) modeling are discussed further in Appendices A and B, respectively, and the effects of systematics are discussed in Appendix C.

II. DATA ANALYSIS
The LAT is a pair-conversion, imaging gamma-ray detector on board the Fermi Gamma-ray Space Telescope (Fermi ), which measures gamma rays from 30 MeV up to > 300 GeV energies [41].Our aim is to look for oscillations in Fermi -LAT gamma-ray spectra caused by ALP-photon mixing.We target the three sources with the brightest flares over the Fermi lifetime: 3C454.3,CTA 102, and 3C279 [42].We use FERMIPY v1.0.11[43] and Fermi Science Tools v2.0.8 2 for the analysis.

A. Data selection
Initially, we analyze each source over a significant fraction of the Fermi -LAT lifetime (11.7 years between August 4, 2008 and April 1, 2020) to get an average model for each region of interest (ROI).We choose an energy range of 100 MeV to 500 GeV.This long-term ROI model can then be used as an initial condition for fitting the flare observations.Each ROI is centered on the respective source and has a size of 15 • × 15 • .To avoid including gamma-rays produced from the Earth limb, we only use events with a zenith angle θ z ≤ 90 • .We choose a spatial binning of 0.1 • pixel −1 .We use the P8R3_SOURCE_V2 instrument response functions3 (IRFs) and only use events that pass the P8R3 SOURCE event selection.Because we are looking for spectral oscillations, we make use of the EDISP event classes available with the Pass 8 IRFs [44].Events are classified into four classes, EDISP0 to EDISP3, depending on the quality of their energy reconstruction (worst to best respectively).These classes each contain a similar number of events and are analyzed separately with their corresponding IRFs.This allows us to extract the best spectral information from the data possible.For the long-term analysis, we use eight energy bins per decade.Then, for the flare analyzes, we choose the binning so as to reach the smallest resolvable energy scale.This is done by extracting the detector response matrices for our observations and choosing the bin width to match the minimum ∆E/E for the best energy dispersion class (EDISP3).This gives 65 (3C454.3),67 (CTA 102), and 61 (3C279) bins per decade for our sources.

B. ROI fitting
First, we optimize the ROI model for each of our sources, for all the event types combined, over the entire 11.7 year time range defined above.The initial model includes every point source in the 4FGL catalogue (Data Release 1) [26] and the standard diffuse isotropic and galactic background templates 4 .We then free the normalization of all sources, including the diffuse backgrounds.Point sources within 5 • of the ROI center or with test statistic TS > 10 have the rest of their spectral parameters freed too (TS is the loglikelihood ratio of the likelihoods with and without the source).We free the spectral index of the galactic background as well.These free model parameters are then fitted to the data.Within this fitted ROI, we search for new point sources to add to the model by calculating a TS map.This is done by adding a potential point source (with a power-law index, Γ = 2) at each pixel of the ROI and calculating its TS.Sources with √ TS ≥ 5 are then added to the overall ROI model at the position which gives the highest TS.We then reoptimize the entire ROI, and repeat the process until no more sources are found; overall we find four new sources each for 3C454.3 and 3C279, and three for CTA 102.This gives the final best-fit model for each of our ROIs over the long-term time period.The time ranges used for our flares are taken from the light-curve analysis of [42] and are listed in Table I.For each of these ranges, we redo the above analysis using these final best-fit ROI models, including the new sources, as the initial conditions-this time only freeing the galactic background and sources that still have TS > 10.We fit each event type separately, treating them as separate measurements (as in [21]).Once fitted ROI models have been found for each of our flares, we calculate SEDs for our sources of interest.Following the 4FGL catalogue, the spectra of 3C454.3 and CTA 102 are both best fitted by a power law with a superexponential cutoff: whereas the 3C279 spectrum is best fitted with a logparabola, N is the number of photons received per unit area per unit time at photon energy E, N 0 is the spectral normalization, E 0 is the reference energy, E c is the cutoff energy, and Γ 1 , Γ 2 and κ are indices.Each event type will have different best-fit spectral parameters; those for a combined event-type analysis are shown in Table I, and the corresponding SEDs (E 2 dN/dE) are shown in Fig. 1.For clarity, only every other energy bin is plotted, but we utilize the full energy resolution in our analysis steps.
For each event type k, we then extract likelihood curves5 L k (µ i ), in each energy bin i, as a function of expected counts µ i , from these best-fit SEDs6 (shown as blue bands in Fig. 1).As can be seen, the best statistics are at low energies, and no detected emission is seen at energies above about 80 GeV.These curves can then be used to perform a likelihood ratio test between models with and without ALPs (see Sec. IV).For each event type, the total likelihood for the no-ALPs model is then for each source, where μi are the expected counts from the best-fit spectral models, including all photon absorption (see Fig. 3 below), but without ALPs (best-fits with ALPs will later be denoted with a hat as opposed to a bar).

III. ALP-PHOTON OSCILLATIONS
To test whether an ALP signature is present in the data, we must model the spectral oscillations caused by ALP-photon mixing.
In general, a photon of energy E propagating in a homogeneous field B, (with a component B T , transverse to the photon direction of motion, and parallel to one of the photon polarization states) will oscillate into an ALP with mass m a and coupling g aγ , with a wave number [33,45,46] where m T is the effective mass of the photon (see Ref. [40] for the calculation within the jet).χ and Γ γγ are the total dispersion and absorption terms for the surrounding photon fields respectively, and b = 7α 90π is the vacuum QED term describing dispersion off the magnetic field, with B cr the critical magnetic field B cr = m 2 e /|e| ∼ 4.4 × 10 13 G, where e is the electric charge.Assuming absorption is small, this means there are two so-called "critical energies", around which the oscillation length depends strongly on energy (and so ALP-photon mixing could lead to oscillations in energy spectra): which depends on the effective mass difference between the ALP and the photon, and which depends on the dispersion terms.For astrophysical plasma environments, these energies can be in the gamma-ray energy range for interesting ALP parameters.This has been the basis of previous searches, and is the basis of ours.
In order to model these spectral oscillations, then, we need a model of the field strength and orientation of the magnetic fields along the sight between us and the gamma-ray sources.This allows us to calculate the photon survival probability, P γγ , i.e., the probability that the emitted photon arrives as a photon at Earth as a function of photon energy (taking into account both photon-ALP conversion and absorption via pair production).The magnetic fields we include along the line of sight are the jet field and the galactic magnetic field (GMF) of the Milky Way, as we choose a mass range where the intergalactic magnetic field does not contribute strongly (see Sec. IV below) and these sources are not thought to be in highly magnetized clusters.For the GMF, we use the model of Ref. [47], as used in, e.g.[21,48].We also include extragalactic background light (EBL) absorption for propagation through intergalactic space, using the model of Ref. [49].The dominant mixing region we are using, however, is the jet field.We use the gammaALPs PYTHON package7 to solve the ALP-photon mixing equations throughout-see [52] for an overview.

A. Jet modeling
For mixing within a jet, the detailed structure of the jet field needs to be taken into account, as well as dispersion and absorption from the various photon fields within the jet [39,40,53].
We use the Potter & Cotter jet framework (PC, see [54]) for the overall jet properties (shape of the field strength, bulk Lorentz factors, and electron density) 8 , as discussed in the context of ALP-photon mixing in [40].The structure of the PC jet model is a parabolic, magnetically dominated accelerating jet base, which transitions to a decelerating ballistic conical jet in rough energy equipartition at r tr ∼ 10 5 r g from the black hole, where r g is the gravitational radius, which depends on the black hole mass as r g = 2GM/c 2 .In the accelerating region (r ≤ r tr ) the bulk Lorentz factor is Γ ∝ r 1/2 , and it is Γ ∝ log(r) in the decelerating region (r tr < r ≤ r jet ), where r jet is the jet length.This leptonic jet framework is consistent with theory, observation, and simulations, and is capable of reproducing broadband steady-state SEDs for many blazars [54].
For the location of the gamma-ray emitting regions during the flares r em , we use the lower limits found in Ref. [42], derived from the absence of attenuation due to pair production with broad line region (BLR) photons in the gamma-ray spectra.We use the B(1pc) values found in [54] from fits with the PC model to set the initial value of the magnetic field strength B 0 , which is then left free in the fitting (see Sec. IV below).These initial values are slightly lower than those derived from very-long-baseline interferometry core-shift measurements for each of our sources in Ref. [55], where they assume a conical jet throughout 9 .The electron density varies as n e ∝ R −2 , where R is the jet width, with the value at r tr derived from energetic equipartition.Values for the jet parameters used are listed in Table II.
We then model the detailed field structure as in Ref. [40], with a tangled component (B t ) and a helical component (B h ) that transitions from poloidal to toroidal as r increases down the jet.A constant fraction, f , of the total field energy density is in the tangled component: The radius at which the helical field component becomes toroidal is r T ; the transverse component of the helical field varies as B T ∝ r −α for r < r T .The three  parameters f , r T , and α therefore govern the detailed field structure (along with a treatment of the coherence length of the tangled field).Ideally, these parameters would be allowed to vary in the fit in the same way B 0 is.However, because of computational constraints (P γγ has to be recalculated every time one of them changes), it is necessary to fix them.In Appendix A, we motivate our choices for these parameters from observation and simulations, and show that varying them would be unlikely to strongly affect our final results.For all our sources, we use α = 1 and r T = r tr .Figure 2 shows one example field realization for 3C454.3 with this set of parameters.a Values taken from Ref. [56], unless marked otherwise.b From Ref. [42] unless marked otherwise.c From Ref. [57] unless marked otherwise.For R Hβ for 3C279 we use the relation of R Hβ and RLyα in Ref. [56] to convert from RLyα.d From Ref. [53].

B. Photon fields
As well as the magnetic field structure, background photon fields can also affect ALP-photon mixing [46,53] (see χ and Γ γγ in Eqs. ( 4) and ( 7)).This is because the oscillations are sensitive to slight differences in propagation between the ALP and photon states.Specifically, gamma rays will be affected by photon-photon dispersion and absorption via pair production from background photon fields, whereas ALPs will not.The fields we would expect within FSRQs are those from the central AGN (accretion disk, BLR, dust torus), starlight (extragalactic and from the host galaxy), the cosmic microwave background (CMB), and synchrotron photons from the jet plasma itself.Reference [53] investigated the effects of all these fields on mixing within 3C454.3.They found that for emission regions on the scale of the AGN fields, dispersion off of them will dominate and should be included in the calculations.In particular, for our r em ∼ 0.1 pc, we expect the BLR and torus fields to be the most important, as the disk is only relevant at much smaller scales.Dispersion from the CMB can play a large role within the jet at energies above 100 GeV, but, as can be seen from Fig. 1, we are only interested in lower energies.(In fact, gamma-rays at these energies would likely be absorbed by BLR photons in our sources anyway, see Fig. 3 below).modeling of the starlight and synchrotron fields therefore does not have to be extremely precise.Nevertheless, we include all the photon fields for each of our sources, using the same method and models as Ref. [53].The various parameters we use for the AGN fields, along with their sources, are given in Table III.
The χ and Γ calculations depend on the geometry as well as the photon energies and energy densities of the background fields.The disk is modeled as flat, extending radially in the plane perpendicular to the jet between R in and R out , with each radius between the two emitting at only one energy (as in Ref. [56]).We use R in = 6r g for all sources, the expected inner disk radius for a Schwarzschild black hole.
The BLR is modeled as a series of concentric rings, each corresponding to an emission line, and also emitting at only one energy.The radii and luminosities of the lines can be derived from those of the Hβ line for each source (we use all the lines in the Appendix of [56]).
We use the torus model described in Ref. [53]-an extension of the flat model of Ref. [56] to include an elliptical torus cross section.Each torus emits at a single energy, depending on its temperature, Θ, and the fraction of disk radiation reemitted in each case is assumed to be ξ dt = 0.1 (as in [56]).All tori for our sources are given the same size and shape.They extend radially between R 1 and R 2 and have a height so as to give a covering fraction-the fraction of the sky obscured by the torus from the point of view of the black hole-of f c = 0.6, which is considered typical (e.g., [58]).The cloud number density within the torus decreases with R from the black hole ∝ R −1 for all our sources (see [53,56] for details).
We also include the EBL, starlight, and CMB fields exactly as described in Ref. [53], though, as mentioned above, they are subdominant.In order to model the (also subdominant) synchrotron photon field within the jet, we then follow Ref. [53] in modeling broadband SEDs for each of our sources (in both flaring and steady states) and compare them to observations (see Appendix B).This also enables us to confirm the overall self-consistency of our jet and photon-field models.
Figure 3 shows photon survival probability P γγ as a function of observed energy for each of our sources, displaying the total absorption from all the photon fields (including intergalactic EBL absorption).The absorption rates are calculated as in Ref. [53], and are included in every calculation, both with and without ALPs.We note that, even though B 0 changing in the fit would, in principle, change the synchrotron field, we keep all the fields fixed.This is a good approximation because the synchrotron field hardly affects the dispersion, and never affects the absorption (see [53] and Appendix B).

IV. STATISTICAL METHODS
We are now in a position to compute photon survival probabilities, P γγ (m a , g aγ , B j ) for ALP-photon beams propagated through our sources, where B = (B 0 , f, α, r T ) and j denotes a realization of the random magnetic field.Figure 4 shows an example P γγ (E obs ) for one pair of ALP parameters for 3C454.3,including both dispersion and absorption from the background fields.Our aim is to compare models with ALPs to the observed Fermi data.We follow the statistical methods of, e.g., Refs.[20,21,48] closely.For a random B-field realization j, and spectral parameters θ, the expected counts including ALPs are then where P γγ i denotes the average over energy bin i, which is necessary because P γγ can vary on energy scales much smaller than the bin width.It is worth noting that possible uncertainties in the instrument response functions used within FERMIPY could possibly slightly affect this expression.In Appendix C, we show that the inclusion of an extra shift or smear in the energy reconstruction or dispersion would not greatly affect our overall results.We calculate P γγ at 500 fine energy bins, logarithmically spaced across our energy range, before averaging.For each source, one set (m a , g aγ , B j , θ) corresponds to a likelihood, where L(µ i ) are the likelihood curves extracted from the Fermi SEDs, evaluated at the expected ALP counts, and is the prior on B 0 , which takes the form of a Gaussian.B0 is the initial value used (see Table III) and σ B is the error derived for the magnetic field strength in Ref.
[55] 10 .In each case, σ B is around 20% of B0 .For each field realization, we fit the ALP spectrum to the data by varying B 0 and θ in such a way as to maximize L ALP (m a , g aγ , B j , θ).We use the iminuit PYTHON package for the fitting.Note that every time B 0 is changed in the fit, P γγ has to be recalculated completely (as is the case when changing m a or g aγ ).When B 0 changes, only the overall field strength is affected; the random field orientations and domain lengths remain the same for a given j.As shown in Appendix B, changing the synchrotron photon field has a negligible effect on dispersion, so we can keep it constant as well.Large variations of B 0 , such as removing the field completely, are discouraged by the prior term in the likelihood.Best-fit values of the field strength and spectral parameters are denoted by Bj and θ.
We scan an 8 × 7 logarithmically spaced grid in (m a , g aγ )-space, with m a ∈ [5, 5000] neV and g aγ ∈ [0.1, 10] × 10 −11 GeV −1 .This region is where we might expect mixing in the jets (see Ref. [40]), with the critical energies (Eqs.( 6) and ( 7)) lying around the Fermi energy range.For masses m a < 5 neV, the precise jet length becomes important, as does conversion in the IGMF.Higher couplings are ruled out by experiment [59], and lower couplings would lead to oscillations too small to be detectable.
In order to treat the random field statistically (following, e.g., Refs.[21,48]), for each (m a , g aγ ) we perform the fits for 100 magnetic field realizations, then sort them by L ALP and choose j = 95 corresponding to the 0.95 magnetic field realization quantile.Each point on the ALP grid then corresponds to a likelihood value, L k ALP (m a , g aγ , B95 , θ), for each event type.
For each source, the overall ALP and no-ALP hypotheses can be compared with the test statistic, defined in the standard way for a likelihood ratio test, where L k 0 are the maximum likelihoods for the no-ALP model, found in Eq. ( 3), with best-fit spectral parameters, θ, and the additional hats in the denominator denote the maximum likelihood over the whole ALP grid [60].A high value of TS would mean that the ALP hypothesis is more likely than the no-ALP hypothesis.To quantify how confidently we could reject the no-ALP hypothesis for a given TS, we use Monte Carlo simulations to find the null TS distribution.For this, we use the simulate_roi function in FERMIPY to generate 100 simulated ROIs for each event type (for each of our sources), every time removing the source and injecting a new one from the best-fit spectral model, including photon absorption.The whole analysis is then repeated on each simulated ROI to get a distribution in TS (shown as the solid lines in Fig. 5).The TS thresholds can then be read from these distributions.Specifically, because we only have 100 simulations, we fit gamma distributions to the TS distributions (dashed lines), from which we can take the 0.95 threshold values (dot-dashed lines in Fig. 5).As can be seen, TS > 10.5 is required to reject the no-ALP hypothesis with 95% confidence for 3C454.3;TS > 12.52 for CTA 102; and TS > 8.29 for 3C279.
Regardless of whether we can claim an ALP detection, we would also like to find which ALP parameters a given observation is inconsistent with-i.e., which values of m a and g aγ can be excluded.This can be done with a different test statistic: which compares the best fit at each point (m a , g aγ ) with the overall best fit of the whole grid.A large λ means that the best fit at that point is significantly worse than the best fit overall and so can be rejected by the data.The underlying distribution of λ(m a , g aγ ) can be found in the same way as the TS distribution, except this time an ALP spectrum is injected into the simulated ROIs.This distribution could, in principle, be different for each point in ALP parameter space, as the oscillations do not depend trivially on m a and g aγ .Doing the simulations for every point (m a , g aγ ) is not computationally feasible, however.Therefore we calculate the λ distribution at various points and linearly interpolate between them to get the 95% λ thresholds across the grid, λ 95 (m a , g aγ ).For CTA 102, we use eight points, spread over the region of parameter space where we might expect exclusions; the top panel of Fig. 6 shows λ 95 for CTA 102.The black points show the injected (m a , g aγ ) pairs.Within the region bounded by these points, λ 95 is interpolated, and it can be seen that λ 95 is not constant, but varies smoothly across the grid-generally lowering for decreasing coupling and increasing mass, i.e., for weaker oscillations.Outside the region bounded by the injected points, nearest-neighbor λ 95 values are used.This is generally a conservative estimate, as λ 95 would continue to decrease into the unprobed parameter space beyond the lower right edge of the interpolation region, where it would approach the TS threshold.Because of the smoothness of the CTA 102 λ 95 distribution, we use fewer injected points for 3C454.3 and 3C279-seven and three respectively-to save on computing time.
So far, we have only discussed individual source analyzes.It is possible to combine the likelihoods from the different sources in the same way as those from the different event types within one source (or even different energy bins within one event type).The total ALP and no-ALP likelihoods are just the product of the individual source likelihoods, L = Π s L s , where s indexes the different sources.This means the final TS tot and λ tot (m a , g aγ ) formulae are11 and where the minima are found after the product over the sources is taken.Of course, because the intrinsic parameters of each source are different, the different sources will be capable of probing slightly different regions of parameter space to greater or lesser degrees.It is important that the likelihoods are combined in this way so that each source contributes proportionally to the overall likelihood.The distributions for these two test statistics can be found in the same way as those for the individual sources.Figure 5 also shows the TS tot distribution; a value TS tot > 18.2 would be required to reject the no-ALP hypothesis with 95% confidence for all the sources combined.The lower panel of Fig. 6 shows the λ 95 tot thresholds across the ALP parameter space.The same interpolation method is used as before, and, again, the overall distribution varies smoothly.For those points where either one or both of 3C454.3 and 3C279 is missing injected simulations, we use the sum of the individual λ 95 thresholds.This again is a conservative (i.e., over-) estimate of λ 95 tot , as by definition λ 95 tot ≤ s λ s 95 everywhere (the two are only equal if the minima of the likelihood profiles for all the sources lie at the same point).Only the point at 1000 neV does not have injected simulations for both CTA 102 and 3C454.3,which together should dominate the overall thresholds, so in the relevant region of parameter space this approximation is small.

Figure 5 also shows (dotted vertical lines) the data TS
values for all the sources both individually and in combination.The TS values for 3C454.3 and 3C279 are below their respective TS thresholds: TS 3C454.3= 4.97, and TS 3C279 = 3.88.For CTA 102, however, their is a slight preference (2σ) for the ALP case12 : TS CTA102 = 13.37, which is over the threshold of 12.52.TS = 13.37 is in the 97% quantile of the gamma-function fitted to the CTA 102 TS distribution, but falls to the 91% quantile if the distribution is simply read from the simulations.Also, this local significance of ∼ 2σ for an ALP signal in the CTA 102 data would be further reduced by a trial factor of 3 when considering the fact we looked at three sources.Therefore, this is not a very significant preference for the ALP case, and indeed, it disappears in the combined analysis: TS tot = 16.03.Overall then, we cannot rule out the no-ALP hypothesis, or in other words, we have not found an ALP signal in the data.
Nonetheless, we are able to place limits on the parameters m a and g aγ .Figure 7 shows λ for each of the sources.The white contours enclose regions where λ ≥ λ 95 , and so show the 95% exclusion contours for each individual source.For CTA 102, the regions of 1σ (68%) and 2σ (95%) preference over the null hypothesis are also shown as green dot-dashed and red solid contours respectively.The best-fit point (m a = 100.8neV and g aγ = 4.64 × 10 −12 GeV −1 ) is also plotted as a gold cross, along with its significance (0.97; 2.17σ).As can be seen, because of this slight preference for the ALP case, the 95% exclusions contour from the CTA 102 data extend to high masses and low couplings (which approximates the no-ALP case).
Aside from CTA102, the constraints from 3C454.3 are the strongest, as would be expected from the comparatively good statistics of the 3C454.3observations (see Fig. 1).Constraints from 3C279 data are much weaker than the other sources not only because its statistics are not quite as good, but also because the configuration of the field parameters means that, with B 0 free, good fits are generally able to be found to the data; for 3C279, λ is smaller for much of the region that is excluded by the other sources than it is in the high-mass-low-coupling region.This highlights the importance of leaving the magnetic field strength free in the fitting.
As was shown in Fig. 5, the preference for the ALP case shown in the CTA 102 data disappears in the combined analysis; we would therefore expect 3C454.3 to contribute most strongly to the combined exclusions.Indeed, Fig. 8 shows λ tot for the whole scanned parameter space, and the 95% exclusions (again shown by the white contour) are only marginally better than the 3C454.3exclusions.In particular, the high-masslow-coupling region is not excluded.This highlights the benefits of using a combined analysis to derive robust Figure 9. Overall 95% exclusion contours for each source and for the combined analysis.The black dotted contour shows constraints from magnetic white dwarf radio polarization [62].Black dashed contours show previous gammaray constraints [20,21].The red dot-dash contour shows the CERN Axion Solar Telescope (CAST) experimental constraints [59].The dark matter line is shown in grey dot-dash.
exclusions. Figure 9 shows how our 95% exclusion contours compare with current constraints, shown by the black and red dashed contours.The dark matter line is shown as a grey dot-dashed line, below which ALPs could make up all of dark matter [15].As can be seen from the figure, the combined analysis performed here allows the previous gamma-ray constraints to be extended.Overall, we can exclude the parameter space 5neV m a 200 neV and g aγ 5 × 10 −12 GeV −1 with 95% confidence.

VI. CONCLUSIONS
Searches for ALP signals in the high energy spectra of AGN have provided some of the strongest constraints on ALP m a and g aγ so far [20][21][22].These searches have all been analyzes of individual sources and have generally used the turbulent magnetic field of the host cluster as their main mixing region; similar searches are also planned in the future (e.g, [48]).
Here, we have performed, for the first time, a combined analysis on Fermi -LAT data of three bright, flaring FSRQs (3C454.3,CTA 102, and 3C279), with the blazar jets themselves as the dominant mixing region.These sources were chosen because they displayed the brightest flaring periods over the Fermi lifetime.
We analyze each of the sources using the FERMIPY PYTHON package, first over a significant fraction of the Fermi lifetime to get average ROI models which are then used as initial conditions for detailed SED analysis of the flaring time periods.This enables us to extract likelihood curves from the resulting flare SEDs, which can be used to compare ALP spectral models to the data with a log-likelihood ratio test.
To find the ALP spectra, we model the jets within the PC framework, with a helical and a tangled field component as outlined in Ref. [40].In particular, based on observed polarization fractions of the sources, we use jets with 30% of the magnetic energy density in the tangled component.Also, for the first time, we include a full treatment of photon-photon dispersion within the jet, following Ref.[53].This requires the modeling of the disk, BLR, torus, synchrotron, starlight, CMB and EBL photon fields within the jets.We have performed SED modeling with our combined jet and photon-field models to ensure they are consistent with both each other and observations.These jet models then allow us to compute ALP spectra for each of the sources and fit them to the Fermi observations.We treat both the tangled field component and the errors in the data statistically, by running the analysis for 100 field simulations on 100 simulated Fermi data sets.Also, unlike previous work, we account for the uncertainty in our B-field model by leaving the field strength free in the fits, including a prior term in the likelihood function based on core-shift estimates of the field strength.
To find the underlying distributions of the test statistics used to place limits, we performed the analysis on simulated data with various ALP-spectra injected into it.This was done across the ALP parameter space, enabling a 2D test-statistic threshold to be constructed, as opposed to using a single value everywhere.
In the CTA 102 data, we find a marginal (2σ) preference for ALPs, with the best fit occurring at m a = 100.8neV and g aγ = 4.64 × 10 −12 GeV −1 (below the dark matter line).This slight preference disappears in the combined analysis, however, highlighting the benefits of using multiple sources.Overall then, we find no evidence for ALPs, but are able to exclude the parameter space 5neV m a 200 neV and g aγ 5×10 −12 GeV −1 with 95% confidence.This is an improvement on previous gamma-ray searches in this mass range, though it is almost completely contained within the magnetic white dwarf polarization constraints of [62].Our constraints do not quite reach the dark matter line, but are limited in coupling to similar g aγ values as previous Fermi limits (see [21]), which is to be expected.Nonetheless, we reach lower couplings than those projected to be reached by the future ALPS II experiment in the same mass range [63], and comparable couplings to the projected limits of the future IAXO experiment [64].
Future searches, with CTA for instance (like those outlined in [48]), could likely take advantage of greater instrumental sensitivity to probe lower couplings using this same method, with the blazar jets as the dominant mixing region.It would also be interesting to see how these limits could be extended in the event of another flare as bright as the one from 3C454.3 used here, to fully take advantage of the combined analysis method.butions; differences between all the other distributions look similar.This means that the differences between separate realizations of the tangled field outweigh the differences between specific distributions of tangled coherence lengths.We therefore do not need to model the l c distribution in detail, and can fix it to a uniform distribution with l c < R.
The precise values of r T and α are hard to constrain observationally, but we can test their effects on the oscillations with the same average P γγ method.When α = 1 is fixed, the differences between r T = 0.3 pc and r T = r tr (again for 3C454.3)are below 1.5%.It seems that the differences in the ordered component of the field produced by varying r T are swamped by the differences between separate realizations of the tangled component, and also, B 0 can vary in the fit to compen-Figure 12. Overall 95% exclusion contours for the combined analysis (black), and the same when α = 0.5 (red).The black dotted contour shows constraints from magnetic white dwarf radio polarization [62].Black dashed contours show previous gamma-ray constraints [20,21].The red dot-dash contour shows the CAST experimental constraints [59].The dark matter line is shown in grey dot-dash.
sate for any changes.Therefore it is reasonable to take r T = r tr for each of the three sources.Figure 11 shows δ FW for α = 0.5 and α = 1 with r T = r tr , when B 0 is left free in the fit, again for 3C454.3.The differences in this case can be larger (∼ 10%), but only in a few isolated regions.Note that α has a larger effect than r T because varying it can produce a larger change in the transverse field strength at the emission region.It is, of course, possible that these relatively large percentage differences in the P γγ s will not greatly affect the final results because they do not occur in important regions of parameter space.We perform the analysis using α = 1 throughout.In order to test the effects of changing α on our final results, we perform a single analysis of 3C454.3 and CTA 102 with α = 0.5, calculating λ(m a , g aγ ) in the same way as described in Sec.IV. Figure 12 shows how the total (combined) exclusions would vary in this case, using the same λ 95 tot threshold calculated with the ordinary analysis.As can be seen, the overall 95% exclusions are only slightly changed by changing α to 0.5; in some places the exclusions are slightly better, in some places they are slightly worse.This is because the regions where α can make a large difference to the P γγ s are generally beyond, or at the edge of, our exclusion region.This means that, particularly at lower masses (m a 200 neV), our new exclusions are robust despite the approximations made concerning the magnetic field structure.
Appendix B: Self-consistency of field and jet models In order to check the self-consistency of our jet and background field models, we calculate steady-state and flaring SEDs for each of our sources with the agnpy PYTHON package 13 , and compare them to broadband observations.We follow the same method as Ref. [53].This is not supposed to be a detailed SED-modeling of our sources (indeed, the spectral parameters and magnetic field strength will be left free in the actual fits), but rather a check that our overall source models are reasonable.To calculate these SEDs, we line up spherical plasma blobs down the jet, each with a field strength (within σ B , the errors derived from Ref. [55]), electron density, and bulk Lorentz factor taken from our global jet models (see Table II).Every blob contains a population of electrons with a power-law distribution function 13 https://doi.org/10.5281/zenodo.4687123The synchrotron emission from all these blobs can be calculated.We can also accelerate electrons (by increasing E c and adjusting β) within individual blobs and calculate their synchrotron and inverse-Compton emission using our field models14 to simulate localized gamma-ray emission regions.For the steady-state emission, we accelerate electrons in the blob at r ss = r tr , as expected from the PC framework, where acceleration is due to a permanent large-scale feature of the jet (e.g., a standing shock, though the detailed acceleration mechanism is not modeled here or in the PC framework).
For the flare emission regions, we use a blob located within the jet at r em , with a radius R em = R/η em , where R is the jet width and η em > 1.This roughly simulates, e.g., a reconnection or magnetoluminescence event within the highly-magnetized region of the jet, as opposed to a large-scale change in jet structure for the flares, which is disfavoured because of the small flaring timescales.Table IV shows the parameters used for the various blobs.Figure 13 shows our model SEDs; they are largely consistent with observations in both the flaring and steady-state cases, so we can have confidence in

Figure 1 .
Figure 1.SEDs for our sources during the flares: 3C454.3(top), CTA 102 (middle), 3C279 (bottom).Best-fit spectral parameters (corresponding to the red lines) and time ranges used are listed in Table I.Black points show detections and triangles show 95% upper limits.Shaded regions show the likelihood curves for each energy bin.For clarity, only points for every other energy bin are plotted.

Figure 2 .
Figure 2. One realization of the transverse component of the magnetic field, BT , for 3C454.3,using the parameters f = 0.3, α = 1, and rT = rtr.Vertical dashed line shows rtr.

Figure 3 .
Figure 3. Photon survival probability Pγγ as a function of observed energy for each of our sources, displaying the total absorption for each of our sources, including all photon fields.

Figure 4 .
Figure 4. Example photon survival probability Pγγ for one realization of the 3C454.3field, using f = 0.3, α = 1, and rT = rtr.Twenty more realizations are shown in the background.Dispersion and absorption off of the background photon fields are included.

Figure 5 .
Figure 5. Cumulative distribution functions (CDF) for the TS values for the individual sources (solid), and the total TStot distribution (black).Dashed lines show the best-fit gamma distributions.Dot-dashed vertical lines show the 95% thresholds and dotted vertical lines show the TS values of the data.

Figure 6 .
Figure 6.λ95 thresholds for CTA 102 (top) and all sources (λtot) (bottom).Points show injected (ma,gaγ) points, which are interpolated between.Outside the bounds of the injected points, the nearest neighbor is used for the threshold.

Figure 7 .
Figure 7. λ(ma, gaγ) for 3C454.3(top), CTA 102 (middle), and 3C279 (bottom).White contours show the 95% exclusions, i.e., they enclose regions where λ ≥ λ95.The green dot-dashed and the red solid contours show the 1σ and 2σ preference regions for CTA 102 respectively, and the gold cross shows the location of the best-fit point.
ACKNOWLEDGEMENTS M. M. acknowledges support from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program Grant agreement No. 948689 (AxionDM) and from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy -EXC 2121 "Quantum Universe" -390833306.The Fermi LAT Collaboration acknowledges generous ongoing support from a number of agencies and institutes that have supported both the development and the operation of the LAT as well as scientific data analysis.These include the National Aeronautics and Space Administration and the Department of Energy in the United States, the Commissariat à l'Energie Atomique and the center National de la Recherche Scientifique / Institut National de Physique Nucléaire et de Physique des Particules in France, the Agenzia Spaziale Italiana and the Istituto Nazionale di Fisica Nucleare in Italy, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), High Energy Accelerator Research Organization (KEK) and Japan Aerospace Exploration Agency (JAXA) in Japan, and the K. A. Wallenberg Foundation, the Swedish Research Council and the Swedish National Space Board in Sweden.Additional support for science analysis during the operations phase is gratefully acknowledged from the Istituto Nazionale di Astrofisica in Italy and the center National d'Études Spatiales in France.This work performed in part under DOE Contract DE-AC02-76SF00515.

Table I .
Time ranges (tstart to t end ) and best-fit spectral parameters (see Eqs. 1 and 2) for a combined event-type analysis of the flares.N0 is the spectral normalization, E0 is the reference energy, Ec is the cutoff energy, and Γ1, Γ2 and κ are indices.

Table II .
Jet properties for our sources: rem and rtr are the locations of the emission region and jet-base transition region respectively; rjet is the jet length; B0 is the field strength, ne is the electron density; Γ is the bulk Lorentz factor; and rT , α, and f are the field structure parameters.

Table III .
Parameters used for the AGN fields.From the top: disk luminosity, black hole mass (for setting the gravitational radius rg), inner disk radius, outer disk radius, Hβ line luminosity, Hβ line radius, torus temperature, inner torus radius, outer torus radius, semiminor to semimajor axis ratio of torus cross section.

Table IV .
Parameters used for the blobs down the jet.The steady-state and flaring gamma-ray emission is produced from blobs at r vhe and rem respectively.