© The Authors 2025

1 Introduction

Blazars are among the most enigmatic and intriguing objects in the Universe. They are a subclass of active galactic nuclei (AGNs) and are characterised by an intense and highly variable emission across the whole electromagnetic spectrum. According to a well-established picture, at their cores lies a supermassive black hole surrounded by an accretion disk of swirling gas and dust. The black hole powers a relativistic jet that is closely aligned with the observer’s line of sight (Blandford & Königl 1979). Because of the bulk relativistic motion of the emitting plasma towards the observer, the apparent luminosity of blazars is affected by relativistic beaming, which can dramatically increase the amplitude and reduce the timescale of the observed flux density variations. Doppler boosting also results in an apparent superluminal motion of the plasma itself. The blazar class is traditionally divided into two subclasses: flat-spectrum radio quasars (FSRQs) and BL Lacertae objects (BL Lacs). While both exhibit intense variability and possess relativistic jets, they differ in their spectral properties and the strength of their emission lines.

Variability studies play a pivotal role in unravelling the mysteries of blazars. Blazars emit radiation across a broad range of wavelengths, from radio to gamma rays, with the radio regime offering unique insights into the nature and behaviour of these sources. The extreme variability of blazars challenges our understanding of the underlying mechanisms driving the intense emissions from their cores. Multi-wavelength variability studies, such those carried out by the Whole Earth Blazar Collaboration¹ (WEBT), can provide invaluable information about the physical properties of the blazars’ emitting region (see e.g. Villata et al. 2002, 2006; Raiteri et al. 2017, 2024). Founded in 1997, the WEBT is an international collaboration of astronomers whose aim is to monitor blazars in chiefly optical bands, but also the radio and near-infrared, to study blazar variability on various timescales. The insights provided by the variability characterisation are essential for discriminating among different emission models, estimating fundamental physical parameters (such as the size of the emitting regions and Doppler boosting factors; see e.g. Liodakis et al. 2017a), and casting light onto the processes that power this class of AGN. Since radiation at different wavelengths is emitted in different parts of the jet, multi-frequency observations provide us with a virtual view of the structure of the jet on different scales. Through the analysis of time delays between variations occurring at different wavelengths, it is possible to infer how perturbations propagate through the emitting regions located in different parts of the jet.

In the realm of radio observations, blazar variability is particularly intense, with fluctuations in radio flux densities occurring on timescales ranging from days to decades. Studying the time variability of blazars in the radio band sheds light on the dynamic processes occurring in their jets, providing clues about their formation, acceleration, and collimation. As confirmed by the Fermi-GST AGN Multi-frequency Monitoring Alliance (F-GAMMA) project, conducted at the Effelsberg telescope in Germany (Fuhrmann et al. 2016; Angelakis et al. 2019), long-term radio monitoring programmes with high spectral coverage are of the utmost importance for disentangling the wide variety of variability manifestations observed in blazars. The study of samples with complementary selection criteria, such as in the case of the TeV Effelsberg Long-term AGN Monitoring (TELAMON) programme (Eppel et al. 2024) can increase the richness of the input data. However, this goal can only be achieved by perpetuating monitoring programmes for several decades, as recently shown by Kankkunen et al. (2025).

In 2004, a programme for the monitoring of blazars, including both FSRQs and BL Lacs, started at the 32 m radio antennas of the INAF-Institute of Radioastronomy (IRA) in Medicina and Noto (see Bach et al. 2007). A short-term blazar monitoring programme had previously been run with the same telescopes from 1996 to 1999 (Venturi et al. 2001). Observations were initially performed at frequencies of 5, 8, 22, and 43 GHz with a monthly cadence. Despite some discontinuities in frequency and time coverage (the 22 GHz observations were soon replaced by observations at 24 GHz; 5 and 43 GHz observations have been dropped because of radio-frequency interference and problems with the receiver, respectively; and periods of prolonged extraordinary maintenance of the antennas have created some gaps in the light curves), the ROBIN² (Radio Observations of Blazars with INaf telescopes) programme has been perpetuated across the years and is still operational, having reached 20 years of activity. The data collected within the programme have been used in several studies and publications, mostly in correspondence with WEBT multi-frequency campaigns (see e.g. Otero-Santos et al. 2025; Raiteri et al. 2024; Villata et al. 2006; Raiteri et al. 2007, 2021, 2009; Larionov et al. 2008; Carnerero et al. 2015) and in collaboration with the gamma-ray community (MAGIC Collaboration 2024; Abe et al. 2023; Abdo et al. 2010; Vercellone et al. 2010; Hayashida et al. 2012; Aleksic et al. 2015).

