© ESO 2019

1. Introduction

Blazars are a class of active galactic nuclei (AGN) characterised by rapid and large amplitude spectral variability, assumed to be due to the presence of a relativistic jet pointing very close to the line of sight (Blandford & Rees 1978; Antonucci 1993; Urry & Padovani 1995). The emission of these objects is non-thermal over most of the electromagnetic spectrum, from radio frequencies to hard γ-rays. The observed radiation shows extreme properties, mostly coming from relativistic amplification effects. The observed spectral energy distribution (SED) has a typical shape composed of two broad humps, one peaking between the far-infrared and the soft X-ray band, due to synchrotron emission, and the other peaking in the hard X-ray to γ-ray bands.

If the peak frequency of the synchrotron component (ν_peak) in ν − νf_ν space is higher than 10¹⁵ Hz, a blazar is usually called high-synchrotron peaked blazar (HBL or HSP; Padovani & Giommi 1995; Abdo et al. 2010); in particular, an HSP blazar with a synchrotron peak frequency above 10¹⁷ Hz is known as extreme high-energy synchrotron peaked (EHSP; Giommi et al. 1999; Costamante et al. 2001). Although dozens of objects clearly peaking at frequencies as high as 10^18.5 Hz (∼5–10 keV) have been found, evidence for the synchrotron peak reaching the MeV range is still under debate (Tanaka et al. 2014; Kaufmann et al. 2011; Tavecchio et al. 2011; Giommi et al. 2001).

Observations have shown that both EHSPs and HSPs are bright and extremely variable sources of high-energy γ-ray and TeV photons (TeVCat¹) and that they may be the dominant component of a putative extragalactic TeV background (Padovani et al. 1993; Giommi et al. 2006; Di Mauro et al. 2014; Giommi & Padovani 2015; Ajello et al. 2015). Given that most of the extragalactic objects detected so far above a few TeV are HSPs (Giommi et al. 2009; Padovani & Giommi 2015; Arsioli et al. 2015; The Fermi-LAT Collaboration 2017, see also TeVCat), HSP blazars are the main targets for future γ-ray and VHE observations (see examples in Chang et al. 2017; Arsioli & Chang 2017).

Arsioli et al. (2015, hereafter Paper I) built a catalogue of HSP blazars based on data from the Wide-field Infrared Survey Explorer mission (WISE); it was called 1 WISE HSP (1WHSP) and only selected sources inside a specific area of the colour–colour diagram called Sedentary WISE colour domain (SWCD), which was defined as an extension of the WISE blazar strip (Massaro et al. 2011, 2012; D’Abrusco et al. 2012) in order to include all sources from the Sedentary survey (Giommi et al. 1999, 2005; Piranomonte et al. 2007). They cross-matched the sources in the AllWISE catalogue (Cutri et al. 2013) that are inside the SWCD with different radio and X-ray catalogues using TOPCAT², applied three spectral slope criteria, and selected sources with ν_peak >  10¹⁵ Hz and Galactic latitude b >  |20° |. The slope criteria applied in Paper I are the radio to IR spectral slope, the IR to X-ray spectral slope, and the AllWISE channels W1 to W3 spectral slope³; the criteria are obtained from normalised and rescaled SEDs of three well-known HSP blazars.

Chang et al. (2017, hereafter Paper II) assembled the 2WHSP catalogue, an updated version of 1WHSP extended to lower Galactic latitudes (b >  |10° |) and including bright HSPs in the region close to the Galactic plane. Similarly to Paper I, the 2WHSP catalogue was constructed by cross-matching three radio catalogues (NVSS, FIRST, and SUMSS: Condon et al. 1998; White et al. 1997; Manch et al. 2003) with the AllWISE IR catalogue and then with various X-ray catalogues (RASS BSC and FSC, 1SWXRT and deep XRT GRB, 3XMM, XMM slew, Einstein IPC, IPC slew, WGACAT, Chandra, and BMW: Voges et al. 1999, 2000; D’Elia et al. 2013; Puccetti et al. 2011; Watson et al. 2009; Saxton et al. 2008; Harris et al. 1990; Munz 1992; White et al. 2000; Evans et al. 2010; Panzera et al. 2003). However, the 2WHSP catalogue does not apply the WISE colour–colour diagram, and any IR slope criteria (e.g. the W1 − W3 slope as measured from AllWise). This was done to avoid missing several HSPs where the IR and optical bands are dominated by the host galaxy thermal radiation, as many of the blazars that are classified as 5BZG in the latest BZcat catalogue Massaro et al. (2015).

Chang et al. (2017) used the SSDC SED tool⁴ to examine and fit the synchrotron component with a third-degree polynomial to get ν_peak and synchrotron peak flux (ν_peakf_νpeak) values for each candidate. The 2WHSP catalogue includes a total of 1691 sources with 540 previously known HSPs, 288 new HSPs, and 814 HSP-candidates. The name WHSP, which stands for WISE high-synchrotron peaked blazars, indicates that all sources in 2WHSP (except for one, 2WHSP J135340.2−663958.0) have WISE counterparts. For each 2WHSP source, we adopted as the best coordinates those taken from the WISE catalogue.

The 2WHSP catalogue has been applied as a seed to HE and VHE observations to find new VHE detections or counterparts of VHE catalogues. Arsioli & Chang (2017) analysed bright 2WHSP sources using archival Fermi-LAT Pass 8 data integrated over 7.2-year observations. By using the positions of 2WHSP sources as seeds for the data analysis, they found 150 new γ-ray detections not yet reported in any of previous γ-ray catalogues. The 150 new detections are collected in the First Brazil ICRANet γ-ray blazar catalogue (1BIGB).

Moreover, Padovani et al. (2016) cross-matched the 2WHSP and the HSP subsample of the second catalogue of Hard Fermi-LAT Sources (2FHL; Ackermann et al. 2016) with IceCube neutrinos (IceCube Collaboration 2015) suggesting that, among the blazar family, HSPs are the most probable counterparts for astrophysical neutrinos. Resconi et al. (2017) have presented further evidence of a connection between 2FHL HSPs, with very high-energy neutrinos and ultra high-energy cosmic rays (UHECRs) when cross-matching the 2FHL HBL subsample with UHECRs from the Pierre Auger Observatory (Abraham et al. 2004; The Pierre Auger Collaboration 2015) and the Telescope Array (Abu-zayyad et al. 2012). In a nutshell, HSP catalogues are important for HE, VHE, and multi-messenger astronomy. Their statistical properties, such as completeness, evolution, and possible bias associated with the building of HSP catalogues need to be analysed carefully.

Throughout the paper we adopt a Flat-ΛCDM cosmology with the following parameters (Carroll et al. 1992): Ω_m = 0.3 and H₀ = 70 km s⁻¹ Mpc⁻¹. Poisson errors of 1σ (Gehrels 1986) were applied when the numbers were ≤50 in the log N − log S.

2. Missing sources in the 2WHSP catalogue

There are a number of known HSP blazars that are not included in the 2WHSP catalogue, in particular all those that have not been detected in WISE surveys. For example, the blazar 5BZB J0403−2429 (Fig. 1, top) (Massaro et al. 2015) does not have an IR counterpart but is an extreme HSP with γ-ray counterparts in Fermi Third Catalog of Hard Fermi-LAT Sources (3FHL J0403.2−2428: The Fermi-LAT Collaboration 2017) and Fermi-LAT 8-year Source Catalog (4FGL⁵) catalogues. After checking all the HSP-candidates without IR detection from WISE, we realised that all of them were relatively close to another bright source and it was likely the reason why they were not included in the WISE source list. However, all of them have optical counterparts from the Panoramic Survey Telescope and Rapid Response System (PANSTARRs) or The United States Naval Observatory (USNO) catalogue, and may be robust HSP-candidates or have been catalogued as blazars, like 5BZB J0403−2429. This suggest that the presence of a nearby bright IR source affects the IR detections around it, and therefore offers a good explanation for the reason why several promising HSP-candidates have no WISE counterpart.

Fig. 1. Spectral energy distributions of 5BZB J0403−2429 (top) and 5BZG J0903+4055 (bottom). The sources were not selected for 2WHSP due to the absence of IR data (top) and the slope criterion (bottom).

Moreover, sources where a bright host galaxy dominates the IR and optical bands have low α_{1.4 GHz − 3.6 μm} values and therefore do not satisfy the radio-IR slope criterion adopted for the selection of 2WHSP blazars. This applies for instance to the case of 5BZG J0903+4055 (Fig. 1, lower panel), which is not in the 2WHSP catalogue since it does not satisfy the slope criterion (0.05 <  α_{1.4 GHz − 3.4 μm} <  0.45) owing to host galaxy contamination in the IR band, even thought it is an HSP with a highly variable X-ray spectrum and a luminosity that ranges between 2.4 and 7.3 × 10⁴⁴ erg s⁻¹ in the 0.3–10.0 keV band (Giommi et al. 2019), reflecting changes in ν_peak.

