New Identifications and Multi-wavelength Properties of Extragalactic Fermi Gamma-Ray Sources in the SPT-SZ Survey Field

The fourth Fermi Large Area Telescope (LAT) catalog (4FGL) contains 5064 $\gamma$-ray sources detected at high significance, but 26% of them still lack associations at other wavelengths. The SPT-SZ survey, conducted between 2008 and 2011 with the South Pole Telescope (SPT), covers 2500 $\mathrm{deg^2}$ of the Southern sky in three millimeter-wavelength (mm) bands and was used to construct a catalog of nearly 5000 emissive sources. In this study, we introduce a new cross-matching scheme to search for multi-wavelength counterparts of extragalactic $\gamma$-ray sources using a mm catalog. We apply a Poissonian probability to evaluate the rate of spurious false associations and compare the multi-wavelength associations from the radio, mm, near-infrared, and X-ray with 4FGL $\gamma$-ray sources. In the SPT-SZ survey field, 85% of 4FGL sources are associated with mm counterparts. These mm sources include 94% of previously associated 4FGL sources and 56% of previously unassociated 4FGL sources. The latter group contains 40 4FGL sources for which SPT has provided the first identified counterparts. Nearly all of the SPT-associated 4FGL sources can be described as flat-spectrum radio quasars or blazars. We find that the mm band is the most efficient wavelength for detecting $\gamma$-ray blazars when considering both completeness and purity. We also demonstrate that the mm band correlates better to the $\gamma$-ray band than the radio or X-ray bands. With the next generation of CMB experiments, this technique can be extended to greater sensitivities and more sky area to further complete the identifications of the remaining unknown $\gamma$-ray blazars.


INTRODUCTION
The γ-ray background holds important clues to the nature of galaxy evolution, the cosmic history of black hole accretion, and possibly the nature of dark matter (Ajello et al. 2015;Funk 2015). Unveiling the nature of the unassociated F ermi Large Area Telescope (LAT) sources is one of the biggest challenges in γ-ray astronomy, due to the relatively large point spread function, and is necessary to understand the contribution of various source classes to the γ-ray background.
Blazars are a subclass of active galactic nuclei (AGN) with relativistic jets of high-energy particles pointing near our line of sight (e.g., Urry & Padovani 1995). Their non-thermal emission is generally detected across the entire electromagnetic spectrum from radio to γ-ray bands. Blazars are sub-classified into flat spectrum radio quasars (FSRQs) and BL Lac objects (BL Lacs), according to the equivalent width of the emission lines in their optical spectrum (Stickel et al. 1991;Stocke et al. 1991;Marcha et al. 1996). These two subclasses of blazars are thought to be intrinsically different, perhaps based on their accretion mode (Dermer & Giebels 2016). FS-RQs have high luminosity and a thin and radiatively efficient black hole accretion disk (Malkan & Moore 1986), while BL Lacs are powered by an advection-dominated, low-radiative-efficiency accretion flow (Dermer & Giebels 2016;Blandford et al. 2019). The jet emission is relativistically beamed (Ghisellini 2019), with a Doppler boosting factor corresponding to a bulk Lorentz factor of several to greater than ten (Pushkarev et al. 2009). In both cases, the broad-band spectra consist of two broad humps, one peaking in the IR-to-X-ray regime, and the other peaking in the γ-ray regime. The low energy peak is believed to be due to synchrotron emission, while the high energy peak is likely due to inverse Compton scattering of lowenergy photons of either the same synchrotron photons (for BL Lacs) or external photons from the disk/BLR (for FSRQs) (e.g. Dutka et al. 2017). However, some blazars might not necessarily be detected in γ-rays (e.g. Paliya et al. 2017). Indeed, a recent study showed that blazars undetected in γ-rays are likely to have relatively smaller Doppler factors and more disk dominance (Paliya et al. 2017). In the case of strong Compton scattering, the beaming of γ-rays could be larger than e.g. that seen in the radio (Dermer 1995) leading to the possible nondetection (or reduced detection efficiency) of γ-rays from sources not seen exactly pole-on.  (Everett et al. 2020). In this region, grey points mark all SPT sources in Everett et al. (2020); green circles represent the position of the previously unassociated 4FGL sources which have SPT counterparts (40); blue circles represent the already associated 4FGL sources (211); red circles represent the remaining still-unidentified 4FGL sources (31).
The F ermi Gamma-ray Space Telescope has been observing the high-energy sky from 50 MeV to 1 TeV since 2008 (Atwood et al. 2009). The fourth F ermi LAT catalog (4FGL) of γ-ray sources is based on observations for eight years, and is a deep all-sky survey in γ-ray bands (Abdollahi et al. 2020). The F ermi 4FGL catalog contains 5064 γ-ray sources detected with at least 4σ significance. More than 62% of these sources are associated with AGN. Currently, 26% of F ermi 4FGL sources are unassociated. The high-galactic latitude sources may fall in the following categories: blazars, radio galaxies, and milli-second pulsars (e.g. Abdollahi et al. 2020;Schinzel et al. 2017), with the largest fraction of high-Galacticlatitude (and presumably extra-galactic) sources being blazars, where millimeter-wave (mm) and γ-ray emission are likely to be produced co-spatially in the extremely compact emission regions within the jet (Meyer et al. 2019). This makes the mm regime a very efficient way to identify both blazars and previously unassociated F ermi sources.
Previous studies have targeted monitoring of F ermi γ-ray blazars at various wavelengths, based on, for example, follow-up at radio wavelengths (Schinzel et al. 2015;Schinzel et al. 2017) or the infrared behavior of γ-ray sources with WISE (D' Abrusco et al. 2012;Massaro & D'Abrusco 2016). From these studies, blazars are known to reside at high redshifts (z > 0.1) and exhibit extreme apparent luminosities with strong variability.
The South Pole Telescope (SPT, Carlstrom et al. 2011) is a 10-meter telescope dedicated to studying the cosmic microwave background (CMB). The SPT has been used to survey thousands of square degrees of the southern extragalactic sky at 1.4, 2.0, and 3mm with arcminute reso-lution down to milli-Jansky noise levels. In this work, we focus on the 2500 deg 2 SPT-SZ survey, conducted with the SPT from 2008-2011. The SPT-SZ field is at high Galactic latitude (|b| > 15), and thus, most of the sources are expected to be of extra-galactic origin (Everett et al. 2020). Roughly 3500 synchrotron-dominated sources are detected at high significance in the SPT-SZ maps, providing a powerful tool to identify unassociated 4FGL sources. This wavelength regime is particularly suited to uncovering the extragalactic unidentified F ermi 4FGL sources, because the mm sources are almost exclusively (>80%) blazars, and blazars are the dominant population in the γ-ray regime. The instrument sensitivity, survey area, and resolution of the SPT are well suited to this task.
In this paper, we establish the methodology of associating mm point sources from CMB surveys with γ-ray sources. In Section 2 we present the data used in this study, while Section 3 describes the association method. Section 4 discusses the implications for the associations, the multi-wavelength properties of the new associations, and the comparisons to previous work. In this analysis, we use a flat ΛCDM cosmology with Ω Λ = 0.73, Ω M = 0.27, and H 0 = 71 km s −1 Mpc −1 .