In this paper, we offer an overview of the ROBIN programme, providing the reader with basic information about the monitored sources and the observations. The collected flux density measurements within the time interval that goes from 2004 to 2024 are available at the CDS; the same data, along with more recent ones, are also stored on an online archive³, which is accessible upon request. We want to stress that an in-depth analysis of the variability of the sources is beyond the scope of this article, which aims instead to properly present and describe this important database. A thorough time series analysis of the light curves will be presented in a subsequent publication.

The paper is organised as follows: The source sample and the observations are briefly discussed in Sects. 2 and 3. The data reduction is described in Sect. 4, together with the results of a basic analysis of the data. In Sect. 5 we present our conclusions.

2 The sources

The database of the ROBIN programme comprises flux density measurements of 47 sources, listed in Table A.1; following the classification from Massaro et al. (2015), the sample includes 16 BL Lacs and 27 FSRQs, while for 4 objects no certain classification is given. For 0859+0454, the FSRQ classification is derived from Healey et al. (2007). Currently, all the listed sources are still monitored, except 0323+342, which has not been observed since the end of 2018, because its variability is comparable to the uncertainty in our flux density measurements, and therefore, it is impossible to derive reliable variability characteristics for it.

The sources were selected according to a few broad criteria: a high degree of variability in all the covered radio bands; observability from Medicina and Noto at reasonably high elevations; and flux densities high enough to enable sufficiently precise measurements for a variability study. Some basic information about the sources is provided in Table A.1; the redshift values have been extracted from the NASA/IPAC Extragalactic Database⁴ (NED).

3 Observations

Observations were almost entirely performed using the 32 m dishes located in Medicina and Noto, which are both owned and managed by the National Institute for Astrophysics (INAF; Italy). In the first few months of the project, observations at 5, 8, and 22 GHz were obtained at both sites. Looking at the performances of the antennas, it was successively decided to split the observations, using the Medicina antenna for the monitoring at the higher frequencies, and Noto for the 5 GHz observations, since at this frequency the Medicina antenna was much more affected by radio-frequency interference (RFI). In early 2006, a new 43 GHz receiver became available in Noto (see Leto et al. 2009 for detailed information about observations performed using this receiver), allowing us to extend our monitoring to higher frequencies, where blazar variability is generally stronger. One epoch of observation was performed at the Sardinia Radio Telescope, a 64 m dish equally managed by INAF, when both 32 m dishes were unavailable.

Radio continuum acquisitions were carried out by means of on-the-fly cross-scans in equatorial coordinates. Table 1 lists the telescope main features and employed configurations, while Table 2 provides the scanning parameters used during the observations. The number of cross-scans performed on a given source varied according to its expected flux density. Typical on-source integration times were 37.5 s at 5 GHz, 40.0 s at 8 GHz, 52.5 s at 24 GHz, and 20.0 s for 43 GHz acquisitions. Occasionally, some sources were observed in the same epoch and at the same frequency more than once. For the purposes of the present work, the flux density measurements from redundant observations have been averaged to provide a single measurement per epoch per frequency. Overall, the ROBIN database comprises 20 915 data points, which translates, after the averaging, into the 12 675 measurements that are presented in this paper.

Flux density calibration was carried out observing 3C 123, 3C 286, 3C 295, and NGC 7027, whose reference flux densities were computed, for the observed band central frequency, according to Perley & Butler (2013). For 24 GHz observations, the atmospheric contribution was also taken into account in the calibration procedure; zenithal opacity was estimated by means of skydip acquisitions and then applied to calibrate the signal amplitude.

