TELAMON: Effelsberg monitoring of AGN jets with very-high-energy astroparticle emission

Aims. We introduce the TELAMON program which is using the E ﬀ elsberg 100-m telescope to monitor the radio spectra of active galactic nuclei (AGN) under scrutiny in astroparticle physics, speciﬁcally TeV blazars and candidate neutrino-associated AGN. Here, we present and characterize our main sample of TeV-detected blazars. Methods. We analyzed the data sample from the ﬁrst ∼ 2.5yr of observations between August 2020 and February 2023 in the range from 14GHz to 45GHz. During this pilot phase, we observed all 59TeV-detected blazars in the Northern Hemisphere (i.e., Dec > 0 ◦ ) known at the time of observation. We discuss the basic data reduction and calibration procedures used for all TELAMON data and introduce a sub-band averaging method used to calculate average light curves for the sources in our sample. Results. The TeV-selected sources in our sample exhibit a median ﬂux density of 0.12Jy at 20mm, 0.20Jy at 14mm, and 0.60Jy at 7mm. The spectrum for most of the sources is consistent with a ﬂat radio spectrum and we found a median spectral index ( S ( ν ) ∝ ν α ) of α = − 0 . 11. Our results on ﬂux density and spectral index are consistent with previous studies of TeV-selected blazars. Compared to the GeV-selected F-GAMMA sample, TELAMON sources are signiﬁcantly fainter in the radio band. This is consistent with the double-humped spectrum of blazars being shifted towards higher frequencies for TeV-emitters (in particular for high-synchrotron peaked BL Lac type objects), which results in a lower radio ﬂux density. The spectral index distribution of our TeV-selected blazar sample is not signiﬁcantly di ﬀ erent from the GeV-selected F-GAMMA sample. Moreover, we present a strategy to track the light curve evolution of sources in our sample for future variability and correlation analysis.


Introduction
Blazars are radio-loud active galactic nuclei (AGN) hosting relativistic jets pointed close to our line of sight.Their emission is highly beamed and Doppler-boosted, which makes them variable broadband emitters from radio to γ-ray energies.Their spectral energy distribution (SED) typically shows a two-peaked (double-humped) spectrum.The first peak corresponds to synchrotron emission while the second peak is often attributed to inverse Compton scattering.High-synchrotron peaked BL Lac type objects (HBLs) are defined as sources whose primary (synchrotron) emission hump peaks above 10 15 Hz in νF ν scale (Padovani & Giommi 1995).In the most extreme cases, they peak at even higher frequencies by up to two orders of magnitude for the so called extreme HBLs (EHBLs, Ghisellini 1999;Biteau et al. 2020).The second, high-energy peak of blazar SEDs can stretch into the very-high-energy (VHE) regime at TeV γ-rays.Imaging Air Cherenkov Telescopes have been able to detect VHE emission from ∼80 blazars, with the majority of them being HBLs1 .Most of these sources are faint radio sources, which makes TeV-emitting blazars difficult to study in the radio band.Blazars are of utmost interest for astroparticle physics as potential sources of ultrahigh-energy cosmic rays and neutrinos (e.g., Hillas 1984;Mannheim 1995).In particular, HBLs and EHBLs have been considered in some recent theoretical works as relevant neutrino sources (Tavecchio et al. 2014;Padovani et al. 2015;Giommi et al. 2020).We have established the TeV Effelsberg Long-term AGN Monitoring (TELAMON) program in August 2020, a radio monitoring program which uses the Effelsberg 100-m telescope to investigate the radio properties of TeVemitting blazars.We perform radio observations of a sample of TeV-selected blazars at high frequencies up to 44 GHz to trace dynamical process in these objects.Our program is designed to monitor the radio spectra of TeV-blazars and candidate neutrino AGN.A first study of TeV-selected HBLs and EHBLs has been presented by Lindfors et al. (2016) at 15 GHz.We extend their observations with radio coverage at multiple frequencies, and, for the first time, a spectral characterization of these objects.In Section 2, we explain in detail the sample selection and observing setup used with the Effelsberg 100-m telescope.On top of that, we present the TELAMON analysis pipeline including cross-scan setup and calibration procedures.In Section 3, we present first results for sources in our sample and present average flux densities and spectral indices for all observed sources.These results are discussed and compared with previous studies, e.g., the Lindfors et al. (2016) study, in Section 4.Moreover, we compare source properties of the TeV-selected sources in our sample with the GeV-selected F-GAMMA (Fuhrmann et al. 2016;Angelakis et al. 2019) sample and discuss similarities and differences.In Section 5, we give an overview about planned future publications and the overall relevance of the TELAMON project within current developments in astroparticle physics.

Program description & sample
Our sample consists of all 59 TeV-detected blazars 1 in the Northern Hemisphere (i.e., Dec. > 0 • ).In addition, we include 5 TeV-detected blazars from the Southern Hemisphere and a TeV-detected radio galaxy (3C 264) as sources of special interest.This TeV-selected sample, consisting of 65 sources, is presented in Table 1.It includes 16 EHBLs (i.e., ν peak > 10 17.38 Hz = 1 keV), 32 HBLs (i.e., 10 15 Hz < ν peak ≤ 10 17.38 Hz), 4 intermediate-peaked BL Lac type objects (IBL, i.e., 10 14 Hz < ν peak ≤ 10 15 Hz), 4 low-peaked BL Lac type objects (LBL, i.e., ν peak ≤ 10 14 GHz), 8 flat-spectrum radio quasars (FSRQ) and 1 Radio Galaxy (RG).For J0316+4119 (IC 310), a clear classification is complicated since the source shows some blazar-like properties (Kadler et al. 2012) but its jet might be misaligned by 10−20 • from the line of sight (Aleksić et al. 2014b).We include the source in our HBL count, according to its peak frequency found in the literature.Note that our classification for BL Lac type objects solely depends on the synchrotron peak frequency, as per the references listed in Table 1.We adapt the extreme synchrotron limit (i.e., ν peak > 1 keV) from Biteau et al. (2020), but do not account for the other types of extreme behaviours introduced in their work.Therefore, in our analysis, some of the blazars labeled as 'extreme' by Biteau et al. (2020) fall into the HBL category in our work.On the other hand, not all EHBL sources from our analysis are covered in Biteau et al. (2020).

