The local-filament pattern in the anomalous transparency of the Universe for energetic gamma rays

The propagation length of high-energy photons through the Universe is limited by the absorption due to $e^+e^-$ pair production on extragalactic background radiation. Previous studies reported some discrepancies between predicted and observed absorption, suggesting explanations in terms of new physics. However, these effects are dominated by a limited number of observed sources, while many do not show any discrepancy. Here, we consider the distribution in the sky of these apparently anomalous objects, selected in two very different approaches: the study of unphysical hardenings at distance-dependent energies in deabsorbed spectra of TeV blazars, and the observation of ultra-high-energy air showers from the directions of BL Lac type objects. In both cases, directions to the anomalous sources follow the projected local distribution of galaxies, meaning that the distant sources, contributing to the anomalies, are seen through the local filament. This is in line with the proposed earlier explanation of the anomalies based on mixing of photons with axion-like particles in the filament's magnetic field.


I. INTRODUCTION
Since the beginning of very-high-energy (VHE; above ∼ 100 GeV) gamma-ray astronomy in 1990s, it has been clear that very distant active galaxies, blazars, are among the strongest sources in this band. They were detected even with the first modest instruments, despite the predicted [1] attenuation of the flux because of interaction of energetic gamma rays from distant sources with infrared and optical photons constituting the extragalactic background light (EBL), producing e + e − pairs. Though direct observations of EBL suffer from serious uncertainties, models of the EBL intensity and evoltion, see e.g. Refs. [2,3] for reviews, allow one to take the absorption into account. It has been pointed out that gamma-ray spectra of some distant blazars, deabsorbed with these models, are too hard compared to those of similar nearby objects. This problem, in its early formulation, was called "the infrared/TeV crisis" [4] in 2000, and no successful astrophysical solution to it has been proposed. Since the problem was pointed out for a few sources only, and new, more conservative EBL models have been subsequently developed, the crisis attracted a limited attention only.
Further development of the field involved modern atmospheric Cherenkov telescopes which enlarged the amount of blazars detected in VHE to several dozens, opening the way to analyzing statistical samples. Starting from the first sample of 7 objects in Ref. [5], and subsequently with a much larger sample in Ref. [6], the anomalous behavior of a VHE blazar has been identified as a hardening, that is an upward change of the spectral slope, in the spectrum deabsorbed with the ultimate lower-limit EBL model. For blazars located at different distances from the observer D, the absorption on EBL becomes important at different energies E 0 (D), and the hardenings are observed precisely at these energies E 0 , indicating an incorrect account of the absorption.
The observed phenomenon has been termed "anomalous transparency of the Universe".
However, with more and more blazars discovered in VHE, and notably with improving the quality of the measurement of their distances from the Earth, the overall significance of the anomalous transparency seen in VHE blazars became less impressive, see e.g. Refs. [7][8][9] for recent works 1 . While spectra of many newly discovered blazars are consistent with the pair-production attenuation for low-EBL models, several bright objects observed previously continue to demonstrate the anomaly. In addition, new sources with anomalous hardenings have been discovered, see Ref. [7] for details and lists.
On the same timescale of decades, another puzzling effect has been observed and widely discussed, the directional correlation of some of ultra-high-energy (UHE; above ∼ 10 18 eV) cosmic-ray air showers with BL Lac type objects, a subclass of blazars. The correlation was found in Ref. [14] for the published set of cosmic-ray events with the primary energy E > 10 19 eV detected by the High Resolution Fly's Eye (HiRes) fluorescence airshower detector in the stereoscopic mode [15]. Given the deflection of charged cosmic-ray particles in astrophysical magnetic fields, the observed directional coincidence implies neutral primaries. No neutral particle with that high energy can however reach us from the distances at which BL Lacs are located [16]; in particular, photons of E ∼ 10 19 eV produce e + e − pairs on the extragalactic radio background and have the mean free path of a few Megaparsecs. Subsequently, the HiRes collaboration has confirmed [17] the observation and, making use of internal data, has extended it to lower energies, E 10 18 eV.
Unfortunately, the angular resolution of the HiRes experiment in the stereo mode, 0.6 • , remains unsurpassed, and the correlation has not been tested with data of modern fluorescence cosmic-ray detectors yet.
From the theoretical side, explanations of the anomalous transparency effects at VHE require either the suppression of the pair production or the models in which the observed photons do not come from the source but are born much closer to the observer. The former option is possible only in models with the Lorentz-invariance violation which would unfortunately suppress the development of air showers in the atmosphere and therefore make the VHE photons invisible for Cherenkov telescopes [18,19], contrary to observations (see, however, Ref. [20]). The latter option may in principle be realized in two approaches, "cascade" and "conservation". Cascades can develop in the cosmic background fields and radiation, starting either from a VHE photon [21] or from an accompanying UHE proton from the same source [22]; what we detect at the Earth may be secondary particles born at relatively short distances from us. Conservation means that the photon converts near the source to some particle, which does not produce pairs on EBL, then this particle travels unattenuated and reconverts back to a photon close to us. A working example of this mechanism involves oscillations of photons to hypothetical axion-like particles (ALPs) in the external magnetic fields, see Ref. [23] for the general theory, Ref. [24] for an early astrophysical proposal, Refs. [25,26] for two approaches to the problem under discussion and Ref. [27] for a recent short review. Amazingly, the very same ALP-based mechanism can [28] explain in a consistent way the correlations of UHE showers with BL Lacs, which otherwise remain misterious.
The two approaches to the explanation of the anomalous transparency are opposite to each other in the requirements of the presence of ambient magnetic fields along the line of sight. Indeed, while the ALP-photon mixing requires magnetic fields, the cascades get distorted by the fields so that secondary charged particles, and hence photons born in their interactions, no longer point to the original source of the high-energy emission. Therefore, the two scenarios predict opposite patterns of anisotropy in the distribution of the effect strength over the sky: in the cascade case, anomalies are expected to follow low-field regions, while in the ALP case, regions of larger field are favoured. In particular, it has been pointed out that the VHE anomalies may be explained by the ALP-gamma conversion on the Galactic magnetic field and the corresponding anisotropy might be related to the Galactic plane [26]. For UHE photons, the conversion on the Galactic field is suppressed, but local extragalactic structures, small-scale filaments, provide for the required conditions; conversion of VHE photons on the filament field is also possible for certain ALP parameters [28].
Several attempts to find deviations from isotropy in the anomalous transparency have been made. In Ref. [26], the distribution in the sky of a few distant (z > 0.1) VHE blazars, known at that time, was compared to the maps of the photon-ALP conversion probability calculated with three models of the Galactic magnetic field. Visually, a possible correlation with the high-probability regions was pointed out for one of the field models, that of Ref. [29]. No selection of the anomalously hard spectra was done and no statistical study was performed. Subsequently, in Ref. [13], the very same study has been repeated with the same field model and an enlarged sample of VHE blazars, qualitatively confirming the trend. Unfortunately, as it was recognized in Ref. [13], the model of the field is outdated and the test is not perfect because it does not distinguish truly anomalous objects, nor the selection of the sample is isotropic. More recent magneticfield models predict different patterns of the conversion probability in the Milky Way, see e.g. Refs. [30,31], so the anisotropic pattern seen in Refs. [13,26], if real, may be caused by something else. In Ref. [28], it was shown that the cosmic rays correlated to BL Lac type objects do not follow the distribution in the sky of other events, determined by the HiRes exposure, but no particular patterns were tested. In this work, we revisit the question of anisotropy in the anomalous transparency for both VHE and UHE cases and concentrate on the pattern associated with local filaments, as suggested in Ref. [28].