4 Data reduction and analysis

Our group developed an IDL pipeline - the cross-scan analysis pipeline (CAP) - to reduce and analyse the cross-scans. It is fully described in Giroletti & Righini (2020). The pipeline performs the following operations:

• Data rearrangement: The dataset is organised by frequency and the purpose of the observed sources (calibrators, skydips, and targets).

• Data flagging: Users can optionally use a graphical user interface (GUI) to inspect each sub-scan and assign flags, determining which data to include in the analysis. Since FITS files include both left circular polarisation (LCP) and right circular polarisation (RCP) data streams, which can be impacted differently by interferences (RFI), users can choose to include or exclude specific polarisations.

• Skydip reduction: For frequencies above 10 GHz, skydip acquisitions are used to estimate atmospheric opacity. These are fitted using atmospheric models, telescope parameters (e.g. receiver temperature), and weather data from the acquisition period, producing a table of zenithal opacity estimates.

• Calibrator reduction: Integrating acquisitions from flux density calibration sources involves aligning and integrating flagged sub-scans, correcting for antenna gain curve and atmospheric opacity as needed. Each sub-scan is handled individually due to elevation dependence. A Gaussian fitting is performed to measure signal amplitude in raw counts. The measured amplitude is corrected for pointing errors, and the theoretical flux density of the source is calculated based on the observed band, providing a counts-to-jansky conversion factor.

• Conversion timeline: After processing data from flux density calibrators, users must decide how to interpolate the counts-to-jansky conversion factors. These factors, ideally constant after correcting for antenna and weather effects, can still fluctuate due to weather or instrument instabilities. Users can choose to average or linearly interpolate the conversion factors within a chosen time window.

• Target reduction: This process mirrors the ‘calibrators reduction’ phase, where the raw amplitude from Gaussian fitting is multiplied by the appropriate calibration factor from the counts-to-jansky timeline. This results in a table of measured flux densities for each scan. LCP and RCP data streams are processed separately, producing independent flux density measurements. For each polarisation, the RA and Dec integrated semi-scans are averaged, yielding a single flux density. Partial results are documented in a separate table, and users can inspect the results of each original sub-scan if desired.

It must be noted that, as the operating and acquisition systems at the telescopes were still being developed during the many years of observations, the FITS file content details varied along time. CAP was conceived to manage this range of input formats.

4.1 System performance and statistics

At 24 GHz, even though the reduction pipeline is able to compensate flux density measurements for the pointing error, observations included pre-scan corrections of the antenna pointing. This was achieved by means of dedicated cross-scans on the targets themselves - when sufficiently bright - or on close bright sources. This allowed the subsequent acquisitions to be carried out with a pointing error mostly smaller than one-tenth of the half-power beam width (HPBW). Errors larger than about 0.3 times HPBW were usually associated with bad weather conditions, also leading to poor quality data, so the associated acquisitions were discarded from the analysis. Individual measurements, especially in X band - the cleaner frequency band, less affected by RFI and weather effects - show flux density errors as low as 1%. To accommodate for flagging inaccuracies and fluctuations caused by the weather or by residual errors in the calibration, we applied a conservative minimum value to the flux density errors: 3% for C- and X-band measurements, 8% for K-band flux densities, and 10-15% for Q-band flux densities. Refined results are expected to be provided by a future version of the CAP pipeline, involving artificial intelligence procedures.

An example of the multi-frequency light curves collected within the ROBIN programme is shown in Fig. 1. The monthly cadence provides for adequate coverage of the fast variability of the sources: both the rising and falling part of the flares can be reconstructed with reasonable accuracy. The full set of the multi-frequency light curves is shown in Appendix B.

Fig. 1 Light curves of 2251+158 (i.e. 3C454.3) from the ROBIN programme. The 5, 8, 24, and 43 GHz data are plotted as blue, orange, green, and magenta dots, respectively.

4.2 Basic variability characteristics

Some simple analysis techniques have been applied to the data collected within the ROBIN programme to extract basic information concerning their variability characteristics.

