Search for point sources of ultra-high energy photons with the Telescope Array surface detector

TL;DR

This study used the Telescope Array surface detector to search for point sources of ultra-high energy photons in the Northern sky, setting new upper limits on photon fluxes at energies above 10^18 eV, with no significant sources detected.

Contribution

First-time setting of photon point-source flux limits for energies above 10^18.5 eV using the Telescope Array data with a multivariate classification method.

Findings

No significant photon point sources found.

Established upper limits on photon fluxes at various energies.

Provided detailed numerical results for each sky direction.

Abstract

The surface detector (SD) of the Telescope Array (TA) experiment allows one to indirectly detect photons with energies of order $10^{18}$ eV and higher and to separate photons from the cosmic-ray background. In this paper we present the results of a blind search for point sources of ultra-high energy (UHE) photons in the Northern sky using the TA SD data. The photon-induced extensive air showers (EAS) are separated from the hadron-induced EAS background by means of a multivariate classifier based upon 16 parameters that characterize the air shower events. No significant evidence for the photon point sources is found. The upper limits are set on the flux of photons from each particular direction in the sky within the TA field of view, according to the experiment's angular resolution for photons. Average 95% C.L. upper limits for the point-source flux of photons with energies greater than $10^{18}$, $10^{18.5}$, $10^{19}$, $10^{19.5}$ and $10^{20}$ eV are $0.094$, $0.029$, $0.010$, $0.0073$ and $0.0058$ km$^{-2}$yr$^{-1}$ respectively. For the energies higher than $10^{18.5}$ eV, the photon point-source limits are set for the first time. Numerical results for each given direction in each energy range are provided as a supplement to this paper.

Tables

Table 1. Table 1: Bias in the reconstruction of the zenith angle and angular resolution for the photon primaries at various energies.

$E_{\gamma }$, eV	$⟨{\theta }_{\mathrm{rec}.}-{\theta }_{\mathrm{true}}⟩$	ang. resolution
$>{10}^{18.0}$	$-{2.25}^{\circ }$	${3.00}^{\circ }$
$>{10}^{18.5}$	$-{2.24}^{\circ }$	${2.92}^{\circ }$
$>{10}^{19.0}$	$-{2.16}^{\circ }$	${2.64}^{\circ }$
$>{10}^{19.5}$	$-{2.06}^{\circ }$	${2.21}^{\circ }$
$>{10}^{20.0}$	$-{1.72}^{\circ }$	${2.06}^{\circ }$

Table 2. Table 2: Point-source photon-flux upper limits and proton backgrounds averaged over all pixels together with the maximum pre-trial significance of the photon excess over proton background.

$E_{\gamma }$, eV	$⟨F_{\gamma }⟩\le$, ${\mathrm{km}}^{-2}{\mathrm{yr}}^{-1}$ (zero bg.)	$⟨F_{\gamma }⟩\le$, ${\mathrm{km}}^{-2}{\mathrm{yr}}^{-1}$ (“real” bg.)	$⟨N_{\mathrm{bg}}⟩$	max. $\gamma$ signif. (pre-trial)
$>{10}^{18.0}$	0.094	0.069	0.49	$2.72\sigma$
$>{10}^{18.5}$	0.029	0.021	0.52	$2.71\sigma$
$>{10}^{19.0}$	0.010	0.0074	0.34	$2.89\sigma$
$>{10}^{19.5}$	0.0071	0.0055	0.10	$2.76\sigma$
$>{10}^{20.0}$	0.0058	0.0045	0.029	$3.43\sigma$

Equations

$$S_b=\sum _i\left[S_i\times {\left(r_i/r_0\right)}^b\right],$$

$$\mathcal{A}=\frac{\sum _{i,\alpha }|s_{i,\alpha }^{upper}-s_{i,\alpha }^{lower}|}{\sum _{i,\alpha }|s_{i,\alpha }^{upper}+s_{i,\alpha }^{lower}|},$$

$$F_{\mathrm{UL}}=\frac{{\mu }_{\mathrm{FC}}(N_{\mathrm{obs}},N_{\mathrm{bg}})}{A_{\mathrm{eff}}}$$

$$A_{\mathrm{eff}}^i=S\cdot T\cdot \overline{\cos {\theta }_i}\frac{N_{MC,\gamma }^i(\xi >{\xi }_0)}{N_{MC,\gamma }^i}$$

$$\overline{cos\theta }=\cos {\lambda }_0\cos \delta \sin {\alpha }_m+{\alpha }_m\sin {\lambda }_0\sin \delta ,$$