Finally, some sources without radio data could be suitable HSP-candidates with extreme ν_peak values. For example, the source 6dF J0213586−695137 (Fig. 2) is clearly an HSP with strong X-ray γ-ray emission, but with no radio counterpart in the Sydney University Molonglo Sky Survey (SUMSS) catalogue. The Fermi Pass 8 analysis (triangle points), 3FHL, and 4FGL data show that the γ-ray fluxes are consistent with data from the lower energy bands, suggesting that its inverse Compton emission might peak in the TeV band. This type of very faint radio source clearly could not be selected when using the criteria adopted for 2WHSP.

Fig. 2. Spectral energy distribution of 6dF J0213586−695137. This source was missed by 2WHSP because there is no radio data available for the source.

3. Building a more complete catalogue

As discussed above, the 2WHSP catalogue still misses some relevant HSPs. Among the factors that cause this deficiency, the slope criterion is the most relevant; also, as Fig. 1 shows, an HSP is not necessarily detectable by WISE. To recover these missing sources in this work we do not demand that all candidates have an IR counterpart; instead, we only apply the X-ray-to-radio flux ratio criterion. We note that for the new catalogue, all criteria are applied to find more blazars and not to eliminate sources that were already in previous WHSP catalogues. This was done only when we found new data implying that the sources are not HSPs.

3.1. Cross-matching and radio and X-ray slope criteria

The first step we took to increase the 2WHSP completeness was to cross-match the RASS and NVSS catalogues using a fixed 0.8 arcmin radius (which is larger than the positional uncertainty of more than 99% of the RASS sources) and to choose those matches for which the X-ray to radio flux ratio f_X/f_r >  9 × 10⁻¹¹ erg cm⁻² s⁻¹ Jy⁻¹. The X-ray RASS flux used in this work were corrected for Galactic absorption. We note that this value is obtained from the average radio-to-X-ray flux ratio of 2WHSP sources with ν_peak close to 10¹⁵ Hz. Sources that are already part of the 2WHSP catalogue and those that are close to the Galactic plane (|b|≤10°) were excluded. This procedure led to a list of 3011 additional pre-selected candidates.

We then cross-matched this list with the 5BZCat (Massaro et al. 2015), the XMM-Newton Optical Monitor Serendipitous UV Source Survey (XMMOM, Page et al. 2012), and the Fermi 3FHL (The Fermi-LAT Collaboration 2017) catalogues, using radii of 0.3, 0.3, and 20 arcmin, respectively (i.e. values that are somewhat larger than the largest positional uncertainties of these catalogues) to identify HSPs that have already been catalogued, and to check on all candidates having UV or γ-ray non-thermal emission. We used XMMOM and Fermi data to reduce the pre-selected sample since a source with UV or γ-ray detection is more likely to be a blazar. These cross-matches were not performed to exclude blazars that were already listed in other catalogues, on the contrary, they were done to help us find more HSP blazars, especially those already in 5BZCat or with a γ-ray (Fermi 3FHL) or UV (XMMOM) counterpart. This process led to 254 pre-selected sources, of which 30 are in 5BZCat, 149 in XMMOM, and 75 in Fermi 3FHL. Ten sources are in both 5BZCat and 3FHL, while three sources are in both XMMOM and 3FHL.

We note that all the cross-matchings described above use fixed radii that are relatively large compared to the average positional error of each catalogue. This was done to avoid loosing candidates in the first step of the selection process. Further refinement of all the tentative associations is carried out through the examination of the error circle map of each candidate. This map is produced by the VOU-Blazar tool (available at the Open Universe website⁶ and on GitHub⁷), which generates a spatial-plot where all the sources in the field and their position uncertainties (error regions) are visualised. After thoroughly checking this map for error circle consistencies, we inspected the SEDs of all the pre-selected candidates before making a final decision. Along with the SED, we also checked the optical spectrum, the bibliography on NED, and the radio, IR, and optical images for every source in order to classify each of them. Moreover, we excluded sources close to clusters of galaxies or members of a galaxy cluster, and removed all those with extended X-ray emission. For this purpose we cross-matched our list with a number of catalogues of clusters of galaxies (Zwicky, PLANCK, and MCXC: Zwicky et al. 1968; Planck Collaboration XXVII 2016; Piffaretti et al. 2011), and carefully checked the X-ray image for every 3HSP source, when available (see Sect. 2.4 in Paper II). Overall, we expect to have very few spurious sources due to contamination of X-rays from the hot gas in clusters. After the inspection of the 254 pre-selected sources, 58 of them were determined to be classifiable as HSPs or HSP-candidates and were added to the current HSP catalogue.

There are still 2757 radio–X-ray matching sources to be checked. In this version of the catalogue we added the 58 sources that have a γ-ray detection or that have already been listed in 5BZCat. Given that nearly half of the sources in our HSP catalogue have been detected in the γ-ray band, we expect that there should be approximately 60 HSPs with γ-ray flux below the Fermi-LAT sensitivity among the other 2757 sources. In a future update of the catalogue we will use the VOU-Blazars tool (see Sect. 3.3) to carefully examine the remaining 2757 radio–X-ray matching sources.

3.2. Searching for extra sources using Fermiγ-ray catalogues

We note that robust HSP-candidates, especially those with γ-ray detection, do not necessarily require both radio and X-ray data to be present in the currently available archives, as shown by Fig. 2. There are still several sources that cannot be selected on the basis of the radio/X-ray flux ratio, but that are detected in γ-rays. To identify these HSPs, a careful examination of all the sources in the 3FHL catalogue and their possible HSP counterparts was performed. The recently released Fermi 4FGL catalogue⁸ (The Fermi-LAT Collaboration 2019) was also searched in an effort to find new HSP sources. HSPs typically have γ-ray photon index Γ <  2.0 (see e.g. the Fermi-3LAC catalogue; Ackermann et al. 2015), thus we checked those 4FGL γ-ray sources with hard γ-ray slopes.

To search for additional HSPs among Fermi sources, we assumed that a Fermi γ-ray detected blazar has a counterpart in existing radio or X-ray catalogues, but not necessarily in both. Thus, we checked every radio and X-ray source around the Fermi detections with the inspection tools available via the Open Universe portal. Then we searched for possible optical, IR, and UV counterparts for them. Even though there might be more than one optical source within the Fermi error region, we inspected those that have a radio or X-ray counterpart. After that we carefully examined every possible counterpart with a multi-frequency approach (see Sect. 3.1). This lead to a total of 168 HSPs that were not in the 2WHSP sample, but were in 3FHL. We note that here the number is approximately three times larger than in the last step (168 vs. 58 sources), suggesting that most of the sources with radio and X-ray matches had already been selected in previous steps. Moreover, out of the 389 4FGL sources with a hard spectrum that are still out of our selection, 121 HSP blazars or candidate blazars were identified and added to our catalogue after carefully checking their SEDs.

As expected, approximately half of the sources added from the 3FHL or the 4FGL catalogue do not have an associated radio or X-ray detection. These sources are either without radio or X-ray data or have an X-ray-to-radio flux ratio that is not typical of previously known HSPs (i.e. f_X/f_r <  9 × 10⁻¹¹). We discuss the selection of sources without a radio counterpart in the Appendix A. In conclusion, we only select HSP-candidates after a careful review of each field, using all available information from multi-frequency catalogues.

3.3. The VOU-Blazar tool

To increase the efficiency of our search for new HSPs, we developed a tool called VOU-Blazars⁹ that uses Virtual Observatory (VO) protocols to retrieve multi-frequency information from a large number of services distributed worldwide and combines the retrieved data to find sources with SEDs similar to those of blazars. This tool allowed us to select 48 additional HSP-candidates. Most of these new sources have been detected only in recent Swift X-ray telescope (XRT) or X-ray Multi-Mirror Mission (XMM-Newton) observations, and therefore are not included in the catalogues used by the SSDC-SED builder tool provided by the Italian Space Agency’s Space Science Data Center (ASI-SSDC). With the VOU-Blazars, it is also possible to retrieve results from recent X-ray observations, and this allowed us to significantly improve our search with respect to the 2WHSP catalogue. The VOU-Blazars tool has already been used to locate all the probable and confirmed blazars within the uncertainly region of the astrophysical neutrino IceCube-170922A (Padovani et al. 2018).

3.4. SED selection criteria and the estimation of synchrotron peak

The process discussed above led to the addition of a total of 395 new sources to the updated HSP catalogue. Only objects having non-thermal data in at least three different bands were considered.

