High–Redshift Quasars at z ≥ 3: Radio Variability and MPS/GPS Candidates

1. Introduction

Quasars [1] are highly variable celestial objects [2] that emit radiation across the electromagnetic spectrum. In the radio band, their variability [3,4] is often attributed to the presence of outflows or jets that emanate from the central region of a galaxy [5]. These jets can exhibit a range of behaviors, including propagating shocks (e.g., [6,7]), precession (e.g., [8,9]), bending (e.g., [10]), and interaction with gas clouds (e.g., [11]). In addition, several other phenomena may contribute to the observed flux density fluctuations in quasars, such as accretion disk instabilities (e.g., [12,13]), interstellar scintillation (e.g., [14,15]), and lensing by the ionized structures (e.g., [16,17]).

Quasars at high redshifts play a crucial role in understanding the evolution of galaxies and the influence of active galactic nuclei (AGNs) in shaping the Universe. However, they are rarely detected in the radio band due to their weak flux densities. These quasars typically exhibit either a steep, flat, or peaked radio spectrum [18,19,20]. The peaked-spectrum sources are classified as high-frequency peakers (HFPs) as well as megahertz and gigahertz peaked-spectrum (MPS and GPS) sources [21,22,23] based on the location of the peak frequency in their radio spectra. The peaked-spectrum (PS) sources represent an early stage in the radio source evolution, where HFPs may evolve to GPSs, and some GPSs may further develop into extended FRI/II radio galaxies [24,25,26,27]. Some other GPS sources exhibit intermittent and short-lived radio emission, offering insights into the dynamic nature of quasars and AGN evolution [28]. Coppejans et al. [29] proposed to search for high-redshift ($z>2$) sources by identifying MPS sources that are compact on scales of tens of milliarcseconds; though later Coppejans et al. [30] found a roughly equal proportion of flat, steep, and peaked radio spectra among 30 quasars at $z>4.5$.

In a study of the radio spectra of over a hundred radio-loud $z\ge 3$ quasars, we found that almost half of the sample can be classified as PS sources, either GPS or MPS [31]. A significant fraction of PS sources in this sample are identified as blazars which are more easily detectable due to the beamed emission from their jets compared to unbeamed radio galaxies and AGNs. Variability serves as an important indicator to distinguish between blazars and unbeamed PS sources. Genuine GPS sources exhibit low levels of variability (e.g., [32,33,34,35,36,37,38,39]), which may be attributed to the expansion of the compact lobes. This expansion usually does not change the shape of the spectrum, although the peak frequency might shift [40]. In contrast, blazars exhibit strong and rapid flux variations in a wide range of timescales and frequencies, often accompanied by changes in the spectrum shape, such as a peaked spectrum during flaring activity (e.g., [39,41,42]). By classifying GPS/MPS sources and blazars among high-redshift radio sources, we can gain insights into the underlying physical mechanisms that govern their formation and evolution, which may have profound implications for our understanding of the early Universe. In this study we determine the incidence of GPS and MPS sources in the early Universe and assess the degree to which they are contaminated by blazars.

VLBI studies of high-redshift quasars reveal two major morphological patterns: a core–jet structure, indicative of relativistic beaming of the jet, and a compact double-lobed morphology, the so-called compact symmetric objects (CSOs) with their jets lying on the plane of the sky and relativistic beaming is not significant. Typically, quasars with the core–jet morphology have flat radio spectra at gigahertz frequencies, while CSOs have peaked radio spectra, small sizes ($<1$ kpc) and are classified as MPS/GPS sources [21,22,23]).

To identify blazars, GPS, and MPS sources, multi-frequency radio measurements provide fundamental insights into the amplitude, time scale, and frequency evolution of the flux density variability. The features of variability can be used to identify blazar candidates as well as GPS and MPS sources. The aim of this paper is to study the radio variability properties of bright quasars at $z\ge 3$ on historical timescales of 30–40 yrs. We used the multi-frequency RATAN-600 measurements obtained from 2017 to 2020 quasi-simultaneously at frequencies of 2.3–22.3 GHz, complemented by available literature data. Utilising the dataset compiled by Sotnikova et al. [31], we present the average variability parameters for the quasars, discuss flux density variations for individual objects, and classify GPS/MPS sources according to their long-term spectral behaviour.

2. Sample

