The Remarkable Predictive Power of Infrared Data of Blazars

Blazars are the brightest and most abundant persistent sources in the extragalactic γ-ray sky. Due to their significance, they are often observed across various energy bands, where the data of which can be used to explore potential correlations between emission at different energies, yielding valuable insights into the emission processes of their powerful jets. In this study we utilized IR data at 3.4 and 4.6 μm from the Near-Earth Object Wide-field Infrared Survey Explorer Reactivation Mission, spanning 8 yr of observations, X-ray data from the Neil Gehrels Swift Observatory collected throughout the satellite’s lifetime, and 12 years of γ-ray measurements from the Fermi Large Area Telescope’s all-sky survey. Our analysis reveals that the IR spectral slope reliably predicts the peak frequency and maximum intensity of the synchrotron component of blazar spectral energy distributions, provided it is uncontaminated by radiation unrelated to the jet. A notable correlation between the IR and γ-ray fluxes was observed, with the BL Lacertae subclass of blazars displaying a strong correlation coefficient of r = 0.80. IR band variability is more pronounced in flat spectrum radio quasars than in BL Lacertae objects, with mean fractional variability values of 0.65 and 0.35, respectively. We also observed that the synchrotron peak intensity of intermediate-high-energy-peaked objects can forecast their detectability at very high γ-ray energies. We used this predicting power to identify objects in current catalogs that could meet the detection threshold of the Cerenkov Telescope Array extragalactic survey, which should encompass approximately 180 blazars.


