BLAZAR RADIO AND OPTICAL SURVEY (BROS): A CATALOG OF BLAZAR CANDIDATES SHOWING FLAT RADIO SPECTRUM AND THEIR OPTICAL IDENTIFICATION IN PAN-STARRS1 SURVEYS

To select flat-spectrum radio sources, we use the TGSS-NVSS spectral index maps and catalog (TGSS-NVSS spindx/catalog, de Gasperin et al., 2018), which contains about 0.5 million sources. Fluxes and the spectral index are taken from 0.15 GHz TGSS and 1.4 GHz NVSS catalogs over the wide region of $\delta>-40^{\circ}$. In the TGSS-NVSS catalog, the images of the NVSS and TGSS are re-processed, combined, and then the sources are extracted. Therefore, the TGSS-NVSS catalog differs slightly from the original TGSS and NVSS catalogs (see de Gasperin et al., 2018, for details). It is estimated that the typical background root mean square (RMS) noise of the NVSS-TGSS catalog is $\sim 0.45$ mJy beam^-1 at 1.4 GHz and $\sim 8$ mJy beam^-1 at 0.15 GHz. The RMS noise value at 0.15 GHz is slightly higher than the value in the original TGSS catalog, as TGSS images are convolved to adjust the resolution for NVSS images in the TGSS-NVSS catalog. Mean estimated astrometric error is 3.5 arcsec for the TGSS-NVSS catalog.

From the TGSS-NVSS catalog, we excluded sources located within the galactic plane of $|b|<10^{\circ}$ where galactic objects are dominant. Then, we selected sources with the source classification flag of ’S’ in the TGSS-NVSS catalog, corresponding to the sources detected both in TGSS and NVSS with no other sources in the search region (see Section 2, de Gasperin et al., 2018, for details).

Among the rest, we selected flat-spectrum and compact objects, i.e., blazar candidates. To define the spectral flatness and the compactness of blazars in this research, we examined known blazars detected by Fermi/LAT. Figure 1 shows a scatter plot of the compactness and spectral index of all TGSS-NVSS sources and gamma-ray blazars (4LAC; The Fermi-LAT collaboration, 2019b), as shown in the TGSS-NVSS catalog. The spectral index $\alpha$ is defined as $F_{\nu}\propto\nu^{\alpha}$ between the 0.15 and 1.4 GHz. The compactness $C$ is defined as $C\equiv F_{{\rm integral}}/F_{{\rm peak}}$, where $F_{{\rm integral}}$ and $F_{{\rm peak}}$ are the integrated and peak flux measured using TGSS, respectively. We introduce $C$ to avoid contamination of radio galaxies (misaligned blazars). As shown in the figure, 94% of Fermi-LAT gamma-ray blazars are located in $C\leq 1.6$ and $\alpha\geq-0.6$. Thus, we selected sources with those criteria for this research.

After applying the selection above, 88,211 BROS sources were left. From Figure 1, 58 Fermi-LAT gamma-ray blazars which include several famous blazars such as 1ES 1215+303 ($z=0.131,C=2.0$) and OQ 530 ($z=0.151,C=1.9$) were excluded. Twenty-four of these blazars have $C>1.6$ and mostly reside nearby ($z<0.5$), compared to samples with $C<1.6$. These sources are actually resolved in TGSS images. We also note that a small fraction $\lesssim 0.5$% of gamma-ray blazars have steep spectrum lobes at low frequency (Giroletti et al., 2016). These populations are missed in the BROS catalog.

To search for the optical counterparts of the selected flat-spectrum radio sources, we utilized the Pan-STARRS photometric catalog²²2available at http://mastweb.stsci.edu/ps1casjobs/home.aspx (Chambers et al., 2016). The covered region of the PS1 data ($3\pi$ sky above $\delta>-30^{\circ}$ ) matches that of the TGSS-NVSS catalog. The PS1 data were taken between 2010 and 2012 with a 1.8 m telescope near the summit of Haleakala in Hawaii.