We determine the ν_peak for sources without a radio or X-ray counterpart as follows: for sources without an X-ray measurement we verified that ν_peak was higher than 10¹⁵ Hz using radio, IR, and UV data. For sources without a counterpart in current radio catalogues we used the (non-thermal) IR, UV, and X-ray data (see the Appendix for more details). If we were able to identify a non-thermal component in the IR, UV, or X-ray bands, we estimated the ν_peak value by fitting those data points. This was done only when the source had at least one data point in the X-ray band, otherwise a lower limit or an uncertain value was derived. Only a very few candidates without radio data were found; in these rare cases we used the UV and X-ray data to estimate the synchrotron peak position. In cases where no UV data were available, the synchrotron peak was estimated only if a good quality X-ray spectrum (producing several SED good signal-to-noise points) could be found in the Swift XRT or XMM-Newton public archives.

The criteria used for source selection are somewhat inhomogeneous given that different data are available for each source. In practice, we selected a source based on the following criteria:

• the SED includes at least three data points that can be attributed to non-thermal emission from, for example, radio, IR/optical, X-ray; radio, optical/UV, X-ray; IR, optical/UV, X-ray; or radio, IR, UV;

• the optical spectrum and the SED is blazar-like (see e.g. Fig. 1 of Padovani et al. 2017) and the non-thermal data are sufficient to reliably establish that ν_peak ≥ 10¹⁵;

• the source is included in a blazar catalogue such as 5BZCat or Sedentary;

• the radio-to-synchrotron peak flux ratio is similar to that of standard HSPs, i.e. ≈4 − 5 orders of magnitude;

• the radio emission is within the X-ray position error, and the optical, IR, UV counterparts coincide with the radio emission;

• the source was detected in a very high-energy survey and the SED built with the γ-ray, and the lower energy data is consistent with that of an HSP blazar.

3.5. Estimation of spurious sources

The contamination of the 3HSP sample with other types of objects is very much dependent on the position uncertainties of each individual multi-frequency detection. For the γ-ray detected sources, the number of expected spurious associations is less than three, while for other sources it depends on the quality of the available X-ray data, which is heterogeneous. The uncertainties are mostly dominated by X-ray sources that have only been detected in the RASS survey, given that for those cases the positional error can be relatively large, up to ∼40 arcsec. For sources with a more precise positioning (1–5 arcsec), like those detected by XMM-Newton or Swift, the expected number of spurious sources is very small, likely less than one.

The number of spurious associations between 3HSP objects and γ-ray sources can be estimated from the Fermi positional uncertainty, the number of Fermi sources, and the number density of HSPs. For the 4FGL sources above the Galactic plane (b ≥ |10|°), the average error ellipse major axis is 4.46 arcmin, and the average error ellipse minor axis is 3.59 arcmin. Thus, the average area associated with the positional uncertainty of 4FGL sources at b ≥ |10|° is ∼0.014 deg². There are 3663 4FGL sources with b ≥ |10|°, so the total area covered by 4FGL error regions is 0.014 × 3663 = 51.3 deg². The radio log N − log S of 3HSP gives a surface density of 3HSP objects with radio flux > 3.5 mJy of ≈0.05 deg⁻². Multiplying the total area 51.3 deg² by the number density 0.05, we get 2.6. This should be seen as an upper limit since this calculation assumes that all 3HSP sources have a radio flux density of 3.5 mJy. Obviously, the radio flux density of the large majority of our sources is significantly higher than the limit of current radio surveys. In conclusion, the expected number of spurious association of 3HSP objects with (high b) Fermi 4FGL sources is less than 2.6 sources, that is < 2.6/2011 or < 0.15%.

The expected number of spurious associations related to the radio–X-ray matching sources can be estimated as follows. The average uncertainty radius in the RASS survey is 18.8 arcsec, corresponding to an area of 8.6 × 10⁻⁵ deg². This value, multiplied by the number of sources in the RASS catalogue at b ≥ |10|° and Dec >  − 40°, gives an area of 7.6 deg², corresponding to the surface covered by all RASS error regions in the part of the sky covered by the NVSS survey. Given that in the NVSS catalogue at b ≥ |10|° there are ∼43.5 sources per square degree, the expected number of random NVSS/RASS matches is ∼330. Since only 64% of sources in the NVSS catalogue matches a WISE source, and only 14.5% of the WISE sources are located in the SWCD area, the expected fraction of NVSS/RASS matches including a WISE source in the SWCD area is 0.09. The total number of expected spurious associations that would meet our selection criteria is therefore 30.6 = 330 × 0.09, equivalent to less than 2% of the 3HSP sources at Dec >  − 40°.

The equivalent calculation applied to the subsample of 3HSP objects with no infrared counterpart but with UV information, with WISE replaced by Galaxy Evolution Explorer (GALEX) and the SWCD replaced by UV to X-ray spectral slope constraint, gives an expected value of about nine spurious sources. These numbers should be considered conservative upper limits as in most cases, in addition to the basic multi-frequency data considered above, we have information such as radio measurements at different frequencies, X-ray data with much better positioning than that of the RASS survey and γ-ray data, for example, that is fully consistent with the assumption of an HSP blazar, both in terms of positional uncertainties and SED constraints. All sources with poor multi-frequency coverage and large X-ray positional uncertainty have been flagged as “candidate”.

A careful inspection, based on newly available multi-frequency data, has been carried out on the sources of the 2WHSP catalogue, resulting in the elimination of 73 objects (4.32% of the total) as spurious associations or blazars with intermediate ν_peak energy. As explained above we expect the fraction of incorrect 3HSP associations to be significantly lower than that of the 2WHSP catalogue.

3.6. Comparison between new sources and the 2WHSP catalogue

A comparison between the new sample and the 2WHSP catalogue is shown in Figs. 3 and 4, normalising on the total number of objects in each sample. Figure 3 shows that many of the new sources have ν_peak close to the threshold of 10¹⁵ Hz and have higher radio flux compared to the 2WHSP sample. Figure 4 indicates that there are more high α_IR − X sources in the new sample, suggesting that relaxing the slope criteria allows the selection of new HSPs with relatively faint X-ray emission. The figure also shows that the α_r − IR distribution for the subsample of new sources is similar to that of the 2WHSP sample. This is a strong indication that 2WHSP is close to complete regarding galaxy dominated HSP sources. We have selected almost all the nearby and bright HSPs. Around ten sources have α_r − IR higher than 0.5, and we examined all of them thoroughly. According to their radio morphology and data (1.4 GHz and 5 GHz Planck data), most of them are moderately radio extended or have misaligned jets. Their IR emission is on average lower than that of other sources and therefore have relative low IR-to-radio flux ratios.

Fig. 3. Distribution of the new sources and those in the 2WHSP catalogue (νpeak, top; radio flux density, bottom) normalised according to the number of objects in each sample.

Fig. 4. Distribution of the new sources and the 2WHSP (X-ray–IR slope, top; radio–IR slope, bottom) normalised based on the total number of objects in each sample.

In conclusion, compared to 2WHSP, the new catalogue includes more sources with ν_peak close to the 10¹⁵ Hz limit with brighter radio fluxes. By examining all possible blazar counterparts in γ-ray catalogues, we were able to find additional HSP sources with γ-ray emission, which we might have missed if we had only relied on our radio and X-ray selection methods (Sect. 3). We called the new catalogue 3HSP, removing the “W” as the WISE counterpart is not a requirement in the new catalogue.

3.7. Redshift estimations

It is well known that many BL Lacs have no redshift determination because of the lack of any detectable feature in their optical spectra. However, for a good fraction of our HSPs, the signature of thermal emission from the host galaxy can be recognised in their broad-band SEDs, especially in the IR band, and this can be used to obtain a photometric redshift (photo-z) estimation.

Urry et al. (2000) showed that the host galaxies of BL Lac objects are giant ellipticals with almost constant absolute magnitude of M_R = −23.5 (see their Fig. 5 for details). By fitting the host galaxy contribution to the SED using a giant elliptical galaxy template (Coleman et al. 1980) and assuming that the elliptical galaxy template is approximately a standard candle proxy, with luminosity equal to M_R = −23.5 a photometric redshift can be estimated. Therefore, whenever the host galaxy contribution could be distinguished from the non-thermal emission in the SED, we applied this method to estimate the redshifts of all sources with a featureless optical spectrum and those for which no optical spectrum is available.