INTRODUCTION
Blazars are a unique category of active galactic nuclei (AGNs) with a distinctive observational perspective due to the alignment of their jets with the observer's line of sight Urry & Padovani (1995).This alignment results in a Doppler boost of the emission from the relativistic jet, which often outshines the emissions from other components (Padovani et al. 2017).This unique characteristic makes blazars a special laboratory for investigating the formation, propagation, and processes within relativistic jets.
One of the distinctive features of blazars is that their emission spans from radio to high energy (HE; > 100 MeV) or very high energy (VHE; > 100 GeV) γ-ray bands (Padovani et al. 2017), a very broad multiwavelength emission that is known to vary on short timescales and often with high amplitude (see e.g., Aleksić et al. 2014;Ackermann et al. 2016;Shukla et al. 2018;Hood et al. 2023).The spectral energy distribution (SED) of blazars is typically characterized by a double bump structure.The first peak is observed in the infra-red, optical or X-ray bands, while the second peak is typically in the γ-ray band.The first component of the SED is believed to be due to the synchrotron emission of electrons accelerated within the jet.The origin of the second peak can be attributed to the interaction of either electrons (in the leptonic scenario), protons (in the hadronic scenario), or electrons and protons (in the lepto-hadronic scenario).In the leptonic case, the second peak is due to inverse Compton scattering of electrons, which can be either internal (synchrotron-self Compton model, SSC; Ghisellini et al. 1985;Maraschi et al. 1992;Bloom & Marscher 1996) or external, as in the external inverse Compton model (EIC, e.g., Sikora et al. 1994;Dermer et al. 1992;Dermer & Schlickeiser 1994;Sikora et al. 1994;B lażejowski et al. 2000).In the case of hadronic or lepto-hadronic models, the accelerated protons can significantly contribute to the formation of the second peak.This contribution can be either by direct synchrotron emission (Mücke & Protheroe 2001) or by the secondaries generated in photo-pion and photo-pair interactions (Mannheim 1993;Mannheim & Biermann 1989;Mücke & Protheroe 2001;Mücke et al. 2003;Böttcher et al. 2013;Petropoulou & Mastichiadis 2015;Gasparyan et al. 2022).The hadronic or lepto-hadronic models are particularly interesting as they can account for the neutrinos observed from the direction of blazars; as an example of TXS 0506+056 (IceCube Collaboration et al. 2018a,b;Padovani et al. 2018) and PKS 0735+178 (Sahakyan et al. 2023).
Blazar populations have been traditionally subdivided into two large classes based on their optical emission lines: Flat Spectrum Radio Quasars (FSRQs), which exhibit strong, quasar-like line emissions, and BL Lacertae objects (BL Lacs), which show either no or weak emission lines.The observations of a peculiar behavior in the γ-ray blazar 4FGL J1544.3-0649showed the possible existence of transient-like blazars (Sahakyan & Giommi 2021).This source remained undetected in the X-ray and γ-ray bands until May 2017, then it transformed into a highly luminous source for several months and was detected by the Fermi Large Area Telescope (Fermi-LAT) and the All-sky X-ray Image (MAXI) sky monitor.This observation implies the potential existence of a yet-to-be-discovered population of blazars that may sporadically exhibit very large flaring activity.
Blazars are further classified based on their broad-band emission, depending on where the peak of their synchrotron component (ν s p ) is located: they are called low-energy peaked (LSP/LBL) when ν s p < 10 14 Hz, intermediate-energy peaked (ISP/IBL) when 10 14 Hz < ν s p < 10 15 Hz, and high-energy peaked (HSP/HBL) when ν s p > 10 15 Hz (Padovani & Giommi 1995;Abdo et al. 2010).More recently Giommi & Padovani (2021) have shown that IBLs and HBLs, collectively referred to as intermediate-high-energy-peaked objects (IHBLs), share common properties that differ from those of LBLs.As a result, blazar classification can be simplified to just two categories: LBLs, which have a peak ν s p < 10 13.5 Hz, and IHBLs, which have ν s p > 10 13.5 .Throughout this paper, we will use the LBL and IHBL classification.
The two classifications largely, but not completely, overlap since FSRQs are nearly always LBLs and BL Lacs are both LBLs and, more frequently, HBLs.This could be the result of selection effects induced by the strong non-thermal radiation from the jet that occasionally outshines the broad line emission, or to deeper physical differences between the two classes, like, e.g.strong differences in the accretion disk properties with FSRQs having more radiatively efficient disks than BL Lacs or, where the higher energy emission in IHBLs results from hadronic processing associated to neutrino emission, like in the hypothesis put forward by Giommi & Padovani (2021).
Due to their prominent emission in all bands of the electromagnetic spectrum, blazars are often targets of multiwavelength and multimessenger observations.These observational campaigns not only generate large amount of data, which can be used to build different parts of their multiwavelength SEDs, but also enable correlation studies between the emission in different bands.The correlation between the emission in the HE γ-ray band with those at lower energies, e.g.infrared (IR) or optical, is one of the keys to understand the emission processes in these objects.Previous studies have already compared the γ-ray and infrared properties of blazars using a smaller number of sources and data from WISE.Specifically, Massaro et al. (2011) demonstrated that blazars cover a distinct region in the 3.4-4.6-12µm color-color diagram, referred to as the WISE blazar Strip.Additionally, Massaro & D'Abrusco (2016) found a tight correlation between the mid-IR colors and the γ-ray index.
In this work, we present a comprehensive study on the properties of blazars as observed in the γ-ray and IR bands.Since 2008, a large number of blazars are continuously monitored in the HE γ-ray band by Fermi-LAT, offering deep insight into their γ-ray emission properties.Fermi-LAT is sensitive to γ-rays with energies from 100 MeV to over 300 GeV, and it systematically scans the entire sky every three hours (Atwood et al. 2009).The Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010) was launched in 2009, and during the nominal mission conducted a comprehensive survey of the IR sky, specifically in the 3.4, 4.6, 12, and 22 µm bands.After the conclusion of its primary mission in 2011, the WISE spacecraft was placed in hibernation mode until 2013 when it was reactivated under the name NEOWISE.Its primary goal was to map the entire sky in the 3.4 µm and 4.6 µm bands every six months (Mainzer et al. 2014).Both WISE and NEOWISE detected a large number of blazars thus yielding a very rich data archive.
Motivated by the substantial number of blazars observed by Fermi-LAT that have also been surveyed in the IR band by WISE/NEOWISE, we performed in an in-depth comparison of blazar emission features in the IR and γ-ray bands.The organization of this paper is as follows: Section 2 the source sample selection is described.The estimation of the BLs IHBLs synchrotron peak from IR data is presented in Section 3. Section 4 discusses the relationship between IR and γ-ray emissions.Finally, Section 5 encompasses the discussion and conclusions.