An estimation of the amplitude of the variability has been made by means of the intrinsic modulation index algorithm developed by Richards et al. (2011): $$\overline{m}=\frac{{\sigma }_0}{\overline{S_0}},$$(1)

where σ₀ is the intrinsic standard deviation of the distribution of source flux densities in time, measured in units of the intrinsic source mean flux density, $\overline{S_0}$.

While the standard definition of the modulation index expresses the variability in terms of the standard deviation of the flux densities of the source, σ, and the average flux density, S̄, the intrinsic modulation index uses the flux densities and variations as would be observed with uniform sampling of adequate cadence and zero observational error. To estimate these quantities, Richards et al. (2011) developed a likelihood analysis starting from the assumption that the true flux densities are normally distributed, and therefore both $\overline{S_0}$ and σ₀ are well defined, and can be inferred from the available data.

Despite the improvements introduced by the intrinsic modulation index with respect to the standard definition, some problems still persist. The assumption of stationarity, implicit in the definition of a true flux and a well-defined standard deviation, cannot be regarded as generally valid in the case of the radio emission of blazars. The radio variability of blazars often appears as a superposition of several components with different timescales: the longest timescales can be as long as decades (see e.g. Hughes et al. 1992, and, more recently, Kankkunen et al. 2025). This means that, despite the remarkable time span covered by the ROBIN programme, the light curves obtained within the project are still not long enough to assume that $\overline{S_0}$ and σ₀ are always well defined. For this reason, the intrinsic modulation index values, which depend on $\overline{S_0}$ and σ₀, are sometimes almost as sensitive to the duration of the observations as the standard modulation index values, showing a tendency to increase with the length of the light curve.

Since the duration of the observations varies with the source and the observing frequency, the estimated m̄ values, reported in Table A.2, for the different bands, may not always be comparable to each other. We obtained an alternative estimation of the variability amplitude, not directly related to the temporal extension of the observations, via a first-order structure function (Simonetti et al. 1985; henceforth, simply ‘structure function’) that we calculated at a time lag of ≈1.5 years, $$S{F}_{1.5}=\frac{1}{N}\sum _{i,j}(S(t_i)-S(t_j){)}^2,$$(2)

where the sum extends to all N pairs (t_i, t_j) for which 1.4 yr < (t_i - t_j) < 1.6 yr.

S F_1.5 corresponds, roughly speaking, to twice the variance of the signal estimated using pairs of data with 1.5 years of separation. Starting from this consideration, we defined a new function, S F′_1.5, as follows: $$S{F}_{1.5}^{\mathrm{\prime }}=\frac{\sqrt{S{F}_{1.5}/2}}{\overline{S}}.$$(3)

Since this parameter, which is adimensional as the modulation index, provides a measure of the variability over a fixed interval of time, it is not affected by the dependence on the duration of the observations. The estimated S F′_1.5 values for the different bands are reported in Table A.2.

The time interval of 1.5 years we set for the SF′ calculation was defined looking at Hovatta et al. (2007); by analysing decades-long light curves of a large sample of AGNs at gigahertz frequencies, the authors found characteristics timescales (which, in the case of multiple variability components, can be assumed to be the dominant ones) mostly between 1 and 2 years (see also Hughes et al. 1992). A time interval of 1.5 years should therefore be sufficient to sample a large fraction of the dominant variability component for most of the sources. It should be noted that S F′_1.5 does not take into account the uncertainty in the measurements.

In Table A.2, we also report the values of S F′_1.5/SF′_3.0, where S F′_3.0 is calculated by summing the contribution of the pairs (t_i, t_j) for which 2.9 yr < (t_i - t_j) < 3.1 yr. Since the structure function provides a way to evaluate the variability amplitude at different timescales, the ratio S F′_1.5/S F′_3.0 assesses how strong the variability on a short timescale is compared to the variability on longer timescales. In this sense, it can be considered as a rough indicator of the speed of the flux density variations.

The average spectral index between 8 and 24 GHz, α (Table A.2, Col. 14), were calculated over all the observing sessions for which both 8 and 24 GHz data are available. As a convention, we used S (ν) = S′ν^α, which means that positive α values denote an increase in flux density with frequency.