II. ANALYSIS AND RESULTS
We use the weighted density of galaxies along the line of sight as a tracer of the local large-scale structure. The method was used many times in the studies of cosmic-ray anisotropies at UHE, see e.g. Refs. [32][33][34]. The starting point is a flux-limited catalog of galaxies for which we use the 2MASS Redshift Survey (2MRS), Refs. [35,36]. Weights [37,38] are chosen in such a way that nearby galaxies contribute more, so the weighted density is higher for the directions parallel to the local filament in which the Milky Way galaxy resides, a ∼Mpc thick sausage extending from the Virgo cluster to the Fornax cluster, see e.g. Ref. [39]. The weighted density f (l, b) is determined as a function of the direction in the sky given by the Galactic coordinates (l, b) as described in Appendix A. Following Ref. [32], we calculate two sets of values, S = {f (l i , b i )} S for the set S of directions which are associated with the anomalous transparency effects, and B = {f (l i , b i )} B for the control sample B of directions in the sky. The two sets of numbers, S and B, are then compared by means of the Kolmogorov-Smirnov test which gives the probability that they are derived from one and the same distribution.
This method requires to fix the angular smoothing scale when the function f (l, b) is calculated, see Appendix A. While in cosmic-ray studies the scale was determined by expected deflections of charged particles in cosmic magnetic fields, here we do not have such a guidance; we need a particular quantitative hypothesis to fix the scale apriory. While we are interested in the interpretation involving the filament's magnetic fieds, they are poorly measured and we do not know to which extent the number density of galaxies traces the field. Hence we choose to treat the smoothing as a free parameter of the model, see Appendix A, and account for this freedom in a standard statistical approach [40] described in Appendix C.
For the VHE blazars, we start with the data sets of Ref. [7]. There, a sample of blazars with known redshifts z was constructed starting from the TeVCat catalog [41] for the sources observed by Imaging Atmospheric Cherenkov Telescopes (IACTs; 66 sources) and from the 3FHL catalog [42] for those observed with Fermi LAT (307 sources at z > 0.2). Only 26 and 5 of them, respectively, were detected at considerable opacities, and 8 and 2 have hardenings (upward spectral breaks inconsistent with zero at 68% CL) at the appropriate energies in the spectra deabsorbed with the most recent EBL model [43]. These latter 10 objects, see Table I in Appendix B, are considered here as the signal sample S V of objects demonstrating anomalies, while the control sample B V is constructed of 1000 randomly selected sets of 8 of 66 IACT and 2 of 307 Fermi-LAT sources. Note that, because of small fields of view of IACTs, no survey of a significant part of the sky is available at TeV energies. The TeVCat sample of 66 objects is constructed on the base of individual observations while non-detections are normally not published and the sky coverage cannot be quantitatively determined. Contrary, the 3FHL sample is almost flux-limited and covers the full sky. The distribution of objects from the VHE samples in the sky are shown in Fig. 1. The pre-trial Kolmogorov-Smirnov probability that the control and signal samples are derived from one and the same distribution is 0.013 and corresponds to the disk smoothing of 4 • (see Appendix A for details and Fig. 5 in Appendix C for the dependence of the local p-value on the smoothing). The post-trial probability is 0.028.
We now turn to the UHE case. Here, we start with Ref. [14] where a certain number of HiRes air showers with reconstructed primary energies above 10 19 eV where found to correlate with BL Lac type objects. Of 156 sources selected in previous studies, 11 were found to be within 0.8 • from the shower arrival directions, while only 3 were expected assuming isotropy. This angle was determined from the angular resolution of the experiment, 0.6 • , by a Monte-Carlo simulation as maximizing the signal-to-noize ratio, assuming neutral primary particles from BL Lac's. We consider these 11 directions as the signal sample S U , see Table II from Appendix B for the list. The control sample B U is given by arrival directions of all 271 HiRes air showers in the sample [15], see Fig. 2. The pre-trial Kolmogorov-Smirnov p-value for these two samples is 5.2×10 −4 , achieved for the Gaussian smoothing at the 2 • scale. The post-trial probability is 1.1 × 10 −3 .
We now turn to the combination of the two results representing completely independent tests of the anomalous transparency of the Universe at very different energies. We follow the same procedure but compare now two combined distributions by the Kolmogorov-Smirnov test: the signal distribution S V +U consists of directions of all 21 anomalous sources while the control sample B V +U is obtained by taking 1000 random subsamples, each consisting of 8 of 66 TeVCat blazars, 2 of 307 3FHL blazars and 11 of 271 HiRes arrival directions. This results in the pre-trial p 1 = 5.2×10 −4 , again achieved for the Gaussian smoothing at the 2 • scale. The post-trial probability is p = 7.5 × 10 −5 . Were the statistics Gaussian, this would correspond to the significance of 4.0 standard deviations (post trial). Figure 3 shows the distribution of anomalous directions in the sky superimposed on the density plot of the projected weighted galaxy distribution f (l, b).