OBSERVATIONAL DATA
We use multi-wavelength data to associate and characterize the 282 4FGL sources within the 2500 deg 2 SPT-SZ survey field ( Figure 1). We note that the SPT-SZ field is at high Galactic latitude and thus it is assumed that the vast majority of the sources are of extra-galactic origin. The data sets used in our analysis are summarized below.

Gamma rays
We rely on the γ-ray data from the F ermi 4FGL catalog (Abdollahi et al. 2020). We use the source coordinates and 95% uncertainty ellipse for source crossmatching and assume the beam to be Gaussian. The median effective radius for the sources within the 2500 deg 2 SPT-SZ survey field is 3.2 arcmin. The 0.1-100 GeV energy flux and its uncertainty are used to study the multi-wavelength flux correlation and the spectral properties of the associated sources. The class designation in 4FGL is used to evaluate the multi-wavelength associations.
2.2. X-ray X-ray data are from the ROSAT All-Sky Survey (RASS) Bright Source Catalog (Voges et al. 1999) and the Faint Source Catalog (Voges et al. 2000) at X-ray energies 0.1-2.4 keV. The source coordinates and count rate are both used for 4FGL association. We also use source count rate to analyze potential flux correlation. The statistical signal-to-noise ratio reported in the catalog is used when we plot the uncertainties of count rates. Because of the high Galactic latitude of the SPT-SZ field (|b| > 15), the soft X-ray ROSAT data are not affected by photoelectric absorption from our own Galaxy. In addition, blazars generally don't exhibit significant intrinsic soft X-ray absorption (Perlman et al. 2005). Therefore, the ROSAT X-ray measurement should be a reliable measure of the X-ray brightness of the blazars, although dependent on the absorption along the line of sight and the different components sampled in blazars of different classes (e.g. synchrotron emission in BL Lacs and inverse Compton emission in FSRQs). However, the X-ray brightness of a blazar may depend on the blazar sub-class and which component (synchrotron or inverse Compton scattering) is being sampled by ROSAT.

Infrared
Infrared data are taken from the Wide-field Infrared Survey Explorer (WISE ) AllWISE Source Catalog at 3.4, 4.6, 12, and 22 µm (Wright et al. 2010;Cutri & et al. 2013). The angular resolution of WISE ranges from 6 to 12 from short to long wavelengths. The source coordinates and flux at 22 µm (W4) are both used for 4FGL association. The flux at 22 µm is also used to study potential flux correlation. The 4-band magnitudes are used to perform the WISE color analysis on 4FGL blazars. Both the statistical noise reported in the catalog and and additional 10% uncertainty from calibration in W4 (Wright et al. 2010) are included when we plot the flux uncertainty.

Millimeter
Millimeter-wave point sources are from the 2500 deg 2 SPT-SZ survey (Everett et al. 2020) which has a spatial resolution of 1 .15 at 150 GHz (2 mm) and an absolute astrometric uncertainty of 2 (Vieira et al. 2010). The source coordinates and flux at 150 GHz are both used for 4FGL association. The flux at 150 GHz is also used to study the flux correlation and spectral classification. The flux uncertainty includes the statistical noise reported in the catalog, 1.15% uncertainty for the absolute calibration at 150 GHz and 20% uncertainty from The logic sequence of multi-wavelength cross-matching of 4FGL sources. We initially conduct the preliminary 4FGL crossmatching with SUMSS, SPT, and RASS to evaluate the most efficient wavelength for the γ-ray source identification, where the SPT+SUMSS association turns out to maximize the completeness while minimizing the impurity. The true counterparts are selected based on the statistics of the Poissonian probability from the preliminary results. We then take the advantage of mm-radio (SPT+SUMSS) counterparts to refine the position and enable a reliable cross-match with WISE. mm source variability. The 20% variability is the median value of the variance in the light curves of the brightest 200 blazars in the SPT-SZ field. The SPT data used in this paper were collected between 2008 and 2011. Note that the published SPT-SZ catalog reaches down to the 4.5σ significance level, but in this work we extended our search down below 4.5σ to 1σ using forced photometry directly from the 150 GHz maps using the SUMSS positions as priors. This combination of SUMSS (radio) and SPT (mm) associations will hereafter be referred to as SPT+SUMSS counterparts.
2.5. Radio Radio data are from the Sydney University Molonglo Sky Survey (SUMSS, Mauch et al. 2003) at 843 MHz with 45 angular resolution and a detection threshold of 6 mJy. The data were taken between 1997 and 2003. The source coordinates and integrated radio flux density are both used for 4FGL association. The flux density is also used to study the flux correlation and spectral classification. In addition to the statistical noise reported in the catalog, we also adopt a 20% uncertainty to account for source variability. This uncertainty is what we observe for sources in the mm on roughly year-long time scales, but note that this is most likely an underestimate for the source variability on the time scales we are comparing fluxes for this work.