Using the PS1 photometric catalog, we identified the optical counterparts of BROS sources by finding the closest PS1 neighbor within a tolerance radius $r_{\rm tol}$. First, we collected all PS1 sources detected in $r$ band within 20 arcsec of the selected flat-spectrum radio sources, which is large enough compared to positional uncertainty in the radio band (typically a few arcseconds).Figure 2 shows the number count of all PS1 neighbors to BROS sources in an annulus between $r_{PB}$ and $r_{PB}+dr_{PB}$ (we set $dr_{PB}=0.066$ arcsec to reduce uncertainties of each bin less than 5%). The number count peaks at 0.5 arcsec and increases again toward larger radii. The first peak is mostly dominated by a true-matched source, while the latter increse is caused by unrelated sources randomly coincident, $\propto\pi nr_{PB}dr_{PB}$, where $n$ is the number density of PS1 sources and $r_{PB}$ is the separation of neighbors with respect to a BROS source. The green line in Figure 2 indicates the linear line representing the association by chance between random positions and PS1 sources. This line is a measure of ”the background”. The tolerance radius $r_{\rm tol}$ is determined empirically to minimize contamination by chance as well as to maximize the true-matched sources. We adopted $r_{\rm tol}$ = 2 arcsec. Finally, we obtained 42,757 BROS optical counterparts. From Figure 2, we estimated the contamination ratio $\sim$ 10% and overlooked ratio $\sim$ 19% for our data selection.

Gigahertz peaked-spectrum (GPS) and compact steep-spectrum (CSS) radio sources are a kind of compact AGN. Given that we selected sources using only two bands (0.15 and 1.4 GHz), the sources with a spectral peak at $\sim$ GHz can be misidentified as flat spectrum objects in the BROS.

To investigate the contamination fraction of GPS and CSS radio sources, we checked the multi-band spectrum in the radio band, taking the cross-match between BROS and the multi-frequency radio source catalog by Kimball & Ivezic (2014, hereafter, Kimball catalog), which contains flux values from GB6 at 4.85 GHz and WENSS at 325 MHz (Rengelink et al., 1997; Gregory et al., 1996). This enables the blazars and GPS/CSS sources to be distinguished by checking the existence of a peak at $\sim$ GHz.

Figure 3 shows examples of GPSs and blazars. We estimated the expected GB6 and WENSS flux densities using only NVSS flux and TGSS-NVSS spectral index (i.e., assuming a single power-law spectrum) and comparing it with the flux reported in the Kimball catalog. $\sim 10,000$ BROS sources have a counterpart in the GB6 and/or WENSS catalogs. About 25% of them are identified as GPSs whose source flux is deviates from the flat spectrum at the $\geq$5 $\sigma$ confidence level. For those objects, we added a flag of ”G” to the ”SourceFlag” column. However, the GB6 and WENSS catalogs do not cover all-sky condition and have different sensitivity compared to the TGSS-NVSS catalog. Therefore most of the BROS samples ($\sim 70,000$ sources) do not have counterparts in the Kimball catalog. These sources are treated as ”blazars” in the BROS catalog; however, the GPS/CSS sources may be contaminated.

A description of the BROS catalog is summarized in Table 1. The catalog is available at http://www.bao.city.ibara.okayama.jp/BROS/ and the Centre de Données astronomiques de Strasbourg https://cds.u-strasbg.fr/.

Figure 4 shows the sky distribution of BROS sources. The distribution is not perfectly uniform. A higher source density region (strip) is present, for example, at (l, b) = ($0^{\circ},90^{\circ}$) to (l, b) = ($60^{\circ},10^{\circ}$). We think that this is not intrinsic, but instead due to the nonuniform sensitivity of the TGSS (i.e., the nonuniform distribution of TGSS sources), as evidenced by the spatial distribution of TGSS sources in Fig. B.2 of Intema et al. (2017b).

CRATES is a well-known blazar candidate catalog (Healey et al., 2007); here, they selected radio sources of $S_{4.8{\rm GHz}}\geq 65$ mJy showing flat radio spectra of $\alpha\geq-0.5$ between 1.4 GHz (or 0.84 GHz for the sources located at $\delta<-40^{\circ}$) and 4.8 GHz. CRATES has 11,131 sources and is, thus, a quarter size of BROS. One major difference with respect to BROS is the radio frequencies. Given that CRATES utilizes the NVSS catalog for 1.4 GHz and GB6, PMN, and S5 catalogs for 4.8 GHz, the detection limit of CRATES depends on its declination with a typical detection limit at $\delta<+75^{\circ}\sim 65$ mJy at 4.8 GHz. Thus, the CRATES catalog is shallower than BROS.