$${\alpha }_m=\{\begin{array}{ll}0& ;\zeta >1,\\ \pi & ;\zeta <-1,\\ \arccos \zeta & ;-1<\zeta <1;\end{array}$$

$$\zeta =\frac{(\cos {\theta }_{\text{max}}-\sin {\lambda }_0\sin \delta )}{\cos {\lambda }_0\cos \delta }.$$

Full text

††pagerange: Search for point sources of ultra-high energy photons with the Telescope Array surface detector–Search for point sources of ultra-high energy photons with the Telescope Array surface detector††pubyear: 2015

1 Introduction

Ultra-high energy photons are an important tool for studying the high-energy Universe. A plausible source of EeV-energy photons is provided by ultra-high energy cosmic rays (UHECR) undergoing the Greisen-Zatsepin-Kuzmin process (Greisen, 1966; Zatsepin & Kuzmin, 1966) or pair production process (Blumenthal, 1970) on a cosmic background radiation. In this context, the EeV photons can be a probe of UHECR mass composition as well as of the distribution of their sources (Gelmini et al., 2008; Hooper et al., 2011). At the same time, the possible flux of photons produced by UHE protons in the vicinity of their sources by pion photoproduction or inelastic nuclear collisions would be noticeable only for relatively near sources, as the UHE photons attenuation length is smaller than that of UHE protons (see e.g. Bhattacharjee & Sigl (2000) for a review). There also exists a class of so-called top-down models of UHECR generation that efficiently produce the UHE photons, for instance by the decay of heavy dark-matter particles (Berezinsky et al., 1997; Kuzmin & Rubakov, 1998) or by the radiation from cosmic strings (Berezinsky et al., 1998). The search for the UHE photons was shown to be the most sensitive method of indirect detection of heavy dark matter (Kalashev & Kuznetsov, 2016, 2017; Kuznetsov, 2017; Kachelriess et al., 2018; Alcantara et al., 2019). Another fundamental physics scenario that could be tested with UHE photons (Fairbairn et al., 2011) is the photon mixing with the axion-like particles (Raffelt & Stodolsky, 1988) that could be responsible for the correlation of UHECR events with BL Lac type objects observed by the HiRes experiment (Gorbunov et al., 2004; Abbasi et al., 2006). In most of these scenarios, clustering of photon arrival directions rather than diffuse distribution is expected, therefore point-source searches can be a suitable test for them. Finally, the UHE photons could also be used as a probe for the models of Lorentz-invariance violation (Coleman & Glashow, 1999; Galaverni & Sigl, 2008; Maccione et al., 2010; Rubtsov et al., 2012, 2014).

Telescope Array (Abu-Zayyad et al., 2013c; Tokuno et al., 2012) is the largest cosmic-ray experiment in the Northern Hemisphere. It is located at $39.3^{\circ}$ N, $112.9^{\circ}$ W in Utah, USA. The observatory includes a surface detector array (SD) and 38 fluorescence telescopes grouped into three stations. The SD consists of 507 stations that contain plastic scintillators each of 3 m2 area (SD stations). The stations are placed in the square grid with the 1.2 km spacing and covers the area of $\sim 700\;{\rm km}^{2}$. The TA SD is capable of detecting EAS in the atmosphere caused by cosmic particles of EeV and higher energies. The TA SD operates since May 2008.

A hadron-induced extensive air shower (EAS) significantly differs from an EAS induced by a photon: the depth of the shower maximum $X_{\rm max}$ for a photon shower is larger, a photon shower contains less muons and has more curved front (see (Risse & Homola, 2007) for the review). The TA SD stations are sensitive to both muon and electromagnetic component of the shower and therefore may be triggered by both hadron-induced and photon-induced EAS.

In the present study, we use 9 years of TA SD data for a blind search for point sources of UHE photons. We utilize the statistics of the SD data, which benefits from high duty cycle. The full Monte-Carlo (MC) simulation of proton-induced and photon-induces EAS events allows us to perform the photon search up to the highest accessible energies, $E\gtrsim 10^{20}$ eV. As the main tool for the present photon search we use a multivariate analysis based on a number of the SD parameters that make possible to distinguish between photon and hadron primaries.

While searches for diffuse UHE photons were performed by several EAS experiments, including Haverah Park (Ave et al., 2000), AGASA (Shinozaki et al., 2002; Risse et al., 2005), Yakutsk (Glushkov et al., 2007, 2010; Rubtsov et al., 2006), Pierre Auger (Abraham et al., 2007, 2008a; Bleve, 2016; Aab et al., 2017c) and TA (Abu-Zayyad et al., 2013b; Abbasi et al., 2019a), the search for point sources of photons at ultra-high energies has been done only by the Pierre Auger Observatory (Aab et al., 2014, 2017a). The latter searches were based on the hybrid data and were limited to $10^{17.3}<E<10^{18.5}$ eV energy range. In the present paper we use the TA SD data alone. We perform the searches in five energy ranges, namely $E>10^{18}$, $E>10^{18.5}$, $E>10^{19}$, $E>10^{19.5}$ and $E>10^{20}$ eV. We find no significant evidence of photon point sources in all energy ranges and set the point-source flux upper limits from each direction in TA field of view. The search for unspecified neutral particles was also previously performed by the Telescope Array (Abbasi et al., 2015). The limit on the the neutral particles point-source flux obtained in that work is close to the present photon point-source flux limits.

2 TA SD data and reconstruction

2.1 Data set and Monte-Carlo

The data and Monte-Carlo sets used in this study are the same as in the recent TA search for diffuse photons (Abbasi et al., 2019a). We use the TA SD data set obtained in 9 years of observation, from May 11, 2008 to May 10, 2017. During this period, the duty cycle of the SD was about 95% (Abu-Zayyad et al., 2013a; Matthews, 2018).

Monte-Carlo simulations used in this study reproduce 9 years of TA SD observations, as it was shown in (Matthews, 2018). We simulate separately showers induced by photon and proton primaries for the signal and background estimation respectively1)††1) We justify the proton background assumption in the Sec. 3.1., using the CORSIKA code (Heck et al., 1998). The high energy nuclear interactions are simulated with QGSJET-II-03 model (Ostapchenko, 2006), the low energy nuclear reactions with FLUKA package (Ferrari et al., 2005) and the electromagnetic shower component with EGS4 model (Nelson et al., 1985). The usage of the PRESHOWER package (Homola et al., 2005) that takes into account the splitting of the UHE photon primaries into the Earth’s magnetic field allows us to correctly simulate photon-induced EAS up to the 100 EeV primary energy and higher. The thinning and dethinnig procedures with parameters described in (Stokes et al., 2012) are used to reduce the calculation time.