Observations
For our observations, we are using the Effelsberg 100-m telescope operated by the Max-Planck-Institute for Radio Astronomy in Bonn, Germany.The observations presented here are conducted with the 20 mm, 14 mm and 7 mm receivers mounted in the secondary focus, in continuum observing mode with the backend "dual-spec-OPTOCBE".This results in four sub-bands (centered at 19.25 GHz,21.15 GHz,22.85 GHz and 24.75 GHz) for the 14 mm receiver, and four sub-bands (centered at 36.25 GHz,38.75 GHz,41.25 GHz,43.75 GHz) for the 7 mm receiver, which we are using since August 2020.Since spring 2021, the 20 mm receiver has been added yielding two additional frequency bands (centered at 14.25 GHz and 16.75 GHz).Note that we also started using the "SpecPol" receiver backend to record polarization information at the same time (cf.Heßdörfer et al. 2023).The polarization analysis and first results will be assessed in a separate publication.All receivers are equipped with two horns.The first horn is pointed directly at the target and passes the signal directly to the receiver.The second horn is pointed at the atmosphere off-source and is used to subtract weather effects from the first horn.We apply this weather subtraction to all of our observations, because at the used frequencies Earth's atmosphere (especially on cloudy days) emits thermal radiation that needs to be accounted for.In order to measure the flux density of a source, "cross-scans" are performed on the targets, consisting of typically 8 sub-scans (4 in azimuth and 4 in elevation) at 20 mm, 14 mm and 7 mm (36.25 GHz,38.75 GHz).For the higher 7 mm frequencies (41.25 GHz,43.75 GHz), 16 sub-scans (8 in azimuth and 8 in elevation) are used.During a cross-scan the telescope slews over the point-like source region in azimuth-and elevation-direction multiple times while measuring the antenna temperature of the receiver.About every four hours, a calibrator source is observed in order to focus the telescope and to extract calibration factors (cf.Sect.2.4).
For this work, we consider all observations from within the first 2.5 years of the program, i.e., from August 2020 to February 2023.This comprises data taken in 95 observing sessions (epochs) with a total observing time of 1160.3 hours.Note that part of the observing time was also used for observations of neutrino-candidate sources (cf.Kadler et al. 2021), which will be discussed in a separate publication.During this time, we have observed 65 distinct TeV-sources in total, 53 at 20 mm, 61 at 14 mm and 45 at 7 mm.The number of observations per source and per receiver varies within the sample since we have tested different observing strategies depending on the source flux density.Our optimized monitoring strategy has been implemented since mid-2022 (cf. Sec. 4.3).In order to select the statistically most complete sample from our observed targets, we will from here on only address the 59 Northern (i.e., Dec. > 0 • ) TeVblazars.Results for the 5 Southern sources and the radio galaxy 3C 264 (J1145+1936) are presented in Appendix B. Plots of the latest available light curves and spectra are publicly available on the dedicated TELAMON website2 .b For some sources no conclusive redshift value is available in the literature.See also Foschini et al. (2022) for a more detailed collection of redshifts for most sources including information about their reliability.*Sources overlapping with Biteau et al. (2020).† Observations are coordinated with the MOMO program (Komossa et al. 2023).‡ The adopted redshift is a statistical estimate, which may deviate from the z = 0.31, previously suggested based on marginal detection of the host galaxy by Nilsson et al. (2008).

Data acquisition & reduction
The general data analysis for pointed flux density measurements with the Effelsberg 100-m telescope has been described in detail by Angelakis et al. (2019) in context of the F-GAMMA program.We use a very similar data reduction procedure and will therefore restrict this section to the changes and improvements compared to the analysis of Angelakis et al. (2019).The major improvements presented here are a semi-automated flagging algorithm, a new calibration procedure, and an in-depth measurement uncertainty discussion.In principle, this section explains how the raw data output of every scan, namely antenna temperatures, is converted into astrophysical units of jansky.The antenna temperatures are calculated by using a noise diode switching system (cf.Müller et al. 2017).Left-handed circular polarization (LCP) and right-handed circular polarization (RCP) are averaged.

Sub-scan fitting
As mentioned in Sect.2.2, every scan of a source at a specific frequency consists of multiple sub-scans.These sub-scans are used to calculate one flux density value for each frequency per scan.For every scan, all azimuth sub-scans are used to generate a single average azimuth scan and all elevation sub-scans are used to generate a single average elevation scan.Since all of the sample sources are assumed to be point-like and the primary beam/point spread function (PSF) of the telescope is well described by a 2D Gaussian, a Gaussian curve is fitted to both average scans.This fitting process is carried out with the toolbox software (Kraus et al. 2003).At this point, it is necessary to perform a data quality check to see if all sub-scans can be considered clean scans.This is important, since some sub-scans might be corrupted by radio frequency interference (RFI), telescope errors, or atmospheric effects, and therefore might lead to errors in the derived flux density values.
In order to filter out corrupted sub-scans, a flagging system has been developed, which uses multiple criteria to detect corrupted scans.The detailed flagging criteria are listed in Appendix A. If an averaged scan is flagged according to these criteria, it is not trivial to judge whether the source was below the detection limit (not visible in all sub-scans) or if there were simply a few corrupted sub-scans without checking every single subscan for every source.As done by Angelakis et al. (2019) a manual analysis of individual sub-scans is possible but very tedious, time-consuming, and also not clearly reproducible since it depends on the analyzer's individual eyes.Therefore, we have developed a semi-automated analysis tool which takes care of detecting and sorting-out corrupted sub-scans.The general principle of the algorithm is presented in Appendix A. The automated analysis quality is superior to the manual analysis since every sub-scan deletion can be tried out, instead of simply judging by the human eye what sub-scan deletion might improve the overall fit quality.On top of that, due to the appointed flagging criteria, it is completely reproducible.As a final check, the data are again inspected manually to sort out any outliers and left-over corrupted scans, which the algorithm was not able to detect.