SOURCE SAMPLE
To identify γ-ray-emitting blazars that have been also detected in the IR band we cross-matched the Fermi Large Area Telescope Fourth Source Catalog data release 3 (4FGL-DR3 Abdollahi et al. 2022) with the catalog of NEOWISE detected sources in the 2022 Data Release (Mainzer et al. 2014).The latter provides a single file containing data from eight years of full-sky monitoring1 .The 4FGL-DR3 catalogue, which spans the period from 2008 to 2020, is a comprehensive survey of the γ-ray sky including 6,658 sources, 3,743 of which are blazars of various types.The NEOWISE Data Release products offer flux measured at 3.4 and 4.6 µm, taken from eight years of the survey.During this period, the entire sky was scanned nearly sixteen times.
Out of all the sources observed by Fermi-LAT, we focused our study on the high Galactic latitude blazars (|Gal.latitude| > 10 • ) that were also observed in the NEOWISE surveys, totaling 2,279 objects.In order to check the effects of source confusion, we cross-matched the position of our target source with entries in the Hubble Space Telescope Guide Star Catalog (HSTGSC, Lasker et al. 2008).If another object was found within 9.0 arcsec of the target source, the target source itself was excluded from our analysis.For each of the selected sources, we gathered NEOWISE and other archival multiwavelength data using the VOU-Blazars tool (Chang et al. 2020).Developed under the Open Universe initiative, this tool can identify blazars and construct their multiwavelength spectral energy distributions (SEDs).It accesses a variety of catalogs, thereby providing blazar characteristics across multiple wavelengths.For this work we used VOU-Blazars V2.00, an evolution of the original version offering several new datasets, which has been implemented in the Firmamento platform2 (Tripathi et al. 2023).NEOWISE individual observations are typically short and span one to a few days period approximately every six months.To increase the signal to noise ratio we averaged fluxes that were separated in time by less than 10 days.From the resulting IR flux measurements, we computed the overall average flux, variability amplitude, and fractional variability.
We implemented a series of exclusion criteria to refine our dataset, focusing only on the most reliable sources.We excluded sources based on the following criteria: • Sources with a γ-ray band detection significance less than 5σ, ensuring we only consider sources with reliable γ-ray detection and good signal-to-noise ratio.
• Sources with the IR flux that is close to the NEOWISE flux limit, specifically those with an averaged νF ν flux less than 1.8 × 10 −13 erg cm −2 s −1 or a maximum flux less than 3.0 × 10 −13 erg cm −2 s −1 .
• Sources potentially contaminated by the host galaxy in the IR band, as indicated by an IR slope greater than or equal to 1.5, and fractional variability less than 0.12.
• Sources where the IR band might be dominated by IR torus emission, particularly those with both an IR slope and fractional variability less than 0.1.
By applying these exclusion criteria the final sample consists of 1,109 sources.Within this sample, 636 sources are identified as LBLs and 473 as IHBLs.Furthermore, the sample includes 348 FSRQs, 502 BL Lacs, and 252 blazars of an uncertain type.The sky distribution of these blazars, in Galactic coordinates, is depicted in Fig. 1, where the positions of LBLs and IHBLs are respectively indicated with blue squares and orange circles.