We simulated 2100 CORSIKA showers for photon primaries and 9800 for proton primaries in $10^{17.5}-10^{20.5}$ eV primary energy range. The power spectrum for CORSIKA photon events is set to $E^{-1}$. The showers from the photon and the proton libraries are processed by the code simulating the real time calibration SD response by means of GEANT4 package (Agostinelli et al., 2003). Each CORSIKA event is thrown to the random locations within the SD area multiple times. For photons, these procedures also include reweighting of the events to the $E^{-2}$ differential spectrum, which is assumed for primary photons in this work. As a result, a set of 57 million photon events with $E^{-2}$ spectrum was obtained. The proton Monte-Carlo set used in this study contains approximately 210 million of events. Details of proton Monte-Carlo simulations are described in Refs. (Abu-Zayyad et al., 2013a; Abu-Zayyad, T. and others, 2014; Matthews, 2018). The format of the Monte-Carlo events is the same as the one used for real events, therefore both data and Monte-Carlo are processed by one and the same reconstruction procedure (Abu-Zayyad, T. and others, 2014) described below.

2.2 Reconstruction

In this paper, the same procedure to reconstruct shower parameters is used as in the previous TA photon searches (Abu-Zayyad et al., 2013b; Abbasi et al., 2019a). Each event real or simulated, is reconstructed by a joint fit of the shower-front geometry and the lateral distribution function (LDF) that allows us to determine the shower parameters, including the arrival direction, the core location, the signal density at the fixed distance from the core and the shower front curvature parameter (see (Abu-Zayyad et al., 2013b) for details).

We apply the following set of the quality cuts for both MC and data events:

zenith angle cut: $0^{\circ}<\theta<60^{\circ}\,$, 2. 2.

the number of stations triggered is 7 or more, 3. 3.

the shower core is inside the array boundary with the distance to the boundary larger than 1200 m, 4. 4.

joint fit quality cut: $\chi^{2}/{\rm d.o.f.}<5$.

We also use an additional cut to eliminate the events induced by lightnings. It was previously found by the TA collaboration that lightning strikes could cause events mimicking EAS events, the so-called terrestrial gamma-ray flashes (TGF) (Abbasi et al., 2017, 2018a). Moreover, as the lightning events are expected to be electromagnetic, they resemble photon-induced showers. Therefore, the rejection of these events is crucial for photon search. To make this rejection, we use the Vaisala lightning database from the U.S. National Lightning Detection Network (NLDN) (Cummins & Murphy, 2009; Nag et al., 2011; NLDN, ). From this database we extract the list of the NLDN lightning events detected within a 15–mile radius circle from the Central Laser Facility of the TA, that contains all the TA SD stations, in a time range from 2008-05-11 to 2017-05-10. The list contains 31622 events grouped in time in such a way that a total of 910 astronomical hours contain one or more lightnings. To clean up all possible lighting-induced events from the data set we remove all the events that occur within 10 minutes time intervals before or after the NLDN lightnings. This cut removes the events known to be related to the TGFs reducing the total exposure only by 0.66% and the total number of data events by 0.77%.