Data corrections
All averaged scans that passed our semi-automated flagging algorithm undergo further data reduction processes.These data reduction are composed of the pointing offset correction, atmospheric opacity correction, and elevation-dependent gain correction presented by Angelakis et al. (2019) in their Sect.4. The Table 2: Gain curve parameters used for the different receivers of the Effelsberg 100-m telescope (second order polynom, cf.Eq. ( 5) of Angelakis et al. 2019).only difference is that we determine the zenith opacity τ using a water-vapor radiometer located at the focus cabin of the telescope by measuring the strength of the 22 GHz water-vapor line (Roy et al. 2004).On top of that, we use updated gain curve parameters A i as presented in Table 2.

Calibrator sources
In this section, it is discussed how the calibration factor for the analysis is determined.Baars et al. (1977) have introduced a set of secondary calibrators that can be used on a day-to-day basis for calibration purposes.Their flux density is very constant over long periods of time (several decades) and they are also very compact sources.These secondary calibrators have been further investigated and monitored by Ott et al. (1994) and most lately by Perley & Butler (2013, 2017).According to the latest publication (Perley & Butler 2017), the best suited calibrators in the TELAMON frequency range (14 GHz−44 GHz) are 3C 286 and 3C 295.They are therefore used as calibrator sources in this work.Perley & Butler (2017) provide parametrized spectra for these sources that are used to calculate the frequency dependent calibrator flux densities for our analysis.Since 3C 295 is quite faint (≲ 1 Jy) at higher frequencies (ν ≳ 35 GHz) and therefore not always detectable (especially during bad weather sessions), the sources NGC 7027 and W3(OH) are also included as calibrator sources.NGC 7027 is a planetary nebula and has been proposed as a secondary calibrator by Baars et al. (1977).Even though its flux density is gradually fading with time, this source can be used as a calibrator since its behavior has been well characterized by Zijlstra et al. (2008).The spectral model provided in this paper is used to calculate time-dependent calibration flux densities for NGC 7027.W3(OH) is a star forming region that exhibits a strong water maser in the 14 mm band (Hachisuka et al. 2006).Therefore, it is excluded for 14 mm calibration, but since the source is very bright at 7 mm (≈ 3 Jy) it is useful as a substitute for 3C 295 in this band.For W3(OH) we use our own calibration model which was created using Effelsberg archival data and assuming a free-free emission model.Due to their brightness, the two main calibrators observed with highest priority are 3C 286 and NGC 7027.3C 295 is used as a backup calibrator at 20 mm and 14 mm, and W3(OH) as a backup calibrator at 20 mm and 7 mm.For one epoch in 2021, we used 3C 138 as a calibrator for all bands, using the Perley & Butler (2017) model.

Calibration process
As mentioned in Sect.2.2, usually one of the secondary calibrator sources is observed every four hours during an observing session to determine the calibration factor.This calibration measurement typically consists of up to three scans per frequency on  1)), as used in the final analysis.The dark gray area represents the calibrator scatter σ sc,ν , while the light gray hashed area represents the total calibration uncertainty σ cal , including the model uncertainty σ model and σ sc,ν .the calibrator.Usually, more than one scan of the calibrator is taken during a calibration measurement, therefore the mean and standard deviation are calculated to get a value with uncertainty for the calibration factor at the given time.In the case of only one calibrator measurement per calibration measurement the uncertainty is estimated depending on the frequency (see Sect. 2.5).For each observing session, one therefore gets a calibration factor every ∼4 hours, in the ideal case.An example is given in Fig. 1, where all calibration factor data points (black dots) are averages of the underlying (multiple) calibrator scans at the time.When talking about calibration factors in the following, it is referred to these, already (sub-)averaged calibration factors like the (black) ones presented in Fig. 1.
In principle, the calibration factor only depends on the sensitivity of the telescope at a given frequency and should therefore be constant throughout an observing session (for a given frequency).As one can see in Fig. 1 this is not always the case.The calibration factor does change throughout the observing session, mostly due to temperature changes, which affect the ideal focus position.If the telescope is out of focus, the received signal is weaker and therefore the calibration factor rises.The telescope is being kept in focus by readjusting the focus every four hours (or more often in the case of significant temperature changes) to keep the calibration factor constant, but still the calibration factor varies throughout the observation.Modelling the behaviour of the focus is non-trivial, since it is unknown how much of the calibration factor fluctuation is statistical and how much is systematic due to, e.g., shifting focus.If the fluctuation was solely statistical, one would have to take the mean calibration factor throughout the entire session.If the fluctuation was solely systematic it would be best to linearly interpolate between the calibration factors.In the present case, however, interpolation between all data points would be an over-interpretation of the data and taking the mean cannot account for the fluctuation of the calibration factor.In order to take this fluctuation into account and also not to over-interpret the data points, the calibration factor is modelled by using a simple moving average (SMA) interpolation (Chou 1975).Essentially, this approach is a combination of both methods, since first, we take the mean between two adjacent calibration factor values and then interpolate between these mean values.If the calibration factor at time t i is Γ c,i , the interpolated values are calculated via where Γ c,i+1 is the calibration factor adjacent to Γ c,i .The same procedure is used for the time interpolation If initially there are n calibration factors Γ c,i , this interpolation will result in n − 1 interpolated calibration factors Γ SMA,i .As indicated by the blue line in Fig. 1, we interpolate linearly between the (t SMA,i , Γ SMA,i ) values to get a general expression for the calibration factor at any given time.For times t < t SMA,1 , the calibration factor is modelled as constant Γ SMA,1 .Analogously, for times t > t SMA,n−1 , the calibration factor is modelled as constant Γ SMA,n−1 .If there is only one or two calibration factors available (n = 1, 2), a constant calibration factor is assumed for the entire epoch.For all other cases, the interpolation (blue line) is used to determine the calibration factor at any given time during the observing session.This means, for every source scan, one calculates the corresponding calibration factor Γ c using the SMAinterpolation and then uses this calibration factor to calculate the flux density of the source S source via where T src is the observed source temperature.