An example of the application of this method is shown in Fig. 5 (top left) where we show the optical spectrum of the source 3HSP J154433.1+322148, taken from the Sloan Digital Sky Survey (SDSS) Dr14 (Blanton et al. 2017; Abolfathi et al. 2018). The spectrum is clearly featureless and no spectral redshift can be estimated. In the bottom pane we show the radio to X-ray SED of the same source. The top panel of Fig. 5 shows that the slope of the SDSS spectrum changes above ≈7600 Å. According to the SED and the optical spectrum, the SDSS i and z bands are clearly from the host galaxy. We then used only the r, g, and u bands in the optical to UV to fit ν_peak. In addition, the data extracted directly from the SDSS DR14 optical spectrum fits well with the optical SED and indicates the same thing. By fitting the Two Micron All-Sky Survey (2MASS), WISE W1 and W2, and part of the SDSS data using the giant elliptical template available in the SSDC SED builder, we obtained a photo-z estimation of 0.32.

Fig. 5. Example of photometric redshift (photo-z) estimation. The SDSS Dr14 optical spectrum (top) and SED (bottom) of the source 3HSP J154433.1+322148. The red line represents the non-thermal component, the green line the giant elliptical template fitted to the IR (and partly optical) data at z = 0.32, and the blue line the data extract from the SDSS DR14 optical spectrum.

For sources with a featureless spectrum and where non-thermal emission was dominating at all wavelengths, we estimated a redshift lower limit. This was done, as in Paper I (1WHSP), by assuming that in the optical band the host galaxy is swamped by the non-thermal emission and leaves no imprint on the optical spectrum when the observed non-thermal flux is at least ten times higher than the host galaxy flux (for details, see Eq. (5) of Paper I).

To estimate the uncertainty of the photometric evaluation of z, we randomly selected 50 sources with firm z measurements, and plotted in Fig. 6 the photometric estimation vs. the optical spectra estimation (from SDSS Dr14) of z. The figure shows that the differences between the two redshift estimations is at most 0.1 for sources with z <  0.6, while for z >  0.6, the error is at most 0.3. However, the estimation differences for the majority of these high-z sources is around 0.2.

Fig. 6. Comparison between photo-z and spectral-z. The central solid line has slope of one and represents the equality between the two redshift estimations, while the red, green, and blue dashed lines indicate a scatter of 0.1, 0.2, 0.3, respectively.

There are 80 sources with photo-z higher than 0.6, and we double checked all of them. The error might be large given that the contrast between the host galaxy features and the non-thermal emission is low and not easy to resolve. The amplitude of the error also depends on the availability of IR–optical data from host galaxy. For some sources it is more difficult to tell the origin of the IR–optical emission; therefore, they are marked with appropriate source and redshift flags.

Figure 6 illustrates that the photo-z estimation for high-z are underestimated when comparing with spectral-z. One of the most likely reasons for this bias is due to the difficulty in finding the correct position of the 4000 angstrom break in the optical spectrum. For a high-redshift/high-luminosity source the non-thermal flux may be higher than that of the host galaxy at the position of the break. It could be that the “real” break is at a lower frequency than where the non-thermal flux merges with the emission from the host galaxy, which is where a slope change can be detected in the optical spectrum. That is why we underestimated the photo-z from the SED photometric data.

Recently, a new spectroscopic redshift estimation has become available for 3HSP J062753.2−151956 with spectral-z = 0.31(Paiano et al. 2018), which is very close to our photo-z = 0.29.

The breakdown of redshift determinations is as follows:

• 31.8% of 3HSPs with firm redshift (flag 1);

• 5.3% of 3HSPs with uncertain redshift (flag 2);

• 4.7% of 3HSPs with lower limit redshift (flag 3);

• 7.2% of 3HSPs with photometric redshift and featureless optical spectrum (flag 4);

• 39.0% of 3HSPs with photometric redshift and without optical spectrum (flag 5);

• 11.9% of 3HSPs without any redshift estimation or measurements.

In Fig. 7, we show the redshift distribution for the entire 3HSP catalogue. The figure suggests that the redshift of 3HSP sources is centred around z = 0.35, which is similar to previous results (e.g. 1/2WHSP, 5BZCat, Sedentary Survey).

Fig. 7. Redshift distribution of the 3HSP sample. The red line represent all the 3HSP sources with a redshift estimation, while the blue line represents sources with an estimation, but not lower limits.

The redshift of the sources without redshift might be significantly higher than the average value. According to the Monte Carlo simulations of Giommi et al. (2012a), the average predicted redshift of BL Lacs with featureless spectra is around 1.2. Figure 9 of the 2WHSP paper also suggests that HSPs with featureless spectra might be much more distant than we think.

3.8. XRT data analysis with the Swift Deepsky pipeline

During the selection of 3HSP objects, we proposed Swift observations for 210 sources without a good X-ray spectrum (50 from 1WHSP, 80 from 2WHSP, and 80 from 3HSP), 190 of which were kept in the 3HSP final version. As of March 2018, 151 of these sources have been observed by Swift. We analysed all the XRT data using the Swift DeepSky software¹⁰ (Giommi et al. 2019), a pipeline tool based on HEASoft6.25 and the XIMAGE package that we assembled and that can run on a Docker container¹¹ (Morris et al. 2017). This software detects sources and estimates fluxes (or upper limits) in four energy bands. For the sources detected with at least 100 photons we used XSPEC to estimate a best fit spectrum using power law and log parabola models (see Fig. 8).

Fig. 8. Aitoff projection of 3HSP sources in Galactic coordinates. Red, blue, and green points represents FOM >  1.0, 0.2 <  FOM <  1, and FOM <  0.2, respectively (see text for details).

Our analysis resulted in 147 new X-ray detections in association with 3HSP sources. All new X-ray data, together with optical and UV measurements from Swift Ultraviolet/Optical Telescope (UVOT), were used to estimate or update ν_peak (analysed using the SSDC on-line interactive analysis tool). The high Swift detection rate clearly shows that the use of multi-frequency information for selecting X-ray targets is very effective for the detection of blazars in that band; the only four non-detections (< 1%) are objects observed with very short exposures.

4. 3HSP: the largest and most complete catalogue of HSP blazars to date

In total, the new version of the catalogue includes 2013 sources, 88 of which are close to the Galactic plane (|b|< 10°). This catalogue has been named 3HSP, which means the third catalogue for HSP blazars. The W in the acronym has been removed in this edition of the catalogue as sources are no longer required to be detected in the WISE all sky survey. Of the 2013 3HSP sources, 1618 are also in 2WHSP, 657 are in 5BZCat, and 1007 have a counterpart in one of the Fermi-LAT catalogues. Figure 8 illustrates the Aitoff projection map of the 3HSP catalogue in Galactic coordinates, and most of the 3HSP sources are located out of the Galactic plane.

The content of the 3HSP catalogue presented in this paper differs from 2WHSP as detailed below:

• it includes 395 new HSP blazars/candidates;

• 73 sources that were in the 2WHSP catalogue have been removed, due to reclassification as intermediate synchrotron peaked (ISP), low-energy synchrotron peaked (LSP) objects¹², or spurious associations based on new optical spectra and XRT data;

• we added sources mainly based on the radio/X-ray flux ratio and γ-ray catalogues;

• photo-z values were estimated for 930 sources;

• appropriate flags were assigned to uncertain sources.

We summarise the steps followed to build the 3HSP catalogue here:

1. Cutting the NVSS-RASS matched sources with flux ratio ≥9 × 10⁻¹¹ erg s⁻¹ cm⁻² Jy and excluding those already in 1/2WHSP selection process, leading to 3011 pre-selected sources;

2. Cross-matching the pre-selected sources with 5BZCat, XMMOM, and Fermi 3FHL, reducing the pre-selected sample to 254 sources;

3. Examining the pre-selected candidates one by one, adding new 58 HSPs and HSP-candidates to the catalogue;

4. Checking additional blazars from the Fermi 3FHL catalogue (168 sources) and Fermi 4FGL catalogue (121 sources);

5. Cleaning non-confirmed sources, especially those with no Fermi counterpart, removing 73 sources previously listed in 2WHSP;

6. Adding 48 new sources identified using the VOU-Blazars tool.

Table 1 lists the ν_peak, redshift, γ-ray counterpart, 2WHSP counterpart, BZCat counterpart, and a figure of Merit (FOM) of a representative subsample. The full catalogue is available in electronic form¹³. The FOM was defined in Arsioli et al. (2015) as the ratio of the flux at the synchrotron peak (ν_peakf_νpeak) of a given source to the peak flux of the faintest 1WHSP blazar that had been detected in the TeV Band. Here we re-define the FOM in units of 2.5 × 10⁻¹² erg s⁻¹ cm⁻² (Log ν_peakf_νpeak = −11.6), which is the peak flux of the faintest 3HSP source in the current version of TeVCat. This parameter was introduced to provide a simple quantitative measure of potential detectability of HSPs by TeV instruments (see Paper I for more details).