The basic observables such as zenith angle, calculated in the reconstruction procedure together with several additional parameters (see below), are used to distinguish photon and proton events by means of a multivariate analysis. Some of the observables are utilizing the features of the experiment’s SD technical design, such as the double-layered scintillators. The detailed description of these technical parameters is given in (Abu-Zayyad et al., 2013c). The full list of 16 parameters used in the present photon search is the same as in the TA SD search for diffuse photons (Abbasi et al., 2019a) and the TA SD composition study (Abbasi et al., 2019b). These parameters are:

Zenith angle, $\theta$. 2. 2.

Signal density at 800 m from the shower core, $S_{800}$. 3. 3.

Linsley front curvature parameter, $a$ obtained front the fit of the shower front with the AGASA-modified Linsley time delay function (Teshima et al., 1986; Abu-Zayyad et al., 2013b). 4. 4.

Area-over-peak (AoP) of the signal at 1200 m (Abraham et al., 2008b). 5. 5.

AoP slope parameter (Rubtsov & Troitsky, 2015). 6. 6.

Number of stations with Level-0 trigger (Abu-Zayyad et al., 2013c) (triggered stations). 7. 7.

Number of stations excluded from the fit of the shower front due to large contribution to $\chi^{2}$. 8. 8.

$\chi^{2}/d.o.f.$ of the shower front fit. 9. 9.

$S_{b}$ parameter for $b=3$; $S_{b}$ is defined as $b$-th moment of the LDF:

[TABLE]

where $S_{i}$ is the signal of $i$-th station, $r_{i}$ is the distance from the shower core to a given station, $r_{0}=1000$ m. The sum is calculated over all triggered non-saturated stations. The $S_{b}$ is proposed as a composition-sensitive parameter in (Ros et al., 2013). 10. 10.

$S_{b}$ parameter for $b=4.5$. 11. 11.

The sum of signals of all triggered stations of the event. 12. 12.

An average asymmetry of signal at upper and lower layers of the stations defined as:

[TABLE]

where $s^{upper|lower}_{i,\alpha}$ is the FADC value of upper or lower layer of $i$-th station at $\alpha$-th time bin. The sum is calculated over all triggered non-saturated stations over all time bins of the corresponding FADC traces. 13. 13.

Total number of peaks of FADC trace summed over upper and lower layers of all triggered stations of the event. To suppress accidental peaks as a result of FADC noise, we define a peak as a time bin with a signal above 0.2 Vertical equivalent muons (VEM) which is higher than a signal of the 3 preceding and 3 consequent time bins. 14. 14.

Number of peaks for the station with the largest signal. 15. 15.

Total number of peaks present in the upper layer and not in the lower one, summed over all triggered stations of the event. 16. 16.

Total number of peaks present in the lower layer and not in the upper one, summed over all triggered stations of the event.

For each MC and data event we also define the “photon energy” parameter $E_{\gamma}$ which is the expected energy of the primary particle assuming it is a photon. This energy parameter is calculated as the function of the zenith angle and the $S_{800}$ parameter from the photon MC simulations (Abu-Zayyad et al., 2013b). For proton MC events, as well as for the majority of data events, the $E_{\gamma}$ parameter is not the actual primary energy but merely a parameter needed for the consistent comparison of proton events and possible photon events. It is important to note that for majority of proton-induced events, the reconstructed $E_{\gamma}$ parameter is systematically higher than that of photon-induced events of the same primary Monte-Carlo energy. For instance, at $\sim 10$ EeV Monte-Carlo energy the mean $E_{\gamma}$ for protons is $\sim 40\%$ higher than that for photons, if we assume the averaging over zenith angle. Due to this fact the proton background for SD photon search is higher with respect to the hypothetical ideal situation when the energy reconstruction bias is independent of the primary particle type. All the energy values considered in this work is assumed to be $E_{\gamma}$ values, unless the other meaning is specified.