III. DISCUSSION AND CONCLUSIONS
Recall that of the three lists we used for construction of control samples, the 3FHL catalog and the list of HiRes stereo events represent complete samples (that is, fol-  FIG. 3. Objects with anomalous hardenings in VHE (red circles), HiRes stereo UHE air showers correlated with BL Lac type objects (blue diamonds) and gamma-ray bursts detected in VHE (black triangles) together with the weighted galaxy distribution (density plot) in the celestial sphere (supergalactic coordinates).
low a known sky coverage), while the TeVCat sample is not complete in this sense and may be biased. For the present analysis, however, this possible bias is conservative. Indeed, the anomalous blazars are defined in our VHE sample as those whose deabsorbed spectra harden right at the energy where the deabsorption correction becomes essential. We assume that the true intrinsic spectra of blazars are either power-law or concave since no physical reason for convexity at these energies is known. Suppose now that there is a region in the sky in which direction the Universe is anomalously transparent at high energies, that is the standard deabsorption results in too high fluxes at high energies (hardenings). Then those spectra which are power-law would be reconstructed with upward breaks; this gives our signal sample S V . Intrinsically concave spectra would not exhibit hardening after deabsorption, but these objects would look brighter than physically similar sources observed in other directions in the sky. These stronger high-energy fluxes would make them more probable to enter the sample which consists of detected objects only. Interestingly, this is right what we observe: in Fig. 1, objects from the TeVCat sample without anomalies, represented by black dots, also have a tendency to concentrate towards directions with larger values of f (l, b). In this context it is interesting to note that all three gamma-ray bursts detected at VHE and listed in TeVCat are also seen in the direction of the local filament, see Fig. 3 (we did not study whether they exhibit any kind of anomalies in their spectra). While the tendency is much weaker for the objects without hardenings than the effect we find for S V , it explains the conservative character of the TeVCat selection bias: the control sample contains many "signal" objects. We believe that this is the reason for the p-value for S V to be considerably larger than that for S U . This may be cured when a more complete sample of VHE blazars, incuding non-detections, will be available, for instance, from the Cherenkov Telescope Array [44]. Another option is to develop a more sensitive method for the Fermi-LAT data. This instrument has a small effective area, as compared to IACTs, and presently only five blazars are significantly detected by LAT at the energies where the absorption is important [7]. Other statistical methods than used here should be applied in order to benefit from the Fermi-LAT full sky coverage. In future, the sensitivity in the Fermi-LAT band might be improved with low-threshold high-altitude IACTs [45], e.g. ALEGRO [46].
Further tests of the BL Lac/UHE-shower correlations are necessary before a study of the anisotropy of the effect can be performed with new data. Poor angular resolution of modern UHE cosmic-ray experiments 2 makes it hard to test the correlation found in the HiRes data, though this may be partly compensated [47] by huge statistics accumulated by now. An attempt reported by Auger in 2007 [48] is not conclusive because of low statistics and, more importantly, low sensitivity of the experiment to photon-induced showers [49,50]; see Ref. [28] for a discussion. Large-statistics tests of the BL Lac correlations with the TA data are a prerequisite for further studies of the effect reported here at UHE.
We turn now to the interpretation of the effect in terms of the ALP-photon mixing in filaments, see Ref. [28] for details. The mechanism involves the conversion of a part of emitted photons to ALPs in the magnetic field of a filament containing the source and their reconversion back to photons in the filament containing the observer. The maximal-mixing conditions, see Ref. [28], depend on two ALP parameters, -the mass m and the ALP-photon coupling g, -the photon energy E, the magnetic field B and the size of the field-filled region L. The observed anisotropy suggests that the conversion happens in the directions along the filament and not in the transverse direction, as it is graphically shown in Fig. 6 in Appendix D. This means that 0.5 Mpc L 20 Mpc. The magnetic field in filaments exists [51] but is very poorly known. Computer simulations, see e.g. Ref. [52], indicate that B ∼ 10 −8 G is a reasonable order-of-magnitude estimate for cluster outskirts. It is quite nontrivial to satisfy the maximal-mixing conditions for very different energies, E ∼ 10 12 eV (VHE) and E ∼ 10 19 eV (UHE), simultaneously, especially given that the ALP-photon mixing is suppressed by Quantum-Electrodynamics effects for large values of the product EB, see Ref. [28] for details. It is therefore remarkable that these conditions are indeed satisfied, for both energy bands, in the local filament for g ∼ (a few) × 10 −11 GeV −1 and m ∼ (a few) × 10 −9 eV. More precise estimates of the required ALP parameters are hardly possible because of the lack of a firm magneticfield model. This parameter range is allowed (see e.g. Ref. [27]) by present experimental and astrophysical limits 3 and is a subset of that previously invoked for the explanation of the anomalous transparency of the Universe, see e.g. Ref [13]. At the same time, it is within the reach of experiments of the near future, including solar axion telescopes TASTE [54] and IAXO [55], as well as the laboratory experiment ALPS-IIc [56].