The Candidate Gamma-Ray Blazar Survey catalog (CGRaBS Healey et al., 2008) is essentially based on the CRATES and contains 1,625 sources with radio and X-ray data. CGRaBS also contains red-shift and blazar-type classifications (FSRQ or BL Lac) based on optical spectroscopy. The Roma-BZCAT catalog (Massaro et al., 2015) contains 3,561 blazars selected by multi-frequency surveys. In the Roma-BZCAT, most of the blazars are classified into two types: flat-spectrum radio quasars (FSRQs) and BL Lac objects (BL Lacs), by optical spectroscopic data. 3HSP catalog (Chang et al., 2019) contains 2,013 high synchrotron peaked blazar candidates made for high energy astrophysics and neutrino astronomy. CGRaBS, Roma-BZCAT and 3HSP focus on multi-wavelength observations and are useful for identifying blazar properties; however, the number of sources is smaller than that of BROS.

Table 2 shows a summary of blazar catalogs. Hereinafter, we focused on a comparison with CRATES, which adopts similar selection criteria to BROS. We cross-matched the BROS and CRATES catalogs. First, we selected CRATES sources located at $\delta>-40^{\circ}$; the resulting CRATES source number was 9,071. For these sources, we searched for the TGSS counterpart. Then, we identified 7,971 sources (88%) with TGSS-NVSS counterparts; additionally, 5,178 sources (57%) have BROS counterparts. The sources that have TGSS-NVSS counterparts but are not listed in the BROS catalog, are due to the soft radio spectra of $\alpha<-0.6$ or have a non source classification flag of ’S’ or a radio compactness of $C>1.6$.

We also checked the reason why about 12% of the CRATES sources are not detected by the TGSS-NVSS catalog. One reason would be due to the sensitivity of the TGSS-NVSS catalog at 0.15 GHz. The typical RMS noise is 8 mJy beam^-1 at 0.15 GHz (de Gasperin et al., 2018); sensitivity is not uniform. In the CRATES catalog, 4.8 GHz fluxes of the sources are larger than 65 mJy. Thus, if a source has a very hard spectrum of say 0.6, the source cannot be detected by the TGSS. To check this possibility, we made a histogram of $\alpha_{1.4-4.8{\rm GHz}}$ for CRATES sources that are listed or dropped in the TGSS catalog (Figure 5). We find that the CRATES sources not detected by the TGSS-NVSS show a relatively hard spectral index compared to the TGSS-CRATES sources. In addition, contamination of GPS/CSS sources in the CRATES catalog may also be a reason for non-detection sources at 0.15 GHz, because it has a spectral break peaking for the GHz band.

1861 Roma-BZCAT sources have BROS counterparts within the BROS sky region, in other words, 1253 ($\sim 40$%) of the Roma-BZCAT objects are missed in the our catalog. Among the missing Roma-BZCAT sources, 67% are recorded as either of ’BL Lac’, ’BL Lac Candidate’, ’BL Lac-galaxy dominated’, ’Blazar Uncertain type’, while the remaining objects are classified as ’QSO RLoud flat radio sp’. Figure 6 shows the redshift distribution of the missing ratio of the Roma-BZCAT objects in BROS. The bump at low redshift ($z\lesssim 0.5$) in the Figure 6 is mostly dominated by the BL Lac sources, while at $z>0.5$ the missing population is dominated by FSRQs. The redshift measurements of BL Lac sources are known to be biased towards low redshifts $\sim 0.5$ because of spectroscopic limitations (see Ajello et al., 2014, Figure 7 for detail). We further investigate the reason for the differences between BROS and Roma-BZCAT. About 70% sources of the missed Roma-BZCAT sources are not listed or having source classification flag of ”L”, i.e., non-TGSS detection, (see de Gasperin et al., 2018, for detail) in the TGSS-NVSS catalog. By estimating the 0.15 GHz flux, we confirm that the corresponding flux of the missed Roma-BZCAT sources are below the detection limit of TGSS, where we extrapolate from the NVSS flux assuming a spectral index of 0.5. Therefore, we conclude the TGSS flux limit causes this mismatch. The remaining 30% are excluded by its radio morphology.