The reconstructed values of shower zenith angle, $\theta_{\rm rec}$, for photon primaries are systematically underestimated. The possible reason for this is the azimuthal asymmetry of the shower front, that originates from the fact that the shower arrives younger to the front-side stations and older to the back-side ones. The reconstruction bias is defined as a deviation of the event $\theta_{\rm rec}$ from a real Monte-Carlo zenith angle of this event, $\theta_{\rm true}$. The average values of this bias for various energies $E_{\gamma}$ are given in Table 1. In this study we correct both proton and photon Monte-Carlo events and data events by these average bias values. This correction allows us to restore arrival directions of possible photon-induced events more accurately, while not affecting the background of hadron-induced events, which is known to be highly isotropic (Deligny et al., 2017). Another crucial parameter for the point-source search is the angular resolution of the experiment. It is defined as a $0.68$ percentile of a distribution of Monte-Carlo events over opening angle between event reconstructed arrival direction and real Monte-Carlo arrival direction. The angular resolution of TA SD for proton primaries at “proton energy”, $E_{\rm p}=10^{19}$ eV, was estimated to be $1.5^{\circ}$ (Abu-Zayyad et al., 2012). As it was mentioned above, in the present study we use the reconstruction of (Abu-Zayyad et al., 2013b) for both data and Monte-Carlo events. Using the photon Monte-Carlo set, after applying the zenith-angle bias correction described above, we estimate the angular resolution for photon primaries at various energies $E_{\gamma}$. The results are shown in Table 1.

3 Analysis

3.1 Multivariate analysis

The analysis method used in this study to distinguish between photon and proton events is a boosted decision tree (BDT) classifier built with the 16 observable parameters listed in the previous section. As an implementation of this method, we use the AdaBoost algorithm (Freund & Schapire, 1997) from the TMVA package (Hocker et al., 2007) for ROOT (Brun & Rademakers, 1997), in the same way as in the recent TA studies (Abbasi et al., 2019a, b).

The BDT is trained to separate proton MC events from photon MC events. Both proton and photon MC sets are split into three parts with equal amount of events in each: one for training the classifier, the second one for testing the classifier and the last one for the calculation of proton background and photon effective exposure, respectively. We train the classifier separately in five photon energy ranges: $E_{\gamma}>10^{18}$ eV, $E_{\gamma}>10^{18.5}$ eV, $E_{\gamma}>10^{19}$ eV, $E_{\gamma}>10^{19.5}$ eV and $E_{\gamma}>10^{20}$ eV. As a result of the BDT procedure, the single multivariate analysis (MVA) parameter $\xi$ is assigned to each MC and data event. $\xi$ is defined to take values in the range $-1<\xi<1$, where proton-induced events tend to have negative $\xi$ values, and photon-induced events — positive $\xi$ values. The resulting $\xi$ distributions of the MC events from the testing sets and the data events for all considered energy ranges are shown in Fig. 1.

From Fig. 1, which shows the distributions of data and Monte-Carlo irrespective of the direction in the sky, one can observe no deviation from the proton distribution in the expected photon signal region. However, possible excesses in one or several separate directions in the sky could be overlooked if we analyze the all-sky averaged $\xi$ distribution. Hereafter we discuss method to set photon–flux upper limit and to search for photon excesses from separate directions on the sky and present respective results.

It is important to note that at primary energies of order EeV and higher there is a potential systematic uncertainty in the estimation of the hadron background for the photon signal. The bulk of the events are induced by protons and/or nuclei, but their mass composition is not known precisely (Aab et al., 2017b; Abbasi et al., 2018b, 2019b). We have examined $\xi$ distributions of the iron nucleus-induced events and found that in average the iron-induced events are less ”photon-like” than proton-induced events. The results of TA work (Abbasi et al., 2019b), where the similar BDT-classifier was used, implies that a mixed nuclei $\xi$ distribution would also deviate from the photon $\xi$ distribution stronger than the proton $\xi$ distribution. Therefore, the assumption of the proton background for the photon search is conservative. However, we also perform independent photon search assuming more realistic mixed nuclei background inferred form the TA SD data in our study (Abbasi et al., 2019b).

3.2 Photon-flux upper limit

In general, the flux upper limit for the particular type of primaries is defined as:

[TABLE]

where $N_{\rm obs}$ is the number of detected events of a given type in a given energy range, $N_{\rm bg}$ is the estimated number of background events in the same energy range, $\mu_{\rm FC}$ is the upper bound of the respective Poisson mean for the given confidence level, defined according to  (Feldman & Cousins, 1998), and $A_{\rm eff}$ is the effective exposure of the experiment for the given type of primaries in the same energy range.

In the present upper-limit calculation we assume the “null hypothesis”, i.e. that there is actually no photons and any excess counts from the expected background, $N_{\rm obs}-N_{\rm bg}$, is considered as a fluctuation of background.

