1 Introduction

In November 2012, the IceCube Collaboration announced the detection of two showerlike events with energies slightly above 1 PeV by analyzing the data taken during May 2010–May 2012 [1].

A follow-up analysis of the same data published in November 2013 revealed additional 26 events in the energy range \(\sim \)30 to 250 TeV [2]. Reconstruction of these events shows that 21 events are showerlike, mostly caused by \(\nu _e\) and \(\nu _{\tau }\) and seven are muon track events. These 28 events have flavors, directions, and energies inconsistent with those expected from the atmospheric muon and neutrino backgrounds and probably this is the first indication of extraterrestrial origin of high energy neutrinos. The track events have uncertainty of order one degree in their arrival directions and the angular resolution for 21 shower events is poor, ranging from\(\sim \)10\(^{\circ }\) to \(\sim \) \(50^{\circ }\). The IceCube analysis ruled out any spatial clustering of the events. The third year (2012–2013) data analysis revealed additionally nine events, of which two are track events and the rest are shower events [3]. The event 35 is the most energetic one so far observed. In the full 988-day data, the muon background is expected to be \(8.4\pm 4.2\) and for the atmospheric neutrinos we have \(6.6^{+5.9}_{-1.6}\). Five events are down going muons and are consistent with the expected background muon events. This shows that the IceCube events are predominantly shower events. For a \(E^{-2}_{\nu }\) spectrum the best fit diffuse flux obtained by IceCube per flavor is \(F_{\nu }=(0.95\pm 0.3)\times 10^{-8} \,\mathrm{GeV}\,\mathrm{cm}^{-2}\,\mathrm{s}^{-1}\,\mathrm{sr}^{-1}\), which is consistent with the Waxman–Bahcall bound [4]. The observation of these neutrinos triggered a lot of excitement as regards understanding of their origin and production mechanism. While interpreting these events in terms of astrophysical models seems challenging, several possible galactic and extra galactic sources have been discussed, which include the Galactic center [58], \(\gamma \)-ray bursts (GRBs) [9, 10], active galactic nuclei (AGN) [1114], high energy peaked blazars (HBLs) [1517], starburst galaxies [18, 19] etc. In Ref. [15] many positional correlations of BL Lac objects and galactic pulsar wind nebulae with the IceCube events are shown. It is also very natural to expect that these neutrinos might come from diverse sources having different production mechanisms and different power laws and this information can probably be extracted from the directionality of the observed neutrino events. The largest concentration of seven events is around the Galactic center and also clustering of the events could be associated to the Norma arm of the Galaxy [20]. As the statistics is too sparse, it is premature to draw any conclusion regarding the galactic origin of these events. There are also nonstandard physics interpretations of these events [2127].

We found coincidence of 12 TeV emitting HBL positions and the FR-I galaxy Centaurus A (Cen A) within the error circles of ten IceCube events from the online catalog for TeV astronomy (TeVCat) [28]. Due to the observed multi-TeV emission, these objects have long been believed to be sources of ultra high energy cosmic rays (UHECRs). A few years ago, the Pierre Auger (PA) collaboration reported two UHECR events within \(3.1^{\circ }\) around Cen A.

Therefore, in this work we focus our analysis on these candidate sources to find how the IceCube events with the desired energies can be produced through photohadronic interaction within the core region of the emanating jets.

2 Hadronic model

In the framework of the unification scheme of AGN, blazars, and radio galaxies [2931], all are intrinsically the same objects, viewed at different angles with respect to the jet axis. The blazars have jets pointing toward us. The double-peak spectral energy distribution (SED) structure is common to all these objects. This model is successful in explaining the multi-wavelength emission from BL Lac objects and FR-I galaxies [3236]. However, multi-TeV emissions during flaring and non-flaring events from these objects are difficult to reconcile. Also the most important challenge for the leptonic model is to explain the orphan flaring from the blazars 1ES1959\(+\)650 and Mrk 421. So variants of the hadronic models or the lepto-hadronic models are proposed to explain these multi-TeV emissions.