Discussion of uncertainties
In this section, we discuss the determination of the total flux density uncertainties σ tot .The final uncertainty has to include the main uncertainties due to the sub-scan fitting process with the data corrections σ fit , and the uncertainty of the calibration factor, σ cal .First, we focus on the uncertainty due to sub-scan fitting and data corrections, σ fit .Its value is calculated through Gaussian error propagation from the fitting and data correction process explained in Sect.2.3.2.Note that the gain curve is assumed to be free of uncertainty here, since the accuracy of the gain curve is also reflected in the fluctuation of the calibration factors and therefore included in the calibration uncertainty σ cal .Usually, atmospheric corrections have the biggest impact on the data correction uncertainty, especially in sessions affected by bad weather.In total, σ fit is on average on the order of ∼ 1 %.
The main contribution of the flux density uncertainties comes from the calibration uncertainty, σ cal .There are two sources of uncertainty that determine the total calibration uncertainty.First, one needs to consider the fluctuation of the calibration factor σ sc,ν for each frequency and observing epoch.Secondly, one needs to account for the uncertainty of the underlying calibrator model σ model .As explained in the previous section, an SMA-interpolation is used to calculate interpolated calibration factors.To be conservative, the uncertainty of the interpolated calibration factors is assumed to be constant.Therefore, the rootmean-square deviation (RMSD) of all (non-interpolated) calibration factors for each frequency and observing epoch is used as the uncertainty Fig. 2: Illustrative calibration factor evolution at 41.25 GHz throughout the entire program.We determine the calibration factor scattering uncertainty σ sc,ν by taking the standard deviation of these values at each frequency.In this example, σ sc,ν is on the order of 10 %.Significant outliers can be explained by bad weather epochs.
Here, Γ c is the mean over the non-interpolated calibration factors Γ c,i .This uncertainty is also illustrated in Fig. 1 by the dark gray background.In the case where only one calibration factor was measured during an observing epoch, it is non-trivial to determine a sensible value for this uncertainty.In order to deal with this special case, the calibration factor scattering for each frequency from all observing epochs has been analyzed.An illustrative plot of the calibration factor evolution at 41.25 GHz throughout the entire program since August 2020 is shown in Fig. 2. In order to derive a sensible uncertainty for these epochs, the standard deviation of all calibration factors throughout the program is used to get an estimate of the average calibration factor scatter.This analysis is performed for every frequency band.Following this procedure, a conservative calibration uncertainty of 5 % is used at 20 mm and 14 mm and 10 % at 7 mm for epochs with only one calibration factor.In Fig. 2, one can also see that for some epochs the calibration factor varies more than for others.This is due to varying focus because of significant temperature and weather fluctuations (on shorter time scales than the usual focus adjustment interval every ∼4 hours) in some observing sessions.
In addition to the calibration factor scattering σ sc,ν , one also needs to take into account the accuracy of the calibrator models σ model .The models of 3C 295, 3C 138 and 3C 286 by Perley & Butler (2017) have an estimated accuracy of 3 % to 5 %, with the larger uncertainty at the lowest (∼ 50 MHz) and highest (∼ 50 GHz) ends.Since the observations take place from 14 GHz−44 GHz, which is at the higher end of their scale, an uncertainty of 5 % for the 3C 295, 3C 138 and 3C 286 models is assumed.For the timely variable model of NGC 7027, Zijlstra et al. (2008) provide an uncertainty of 6 %.It is estimated that the uncertainty for the W3(OH) model is also in the same range, since it has a similar underlying free-free emission model.As a conservative estimate, a general accuracy of the calibrator models of σ model /Γ c = 6 % is assumed, which is the maximum uncertainty out of the models.One must not use the Gaussian law of error propagation to combine the uncertainties of the different models, since they are not statistically independent and all based

Flux Density [Jy]
Raw Data Averaged Data Spectral Fit Fig. 3: Example of the power-law spectral fitting process (J1653+3945, also known as Mrk 501, on July 28, 2021).The magenta squares indicate the raw data, uncertainties are not shown to enhance readability.For every frequency, where more than one raw flux density value is available, the flux densities are averaged (black dots).The averaged flux densities are then fitted with a power-law spectrum (blue line).In the presented case, one finds α = −0.262± 0.021.on the same flux density scale by Baars et al. (1977).To calculate the total calibration uncertainty σ cal , the estimated model uncertainty σ model and the calibration scatter σ sc,ν , which is individual for every epoch and frequency, are added quadratically.The total calibration uncertainty is therefore given by In order to obtain the total flux density uncertainty σ tot , one needs to combine σ cal with the fitting uncertainty σ fit .This is done by Gaussian error propagation.Following Equation (3), one finds for the total flux density uncertainty.Note that for the following analysis the model uncertainty is considered to be zero, since it is purely systematic throughout the entire program.Using it would therefore lead to an overestimation of the statistical flux density uncertainty.If the presented flux density values were to be combined with other radio data (using a different flux density scale), or if the absolute flux density values are of interest, the model uncertainty has to be taken into account.

Median spectral indices
In order to characterize the TELAMON TeV-sample in terms of the spectral index, we calculate spectral indices for every observation and present their medians in Table 3.For some sources there is no sufficient frequency coverage available (i.e., less then three individual frequencies detected).Therefore, we could determine a sensible spectral index for 43 of our sample sources.