An example of the spectral indices calculated for each epoch is shown in Fig. 2 for BL Lacertae. Given the large gap in the 24 GHz, the data have been complemented with the 20 GHz observations from the F-GAMMA programme (see Fuhrmann et al. 2016 and Angelakis et al. 2019), re-scaled by a factor of 1.05, which accounts for the slight difference in the flux densities measured at non-coincident frequencies and has been obtained by calculating the average ratio between quasi-simultaneous (i.e. separated by no more than 10 days) ROBIN and F-GAMMA flux density measurements. The spectral indices shown in the lower panel have been calculated after merging the ROBIN and the F-GAMMA 20-24 GHz data: the clear trend of an increasing spectral change between 24 and 8 GHz data is confirmed by the trend of α, which may have started before MJD 58000 and seems about to reach a plateau.

Fig. 2 Upper panel: ROBIN light curves of BL Lacertae at 24 and 8 GHz (light green and orange dots, respectively) shown together with the re-scaled (see the main text for the details) 20 GHz ones from the F-GAMMA programme (dark green dots). Lower panel: spectral indices calculated using the combined ROBIN and F-GAMMA data at 20-24 GHz and the ROBIN data at 8 GHz. The clear increasing trend of the spectral index (magenta line) reflects the increasing discrepancy between the light curves in the last 2000 days of observations.

4.3 Comparison between variability amplitude estimators

A strong correlation is expected between m̄ and S F′_1.5. The results of our analysis, shown in Fig. 3 and Table 3, confirm this expectation: the correlation coefficients calculated for the four observed bands range between 0.74 and 0.91.

It should be noted that the 8 GHz light curve of 2344+514 has not been taken into account for the comparison between m̄ and S F′_1.5, because, according to m̄, the source, at this frequency, does not show variability. Most likely the evaluation of this dataset as non-variable (which might be due to a combination of poor sampling and high uncertainties in the flux density measurements corresponding to possible flares) is not correct, as the variability of the source at 8 GHz has a timescale of months and is not consistent with white noise; the displayed variations, moreover, are correlated to the ones observed at 24 GHz, which are certainly source intrinsic.

Intuitively, we expect the correlation coefficient r and the slope B of the linear regression between m̄ and S F′_1.5 to increase with frequency; in fact, the higher the frequency, the stronger and faster should be the intrinsic variability of the source. As a consequence, the m̄ and the S F′_1.5 values should tend to converge, and B and r should both tend to 1. In this regard, the behaviour of the variability parameters at 24 GHz is puzzling. The slope B is considerably lower than at other frequencies, while the correlation coefficient aligns with the values estimated at 5 and 8 GHz, very far from the almost perfect correlation at 43 GHz. This behaviour seems to indicate that either S F′_1.5 overestimates the variability (possibly because of a higher level of noise in the data at this frequency) or m̄ underestimates it (this could be the consequence of an overestimation of the uncertainty in the data, which reduces estimated intrinsic variability).

The high correlation between m̄ and S F′_1.5 at all frequencies shows that both tools are capable of characterising the variability amplitude of the light curves reasonably well. However, as we mention in Sect. 4.2, they both suffer from some bias. In the case of S F′_1.5, it is important to take into account that this estimator does not take into account the observational uncertainties, which causes an overestimation of the variability amplitude. Concerning m̄, it does not distinguish variations according to the timescale on which they occur, resulting in a substantial dependence on the duration of the observations. By combining the two methods, however, we can reach a better understanding of the variability properties of the light curves.

Fig. 3 Comparison of the variability amplitude parameters at 5 (blue dots), 8 (orange squares), 24 (green triangles), and 43 GHz (magenta diamonds). The linear regression results are plotted as lines with the same colours.

Fig. 4 Minimum value of the spectral index calculated over all the available epochs versus redshift. Both a linear regression of the data (green line) and the average calculated over increasingly larger bins (in order to keep the number of contributing points comparable; orange diamonds) indicate a moderate anti-correlation between the two parameters.

4.4 Other parameters