Spectroscopic Redshifts
For each associated 4FGL source, we gather all multiwavelength data and adopt the most accurate position from radio/mm/infrared/optical counterpart if available. We obtained spectroscopic redshifts from the NASA/IPAC Extragalactic Database (NED) for any sources for which they were available. Using the archival redshift data, we study the luminosity distribution, multi-wavelength flux evolution, and classification distribution of the 4FGL blazars.

METHOD
In this section, we describe the method we use for identifying and associating multi-wavelength counterparts to the F ermi 4FGL sources. There are 282 4FGL sources -Selected samples illustrating various kinds of 4FGL associations with mm and radio sources. 4FGL source name and thumbnail index (see Appendix for all the sources studied in this work) are shown at the top-left of each panel. Each 0.7 • × 0.7 • thumbnail in the grey background is the high-pass filtered SPT 150 GHz image. The blue ellipse at the center show the 4FGL 95% uncertainty position and the dashed blue ellipse represents the 4σ uncertainty. The green diamond marks the position of the SPT point source. The red circle shows the position of SUMSS point sources. The first two examples are of previously associated 4FGL sources, which already had known counterparts. In both cases, these known sources are also uniquely identified by SPT (the first is very bright in the mm waveband, while the second is fainter). The last two panels show previously unassociated 4FGL sources. The third panel shows a new unique association with an SPT source. The last panel shows an unusual case where two SPT sources fall within the 4FGL error circle, either or both of which could be producing the γ-ray emission. The grey dots represent all the point sources within 10σ 4FGL positional uncertainty. The vertical and horizontal dashed lines indicate our adopted r < 4σ 4FGL positional separation cut and p-value < 0.1 cut, respectively. The black dots represent the source with the lowest p-value associated to a 4FGL source. orange diamonds represent faint mm sources with fluxes assigned by forced photometry using positional priors from the SUMSS radio catalog. Note that σ 4FGL refers to the standard deviation of the 4FGL position, which varies source-by-source. Left Panel: The SUMSS radio associations are heavily contaminated due to the excessively high source number density; thus, the region of spurious association extends from the top-right region to roughly 4σ 4FGL in separation and 0.1 in probability of false association. Middle Panel: The SPT mm association shows more distinct separation between the two groups of the associations in p-r space, while for radio association, the contamination of the spurious cross-matching causes the large overlap in the spurious region. With both the high completeness and purity, mm provides the best combined performance in completeness and reliability. Right Panel: The RASS X-ray association shows good purity according to the low probability of false association (most grey dots within 4σ 4FGL have p-value less than 0.1) but the completeness is the lowest compared with the others (see Table 1). Multi-wavelength comparisons of sources within 4σ 4FGL in single parameter space are histogrammed in Figure 5. Note that the sources shown at low p-value and separation < 4σ 4FGL which are shown as grey are where there are multiple counterparts for a given 4FGL source and we have adopted the counterpart with the lowest p-value as the association.
within the 2500 deg 2 SPT-SZ survey field, 71 (25%) of which are previously unassociated with any counterpart. Attempting to cross-match the 4FGL sources to external catalogs by selecting all potential counterparts within a given search radius will often yield multiple potential counterparts, particularly when the external catalog has a high surface number density. There are two quantities associated with every external source that we use for cross-matching. The first is the separation between the potential counterpart and the 4FGL position in units of positional uncertainty (σ 4FGL ). The second is the probability of false association for a given external source, which is a function of separation, flux density, and catalog density.
In Figure 2, we show the logical flow diagram we use to identify the multi-wavelength counterparts, which is described in more detail later in this section. Figure 3 shows the 4FGL positional uncertainty ellipses on top of representative thumbnail images from the SPT-SZ survey and also the positions of radio sources from the SUMSS These histograms indicate that SPT+SUMSS is the most efficient identifier of γ-ray counterparts with high completeness and reliability. Left Panel: Histogram of angular separation of 4FGL sources from SUMSS (radio), SPT (mm), SPT+SUMSS (mm), and RASS (X-ray) positions. mm and radio associations are highly complete, while X-ray only has around half of their completeness. Right Panel: Histogram of the probability of false associations of the 4FGL sources with SUMSS, SPT, SPT+SUMSS, and RASS sources. A low probability of false association corresponds to a high certainty of a real counterpart. Thus, the mm and X-ray counterparts are more secure because most of them have the probability of false association less than 0.1, while only around 2/3 radio potential counterparts are within 0.1.
catalog. (see Appendix C for the thumbnails of the entire sample.)

Probability of False Association
To start our cross-matching analysis, we search for external candidates for each source in the 4FGL catalog within its 4σ positional error ellipse. The choice of 4σ was a qualitative choice, designed to be inclusive of the relatively large positional uncertainties of both F ermi and the external catalogs.
In order to robustly identify the most probable counterpart when there are multiple candidates, and to statistically evaluate the odds of spurious associations, we calculate a simple Poissonian probability (Browne & Cohen 1978;Downes et al. 1986;Biggs et al. 2011), which gives the probability of a matched source being randomly associated within a given area. The Poissonian probability takes into account the angular separation, flux densities, and the number density of potential counterparts. The expected number of random associations is: where a and b are semi-major and semi-minor axes, respectively, and n s (> S) is the surface number density of sources brighter than the candidate counterpart. The Poissonian probability (hereafter "p-value"), is thus defined as: The p-value constructed in this way is the probability that an association is false, due to spurious coincidence. Thus, if p 1, the association is highly likely to be real. We note that there are small corrections that can be made to this probability (see, e.g. Downes et al. 1986;Biggs et al. 2011), which we neglect in this work as the density of sources is relatively low and our sources of interest are rare. In this work, we accept all sources as potential counterparts with p-value < 0.1 within the 4σ association area. To visualize our cross-matching selection procedure, Figure 4 shows all the sources within 10σ of the 4FGL position with the 4σ 4FGL and p-value < 0.1 cuts indicated by dashed lines. We plot the p-value versus the angular separation in units of positional error (σ 4FGL ) for each of the external catalogs. In some cases there are multiple counterparts that meet this criteria. When this is the case, we adopt the counterpart with the lowest p-value as the most probable counterpart and show that source in Figure 4 with a black point.