The similar situations are also realized when comparing the BROS catalog with the 3HSP catalog (363 missing sources, 21.5%) and 4LAC (784 missing sources, 32%). The BROS catalog is uniform over the footprint, but tends to miss the sources that have hard radio spectrum limited by TGSS sensitivity at 0.15 GHz.

Here we discuss the optical properties of BROS sources. Figure 7 upper and lower panel shows the $g-r$ vs. $r$ color-magnitude and $g-r$ vs. $r-i$ color-color diagrams for the BROS sources, respectively. Magnitudes are in Kron magnitude, which is photon flux measured with the Kron technique for an extended source (Kron, 1980), because many BROS sources show extended optical emission, i.e., host galaxies. We exclude sources with large photometric uncertainties (g-error, r-error or i-error $>$ 0.3 mag.) from the plot.

The upper panel in Fig. 7 clearly shows the presence of two populations: one cluster at $(g-r,r)\sim(0.2,18.0)$ to $(0.2,22.0)$ and the other at $(g-r,r)\sim(0.8,18.0)$ to $(1.6,22.0)$. The latter population is not seen in CRATES catalog (see Appendix A for details about the CRATES-PS1 sources and their comparison with BROS-PS1 sources). The horizontally elongated distribution at $(g-r,r)\sim(0-1.8,22.0)$ shown in the upper panel is not intrinsic, but is due, instead, to the large error associated with the detection limit of the PS1 catalog.

To understand the nature of these two populations, we calculated the expected gri magnitudes and colors, utilizing phenomenological spectral energy distribution (SED) templates of a composite of blazars and elliptical galaxies. Blazar’s broadband SED consists of two spectral peaks (synchrotron and, most likely, inverse Compton emissions) and that the fainter blazars show higher peak frequency due to inefficient radiative cooling, so-called blazar sequence (e.g., Fossati et al., 1998; Kubo et al., 1998; Ghisellini et al., 2017). We adopted the latest blazar sequence model in Ghisellini et al. (2017). As described later, this model well explained its magnitude-color diagrams for the BROS sources. However, we note that the existence of the blazar sequence is still under debate because of selection effect (e.g., Caccianiga & Marchã, 2004; Padovani et al., 2007; Giommi et al., 2012; Giommi & Padovani, 2015; Padovani et al., 2019). For elliptical galaxies, we used the stellar population synthesis model by Bruzual & Charlot (2003) assuming the Salpeter initial mass function, instantaneous star formation, age of 10 Gyr old, and solar metallicity. Figure 8 shows the overall SED with bright host galaxy SED templates of various stellar mass.

Figure 9 shows the $g-r$ vs $r$ color-magnitude and $g-r$ vs $r-i$ color-color diagrams of blazars with several gamma-ray luminosity $(L_{g}$ and host galaxy stellar mass $M_{\rm host})$, at various redshifts. BROS sources are also shown. The model curves based on our simple calculation closely reproduce the observed population, particularly $(L_{g},M_{\rm host})=(10^{43}~{}{\rm erg\ s^{-1}},10^{12}M_{\odot})$ and $(10^{44}~{}{\rm erg\ s^{-1}},10^{10}M_{\odot})$. The optical emission is dominated by bright hosts in the former model, whereas in the latter case, the emission is dominated by blazars. Thus, the right population is dominated by bright elliptical galaxy ($z<0.5$) emission, whereas the left population is dominated by blazar emission. Large samples (13,153 sources, 52%) which have bright host galaxy are newly discovered by BROS, while hundreds of similar samples are reported in Roma-BZCAT, 3HSP and other catalogs (e.g., 1,207 sources, 19% are classified as latter population in CRATES catalog, see Appendix A). We also note that ’big blue bump’ component originated in the disk is confirmed in several FSRQs in the optical band. For reference, the colors-magnitude and color-color diagram of normal SDSS quasar located at low redshift ($z<1$) are shown in Figure 9 as yellow contour (Pâris et al., 2017). As we see in the figure, the observed low redshift SDSS quasars which do not deviate from the simple Blazar model confirm that the contribution from the disk component is indistinguishable in our distinction.