The AGNs are efficient accelerators of particles through shock or diffusive Fermi acceleration processes with a power-law spectrum given as \({\mathrm{d}N}/{\mathrm{d}E} \propto E^{-\kappa }\), with the power index \(\kappa \ge 2\) [37]. Protons can reach an ultra high energy (UHE) through the above acceleration mechanisms. Fractions of these particles escaping from the source can constitute the UHECRs arriving on Earth. These objects also produce high energy \(\gamma \)-rays and neutrinos through pp and/or \(p\gamma \) interactions [38, 39]. The multi-TeV flaring events in blazars can be well explained by invoking hadronic model through \(p\gamma \) interactions [4042]. Here it is assumed that the multi-TeV flaring in blazar occurs within a compact and confined region with a comoving radius \(R'_f\) inside the blob of radius \(R'_b\) [40] (henceforth \('\) implies the jet’s comoving frame). In the context of leptonic model, the SED of the HBLs (the synchrotron and the IC peaks) are fitted by taking into account different parameters (the blob radius \(R'_b\), magnetic field \(B'\), Doppler factor \(\delta \), bulk Lorentz factor \(\Gamma \) etc.). For the present work, instead of discussing in detail the SED of the individual HBLs, we use the best fit parameters from the leptonic models; these are shown in Table 1 and the references are given for these objects. As discussed earlier, in the inner region, the photon density \(n'_{\gamma ,f}\) is very high compared to the photon density \(n'_{\gamma }\) in the outer region.

Table 1 The objects HBLs/AGN are shown in first column which are in the error circles of the IceCube events ID (second column). After each object we also put their coordinates, Declination, and Right Ascension (Dec, RA) in degrees, and the redshift (z) and the Doppler factor (\(\delta \)). In the third and the fourth columns the observed neutrino energy \(E_{\nu }/\mathrm{TeV}\) and the corresponding seed photon energy \(\varepsilon _{\gamma }/\mathrm{keV}\) are given. In the fifth and the sixth columns the radius of the inner blob \(R'_{f}\) and the outer blob \(R'_{b}\) are given in units of \(R'=10^{15} R'_{15}\, \mathrm{cm}\). The seed photon density in the inner blob \(n'_{\gamma ,f}\) in units of \(n'_{\gamma ,f}=10^{10}\,n'_{\gamma ,f,10}\,\mathrm{cm}^{-3}\) is given in the seventh column and the diffuse neutrino flux \(F_{\nu }\) in units of \(F_{\nu }=10^{-9}\, F_{\nu ,-9}\, \mathrm{GeV}\,\mathrm{cm}^{-2}\,\mathrm{s}^{-1}\,\mathrm{sr}^{-1}\) is given in the eighth column. In the last column we show the \(\delta \chi ^2\) value for each event. The reference to each object is given in the first column

The UHE protons undergo photohadronic interaction with the seed photons in the inner region in the self-synchrotron Compton (SSC) regime through the intermediate \(\Delta \)-resonance. In a normal blazar jet, however, the photohadronic process is not an efficient mechanism to produce multi-TeV \(\gamma \)-rays and neutrinos because \(n'_{\gamma }\) is low, which makes the optical depth \(\tau _{p\gamma }\ll 1\). However, the assumption of the compact inner jet region overcomes this problem. The pion production in \(p\gamma \) collision through the \(\Delta \)-resonance is

$$\begin{aligned} p+\gamma \rightarrow \Delta ^+\rightarrow \left\{ \begin{array}{l l} p\,\pi ^0, &{} \quad {\mathrm{fraction}~ 2/3}\\ n\,\pi ^+, &{} \quad {\mathrm{fraction}~ 1/3}\\ \end{array} \right. . \end{aligned}$$
(1)

The \(\pi ^+\) and \(\pi ^0\) will decay to GeV–TeV neutrinos and \(\gamma \)-rays, respectively. The optical depth of the \(\Delta \)-resonance process in the inner compact region is \(\tau _{p\gamma }=n'_{\gamma ,f} \sigma _{\Delta } R'_f\), where \(n'_{\gamma ,f} \) is not known. By assuming that the Eddington luminosity is equally shared by the jet and the counter jet in the blazar, for a given comoving photon energy \(\varepsilon '_{\gamma }\) in the synchrotron/SSC regime we can get the upper limit on the photon density as \(n'_{\gamma ,f} \ll L_{Edd}/(8\pi R'^2_f \varepsilon '_{\gamma })\). We can also compare the proton energy loss time scale \(t'_{p\gamma }\simeq (0.5\,n'_{\gamma ,f}\sigma _{\Delta })^{-1}\) and the dynamical time scale \(t'_{d}=R'_f/c\) in this region to estimate \(n'_{\gamma ,f}\), so that the production of multi-TeV \(\gamma \)-rays and neutrinos take place. Not to have overproduction of neutrinos and \(\gamma \)-rays, we can assume a moderate efficiency (a few percents) by taking \(\tau _{p\gamma } < 1\), which gives \(n'_{\gamma ,f} < (\sigma _{\Delta } R'_f)^{-1}\). The kinematical condition for the production of \(\Delta \)-resonance in the observer’s frame is

$$\begin{aligned} E_p \varepsilon _{\gamma }=0.32 \frac{\Gamma {\delta }}{(1+z)^2}\, \mathrm{GeV}^2, \end{aligned}$$
(2)

where \(E_p\) and \(\varepsilon _{\gamma }\) are the proton and the seed photon energies, respectively.

In the decay of the \(\Delta \)-resonance to nucleons and pions, each pion carries \(\sim \) \(0.2\) of the proton energy and from the pion decay each neutrino and photon carries 1 / 4 and 1 / 2 of the pion energy, respectively. So the individual neutrino and photon energies are, respectively, \(E_{\nu }=E_p/20\) and \(E_{\gamma }=E_p/10\). This gives

$$\begin{aligned} E_{\nu } \varepsilon _{\gamma }=0.016 \frac{\Gamma {\delta }}{(1+z)^2}\, \mathrm{GeV}^2. \end{aligned}$$
(3)

In a HBL, \(\varepsilon _{\gamma }\) can be calculated from the given neutrino energy if \(\Gamma \) and \(\delta \) are known.

3 Results

We found coincidence of the positions of 12 HBLs and one radio galaxy, Cen A, within the error circles of ten IceCube events. These HBLs and AGN are taken from the online catalog TeVCat [28] and are observed in multi-TeV \(\gamma \)-rays. However, the redshift, Lorentz factor, and Doppler factor of some of these HBLs are not yet known. So whichever HBL has the known z, \(\Gamma \), \(\delta \), and SED and lies within the error circle of the IceCube event we calculate the seed photon energy \(\varepsilon _{\gamma }\) necessary to produce the desired neutrino energy \(E_{\nu }\) through the photohadronic interaction. The events 25 and 34 have very large errors, \(>\) \(40^{\circ }\), so we neglect these two events in our analysis. For the calculation of \(n'_{\gamma ,f}\), first we estimate the radius of the inner blob \(R'_f\), which will satisfy the restriction \(R_s < R'_f < R'_b\), where \(R_s=2 G_N M_{BH}/c^2\) is the Schwarzschild radius of the central black hole of mass \(M_{BH}\). The \(R'_b\) is obtained from the leptonic model fit to the SED of the object. The values of \(R'_f\) and \(R'_b\) for the objects are shown in Table 1. We assume a very conservative \(1~\%\) energy loss of the UHE protons in the inner blob on the dynamical time scale \(t'_{d}\), which corresponds to a optical depth of \(\tau _{p\gamma }\sim 0.01\) and \(n'_{\gamma ,f}\sim 2\times 10^{10}\, R'^{-1}_{f,15}\, \mathrm{cm}^{-3}\). The proton in the inner jet region has maximum energy \(E_{p,\mathrm{max}}\sim 3\times 10^{17}(B'_f/G) R'_{f,15}\, \mathrm{eV}\), where \(B'_f\) is the comoving magnetic field, which is higher than the outer region. For all neutrino flavors \(\alpha \) we assume a power-law spectrum of the form

$$\begin{aligned} J_{\nu _{\alpha }}(E_{\nu }) = A_{\nu _{\alpha } }\left( \frac{{E_{\nu }}}{{100\, \mathrm{TeV}}}\right) ^{-\kappa }, \end{aligned}$$
(4)

and the neutrino flux can be given as

$$\begin{aligned} F_{\nu }=\sum _{\alpha } \int _{E_{\nu 1}(1+z)}^{E_{\nu 2}(1+z)}\mathrm{d}E_{\nu } E_{\nu } J_{\nu _{\alpha }} (E_{\nu }). \end{aligned}$$
(5)

The normalization factor \(A_{\nu _{\alpha }}\) is calculated by using the 988 days IceCube data [2]. The integration limit is from 25 TeV to 2.2 PeV [57] and \(\kappa \) is the spectral index. For the luminosity distance calculation we take the Hubble constant \(H_0=69.6\,\mathrm{km}\,\mathrm{s}^{-1}\,\mathrm{Mpc}^{-1}\), \(\Omega _{\Lambda }=0.714\), and \(\Omega _m=0.286\). All the 37 IceCube events with their individual error circles in equatorial coordinates are shown in the sky map in Fig. 1. The 12 HBLs and Cen A are within the error circles of ten IceCube events, which are also shown in the sky map. In Table 1, we summarize all the relevant parameters of these 13 objects. All the correlated IceCube events are shower events with sub-PeV energies and the event 35 is the only PeV event with \(E_{\nu } \simeq 2\) PeV. Except the HBL, KUV00311-1993 [28], all other ones have their z, \(\Gamma \), and \(\delta \) measured/fitted and SEDs are calculated from the leptonic model. For most of the objects \(\varepsilon _{\gamma }\) lies between the synchrotron peak energy and the forward falling tail of synchrotron energy with the exception of RGBJ0192\(+\)017 [43] and 1ES1011\(+\)496 [46]. In these two HBLs \(\varepsilon _{\gamma }\) lies in the beginning of the SSC spectrum and the values are 179 and 69 keV, respectively. The corresponding photon densities and the neutrino fluxes are shown in Table 1. Our estimate of \(n'_{\gamma ,f}\) is based on the assumption of \(1~\%\) energy loss of the UHECR proton for all the HBLs/AGN. We observed that by varying \(\kappa \) between 2.2 and 3.08 we found a small variation is the neutrino flux. So here we fix its value to 2.2.

Fig. 1
figure 1

The sky map is shown in the equatorial coordinates with the 37 IceCube events and their individual errors (only for shower events). Here plus signs are for are shower events and the times sign are for track events with their corresponding event ID. We have also shown the positions of the HBLs with their names, which are within the error circle of the IceCube events. The TA hotspot is shown as a shaded closed contour and the galactic plane is shown as a dashed line

The HBL H2356-309 [44] is within the error circles of the three IceCube events 7, 10, and 21 and their corresponding synchrotron energies, \(n'_{\gamma , f}\), and neutrino flux are shown in Table 1. Two other HBLs, SHBLJ001355.9 [45] and KUV00311-1938, are also within the error circle of event 21 and SHBLJ001355.9 has the corresponding synchrotron energy \(\varepsilon _{\gamma }\simeq 45\) keV. The blazar PKS2005-489 [47] is in the error circles of the events 12 and 15 and to produce these neutrino events the photon energy is in the range 30–53 keV, which is near the synchrotron peak, and the corresponding proton energy is in the range \(1.2\, \le E_p \le 2.1\,{\mathrm{PeV}}\). These two events are also spatially correlated with the Fermi bubble. The event 17 has a mean energy of 200 TeV, correlated with the HBL, PG1553\(+\)113 [48] and is the farthest one in our list with a redshift of \(z=0.4\). The \(n'_{\gamma , f}\) and neutrino fluxes for PKS2005-489 and PG1553\(+\)113 are shown in Table 1.

Very recently the Telescope Array (TA) observed an UHECR hotspot above 57 EeV in a region within \(20^{\circ }\) radius circle centered at RA \(=\) \(146.7^{\circ }\) and Dec. \(=\) \(43.2^{\circ }\) [53], the shaded closed counter in the sky map in Fig. 1. This region correlates with the three neutrino events 9, 26, and 31. We found three HBLs: Mrk 180, 1ES0502\(+\)675, and RGBJ0710\(+\)591 within the error circle of IceCube event 31. Interestingly, the positions of two blazars, Mrk 421 [36] and 1ES1011\(+\)496 [46], are also simultaneously within the error circle of IceCube event 9 and within the TA hotspot [15, 54]. The required \(E_p\) and \(\varepsilon _{\gamma }\) for Mrk 421 are 1.3 PeV and 46 keV, respectively. The photon density and \(F_{\nu }\) are shown in Table 1. Similarly for 1ES1011\(+\)496 also we have shown the \(n'_{\gamma ,f}\) and \(F_{\nu }\) in Table 1.

Cen A is the nearest active radio galaxy and has since long ago been proposed as the source of UHECRs. A few years ago the Pierre Auger (PA) Collaboration reported two UHECR events above 57 EeV within \(3.1^{\circ }\) around Cen A [55]. Its position coincides within the error circle of IceCube event 35, having the highest neutrino energy so far observed by IceCube: 2 PeV. In terms of the hadronic model discussed above the 2 PeV neutrino energy corresponds to a proton energy of \(\sim \) \(40\) PeV and the seed photons energy is \(\varepsilon _{\gamma }\sim 56\) eV in the valley formed by the synchrotron and the SSC photons. The seed photon density \(n'_{\gamma ,f} \sim \times 10^{10}\, {\mathrm{cm}}^{-3}\) around \(\varepsilon _{\gamma }\sim 56\) eV is also high. For \(\varepsilon _{\gamma }\, <\, 56\) eV, synchrotron emission dominates and the low energy seed photon density increases rapidly [56]. So in principle \(E_{\nu } > 2\) PeV can be produced more efficiently. But non-observation of neutrinos above 2 PeV from Cen A can be due to: (i) a low flux of UHECR above 40 PeV and/or (ii) a cut-off energy around 40 PeV, beyond which the relativistic jet is unable to accelerate protons. Probably many more years of data taking are necessary to shed more light on this possible correlation between the IceCube event and the position of Cen A. The position of another HBL 1ES1312-423 also almost coincides with the position of the Cen A and thus falls within the error circle of IceCube event 35. For this HBL, \(\varepsilon _{\gamma }=0.32\) keV, and the corresponding observed photon flux is \(F_{\gamma }\sim 6\times 10^{-12}\,\mathrm{erg}\,\mathrm{cm}^{-2}\,\mathrm{s}^{-1}\), which is close to the synchrotron peak [52].

The multi-TeV flarings of the objects 1ES1959\(+\)650, Mrk 421, and M87 are interpreted through the photohadronic interaction as discussed in Sect. 2 [4042]. The maximum energy of these high energy \(\gamma \)-rays is less than 20 TeV (Mrk 421 [42]), which corresponds to a proton energy \(E_p < 200\) TeV and a neutrino energy \(E_{\nu }< 10\) TeV. But for the interpretation of the IceCube events the necessary proton energy will be \(E_p=20\times E_{\nu }\). For \(30\, {\mathrm{TeV}} \le E_{\nu }\le 2\, {\mathrm{PeV}}\), the proton energy will be in the range \(600 \,{\mathrm{TeV}} \le E_p \le 40\, {\mathrm{PeV}}\). So the neutrino flux from the interaction of these very high energy protons with the background photons can be small from an individual HBL. Apart from this, we have only observed flaring episodes of very few HBLs. So it is very hard to justify the temporal correlation of IceCube events during a flaring episode of a HBL. We have to wait a longer period and have sufficient data to comment on the correlation between the IceCube events and the flaring episode of the object.

In the photohadronic scenario the TeV–PeV neutrinos and the TeV–PeV \(\gamma \)-rays are correlated as both are produced from the decay of charged and neutral pions, respectively, as shown in Eq. (1). The background seed photons responsible for the production of these high energy neutrinos and \(\gamma \)-rays have energies above few keV. These photons have energies in between the synchrotron peak and the low energy tail of the SSC spectrum. The TeV–PeV photons produced from the \(\pi ^0\) decay will interact mostly with the same \(\sim \) keV seed photons in the inner blob region to produce \(e^+e^-\) pairs. The required threshold energy for the seed photon to produce the pair is \(\varepsilon _{\gamma ,\mathrm{th}}\ge 2 m_e^2/E_{\gamma }\), which is mostly in the microwave range. Also the \(\sigma _{\gamma \gamma }\sim 1.7\times 10^{-25}\, \mathrm{cm}^{-2}\) is the maximum in the microwave range and the pair creation cross section for the keV background photon is very small \(\sigma _{\gamma \gamma }\le 10^{-29}\,\mathrm{cm}^2\). In the region where the TeV–PeV photons and neutrinos are produced, the microwave photon density is very low. So even if the seed photon density is high (in the keV range), the mean free paths for the TeV–PeV photons satisfy \(\lambda _{\gamma \gamma }\gg R'_f\), hence, there will be negligible attenuation of these photons in the inner blob region. Again, in the outer blob, the low energy photon density is an order of magnitude smaller than the inner blob, so there is no attenuation in the outer region. However, on their way to the Earth, these TeV–PeV photons can interact with the low energy photons to produce pairs.

We have also done a statistical analysis to look for the correlation between the IceCube events and the 42 TeV emitting HBL/AGN from the TeVCat [28]. Here we adopt the method used in Ref. [57] and convert the coordinates (RA and Dec) into unit vectors on a sphere by

$$\begin{aligned} {\tilde{x}}=(\sin \theta \, \cos \phi , \sin \theta \, \sin \phi , \cos \theta ), \end{aligned}$$
(6)

with \(\phi =\mathrm{RA}\) and \(\theta =\pi /2-\mathrm{Dec}\), where i and j correspond to the event coordinates and object coordinates, respectively. The angle between the two unit vectors \({\tilde{x}}_i\) and \({\tilde{x}}_j\) is given as \(\gamma =\cos ^{-1} ({\tilde{x}}_j\cdot {\tilde{x}}_i)\), which is independent of the coordinate system and is a measure of the correlation between the events and the objects. Then one makes use of the quantity

$$\begin{aligned} \delta \chi ^2_i=\min \nolimits _j ( \gamma ^2_{ij}/\delta \gamma ^2_i ) \end{aligned}$$
(7)

where \(\delta \gamma ^2_i\) is the error on the ith coordinate. Only ten events meet the condition that \(\delta \chi ^2\le 1\) with the 13 objects which are shown in the sky map and also in Table 1. The \(\delta \chi ^2\) values of these events are given in the last column of Table 1. From the Monte Carlo simulation we estimate the significance of any correlation with IceCube events by randomizing the RA of the 42 objects within their allowed ranges. One has to remember that the value of \(\delta \chi ^2\) for the object closest to the neutrino event is chosen in this method. The distribution of \(\delta \chi ^2_i\) is realized by repeating this process one million times and the p value is calculated by counting the number of times 10 or more IceCube events satisfy \(\delta \chi ^2\le 1\) divided by the total number of realizations. In Fig. 2, the shaded histograms correspond to the number of correlated neutrino events with the 42 objects of the TeVCat in different ranges of \(\delta \chi ^2\) value. The open histograms correspond to the expected number of correlated neutrino events from the simulations (continuous line for the randomized RA) with their corresponding p-value, which is 0.647, which corresponds to a confidence level (CL) of \(\sim \) \(35~\%\). In another simulation we select the IceCube events which have angular errors \(\le \) \(20^{\circ }\). In this case IceCube events 7, 21, and 31 will not contribute. So with this constraint in the angular resolution, we have only seven events instead of ten events as considered earlier. In this simulation we found the CL \(\sim \) \(42~\%\); this is shown in Fig. 3. Both of these analyses show that there is no significant correlation between the IceCube events and the HBLs positions. As we have shown by increasing the angular resolution from \(40^{\circ }\) to \(20^{\circ }\) the CL increases by \(\sim \) \(7~\%\). Also, we believe that 42 objects from the TeVCat are not enough to give a better statistics when the events are isotropic. Apart from these objects there may be other types of sources which will contribute but are not included in our list. In the future we would like to consider more sources for our analysis.

Fig. 2
figure 2

The observed IceCube events (shaded histograms) and the simulated events (open histograms with continuous lines are for random RA) for different \(\delta \chi ^2\) distributions are shown for an angular resolution of the IceCube events \(\le \) \(40^{\circ }\). The p values for the open histograms are also given

Fig. 3
figure 3

Same as Fig. 2, but for an angular resolution of the IceCube events \(\le \) \(20^{\circ }\)

4 Conclusions

The astrophysical interpretation of the 37 TeV–PeV neutrino events by IceCube is challenging and several viable candidates have been proposed. HBL is one of them. The HBLs are the sources capable of producing multi-TeV \(\gamma \)-rays. In the photohadronic scenario, TeV \(\gamma \)-rays are accompanied by multi-TeV neutrinos from the decay of charged pions and kaons. By analyzing the online catalog TeVCat [28] we found coincidence of 12 HBLs and one FR-I galaxy Cen A position within the error circles of ten IceCube events. All these events are found to be shower events. The position of the HBL H2356-309 coincides with three IceCube events. We found that the positions of Mrk 421 and 1ES1011\(+\)496 are within the error circle of the IceCube event 9 as well as within the error circle of the TA hotspot. The observed highest energy PeV event coincides with the positions of Cen A and the HBL 1ES1312-423. Although from the statistical analysis we found no significant correlation between the IceCube events and the 42 objects in the TeV Catalog, it does not necessarily discard the photohadronic model interpretation for some of the IceCube events. Many more years of data taking are necessary to confirm or refute the positional correlations of the HBLs/AGN with the IceCube events. Also these possible candidate sources should be constantly monitored and studied in greater detail to have a better understanding of their properties and emission mechanisms.