Completeness and Purity of Individual Catalogs
In Figure 5, we quantify the behavior of the multiwavelength associations in the space of angular separation and the p-value by histograming potential multiwavelength counterparts of 4FGL point sources within their 4σ-beam before cutting on p-value.
In the left panel of Figure 5, we show the histogram of angular separation in units of positional error for each of the external catalogs within the adopted cutoff of 4σ 4FGL . What is immediately apparent is that the curves are mostly flat above 1σ 4FGL for SUMSS and RASS, but peaks around 1σ 4FGL for SPT (or SPT+SUMSS), as the separations increase. In the case of SUMSS, the sources Note.
-Three wavelengths -X-ray (RASS), mm (SPT and SPT+SUMSS), radio (SUMSS) -are used to study the performance of the 4FGL association within the 2500 deg 2 sky covered by the SPT-SZ survey. The second column refers to the surface number density of each survey. The third to fifth column represent the total number of potential counterparts for all 4FGL sources, previously identified 4FGL sources, and previously unassociated 4FGL sources, respectively. The last two columns are completeness and purity for the multi-wavelength association, where the error of completeness is Poissonian and error of purity is evaluated from bootstrapping. density is high (∼ 30 sources per square degree) and so there are always potential counterparts to match with, and there are more the further out you include, and so the problem becomes in deciding which is the true counterpart. For RASS, there are simply not many counterparts to match with. SPT (or SPT+SUMSS), however, often has a counterpart to 4FGL, and the density of background sources is lower, so the association histogram peaks and then decreases at larger separations.
In the right panel of Figure 5, we plot the histogram of p-value for all sources within the 4σ 4FGL radius. The distribution of SUMSS counterparts is bimodal, with the counterparts split into two groups near 0 and 1 respectively, which indicates that nearly 1/3 of SUMSS sources inside 4σ 4FGL beams are likely to be spurious associations. The SPT (or SPT+SUMSS) and RASS catalogs, however, with much lower densities, have just one peak at low p-value, indicating a high certainty of association.
We adopt the critical p-value (p cri ) to best separate the real and false associations. We chose 10% as an acceptable false association rate and indicate this cut by the horizontal dashed lines in Figure 4.
For the purposes of defining the efficiency of an external catalog associating with 4FGL sources, we define completeness as the probability of a catalog providing at least one viable counterpart (separation < 4σ 4FGL and p-value > 0.1) to a source in the 4FGL catalog, and the purity as the probability that the association is false (separation < 4σ 4FGL ). To calculate the purity of each catalog we use the equation: where N and r refer to number of samples and angular separation, respectively. The X-ray (RASS) associations have a high purity ( Figure 5 right panel) but only 54% 4FGL sources have RASS sources inside their 4σ uncertainty regions, and is thus largely incomplete.
Matching the 4FGL sources to the SUMSS catalog we find that 95% of the 4FGL sources have at least one SUMSS source falling inside their 4σ uncertainty region. However, as Figure 4 demonstrates, many SUMSS sources which are within the 4σ 4FGL uncertainty regions also have a high probability of being spurious associations, i.e. they lie near p = 1 and r < 4σ 4FGL in p-r space. Figure 5 also demonstrates that ∼ 1/3 of SUMSS associations have high probability of false association (p-value > 0.1) and there are often multiple possible radio counterparts within the 4σ uncertainty region. We find that ∼55% of 4FGL sources have more than one potential SUMSS counterpart within the 4FGL positional uncertainty. This demonstrates that while the radio catalog has a high completeness, it also has a low purity (i.e. low confidence of a true association). When we cut the sources within 4σ 4FGL , we are left with 232 associations and a completeness of 82%.
Matching the 4FGL sources to the SPT-SZ catalog we find that 204 (72%) of the 4FGL sources have at least one SPT source falling inside their 4σ uncertainty regions and p-value < 0.1. Each SPT source was also crossmatched with SUMSS, within 1 arcmin of the SPT position at 150 GHz and with p-value less than 0.1. Nearly all (99%) of the SPT sources that are located inside 4σ uncertainty regions of 4FGL have SUMSS counterparts. Only five 4FGL sources (#80, #127, #177, #211, and #225 in Appendix A. For details see thumbnails in Appendix C) have SPT counterparts but no SUMSS counterparts. When we include the faint SPT fluxes at > 1σ derived from the SUMSS positions (SPT+SUMSS) we find an additional 35 have mm counterparts within the 4σ 4FGL beams and p-value < 0.1, bringing the completeness to 85%. The mm associations also have high completeness, but because of the far lower surface density of SPT sources, the associations retain a high purity (see Figure 5, right panel). Among the 71 previously unassociated 4FGL sources in this region, 40 (56%) have mm counterparts in the SPT-SZ data.
In Table 1, we show the completeness and purity of each of the catalogs we associate to the 4FGL sources after the p-value cut. We estimate the purity for each association as discussed earlier. The completeness is defined as the fraction of 4FGL sources with a probable counterpart. In the results of multi-wavelength associations, the mm (including the faint sources from the forced photometry using SUMSS priors) and radio bands have high completeness of association (85 % and 82 %, respectively), while the X-ray band is less than half complete (48 %).
We thus conclude that the SPT+SUMSS is the most efficient means to identify 4FGL sources, as it maximizes both completeness and purity. Almost all of the SPT identified sources (99%) have corresponding SUMSS counterparts which appear to be flat-spectrum at mm and radio wavelengths (Figure 7).
The results of multi-wavelength associations demonstrate the promise of mm association to identify γ-ray sources. Thus, the combination of SUMSS (radio) and SPT (mm) associations (SPT+SUMSS) can best characterize the 4FGL sources because of the high certainty and high completeness of γ-ray association with joint mm and radio counterparts. For the previously identified 4FGL sources without SPT+SUMSS counterparts, we find that three of them are pulsars, one is a faint BL Lac, and the rest are undetermined blazar candidates with X-ray emission.

The Construction of Our Combined 4FGL
Multi-Wavelength Catalog To summarize the previous sections, we adopt the following criteria to select multi-wavelength counterparts of 4FGL sources for the analyses in this work (see Figure 2): We select SUMSS/SPT(+SUMSS)/RASS counterparts within 4σ 4FGL and p-value < 0.1 as the 4FGL associated counterparts.
The 4FGL-SUMSS cross-matching, while the most complete, is often degenerate with multiple possible associations. Recall that to help break the degeneracy of multiple possible SUMSS counterparts where there is no SPT source at ≥ 4.5σ, we perform forced photometry at the SUMSS position in the SPT 150 GHz map. In this way, we are able to dig deeper into the SPT map and associate the most probable radio counterpart based on the mm flux. When there are multiple sources which meet the criteria of 4σ-beams and p-value < 0.1, we adopt the source with the lowest p-value in the mm. Some examples are shown in Figure 3 to illustrate how the 4FGL sources are associated by SPT and SUMSS. Of the 239 4FGL-SPT matches, 204 sources are detected at > 4.5σ in Everett et al. (2020), while the remaining 35 have had mm-fluxes assigned using forced photometry using the SUMSS position as a prior. For clarity, these sources are highlighted in Figure 7 and Figure 8 labelled as FSPT for "faint". Now, with this catalog of multi-wavelength associations in hand, we can associate the 4FGL sources to the WISE catalog to characterize their infrared emission. Naive cross-matching 4FGL with WISE is difficult because the high source density (∼10 4 deg −2 ) results in high contamination from spurious associations. Each 4FGL source on average has ∼200 potential WISE counterparts within the 95% (2σ) uncertainty ellipse. However, SPT-selected blazars with a SUMSS counterpart can provide positional accuracy better than 15 (Vieira et al. 2010) when the signal-to-noise ratio is greater than 5 (Ivison et al. 2007). Using the 4FGL-SPT-SUMSS associations, we performed the 4FGL association with WISE based on the refined position, and selected the most likely candidates according to their lowest p-value. When cross-matching with WISE, the p-value is especially useful to eliminate spurious false associations. We adopt WISE associations using p-value < 0.03 and separation within 9 arcsec of the combined SPT+SUMSS position. The characterization of the infrared emission of the 4FGL sources is discussed in Section 4.3.
Nearly all the SPT counterparts (99%) of 4FGL sources also have a SUMSS counterpart allowing us to characterize each 4FGL source with both mm and radio fluxes. In the Figure 7, mm sources with radio counterparts (grey dots) are roughly separated into two classes -steep-spectrum AGN (yellow oval) and flat-spectrum AGN (blue oval). The majority of γ-ray sources (blue and green dots) have comparable mm and radio emission, which indicates that they are either BL Lacertae objects (BL Lacs) or flat-spectrum radio quasars (FSRQs).

RESULTS & DISCUSSION
Robust multi-wavelength counterpart identification to the sources in the 4FGL γ-ray catalog is crucial for understanding the γ-ray source population and the diffuse γ-ray background. In this section we study and discuss the multi-wavelength properties of F ermi γ-ray blazars, including the multi-wavelength flux correlation, multiwavelength color analysis, redshift dependency, and implications for future surveys.
We adopt two approaches to study the multiwavelength properties of the 4FGL sources. The first approach directly adopts the multi-wavelength associations derived from SPT+SUMSS and described in Section 3. For the second approach, we divide the 4FGL sources into two groups -previously associated 4FGL sources (ID-4FGL) and previously unassociated 4FGL sources with new mm identifications (UID-4FGL-SPT). For ID-4FGL sources (which have well-measured ∼arcsecond associated astrometry) we apply a simple angular-separation-based cross-matching with the external multi-wavelength catalogs. For UID-4FGL sources, i.e. the ones without a previous counterpart association, we use the counterpart associations found in this work using SPT+SUMSS. Both methods -either using the provided cross-matching from the 4FGL catalog or our own independent cross-matching using SPT+SUMSSproduce qualitatively similar results. For the rest of the figures in this paper, for the ID-4FGL sources we use the supplied position for multi-wavelength associations, and for the UID-4FGL sources we use the UID-4FGL-SPT position derived in this work.
The spectroscopic redshift data were acquired from the NED based on the refined source positions from the multi-wavelength associations. We found 75 sources (out of 239) with redshift measurements.

The γ-ray Flux Correlation
With the cross-matched multi-wavelength catalog in hand, we can study the flux correlation across bands for the 4FGL catalog. For each band, we calculate the Spearman's rank correlation coefficient (r SC ) of the flux correlation and the results are shown in Figure 6. The UID-4FGL-SPT sources (black dots) are typically faint sources at all wavelengths studied here. The inclusion of these faint sources increases the completeness the 4FGL associations.
Starting in the radio, we see that there is a slight correlation between the radio and γ-ray flux, with r SC ∼ 0.3.
The mm flux is well-correlated (r SC ∼ 0.5) to the γ-ray flux and can be parameterized by the following relation: Most of these sources are flat spectrum AGN, but a few sources have stronger radio emission (see 4FGL sources with S 150GHz > 100 mJy outside the blue oval in Figure 7). This may be partially related to source variability and may contribute to the radio flux correlating less well than the mm flux. -Flux correlations between 4FGL and other surveys (Xray: blue; Infrared: light blue; mm: green; radio: red), where black dots indicate previously unassociated 4FGL sources with mm counterparts. The Spearman's rank correlation coefficient (r SC ) of the flux correlation is shown at the top-right of each plot. Note that r SC,ID refers to previously associated 4FGL sources and r SC,UID+ID counts both previously associated and extra mmidentified 4FGL sources. Panel (a): As most of the 4FGL-SPT blazars are flat spectrum AGN, the r SC in the radio-band should be similar that of the mm-band. However, r SC in the radio-band is significantly lower than the mm-band because some 4FGL-SPT blazars have excess radio emission compared with most flat spectrum sources (see Figure 7). Panel (b): the mm and γ-ray bands are highly correlated, which is consistent with previous studies of multi-wavelength associations of 4FGL sources. The dashed line represents least squares fitting between the fluxes in the two bands. Panel (c): The infrared-band has comparable correlation to the γ-ray band as the mm-band. However, the association between 4FGL and WISE has been refined by first matching to SPT because normal cross-matching would be heavily contaminated by the spurious associations due to the high source number density of WISE (∼ 10 4 deg −2 ). Panel (d): Unlike any other wavelengths, the X-ray fluxes are uncorrelated with γ-ray fluxes. This may be because the γ-ray traces the jet while the X-ray traces the coronal region of the AGN (Ghisellini et al. 2017). The grey dots represent all the mm sources with radio counterparts. The blue/green dots represent those with previous known/unknown γ-ray detection. The orange diamonds represent faint mm sources with fluxes assigned by forced photometry using positional priors from the SUMSS radio catalog. Note that the steep-spectrum AGN roughly lie within the yellow oval and run up the vertical axis clustered around S 150GHz =10 mJy, while the flat spectrum AGN lie within the blue oval which is nearly diagonal with slope=1. 239 4FGL sources are selected if they have both SPT and SUMSS counterparts. These sources are shown in green and blue. Note the total number of both blue and green dots in the plot is 265. This is because some 4FGL sources have multiple mm-radio counterparts as illustrated in Figure 3. Similar to Figure 1, the blue points are previously associated 4FGL sources, while the green points are previously unassociated 4FGL sources that have SPT counterparts. The majority of 4FGL sources are flat spectrum AGN.
The infrared flux is statistically as correlated as the mm (r SC ∼ 0.5). However, given the high density of the infrared WISE catalog, we first needed to match the 4FGL to a SPT+SUMSS source, and then to WISE in order to refine the positional uncertainty. Thus, this correlation may be biased towards the mm-detected sources.
The X-ray flux is largely uncorrelated to the γ-ray flux (r SC ∼ 0). Naively, the lack of correlation between the γ-ray and X-ray is a surprise, given their proximity along the electromagnetic spectrum, particularly relative to the radio. This is likely due to the fact, that while gamma rays in the jet are always produce via inverse Compton scattering (of synchrotron photons in BL Lac, or external photons in FSRQs), the X-rays are generated via synchrotron in BL Lacs and via IC scattering in FSRQs.
Our study demonstrates that the mm band is the most efficient band to associate γ-ray blazars to multiwavelength counterparts. As shown in the Method Section 3, RASS's X-ray sensitivity is insufficient to detect all the 4FGL counterparts and the X-ray flux of 4FGL sources is observed to be uncorrelated to their γ-ray fluxes ( Figure 6). Thus, the RASS's X-ray catalog is highly incomplete in terms of associating 4FGL sources orange diamonds). The previously identified 4FGL sources without mm detections (ID-4FGL-NoSPT) and the remaining unidentified 4FGL sources (UID-4FGL) are indicated by blue and orange arrows, respectively, as 3σ upper limits from SPT-SZ. Note that three bright ID-4FGL-NoSPT sources (S 0.1−100GeV > 10 −11 erg cm −2 s −1 ) in blue arrows are previously identified pulsars and therefore fails to be identified by SPT. Within associated 4FGL sources, FSRQs are more likely to be brighter than BL Lacs in mm wavelength. Most unassociated 4FGL sources are faint in both mm and γ-ray bands, which indicates that source identification might be limited by the sensitivity of the current generation of surveys. Lower Panel: The mm-radio spectral index (α 843MHz 150GHz ) as the function of γ-ray flux. Sources are labeled the same way as the left panel. Cyan dots denote the median and standard deviations of the spectral index in each logarithmic bin of γ-ray flux. As shown in the plot, sources with brighter γ-ray fluxes are more likely to have flatter spectra except a few sources with γ-ray flux around 10 −11 erg cm −2 s −1 . These outliers have excess radio emission compared with normal flat-spectrum AGN. Among these sources, some might have extra radio emission from nearby radio lobes; some are just ambiguous cross-matching that blends several radio counterparts; some have valid multi-band counterparts and need further investigation.
to multi-wavelength counterparts. The SUMSS radio associations have a high completeness, but often have multiple counterparts and are thus confused. WISE, due to its high source density, produces a large number of spurious associations, and thus needs an accurate prior on the position to enable accurate source association. The mm catalog with arcminute resolution, such as the SPT catalog, provides 4FGL associations with both high completeness and purity (see Table 1 and Figure 5).