ACKNOWLEDGMENTS
The author is indebted to A. Korochkin, M. Kuznetsov, M. Libanov, A. Plavin, M. Pshirkov, V. Rubakov, G. Rubtsov and P. Tinyakov for illuminating and helpful discussions of various aspects of this study. This work was supported by the Russian Science Foundation, grant 18-12-00258.
Appendix A: Construction of the weighted density of galaxies For the construction of the weighted density of galaxies in the local Universe, we follow the approaches developed and applied previously in the UHECR context, see e.g. Refs. [32][33][34]. We start from the flux-limited (K S < 11.25) complete catalog of galaxies with coordinates and radial velocities, the 2MASS Redshift Survey (2MRS) in its 2019 release, Refs. [35,36]. Like in Ref. [34], we select galaxies with luminosity distances 4 D L between 5 Mpc and 250 Mpc and attribute a weight to each galaxy. This weight is a product of 1/D 2 L and from the lack of irregularities in gamma-ray spectra of sources embedded in the magnetic field of galaxy clusters have been recently shown to suffer from orders-of-magnitude systematic uncertainties [53]. 4 The catalog gives radial velocities which we convert to distances assuming the flat ΛCDM cosmological model with Ω M = 0.308, Ω Λ = 0.692, H 0 = 67.8 km/s/Mpc. an additional factor which accounts for progressive incompleteness of the flux-limited catalog at large distances. This latter factor is calculated with the help of the "sliding-box" method described in detail in Ref. [38], adopting in turn Ref. [37]. Figure 4 presents this factor as a function of D L . In this way, we obtain a list of galaxies with their Galactic coordinates (l i , b i ) in the sky and their weights, . . , N tot , where N tot = 41706 is the total number of galaxies with 5 Mpc≤ D L ≤ 250 Mpc and K S ≤ 11.25 in the 2MRS catalog. The projected weighted density used in our analysis is a function of coordinates (l, b) determined as a sum of these weights over all galaxies in a (smoothed) given direction. We consider two options of smoothing: the disk smoothing, equivalent to simple summing of w i within a cone of the opening angle θ c centered on (l, b), so that the additional weight is where θ is the angle between the directions (l, b) and (l i , b i ); and the Gaussian-like smoothing with the additional weight function Finally, we need to account for the zone of avoidance: the 2MRS catalog is complete for Galactic latitudes |b| ≥ 8 • in the direction to the Galactic center, |l| ≤ 30 • , and for |b| ≥ 5 • for other Galactic longuitudes l. We introduce an additional mask-related weight for every direction (l, b) equal to the inverse integral ofw d orw G over the sky. If this integral is zero, which may happen for smal values of θ c and a narrow band near the Galactic plane, the direction is dropped from the analysis. Various values of θ c were used for this study, see Appendix C.
Appendix B: Coordinates of anomalous sources Tables I and II give the lists of anomalous sources used in the study.