Unsurprisingly, there is a strong correlation between the spectral index of the sources and their variability amplitude; higher (i.e. more inverted) spectral indices correspond to more variable sources, independently of the estimator used for the variability amplitude and the frequency at which this is quantified. This is due to the fact that more inverted spectra are supposed to be associated with more compact emitting regions, in which flux density variations are more extreme. On average, the correlation of α with S F′_1.5 provides higher correlation coefficients than the one between α and m̄, but the difference is moderate (on average 0.51 for S F′_1.5 and 0.45 for m̄). The highest correlation, 0.56, is found with S F′_1.5(24 GHz), while the correlation with m̄(24 GHz) is only 0.44. This is an interesting indication with respect to the issue discussed in Sect. 4.3 concerning the linear regression coefficients B and r at 24 GHz. The higher correlation between α and S F′_1.5 suggests that, in this case, the latter provides us with a more faithful evaluation of the variability properties compared to m̄.

No correlation is found between m̄ and S F′_1.5/S F′_3.0 any observing frequency, which suggests that higher-variability amplitudes are not systematically associated with an increase or a decrease in the variability timescale. No obvious correlation is seen between α and S F′_1.5/S F′_3.0 either, which, again, can be interpreted as evidence that the spectral index is independent of the speed of the variations.

All the parameters reported in Table A.2 have been checked for possible correlations with redshift. None was found; however, the minimum values of the spectral index calculated over all epochs show a weak anti-correlation with redshift at a level of −0.30 (see Fig. 4), which implies a probability of a chance correlation of ~4%. In particular, while at low redshift the minimum spectral index can go from very steep (−1.19) to almost flat (−0.21), with an average value of about −0.5, at high redshift it does not reach higher values than −0.56. Given the low number of data points, this result should be considered cautiously.

4.5 BL Lacs versus FSRQs

The mean variability parameters for BL Lacs and FSRQs are plotted separately on Fig. 5. Given the relatively low number of objects in the sample and the large spread in their characteristics, the uncertainties in the mean values are large, and the results discussed below should be considered as purely indicative.

Both m̄ and S F′_1.5 show that the variability of our sample of BL Lacs is higher than that of FSRQs; the difference, however, is remarkable only at 43 GHz. This result is in agreement with the one documented by Liodakis et al. (2017b); a higher variability for FSRQs, concerning the first 2.5 years of monitoring of the F-GAMMA programme, is reported by Fuhrmann et al. (2016), which though estimates the variability amplitude in terms of flux density variance instead of modulation index.

The light curves show a general increase in the variability with frequency, as expected from previous studies (see e.g. Fuhrmann et al. 2016), partially contradicted by the sudden drop of m̄ (43 GHz) for FSRQs. This drop is not seen in S F′_1.5(43 GHz), which may be surprising given the high level of correlations between the two parameters at this frequency (see the discussion in Sect. 4.3). The difference between the two is illustrated by the high value of the intercept A (see Table 3), which reflects a systematically higher value of S F′_1.5 with respect to m̄; this seems a further confirmation of the tendency of m̄ to underestimate the intrinsic variability of the sources.

For the mean values of S F′_1.5/S F′_3.0 we see systematically higher values for BL Lacs than for FSRQs (with a tendency to increase with frequency), suggesting that the former show generally faster variability than the latter and that variability timescales decrease with frequency. Note that values higher than 1, as obtained at 43 GHz, do not have an obvious physical interpretation, as S F′_1.5/S F′_3.0 should oscillate around 1 when the variability timescale is lower than 1.5 years, unless the light curves display quasi-periodic oscillations. The behaviour here observed is certainly due to the large uncertainty in the data. Nevertheless the increasing trend of speed with frequency is consistent with our expectations.

We evaluated the presence of statistical difference in the variability parameters of BL Lacs and FSRQs through a Kolmogorov-Smirnov test; no significant difference was found.

Fig. 5 Mean variability parameters for BL Lacs and FSRQs plotted separately as a function of frequency. Upper panel: m̄ and S F′1.5 values. Lower panel: speed of the variability calculated through S F′1.5/S F′3.0.