Previously Unidentified and Faint γ-ray Population
As shown in Figure 6, γ-ray emission is correlated to the mm emission. Thus, the brighter 4FGL sources have a higher completeness of mm associations. The previously associated 4FGL sources (ID) typically have brighter multi-wavelength fluxes (see non-black sources in Figure 6), where 94% of them are also SPT-identified. The previously unassociated 4FGL sources (UID) tend to be fainter (S 0.1−100GeV <10 −11 erg cm −2 s −1 ) at longer wavelengths (see black sources in Figure 6), where only On the background color map, each dot marks the SPT-selected source. These sources are also labeled by their identifications (previously associated 4FGL source: blue dot; previously unassociated 4FGL source with mm counterpart: green triangle; BL Lac: black diamond; FSRQ: red square). Red dashed lines outline the WISE Gamma-ray Strip (WGS, Massaro et al. 2012). Most of the mm associations lie inside the WGS, which indicates the consistency of mm correlation and WGS parameterization. Besides, the mm cross-matching provides more intrinsic associations and pinpoints the valid outliers outside WGS in WISE color space.
56% of them are SPT-identified. (See Table 1 for a summary of these associations). In addition, as can be seen in Figure 8 (left), the remaining 4FGL sources without SPT-identification are mostly faint γ-ray sources (S 0.1−100GeV 5 × 10 −12 erg cm −2 s −1 ), which might explain why they have remained unassociated thus far. Therefore, the identification of the remaining unassociated 4FGL blazars will be further completed by either deeper catalogs from upcoming surveys in mm wavelengths (e.g. SPT-3G, Simons Observatory, CMB-S4), or by dedicated pointed observations with e.g. ALMA, SMA, or NOEMA.
In Figure 7 we plot the radio versus mm flux for all sources in the SPT catalog with a radio counterpart and highlight the sources with γ-ray counterparts. Some flat spectrum AGN with strong mm emission (S 150GHz > 100 mJy) are still undetected in γ-ray. Roughly half of them are γ-ray-quiet blazars and can also be found in CGRaBS (see Healey et al. 2008;Paliya et al. 2017) or ROMA-BZCAT (Massaro et al. 2015). The rest of them are also blazar-like, but lack γ-ray detection. The γray-quiet blazars are often associated with small Doppler factors and high disk dominance (Paliya et al. 2017). As noted by several authors (Dermer 1995;Paliya et al. 2017;Ghisellini et al. 2017), there are two important selection effects which could make a luminous but highredshift blazar not detected by F ermi: (1) Luminous sources usually have the high-energy peaks of their SEDs Most wavelengths (radio, mm, infrared, and γ-ray) show no apparent correlations between flux density and redshift, while in X-rays the upper bound of count rate declines for the fainter sources. Therefore, at most wavelengths, flat-spectrum AGN are equally likely to be detected given the flux above the threshold of flux detection. In the X-rays, the extra flux dependence constraining the upper bound is due to the relative shallow depth of RASS. The bright X-ray sources at high redshift are either less likely to be detected by RASS, or the origin of X-ray emission can be different from other wavelengths.
shifted to lower energies, even as measured in the source frame, which is in addition to the effect of redshift itself.
(2) These sources may be highly beamed at gamma rays and slightly misaligned making them currently undetected by Fermi-LAT. Unfortunately we do not have the data necessary to resolve this issue at this time, but, for instance, a stacking analysis could be done in the Grey dots represent all the SPT sources with redshift measurements. Those with 4FGL counterparts (4FGL-SPT) are covered by blue dots. We also label the blazar type for each 4FGL-SPT source if it is classified in the catalog. Given the same redshift, mm-identified 4FGL sources (flat spectrum AGN) are brightest, where the specific luminosity is above the threshold of mm detection by orders of magnitude. This is resulting from the mixing of two types of sources (flat-spectrum and steep-spectrum AGN) and the selection bias of redshift data. Details are discussed in the discussion section (4). We can roughly distinguish the population of FSRQs and BL Lacs by their mm specific luminosity. Right Panel: Histogram of redshift distribution. The label colors are consistent with the left plot and redshifts are binned logarithmically. This histogram shows that the population of FSRQs and BL Lacs can be distinguished by their redshift measurements as well. FSRQs are mostly observed with higher redshift than BL Lacs. This distinction should come from the detection bias of optical spectrum measurement. In general, FSRQs have strong emission lines in optical broadband, while the broadband spectra of BL Lacs have weak emission lines or even featureless.