Appendix C: Account of multiple trials
In the present study, we do not assume any particular quantitative model of the anomalous transparency. While our results support the photon-ALP conversion on the magnetic field in the local filament, this field is poorly known. To which extent the number density of galaxies is a tracer of the filament field is also unknown. Quantitatively, this lack of direct relation is parametrized by the angular smoothing we introduce in the calculation of the weighted density. We choose to scan over the smoothing angular scale, 1 • ≤ θ c ≤ 25 • in steps of 1 • , and to try both ways of smoothing described in Appendix A, treating the scan on the unknown parameter as multiple trials in our statistical study. As it is customary for the account of multiple trials, we first calculate the local p-value for each of the 50 variants of smoothing (two functionsw d,G and 25 values of θ c ), see Fig. 5. The minimal p-value, obtained for a certain variant of smoothing, is called the pre-trial p-value, p 1 . The next step is necessary to directly estimate how often this or lower p-value can appear as a fluctuation because of the large number of trials. Were the trials statistically independent, this step would result in the multiplication of p 1 by the number of trials, 50 in our case. However, this is not the case here: different trials correspond to different versions of smoothing of one and the same distribution and hence are strongly statistically interdependent. We use the standard method of treating the multiple comparisons issue for non-independent data, see e.g. Ref. [40]. We generate a large number, N = 10 5 , Monte-Carlo samples which imitate the sample of anomalous sources, and repeat the procedure of calculating the local p-values and taking the minimal of them, p 1,k , for each sample k. The number M of random samples for which p 1,k ≤ p 1 determines the global, or post-trial, pvalue, p ≈ M/N . This p is interpreted as the probability that the distribution of weighted densities of galaxies in the directions to anomalous sources differs from that for all sources due to a random fluctuation for any assumed smoothing function. By definition, this is the probability that determines the significance of the rejection of the null hypothesis of isotropy in our analysis.
Appendix D: Illustration of the interpretation of the main result of the paper Figure 6 illustrates the interpretation of our main result in terms of the photon-axion mixing.