5 Summary

We have presented the data collected as part of the ROBIN programme. Launched in 2004, the project has so far gathered about 21 000 flux density measurements at frequencies of 5, 8, 24, and 43 GHz for a sample of 47 blazars. The light curves of all the sources in the sample, which comprises 27 FSRQs and 16 BL Lacs, plus 4 objects with uncertain classification, are shown in Appendix B.

We applied some basic analysis tools to the data to estimate the variability amplitude (through both the intrinsic modulation index and a structure-function-based algorithm), the speed of the variability, and the spectral index between 8 and 24 GHz. A linear regression was applied to the estimated parameters in a search for correlations among them. A significant correlation has been found only between the spectral index and variability amplitude of the sources, independently of the way the latter is estimated: the most variable sources are those with the most inverted spectra.

No difference is found between the variability parameters of BL Lacs and FSRQs on a significant statistical level, estimated through a Kolmogorov-Smirnov test. In general, the former are more variable at all wavelengths than the latter; also, their variability occurs on faster timescales and their spectral index is higher, although the differences are not significant.

After more than 20 years of activity, the ROBIN programme has been renewed for five more years, which will strengthen an archive that is among the most longstanding in this field of research. Up-to-date light curves of single sources from the ROBIN programme are available from the corresponding author upon reasonable request.

Data availability

The flux density measurements collected within the ROBIN programme between 2004 and 2024 are available at the CDS via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/702/A140.

Acknowledgements

The Medicina and Noto radio telescopes are funded by the Ministry of University and Research (MUR) and are operated as National Facilities by the National Institute for Astrophysics (INAF). The Sardinia Radio Telescope is funded by the Ministry of University and Research (MUR), Italian Space Agency (ASI), and the Autonomous Region of Sardinia (RAS) and is operated as National Facility by the National Institute for Astrophysics (INAF). We acknowledge financial support from the INAF Fundamental Research Funding Call 2023 (PI: Raiteri).

Appendix A Tables

Appendix B Light curves

Here we present the light curves of the 47 sources monitored within the ROBIN programme.

Fig. B.1 Multi-frequency light curves for all the sources monitored within the ROBIN programme. In all plots, blue dots represent 5 GHz data, while orange, green, and magenta dots respectively represent 8, 24, and 43 GHz data.

References

All Tables

All Figures

	Fig. 1 Light curves of 2251+158 (i.e. 3C454.3) from the ROBIN programme. The 5, 8, 24, and 43 GHz data are plotted as blue, orange, green, and magenta dots, respectively.
In the text

Fig. 2 Upper panel: ROBIN light curves of BL Lacertae at 24 and 8 GHz (light green and orange dots, respectively) shown together with the re-scaled (see the main text for the details) 20 GHz ones from the F-GAMMA programme (dark green dots). Lower panel: spectral indices calculated using the combined ROBIN and F-GAMMA data at 20-24 GHz and the ROBIN data at 8 GHz. The clear increasing trend of the spectral index (magenta line) reflects the increasing discrepancy between the light curves in the last 2000 days of observations.

In the text

	Fig. 3 Comparison of the variability amplitude parameters at 5 (blue dots), 8 (orange squares), 24 (green triangles), and 43 GHz (magenta diamonds). The linear regression results are plotted as lines with the same colours.
In the text

	Fig. 4 Minimum value of the spectral index calculated over all the available epochs versus redshift. Both a linear regression of the data (green line) and the average calculated over increasingly larger bins (in order to keep the number of contributing points comparable; orange diamonds) indicate a moderate anti-correlation between the two parameters.
In the text

	Fig. 5 Mean variability parameters for BL Lacs and FSRQs plotted separately as a function of frequency. Upper panel: m̄ and S F′1.5 values. Lower panel: speed of the variability calculated through S F′1.5/S F′3.0.
In the text

	Fig. B.1 Multi-frequency light curves for all the sources monitored within the ROBIN programme. In all plots, blue dots represent 5 GHz data, while orange, green, and magenta dots respectively represent 8, 24, and 43 GHz data.
In the text