Spectral index calculation
As established earlier, there can be two measurements of a source per epoch and frequency, which means that there exists more than one flux density value for the same frequency.This can be seen in the example presented in Fig. 3 (magenta squares).In this case, the average flux density per frequency is calculated by taking the mean.If there is only one flux density value per frequency available, it is used as-is.The averaged values are depicted as black dots in Fig. 3.The uncertainty is determined by Gaussian error propagation, since distinct scans are considered independent measurements.After this first averaging process, a spectral power-law fit is performed to the data for every source at every observed epoch.This power-law is defined via where S is the source flux density, ν the observed frequency and α the spectral index.The fit is performed using a Levenberg-Marquardt fitting algorithm (Moré 1978) which ensures a robust least-squares optimization.Moreover, it allows for bounds on the fitting parameters, which is important for the sub-band averaging (cf.Sec.3.2.1).An example of such a fit for the source J1653+3945 (Mrk 501) on July 28, 2021 is shown in Fig. 3, which results in a spectral index α = −0.262±0.021.This analysis is applied to the entire data set, where at least three sub-band detections across all receivers (i.e., 14 GHz−45 GHz) are available for the same epoch and source. 3 In order to get a typical value for the spectral index for every source, we list their median spectral indices in Table 3, taken over all epochs.This is done only for sources for which spectral indices could be determined in three or more independent epochs to avoid outliers.

Spectral indices in the TELAMON sample
All median spectral indices across all bands (i.e., 14 GHz−45 GHz) are presented in Table 3, and in Fig. 4 as a histogram plot binned according to the source type.We find an overall average spectral index of −0.07, with a standard deviation of 0.20 and a median of −0.11.This is consistent with the expectation of a flat spectrum.We therefore conclude that the TeV-selected sub-sample of blazars presented here does not show any unexpected or special spectral features.In order to test if the spectral index distribution of HBL and EHBL is statistically different, we perform a two-sample Kolmogorov-Smirnov-test (KS-Test) which leads to a p-value of p = 0.22.This indicates that the distributions cannot be clearly distinguished and are consistent with having the same underlying statistics.For all other source classes, the sample size is not sufficient to perform significant statistical studies. 3In some cases (especially at 7 mm) a source was occasionally observed but not detected (e.g., due to poor weather conditions).In principle, this non-detection yields information on an upper limit to the flux density of the source which could be used to put additional constraints on the spectral index.However, due to oftentimes rapidly changing weather conditions (especially cloud coverage) the determination of such upper limits is not reliably reproducible.We discuss the detectability of sources further in Sec.4.3 where we derive flux density thresholds for which one can expect a significant detection for each receiver.Using these thresholds as estimates for the upper flux-density limits in case of non-detections, it becomes clear that they are too high to have any significant impact on the spectral index fit.

Average flux densities
In order to simplify the analysis, the observed flux densities of every source and epoch are averaged over the sub-frequencies of the 20 mm, 14 mm and 7 mm receivers, respectively.This means for every source we derive an average flux density value at 20 mm, 14 mm, and 7 mm for each epoch, given that the source was observed and detected in these bands.The main reason for this is that in some cases not all of the four (two for 20 mm) sub-bands of each receiver show a significant source detection.This can be due to RFI, because the source is too weak at the highest frequencies, or due to background noise.Averaging over the sub-bands of each receiver will make it possible to compare all epochs with each other, even if a sub-band flux density value might be vacant in one or more epochs.However, taking the mean or the weighted mean of all measured values does not suffice to ensure the comparability of epochs.For example, if a source was detected at 19.25 GHz and 21.15 GHz at Epoch A and at another Epoch B at 22.85 GHz and 24.75 GHz, taking the mean will shift the mean frequency of this average.Since the spectrum of the sample sources is not always flat (S (ν) const.), this could turn out to be problematic when comparing flux densities from two distinct epochs.We therefore introduce a new method of obtaining average flux densities from the receiver subbands using spectral fits in the following section.

Sub-band averaging
For each receiver (7 mm, 14 mm and 20 mm), at first a powerlaw spectral fit (cf.Eq. ( 7)) is performed to the raw data within the receiver band width.Similar to the spectral index calculation in Sect.3.1.1,a Levenberg-Marquart fitting algorithm is used to calculate the fit.Note that here, we set a bound of |α| < 0.5 to the spectral index, since the analysis presented here is more sensitive to outliers.From the fit, the average flux density S in the range of the receiver bandwidth is calculated: −0.29 a Average flux density and standard deviation at the given wavelength taken over all observed epochs.b Detection rate (i.e., number of epochs with detection divided by the number of epochs where the source was observed) at the given wavelength.c Number of epochs where the source was observed at the given wavelength.d Median spectral index across all observed epochs with at least three sub-band detections.e Source was observed but not detected at 20 mm and 14 mm.We have confirmed with independent 45 mm observations that the source is indeed very faint, i.e., a successful detection at wavelengths < 20 mm is unlikely.For these faint sources we have defined a new observing category and are now monitoring them with high detection rates at 45mm.f Affected by source confusion with a brighter neighbouring source, therefore no detection was possible.

Flux Density [Jy]
Best Fit 0506+056 min/max Fit Fig. 5: Example of the sub-band averaging process (J0509+0541, known as TXS 0506+056 on Jan 2, 2021 with two distinct scans at 14 mm).The average flux density is calculated by the integral of the best fit between 19 GHz and 25 GHz.The uncertainties are determined by integrating the min/max fits.
For the 20 mm receiver, the integration limits are ν 1 = 14 GHz and ν 2 = 17 GHz, for the 14 mm receiver ν 1 = 19 GHz and ν 2 = 25 GHz, and for the 7 mm receiver ν 1 = 36 GHz and ν 2 = 44 GHz.This method ensures that the average flux density is always calculated from the same frequency range, with best knowledge about the intrinsic spectrum of the source for each epoch.In the case where only data for one sub-frequency (e.g., only for 19.25 GHz) are available, we use the data as-is.This is equivalent to assuming a flat spectrum (α = 0) over the receiver bandwidth, which is typical for compact radio jets (e.g., Zensus 1997).
In order to define an uncertainty for the average flux density values, different approaches are used according to how many sub-bands are detected.In the case of two, three or four (i.e., all) sub-frequencies detected, two additional (uncertainty-)fits are performed: One with the uncertainties subtracted from the measured flux density values and the other one with the uncertainties added to the measured values (min-and max-fit).This is justified since in this case, the uncertainty mostly consists of systematic calibration uncertainty, i.e., the measured values may all together be higher (or lower) than the best values but relative to each other, they are known much more precisely.This procedure is illustrated in Fig. 5.The uncertainty is determined by the difference between the best fit to the minimum and maximum flux density (min/max-fit), and the best fit to the flux density value.In the case of only one sub-frequency detected, this is vastly different, since no sensible fit can be performed.The uncertainty of this value is estimated in a similar manner to the min/max-fits mentioned earlier.Again, the measurement uncertainties are subtracted/added to the flux density value.Then, different spectral distributions with α ∈ {−0.5; 0; 0.5}, originating at the min/max flux densities, are considered.This means, one gets three alternate spectral distributions for minimum and maximum flux density, respectively.Analogous to the previous calculations, the difference between the integrated alternate spectral distributions and the integrated best fit (in this case the flat spectrum) is calculated.The maximum difference is then used as the average flux density uncertainty.