SYNCHROTRON PEAK CHARACTERIZATION USING WISE/NEOWISE DATA
The data from the WISE/NEOWISE surveys, enable the determination of the spectral slope of the IR emission.We show here that when this slope is primarily due to synchrotron radiation it can provide crucial constraints on our understanding of blazar properties.To compute the spectral slope in the IR band (slope IR ) in the ν-νF ν space, we use the fluxes at 3.4µm and 4.6µm from the NEOWISE surveys and apply the following formula: computing the errors in the slope using the error propagation method.As noted by Massaro et al. (2011) in their analysis of WISE colors, the value of the spectral slope (which is equivalent to a color) varies depending on the blazar type and defines distinct regions of the synchrotron spectrum.For instance, in FSRQs, the slope slope IR is typically less than zero, indicating that the IR band is located after the synchrotron peak, within the high energy cut-off.Conversely, in high-peaked BL Lacs, the slope is positive implying that the IR band is within the still rising part of the synchrotron component.As we will see in the following, this slope provides a reliable quantitative estimate of the peak of the synchrotron component, referred to as the W-peak (ν p W ). To explore the relationship between the synchrotron peak and the IR spectral slope, we examined the sample of bright blazars that have been frequently observed by Swift (Giommi et al. 2021), which consists of 43 well-known sources with many multi-frequency measurements and IR emissions predominantly from synchrotron radiation originating in the jet, therefore minimizing contributions from the host galaxy, infrared torus, or other jet-unrelated components.We computed the synchrotron peak of these selected blazars using BlaST (Glauch et al. 2022), a machine-learning estimator of the synchrotron peak of blazars based on multiwavelength data.We then divided the sample into two groups: blazars with a synchrotron peak below the infrared energy band and those with a peak above this threshold.
Fig. 2 plots IR slopes of the blazars in our sample against the synchrotron peak frequency, revealing a clear relationship between these two variables.We estimated the best fit to a linear relationship using only the sample of sources frequently observed by swift, which ensures that the SEDs are well populated both in terms of frequency and time.
These relationships can be expressed as: when log(ν peak ) < 13. when log(ν peak ) > 13.2.The linear relations between ν p W and IR slope are plotted in the left panel of Fig. 2 together with the data of the full sample (light gray points), and that of the sub-sample of blazars that have been frequently observed by Swift (red points).To enhance clarity and reliability, from the complete source sample, we only include sources for which the average IR band flux exceeds 5.0 × 10 −13 erg/cm 2 /s, comfortably above the flux limit, to assure both accurate slope estimates and reliable measurements.It is noticeable that although the overall trends in the source sample are consistent with the predictions made by Equations 2 and 3, the scatter in the gray points is somewhat larger than that in the sample of Swift sources.This is most likely due to the fact that the SEDs used to estimate the ν p W of gray sources often include sparse multi-frequency data taken from single or very few observations carried out when the source happened to be in a random state (including flares or temporary faint periods), whereas the red sources are more representative of their mean long-term state, since their NEOWISE, X-ray and optical fluxes are averaged over many observations carried out over a long period of time.
In addition to the peak frequency of the synchrotron component, another crucial parameter is the flux at the synchrotron peak.In the fourth catalog of AGN detected by Fermi-LAT (Ajello et al. 2022) the synchrotron flux at the peak (νF ν syn ) was estimated by fitting the data with a third-degree polynomial function.Data collected during flaring episodes or those dominated by thermal emission from the accretion disk or the host galaxy were excluded from this analysis.To determine a correlation between the ν p W and the synchrotron peak flux we fit the following functions to the data of our sample.log(νF For sources where ν p > 6.3 × 10 13 Hz, and log(νF in sources where ν p < 6.3 × 10 13 Hz. The best fit gives θ h = 10.79 • and θ l = 7.97 • , respectively.This functional relationship facilitates the calculation of the anticipated value (or range) of νF ν syn , provided the averaged flux in the IR range and the peak frequency are known.In the right panel of Fig. 2, we compare two different estimates of peak flux: one derived from polynomial fitting and the other obtained from the average IR flux.Sources frequently observed by the Swift observatory are highlighted in red, while the entire sample is shown in light gray.A linear correlation is evident between the two flux estimates: the correlation coefficient below and above the infrared threshold of 6.3 × 10 13 Hz being r = 0.83 and r = 0.87.Additionally, an inset graph plots the ratio of the peak flux estimated through polynomial fitting to the flux deduced from the average IR band.This ratio is plotted against the peak synchrotron frequency, ν p W . Remarkably, this ratio approaches unity across the full range of observed synchrotron frequencies.This finding suggests that one can reliably estimate both the peak frequency and the peak flux of the synchrotron component by simply measuring the spectral slope and average flux in the IR band.To illustrate some of the results of measuring the peak frequency and flux derived from the IR data, the upper panel of Figure 3 presents the broadband SEDs of 1H 1013+498 and ON 231 (in blue) together with the predicted νF ν syn and ν p W (depicted in red).For νF ν syn , the expected range is calculated considering the minimum and maximum IR flux values (red arrows).For 1H 1013+498 and ON 231, the peak synchrotron frequency and flux are predicted to be at 10 16 Hz and 1.19 × 10 −11 erg cm −2 s −1 , and at 2.51 × 10 14 Hz and 1.47 × 10 −11 erg cm −2 s −1 , respectively.These values visually correspond to the anticipated shape of the synchrotron component.Conversely, for the cases of PMN J0151-3605 and 4C +20.25 (lower panel of Figure 3), the predicted frequency and flux are noticeably shifted towards higher frequencies, likely because part of the IR flux is due to the host galaxy.For both sources, ν p W ≃ 2 × 10 15 Hz but a synchrotron component peaking at this frequency would exceed the observed X-ray data.Nevertheless, for a large number of sources, as demonstrated in Figure 2, the non-thermal IR data enables a robust prediction.
The formulae presented in Equations 2 and 3, which link the synchrotron peak frequency to the IR slope and the correlation between the average IR flux to the synchrotron peak frequency (Equations 4 and 5), provide a valuable and accessible method for estimating the synchrotron component of blazar jets.This is especially useful when studying large samples of blazars, where obtaining high-quality data for all sources may not be practical, such as in the X-ray or γ-ray bands.However, it is crucial to note that the ν p W and peak flux estimates are based on long-term averaged infrared data and represent the average frequency and intensity of the synchrotron peak of an object.As a result, this method cannot be used to determine the peak value during flaring events or low-level states.

IR-GAMMA-RAY RELATIONSHIP
With the availability of γ-ray and IR data for extensive blazar samples, it is now feasible to investigate potential correlations among various parameters.In this section, we explore the relationships between different quantities estimated in the IR and γ-ray bands within the blazar sample that meets our selection criteria.

IR flux versus the γ-ray flux
In Figure 4, we study the relationship between IR and γ-ray emissions using the plane defined by the energy flux (F γ ) in the 0.1-100 GeV range (from 4FGL-DR3) and the mean IR flux at 4.6 µm (F 4.6µm ), calculated by averaging observations performed within a time separation of less than 10 days.The left panel shows the properties of FSRQs and BL Lacs, while the right panel of LBLs and IHBLs.For FSRQs and BL Lacs, there is a strong correlation between log10(F γ ) and log10(F 4.6µm ), with coefficients of r = 0.66 and r = 0.80 respectively, and a probability of p < 10 −5 .Fitting a linear function to the BL Lac sample in the form: yields κ = 0.81 and F 0 = −1.51.A similar fit for the FSRQ sample gives κ = 0.78 and F 0 = −1.72.The positive index points to a direct correlation between the average IR flux and the γ-ray flux for both classes, but it is more pronounced in BL Lacs than in FSRQs.The same relationship between LBLs and IHBLs is depicted in the right panel of Fig. 4. For LBLs, the correlation coefficient is r = 0.65, whereas for IHBLs, it is r = 0.82.These correlation coefficients mimic the results presented for FSRQs and BL Lacs.Considering the entire sample of sources, the correlation remains significant with a coefficient of r = 0.70 implying a statistically significant correlation between IR and γ-ray fluxes of blazars.Similar results have been obtained in the past by e.g.D'Abrusco et al. ( 2012) in smaller samples of blazars.

IR slope versus γ-ray slope
It is well-known that distinct spectral properties characterize different types of blazars.To investigate the correlation between the IR and γ-ray bands, we performed a statistical analysis of their slopes.The γ-ray slope, derived from the power-law index, is defined as slope γ = 2 − Γ, where N (E γ ) ∼ E Γ γ .We restricted our discussion to sources with a relative uncertainty on the IR slope of less than 50% (i.e., slope IR /slope IR,err > 2.0), to exclude sources with large uncertainties.Fig. 5 presents the scatter plot of slope IR − slope γ .The distributions of BL Lacs and FSRQs (left panel) extend towards opposite ends of the scatter plot, indicating distinct trends in their spectra.However, a notable overlap in the central region suggests that the two populations share common properties within a specific range.This is expected since some BL Lacs are LBLs, but they exhibit slightly harder γ-ray photon indices.
The opposite ends of the BL Lac and FSRQ distributions suggest distinct spectral indices for each.Specifically, BL Lacs are primarily characterized by indices slope IR > 0.1 and slope γ > −0.1, while FSRQs predominantly exhibit slope IR < 0 and slope γ < 0. To compute the density of points in the plot, we employed the Kernel Density Estimation (KDE).This analysis revealed that the region of highest density for BL Lacs lies at slope IR = 0.27 and slope γ = −0.06,whereas for FSRQs, the values are slope IR = −0.53 and slope γ = −0.43.The contrast between these values underscores the differences between the BL Lac and FSRQ distributions.Furthermore, a correlation analysis of slope IR − slope γ shows a stronger overall correlation for BL Lacs (0.62) compared to FSRQs (0.38).Results in line with these findings have been previously reported by D' Abrusco et al. (2012).A linear function fit of the form slope γ = κ × slope IR + α 0 to the BL Lacs distribution yields κ = 0.33 and α 0 = −0.16.
In the slope IR − slope γ plane (shown in Fig. 5, right panel), the distinctions between LBLs and IHBLs is more pronounced.The region of highest density for LBLs is centered at slope IR = −0.48 and slope γ = −0.39,while for IHBLs, it is at slope IR = 0.25 and slope γ = −0.07.Both LBLs and IHBLs display moderate correlations between the IR and γ-ray slopes, with coefficients of 0.43 and 0.48, respectively.Notably, the correlation coefficient for IHBLs is lower than that of BL Lacs.This discrepancy is due to the inclusion of blazars of uncertain type (bcu) in the IHBL classification.These bcu blazars possess distinct properties, influencing the overall correlation.Fitting the slope IRslope γ relation yields parameters κ = 0.40 and α 0 = −0.23 for LBLs, and κ = 0.39 and α 0 = −0.20 for IHBLs, which reveal a positive correlation between the slope IR and α γ slopes.

Fractional variability in the IR and γ-ray bands
In this subsection, we explore and quantify the variability of selected blazars in the IR and γ-ray band.The fractional variability (F var ) is computed using the standard formula as described in, for instance, Schleicher et al. (2019).This computation considers IR flux measurements at 4.6µm and its associated uncertainty, taken with a six-month time separation between NEOWISE surveys.
The comparison of fractional variability between BL Lac objects and FSRQs is illustrated in Fig. 6 (left panel).The FSRQ distribution is noticeably shifted toward higher values compared to the BL Lac distribution, indicating that FSRQs, on average, exhibit more pronounced IR variability than BL Lacs.The mean value for the FSRQ distribution is 0.66, while that for BL Lacs is 0.35.A comparable trend is observed between LBLs and IHBLs as presented in Fig. 6 (right panel).Specifically, the mean value for the IHBL distribution is 0.30, while for LBLs, it is 0.58.
It is interesting to examine the relationship between the fractional variability in the IR and γ-ray bands.In the 4FGL catalog of AGNs (4LAC-DR3; Ajello et al. 2022) the fractional variability for blazars is calculated using fluxes estimated annually from 2008 to 2020.In this work, we consider only γ-ray fractional variability measurements with a relative uncertainty of less than 50 percent.Fig. 7 (left panel) shows the IR and γ-ray fractional variability for both BL Lacs and FSRQs.The correlation coefficient for BL Lacs is 0.38, pointing to a moderate correlation between the fractional variability in the IR and γ-ray bands.In contrast, FSRQs display a weaker correlation with a coefficient of 0.30.Further, Fig. 7   frequencies are related to the properties of the relativistic electrons in the jet and the surrounding environment, and their comparison can reveal essential information about the physical processes and the dominant emission mechanisms in these sources.ν p IC is a new parameter introduced in the 4LAC-DR3 catalog of AGNs (Ajello et al. 2022) that is estimated directly from the curvature of the γ-ray spectrum.For our analysis, we have selected only those blazars for which the relative uncertainty in the ν p IC estimation is below 50%.We note that if the SSC mechanism is responsible for the emission, then in the Thomson regime, the peak frequencies of synchrotron and inverse Compton are related through the relation γ Giommi et al. 2012).In our case the ν p IC estimation based on the curvature of the γ-ray spectrum could introduce biases due the limited γ-ray band, e.g.difficulties in estimating peaks at or below the Fermi-LAT low-energy threshold.Additionally, issues related to the sensitivity dependence on spectral slope, may contribute to discrepancies between our observed correlation and theoretical expectations.
In Fig. 8 (left panel), we present the comparison of ν p W and ν p IC for FSRQs and BL Lacs.A majority of the FSRQs are found in the region where ν p W < 10 13 Hz and ν p IC < 2 × 10 23 Hz, while BL Lacs display a broader distribution.The highest density region for FSRQs is located at (ν p W = 4.02 × 10 12 Hz, ν p IC = 5.2 × 10 22 Hz), while for BL Lacs it is at (ν p W = 5.6 × 10 15 Hz, ν p IC = 1.8 × 10 24 Hz).We fit a linear function log(ν p W ) = κ × log(ν p IC ) + ν 0 to the data, obtaining κ = 0.30 and ν 0 = 19.6 for BL Lacs, and κ = 0.20 and ν 0 = 20.3 for FSRQs.The correlation coefficient between log(ν p W ) and log(ν p IC ) is 0.74 for BL Lacs, indicating a strong correlation while it is only 0.21 for FSRQs.This is expected as in BL Lacs, ν p IC is determined by the SSC component, which in turn is defined by the synchrotron component.
The comparison between LBLs and IHBLs is displayed in the right panel of Fig. region for IHBLs is located at (ν p W = 5.4 × 10 15 Hz, ν p IC = 1.9 × 10 24 Hz), while the peak of the LBL distribution is at (ν p W = 4.7 × 10 12 Hz, ν p IC = 7.0 × 10 22 Hz).The correlation analysis between log(ν p W ) and log(ν p IC ) yields a coefficient of r = 0.34 for LBLs, corresponding to a weak to low-moderate correlation, while r = 0.58 for IHBLs.The linear fit to the data produces κ = 0.67 and ν 0 = 14.4 for LBLs, and κ = 0.37 and ν 0 = 18.4 for IHBLs.

Prediction of VHE detectability
The ability of estimating the location and intensity of the peak of the synchrotron emission, combined with the correlation between ν p W and ν p IC and the γ-ray spectral slope, can be used to predict the intensity in the VHE band, hence the detectability of known blazars with current and planned TeV telescopes.In the following we build the ν p f (ν p ) distribution of TeV detected IHBL blazars, as reported in the TeVcat catalog3 , as well as that of blazars that have been observed but not detected by the VERITAS γ-ray observatory (Archambault et al. 2016).Both samples are certainly not ideal statistical data sets as they are likely biased in various ways, the first one probably in favour of sources observed during flares or with long exposures, and the second one likely suffering the opposite bias, that is shorter observations in non-ideal conditions or sources pointed during low states.Nevertheless, these are the only samples of VHE observed blazars easily available, they are relatively sizable and the opposite biases tend to partially cancel out.Fig. 9 shows the distributions of peak flux intensities for the subsamples of detected (red histogram) and undetected IHBL blazars (blue histogram), which include 48 and 49 sources respectively.The histogram of detected blazars is notably shifted towards fluxes approximately three times larger than those of undetected objects.This suggests that the observed difference cannot be attributed to variations in observational conditions, thereby indicating that the synchrotron peak intensity of IHBL blazars serves as a valid indicator of VHE detectability.Fig. 9 illustrates that for peak fluxes with logarithms larger than -10.75 all sources are detected.At lower ν p f (ν p ) fluxes the detection percentage gradually decreases as indicated by the numbers above each bin.In the next subsection we use these percentages to estimate the number of blazars that could be detected in an hypothetical Cherenkov Telescope Array (CTA) extragalactic survey starting from the catalog of known blazars.

A preview of a mini CTA extragalactic survey
The histograms shown in Fig. 9 are based on data from the current generation of IACTs.Assuming that the planned CTA extragalactic survey will be approximately three times more sensitive than current observatories, then the detection percentages shown Fig. 9 should be associated to peak fluxes that are three times fainter, or approximately 0.5 decades lower in log space.To define a survey with sensitivity similar to the CTA extragalactic survey that is sufficiently compact to be considered in this work, but still numerically significant, we consider all known blazars that are located in the sky region delimited by 0 • < dec < 5 • and Galactic latitude (-b-) larger than 30 • , corresponding to an area of 1,100 square degrees.This is equivalent to slightly more than one tenth of the planned CTA extragalactic survey.This part of the sky includes 108 IHBL blazars.For each of these sources we constructed the SED using the Firmamento platform (Tripathi et al. 2023) and, whenever possible, we calculated ν p and ν p f (ν p ) using IR data and  the relationships described above.In those cases where the NEOWISE data are significantly contaminated by the host galaxy or are not available we used the BLAST tool to estimate ν p and calculated ν p f (ν p ) by fitting the available multi-frequency data near the peak.To estimate the probability of detection of each of the 108 blazars in our mini survey, we used the distribution of detection fractions shown in Fig. 9 after shifting the Log(ν p f (ν p )) axis by 0.5.The results are reported in table 1 where column1 gives the interval of ν p f (ν p ) values considered, the second column gives the probability of detection, the third column gives the observed number of sources in our sample with ν p f (ν p ) in the bin, the fourth column gives the expected number of detections (equal to the observed number multiplied by the detection probability), and the last column gives the expected number of detection scaled to 10,000 sq degrees, the assumed area of the CTA extragalactic survey.The expectation of a total of 158 detections of IHBL blazars in the CTA survey is somewhat higher but of the same order of an earlier estimation presented in Cherenkov Telescope Array Consortium et al. (2019).The detectability of LBL blazars cannot be estimated with our method since these sources are characterised by very different and extremely variable γ-ray to νF ν peak flux ratios, a parameter also known as the Compton dominance.Since the TeVcat catalog includes approximately 15 percent of LBL blazars, if we assume that the same fraction of LBLs will be detected in the CTA extragalactic survey, the total number of blazars in our simulated CTA extragalactic survey increases to about 180.

DISCUSSION AND CONCLUSION
We have investigated the multi-wavelength properties of a large sample of blazars using data from the NEOWISE extragalactic sky survey, hundreds of Swift X-ray observations, and twelve years of γ-ray data from the Fermi-LAT survey.The γ-ray sample consists of 1,109 blazars also monitored at IR frequencies.Examining the IR and γray emissions of this extensive sample provides important constraints on their emissions and offers valuable insights into the underlying physical processes occurring in their jets.These insights can help improve our understanding of the dominant emission mechanisms in the IR and γ-ray bands, jet properties, and aid in the development of more precise blazar classification schemes.Our analysis shows that the IR slope, calculated between the 3.4µm and 4.6µm wavelengths for the subset of blazars whose IR emissions originate from the jet, exhibits a tight linear relationship with the position of the synchrotron peak.Furthermore, the mean 4.6µm IR flux enables the estimation of the peak intensity of the synchrotron component.In other words, the IR slope and flux can be used to estimate properties such as the power content and the maximum acceleration energy in blazars.This method has been implemented in a tool called W-Peak, which is publically available within the Firmamento platform (Tripathi et al. 2023) and can be applied in various contexts to help derive valuable information suitable for blazar selection as well as providing insights on their underlying physical processes.For instance, the connection between the IR and synchrotron peak enables a more efficient selection and prioritization of blazars for multiwavelength and multimessenger studies.
The relationship between the IR slope and the synchrotron peak frequency, as quantified by Eqs. 2 and 3, can be interpreted within the context of the physical mechanisms responsible for blazar jet emission.The radio to optical/Xrays emission from the jet is primarily governed by synchrotron radiation generated by accelerated electrons spiraling around magnetic field lines.The spectral slope in the IR band is particularly sensitive to the shape of the synchrotron spectrum, which is determined by two main factors: the energy distribution of the radiating particles and the strength of the magnetic field.A flat (close to 0) or positive IR slope in the SED implies an energetically dominant contribution from high-energy particles.Conversely, a steeper (more negative) slope suggests fewer high-energy particles in proportion.The synchrotron peak frequency, which marks the frequency where the synchrotron emission is maximum, is the outcome of the balance between the number of high-energy and low-energy particles (i.e., the break energy).So the IR slope and average flux link the break energy in the electron distribution with the magnetic field strength, thereby constraining both these parameters.
Our analysis also reveals a statistically significant correlation between the fluxes in the IR and γ-ray bands.This correlation is especially pronounced for BL Lacs, which can be explained by considering the underlying physical processes.The emission in the IR band arises from the synchrotron radiation of relativistic electrons, while the emission in the γ-ray band results from inverse Compton scattering of synchrotron photons by the same population of electrons (Ghisellini et al. 1985;Maraschi et al. 1992;Bloom & Marscher 1996).Any change in the energy distribution of the electrons or the conditions within the jet will impact both the synchrotron (IR) and SSC (γ-ray) components.For example, an increase in the density of high-energy electrons or a stronger magnetic field will lead to an increase in the synchrotron radiation, which will be observed as an enhanced IR emission.Concurrently, the increased synchrotron photon density will lead to a higher rate of self-Compton scattering and a stronger γ-ray emission.In contrast, in FSRQs, several external photon fields (e.g., broad-line regions and dusty torus) can serve as a significant source of external photons.These photons can be scattered by the relativistic electrons in the jet, and the γ-ray spectrum will be influenced by a combination of synchrotron and external photons (Sikora et al. 1994;Dermer et al. 1992;Dermer & Schlickeiser 1994;Sikora et al. 1994;B lażejowski et al. 2000).The presence of these external photon fields introduces additional complexities.As a result, the correlation may weaken or even be hidden, explaining the less pronounced relationship observed in FSRQs.This can also explain the correlation between the peak of the synchrotron and high energy components obtained for BL Lacs or IHBLs.The correlations between IR emissions and other energy bands across blazar SEDs provide compelling evidence for the effectiveness of a one-zone SSC radiation model in describing the average broadband SEDs of BL Lacs.Specifically, the observed correlations between between slope IR and ν p W and slope γ , between log 10 (F γ ) and log 10 (F 4.6µm ), and between ν p W and ν p IC strongly suggest that the emissions in the investigated bands originate from a dominant population of electrons, largely confined to a single emission region.In contrast, in the case of emission from multiple independent regions (as in multi-zone models), one would not expect strong correlations.
The amount of variability in the IR band varies among the type of the sources being studied.On average, FS-RQs/LBLs exhibit greater fractional variability in the IR band than BL Lacs/IHBLs.Several factors might account for this observed difference in fractional variability.These include differences in particle acceleration rates between blazar subclasses or variations in the structure and magnetic field configuration of the jets in FSRQs and LBLs.However, the primary factor might arise from the physical process underpinning the emission in the IR band.The IR spectral shape of FSRQs and LBLs is soft, so it characterizes the higher-energy tail of the synchrotron emission, implying that the emission is due to electrons that are more energetic than those that are in balance between emission and acceleration.These high-energy electrons cool quickly and their number may fluctuate due to acceleration instabilities near maximum energy.In contrast, the more moderate amount of IR variability observed in BL Lacs and IHBLs is likely related to their harder IR slope, which characterizes the increasing shape of the synchrotron component, where the emission is due to a more stable population of electrons with energy well below where maximum acceleration occurs.
Finally, we have shown that the synchrotron peak intensity of IHBL blazars is an indicator of detectability at VHE γ-ray energies.Based on this we identified objects in current blazar catalogues that could be detected by the upcoming CTA observatory, and conclude that the CTA extragalactic survey could include approximately 180 blazars.

Figure 1 .
Figure 1.Hammer-Aitoff projection in Galactic coordinates of the sky distribution of Fermi-LAT detected blazars with NEOWISE data.

Figure 2 .
Figure 2. Left panel: NEOWISE slope versus the synchrotron peak estimated using BlaST.The best-fit line is shown in red.All source sample is shown in light gray.Right panel: The prediction of the synchrotron peak flux based on the averaged IR flux, compared with estimations from the polynomial fitting.The inset shows the ratio between the prediction and estimation versus the ν p W .In both plots, the sub-sample of blazars that has been frequently observed by Swift is shown in red.

Figure 6 .
Figure 6.Distribution of fractional variability for BL Lacs and FSRQs (left panel) and for LBLs and IHBLs (right panel).

Figure 7 .
Figure 7.Comparison of fractional variability in the IR band with the γ-ray band for BL Lacs and FSRQs (left panel) and for LBLs and IHBLs (right panel).

Figure 9 .
Figure 9. Synchrotron peak intensity histograms of VHE detected (from the TeVCat catalog, red) and undetected (from the VERITAS list of upper limits, blue) IHBL blazars.The number shown at the top of each bin represents the percentage of VHE detected objects in that bin.
8. It is evident that LBLs and IHBLs occupy distinct regions, with IHBLs populating the area characterized by high ν p W and ν p IC .The highest density Comparison of ν p W and ν p IC .Left panel: FSRQs versus BL Lacs.Right panel: LBLs versus IHBLs.

Table 1 .
Results of the mini survey and extrapolation to the full CTA extragalactic survey