Some of the 3HSPs sources have sparse non-thermal data, poor quality data, or peculiar flux ratios. These sources will need to be carefully followed up in future versions of the HSP catalogue. These objects were assigned a flag value that reflects the reason for the uncertainty. For sources with little non-thermal data or with fairly large X-ray positional errors, we gave a flag value of 1. Flag 2 is for cases with one or two doubtful data points (due to large positional uncertainty or problematic photometry) but that still provide a reasonable HSP-like SED. Sources with a low ratio between synchrotron peak flux and radio flux (possibly due to jet misalignment) were given source flag 3. Source flag 4 is for cases where the observed IR or optical emission cannot be safely attributed to the host galaxy or to synchrotron emission.

One of the main purposes of the 3HSP catalogue is to provide possible seed sources for future HE and VHE observations, so we would like to select as many candidates as possible. There might be a small fraction of sources that still need further data to be confirmed as HSPs. Among those with an uncertain flag, sources with flag 1 are the most ambiguous, but we have no reason to remove them from the current version of the catalogue (some of them even show a blazar-like optical spectra), although we highlight the need of a follow-up to confirm their classifications.

There are 374 sources (18.6%) assigned with flag 1 in the 3HSP catalogue, which implies an efficiency of at least 81.4%. We note that less than 5% of the sources in 2WHSP turned out to be spurious HSPs. Therefore, given that we now have access to more data and better tools, we expect a lower rate of spurious classifications for the 3HSP catalogue. A more detailed estimation of the number of spurious sources is given in Sect. 3.5. In addition, we have a good record of selecting HSP-candidates for further X-ray observations, especially with Swift. We have already had three successful Swift observation campaigns carried out based on 1/2WHSP sources, resulting in 200 dedicated observations to HSPs and HSP-candidates.

Flags are also associated with ν_peak and redshift values with the following meaning:

1. firm estimation;

2. uncertain value;

3. lower limit;

4. photometric redshift of an object with featureless optical spectrum;

5. photometric redshift of an objects for which no optical spectrum is available.

We note that flags 4 and 5 only apply to redshift. Moreover, the uncertain value for the ν_peak estimate means that we are not sure exactly what the synchrotron peak frequency is, due to insufficient non-thermal data, but we could still tell that the frequency is higher than 10¹⁵ Hz. Sources with an uncertain flag for synchrotron peak frequency are not necessarily candidates, as some of them have blazar like optical spectrum or have already been included in 5BZCat. The source flag and the synchrotron peak frequency flag are marked independently.

As shown in Fig. 5, sources with a redshift flag equal to 4 have featureless optical spectra, but the emission from the host galaxy is not completely overwhelmed by the non-thermal radiation, thus a photo-z can still be estimated from IR data or in part of the optical band. Flag 3 (lower limits) sources also have featureless spectra; however, their SEDs are totally dominated by synchrotron radiation. There are still some sources with SEDs that are non-thermal dominated and for which no optical spectrum is available. The redshift in this case remains blank; this applies only to 11.94% of the sources.

5. Completeness and statistical properties of the 3HSP catalogue

The demographic properties and the cosmological evolution of blazars have been extensively debated by the community, and for a long time the existence of the so-called blazar sequence has been a controversial topic. With the largest ever HSP blazar catalogue, the overall properties of the 3HSP sample can be discussed thoroughly, and here we present arguments in tension with the blazar sequence scenario. We start by checking the completeness of the 3HSP catalogue in the radio and the X-ray energy bands.

Figure 9 shows the radio–X-ray flux plane of the 3HSP, comparing it with the Sedentary survey (Giommi et al. 1999, 2005; Piranomonte et al. 2007) and the Deep X-ray Radio Blazar Survey (DXRBS, Padovani et al. 2007). The corresponding approximate synchrotron peak frequencies were estimated using extrapolation, and the relationship between ν_peak and the X-ray-to-radio flux ratio (f_X/f_r) with the equation, ν_peak = (f_X/f_r + 16.068 ± 0.306)/0.377 ± 0.019. According to the error of the two fit parameters, we estimate that the uncertainty on this estimation is around one order of magnitude for ν_peak. We note that this relationship was derived from the 3HSP subsample, so we only convert the radio-to-X-ray flux ratio to synchrotron peak frequency for HSPs (or equivalently ν_peak >  10¹⁵ Hz).

Fig. 9. Sampling of the radio flux density–X-ray flux plane with the 3HSP, Sedentary, and DXRBS samples. The blue dashed and dotted lines are the slope limits for the 3HSP sample. The black dotted line is the flux ratio that separates LBL and HBL. The green and yellow dotted lines denote the forbidden region for BL Lacs and FSRQs. The red dashed lines indicate the incomplete regions of 2WHSP.

The black dotted line indicates the X-ray-to-radio flux ratio that separates low-energy peaked BL Lacs (LBL) and HBL, which is f_X/f_r = 3.2 × 10⁻¹² erg cm⁻² s⁻¹ Jy⁻¹ in Padovani et al. (2003) and Padovani (2007). Three black dashed lines represent the flux ratios, from left to right, f_X/f_r = 8 × 10⁻¹⁴,  3 × 10⁻¹⁰, and 2 × 10⁻⁹ erg cm⁻² s⁻¹ Jy⁻¹, respectively. The first value represents the minimum X-ray-to-radio flux ratio for BL Lac objects and for flat spectrum radio quasars (FSRQs; Padovani 2007) based on the data from Padovani (1997) and Siebert et al. (1998). The second value is the maximum flux ratio for FSRQs, while the third is the maximum flux ratio for BL Lacs (ν_peak ≈ 10²⁰ Hz: Padovani et al. 2003; Padovani 2007).

Among the radio and X-ray catalogues used for the 3HSP selection, NVSS and RASS have the largest sky coverage, and we used the NVSS and RASS catalogues to estimate the radio and X-ray limits for Fig. 9. Given that the minimum radio flux cut applied for the 3HSP–NVSS subsample, the radio limit for this 3HSP subsample in this figure is set to 3.5 mJy, while the X-ray limit for the 3HSP–RASS subsample is set to the minimum RASS flux value in the subsample, ≈2 × 10⁻¹³ erg cm⁻² s⁻¹. However, it should be noted that in some cases radio and X-ray flux limits can be lower as we used several other radio and X-ray catalogues to build the 3HSP sample. In practice, the exact radio and flux limits of 3HSP are lower than the limits set here.

According to Fig. 9 the 3HSP is not complete, neither in radio nor in X-rays. Sources with radio flux brighter than 3.5 mJy but X-ray flux fainter than 2 × 10⁻¹³ erg cm⁻² s⁻¹ (upper red dashed region) are not selected as they are not detected by current large-area X-ray surveys. On the other hand, some bright X-ray but faint radio sources (lower red dashed region) are missed, since these sources are not included in today’s large-area radio catalogues. If we were to increase the radio flux limit to ≈22 mJy, we would define a complete, radio flux-limited HBL sample. Similarly, when setting the X-ray flux limit to a higher value ≈7 × 10⁻¹² erg cm⁻² s⁻¹, the sample would become complete in the X-ray band.

The blue dotted line in the figure represents the radio to X-ray ratio f_X/f_r = 9 × 10⁻¹¹ erg cm⁻² s⁻¹ Jy, which is the flux criteria applied when adding new sources that are not in the 2WHSP catalogue. This value is just slightly lower than the average value we obtained from the old 2WHSP subsample with ν_peak close to 10¹⁵ Hz. The left blue dashed line (f_X/f_r = 1 × 10⁻¹¹) corresponds to the slope criteria applied when building the 1/2WHSP, and yields a relatively low ν_peak value (10^14.2 Hz), while the blue dotted line marks the sources with ν_peak around 10¹⁵ Hz. Some HSPs have a radio-to-X-ray flux ratio lower than f_X/f_r = 9 × 10⁻¹¹, and to select them during the assembling of previous HSP catalogues (1/2WHSP) we applied a selection criterion corresponding to lower ν_peak sources.

The “lost sources” in the upper red dashed region of Fig. 9 are mainly lower ν_peak blazars, while those in the lower region are essentially higher ν_peak blazars. The consequences of the slope criterion applied¹⁴ (Arsioli et al. 2015) are the blue dashed lines displayed, and the estimated lower and upper limits of f_X/f_r with this criterion are 1.490 × 10⁻¹¹ erg cm⁻² s⁻¹ Jy (ν_peak ≈ 10^14.2 Hz) and 2.063 × 10⁻⁹ erg cm⁻² s⁻¹ Jy (ν_peak ≈ 10²⁰ Hz), respectively. The two slope limits are estimated with the assumption that the average X-ray spectral slope is 0.9 and the IR slope is 0.3.