mm
Fig. 6: Distribution of the average flux density found in the TELAMON sample sources for the 20 mm (top), 14 mm (center), and 7 mm (bottom) receivers.The histogram is divided into the source classes EHBL (blue with diagonal stripes), HBL (blue), IBL (yellow with dots), LBL (magenta with circles) and FSRQ (black with vertical stripes).The red vertical lines indicate the flux density above which one can expect a > 50 % detection chance per epoch (see Sect. 4.3).Below these thresholds the statistics are significantly affected due to the limited telescope sensitivity.For 20 mm we do not have enough data on the faintest sources to determine this limit.
In the following sections, when referring to flux densities at 7 mm, 14 mm or 20 mm, we always refer to the sub-bandaveraged flux density in each band, calculated using the method introduced in this section, if not declared otherwise.

Flux densities in the TELAMON sample
The sub-band averaging procedure introduced in the previous section provides us with one flux density value at 20 mm, 14 mm and 7 mm for every source at every epoch, given that it was observed and detected.In order to get an overview of the general source properties, we calculate the average flux densities of every source for 20 mm, 14 mm and 7 mm over time from the sub-band averaged flux densities.This allows us to characterize the observed sample in terms of average source flux density.All average flux density values are presented in Table 3. Their uncertainties correspond to the standard deviation of the sub-band averaged flux densities and therefore also reflect the intrinsic source variability.In the case where a source has only been detected once, the presented uncertainty is equal to the uncertainty of the sub-band averaged flux density.All average flux densities presented in Table 3 are also depicted as histogram plots in Fig. 6, binned according to the source type.Here, one can see that the 7 mm flux densities (Fig. 6, bottom) are by far the smallest sample and limited to flux densities ≳150 mJy.Only 25 out of 41 observed sources show a significant at 7 mm.This can be explained by the detection limit of the 7 mm receiver (see Sect.4.3) and by the fact that most of the sources in our sample are too faint to be detected at 7 mm.This is vastly different at 20 mm and 14 mm.At 14 mm (Fig. 6, center), 51 of the 55 observed sources show a significant detection in at least one epoch.At 20 mm (Fig. 6, top), 44 of the 47 observed sources show a significant detection in at least one epoch.Only three observed sources (J0152+0146, J0847+1133, J2250+3825) have not been detected at all.We have confirmed with independent 45 mm observations that these sources are amongst the faintest in our sample and therefore most likely below the detection limit of the telescope at shorter wavelengths with the current setup.
Following from the values presented in Table 3, we find an average 7 mm flux density of 1.4 Jy with a standard deviation of 1.8 Jy and a median of 0.6 Jy in our sample.This can only be considered as an upper limit to the average 7 mm flux density of the entire sample, since it is heavily biased by the non-detection of sources with flux densities ≲150 mJy due to the sensitivity limit at 7 mm.At 14 mm, the average source flux density is 0.7 Jy with a standard deviation of 1.3 Jy and a median of 0.2 Jy.At 20 mm, we find an average source flux density of 0.24 Jy with a standard deviation of 0.39 Jy and a median of 0.12 Jy.At the latter two wavelengths, the receiver sensitivity also limits the detection of very faint sources, but since almost all sources from our sample were detected in these bands, we consider this bias to be much less significant than for 7 mm.This allows us to perform a sensible statistical comparison of the 20 mm and 14 mm data.For both wavelengths, we find that according to a two sample KS-test one cannot reject the null-hypothesis that EHBL and HBL have the same underlying statistical distribution (14 mm: p = 0.02; 20 mm: p = 0.19).However, the sub-sample of EHBLs seems to populate lower flux densities than the sub-sample of HBLs according to the median average flux densities: At 20 mm, the subsample of EHBLs exhibits a median average flux density of 0.048 Jy (mean: 0.15 Jy, RMSD: 0.30 Jy), while for HBLs one finds a median of 0.13 Jy (mean: 0.16 Jy, RMSD: 0.13 Jy).At 14 mm, the subsample of EHBLs exhibits a median average flux density of 0.045 Jy (mean: 0.15 Jy, RMSD: 0.28 Jy), while for HBLs one finds a median of 0.15 Jy (mean: 0.27 Jy, RMSD: 0.31 Jy).The larger mean values are driven by one single EHBL (Mrk 501) with unusually high flux density.Considering instead the median flux densities, the impact of this single source is reduced and it appears that a majority of sources from the subsample of HBLs shows higher flux densities than the sub-sample of EHBLs, at moderate significance, which is in line with visual inspection of Fig. 6.This is consistent with the per definition higher synchrotron peak frequencies of EHBLs which results in a shift of their SED towards higher frequencies.Consequently, we find that the high-synchrotron peaked sources are fainter in the radio band than the lower-peaked sources in our sample.