We consider two options of the background estimation. First one is $N_{\rm bg}=0$, this assumption is conservative since for the fixed $N_{\rm obs}$ the upper-limit value is higher for a lower value of $N_{\rm bg}$. Second one is a “real” background of mixed nuclei with the mean $\ln A$ following the one derived from the same TA SD data with the same MVA method in our work (Abbasi et al., 2019b). This background is estimated by down-scaling of the proton background to the respective mean $\ln A$, linearly with $\ln A$, taking into account the recaling of SD energy scale used in work (Abbasi et al., 2019b) to the $E_{\gamma}$ energy scale used in this work.

The separation between photon and proton primaries is defined by a cut on MVA-variable $\xi$. The cut is set at some value $\xi_{0}$ so that any proton with $\xi>\xi_{0}$ is considering as a photon candidate and any photon with $\xi>\xi_{0}$ is contributing to the effective exposure.

To find the minimum value of $F_{\rm UL}^{\gamma}$ as a function of $\xi_{0}$ we optimize the cut position assuming: $N_{\rm obs}=N_{\rm p}(\xi>\xi_{0})$,

where $N_{\rm p}(\xi>\xi_{0})$ is the number of protons passing the $\xi$-cut. As one can see from Fig. 1, the number of MC photon events passing the $\xi$-cut is decreasing with the growth of $\xi_{0}$ leading to the respective decrease of the exposure $A_{\rm eff}^{\gamma}$, also the number of photon candidates, $N_{\rm obs}=N_{\rm p}(\xi>\xi_{0})$, is decreasing, but $N_{\rm obs}=0$ yields a constant non-zero value of $\mu_{\rm FC}$. This implies that there indeed should be a non-trivial minimum value of $F_{\rm UL}^{\gamma}$ as a function of $\xi_{0}$.

For the $\xi$-cut optimization we use the proton Monte-Carlo set normalized to the size of the data set and the photon Monte-Carlo set, the latter is used for the calculation of the photon effective exposure $A_{\rm eff}^{\gamma}$. It is important to note that the optimization procedure tends to place $\xi_{0}$ at the right edge of the proton distribution in Fig. 1, therefore the number of candidates $N_{\rm obs}$ and the upper-limit value $F_{\rm UL}$ are subject to fluctuations. These fluctuations become apparent when one considers upper limits for particular directions in the sky with small number of events.

Up to this moment the procedures of upper-limit calculation and cut optimization were similar to those used to search for diffuse photons (Abbasi et al., 2019a). The difference of the analysis procedure used here from that of (Abbasi et al., 2019a) is in the usage of the separate event sets for different directions in the sky. The $\xi$-cut is also optimized separately for every direction studied. We pixelize the sky in equatorial coordinates $\{\alpha,\delta\}$ using the HEALPix package (Gorski et al., 2005) into 12288 pixels ($N_{side}=32$). For the pixel “i” with the center $\{\alpha_{i},\delta_{i}\}$ the corresponding data set contains events located inside a spherical cap region around the pixel center within an angular distance that equals to the experiment’s angular resolution at the respective energy (see Tab. 1)1)††1) The distance between any pixel centers is smaller that experiment’s angular resolution at all considered energies, therefore the experiment FOV is overlapped without gaps, but some events in adjacent pixels could be the same..

The effective exposure of the experiment to photons at the pixel “i” is given by:

[TABLE]

where $S$ is the area of the experiment, $T$ is the period of observation, $\theta_{i}$ is the zenith angle at which the pixel “i” is seen by the experiment, $N^{i}_{MC,\gamma}$ is the total number of photon events simulated in the respective pixel and $N^{i}_{MC,\gamma}(\xi>\xi_{0})$ is the number of these events that pass the $\xi$-cut. The same pixel in equatorial coordinates is seen by the experiment at different $\theta$ depending on time, therefore the diurnal mean value $\overline{cos{\theta}}$ is used. It is given by the expression (Sommers, 2001):

[TABLE]

where $\delta$ is the declination, $\lambda_{0}$ is the geographical latitude of the experiment, $\theta_{max}$ is the maximum zenith angle of the events considered in the particular analysis and $\alpha_{m}$ is given by the expression

[TABLE]

where

[TABLE]