Figure 9 shows that only by setting a double flux limit (both at radio and X-rays) do we manage to select a 3HSP subsample that avoids the red areas, and that is therefore complete. Thus, we built a statistically complete sample by applying both radio and X-ray cuts. Among the radio and X-ray catalogues used for the 3HSP selection, NVSS and RASS have the largest coverage, thus all the statistical tests are based on the subsample of 3HSP sources with both NVSS and RASS counterparts. We built a subsample of sources detected both by NVSS and RASS, which we name 3HSP-NVSS-RASS, to be used for statistical studies. It has a cross-matching radii of 0.3 arcmin and 0.8 arcmin, respectively, and includes a total of 1247 sources.

The sensitivity of a survey, as for the case of the ROSAT survey, may not be the same everywhere on the sky. In this particular case only a small fraction of the sky is actually observed with the highest sensitivity. Here and in the following sections we take this into account by weighing the area of sky available based on the X-ray flux of each source (Voges et al. 1999, 2000).

Figure 10 illustrates the X-ray-weighted synchrotron ν_peak distribution for the 3HSP-NVSS-RASS subsample with different radio flux density cuts. This figure suggests that the 3HSP-NVSS-RASS sample is complete for ν_peak ≳ 10¹⁶ Hz; for all radio flux density cuts the distribution flattens at the low end ν_peak ≲ 10¹⁶ Hz, and for the 3.5 mJy cut the data suggests incompleteness for ν_peak ≲ 10¹⁶ Hz.

Fig. 10. Differential (top) and cumulative (bottom) synchrotron νpeak distribution for the 3HSP-NVSS-RASS sample with different radio flux density cuts. Different colours indicate different radio flux cuts. The dashed lines represent the best fixed-slope linear fitting for each radio-cut subsample.

The cumulative ν_peak distribution, compared with the radio complete Sedentary catalogue (green triangles), are shown in Fig. 10. The dashed lines, representing the standard slope for a non-evolving complete sample, show a likely incompleteness at low ν_peak and low radio flux density values. However, even for a 50 mJy cut the distribution tends to flatten at the low end; there are still several bright radio sources with ν_peak close to the selection threshold that might be missed due to their faint X-ray flux.

Figure 11 shows the distribution of radio flux density with different ν_peak cuts compared to that of the Sedentary survey (Giommi et al. 1999). It looks like the 3HSP-NVSS-RASS subsample is very close to being complete for F_1.4 GHz ≳ 10 mJy; however, compared with the distribution from Sedentary, the slope is flatter. There are still some sources with F_1.4 GHz between 10 and ≈25 mJy lost from the X-ray selected sample, implying that the 3HSP subsamples are complete only for radio flux cuts brighter than 25 mJy. The figure also shows a flattening with higher ν_peak values for faint radio sources; therefore, those relatively high ν_peak sources might be missing due to evolution. Since there was a radio cut already applied for the 3HSP-NVSS-RASS subsample and there was no IR constraint when selecting the 3HSP sources, we do not expect to lose high ν_peak sources.

Fig. 11. Radio flux density distribution with different νpeak cuts.

The average radio flux density with respect to the synchrotron ν_peak is shown in Fig. 12, suggesting that the radio flux density does not depend on the synchrotron peak. Without the radio cut, the average radio flux density is slightly higher on both ends, probably resulting from the incompleteness and evolution of the HSP blazars. For high ν_peak sources, we are probably missing faint radio flux sources because they might evolve negatively, while for sources with ν_peak close to the selection threshold, this is likely due to incompleteness. For the subsample with a radio flux limit of 25 mJy (blue points in Fig. 12) it is clear that the radio flux density is independent of ν_peak, indicating that the fraction of HSP among blazars is independent of radio flux densities.

Fig. 12. Radio flux density vs. νpeak for the entire 3HSP sample (black), and for a complete 3HSP subsample with radio flux cut at 25 mJy (bue).

A good way to check for the completeness of a catalogue is to evaluate the log N − log S. The integral radio log N − log S for the 3HSP-NVSS-RASS subsample with different ν_peak cuts is shown in Fig. 13. Number counts for Sedentary HBL, DXRBS BL Lacs, and DXRBS HBL only are shown as well for comparison. The log N − log S for DXRBS are at 5 GHz; however, given that BL Lacs typically have radio spectra with α_r  ∼  0 (Giommi et al. 2012b; Caccianiga et al. 2001), no conversion between the 5 GHz counts and 1.4 GHz is necessary. The dashed lines in bright bins correspond to a fixed slope of −1.5, the expected value for a complete sample of a non-evolving population in a Euclidean Universe. Since the radio surveys that we use have different sensitivities, here we only considered the subsample of sources included in the NVSS survey with a radio flux density ≥3.5 mJy.

Fig. 13. Radio log N − log S of different samples of blazars. Blue circles represent the 3HSP sample, green triangles the Sedentary sample, red squares the subsample of BL Lacs from DXRBS, red trapezoids the subsample of HSPs from DXRBS, and orange circles the 3HSP subsample with log(νpeak) > 17.

From Fig. 13 we see that the surface density of the 3HSP subsample is approximately a factor of ten higher than that of the Sedentary Survey. This large difference is expected since the latter includes only extreme sources (its ν_peak, distribution peaks at logν_peak ∼ 16.8, compared to logν_peak ∼ 15.5 for the 3HSP sample). Similarly, for the case of DXRBS all BL Lacs outnumber the 3HSPs in every flux bin. For 3HSP subsamples, the higher the ν_peak cut, the lower the density, and the Sedentary number density are consistent with 3HSP subsample with a ν_peak cut at 10^17.2 Hz. The DXRBS HBLs are also in very good agreement with the 3HSP number counts in the region of overlap, which shows that our selection criteria are robust.

Apart from the different normalisation, the log N − log S of every sample or subsample shows a similar trend deviating from the Euclidean slope at radio flux densities lower than ≈25 mJy. The number densities for different BL Lac groups and different ν_peak HBL are almost parallel to each other, implying that the ratio of high ν_peak BL Lacs to low ν_peak BL Lacs remains the same regardless of the radio flux density. It is consistent with Fig. 12 which shows that the average radio flux density does not depend on synchrotron peak frequency and it seems that there is no preference for high ν_peak sources with faint flux and vice versa. This indeed deviates from the prediction of a “blazar sequence” scenario. A clear trend can be seen going from the Sedentary Survey of extreme HSPs to the 3HSP sample, to the entire population of BL Lacs as estimated in the DXRBS survey, with an increase in number by approximately a factor 10 at every step.

The log N − log S indicates that both samples are complete at the bright end, but significantly deviate from the Euclidean slope, −1.5 (dashed line) at the faint end. The Sedentary is a complete sample in the radio band, and the flattening at faint fluxes is due to its cosmological negative evolution, meaning that there are fewer faint radio sources in the young universe. The 3HSP flattening, however, appears to be somewhat stronger than that of the Sedentary Survey, likely because of the onset of some degree of incompleteness at low radio flux densities in addition to the evolutionary effects discussed in the Sedentary Survey paper (Piranomonte et al. 2007).

The faintest bins of the 3HSP-NVSS-RASS subsample have a surface density ≈0.06 deg⁻², indicating that the maximum surface density corresponds to a total of ∼2400 HSP blazars over the whole sky. Given that this number refers only to sources with 1.4 GHz flux density ≥3.5 mJy, and because of the incompleteness discussed above, this can be considered a robust lower limit.

Although the full 3HSP catalogue is not a complete sample in the radio band, we can build a complete subsample from it if we apply specific cuts. For instance, by using a radio and X-ray flux-limited sample and selecting only sources with ν_peak ≳ 10¹⁶ Hz, we built a complete HSP subsample. The subsample with limits in radio flux density, X-ray flux, and ν_peak value, comprising 1096 sources, can be considered highly complete and therefore useful for statistical purposes

6. 3HSP blazars as VHE source candidates

High synchrotron peaked blazars are known to emit in the very high-energy (VHE) γ-ray band. By definition, the synchrotron ν_peak of HSPs are the highest among blazars, sometimes reaching and exceeding 10¹⁸ Hz. Consequently, the peak frequency of the second hump in the SED of these extreme sources goes close to and sometimes well into the TeV band. Synchrotron ν_peak, values ≳10¹⁸ Hz imply that the electrons responsible for the radiation must be accelerated to extremely high energies (Rybicki & Lightman 1986; Costamante et al. 2001).

Considering a simple SSC model where ${\nu }_{\mathrm{peak}}=3.2\times {10}^6{\gamma }_{\mathrm{peak}}^{2}B\delta$ (e.g. Giommi et al. 2012a), assuming B = 0.1 gauss and a Doppler factor δ = 10, HSPs characterised by ν_peak ranging between 10¹⁵ and ≳10¹⁸ Hz require γ_peak to be in the range 2 × 10⁴ − ≳10⁶, where γ_peak corresponds to the Lorentz factor of relativistic electrons emitting at the synchrotron peak frequency.