Detection rates
In order to analyze the impact of the receiver sensitivity on the detectable source flux density, we calculate detection rates for all observed sources and receivers.For every source, we count the observing epochs (dates) during which the source was observed and the epochs where the source was detected in at least one subfrequency of the given receiver.The detection rate is then given by the number of detections divided by the number of observations.The detection rates and number of observations for each receiver and source are presented in Table 3 next to the average flux density values.These detection rates are further discussed in section 4.3.

Discussion
In this section, we discuss and compare the TELAMON source properties with previous programs, namely the F-GAMMA monitoring (Fuhrmann et al. 2016;Liodakis et al. 2017;Angelakis et al. 2019) and a study on TeV-selected blazars by Lindfors et al. (2016).

Spectral index discussion
Since the study by Lindfors et al. (2016) was carried out only at a single radio frequency (15 GHz), we can only compare spectral indices with the F-GAMMA program.Angelakis et al. (2019) provide spectral indices for various frequency windows.Their high-band (i.e., 14.6 GHz−43 GHz) is compatible with the TELAMON frequency range (i.e., 15 GHz−44 GHz).We therefore use the F-GAMMA high-band spectral indices provided by Angelakis et al. (2019) to compare the spectral indices of the TELAMON and F-GAMMA samples.Note that the F-GAMMA data (Angelakis et al. 2019) includes several (nonblazar) sources with a steep spectrum (α ≲ −0.8), which the TELAMON sample does not show.These are mostly calibrator sources which we exclude for further discussion.Figure 7 shows a direct comparison of the distribution of both samples in form of a histogram plot.The overall sample size of F-GAMMA (118 sources with spectral index) is superior to the TELAMON sample (43 sources with spectral index).According to a twosample KS-test, we find that the samples most likely have the same underlying distributions (p = 0.24).This suggests that GeV-emitting blazars and TeV-emitting blazars exhibit similar spectral indices.

Flux density discussion
Similar to the spectral index comparison, we have calculated average flux densities for the F-GAMMA data at 23.05 GHz from Angelakis et al. (2019) to compare with the TELAMON average flux densities at the overlapping wavelength of 14 mm.On top of that, we use the average flux density values provided by Lindfors et al. (2016) at 15 GHz for an additional comparison.All three flux density distributions are depicted in Fig. 8.One can clearly see that the flux densities of the GeV-selected F-GAMMA sources are in general much higher than for the TeVselected Lindfors et al. (2016) and TELAMON samples.This is consistent with the fact that most TeV emitting blazars have their SED shifted to higher energies and are consequently fainter in the radio band.We further performed a two-sample KS-test to compare TELAMON flux densities with the F-GAMMA and Lindfors et al. (2016) sample and find with high significance that the TELAMON and Lindfors et al. (2016)    a similar distribution (p = 0.30) while we have to reject the null-hypothesis of similar underlying distributions for TELA-MON and F-GAMMA (p = 1.4 × 10 −16 ).We therefore conclude that this difference in flux density seems to be characteristic for the target selection of either GeV-(F-GAMMA) or TeV-selected (TELAMON, Lindfors et al. 2016) blazars.TeV-emitting blazars seem to be fainter radio emitters than GeV-emitting blazars.Note that all three programs are limited by similar telescope sensitivity, which means we expect that this has no impact on the statistical comparison of the samples.

Long-term monitoring strategy
In order to address the question of the telescope sensitivity limit and to suggest a TeV-blazar sample well suited for long-term monitoring observations with the Effelsberg 100-m telescope, we discuss the relationship between source flux density and detection rate for all three receivers (20 mm, 14 mm and 7 mm) in this section.One can see in Table 3 that for some very faint sources the detection rate is rather low (i.e., ∼10 %).Observing these sources in a long-term monitoring study would lead to a waste of telescope time that could instead be used to include other (brighter) sources with better detection potential more frequently.On top of that, in order to perform variability and multiwavelength correlation studies, a sufficient sample size (flux density measurements per source) is required, which would take decades for sources with low detection rates and an observing cadence of two to four weeks.Moreover, such a sensitivity study can be useful to derive upper flux density limits in case of a nondetection in future studies.
In principle, the detection rate should only depend on the sensitivity of the receiver and can be calculated by telescope intrinsic parameters.However, given that we are performing a long-term monitoring program, bad weather epochs and intrinsic source variability can influence the detectability of fainter sources, especially at the highest frequencies.To maximize source detections, we try to prioritize the faintest sources during good weather conditions, but this is not always possible since we are trying to keep a cadence of two to four weeks for all sources.In principle, the detection chances for the faintest sources can be increased by using more scans, i.e., spending more observing time per source.However, due to the limited observing time and large number of sources this was usually not feasible.We consider the 2.5 year observations presented in this paper to be a good representation of such effects and therefore well suited to identify a flux density limit above which it is sensible to monitor sources as part of a long-term study with the given telescope and observing setup.For each receiver, we calculate the average detection rate for all sources below a certain flux density threshold.This is done for 10,000 different flux density thresholds.Moreover, this entire procedure is carried out 1,000 times with varying average source flux densities in a Monte-Carlo (MC) way, assuming the average flux densities and their uncertainties in Table 3 correspond to Gaussian distributions.The mean average detection rate and its standard deviation taken from the 1,000 MC-iterations is then plotted in Fig. 9 for 10,000 different flux density threshold values.In order to guarantee statistical convergence, we require a minimum of three sources to be included in the average detection rate count, therefore the average detection rates in Fig. 9 can only be presented above a certain flux density threshold.If we assume that the detection rate is monotonically increasing with flux density, the detection rate curves in Fig. 9 correspond to lower limits of the actual detection rate at any given flux density.For a sensible long term monitoring, we require a minimum detection rate of 50 % (i.e., a source detection at least every other epoch) to optimize the scientific outcome of the observing time and to prevent too many non-detections.Using the data from Fig. 9, we can determine a flux density limit S 0 for each receiver above which we expect a detection rate of at least 50 %.For the 14 mm receiver, we find S 0 = 70.8+8.3 −6.3 mJy, and for 7 mm S 0 = 337 +36 −50 mJy.For the 20 mm receiver, we do not have enough data on the faintest sources to determine a flux density threshold for a detection rate of 50 %.All sources with a 20 mm flux density ≳ 30 mJy have a detection chance > 50 %.
Sources below a 20 mm flux density of ∼30 mJy have a limited detection chance with either of the three introduced receivers.We therefore started observing these sources at even longer wavelengths (i.e., 45 mm, sub-sample I) to have the best chances to detect them in our long-term monitoring program.The second faintest sources, exhibiting a 20 mm flux density ≳ 30 mJy but a 14 mm flux density ≲ 70 mJy, are monitored at 20 mm (sub-sample II) in the current observing setup.Brighter sources (S 14mm ≳ 70 mJy and S 7mm ≲ 350 mJy) can be moni-tored at 20 mm and 14 mm (sub-sample III) in a sensible way, and for the brightest sources (S 7mm ≳ 350 mJy) sensible monitoring is possible at both 14 mm and 7 mm (sub-sample IV).Note that in principle these thresholds can be improved by spending more time on each source, i.e., performing more sub-scans.However, this tactic is not favorable for our program, since we need to establish an ideal trade-off between observing time spent per source and the total number of sources that can be monitored.In Fig. 10, we show example light curves of one representative source per sub-sample.One can see in light curves (a), (b), and (c), that we tried monitoring the sources with higher frequency receivers than suggested by the sub-sample, which resulted in fairly low detection rates.Since the adaption of the sub-samples, the detection rates have significantly improved.Ongoing TELA-MON observation are conducted according to this strategy as it is the best possible setup for our purposes.