The “effective” part of the exposure, $\frac{N^{i}_{MC,\gamma}(\xi>\xi_{0})}{N^{i}_{MC,\gamma}}$, is calculated using the photon Monte-Carlo. To have enough statistics for this calculation, one needs to generate separate Monte-Carlo sets for each sky-map pixel. However, it is technically unreasonable, since the exposure depends only on declination of the given pixel. We use the following method to increase the Monte-Carlo statistics in each pixel: $\xi_{0}$ is optimized over the events belonging to the whole constant–declination band whose width is twice the angular resolution centered in the given pixel. This method resembles the so-called “scrambling technique” (Cassiday et al., 1990), which was used for instance in the Pierre Auger Observatory search for photon point sources (Aab et al., 2014). The additional advantage of the used method is the preservation of relatively large effective statistics of the Monte-Carlo events in each pixel, including the variety over $\xi$–parameter. We have found that in this case, fluctuations of the $\xi_{0}$ position between adjacent pixels are smaller, compared to the standard scrambling technique. It is reasonable to smooth these fluctuations even further by making a least-squares fit of a $\xi_{0}$ position as a function of declination with a smooth function, for which we use a second–order polynomial. As it has been mentioned before, the flux upper limit remains conservative after this operation. The examples of $\xi_{0}$ as a function of declination and its smooth fitting are shown in Fig. 22)††2) The $\xi_{0}$ points for adjacent declination bands are clustered because these bands are overlapping with each other and a part of their MC events is one and the same..

As the $\xi_{0}$ position for the pixel “i” is fixed, the actual upper-limit value is calculated using the definition (3) with $N^{i}_{\rm obs}=N^{i}_{\rm data}(\xi>\xi_{0})$, where $N^{i}_{\rm data}$ is a number of data events belonging to the respective pixel and $N^{i}_{\rm bg}=0$ or $N^{i}_{\rm bg}=N^{i}_{\ln A}(\xi>\xi_{0})$, where $N^{i}_{\ln A}$ is a number of “real” background events.

The Telescope Array field of view for the considered zenith angle cut ($0^{\circ}<\theta<60^{\circ}$) spans from $-20.7^{\circ}$ to $90^{\circ}$ in declination. However, the event statistics is low in the constant–declination bands near the edges of this interval. Therefore we reduce the considered sky region to $-15.7^{\circ}\leq\delta\leq 85^{\circ}$. It contains 7848 pixels.

3.3 Results

The 95% C.L. photon-flux upper limits calculated in zero background assumption for each pixel in the Telescope Array field of view and for various photon energies are shown in Fig. 3. The numerical values of these limits as well as the limits calculated with the “real” background assumption are given in the supplementary material. The values of the limits averaged over all pixels are presented in Table 2.

The null hypothesis assumed for the photon upper-limit calculation is not optimal for the photon search. However the rough estimation of the possible photon signal could be made already in this setup. We optimize $\xi_{0}$ in each declination band with the same assumptions as in the previous section, and estimate the background in each pixel as the appropriately normalized number of protons that pass the cut: $N^{i}_{MC,p}(\xi>\xi_{0})$. For the photon excess calculation, the assumption of proton background is conservative as it should be higher than any mixed nuclei background, as it was discussed in Sec. 3.1. The background maps for various photon energies are shown in Fig. 4. The maxima among all pixels pre-trial photon candidate excesses over the proton background are presented in the Table 2 along with the average values of the proton background. The highest pre-trial excess significance, $3.43\sigma$ ($N_{\rm bg}=0.036$ and $N_{\rm obs}=2$), appears in the highest energy bin $E_{\gamma}>10^{20}$ eV, at $\{\alpha=155.3^{\circ},\delta=60.4^{\circ}\}$ pixel. To make a simple estimation of the post-trial p-value one can use the Bonferroni correction, i.e. to multiply the number of trials by the mimimum pre-trial excess p-value (Miller, 1981). In turn, the number of trials could be estimated as the number of non-overlapping pixel-size regions of the map which is several times smaller than the actual number of pixels. The resulting post-trial significances estimated in this way appear to be below $1\sigma$ level for all points of the sky at all considered energies. Therefore we conclude that, at the present level of point-source photon search sensitivity, there is no evidence for the photon signal. The actual results for each sky-map pixel at various energies are given in the text files supplemented to this paper. The format of the files is described in the Supplementary section.

The main systematic uncertainties for the photon-flux upper limits are related to the overestimation of the $E_{\gamma}$ parameter for hadron-induced events and to the uncertainty of the primary hadron mass-composition. The former uncertainty leads to the overestimation of the hadron background and subsequently to the looser photon-flux upper limit. As for the hadron mass composition uncertainty, the assumption of the proton composition which we used for the $\xi$-cut optimization could only make the photon-flux upper limit looser comparing to a mixed nuclei composition case. Therefore, the limits set are conservative with respect to the both of these uncertainties.

Finally, the last assumption that affects the result is the assumption of the background in (Eq. 3). The most conservative limits are set for zero background assumption, while the “real” mixed nuclei background assumption yields somehow more realistic limits.

4 Discussion and Conclusions