Sources with ν_peak >  10¹⁷ Hz are often called extreme blazars and are particularly relevant for high-energy astrophysics since they may be the counterparts of VHE γ-ray, sources, high-energy astrophysical neutrinos, and ultra high-energy cosmic rays (UHECR). Table 2 gives the list of all objects in the 3HSP catalogue with rest-frame ν_peak ≥ 10¹⁷ Hz, and with ν_peak ≥ 10¹⁷ Hz and no redshift, which implies that their rest-frame ν_peak is at least 10¹⁷ Hz. In total, there are 384 extreme HSPs in our catalogue, a much larger number than in any previous catalogue. Based on their radio number counts (Fig. 13) we expect an all-sky content of about 370 extreme blazars down to the 3.5 mJy flux limit for the 1.4 GHz channel. Given that the number counts are relatively flat, we estimate that the total number of extreme blazars in the sky is ≈400.

One of the main reasons for assembling the 3HSP catalogue was the need to find sources that could be detected in γ-ray surveys and provide promising targets for VHE/TeV observations. About 50% of the 3HSP sources already appear in one of the existing γ-ray catalogues, sources in the other half are still undetected, but a large fraction of them are expected to be above the sensitivity limit of upcoming VHE telescopes like the Cherenkov Telescope Array (CTA).

An example of this is PGC 2402248, a blazar recently detected by the Major Atmospheric Gamma Imaging Cherenkov Telescope (MAGIC) collaboration (Mirzoyan 2018) that is in the 2WHSP catalogue and is listed in 3HSP as a source with extreme ν_peak = 10^17.9 Hz. The 1BIGB catalogue (Arsioli & Chang 2017) also demonstrates the potential of HSP catalogues built on the basis of multi-wavelength data for the detection and identification of γ-ray, and VHE sources, and for the selection of targets for TeV observations.

A significant fraction of Fermi-LAT detected sources still do not have an assigned counterpart. As blazars are the dominant population of extragalactic persistent γ-ray sources, we checked whether some of the objects in the 3HSP sample could be the counterparts of still non-associated Fermi sources. We found that many 3HSPs may be the counterparts of Fermi 3FHL and 3FGL sources that still do not have a counterpart assigned in the current Fermi-LAT catalogues. In Table 3 we list 19 possible counterparts for these non-associated Fermi 3FHL or 3FGL detections. These 19 sources have no association counterpart in the recent 4FGL catalogue either.

We note that Kaur et al. (2019, hereafter K19) presented the results of an identification campaign of unassociated sources from the Fermi 3FHL catalogue. In Table 3, we indicated those already in K19, and there are ten sources that are not in their identification list. Among their identifications, only 10 out of 110 sources are not included in the 3HSP. We checked all of them, and they are either non-HSP blazars or too uncertain to be included in the 3HSP catalogue. The other 91 sources (i.e. 110 minus 10, which are not in the 3HSP, and 9 already in Table 3) either have an association in the 4FGL catalogue or in the 3FHL catalogue.

We applied the FOM value to estimate the potential detectability of HSPs by current and future VHE-TeV telescopes. FOM is defined in units of 2.5 × 10⁻¹² erg s⁻¹ cm⁻² (Log ν_peakf_νpeak = −11.6), which is the peak flux of the faintest HSP blazar that is included in TeVCat (see Paper I for more details). Figure 14 indicates that the TeV detection rate increases with the FOM between 1 and 10 and that the majority of cases with FOM > 10 have already been reported in TeVCAT. Sources with a higher FOM value (i.e. brighter flux) are expected to be more easily detectable by the current generation of VHE-TeV telescopes. We note that there are still 262 sources with FOM ≥ 1 that are not yet listed in TeVCAT, and 251 of them have already been detected by Fermi as reported in the FGL or FHL catalogues. Those sources should be of great interest to CTA, especially for the planning of observation campaigns.

Fig. 14. Fraction of sources detected in the TeV band for FOM values above the one indicated.

7. Conclusions

We presented 3HSP, the largest compilation of high-synchrotron peaked blazars. 3HSP is an evolution of the 2WHSP (Chang et al. 2017) and 1WHSP catalogues (Arsioli et al. 2015), which were assembled starting from Wise IR data, but were largely incomplete in the radio, X-ray, and γ-ray band. The 3HSP catalogue contains 2013 HSPs and HSP-candidates, with 1618 of them also in the 2WHSP catalogue and 1007 having Fermiγ-ray counterparts. Only 657 3HSP sources are in the 5BZCat (Massaro et al. 2015), implying that the number of known HSP blazars has tripled compared to 2015 when BZCAT was the most complete list of blazars available. Another distinctive aspect of the 3HSP catalogue is that it provides redshift estimates for 88% of the sources, a much higher percentage than in any previous catalogue.

Providing robust candidates for high-energy, VHE/TeV, and γ-ray observations was one of the main motivations for building the 3HSP catalogue. Previous versions of the catalogue have already been used as a seed for VHE or TeV observations, and several new detections from HE or VHE has been secured based on the positions of 2WHSP sources. For example, a new VHE counterpart, PGC 2402248 (2WHSP J073326.7+515355), has recently been detected by MAGIC (Mirzoyan 2018), showing that the 3HSP can contribute potential VHE candidates for future surveys. The 1BIGB catalogue (Arsioli & Chang 2017) is another successful example that demonstrated the presence of new γ-ray sources found based on 2WHSP sources. Therefore, with multi-frequency-based catalogues the search for new VHE sources might become more efficient as 3HSP could point out the possible location of relevant VHE sources that were not known before. We note that 26 3HSP sources have been proposed as counterparts of non-associated Fermi sources (see Table 3).

We also presented several tests to investigate the completeness and the general properties of the catalogue. Our results suggest that the radio and X-ray selected 3HSP subsample is complete for ν_peak >  10¹⁶ Hz. This complete and large (1096 sources) subsample is suitable for the detailed investigation of intrinsic statistical properties associated with HSP and extreme blazars. This will be presented in future publications.

We note that after this paper was completed Paliya et al. (2019) applied a new γ-ray data analysis technique to 337 3HSP sources with ν_peak >  10¹⁷ Hz. They detected 165 objects and report a cumulative signal at > 32σ confidence for the remaining 172 Fermiγ-ray undetected 3HSP sources. Their average spectral slope is very flat (Γ <  2). This justifies the usefulness of the 3HSP catalogue in finding new VHE γ-ray sources.

Acknowledgments

YLC is supported by the Government of the Republic of China (Taiwan), BA is supported by São Paulo Research Foundation (FAPESP) with grant no. 2017/00517-4. PP thanks the SSDC for the hospitality and partial financial support of his visits. CHB would like to thank the Brazilian government and the CAPES foundation for supporting this work under the project BEX 15113-13-2. PG acknowledges the support of the Technische Universität München – Institute for Advanced Study, funded by the German Excellence Initiative (and the European Union Seventh Framework Programme under grant agreement no. 291763) This work was supported by the SSDC, Agenzia Spaziale Italiana Science Data Center; and University Sapienza of Rome, Department of Physics. We made use of archival data and bibliographic information obtained from the NASA/IPAC Extragalactic Database (NED), data and software facilities from the SSDC managed by the Italian Space Agency (ASI).

References

Appendix A: Association of no-radio sources around Fermi detections

There are several sources within the 3HSP catalogue without radio detections. They were found by directly examining the area around Fermi detections reported in the 3FHL catalogue. Section 3 shows that checking γ-ray detected regions is a powerful way to search for blazar candidates. Here we show an example to illustrate how we estimated the synchrotron ν_peak and why these sources are robust HSP-candidates.

First, we discuss the case of 6dF J0213586−695137 (see Sect. 2 and Fig. 2). We carefully checked the error region of the Fermi source detected as part of the 3FHL, 8-year catalogue (4FGL) and 9 yr MST (Campana et al. 2018) lists and we confirm that 6dF J0213586−695137 is the only plausible counterpart (see Fig. A.1). No Fermi blazars emit in the radio or X-rays, so finding a Fermi counterpart can begin by checking radio sources and X-ray sources around it. Figure A.1 shows the positions of the γ-ray detections (from 3FGL, 3FHL, and 4FGL, with the error circle referring to the last: purple triangles) and of the X-ray detections of 6dF J021358.6−695137 (from XMM-Newton slew and XRT: blue circles). The corresponding UV (GALEX), optical (GAIA), and IR (AllWISE) counterparts are all within ≈5 − 10 arcsec from the X-ray positions. Based on the available information 6dF J021358.6−695137 is the only source within the γ-ray error ellipse with a blazar-like optical/UV/X-ray SED consistent with the γ-ray data, and therefore the very likely counterpart of the Fermiγ-ray detection.