4.3.
Gamma-ray Blazars in Radio, mm, and Mid-IR By using multi-wavelength profiles and labeling the source type of associated 4FGL sources, we studied how the spectral characteristics influence the patterns in multiple parameter spaces.
First, we looked into the SPT-identified 4FGL sources within the SPT-SZ field in the radio, mm, and γ-ray band fluxes. FSRQs and BL Lacs have already been classified in the 4FGL catalog based on their optical spectra (Abdollahi et al. 2020). As shown in the left panel of Figure 8, FSRQs are generally brighter than BL Lacs in mm wavelengths, but indistinguishable in γ-ray. The right panel demonstrates how the spectral indices vary with γ-ray flux, where the spectral index is defined as A flat spectrum from mm to radio corresponds to α 843MHz 150GHz = 0. As shown in the right panel of Figure 8, the radio spectral index of bright γ-ray sources tend to be more flat, presumably because the viewing angle is more aligned with the center of the jet thus the observed emission is dominated by the central engine of the jet. As γ-ray flux decreases below ∼ 5×10 −12 erg cm −2 s −1 , flatspectrum and steep-spectrum AGN are heavily mixed together (S 843MHz and S 150GHz 100 mJy in Figure 7), so the spectral index exhibits a higher scatter at lower γray fluxes. Most of the sources (< 10 −11 erg cm −2 s −1 ) have a steeper (i.e. less flat) radio spectral index because they are less jet-dominated. The radio emission from jets likely originates from an optically thick regime. Therefore, there could be additional radio emission if the surface is large. However, the jet is optically thin in the mm, where mm and γ-ray emission are likely to be produced co-spatially in the compact emission region in the jets (Meyer et al. 2019). This is consistent with the higher flux correlation between mm and γ-ray than with the radio (see Figure 6).
Next, we studied the 4FGL sources with joint mm and infrared counterparts in WISE color space. The WISE two-color diagram has clear patterns for the classes of WISE objects (see also Wright et al. 2010). As shown in Figure 9, the background color contours represent the source number density of WISE objects, with density increasing from blue to red. Each population of WISE objects is labeled on the diagram. For the SPT sources with WISE counterparts, some of them have infrared colors indicative of being AGN-dominated, while others have a substantial component of starlight from an earlytype host galaxy. Many other SPT sources which are not γ-ray emitters have very blue [3.4]- [4.6] colors which is indicative of redshifts reaching up to around 1. The red dashed lines outline the region of WISE Gammaray Strip (WGS, Massaro et al. 2012) which indicates the location of the known blazars from ROMA-BZCAT in WISE color space. The majority of SPT-identified blazars with γ-ray emission are located within both QSOs population and WGS. We find general consistency between our association results with WISE source classification and WGS parameterization. Moreover, since the mm is efficient at associating γ-ray-loud flat spectrum AGN, it also provides a method to validate outliers from the WGS region and also pinpoint other possible populations that contain γ-ray emission along with infrared emission.