The upper limits are set on the fluxes of photons from each particular direction in the sky in the TA field of view, according to the experiment’s angular resolution with respect to photons. The only results of the ultra-high energy point-source photon flux upper limits were presented so far by the Pierre Auger experiment (Aab et al., 2014, 2017a). The comparison of those results to ours is not straightforward as the photon energy range of the Auger search, $10^{17.3}<E_{\gamma}<10^{18.5}$ eV, does not fully coincide to any of our ranges of search. Regardless of that, the average point-source photon flux upper-limit of Auger $\langle F_{\gamma}\rangle\leq 0.035\;{\rm km}^{-2}{\rm yr}^{-1}\;,$ is from two to three times lower than our average limit for the $E_{\gamma}>10^{18}$ eV. The results for the energies larger than $E>10^{18.5}$ eV are obtained here for the first time.

The point-source photon-flux upper limits derived in the present study can be used to constrain various models of astrophysics and particle physics. One can assume a distribution of photon sources and impose the constraints on their properties using the combination of point-source limits. In principle these constraints could be stronger than those derived from the diffuse photon-flux limits. The models that could be probed with the present photon point-source flux limits include cosmogenic photon generation models as well as top-down models of ultra-high energy photons production such as heavy decaying dark matter.

Supplementary

Photon point-source flux upper limits and photon excess pre-trial significances for all sky-map pixels are summarized in the separate file for each energy bin, named “limit_$[\log(E_{\gamma}/{\rm eV})]$.txt”. The file contains several columns with the following data.

Column 1: HEALpix pixel number (RING, started from 24)

Column 2: pixel $\alpha$, rad.

Column 3: pixel $\delta$, rad.

Column 4: $\xi$-cut value

Column 5: proton background value

Column 6: number of $\gamma$-candidate events

Column 7: 95% C.L. $F_{\gamma}$ upper limit, ${\rm km}^{-2}{\rm yr}^{-1}$ (zero background assumption)

Column 8: 95% C.L. $F_{\gamma}$ upper limit, ${\rm km}^{-2}{\rm yr}^{-1}$ (“real” background assumption)

Column 9: pre-trial $\gamma$ excess p-value

Column 10: pre-trial $\gamma$ excess significance

As it was mentioned in Sec. 3.3 the proton background value is used only for the calculation of photon excess p-value and significance, while the upper limits are calculated in either zero background assumption or “real” mixed nuclei background assumption. For pixels with the number of $\gamma$-candidates less than $p$ background both p-value and significance are set to zero.

Acknowledgments

The Telescope Array experiment is supported by the Japan Society for the Promotion of Science (JSPS) through Grants-in-Aid for Priority Area 431, for Specially Promoted Research JP21000002, for Scientific Research (S) JP19104006, for Specially Promoted Research JP15H05693, for Scientific Research (S) JP15H05741, for Science Research (A) JP18H03705 and for Young Scientists (A) JPH26707011; by the joint research program of the Institute for Cosmic Ray Research (ICRR), The University of Tokyo; by the U.S. National Science Foundation awards PHY-0601915, PHY-1404495, PHY-1404502, and PHY-1607727; by the National Research Foundation of Korea (2016R1A2B4014967, 2016R1A5A1013277, 2017K1A4A3015188, 2017R1A2A1A05071429) and Belgian Science Policy under IUAP VII/37 (ULB). The development and application of the multivariate analysis method is supported by the Russian Science Foundation grant No. 17-72-20291 (INR). The foundations of Dr. Ezekiel R. and Edna Wattis Dumke, Willard L. Eccles, and George S. and Dolores Doré Eccles all helped with generous donations. The State of Utah supported the project through its Economic Development Board, and the University of Utah through the Office of the Vice President for Research. The experimental site became available through the cooperation of the Utah School and Institutional Trust Lands Administration (SITLA), U.S. Bureau of Land Management (BLM), and the U.S. Air Force. We appreciate the assistance of the State of Utah and Fillmore offices of the BLM in crafting the Plan of Development for the site. Patrick Shea assisted the collaboration with valuable advice on a variety of topics. The people and the officials of Millard County, Utah have been a source of steadfast and warm support for our work which we greatly appreciate. We are indebted to the Millard County Road Department for their efforts to maintain and clear the roads which get us to our sites. We gratefully acknowledge the contribution from the technical staffs of our home institutions. An allocation of computer time from the Center for High Performance Computing at the University of Utah is gratefully acknowledged. The cluster of the Theoretical Division of INR RAS was used for the numerical part of the work. The lightning data used in this paper was obtained from Vaisala, Inc. We appreciate Vaisala’s academic research policy.