Distance-dependent hardenings in gamma-ray blazar spectra corrected for the absorption on the extragalactic background light
We consider the ensemble of very-high-energy gamma-ray sources observed at distances and energies where a significant absorption of gamma rays is expected due to pair production on the extragalactic background light. Previous studies indicated that spectra of these sources, upon correction for the absorption, exhibit unusual spectral hardenings which happen precisely at the energies where the correction becomes significant. Here, we address this subject with the most recent clean sample of distant gamma-ray blazars, making use of published results of imaging atmospheric Cerenkov telescopes and of the Pass 8 Fermi-LAT data, supplemented by the newest absorption models and individual measurements of sources’ redshifts. We conclude that the distance dependence of spectral hardenings is detected with the statistical significance of 4.5 standard deviations, which suggests that new physical processes, e.g. those involving new hypothetical light bosons (axion-like particles), should be accounted for in the study of the gamma-ray propagation in the Universe.
a,b]Alexander Korochkin, a]Grigory Rubtsov a]and Sergey Troitsky
Prepared for submission to JCAP
Distance-dependent hardenings in gamma-ray blazar spectra corrected for the absorption on the extragalactic background light
Institute for Nuclear Research of the Russian Academy of Sciences,
60th October Anniversary Prospect 7a, Moscow 117312, Russia
APC, Université Paris Diderot, Bâtiment Condorcet, Case 7020,
75205 Paris Cedex 13, France
- 1 Introduction
- 2 The data
- 3 Data analysis and results
- 4.1 Systematic uncertainties
- 4.2 Comparison to previous studies
- 5 Conclusions
- A Observed and best-fit intrinsic spectra of blazars in the main sample.
It is well known since 1960s  that energetic gamma rays produce electron-positron pairs on the extragalactic background light (EBL). Because of this process, gamma-ray flux from distant sources is attenuated severely and the propagation distance is limited. Very-high-energy (VHE) gamma rays (energy GeV) scatter most efficiently on the infrared and visible background photons, and information on the intergalactic absorption may be derived from VHE observations of emitters located at cosmologically large distances.
There exists some long-standing controversy in the determination of the EBL intensity. Direct observations are very difficult to carry out because of strong foreground contamination by the Zodiacal light and by the Galactic foreground. Nevertheless, numerous attempts gave coherent, though quite uncertain, results. At the same time, there exist several theoretical and semi-empirical models for the EBL which, being consistent between each other, point to the intensity a factor of a few lower than that derived from direct observations. Some of the models are based on the galactic counts and represent lower limits on the EBL intensity, given just by the sum of contributions from guaranteed light sources. Others include more theoretical input from the star formation rate history but point to a similar range of intensities. The situation remains unsettled: two independent observations, published in 2017 and using different sophisticated techniques to get rid of the foreground systematics, point consistently to the EBL intensities a factor of two higher than the most recent theoretical model. For the purposes of the present study, we need a conservative (low-absorption) model and use that of Korochkin and Rubtsov (2018)  as the most recent available benchmark for the low-EBL theoretical models. Variations of the absorption model will be discussed in the context of systematic uncertainties of our results.
Distant blazars were observed in energetic gamma rays from the early days of the VHE astronomy. It has been pointed out long ago that some of them are seen from distances for which the optical depth with respect to the pair production on EBL is significant. This apparent tension with expectations was dubbed “the infrared/TeV crisis” . However, subsequent improvements in EBL models demonstrated that the tension is not that strong. Moreover, apparent hardening of the intrinsic spectrum (corrected for the pair-production attenuation, that is “deabsorbed”), seen in several individual sources, might be related to physical conditions in particular objects, though a successful model for these hardenings is missing .
The situation has been changed when the number of sources observed at large optical depths became sufficient for studies of their ensembles. By comparing deabsorbed spectra of two pairs of similar sources, T. Kneiske (unpublished) has pointed out that spectra of distant sources exhibit unusual hardenings at high energies while those of physically similar nearby objects do not. The study of a sample of 7 sources observed by imaging atmospheric Cerenkov telescopes (IACTs) at optical depths by Horns and Meyer  demonstrated that the energy at which such hardenings happen is correlated with the distance to the source in such a way that the spectra become harder only when the correction for the pair-production absorption is significant. This correlation can hardly be caused by any physical reason and suggests that the optical depth is estimated incorrectly. The statistical significance of this “pair-production anomaly” was estimated at the level111Hereafter, we quote statistical significances in terms of standard deviations, , as it is customary in the high-energy physics. One should always understand the following meaning of these estimates: the probability that the observed, or stronger, effect appears as a random fluctuation is the same as it would happen for a normally distributed random quantity deviating from the mean by this number of standard deviations. The underlying probability distribution is never Gaussian in the analyses we discuss..
The most complete sample of sources observed at large optical depths (, 15 objects observed by IACTs and 5 objects observed by Fermi Large Area Telescope, LAT) has been considered  in 2014. There, we first confirmed the existence of significant hardenings in deabsorbed spectra of a number of objects at the energies where the correction for pair production becomes important. We then considered the full sample of objects, including those with and without statistically significant hardenings, and fit their deabsorbed spectra with the broken power-law function, assuming the break at the energy , at which the pair-production optical depth (these energies are different for sources at different redshifts). We found that the break strength, that is the difference between power-law indices below and above , grows with the source redshift. Statistical significance of this distance dependence, which again suggests an unphysical origin of the hardenings and hence an incorrect account of the absorption, was found to be .
All these results are of great importance not only for the gamma-ray astronomy, but also for other fields of physics and astrophysics. As it will be discussed below, overestimation of the absorption can hardly be understood without invoking new physical or astrophysical phenomena, and the most promising explanation of the data requires the existence of new light particles beyond those described by the Standard Model of particle physics.
Since 2014, several important improvements changed the field. Firstly, many new observations of distant VHE blazars have been performed by IACTs. Secondly, the Fermi-LAT team not only accumulated additional years of statistics, very important for faint distant VHE sources, but also issued the new “Pass 8” reconstruction  improving the data processing. Another related novelty is the 3FHL catalog  of hard-spectrum sources making pre-selection of the sample easier and more uniform than before. Thirdly, from the EBL side, new and improved conservative theoretical models have become available [2, 9, 10], but at the same time new sophisticated observational studies strengthened the tension in the determination of the background radiation intensity [11, 12]. Finally, for a number of objects, new information on their redshifts became available. The main purpose of this work is to revisit the claim of Ref.  with the enlarged high-quality sample of sources, taking into account all the listed improvements. In addition, unlike in the letter-like paper of 2014, here we present all details of the analysis so that the readers could repeat the study themselves in case new data are available. We also present some tests demonstrating independence of our results from potential biases and their stability with respect to systematic uncertainties.
The result of the present study supports previous claims of distance-dependent hardenings in spectra of VHE gamma-ray blazars, corrected for the pair-production attenuation with conservative EBL models. Absence of these hardenings is excluded with the statistical significance of .
The rest of the paper is organized as follows. In Sec. 2, we discuss the gamma-ray data used in our study. Section 2.1 summarizes general criteria for the sample construction. Section 2.2 detalizes how they are implemented for the objects observed by IACTs, while Sec. 2.3 discusses the construction of the Fermi-LAT subsample. We list the resulting sample of 28 sources and discuss its general properties in Sec. 2.4.
Section 3 describes the data analysis procedure and presents the main results of the paper. We describe how the correction to the pair-production opacity is implemented and perform a combined analysis of the full sample of 28 sources. We fit all the spectra with absorbed broken power laws fixing the break energy to ; clearly, the value of varies with the distance to the source. Here, we obtain the main result of the paper: a evidence for the distance dependence of spectral hardenings.
This result, confirming our results of 2014 but obtained with the enlarged and cleaned up most recent data set and a new absorption model, is discussed in detail in Sec. 4. There, we first turn, in Sec. 4.1, to potential biases and systematic uncertainties which might affect our result. In Sec. 4.1.1, we discuss variations in the absorption model: besides the most recent model we use throughout the paper, Korochkin and Rubtsov (2018) , we repeat our analysis with three other models for EBL, namely: Gilmore et al. (2012) fixed , which was used in our 2014 paper; Franceschini et al. (2017) ; and a high-absorption model artificially normalized to the most recent direct observational data [11, 12] on the background-light intensity. The statistical significance of the effect reported in Sec. 3 is always above , remaining, for the low-absorption models, consistent with the result for the baseline model we use. Next, in Sec. 4.1.2, we discuss the possible bias related to the fact that only brightest sources are seen at large distances, the Malmquist bias. We extend our sample to include numerous sources not detected in the highest-energy spectral bins, but satisfying at the same time all other selection criteria. We use flux upper limits for the undetected sources and demonstrate that they are consistent with the distance-dependent spectral hardenings found in Sec. 3. In Sec. 4.1.3, we address another potential bias related to the fact that, generally, blazars detected at higher redshifts are mostly flat-spectrum radio quasars (FSRQs) while those at lower redshifts are mostly BL Lac type objects (BLLs). FSRQs might have intrinsic features in the spectrum which mimic the distance dependence of hardenings just because all observed FSRQs are distant. We remove FSRQs from our sample and consider BLLs separately; the effect is clearly seen in the subsample and the reduction of the statistical significance is fully explained by the decrease of the number of objects in the sample. In Sec. 4.1.4, we discuss the effect of uncertainties in the energy determination, which may have important consequences for the estimation of the spectral shape at the highest energies for steeply falling spectra. We perform various tests and conclude that this effect is unlikely to affect our results. Finally, in Sec. 4.1.5 we attempt to estimate a potential impact of possible wrong determination of distances to the sources. Since 2014, new redshift data removed 2 of 20 sources from the original sample of Ref. . We artificially remove 1, 2 or 3 sources from our new sample of 28 objects and demonstrate that this does not change significantly our results, therefore suggesting that (a reasonable number of) erroneous redshifts cannot cause the effect we observe. Therefore, we conclude that none of the biases or uncertainties we are aware of can change the principal result of the paper, and continue in Sec. 4.2 with a brief account of novelties and differences of this result in comparison with previous studies.
While the outline of our study and a summary of results are given in this Introduction, we present a brief account of our conclusions and discuss possible approaches to the interpretation of our results in Sec. 5. Appendix A presents observed and best-fit deabsorbed spectra of 28 blazars entering our main sample.
2 The data
2.1 Selection criteria
The purpose of the present study implies the use of high-quality gamma-ray spectra of sources located at large, confidently known distances. The obvious candidates are blazars with firmly measured redshifts. For various redshifts , the absorption due to pair production becomes important at different energies of photons: for more distant sources, pair production affects the spectra at lower energies222This is true for energies below TeV which we discuss here.. The benchmark energy, at which the absorption becomes important, corresponds to the optical depth and is denoted as hereafter,
which determines for a given absorption model, see Fig. 1.
To study the behaviour of a spectrum in the region of strong absorption, one needs observations at , which requires the use of different instruments for objects located at different redshifts; for large , is typically dozens of GeV, and the relevant instrument is Fermi LAT, while for less distant objects, is of order of several hundred GeV, and the data of IACTs should be used. Though the data are published in quite different ways, it is essential to treat them on equal footing. Here we discuss our general requirements governing selection of the objects, which will be implemented in the next two subsections for the Fermi-LAT and IACT data.
Redshift criteria. Incorrectly determined or uncertain redshifts may hurt any of blazar studies. This happens quite often, especially for BLLs whose spectra do not possess strong emission lines. We require a firm spectroscopic redshift for the objects included into our sample. More specifically, we start from a preselected sample of gamma-ray blazars and check, for every object, the relevant redshift information and references in the NASA/IPAC Extragalactic Database (NED)333Available at http://ned.ipac.caltech.edu.. The redshifts are accepted if they satisfy (R1) and one of (R2A), (R2B), (R2C) criteria:
(R1) the redshift is spectroscopic (not photometric)
EITHER (R2A) the redshift is determined in a dedicated study (if several dedicated studies give different results, the latest one is used),
OR (R2C) the redshift quoted in NED is determined in the Sloan Digital Sky Survey (SDSS), the result is unique and does not change from one release to another.
These criteria are based on the previous experience in the use of redshift data and are of course ad hoc ones. Their application removes many uncertain redshifts.
Spectral criteria. In our work, we use binned gamma-ray spectra, because only this information is publicly available for the objects observed by IACTs. We impose the following criteria for the spectra:
(S1) information for at least 5 energy bins is available;
(S2) the lowest energy of the last spectral bin satisfies AND the lowest energy of the third from below spectral bin is below .
These criteria are aimed at the reconstruction of spectra at both below and above with reasonable accuracy.
Several comments are in order. Firstly, (S2) depends on the assumed absorption model, and therefore the samples we use in Sec. 4.1.1, when studying the sensitivity of our results to the choice of this model, are slightly different from one to another. Secondly, while we require a nonzero measurement in the last bin for our main study, we allow for an upper limit when studying potential effects of the Malmquist bias in Sec. 4.1.2.
2.2 IACT subsample
To construct the sample of blazars observed by IACTs, we start with the database consolidating results published by different observatories, the TeVCat online source catalog 444Available at http://tevcat.uchicago.edu.. The preselected sample includes objects classified as blazars there (classes “HBL”, “‘IBL”, “LBL”, “BL Lac (class uncertain)”, “FSRQ” and “blazar”), 68 objects in total as of summer 2017. Then, for each object in the sample, we checked the redshift selection rules (R1) and ((R2A) or (R2B) or (R2C)), using references from NED, and availability of gamma-ray spectra satisfying (S1), using references from TeVCat.
Some objects have been observed several times, often by multiple instruments. To avoid double counting of information, which may artificially increase or dilute observed effects, we use only one observation for each source. It is chosen on the basis of better statistics (including the number of available spectral bins) and of better coverage of the energy ranges both below and above . Because of strong variability of many blazars, we never combine any two spectra from observations performed at different epochs to enlarge the energy coverage. We also never combine Fermi-LAT spectra with those obtained from IACT observations because of potential systematic differences between the instruments’ energy scales. These requirements represent a major refinement with respect to some of preceding studies.
The next step requires to fix the absorption model. We use the most recent EBL model  as our baseline low-absorption model; when variations of the model are considered (Sec. 4.1.1), this step is repeated. At this step, we calculate the value of for every source and check the condition (S2). We are left with 22 objects which, together with the sources observed by Fermi LAT, are listed in Table 1.
2.3 Fermi-LAT subsample
For blazars observed by Fermi LAT, our starting point is the 3rd Fermi-LAT Catalog of High-Energy Sources, 3FHL 555Available at http://fermi.gsfc.nasa.gov/ssc/data/access/lat/3FHL/ ., presenting a list of sources detected confidently above 10 GeV. From the catalog, we select 1212 objects classified as blazars (classes “BLL”, “bll”, “FSRQ”, “fsrq”, “bcu”). Since the redshifts presented in the catalog may be doubtful or erroneous for distant sources, we check them manually to satisfy (R1) and ((R2A) or (R2B) or (R2C)) criteria, using references from NED, and apply the additional constraint , because for less distant sources, GeV is beyond the limit of reasonable sensitivity of Fermi LAT. This procedure results in a list of 309 blazars. For each of them, we use publicly available Fermi-LAT  data666Available at https://heasarc.gsfc.nasa.gov/FTP/fermi/data/lat/weekly/photon/ . to construct binned spectra satisfying (S1) automatically. To this end, we use Fermi Pass 8  Release 2 data (version 302). We use “SOURCE” class events with version 6 instrumental response functions “P8R2SOURCEV6”, recorded between August 4, 2008, and February 26, 2018 (Mission Elapsed Time interval 239557417 – 541338374). The data are processed with the standard “gtlike” routine from Fermi Science Tools v10r0p5. We impose standard Fermi-LAT cuts for point-source analysis including the Earth zenith angle cut (). The background models used are “glliemv06.fits” for the Galactic component and “isoP8R2SOURCEV6v06.txt” for the isotropic one. For most of the sources (305 of 309), the coordinates of the source are taken from the 3rd Fermi-LAT Source Catalog, 3FGL . However, not all 3FHL sources have 3FGL counterparts, and for the remaining 4 sources, we used 3FHL coordinates.
The observed spectrum was reconstructed for 8 energy bins, each of 03 dex, for GeV. The lower energy limit of 2 GeV was chosen to cut the inverse-Compton peak energy range in order to improve broken power-law fits used in our subsequent analysis.
At the next step, we fix the absorption model and apply the criterion (S2) to the obtained spectra in the same way as it was done for the IACT subsample. A measurement of the flux in a spectral bin is included if the corresponding number of source photons was nonzero (99% CL). Only 6 of 309 blazars satisfy (S2) and enter our final sample. Other sources, however, are used in the Malmquist-bias tests described in Sec. 4.1.2.
2.4 The sample
For IACT observations, references for binned spectra used in this work are given there. We also list classes of the blazars (BLL or FSRQ) as determined in the catalogs used to compile the samples, TeVCat and 3FHL.
represents the distribution of blazars in our main sample in redshift. One may note that, though BLLs dominate at short and moderate distances and FSRQ are in general farther away, both classes cover large ranges of redshifts. We will return to this subject in Sec. 4.1.3.
3 Data analysis and results
The key point of our study is the use of blazar spectra corrected for the pair production in a minimal-absorption model. Since, for IACT observations, only binned spectra are available, we adopt our analysis procedure for this case even for Fermi LAT, in order to process all observations in a uniform way. In contrast with many previous studies, including our work , we do not use bin-by-bin deabsorption because it can introduce systematic biases in the highest-energy bins for the following reasons. The optical depth is a function of the incoming photon energy and the source redshift . It is important to note that the absorption is not uniform within the bin: photons of higher energy are absorbed stronger777See footnote on page 2.. For the energies and distances at which the absorption is significant, this effect may affect strongly the result. It can be taken into account if the shape of the observed spectrum within the bin is known, as we have done in Ref. ; however, in practice, statistical uncertainties in the highest energy bins often make its determination unreliable.
For the present study, we start from the assumption about the intrinsic spectrum of a blazar, for which we use a broken power law with the break at . Recall that this energy depends on the redshift of the source , cf. Fig. 1, and therefore we assume breaks at different energies for different sources. Then we account for the absorption with the selected EBL model, integrate the obtained spectrum over the energy bins and compare resulted bin-by-bin fluxes with the observed data points. In this way, we fit the data with the three parameters of the intrinsic spectrum, that is the overall normalization and two spectral indices. We denote as the difference between spectral indices at and in this best-fit assumed spectrum. By construction, nonzero would indicate the presence of distance-dependent breaks because the assumed break position is redshift-dependent.
demonstrates the main result of this work: the presence of distance-dependent hardenings in deabsorbed spectra. The assumption of results in for 28 degrees of freedom, corresponding to the probability to occur as a result of a random fluctuation. This corresponds to the exclusion of the hypothesis of the absence of distance-dependent breaks at 4.5 standard deviations (recall the footnote at p. 1 for our conventions).
This statistical study excludes the absence of distance-dependent spectral hardenings. Next, we check that this result is not mimicked by hardenings at a certain energy , common for all sources in the sample. To this end, we perform two additional tests. First, we make the break energy a free parameter of the fit, so that it is adjusted independently for every source. Figure 4
compares the fitted break positions with the value of ; we see that they agree well to each other: averaged over the sample, . For a combined fit of all 28 spectra with the hypothesis of breaks at , we obtain d.o.f., while it is 0.73 for the case when is a free parameter. The second test consists of fitting all the spectra with breaks at a uniformly fixed energy (we consider values 100 GeV TeV). Figure 5
presents d.o.f. for this analysis as a function of the assumed break energy , together with the results of the first test. We see that the fits with distance-independent breaks are worse than the fit with breaks at for any fixed assumed break energy.
In this section, we perform a study of possible systematic effects affecting our results. We will see that the effect we observed cannot be explained by any potential bias or systematics we are aware of. Then we compare our result with those of previous studies.
4.1 Systematic uncertainties
4.1.1 Absorption models
The result of the present study, that is the presence of distance-dependent features in deabsorbed spectra precisely at the energies for which the correction for absorption becomes important, may indicate some problems with the absorption model used. For this study, we used the most recent published model by Korochkin and Rubtsov (2018), Ref. . In this section, we repeat our study with several other representative models to see that the effect we found is present independently of the choice of the model. Besides the baseline model of Ref. , we considered the following three models:
(ii) Franceschini et al. (2017) ;
The latter model was obtained by scaling the model (ii) by a multiplication factor fitted to the data of Refs. [11, 12]. These two independent studies used different techniques to distinguish the extragalactic background light from the foreground contribution: Ref.  used spectral templates, which are different for EBL and for foregrounds, while Ref.  benefited from observations of an absorbing cloud to separate foregrounds experimentally. Figure 6
presents 15 data points of Ref.  and a single data point of Ref.  together. All the points, corresponding to different background-light wavelengths, are fitted nicely by the intensity of the model  multiplied by a wavelength-independent factor . For the toy high-absorption model we therefore take the model of Ref.  uniformly upscaled by 3.35. This may be considered as an upper limit on the opacity.
As it has been discussed above, the energies , at which the absorption becomes important, depend on the EBL model, and therefore the samples of blazars selected for the study are different for different models. For the models [10, 13], some blazars from Table 1 do not enter the sample (marked in Table 1). For the toy high-absorption model, the list of 45 sources in the sample is available from the authors by request.
We present the resulting significance of the observed distance dependence of spectral hardenings for various absorption models in Table 2.
|EBL model||Number of objects||Significance,|
|Korochkin and Rubtsov (2018)||28||4.5|
|Gilmore et al. (2012) fixed||22||3.9|
|Franceschini et al. (2017)||23||4.5|
|High EBL normalized to [11, 12]||45||29.7|
We see that for the low-absorption models, the significance remains very similar to that obtained for the model of Ref. , . Clearly, a 3.35 times increase in the opacity, suggested by Refs. [11, 12], had to result in strengthening the effect, and we quantified this consistently within our approach888The fact that the increase of the absorption up to the values suggested by Ref.  sharpens the problems of interpretation of the blazar gamma-ray observations has been illustrated in Ref.  for two particular sources..
4.1.2 The Malmquist bias
Inclusion of sources to our sample is limited by observational capabilities of the instruments through the (S2) criterion of Sec. 2.1. Therefore, sources intrinsically more and more luminous are included at larger distances. In principle, sources may have spectral features correlated with their luminosity, and the flux selection might mimic the observed effect (a particular model for this is not known). In this hypothetical scenarios, weaker objects not detected above would not have distance-dependent spectral hardenings which we observe for the sources included in our sample. This suggests a way to test this possible bias by a study of objects not detected above .
To this end, we consider the same pre-selected 3FHL sample described in Sec. 2.3 but relax the condition (S2) of Sec. 2.1. Note that we need a flux-limited sample for this study and therefore perform it for Fermi-LAT sources only, not including the sources observed only by IACTs. We obtain a sample of 31 blazars with redshifts between 0.354 and 2.534, which are detected significantly in spectral bin containing , but not above. For each of the objects, we determine upper limits on the strength of assumed breaks at for these sources in the way described in Sec. Ref. 3.
The standard Fermi-LAT routine gtlike does not produce reliable upper limits for weak fluxes. For our purposes, we first use the gtsrcprob routine to calculate the weight of each photon: the probability to originate from a given source. The sum of weights is counted and gives the estimated numbers of photons from the source and from backgrounds. The number of source photons is determined as the sum of weights and is typically very small for undetected sources. Since the background photons are considered separately, this number may be treated as the number of signal events with zero background. We round this number to an integer (0 in most cases) and estimate the Poisson-based 95% CL upper limit on the flux in the last bin by division to the exposure calculated with the gtexpcube routine.
The results of the calculation of these upper limits on are presented in Fig. 7.
If distance-dependent hardenings were observed in our sample because of the Malmquist bias, then the upper limits on would take randomly distributed positive and negative values. We see that this is not the case, and therefore our main result cannot be explained by the Malmquist bias.
4.1.3 Source classes
Gamma-ray sources observed at large distances, the blazars, belong to various classes. In this section, we consider a possibility that distant sources are different from nearby ones, and hence may have systematically different spectral features. In general, broadband spectral energy distributions (SEDs) of blazars have two wide bumps, often associated with the synchrotron and inverse-Compton emissions. Relative positions of the bumps are correlated, which may be explained if the same population of relativistic electrons is responsible for both processes. Therefore, the SEDs are, to the first approximation, characterized by the position of the synchrotron peak so that all blazars form “the blazar sequence” . Flat-spectrum radio quasars (FSRQs) have the synchrotron peak frequency, , in the radio band, while blazars with from infrared to X-ray bands represent various subclasses of BL Lac type objects (BLLs). On average, FSRQs are more powerful and less numerous, which results in a potential bias in flux-limited blazar samples: bright FSRQs are seen at large distances while more numerous weaker BLLs dominate the sample at relatively low redshifts. This bias might indeed affect our observations because intrinsic absorption in (distant) FSRQs might introduce spectral features not present in (nearby) BLLs.
To demonstrate that this is not the case, we remove FSRQs from our sample (as it was previously done in several other studies of the gamma-ray absorption, e.g. Ref. ). The classification of FSRQ or BLL given in Table 1 is based on the information given in the 3FHL catalog for Fermi-LAT sources and in TeVCat for sources observed by IACTs. Since the classification is uncertain and various authors may use different criteria for claiming a source is a FSRQ or a BLL, we obtained approximate values of for objects in the sample, based on SEDs available in NED and also on Refs. [44, 45, 46]. We counted objects with Hz as FSRQs and others as BLLs, and found that this was in accordance with the classification given in Table 1 for all objects. We note that, while FSRQs are indeed concentrated at larger distances, BLLs cover a wide range of redshifts and the sample of BLLs only may be used for the analysis described in Sec. 3. The sample of BLLs contains 19 objects. The significance of distance-dependent spectral hardenings for BLLs only is , and this reduction in significance with respect to the main analysis is in accordance with the statistical depletion of the sample. Fig. 8
demonstrates that the same quantity does not depend on the synchrotron peak frequency. The analyses therefore support the conclusion that potential BLL/FSRQ selection effects described in the beginning of this section are not responsible for the effect we observe. Indirectly, this may suggest that the intrinsic absorption in FSRQs is a subdominant effect with respect to the absorption on EBL.
4.1.4 Tail of the falling spectrum
At the energies under discussion, spectra are determined from event-by-event information about observed individual photons. Determination of the photon energies is subject to statistical and systematic uncertainties. Statistical uncertainties are usually taken into account in the spectrum reconstruction procedures; however, they might be underestimated, while systematics may affect the overall energy scale of an instrument. We are interested in the highest energies, where, because of steeply falling spectra, the flux is often estimated from observations of a few photons. In the case of a systematic shift of their energies upwards, or if the statistical uncertainty is larger than expected, overestimation of the energies of these few photons may affect significantly the shape of the spectrum.
We first address potential systematic errors. These are instrument-dependent, and we take advantage of having objects observed by four different instruments in our sample. Table 3
|Instruments||Number of objects||Significance|
|MAGIC HESS VERITAS||22||4.5|
|MAGIC HESS LAT||22||4.5|
|MAGIC VERITAS LAT||19||3.1|
|HESS VERITAS LAT||21||3.5|
presents significances of the observed distance-dependent spectral hardenings in subsamples of objects with data of one of the instruments dropped. The observed significance is stable and agrees well with the expectations from the sample size. Therefore, only a coherent upward systematic error in the energy estimation by the three IACTs and Fermi LAT may result in the effect we found. We note that this situation is unlikely.
We turn now to the possibility of underestimated statistical errors, which may result in occasional overestimation of energies of particular photons. In the case of the published spectra, this effect is normally already taken into account by the observing collaboration by making use of the “unfolding” procedure. In any case, this potential systematics would be strongly suppressed for the objects where the spectral hardening happens at energies for which a sufficient number of photons is observed. To select those from our sample, we replace the (S2) selection criterion from Sec. 2.1 by a much stronger one, requiring two (instead of one) spectral bins above . Only 18 of 28 objects in our sample satisfy these criterion. Performing our analysis for these 18 objects, we find the significance of for this purified sample. This disfavours the relation of our observation to potentially underestimated errors in the energy determination; moreover, the fact that a purified, though depleted, sample exhibits a stronger effect gives a serious support to all our conclusions.
Our study of distance-dependent spectral features is sensitive to correct determination of distances to the sources. However, redshifts of distant blazars are not always determined unambiguously. In particular (see Sec. 4.2), redshifts of 2 of 20 objects in the sample we used for our previous study  in 2014 were found to be unreliable by 2018, and the two objects left our present sample. Though for the present sample we performed a careful selection of objects based on the quality of their redshifts, see criteria (R1), (R2x) in Sec. 2.1, one might imagine that some of the 28 redshifts would still appear wrong in future analyses. To understand how this may affect our results, we artificially removed 1, 2 or 3 objects in all possible combinations from the sample and repeated our analysis with the reduced samples. The lowest significance obtained with 1 object removed is ; with 2 objects removed it is and with three objects removed it is . Figure 9
presents the distribution of significances among samples with 2 objects removed. This simple procedure, which mimics potential impact of would-be incorrect redshifts, demonstrates that the effect we found is not saturated by a few high-redshift sources and would not be ruined if their redshifts are found to be incorrect.
4.2 Comparison to previous studies
In this section, we give a detailed comparison of our methods and results with those of key previous studies of the gamma-ray opacity of the Universe with ensembles of blazar spectra.
4.2.1 Rubtsov and Troitsky (2014) 
The present work uses the same approach as our previous work  and confirms its results. Advantages of the present work, besides the use of new observational data and a more detailed report on systematics, are:
(i) pre-selection of Fermi-LAT sources with the 3FHL catalog and the use of Pass 8 data processing;
(iii) imposing strict criteria on the redshift quality;
(iv) use of a more conservative method which is not based on the bin-by-bin deabsorption.
In particular, these changes resulted in the removal of 4 sources from the sample of Ref. : redshifts of 2 objects (RGB J1448361 and 1ES 0347121) did not satisfy our new quality criteria; the reported value of the redshift of B3 1307433 changed in such a way that the object no longer satisfies the (S2) criterion; for PKS J07301141, new Pass 8 Fermi-LAT reconstruction does not result in a significant detection at , contrary to the Pass 7REP (V15) used in Ref. . This removed the most distant objects in the sample which, together with a more conservative analysis method, resulted in a considerable reduction of significance of the result (was 12.4, now 4.5). The main result remains qualitatively the same.
4.2.2 Horns and Meyer (2012) 
Ref.  was the first study addressing quantitatively spectral hardenings in the ensemble of distant gamma-ray sources.
(1) Both Ref.  and our study address the distance dependence of upward spectral breaks.
(2) Only objects observed by IACTs () were included in the sample of Ref. . They used a stricter requirement of detection at optical depths , resulting in the sample of 7 objects only, while we require and include objects for which the hardening was not significantly detected.
(3) In Ref. , bin-by-bin deabsorption was used but the correction for the mean energy in the bin for deabsorption was not implemented.
4.2.3 Fermi-LAT collaboration (2012) 
The interesting work  presented an analysis of the ensemble of BLLs detected by Fermi LAT at high redshifts and, for the first time, presented a detection of the EBL absorption in stacked spectra of these distant sources. The ensemble of 150 BLLs with redshifts was split into three distance bins, and intrinsic spectra were extrapolated from Fermi-LAT observations at low energies. The EBL absorption correction was subsequently determined and stacked within the three redshift bins.
As we know from our present study, only a few blazars (6 with our criteria) have been confidently observed by Fermi LAT above . Therefore, only the start of the absorption feature was seen in the stacked sample, and the constraints on the amount of absorption were not very tight. For instance, they were expressed as the dependence of for , see Fig. 1 of Ref. , where the width of the allowed band was at 95% (68%) CL, corresponding to the uncertainty of a factor of a few in . While available absorption models are consistent with the upper part of the band, a factor of 2 to 3 lower absorption is allowed for all energies. This is in a good agreement with our results.
4.2.4 Biteau and Williams (2015) 
Ref.  considered 106 spectra of 38 sources observed by IACTs. Every spectrum was deabsorbed bin-by-bin, without applying the correction for the mean energy of the bin, and the fitted by a model spectrum for which a concave shape was assumed. This study did not reveal any anomaly in the absorption on EBL.
We argue that nonobservation of the absorption anomaly in Ref.  does not contradict to the results of the present work. The difference in conclusions is related both to the sample of source spectra used and to the assumptions on the intrinsic spectra. We point out the following differences:
In the sample of Ref. , a different set of sources is used, dominated by nearby ones, while in our study, the anomaly reveals itself stronger for distant sources. The dominance of nearby sources in Ref.  was enhanced by the use of multiple spectra per source, which gave higher statistical weights to better studied nearby objects, see Figure 10.
The individual-source analysis of Ref.  assumed explicitly a concave shape for the deabsorbed spectrum, thus not allowing for the hardenings we study here by construction, see Figure 11. In fact, in most cases (91 of 106 spectra, in particular for all sources with ), the best-fit model corresponded to a power law, that is to the margin of the allowed concave spectra: convex shapes, which give better fits for many of the spectra in our study, were not allowed.
4.2.5 Galanti et al. (2015) 
It is interesting to note that the distance-dependent hardenings can still be noticed with the analysis based on a single power-law fit. Galanti et al. (2015)  considered a sample of 39 blazars observed by IACTs (), for all tested opacities. Deabsorbed spectra were described by single power laws, and the distance dependence of the spectral index was studied. Though this approach does not constructively include spectral breaks, deabsorbed spectra with upward breaks look harder when fitted by a power-law, and this hardness was found to be redshift-dependent. Therefore, results of Ref.  are in agreement with Ref.  and the present work.
In the present work, we readdress the problem of the anomalous transparency of the Universe for high-energy gamma rays, the most significant indication to which  was based on an unphysical distance dependence of spectral hardenings in deabsorbed spectra for an ensemble of blazars. By making use of a robust analysis procedure which avoids bin-by-bin deabsorption, of new gamma-ray and redshift data and of most recent EBL models, we confirm the presence of the anomalous distance dependence with the statistical significance of . Compared to Ref. , the range of the distances spanned by our sample shrinked considerably because we imposed strict criteria on the redshift quality. Together with other changes discussed above, this explains the reduction of the significance with respect to Ref. . We discuss potential biases and systematic errors and, by performig various tests, conclude that they cannot be responsible for the effect we observe.
Interpretations of the anomalous transparency of the Universe have been widely discussed in the literature, see e.g. Ref.  for a brief review and a list of references. The simplest possible culprit for unphysical distance-dependent hardenings might be an overestimated EBL intensity. However, reduction of the intensity below the levels predicted by the models we use might be dangerous because the model intensities are very close to the lower limits from galaxy counting. With the help of the Markov-chain optimization behind the model , we plan to attempt to change the EBL model in such a way that the effect we observe in the present work is reduced and to understand how much in conflict this would be with underlying astrophysical data. This will be a subject of a forthcoming work.
A number of proposals attempted to explain the apparent anomalous transparency of the Universe by adding secondary emission to that arriving directly from the source. If these secondary photons are born relatively close to the observer, they do not have time to produce pairs and arrive unattenuated. With respect to the origin of this secondary emission, it is useful to distinguish the cases when it is caused by electromagnetic cascades  or from interactions of hadronic cosmic-ray particles [51, 52]. The main characteristic feature of both approaches is that they require extremely low magnetic fields all along the propagation path of the cascade, otherwise secondary photons would not point to the sources because of deflections of their progenitor charged particles in the cascade.
Other proposals require modification of the photon interactions. We note that the pair-production cross section has been well determined experimentally, and the only possibility to change it in the kinematic regimes relevant for our study is to allow for small deviations from the Lorentz invariance. However, this change would affect also the development of photon-induced atmospheric showers [53, 54] making the work of IACTs impossible, so that the scenario, in which results of Sec. 3 are explained by the Lorentz-invariance violation, is in fact excluded by the fact that some photons have been detected by these instruments. The remaining explanation invokes a hypothetical axion-like particle (ALP). In external magnetic fields, ALP mixes with photons  and, since ALP does not attenuate on EBL, this makes it possible to detect photons from more distant sources. In one of the scenarios, this mixing takes place in intergalactic magnetic fields along the path from the source to the observer [56, 57], while in the other, a part of photons is converted to ALPs near the source and reconverted back to gamma rays in our neighbourhood (in the magnetic fields of galaxies, clusters and filaments) [58, 59]. Observational consequences of the two scenarios are compared to each other and to data in Ref. . The account of interactions of energetic photons with the axion-like particle allows one to explain consistently all our results without contradicting to any other experimental or observational data.
Appendix A Observed and best-fit intrinsic spectra of blazars in the main sample.
This research has made use of the NASA/IPAC Extragalactic Database (NED), which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. This research has made use of the TeVCat online source catalog  and of the data and tools provided by the Fermi Science Support Center.
We are indebted to Alberto Franceschini for interesting correspondence and for sharing detailed numerical tables of the absorption model  with us; to Manuel Meyer for numerous interesting discussions and helpful comments; to Timur Dzhatdoev, Oleg Kalashev, Mikhail Kuznetsov, Maxim Libanov, Valery Rubakov and Igor Tkachev for interesting discussions. The work of AK and ST on constraining the intensity of EBL and its applications to astrophysical manifestations of axion-like particles was supported by the Russian Science Foundation, grant 18-12-00258.
-  A.I. Nikishov, Absorption of high-energy photons in the Universe, Sov. Phys. JETP 14 (1962) 393 [ZhETF 41 (1962) 549].
-  A. A. Korochkin and G. I. Rubtsov, Constraining the star formation rate with the extragalactic background light, Mon. Not. Roy. Astron. Soc. 481 (2018) 557 [arXiv:1712.06579 [astro-ph.GA]].
-  R. J. Protheroe and H. Meyer, An Infrared background TeV gamma-ray crisis?, Phys. Lett. B 493 (2000) 1 [astro-ph/0005349].
-  F. Aharonian, D. Khangulyan and L. Costamante, Formation of hard VHE gamma-ray spectra of blazars due to internal photon-photon absorption, Mon. Not. Roy. Astron. Soc. 387 (2008) 1206 [arXiv:0801.3198 [astro-ph]].
-  D. Horns and M. Meyer, Indications for a pair-production anomaly from the propagation of VHE gamma-rays, JCAP 1202 (2012) 033 [arXiv:1201.4711 [astro-ph.CO]].
-  G. I. Rubtsov and S. V. Troitsky, Breaks in gamma-ray spectra of distant blazars and transparency of the Universe, JETP Lett. 100 (2014) 355 [Pisma Zh. Eksp. Teor. Fiz. 100 (2014) 397] [arXiv:1406.0239 [astro-ph.HE]].
-  W. Atwood et al. [Fermi-LAT Collaboration], Pass 8: Toward the Full Realization of the Fermi-LAT Scientific Potential, arXiv:1303.3514 [astro-ph.IM].
-  M. Ajello et al. [Fermi-LAT Collaboration], 3FHL: The Third Catalog of Hard Fermi-LAT Sources, Astrophys. J. Suppl. 232 (2017) 18 [arXiv:1702.00664 [astro-ph.HE]].
-  F. W. Stecker, S. T. Scully and M. A. Malkan, An Empirical Determination of the Intergalactic Background Light from UV to FIR Wavelengths Using FIR Deep Galaxy Surveys and the Gamma-ray Opacity of the Universe, Astrophys. J. 827 (2016) 6 [Erratum: Astrophys. J. 863 (2018) 112] [arXiv:1605.01382 [astro-ph.HE]].
-  A. Franceschini and G. Rodighiero, The extragalactic background light revisited and the cosmic photon-photon opacity, Astron. Astrophys. 603 (2017) A34 [arXiv:1705.10256 [astro-ph.HE]].
-  S. Matsuura et al., New Spectral Evidence of an Unaccounted Component of the Near-infrared Extragalactic Background Light from the CIBER, Astrophys. J. 839 (2017) 7 [arXiv:1704.07166 [astro-ph.GA]].
-  K. Mattila, P. Vaisanen, K. Lehtinen, G. von Appen-Schnur, Ch. Leinert, Extragalactic background Light: a measurement at 400 nm using dark cloud shadow II. Spectroscopic separation of dark cloud’s light, and results, Mon. Not. Roy. Astron. Soc. 470 (2017) 2152 [arXiv:1705.10790 [astro-ph.GA]].
-  R. C. Gilmore, R. S. Somerville, J. R. Primack and A. Dominguez, Semi-analytic modeling of the EBL and consequences for extragalactic gamma-ray spectra, Mon. Not. Roy. Astron. Soc. 422 (2012) 3189 [arXiv:1104.0671 [astro-ph.CO]].
-  M. Colless et al. [2DFGRS Collaboration], The 2dF Galaxy Redshift Survey: Spectra and redshifts, Mon. Not. Roy. Astron. Soc. 328 (2001) 1039 [astro-ph/0106498].
-  D. H. Jones et al., The 6dF Galaxy Survey: Final Redshift Release (DR3) and Southern Large-Scale Structures, Mon. Not. Roy. Astron. Soc. 399 (2009) 683 [arXiv:0903.5451 [astro-ph.CO]].
-  S. P. Wakely and D. Horan, TeVCat: An online catalog for Very High Energy Gamma-Ray Astronomy, Proc. 30th International Cosmic Ray Conference, Mexico (2004) 3, 1341. Available at http://tevcat.uchicago.edu.
-  W. B. Atwood et al. [Fermi-LAT Collaboration], The Large Area Telescope on the Fermi Gamma-ray Space Telescope Mission, Astrophys. J. 697 (2009) 1071 [arXiv:0902.1089 [astro-ph.IM]].
-  F. Acero et al. [Fermi-LAT Collaboration], Fermi Large Area Telescope Third Source Catalog, Astrophys. J. Suppl. 218 (2015) 23 [arXiv:1501.02003 [astro-ph.HE]].
-  M. L. Ahnen et al. [MAGIC and Fermi-LAT Collaborations], Very-high-energy gamma-rays from the Universe’s middle age: detection of the blazar PKS 144125 with MAGIC, Astrophys. J. 815 (2015) L23 [arXiv:1512.04435 [astro-ph.GA]].
-  A. Furniss [VERITAS Collaboration], VERITAS Observations of Very High Energy Blazars and Potential for Cosmological Insight, arXiv:1502.03012 [astro-ph.HE].
-  E. Aliu et al. [MAGIC Collaboration], Very-High-Energy Gamma Rays from a Distant Quasar: How Transparent Is the Universe?, Science 320 (2008) 1752 [arXiv:0807.2822 [astro-ph]].
-  J. Aleksic et al. [MAGIC Collaboration], MAGIC discovery of VHE Emission from the FSRQ PKS 122221, Astrophys. J. 730 (2011) L8 [arXiv:1101.4645 [astro-ph.HE]].
-  A. Abramowski et al. [H.E.S.S. Collaboration], H.E.S.S. discovery of VHE gamma-rays from the quasar PKS 1510-089, Astron. Astrophys. 554 (2013) A107 [arXiv:1304.8071 [astro-ph.HE]].
-  S. O’Brien [VERITAS Collaboration], VERITAS detection of VHE emission from the optically bright quasar OJ 287, PoS ICRC 2017 (2018) 650 [arXiv:1708.02160 [astro-ph.HE]].
-  A. S. Madhavan [VERITAS Collaboration], VERITAS Long-Term Observations of Hard Spectrum Blazars, arXiv:1307.7051 [astro-ph.HE].
-  M. L. Ahnen et al., MAGIC observations of the February 2014 flare of 1ES 1011496 and ensuing constraint of the EBL density, Astron. Astrophys. 590 (2016) A24 [arXiv:1602.05239 [astro-ph.HE]].
-  F. Aharonian [H.E.S.S. Collaboration], Detection of VHE gamma-ray emission from the distant blazar 1ES 1101232 with H.E.S.S. and broadband characterisation, Astron. Astrophys. 470 (2007) 475 [arXiv:0705.2946 [astro-ph]].
-  A. Abramowski et al. [H.E.S.S. Collaboration], Multi-wavelength Observations of H 2356309, Astron. Astrophys. 516 (2010) A56 [arXiv:1004.2089 [astro-ph.HE]].
-  E. Aliu et al., A three-year multi-wavelength study of the very-high-energy -ray blazar 1ES 0229200, Astrophys. J. 782 (2014) 13 [arXiv:1312.6592 [astro-ph.HE]].
-  J. Aleksic et al. [MAGIC Collaboration], MAGIC detection of short-term variability of the high-peaked BL Lac object 1ES 0806524, Mon. Not. Roy. Astron. Soc. 451 (2015) 739 [arXiv:1504.06115 [astro-ph.GA]].
-  J. Aleksic et al. [MAGIC Collaboration], Discovery of VHE gamma-rays from the blazar 1ES 1215303 with the MAGIC Telescopes and simultaneous multi-wavelength observations, Astron. Astrophys. 544 (2012) A142 [arXiv:1203.0490 [astro-ph.HE]].
-  V. A. Acciari et al., The Discovery of gamma-Ray Emission From The Blazar RGB J0710591, Astrophys. J. 715 (2010) L49 [arXiv:1005.0041 [astro-ph.HE]].
-  F. Aharonian [H.E.S.S. Collaboration], Simultaneous multiwavelength observations of the second exceptional gamma-ray flare of PKS 2155304 in July 2006, Astron. Astrophys. 502 (2009) 749 [arXiv:0906.2002 [astro-ph.CO]].
-  A. Abramowski et al. [H.E.S.S. Collaboration], H.E.S.S and Fermi-LAT discovery of gamma rays from the blazar 1ES 1312423, Mon. Not. Roy. Astron. Soc. 434 (2013) 1889 5 doi:10.1093/mnras/stt1081 [arXiv:1306.3186 [astro-ph.HE]].
-  V. A. Acciari et al. [AGILE and VERITAS Collaborations], Multiwavelength observations of a TeV-Flare from W Comae, Astrophys. J. 707 (2009) 612 [arXiv:0910.3750 [astro-ph.HE]].
-  F. Aharonian [H.E.S.S. Collaboration], Discovery of VHE gamma-rays from the high-frequency-peaked BL Lac object RGB J0152017, Astron. Astrophys. 481 (2008) L103 [arXiv:0802.4021 [astro-ph]].
-  A. Abramowski et al. [H.E.S.S. Collaboration], Simultaneous multi-wavelength campaign on PKS 2005489 in a high state, Astron. Astrophys. 533 (2011) A110 5 [arXiv:1111.3331 [astro-ph.HE]].
-  F. Aharonian et al. [H.E.S.S. Collaboration], Discovery of VHE gamma-rays from the BL Lac object PKS 0548322, Astron. Astrophys. 521 (2010) A69 [arXiv:1006.5289 [astro-ph.HE]].
-  E. Aliu et al. [VERITAS Collaboration], Multiwavelength observations and modeling of 1ES 1959650 in a low flux state, Astrophys. J. 775 (2013) 3 [arXiv:1307.6772 [astro-ph.HE]].
-  H. Anderhub et al. [MAGIC Collaboration], Simultaneous Multiwavelength observation of Mkn 501 in a low state in 2006, Astrophys. J. 705 (2009) 1624 [arXiv:0910.2093 [astro-ph.HE]].
-  K. Kohri and H. Kodama, Axion-Like Particles and Recent Observations of the Cosmic Infrared Background Radiation, Phys. Rev. D 96 (2017) 051701 [arXiv:1704.05189 [hep-ph]].
-  G. Fossati, L. Maraschi, A. Celotti, A. Comastri and G. Ghisellini, A Unifying view of the spectral energy distributions of blazars, Mon. Not. Roy. Astron. Soc. 299 (1998) 433 [astro-ph/9804103].
-  M. Ackermann et al. [Fermi-LAT Collaboration], The Imprint of The Extragalactic Background Light in the Gamma-Ray Spectra of Blazars, Science 338 (2012) 1190 [arXiv:1211.1671 [astro-ph.CO]].
-  E. T. Meyer, G. Fossati, M. Georganopoulos and M. L. Lister, From the Blazar Sequence to the Blazar Envelope: Revisiting the Relativistic Jet Dichotomy in Radio-loud AGN, Astrophys. J. 740 (2011) 98 [arXiv:1107.5105 [astro-ph.CO]].
-  G. Ghisellini and F. Tavecchio, Fermi/LAT broad emission line blazars, Mon. Not. Roy. Astron. Soc. 448 (2015) 1060 [arXiv:1501.03504 [astro-ph.HE]].
-  M. L. Ahnen et al. [MAGIC Collaboration], Detection of very high energy gamma-ray emission from the gravitationally-lensed blazar QSO B0218+357 with the MAGIC telescopes, Astron. Astrophys. 595 (2016) A98 [arXiv:1609.01095 [astro-ph.HE]].
-  J. Biteau and D. A. Williams, The extragalactic background light, the Hubble constant, and anomalies: conclusions from 20 years of TeV gamma-ray observations, Astrophys. J. 812 (2015) 60 [arXiv:1502.04166 [astro-ph.CO]].
-  G. Galanti, M. Roncadelli, A. De Angelis and G. F. Bignami, Advantages of axion-like particles for the description of very-high-energy blazar spectra, arXiv:1503.04436 [astro-ph.HE].
-  S. V. Troitsky, Axion-like particles and the propagation of gamma rays over astronomical distances, JETP Lett. 105 (2017) 55 [arXiv:1612.01864 [astro-ph.HE]].
-  T. A. Dzhatdoev, E. V. Khalikov, A. P. Kircheva and A. A. Lyukshin, Electromagnetic cascade masquerade: a way to mimic -axion-like particle mixing effects in blazar spectra, Astron. Astrophys. 603 (2017) A59 [arXiv:1609.01013 [astro-ph.HE]].
-  W. Essey and A. Kusenko, A new interpretation of the gamma-ray observations of active galactic nuclei, Astropart. Phys. 33 (2010) 81 [arXiv:0905.1162 [astro-ph.HE]].
-  W. Essey, O. E. Kalashev, A. Kusenko and J. F. Beacom, Secondary photons and neutrinos from cosmic rays produced by distant blazars, Phys. Rev. Lett. 104 (2010) 141102 [arXiv:0912.3976 [astro-ph.HE]].
-  G. Rubtsov, P. Satunin and S. Sibiryakov, Prospective constraints on Lorentz violation from ultrahigh-energy photon detection, Phys. Rev. D 89 (2014) 123011 [arXiv:1312.4368 [astro-ph.HE]].
-  G. Rubtsov, P. Satunin and S. Sibiryakov, Constraints on violation of Lorentz invariance from atmospheric showers initiated by multi-TeV photons, JCAP 1705 (2017) 049 [arXiv:1611.10125 [astro-ph.HE]].
-  G. Raffelt and L. Stodolsky, Mixing of the Photon with Low Mass Particles, Phys. Rev. D 37, 1237 (1988).
-  C. Csaki, N. Kaloper, M. Peloso and J. Terning, Super GZK photons from photon axion mixing, JCAP 0305 (2003) 005 [hep-ph/0302030].
-  A. De Angelis, M. Roncadelli and O. Mansutti, Evidence for a new light spin-zero boson from cosmological gamma-ray propagation?, Phys. Rev. D 76 (2007) 121301 [arXiv:0707.4312 [astro-ph]].
-  M. Simet, D. Hooper and P. D. Serpico, The Milky Way as a Kiloparsec-Scale Axionscope, Phys. Rev. D 77 (2008) 063001 [arXiv:0712.2825 [astro-ph]].
-  M. Fairbairn, T. Rashba and S. V. Troitsky, Photon-axion mixing and ultra-high-energy cosmic rays from BL Lac type objects - Shining light through the Universe, Phys. Rev. D 84 (2011) 125019 [arXiv:0901.4085 [astro-ph.HE]].
-  S. Troitsky, Towards discrimination between galactic and intergalactic axion-photon mixing, Phys. Rev. D 93 (2016) 045014 [arXiv:1507.08640 [astro-ph.HE]].