The Redshift Distribution of Gamma-ray Blazars
With previous multi-wavelength associations, we obtained the redshifts for 31% (75 out of 239 sources) of associated 4FGL sources by cross-matching NED database. An important caveat to the discussion below is that this sample is highly spectroscopically incomplete, and could be highly susceptible to selection biases. For instance, the sources with spectroscopic redshifts are presumably biased towards being optically bright and having strong emission lines. A future spectroscopic survey, or a dedicated targeted spectroscopic campaign, would be needed before any strong conclusions are drawn from this particular discussion.
In Figure 10 we investigate the redshift dependence on the multi-wavelength detection of 4FGL sources. With the exception of the X-ray band, the majority of multiwavelength fluxes of 4FGL sources have no clear redshift dependence. In the left panel of Figure 11, SPTidentified 4FGL sources (blue) are demonstrably more luminous in the mm-band than the sources not detected in 4FGL. In addition to the 4FGL-SPT sources, we also label the 4FGL classification of FSRQ (red) and BL Lac (black) and histogram the redshift distribution based on their types. In the right panel of Figure 11, BL Lacs are mostly observed with lower redshift than FSRQs. This distinction may be due to an observational bias as BL Lacs generally have weak optical emission lines.

CONCLUSION
1. We have shown that the mm flux detected by SPT from flat-spectrum AGN is the most effective method currently available for identifying extragalactic 4FGL γray sources and predicting the possible γ-ray emission. The SPT detection, combined with an accurate source position, finds 4FGL sources with the highest completeness rate and, at the same time, the lowest contamination rate.
2. The effectiveness of the SPT+SUMSS selection has been demonstrated first by confirming the association of 94% (199 out of 211 sources) of already-known 4FGL sources across the 2500 deg 2 SPT-SZ survey field. It is then applied to identify 40 new sources for which 4FGL did not previously have a counterpart at lower energies (out of 71 previously unidentified sources within the SPT-SZ survey field). Deeper and wider mm surveys will soon be available which will greatly complete these 4FGL associations.
3. We have used multi-wavelength data to explain why SPT is more effective to find associations for γ-ray sources. Our interpretation is that the mm flux is closely correlated with the γ-ray flux (Meyer et al. 2019) because both have a common origin in the jet.
4. SPT has shown its extraordinary ability to track jet-dominant AGN and find blazar-like objects. Therefore, SPT can also complete the sampling of γ-rayquiet blazars. SPT has detected 60 bright mm emitters (S 150GHz > 100 mJy) which do not currently have any γ-ray detection. Roughly half of them are γ-rayquiet blazars and can be found in CGRaBS (Paliya et al. 2017) or ROMA-BZCAT (Massaro et al. 2015). The rest of them are also blazar-like sources and lack multiwavelength detection.
5. SPT-3G (Benson et al. 2014) will survey 1500 deg 2 of southern sky ×10 deeper and with polarization sensitivity. The DOE has recently begun planning a nextgeneration experiment (CMB- Stage 4, Abazajian et al. 2016) to cover the entire extragalactic sky. Thus, this initial study enables us to prepare forecasts for the nextgeneration of SPT surveys and CMB experiments which will extend this technique to greater sensitivities and across the entire sky. With much more powerful data sets from SPT-3G and CMB-S4, we should be able to complete the association of the remaining unassociated 4FGL sources and also study the light curves of 4FGL blazars in both mm and γ-ray. The variability analysis will directly answer if mm and γ-ray emission from 4FGL blazars are intrinsically correlated, which is indicative of the radiation process inside relativistic jets.
While deeper future surveys in the X-ray (e.g. eROSITA), the radio (e.g. VLASS, MeerKAT, ASKAP, SKA), and the mm regime (SPT-3G, CMB-S4) will help advance this work and improve statistics, we believe these general characteristics of the multi-wavelength associations to γ-ray catalogs will remain largely unchanged.
We list all 282 4FGL sources within 2500 deg 2 SPT-SZ survey field. For each 4FGL source, all the SPT counterparts are listed on the right. The index is based on the sequence of the energy flux (100 MeV-100 GeV) from high to low. In this table, 4FGL sources with '*' represents that the source has no associated counterpart in the original 4FGL catalog. Most of SPT sources have detections greater than 4.5σ except those SPT sources whose name starts with 'F' (for "faint"), which represent the detection less than 4.5σ but greater than 1σ. Recall that the uncertainty of mm flux includes mm variability, statistical and systematic noises, while the uncertainty of 4FGL energy flux is just the statistical noise. Source separations are listed in units of arcmin (dr 1 ) and σ of the 4FGL positional uncertainty (dr 2 ).