Hereinafter, for the sake of simplicity, we divided the BROS sources into two populations; namely quasar-like ($r>8\times(g-r)+13$ in $g-r$ vs $r$ plot) and elliptical-like populations ($r<8\times(g-r)+13$ in $g-r$ vs $r$ plot). In the BROS catalog, if the sources have a PS1 candidate and are not classified as GPS sources, we set ’Q’ (’E’) flag for the source which locates in a quasar-like (elliptical-like) population to the ’SourceFlag’ column. The sources that do not have a PS1-candidate and are not classified as GPS sources were assigned a ’B’ flag.

We note that there is also a marginal population in the valley between quasar-like and elliptical-like populations. This marginal population is likely a mixture of the distribution tails of both population.. Further detailed studies such as the convolution of blazar and galaxy luminosity functions will be able to elucidate the properties of the marginal population; however, this is beyond the scope of this paper.

We also investigated the spatial extension of sources in the optical band. We use the $iPSFMag-iKronMag>0.05$ formula, where $iPSFMag$ is the photon flux derived with a point source spread function and $iKronMag$ is the photon flux for the extended source, to select the extended source. As a result, $\sim 93$% elliptical-like sources were classified as extended sources; in contrast, $\sim 22\%$ quasar-like sources were classified as extended sources. This result also supports an ”elliptical-like” branch, determined by its optical color, consisting of elliptical galaxies.

Generally, optical emission in FSRQs is not dominated by a host galaxy, as FSRQs have strong jet emission ($L_{g}>10^{46}$ erg s^-1) compared to that of BL Lac-type objects (see Fig. 8). Therefore, it is expected that the elliptical-like population consists mostly of BL Lacs, whereas the quasar-like population consists of both FSRQs and BL Lacs. We also confirmed this by cross-matching between BROS and classified sources (BL Lac or FSRQ) in 4LAC. As a result, 74% (26%) of elliptical-like sources are classified as BL Lacs (FSRQs) in 4LAC, and 44% (56%) of quasar-like sources are classified as BL Lacs (FSRQs).

Figure 10 shows the $g-r$ vs. NVSS flux and $g-r$ vs. radio loudness $R$ (defined as $F_{{\rm 1.4~{}GHz}}/F_{r}$) of the BROS sources together with the magnitude and color diagram. As shown in the upper right panel of Figure 10, quasar-like sources show slightly larger NVSS fluxes compared to the elliptical-like one.

The quasar-like sources tended to have a high $R$ compared to that of elliptical-like ones. Radio loudness $R$ is a useful tool to distinguish the spectral type of blazar, as a high $R$ value is expected for BL Lac-type objects from the unification model for AGN (e.g., Urry & Padovani, 1995a). However, in our sample, optical emission of elliptical-like sources are likely contaminated by the host galaxy, such that our $R$ value tends to be small. Therefore, we believe that the differences in $R$ between quasar-like sources and elliptical-like ones are not due to the spectral difference originating from jet emission.

Figure 11 shows the $\log N$-$\log S$ distributions for the BROS sources based on TGSS and NVSS fluxes. We fit the $\log N$-$\log S$ data points using a power-law function of $N\propto S^{\beta}$. The NVSS-TGSS catalog is estimated to be complete down to $\sim 5$ mJy for NVSS and $\sim 100$ mJy for TGSS sources (de Gasperin et al., 2018). To avoid the effects of source incompleteness for a lower radio flux, we used sources with $\geq 500$ mJy for the derivation of the slope. We obtain the power-law indices of $-1.62\pm 0.03$ and $-1.74\pm 0.03$ for TGSS and NVSS fluxes, respectively. The power-law index becomes $-1.5$ if the space density of objects is constant in Euclidean space (e.g., Peterson, 1997). The derived indices steeper than -1.5 indicate positive evolution for the selected flat-spectrum radio sources, namely the space density in the past is higher than that at present.

The $\log N$-$\log S$ plots for the quasar-like and elliptical-like BROS sources are displayed in Figure 12. The quasar-like and elliptical-like sources show power-law indices of $-1.78\pm 0.07$ and $-1.49\pm 0.05$, respectively, indicating a significant difference between the two populations. The former population shows positive evolution whereas the latter does not, supporting our interpretation that the former is dominated by quasars and the latter consists of a different type of object relatively nearby elliptical galaxies.