Conclusions & outlook
We have presented the first results of the pilot-phase (August 2020 -February 2023) of the TELAMON AGN monitoring program for a complete sample of all TeV-emitting blazars in the Northern Hemisphere.We have used the Effelsberg 100-m telescope to monitor these sources at high radio frequencies from 14 GHz−44 GHz every two to four weeks.We have developed a semi-automated data reduction pipeline which allows for immediate reduction of the data.From the observations, we derived flux densities and spectral indices from the sources in our sample.Plots of the latest available light curves and source spectra are publicly available on the dedicated TELAMON website4 .
Our results are consistent with the findings of a prior study of the same object class by Lindfors et al. (2016).In comparison to the GeV-selected F-GAMMA sample (Fuhrmann et al. 2016;Angelakis et al. 2019), we find that the spectral indices of both GeV-and TeV-selected sources are consistent with a flat spectrum.The average flux density at 14 mm is significantly lower in the TeV-selected sample than what Angelakis et al. (2019) have found in their study.The comparison between GeVand TeV-emitting blazars is still limited by the small number of TeV-emitting blazars, which is expected to increase within the next years with upcoming observations from the Cherenkov Telescope Array (Cherenkov Telescope Array Consortium et al. 2019).
In future studies, we will analyze the variability of TeVblazars and perform correlation studies with multi-wavelength light curves, similar to previous studies of GeV-emitting AGN (e.g., Paraschos et al. 2023;Rösch et al. 2022;Fuhrmann et al. 2014).For several individual sources we have already demonstrated the multi-wavelength capabilities of our program (e.g., Gokus et al. 2022;Menezes et al. 2021), especially in combination with our complementary Southern Hemisphere monitoring program at the Australia Telescope Compact Array (e.g., Eppel et al. 2023;Satalecka et al. 2021).Moreover, the TELAMON program will continue its observations using an optimized observing setup including measurements at 45 mm for very faint radio sources.In addition, we will provide results on follow-up observations of neutrino-candidate blazars, which have not been discussed in this paper.Another publication reporting on the polarization properties of our TeV-sample is currently in preparation (cf.Heßdörfer et al. 2023)

Fig. 1 :
Fig. 1: An example of the calibration factor evolution at 38.75 GHz during the observing session on Oct 24, 2021.The black dots correspond to the (sub-)averaged calibration factor measurements Γ c,i .The blue line and blue squares indicate the simple moving average (SMA) interpolated calibration factor Γ SMA,i (see Equation (1)), as used in the final analysis.The dark gray area represents the calibrator scatter σ sc,ν , while the light gray hashed area represents the total calibration uncertainty σ cal , including the model uncertainty σ model and σ sc,ν .

Fig. 4 :
Fig.4: Distribution of the median spectral index across all available bands (i.e., 14 GHz−45 GHz), as found in the TELAMON sample sources.The histogram is divided into the source classes EHBL (blue with diagonal stripes), HBL (blue), IBL (yellow with dots), LBL (magenta with circles) and FSRQ (black with vertical stripes).

Fig. 7 :
Fig. 7: Histogram plot of median radio spectral indices for all sources from the F-GAMMA sample (Angelakis et al. 2019, blue with diagonal stripes) and the TELAMON (black with vertical stripes) sample between 14 GHz and 45 GHz.

Fig. 9 :
Fig.9: Average detection rate for sources below a given flux density threshold for the 20 mm, 14 mm and 7 mm receiver as calculated from the first 2.5 years of TELAMON observations.The errors are calculated in a Monte-Carlo way, which takes into account the variability of the observed sources.
Fig. 10: Illustrative light curves (upper panels) of four selected sources from the TELAMON TeV-sample.The flux densities are averaged over the individual receiver sub-bands as described in Sect.3.2.1.The lower panels indicate the times of observation with the given receiver.If there is no matching flux density at the same time in the light curve, the source was not detected.

Table 1 :
The TELAMON sample of TeV-emitting AGN.The classification into EHBL, HBL, IBL and LBL is performed according to the synchrotron peak frequency found in the literature.Synchrotron peak frequency in the source rest frame, i.e., ν peak,rest = ν peak,obs (1 + z).Already corrected values from the literature were adjusted to the redshift values in this work. a )

Table 3 :
Average flux densities for all sources from the TELAMON sample as observed during the first 2.5 years.The median spectral index has been calculated across all receiver bands and epochs.rates and the number of epochs each source was observed are shown. .