The sample contains 101 quasars selected from the NRAO VLA Sky Survey (NVSS, [43]) in the declination range from $-{34}^{°}$ to $+{49}^{°}$ with the flux density $S_{1.4}\ge 100$ mJy and with redshifts $z\ge 3$. The spectroscopic redshifts were obtained by cross-identification in the National Aeronautics and Space Administration (NASA)/Infrared Processing and Analysis Center (IPAC) Extragalactic Database (NED).1 Observations of the sample were conducted with the RATAN-600 radio telescope during 2017–2020 at frequencies of 1.2, 2.3, 4.7, 8.2, 11.2, and 22.3 GHz quasi-simultaneously, within a 5–6 minute interval [44]. The flux densities were measured with a standard error of about 5–20% at 4.7, 8.2, 11.2 GHz and 10–35% at 1.2, 2.3, and 22.3 GHz [Section 4.1 in [31]. Our analysis of the radio spectra, which incorporate both the RATAN-600 measurements and the literature data, revealed that 47 out of 101 quasars (46%) exhibit peaked radio spectra, while 24% have flat spectra,2 and 15% have steep ones. A small number of the sources have complex, upturn, or inverted spectra. In this paper the spectral indices are used in terms of ${\alpha }_{\mathrm{low}}$ and ${\alpha }_{\mathrm{high}}$ as they were defined in [31]. For peaked and upturn spectra, ${\alpha }_{\mathrm{low}}$ and ${\alpha }_{\mathrm{high}}$ are determined lower and higher the frequency where the spectral slope changes its sign from positive to negative or vice versa. For the rest of spectral types, the low-frequency spectral index is calculated using decimetre-wavelength measurements, such as TGSS, GLEAM, or WENSS. The high-frequency spectral index is estimated using centimetre catalogue data points. A spectrum is considered flat if the spectral index is $-0.5\le \alpha \le 0$, and inverted if $\alpha >0$. The spectrum is assumed to be complex for indeterminate shapes with two or more maxima or minima. The spectrum is called steep or ultra-steep for $-1.1<\alpha <-0.5$ and $\alpha \le -1.1$ respectively.

The list of the sources and their characteristics are presented in Table A1: Col. 1 is the source name; Col. 2 is the redshift; Col. 3 is the average RATAN-600 flux density at 4.7 GHz; Cols. 4–5 are the radio morphology and its reference; Cols. 6–7 are the spectral indices, calculated in [31]; Col. 8 is the radio spectrum type; Col. 9 is the blazar type; Col. 10 is the number of cluster. Forty-eight quasars in the source list are classified as blazars according to the “Open Universe for Blazars—Reference list V2.0”.3, which includes 6077 blazars identified through radio, X-ray, and $\gamma$-ray surveys as well as objects discovered based on their multi-frequency SED properties. This database serves as the largest compilation of blazars to date [45]. It contains the blazars from the Roma-BZCAT catalogue [46], the classification of which we used in our research. Most of the quasars of the sample (88 of 101) have been observed using VLBI at 2, 5, 8, or 15 GHz, with the majority of them (66 of 88) exhibiting a core–jet morphology.

Compared to the sample in [31] we have excluded the source J1600+0412 (87GB 155733.8+042120). J1600+0412 has $z=0.79$ in the latest data release of the Sloan Digital Sky Survey [47,48], but previously its redshift was miscalculated as $z=3.11$ [49].

There are six gravitationally lensed quasars in the sample.4 Their structure is unresolved by RATAN-600 because the maximum separations between them are in the range of 0.54″–3.4″ which is less than the angular resolution at any RATAN-600 frequency. The notes for them are given in Appendix C.3.

2.1. The Literature Data

The literature data and their use in our study have been previously described in detail in [31]. The vast majority of the external radio continuum measurements were obtained from the astrophysical CATalogs support System (CATS database,5 [50,51]). These data allow us to estimate radio flux variability over a time period of up to 30-40 years at the most representative observation frequency of 4.7 GHz. The external data were adopted within 10% ranges of the RATAN-600 five frequencies. For example, at $\nu =4.7$ GHz we used all data within $4.23$–$5.17$ GHz. The majority of the quasars (94/101) has more than 30 years of observing coverage, and 78 out of 101 have more than 10 radio measurements at 4.7 GHz (Table A2). The median values of the number of measurements (N${ }_{\mathrm{obs}}$) and the years of observations at 4.7 and 22.3 GHz (t) for the sample are 16/35 and 9/14 respectively (Table 1). In our study we define the timescale of radio monitoring in the rest frame as $t_{\mathrm{rest}}=t_{\mathrm{obs}}/(1+z)$; if not specified, the timescale is used in the observer’s frame of reference.

3. Multi–Band Radio Variability

3.1. Variability Estimates

In order to characterise the variability properties of the quasars at different frequencies, we have computed the variability and modulation indices and the fractional variability. The first and third characteristics take into account measurement uncertainties, while the modulation index and the fractional variability are less sensitive to outliers. The variability index was calculated using the formula from [52]:

$$V_S=\frac{(S_{\max }-{\sigma }_{S_{\max }})-(S_{\min }+{\sigma }_{S_{\min }})}{(S_{\max }-{\sigma }_{S_{\max }})+(S_{\min }+{\sigma }_{S_{\min }})}$$

(1)

where $S_{\max }$ and $S_{\min }$ are the maximum and minimum flux densities over all epochs of observations; ${\sigma }_{S_{\max }}$ and ${\sigma }_{S_{\min }}$ are their errors. This formula prevents one from overestimating the variability when there are observations with large uncertainties in the data. We obtain a negative value of $V_S$ in the case where the flux error is greater than the observed scatter in the data.

The modulation index, defined as the standard deviation of flux density ${\sigma }_S$ divided by the mean flux density $\overline{S}$, was calculated as in [53]:

$$M=\frac{{\sigma }_S}{\overline{S}}$$

(2)

The fractional variability $F_S$ is defined as in [54]:

$$F_S=\sqrt{\frac{V^2-{\overline{\sigma }}_{err}^2}{{\overline{S}}^2}}$$

(3)

where $V^2$ is the variance, $\overline{S}$ is the mean flux density, and ${\sigma }_{err}$ is the root mean square error. The uncertainty of $F_S$ is determined as:

$$△F_S=\sqrt{{\left(\sqrt{\frac{1}{2N}}\frac{{\overline{\sigma }}_{err}^2}{F_S\ast {\overline{S}}^2}\right)}^2+{\left(\sqrt{\frac{{\overline{\sigma }}_{err}^2}{N}}\frac{1}{\overline{S}}\right)}^2}$$

(4)

where N is the number of observations.

The distributions of the variability index, modulation index, and fractional variability at 2.3–22.3 GHz are shown in Figure 1. At 1.2 GHz we did not obtain sufficient number of measurements to estimate variability. Table A2 lists these parameters for each quasar, while their statistics are given in Table 2. The negative variability indices have been excluded from the calculation and distribution diagrams.

Light curves of the quasars obtained with RATAN-600 are presented in Appendix D. For PS candidates we included additional RATAN-600 measurements of 2006–2017 from [36,38]. Their radio spectra, compiled using both the RATAN-600 and CATS radio data, are also presented in Appendix D.

The number of observing epochs N${ }_{\mathrm{obs}}$ and the timescale t play a key role in exploring AGN variability. Previous studies have known that variability tends to increase with a greater number of observations, as sparser sampling may lead to missed peaks of variability [55,56,57]. This finding is verified by our observations, i.e., we found a significant correlation between N${ }_{\mathrm{obs}}$ and the variability index at 4.7 GHz with Pearson’s $r=0.35$, p-val. $<{10}^{-3}$ (Figure 2). The observations of quasars in our sample are not homogeneous and differ in timescales (see the parameters $t_{4.7}$ and N${ }_{4.7}$ in Table A2), affecting the revealed variability level. The largest number of observations at 4.7 GHz is N${ }_{\mathrm{obs}}=187$ ($V_{S_{4.7}}=0.67$) for the blazar J0646+4451.

We note that the distributions of M and $F_S$ are similar (see Figure 1). Both these parameters are suitable for a large number of observations N${ }_{\mathrm{obs}}$; however, unlike the modulation index, the fractional variability takes measurement errors into account. Therefore, the third parameter, $V_S$, can be used for any number of observations. For these reasons, we will use the variability index $V_S$ and the fractional variability $F_S$ in our analysis, as these parameters can be effectively utilized regardless of the number of observations.

3.2. Refractive Interstellar Scintillation

The variability observed in this study can be caused by either intrinsic (related to the source properties) or extrinsic factors arising from the interaction of radio emission from a source with an inhomogeneous propagation medium.

The flickers caused by scattering on the inhomogeneities of the interplanetary medium have a characteristic timescale of the order of a second or less [58], so they are smoothed out during RATAN-600 observations due to the relatively long scan time of 3–5 minutes during the source’s transit. Similarly, for measurements at other radio telescopes, interplanetary flickering is also insignificant, given the longer measurement times for continuum flux density.

The scintillations caused by the propagation of radio waves in the interstellar medium can have a diffractive and refractive character. Diffractive scintillation is usually associated with an extremely compact source, such as a pulsar or a fast radio burst. Additionally, observations must be made in a narrow frequency band in order to resolve the fine-scale structures responsible for the diffractive scintillation [59]. This is because the scintillation pattern can be highly frequency-dependent. The continuum radio measurements are carried out in a wide band (from hundreds of MHz to several GHz); therefore, this type of flickering can not be detected in our data.

Let’s focus on estimating the potential contribution of refractive interstellar scintillations (RISS) to the variability at the RATAN-600 frequencies (1.2, 2.3, 4.7, 8.2, 11.2, and 22.3 GHz) for the sources from the sample (in the declination range from $-34°$ to $+49°$). According to [14], the transition frequency separating strong and weak scattering modes lies in the range from 5 to 15 GHz for most quasars. We will adopt these transition frequency values to estimate the modulation level of the flux density and its timescale. We estimated an angular size as ${\theta }_{\mathrm{s}}=5\times {\theta }_{\min }$, where ${\theta }_{\min }=0.6\sqrt{S}/\nu$ is the minimum angular size for the stationary source of synchrotron radiation [60], S is the flux density (in Jy) at the observer’s frequency (in GHz). We adopted $S=100$ mJy, taking into account our sample selection criterion and characteristic flux densities measured with RATAN-600. By employing the formulas outlined by [14], we calculate the flux density modulations and their typical timescales, as shown in Table 3. We can conclude that for quasars far from the Galactic plane, the modulation level is about 1 per cent, which is significantly lower than the average values of the variability level for our sample. Therefore, refractive scintillations contribute insignificantly to the overall variability for most of these objects, except for a few individual cases where variability reaches $\le 5$ per cent. For objects located close to the Galactic plane (the transitional frequency around 15 GHz), refractive scintillations can make a significant contribution at low frequencies of 1–2 GHz. However, at higher frequencies ($>5$ GHz), this contribution is considerably smaller (again, for individual objects with a level of variability of no more than 10%, the contribution of refractive scintillations is significant but not decisive). Only three quasars in the sample are relatively close to the Galactic plane, within the Galactic latitude $-15°<b<+15°$: J0624+3856, J0733+0456, and J2050+3127. Their variability level at 2.3 GHz is in the range of 19–47%, which may have a contribution from RISS but can not be explained only by it.

Thus, with the exception of a few specific objects, the contribution of interstellar scintillation effects to the total variability of most quasars in our sample is insignificant (see Table 3).

3.3. Variability Characteristics

We computed variability amplitudes at certain observing frequencies. They correspond to different rest frame frequencies (Figure 3, top), since the sources span a large redshift range. Furthermore, the timescale of the observations spans from about half a month to 7 years at 22.3 GHz and from 4 to 11 years at 4.7 GHz in the source’s rest frame (Figure 3, bottom), and consequently the median values of the variability timescale are 3 and 8 years at observing frequencies of 22.3 and 4.7 GHz respectively (Table 1). We found fairly wide scatter of the variability index throughout all the range of 9–130 GHz rest frame emission frequencies (Figure 4). Those emission frequencies are all associated with the dominating core component. There are no trends in the $V_S$ to ${\nu }_{\mathrm{rest}}$ relation, confirmed by Pearson’s r=0.12 (p-val.=0.02).

We obtained the average variability level for the sample to be 0.20–0.30 at five frequencies (Table 2). For all spectral types, the average variability is higher at 22.3 GHz, up to $\sim 0.40$ (Table 4). Maximum variability is observed in sources with flat, peaked, and complex spectra. Kolmogorov-Smirnov test confirms that the $V_S$ and $F_S$ distributions are not significantly different for the flat and peaked spectrum quasars.

We did not find a correlation between redshift and variability $V_S$ at 4.7 and 22.3 GHz. At redshifts greater than 3.5, we obtained the spread of variability indices at 4.7 GHz almost comparable to that at smaller redshifts (Figure 5). But there are scarce statistics on both the number of the highest-redshift quasars and the time scale of their radio monitoring. At 4.7 GHz ten extremely distant ($z>4$) sources exhibit a wide range of variability on a long time scale $t_{\mathrm{rest}}\sim 5$–10 yrs. The variability index varies from 0.12, which was obtained for the bright VLBI quasar [61] J2134−0419 ($z=4.34$, N${ }_{\mathrm{obs}}$=18, $t_{\mathrm{rest}}\sim$ 6 yrs), to 0.71 for the peaked spectrum blazar J0525−3343 ($z=4.41$, N${ }_{\mathrm{obs}}$=12, at $t_{\mathrm{rest}}\sim$ 5.5 yrs) with an unresolved or slightly resolved core–jet structure in VLBI images [61]. The most distant blazar in the sample ($z=5.28$) J1026+2542 with a one-sided jet morphology [62] has a well determined peaked spectrum and a variability index of 0.30 over $t_{\mathrm{rest}}\sim 3.9$ yrs. However, even on a scale of 6–10 yrs in the rest frame there are about 10% of quasars with $V_{S_{4.7}}<0.10$ (Figure 6). Among them, about one third are PS quasars.

For averaged spectra of the quasars measured over the entire historical time period, we did not find a correlation between the variability at 2.3, 4.7, 8.2, 11.2 and 22.3 GHz ($V_S$, $F_S$) and the spectral indices ${\alpha }_{\mathrm{high}}$ and ${\alpha }_{\mathrm{low}}$ in terms of Pearson and Spearman correlation coefficients. Such correlation would be expected for AGNs due to the fact that flat spectrum radio sources ($\alpha >-0.5$) are more variable than the steep ones. Indeed, a strong correlation between the spectral index and variability of bright AGNs has been found: for example, [69,70]. This lack of correlation is possibly due to the prevalence of variable PS quasars with steep averaged spectra in the optically thin emission part. In addition, the blazars make up at least half of the sample and almost half of them have a peaked spectral shape [31]. On a long time scale, these PS blazars demonstrate high radio variability (0.25–0.75) near the peak frequency. In Figure 7 the “$\alpha$–$V_S$” plane is divided into four quadrants: steep-spectra with $V_S>0.25$ and $<0.25$, flat-spectra with $V_S>0.25$ and $<0.25$. It is clearly seen that PS quasars (orange and blue stars) make up the vast majority in the steep-spectrum area with $V_S>0.25$. During a flare, the radio spectra of some sources can exhibit characteristic spectrum evolution. In some cases, the spectra move up while maintaining a peaked shape, while in the others, they first move towards a flatter shape with increasing flux density at high frequencies before moving to the steep-spectrum area after the flare, when emission becomes optically thin.

An example illustrating this behaviour is the spectrum evolution of the PS blazar J0646+4451 with $V_{S_{8.2}}=0.74$ and ${\alpha }_{\mathrm{high}}$=-0.62. During the 17 years of RATAN-600 observations (Figure 8), its spectral index ${\alpha }_{\mathrm{high}}$ changed from $-1.51$ to $+0.35$, which corresponds to a spectral index variability of 0.60. We estimated the variability of the spectral index $V_{\alpha }$ using the same equation [1] in terms of the ${\alpha }_{\nu }$ maximum and minimum values with their uncertainties. The blazar had an ultra-steep spectrum ($\alpha =-1.51$) in March 2011 at the beginning of its flux density decreasing at 2.3–22.3 GHz (orange spectrum), when radio emission had become optically thin at the frequencies above 5 GHz after several years of a relatively high state. There is a strong correlation between ${\alpha }_{\mathrm{high}}$ and the flux density at 11.2 GHz (Pearson’s $r=0.81$, p-val.$<{10}^{-3}$). Similar behavior is observed in other PS blazars such as J2316−3349 and J0525−3343. These spectral evolution patterns are probably caused by changes in the physical conditions in the source during a flare, such as variations in the magnetic field strength or electron density, which affect the observed flux densities at different frequencies.

We found significant correlation (Pearson’s $r>0.5$, p-val.$\ll 0.005$) between the spectral index ${\alpha }_{\mathrm{high}}$ and the flux density variability at 11.2 GHz for 15 of the most frequently observed (N $>10$) quasars with $V_{S_{11.2}}>0.25$. This indicates that for these quasars, changes in the spectral index are predominantly driven by variation in the flux density at and above 11.2 GHz. Further studies are needed to investigate the underlying physical processes driving the observed variability in these quasars.

The estimates of AGN variability at different redshifts can be biased due to the varying cosmological time dilation and different frequencies in the rest frame. The poorly sampled light curves also strongly impact variability estimates. Several factors, such as the time scale of monitoring $t_{\mathrm{rest}}$, the frequency ${\nu }_{\mathrm{rest}}$, and the number of measurements N${ }_{\mathrm{obs}}$ should be taken into account to compare variability at different z [71,72]. In order to search for patterns and tendencies by localizing groups of similar objects in a multidimensional parameter space (feature space) and investigating their statistical properties, we used two independent methods of cluster analysis of the quasars (Appendix B). The first uses PCA dimensionality reduction with k-means clustering, the second is based on the self-organizing maps [73,74]. We chose the variability parameters at an observed frequency of 4.7 GHz, which varies from 18 to 30 GHz in the rest frame. As a result, five quasar groups were determined and their parameters described in Table A3, which allowed us to compare the variability features within the clusters. An extremely highly variable quasar group with median $V_{S_{4.7}}=0.48$ was found at the highest redshifts, $z=4.3$–$5.3$, and they were observed as well at the shortest time scales, the median $t_{\mathrm{rest}}=5.7$ yrs. The clustering algorithm also determined a group of the most observed quasars with $N_{\mathrm{obs}}=24$–187. These quasars are at intermediate $z\le 3.9$ with long $t_{\mathrm{rest}}>7$ yrs. It is interesting that these most monitored quasars do not generally show the highest variability indices, which are nevertheless statistically higher than those for low-variable cluster 2 and intermediate-variable cluster 1. It is clearly seen that two factors impact the clustering result: the evolutionary state of the objects and observational constraints, such as $N_{\mathrm{obs}}$ and $t_{\mathrm{rest}}$.

4. Peaked-Spectrum Sources at High Redshifts

Peaked spectrum quasars make up half of our sample Sotnikova et al. [31]. We consider the source as PS when the spectral index changes its sign from positive to negative in the average spectrum compiled with both the RATAN-600 and literature data. It is well known that only several per cent of genuine GPS sources have been found in complete samples of bright AGNs spreading in a wide range of redshifts: see, e.g., [37,38,75,76]. It is commonly believed that the mechanism which produces the absorption in the optically thick part in these compact luminous objects is due to synchrotron self-absorption caused by the high density of synchrotron-emitting electrons [77]. Additionally, free-free absorption (FFA) has been suggested as a mechanism that can produce the GHz-frequency turnover as well as the observed anticorrelation between peak frequency and size in many of these sources [78]. Several explanations have been put forth to account for the properties of the genuine GPS sources: they might be young radio galaxies in the process of evolving into larger ones [27,79], or their compact nature could be a result of confinement through interaction with a dense environment [80,81,82,83,84], or the sources might be recurrent or intermittent sources [27,85,86,87,88].

While most scenarios are well applicable to the known GPS galaxies, the nature of the GPS quasars has received limited attention in the literature. GPS quasars are thought to be quite rare and challenging to identify due to strong contamination from blazars. Compact beamed sources can exhibit temporarily inverted radio spectra and may be easily disguised as GPS sources [35,89]. Study of the radio variability of PS quasars on historical timescales allows one to detect genuine GPS and MPS quasars, to estimate their fraction in the population of bright distant AGNs, and also to check how the MPS method is effective in searching for high-redshift AGNs [22,29].

To search for GPS/MPS sources, we employ the following common criteria: the observed peak frequency ${\nu }_{\mathrm{obs},\mathrm{peak}}$ is less than 1 GHz for MPS sources and from 1 up to 5 GHz for GPS sources; spectral indices of optically thick ${\alpha }_{\mathrm{thick}}$ and thin ${\alpha }_{\mathrm{thin}}$ emission parts are close to $+0.5$ and $-0.7$, respectively (for the peaked spectra they match with ${\alpha }_{\mathrm{low}}$ and ${\alpha }_{\mathrm{high}}$, respectively); radio variability does not exceed 25% on a long timescale, as GPS/MPS sources are known as the least variable class of compact extragalactic radio sources [21,23,33,90,91,92]. Another PS type is high frequency peakers (HFPs), which are considered to be powerful and intrinsically compact (<1 kpc) extragalactic radio sources with a spectral peak ${\nu }_{\mathrm{peak},\mathrm{obs}}$ above 5 GHz, representing the earliest stage of radio source evolution with typical ages of a few hundred years [79,93,94]. According to the criteria in [95,96], the HFPs lacking flux density variability and with low fractional polarization (<1%) and intrinsically small linear sizes (<1 kpc) are considered good candidates for young radio sources.

Table 5 presents the list of PS sources with their peak frequency, the width of the radio spectra at the half maximum level (FWHM), and spectral type based on [31], radio morphology from the existing literature, timescale of radio measurements, and the number of observing epochs.

The distributions of the parameters calculated for PS sources in this section are shown in Figure 9, and their mean values are listed in Table 6. When PS blazars and non-blazars are considered separately, significant difference in the spectral and variability parameters emerge, emphasizing the importance of a careful approach in distinguishing intrinsically different sources that happen to have a spectral peak in the MHz–GHz frequency range. As shown in Table 6, PS blazars have their spectral peaks at higher frequencies and have higher variations of flux densities. The distributions of variability estimates at 2.3–22.3 GHz for the PS sources and blazars of the sample are shown in Figure 1 for comparison with the whole sample.

4.1. Rest Frame Peak Frequency and FWHM

For the sources with a peaked radio spectrum, we estimated the rest frame peak frequency ${\nu }_{\mathrm{peak},\mathrm{rest}}={\nu }_{\mathrm{peak},\mathrm{obs}}(1+z)$, which ranges from 0.8 to 76 GHz. The redshift versus ${\nu }_{\mathrm{peak},\mathrm{rest}}$ relationship is shown in Figure 10, where we additionally plot PS sources at different redshifts from other studies: the sample from [104] marked by cross symbols, sources from [105] marked by star symbols, and from [38] by grey circles. [91] found a trend towards higher turnover frequencies at higher redshifts, but Figure 10 shows a wide spread in the peak frequencies in our sample at redshifts z=3–5 and no evidence for correlation between them. Thus, sources with synchrotron self-absorbed components around and even below 1 GHz in the rest frame occur at high redshifts and at high radio luminosities.

We calculated the Full Width at Half Maximum (FWHM) for the PS sources in units of decades of frequency (dex). For classical GPS sources, the parameter FWHM $\sim 1$–1.5 dex, as has been found in [21,106] and [91] for a sample mostly represented by $z<3$ GPS sources (with only 6 of 74 sources having $z>3$ in the latter study). [37] compared the samples of GPS sources classified as galaxies and quasars and found that the radio spectrum of GPS quasars is on average wider (FWHM $\sim 1.6$) than that of the GPS galaxies (FWHM $\sim 1.4$). In our PS sample sources, the FWHM varies from 1.3 to 7.7 dex with a mean value of $1.93\pm 0.96$ (Table 6).

4.2. Spectral Variability of Identified HFP/GPS/MPS Quasars

Blazars have the potential to contaminate GPS source samples since they exhibit peaked radio spectra during their flaring states [39,41,107,108]. In this section, we focus on the broadband long-term radio spectra of 23 quasars previously classified as HFP, GPS, or MPS sources [21,33,37,105,109,110,111,112,113]. Their radio spectra and RATAN-600 light curves are presented in Figs. Figure A11–Figure A12. The longest-term measurements, spanning up to 30–43 years, are available at frequencies at 4.7, 11.2, and 22.3 GHz. They are listed in Table 7 with the information about their variability index at 11.2 and 22.3 GHz and the historical period of observations.

Among the known HFP/GPS/MPS sources, six are not classified as blazars. Three of these six sources (J0048+0640, J0751+2716, J0753+4231) have relatively small variability indices $V_{S_{4.7}}=0.04$-$0.25$ and $V_{S_{11.2}}=0.05$-$0.12$ over the timescales of 36–43 and 3–9 years, respectively. Two of them have a compact symmetric morphology [98,101]. J1606+3124 has been classified as a galaxy with a radio structure similar to that of CSOs, and it may be the highest redshift ($z=4.56$) CSO galaxy discovered so far [101]. We measured the variability indices as 0.10–0.16 at 11.2 and 22.3 GHz and 0.28 at 4.7 GHz for this source at a long timescale (37–43 yrs). The other two objects, J1340+3754 and J1538+0019, have the variability indices $V_{S_{4.7}}=0.30$–$0.50$ ($t=36$ and 43 yrs) and $V_{S_{11.2}}=0.03$–$0.27$ ($t=5$ and 23 yrs).

Among seventeen GPS quasars in Table 7, which belong to the blazar population, thirteen demonstrate high variability indices $V_{S_{22.3}}$ up to 0.67. Only four blazars, namely J1424+2256, J1616+0459, J2248−0541, and J2316−3349, have relatively low variability of less than 0.26 at 4.7 GHz and less than 0.27 at 11.2 GHz. However, for J2248−0541 and J2316−3349, there are no enough radio data to establish their variability reliably.

Thus, among previously reported 23 GPS, MPS, and HFP sources, seven quasars demonstrate low radio variability ($V_S<0.27$) at 4.7–22.3 GHz, with four of them being blazars. Only two GPS (J0048+0640, J1606+3124) and two MPS sources (J0751+2716, J0753+4231) have variability attributes of young compact radio sources. The description of 17 individual sources from this section that have significant radio variability is presented in Appendix C.1.

4.3. New GPS/MPS Candidates

We found 12 new GPS/MPS candidates, which are marked with an asterisk symbol in Table 5 and described in Appendix C.2. Seven of them, J0100−2708, J0839+2002, J0933+2845, J1301+1904, J1401+1513, J1521+1756, and J2314+0201, we consider as new GPS candidates. Five quasars, J0214+0157, J0905+0410, J1045+3142, J1418+4250, and J2019+1127, we suggest as new MPS candidates.

4.4. HFP Candidates

There are 6 PS sources previously identified as HFP candidates (Table 5), with a turnover peak at about 5 GHz or higher. If we consider the source’s rest frame frequency, it is above 20 GHz. All of them are classified as blazars of the FSRQ type in BZCAT and have high variability indices from 0.26 to 0.65 at 22.3 GHz according to our estimates. Four quasars, J0646+4451, J1616+0459, J2003−3251, and J2129−1538, have been previously classified as HFPs (Appendix C.1). For the remaining two quasars, J0733+0456 and J1503+0419, we obtained the turnover frequency ${\nu }_{\mathrm{peak},\mathrm{obs}}$ equal to 5.8 and 5.6 GHz respectively [31].

J0733+0456 was not previously considered as an HFP candidate. The multifrequency data from the GLEAM survey at 72–231 MHz in 2013–2014 [114] and the RATAN-600 measurements in 2016–2020 made it possible to classify the quasar as a PS source [31] with a high turnover frequency ${\nu }_{\mathrm{peak},\mathrm{obs}}=5.8$ GHz. The wide spectral peak FWHM $=1.8$ reflects the high radio variability of the quasar at all the considered frequencies ($V_S>0.30$).

J1503+0419, initially identified as an FSRQ in [115], exhibits a core and some extension in the VLA 5 GHz image. Later authors [116] confirmed the presence of a lobe in the 1.4 GHz FIRST survey image. However, when observed using VLBI, no discernible structure was detected, as reported in [97]. Unfortunately, there are very little radio data. We measured a low variability index $V_S=0.10$ at 4.7 GHz spanning a time period of $t_{4.7}=39$ years with a total of 8 observations. The variability at 11.2 and 22.3 GHz can not be evaluated because the flux density errors are larger than the observed scatter of the data, according to the only two RATAN-600 observing epochs. The quasar was classified as a PS source for the first time with an estimated peak frequency ${\nu }_{\mathrm{peak},\mathrm{obs}}=5.6$ GHz and FWHM $=1.5$ [31].

Despite the presence of a peak above 5 GHz, we do not consider three quasars as HFP candidates: J0539−2839, J1405+0415, and J1658−0739, since they have two-component spectra or a strongly variable spectral shape with a wide peak (Table 5, Table A2).

Note that both the quasar J0257+4338 with a core–jet morphology [117] and the BL Lac blazar [118] or blazar of an uncertain type in BZCAT with an unresolved core in VLBA images [119]. J0525−2338 could hypothetically be considered as HFP candidates due to their spectral shape, which has a hint of a possible peak near 8 and 20 GHz [31]. The blazar J0525−2338 exhibits a rising spectrum (${\alpha }_{\mathrm{high}}$=+0.60) and has the variability indices $V_{S_{22.3}}=0.50$ and $V_{S_{11.2}}=0.36$. J0257+4338 is one of the most distant quasars in the sample at $z=4.06$[120] with a variability level $V_{S_{8.2}}=0.45$. Its turnover frequency showed a remarkable change from 2 GHz up to 8 GHz between 1990 and 2017, with a strong increase in flux density from 0.13 Jy [121] up to 0.36 Jy.

5. Discussion

Studying the variability of quasars in the radio band is an approach to understanding the astrophysical processes occurring in the sources. Long-term (months, years, and decades) variability provides insights into the properties of the jet, flaring activity, and interaction between the jet and the surrounding medium, while short-term (from several hours to several days) variability could provide information about the interstellar medium based on the measurements of scintillations [122]. In our study, we consider flux density variations at GHz frequencies over long periods of time of a sample of high-redshift AGNs, which represent a superposition of the aforementioned short-term events and long-term processes.

The redshift values range from 3.0 to 5.3 for the quasars in our sample, which corresponds to a wide range of emission frequencies in the source’s reference frame, from 9 to 130 GHz, and the rest frame time scales of monitoring from several weeks to 11 years. Consequently, we can get the most general picture of the dependence of the radio emission variability on the redshift and frequency and on the time scale of monitoring.

The quasars that have showed a high level of variability but are not listed as blazars could be considered as blazar candidates after checking their morphology and other multi-band properties, thus usefully supplementing the current lists of the known blazars in the “Open Universe for Blazars – Reference list V2.0” and the Roma-BZCAT catalogue.

5.1. Variability Properties

Estimates of the flux density variability show that the averaged variability is higher at higher frequencies ($M_{22.3}=0.28$, $F_{S_{22.3}}=0.31$, and $V_{S_{22.3}}=0.36$) with typical values of about 0.20–0.25 at other frequencies. These values are comparable to those observed in blazars at any redshift.

The distributions of $V_S$ and $F_S$ at 22.3 GHz are significantly different from those at 2.3–11.2 GHz, which may be attributed to intrinsic physical differences. At the source’s rest frequencies greater than 45 GHz (i.e., $11.2\times (1+z)$, where $z\ge 3$), the radio emission from the compact AGN core becomes more dominant, leading to a different variability behavior (the $V_S$ distribution peaks around 0.4) compared to the lower frequencies (the $V_S$ distribution peaks around 0.1–0.2). At frequencies lower than 45 GHz, the variability may be mixed with variations from multiple jet components, which can lead to smearing of the variability index.

Only the brightest peaked-spectrum sources can be detected at high redshifts due to their steep radio spectrum above the peak frequencies. Compared to the other sources in the sample, the PS sources are about two times brighter, while the variability level is of the same order at all frequencies. The PS sources are also found to be more variable than the AGNs with power-law spectra at MHz frequencies in [123] and [124]. Some PS sources even demonstrate high amplitude variability of 100–2500% over 24 years at 1.5 GHz [125]. More than half of PS quasars in the sample are blazars (25 out of 47), which affects the variability level of the PS subsample.

We have not found a correlation between radio variability and redshift or the spectral index; apparently, dense and long-term time series are required to study such relationships more robustly. Comparing the variability of high-z AGNs, the contribution of observational restrictions and different source conditions in such a wide redshift range can not be easily distinguished for our sample. We found several quasar groups with distinctive variability features. The most interesting is a small highly variable group of quasars located at redshifts of 4.3–5.3, notable by a shorter time scale of radio monitoring, $t_{\mathrm{rest}}=4.0$–$6.7$ yrs, compared to the rest of the quasars in the sample. One of the simple explanations is that the youngest quasars are more compact, and as a consequence the variability time scale is shorter. This can be confirmed by the recently reported relationship between variability timescales and core sizes in [126]. The observational factors $t_{\mathrm{rest}}$, ${\nu }_{\mathrm{rest}}$, and N${ }_{\mathrm{obs}}$ impact the obtained variability properties first of all, and we should be careful making far-reaching conclusions.

5.2. Population of Peaked Spectrum Sources at High Redshifts

A comparison of radio spectra features and flux density variations of quasars reveals that the peaked spectrum is a typical feature of blazars at redshifts $z>3$, distinguishing them from the flat spectrum low-redshift blazars. Obviously, only a small fraction of quasars can be associated with the so-called classical GPS/MPS sources. The long-term radio data showed a significant variability at GHz frequencies, up to 70 per cent, in the most known GPS quasars (Section 4.2). Previous studies have mentioned the scarcity of genuine GPS quasars [75,90,103], and we need long-term multifrequency observations to determine their spectra reliably.

Because of the young age of HFPs, one could expect a larger proportion of them at high redshifts compared to the intermediate ones. Previously, in [37], we identified about 2% objects in the flux-density-complete sample of GPS galaxies and quasars at lower redshifts (69% sources are at $z\le 2$) as HFP candidates. In our sample of high-redshift quasars, all 6 HFP candidates are FSRQ blazars, which exhibit high variability levels and show a peaked spectrum shape occasionally, during a flare. These six quasars are, on average, more variable, with the mean value of $V_{S_{22.3}}=0.43\pm 0.14$, compared to the whole sample ($V_{S_{22.3}}=0.35\pm 0.18$). Authors [111] found that only 40% of the objects displayed all the typical characteristics of young radio sources. Furthermore, all quasars in the sample were excluded. So, to prevent contamination from selection and boosting effects, authors [104] considered only galaxies in their statistical analysis of HFP sources.

Three quasars with peaked radio spectra and a CSO-type morphology showed relatively low variability levels at all the five frequencies we considered, providing evidence in favor of them being bona-fide CSOs. These findings are consistent with the estimates in [127,128], where they correspond to $\sim 3$% of the sample.

It should be noted that for 12% of the objects there are no multifrequency and long-term measurements, and our analysis is based on non-uniform sparse time series. Only for about 10 quasars there are flux density measurements with high cadence (N${ }_{\mathrm{obs}}=40$–60) over a long time period. The strong correlation between the number of observations and variability level confirms the influence of the observational selection effect on the result.

6. Summary

In this study, we present the results of a multi-frequency radio variability analysis of 101 bright radio quasars at $z\ge 3$ based on the RATAN-600 observations at 2.3–22.3 GHz conducted in 2017–2020 [31], as well as the information from the literature covering a time period of 30–40 years, which corresponds to a timescale of about ten years in the source’s frame of reference. The main results are summarised as follows.

• The quasars show a wide spread of radio variability on any time scale of the measurements. About half of the objects in the sample show a variability level (25–50%) comparable to that of the blazars at low redshifts. Some quasars demonstrate flux density variations of up to 60–90% at 8.2–22.3 GHz. The median values of modulation index M, fractional variability $F_S$, and variability index $V_S$ are about 0.20 at 2.3–8.2 GHz and are higher at 22.3 GHz (0.28–0.36). The proportion of sources with low variability of $V_S<0.10$ varies from 8 to 20% at different frequencies.
• The distribution of the variability index is significantly different at 22.3 GHz compared to lower frequencies, demonstrating the peak at higher levels of variability. This suggests that the emission from the compact core dominates at higher frequencies, while at lower frequencies, the emission from jet components could blend with the core emission, shifting the distribution peak to lower values. In general, the variability we obtained at the range from dozens to about a hundred GHz in the rest frame is attributed to the dominating core component in quasars.
• Several quasar groups with distinctive variability characteristics were found. Six highest-redshift quasars demonstrate average variability of about 0.50 and at a shorter time scale than the others.
• Among the 23 GPS and MPS sources reported previously, only six objects show low radio variability ($V_S<0.25$). We propose seven new candidates for GPS sources (J0100−2708, J0839+2002, J0933+2845, J1301+1904, J1401+1513, J1521+1756, and J2314+0201) and five new MPS candidates (J0214+0157, J0905+0410, J1045+3142, J1418+4250, and J2019+1127) based on the analysis of their spectral parameters and radio variability. We found that about 18 and 9% of the sample of bright quasars at $z\ge 3$ are GPS and MPS candidates, respectively.
• A peaked spectrum is a typical feature of blazars at high redshifts, and this distinguishes them from the flat spectrum low-redshift blazars. The radio variability exceeding more than 30% in high redshift PS sources is common for blazars.

Appendix C.1. Identified HFP/GPS/MPS Candidates which Revealed Significant Radio Variability

J0121−2806 is an FSRQ classified as PS in [37]. We calculated that the spectral peak is at 0.3 GHz [31]. Due to the absence of radio data at frequencies from 0.2 to 1.4 GHz, the spectral maximum can not be determined reliably. The quasar demonstrates high flux density variability with $V_{S_{4.7}}=0.43$ and $V_{S_{8.2}}=0.32$. The spectral index ${\alpha }_{\mathrm{thin}}$ was changing from -0.70 to -0.08, remaining almost flat in 2014–2020.

J0203+1134 is an FSRQ known as a GPS source with a peak frequency of about 17 GHz in the rest frame [21]. The RATAN-600 observations confirm the peaked spectral shape in the time period of 2006–2020. However, occasionally it shows an inverted spectrum (e.g., in 2001) and has been classified as having chaotic variability [134]. The variability indices $V_{S_{4.7}}$ and $V_{S_{22.3}}$, which have been measured at a historical scale of 20–40 years, are 0.33 and 0.45 respectively. The broadband spectra are well described by two components, linear at the observer frequencies of 0.07–0.5 GHz and convex at higher frequencies.

J0324−2918 was considered as a PS source in Coppejans et al. [30], where authors estimated a possible spectral turnover between 1.4 GHz and 7 GHz. The uncertainty of positioning the peak was caused by the scatter of the observed flux densities around 5 GHz, where the AT20G and PMN surveys contributed with different resolutions and at different epochs of observations. In our study we clarified the 2.3–22.3 GHz part of the spectrum as being inverted using quasi-simultaneous RATAN-600 measurements. Thus we see that the radio spectrum of J0324−2918 has a complex structure: below 1 GHz it has a peaked shape, but above it is inverted. We obtained significant variability (up to 45%) for the source during our monitoring, which is in agreement with the conclusion in Coppejans et al. [20] that its VLBI emission is Doppler-boosted.

J0646+4451 is one of the brightest high-redshift PS blazars. The object is known as a high-frequency peaker [111,112] and has for the first time been considered as a GPS source in [109]. The quasar has a spectral peak around 17–18 GHz in the observer’s frame of reference [31,75]. Flux variations are observed on both sides of the peak (Figure A11), and the variability index reaches 0.70 at 5–8 GHz. We observed variable peaked spectra for the blazar during 2006–2011. In 2017–2021 the quasar had flat and even rising spectra, i.e., the spectrum became optically thick at the high frequencies more than 10 GHz, with ${\alpha }_{\mathrm{thick}}=+0.22$ and $+0.35$ in April and August 2020. Its RATAN light curves show a slow flux density decrease about 1.5–2 times, from 4 to 1.5–2 Jy, at 4.7–22.3 GHz in 2012–2020. The steepest radio spectrum ${\alpha }_{\mathrm{thin}}=-1.58$ was observed in 2011, when the flux density started decreasing at 2.3–22.3 GHz, and the spectrum became optically thin at frequencies above 5 GHz after several years of a relatively high state.

J1230−1139 is an FSRQ and has been classified as a GPS candidate in Mingaliev et al. [37]. Low frequency data are well represented by the GLEAM and TGSS surveys. The VLBA observations show a large jet bending ($\sim 90$) within 2 mas from one bright radio component at 8.4 GHz [98]. We have found the variability indices at 4.7 and 22.3 GHz equal to 0.47 and 0.41 on timescales of >43 and 14 years respectively. The RATAN-600 observational data covers the period of 2015–2020. The overall number of observations N${ }_{\mathrm{obs}}=51$ at 4.7 GHz.

J1340+3754 is an FSRQ [115] at $z=3.11$ with a core–jet morphology [97]. It has been classified as a GPS candidate in Mingaliev et al. [37]. Extended monitoring revealed variability up to 30% at our frequencies.

J1354−0206 is a known GPS source and a blazar [21,109] with a core–jet morphology [99]. It demonstrated the maximum variability at 22.3 GHz $V_{S_{22.3}}=0.50$ on a time scale of 23 yrs.

J1357−1744 is a known GPS [21] with a core–jet morphology [99]. It shows quite high flux density variation both at high and low frequencies ($V_{S_{2.3}}=0.37$ and $V_{S_{22.3}}=0.48$).

J1445+0958 at $z=3.55$ has been classified as a GPS source by O’Dea et al. [21] and is also an FSRQ with the core–jet type of radio morphology [98,100]. The compiled spectrum has a well-constrained absorbed part with ${\alpha }_{\mathrm{thick}}=+0.68$ ($0.01$). A high variability index $V_{S_{22.3}}=0.67$ was calculated using 27 observing epochs during a time period of about 20 years. At lower frequencies the variability has a smaller amplitude, $V_S=0.15$ at 2.3 and 4.7 GHz.

The quasar J1538+0019 was classified as a GPS source in [37]. It has flux density variation at 2.3–22.3 GHz in the range $V_S=0.27$–0.58 on a timescale of up to 40 years [135,136]. Its observed peak frequency shifted from 3 to 1 GHz during 2015–2020.

J1559+0304 at $z=3.89$ with a compact core and a two-sided jet at 5 and 8.4 GHz [100] was classified as a PS source with the peak between 72 MHz and 1.4 GHz by Callingham et al. [105]. We have obtained a greater value of the spectral peak ${\nu }_{\mathrm{peak},\mathrm{obs}}=4.5$ GHz and the variability index $V_{S_{22.3}}=0.36$ on a timescale of eight years. The highest variability is observed near the spectral peak, historical 37-year data show $V_{S_{4.7}}=0.41$ (N $=63$). The spectral shape seems to remain convex during the RATAN observations in 2015–2020. But since the beginning of 2019, the flux density has started to increase slowly at 22.3–4.7 GHz.

J1606+3124 is a radio galaxy at $z=4.56$ with a peak in the radio spectrum at ${\nu }_{\mathrm{peak},\mathrm{obs}}=2.5$ GHz. The radio properties of J1606+3124 studied in [98,101] suggest that it is a young radio source, the highest-redshift GPS galaxy known to date. The peak frequency in the source’s rest frame corresponds to 13.9 GHz. A two-sided jet structure is revealed by VLBI at 8.4 GHz.

J1616+0459, an FSRQ at $z=3.21$, was classified as a GPS and HFP source in [21,109,111]. Analysing its compiled radio spectrum, we estimated ${\nu }_{\mathrm{peak},\mathrm{obs}}=4.4$ GHz.The variability index $V_{S_{4.7}}=0.26$ on a timescale of more than 40 years (N $=64$), and $V_{S_{22.3}}=0.26$ on the 23-year timescale (N $=22$). The highest variability was found at 8.2 GHz, $V_{S_{8.2}}=0.42$.

J1658−0739 is an FSRQ [115] at $z=3.74$ with a core–jet morphology [99]. It has been introduced as a GPS candidate with ${\nu }_{\mathrm{obs},\mathrm{peak}}=4.8$ GHz in [92] based on the VLA and ATCA data, but our study revealed variability of 13–33%, which casts doubt that J1658−0739 is a canonical GPS source. Also its radio spectrum suggests the presence of another extended emission at lower frequencies, turning up below 200 MHz.

The FSRQ blazar J1840+3900 at $z=3.1$ was classified as an HFP source in several papers [110,111,137], but the compiled radio spectra demonstrate a high level of variability ($V_S\sim 0.52$) at 4.7 and 22.3 GHz on a timescale of 36 years. The RATAN radio spectra had rising shapes in July 2006 and July 2014 at 4.7–22.3 GHz (${\alpha }_{\mathrm{thin}}=+0.22$ and $+0.34$). We conclude that J1840+3900 was misclassified as a GPS/HFP source.

The core–jet FSRQ blazar J2003−3251 at $z=3.77$ was classified as a GPS and HFP candidate [21,91]. We determined its spectral peak ${\nu }_{\mathrm{peak},\mathrm{obs}}=5.6$ GHz. The quasar had been demonstrating moderate variability $V_{S_{22.3}}=0.25$ up to 2019. Since 2019 its flux density has started to increase from 0.42 to 0.85 mJy at 22.3 GHz and from 0.50 to 0.61 at 11.2 GHz, and radio variability has reached 42% at 22.3 GHz. Further long-term radio observations are required to determine its variability properties.

The quasar J2129−1538 at $z=3.26$ is a bright FSRQ blazar, one of the first known high-redshift GPS sources [21,109]. The quasar was studied as an HFP candidate in [138], where its variability was estimated over 1.5 years of observations (4–6 epochs) with modulation indices of 9 and 22.5% at 22 and 43 GHz respectively; its radio spectrum was reported as steep (${\alpha }_{22-43\mathrm{GHz}}=-0.52$). [103] reported typical well-aligned core–jet structure for J2129−1538 at 8, 15 and 43 GHz radio images. The RATAN-600 observations were made in 2012–2020 and showed slow flux density rising at 4.7–22.3 GHz since 2017. The spectral index ${\alpha }_{\mathrm{thin}}$ was changing from $-0.04$ up to $-0.83$ (by 88%). We estimated its variability index to be 43% at 22.3 GHz using data of 20-years-long observations.

Thus, as we have shown above, longer time periods of monitoring of blazars are likely to reveal higher flux density variations.

Appendix C.2. New GPS/MPS Candidates

The quasar J0006+1415 has a structure with an extended component of 10 mas from the core [97]. Its radio spectrum at hundreds of MHz is presented by the TGSS and mostly by the GLEAM measurements, which demonstrate quite a scatter of flux densities, e.g., at 130–170 MHz the variations of flux density reaches 100% (110–220 mJy). Considering that, in some states its radio spectrum could be flat or steep. We obtained the low level of variability indices, 0.03–0.25 at 2.3–11.2 GHz. However, the spectral indices ${\alpha }_{\mathrm{thick}}=+0.02$ and ${\alpha }_{\mathrm{thin}}=-0.21$ make the spectrum close to flat, with FWHM $=7.7$, and barely corresponding to the canonical GPS spectrum.

The quasar J0100−2708 with a core–jet angular structure [99] has not much information about its radio continuum properties in the literature. It has a clear peaked shape (FWHM $=1.7$) in the radio spectra, where the optically thin part, with ${\alpha }_{\mathrm{thin}}=-0.47\pm 0.01$, is represented by seven epochs of RATAN-600 observations. The optically thick part mostly consists of the GLEAM data (${\alpha }_{\mathrm{thick}}=+0.41\pm 0.01$). Both modulation and variability indices do not exceed a value of 0.20 at 2.3–11.2 GHz.

For J0214+0157 the GLEAM and TGSS data for the optically thick part together with the RATAN-600 data for the optically thin part show a peaked spectrum (FWHM $=1.9$) with ${\nu }_{\mathrm{peak},\mathrm{obs}}=0.5$ GHz. We note a significant scatter in the measured flux densities (0.01–0.26 Jy) below 100 MHz, and $V_S=0.23$ at 2.3 GHz. At the frequencies 4.7–11.2 the variability index equals 0.03–0.23. Due to moderate radio variability and the spectral shape, we suggest J0214+0157 as a new MPS candidate.

The core–jet [61] PS blazar of an uncertain type J0525−3343 at $z=4.41$ has a spectral peak ${\nu }_{\mathrm{peak},\mathrm{obs}}=0.9$ GHz and could be considered as an MPS candidate with the spectral indices ${\alpha }_{\mathrm{thick}}=+0.43\pm 0.02$ and ${\alpha }_{\mathrm{thin}}=-0.96\pm 0.02$. But its radio variability reaches 29–71 per cent at 2.3, 4.7, and 8.2 GHz on a time scale of 30 years. The highest variability is found at 4.7 GHz, the maximum flux density S${ }_{\max }=0.21$ Jy was measured in 1990 in the Parkes-MIT-NRAO (PMN) survey at 4.85 GHz [139]. In the NVSS survey at 1.4 GHz, carried out several years later [43], the flux density is 0.188 mJy, which gives a rising non-simultaneous spectrum.

The megahertz-peaked spectrum core–jet blazar J0624+3856[99] has a spectral peak ${\nu }_{\mathrm{peak},\mathrm{obs}}=0.3$ GHz and FWHM $=2.5$ dex. The highest variability $V_S=0.43$ is found at 8.2 GHz on a timescale of 30 years.

We consider J0839+2002 as a new GPS candidate despite some deviation from the criteria of genuine GPS sources. The low-frequency spectral part is poorly represented by a single TGSS measurement at 150 MHz [140]. The averaged spectral indices of the optically thick and thin parts are $+0.36\pm 0.01$ and $-0.25\pm 0.01$. The variability indices are 0.05, 0.21, 0.15, and 0.29 at 22.3, 11.2, 8.2, and 4.7 GHz respectively.

The core–jet quasar J0905+0410 was detected in X-rays [141], its radio structure has an extension of 10 mas [97]. Its radio spectrum has a clear peak around ${\nu }_{\mathrm{obs}}\sim 0.5$ GHz with the spectral indices ${\alpha }_{\mathrm{thick}}=+0.57\pm 0.06$ and ${\alpha }_{\mathrm{thin}}=-0.53\pm 0.01$. The variation of flux densities is $V_S=0.04$–$0.17$ at 2.3–11.2 GHz (Table A2). We suggest J0905+0410 as the second new MPS candidate in the sample.

The next new GPS candidate J0933+2845 with narrow absorption lines [142] is detected in X-rays, but its structure is unresolved [143]; in the radio band it has an extended structure [100]. Its radio spectrum has a clear peak at ${\nu }_{\mathrm{peak},\mathrm{obs}}\sim 1.7$ GHz. There is a lack of measurements for the optically thick part (${\alpha }_{\mathrm{thick}}=+0.17\pm 0.02$), we rely on only one measurement point at 150 MHz from the TGSS. At frequencies higher than 2.3 GHz we calculated flux densities variability in the range $V_S=0.21$–$0.34$.

The source J1045+3142 is classified as having an unresolved core morphology [100] and features a clear peak in its radio spectrum at 0.9 GHz with moderate variability at 2.3 and 8.2 GHz. We suggest the quasar as a new MPS candidate, although the optically thick spectral part (${\alpha }_{\mathrm{thick}}=+0.25\pm 0.03$) is represented by only three radio measurements [140,144,145].

For the core–jet quasar J1301+1904 we calculated moderate and negligible variabilities at frequencies of 2.3, 4.7, 8.2, and 11.2 GHz (Table A2) and a peaked spectrum with ${\nu }_{\mathrm{peak},\mathrm{obs}}\sim 1.7$ GHz. The low-frequency spectral part is represented by only single measurement at 150 MHz (TGSS) with ${\alpha }_{\mathrm{thin}}=+0.13\pm 0.01$.

A new GPS candidate J1401+1513 has the optically thick spectral part close to flat, ${\alpha }_{\mathrm{thick}}=+0.05\pm 0.01$, and the variability indices are 0.09 and 0.36 at 11.2 and 4.7 GHz ($t_{4.7}=33$ years). The low-frequency spectral part, poorly represented by a single TGSS measurement, should be determined more precisely.

The PS blazar [115] J1405+0415 has a compact core and a resolved jet extending to 18 mas [146,147,148]. Its radio spectrum is composed of several components: steep below 1 GHz, curved with a peak between 1 and 10 GHz, steep from 10 to 110 GHz, and inverted above 110 GHz. The quasar has moderate variability indices ranging from 0.19 to 0.26 at the radio frequencies.

A new MPS candidate J1418+4250 has ${\nu }_{\mathrm{peak},\mathrm{obs}}\sim 0.9$ GHz, two measurements below the peak (325 and 408 MHz), and measurements at five frequencies above it. More long-term observations would be complementary, both for peak frequency clarification and variability estimation. To estimate the variability for the time being, we used a frequency of 4.7 GHz and obtained $V_{S_{4.7}}=0.02$ on a timescale of 34 years, based however only on three measurements of Green Bank in April and October 1987 (4.85 GHz) and RATAN-600 in July 2020 (4.7 GHz). The flux density varies insignificantly from 0.047 and 0.048 Jy to 0.054 Jy, which is within the measurement uncertainties (0.006–0.009 Jy).

The new GPS candidate J1521+1756 has a spectral peak at 3.5 GHz, the variability indices are 0.10, 0.18, and 0.32 at 11.2, 8.2, and 4.7 GHz on a timescale $t_{4.7}=40$ years (Table A2). The VLA images at 1.4 GHz show an unresolved core morphology [100].

The quasar J2019+1127 should be treated carefully because it is associated with the gravitationally lensed system MG B2016+112 consisting of three components [149]: a doubly imaged core and a counter–jet quadruply imaged at higher resolution [150]. The lensed radio images of the object can not be resolved by RATAN-600. Its radio spectrum has a complex structure below the peak frequency, caused either by the several emission components or by the variability of the source.

J2314+0201 is a new GPS source of X-ray emission [151,152]. We estimated the variability indices $V_{S_{4.7-11.2}}=0.14$–$0.20$. Its radio spectrum has a peak, however there is a lack of measurements in the optically thick spectral part, and we rely only on one measurement point at 150 MHz from the TGSS.

Appendix C.3. Gravitationally Lensed Quasars

The lensing galaxy at $z=1.145$[153] creates four images of the quasar J0131+4358 with a maximum image separation of $0.{54}^{ }$[154]. The quasar J0751+2716 is lensed by a galaxy at $z=0.35$[155] into four images separated by $0.8^{ }$[156]. J1401+1513 is a unique sextuple system with a maximum angular separation between the quasar images of about $1.7^{″}$, the lensing group of three galaxies are located at $z\sim 1$ [157]. The quadruple images of J1424+2256, created by a lensing galaxy at $z=0.34$[158], have a maximum separation of $1.3^{ }$[159,160]. The lensing galaxy located at $z=1.165$[161] creates a triple image of J1632−0033 separated by about $1.5^{ }$, according to the VLA observations at 43 GHz [162,163]. J2019+1127 is a narrow-line Seyfert 1 galaxy, reviewed in Section ?? as a new MPS candidate. The separation between three of its components is $3.4^{″}$ [149]. Based on our results, the radio variability properties of the lensed quasars show no peculiarities compared to other quasars.