Fig. A.1. Region around 4FGL J0213.8−6949. Purple triangles represent γ-ray detections (3FHL, 4FGL, and 9 yr MST) and blue circles indicate X-ray detections. 6dF J0213586−695137 is the only X-ray source inside the Fermi detected region.

Figure A.2 is the SED of 6dF J021358.6−695137. The SED is fully consistent with that of an HSP blazar, such as a variable and flat X-ray spectrum and flat γ-ray spectrum with Compton dominance < 1. If the X-ray flux were to be attributed to a nearby AGN, then, based on the observed distribution of optical-UV/X-ray flux ratio (α_oX) of radio quiet AGN, this hypothetical object should be about 50 times brighter than 6dF J021358.6−695137 in the optical and UV bands. No such object is present in (or near) the X-ray error circles, so we can exclude that the X-ray source is a background AGN.

Fig. A.2. Spectral energy distribution of 6dF J021358.6−695137. The red filled circles represent the average SED of MKN421 scaled down to the flux level of 6dF J021358.6−695137.

We re-scaled the SED of Mrk 421, and plot the re-scaled fluxes of Mrk 421 in Fig. A.2 (brown points), suggesting that the shape of the two SEDs (6dF J021358.6−695137 and Mrk 421) are very similar. The extrapolation of the UV emission in the SED closely matches the X-ray spectrum. Moreover, we do not claim that this is a radio quiet HSP since the radio upper limit in its SED is consistent with the flux predicted by the SED of MKN421 (rescaled to the flux of this source, and shown as a brown dotted line). The lack of radio detection in the 3HSP source is due to insufficient sensitivity of the SUMSS survey rather than to an usually low radio emission in the object.

In conclusion, 6dF J021358.6−695137 (3HSP J021358.6−695137) is the counterpart of the Fermi source and is clearly an extreme HSP blazar. This source is not found by chance or just because it is inside the Fermi detected region.

There are only 24 sources in the 3HSP catalogue that do not have a radio counterpart, and we double checked these source one by one again. We list all of them in Table 4 now. Among these sources, there are only three sources that do not have a UV counterpart yet. All three of these three sources have well-described X-ray spectral data, and we could fit the synchrotron peak from it. The other “no-radio” sources all have IR, UV, optical, and X-ray counterparts, and we could clearly tell their peak frequency. After checking all of them, we suggest that these sources are plausible HSPs, even though they have no radio counterparts.

Here we show another example, namely 3HSP J032852.6−571605 (Figs. A.3 and A.4).

Fig. A.3. Region around 4FGL J0328.8−5715. Purple triangles represent γ-ray detections (3FHL, 4FGL, and 1BIGB), blue open circles indicate X-ray detections, while the red filled circle denotes an unrelated radio source outside the γ-ray error circle. 3HSP J032852.6−571605 is the only X-ray source inside the Fermi region. There are three γ-ray detections around this 3HSP sources, and 3HSP J032852.6−571605 is the only plausible counterpart for them.

Fig. A.4. Spectral energy distribution of 3HSP J032852.6−571605. The Fermi spectrum is fully consistent (in intensity and slope) with the synchrotron part of the SED. The X-ray and the UV data imply a νpeak ≈ 1017.5 Hz.

We conclude that those sources without a radio counterpart in our catalogues and inside Fermi-LAT error ellipses are genuine HSP-candidates. We did not select them because they are inside the Fermi error circles or with X-ray detections, and every source in the 3HSPs was carefully examined. Keeping the focus on the aim of this catalogue, which is to provide seed sources for VHE observations, we include all promising HSP and HSP-candidates based on the currently available information. Thus, we could not exclude those sources.

In addition, these sources were not found by chance. According to the X-ray logN-logS, there might be a surface density of 0.06–0.07 deg⁻² HSPs. As the log N − log S suggests, we might find by chance one HSP within 14–15 deg⁻², which is much larger than the area covered by the Fermi error circle. The largest position error area in the Fermi catalogues is around 1°. In summary, the reported HSPs are reliable counterparts for the γ-ray signatures, in agreement with their γ-ray error circles, as reported by the Fermi team. This is so for all the γ-ray counterparts reported in this work.

All Tables

All Figures

	Fig. 1. Spectral energy distributions of 5BZB J0403−2429 (top) and 5BZG J0903+4055 (bottom). The sources were not selected for 2WHSP due to the absence of IR data (top) and the slope criterion (bottom).
In the text

	Fig. 2. Spectral energy distribution of 6dF J0213586−695137. This source was missed by 2WHSP because there is no radio data available for the source.
In the text

	Fig. 3. Distribution of the new sources and those in the 2WHSP catalogue (νpeak, top; radio flux density, bottom) normalised according to the number of objects in each sample.
In the text

	Fig. 4. Distribution of the new sources and the 2WHSP (X-ray–IR slope, top; radio–IR slope, bottom) normalised based on the total number of objects in each sample.
In the text

	Fig. 5. Example of photometric redshift (photo-z) estimation. The SDSS Dr14 optical spectrum (top) and SED (bottom) of the source 3HSP J154433.1+322148. The red line represents the non-thermal component, the green line the giant elliptical template fitted to the IR (and partly optical) data at z = 0.32, and the blue line the data extract from the SDSS DR14 optical spectrum.
In the text

	Fig. 6. Comparison between photo-z and spectral-z. The central solid line has slope of one and represents the equality between the two redshift estimations, while the red, green, and blue dashed lines indicate a scatter of 0.1, 0.2, 0.3, respectively.
In the text

	Fig. 7. Redshift distribution of the 3HSP sample. The red line represent all the 3HSP sources with a redshift estimation, while the blue line represents sources with an estimation, but not lower limits.
In the text

	Fig. 8. Aitoff projection of 3HSP sources in Galactic coordinates. Red, blue, and green points represents FOM >  1.0, 0.2 <  FOM <  1, and FOM <  0.2, respectively (see text for details).
In the text

	Fig. 9. Sampling of the radio flux density–X-ray flux plane with the 3HSP, Sedentary, and DXRBS samples. The blue dashed and dotted lines are the slope limits for the 3HSP sample. The black dotted line is the flux ratio that separates LBL and HBL. The green and yellow dotted lines denote the forbidden region for BL Lacs and FSRQs. The red dashed lines indicate the incomplete regions of 2WHSP.
In the text

	Fig. 10. Differential (top) and cumulative (bottom) synchrotron νpeak distribution for the 3HSP-NVSS-RASS sample with different radio flux density cuts. Different colours indicate different radio flux cuts. The dashed lines represent the best fixed-slope linear fitting for each radio-cut subsample.
In the text

	Fig. 11. Radio flux density distribution with different νpeak cuts.
In the text

	Fig. 12. Radio flux density vs. νpeak for the entire 3HSP sample (black), and for a complete 3HSP subsample with radio flux cut at 25 mJy (bue).
In the text

	Fig. 13. Radio log N − log S of different samples of blazars. Blue circles represent the 3HSP sample, green triangles the Sedentary sample, red squares the subsample of BL Lacs from DXRBS, red trapezoids the subsample of HSPs from DXRBS, and orange circles the 3HSP subsample with log(νpeak) > 17.
In the text

	Fig. 14. Fraction of sources detected in the TeV band for FOM values above the one indicated.
In the text

	Fig. A.1. Region around 4FGL J0213.8−6949. Purple triangles represent γ-ray detections (3FHL, 4FGL, and 9 yr MST) and blue circles indicate X-ray detections. 6dF J0213586−695137 is the only X-ray source inside the Fermi detected region.
In the text

	Fig. A.2. Spectral energy distribution of 6dF J021358.6−695137. The red filled circles represent the average SED of MKN421 scaled down to the flux level of 6dF J021358.6−695137.
In the text

Fig. A.3. Region around 4FGL J0328.8−5715. Purple triangles represent γ-ray detections (3FHL, 4FGL, and 1BIGB), blue open circles indicate X-ray detections, while the red filled circle denotes an unrelated radio source outside the γ-ray error circle. 3HSP J032852.6−571605 is the only X-ray source inside the Fermi region. There are three γ-ray detections around this 3HSP sources, and 3HSP J032852.6−571605 is the only plausible counterpart for them.

In the text

	Fig. A.4. Spectral energy distribution of 3HSP J032852.6−571605. The Fermi spectrum is fully consistent (in intensity and slope) with the synchrotron part of the SED. The X-ray and the UV data imply a νpeak ≈ 1017.5 Hz.
In the text