© The Authors 2024

1. Introduction

Peaked spectrum (PS) sources are sources that show a peak in their radio spectrum and are observed to be compact. They have been hypothesized to be the young progenitors of radio-loud active galactic nuclei (AGN). These sources can be found with a range of intrinsic turnover frequencies (ν_rest, turn). PS sources are typically classified into the sub-classes: gigahertz-peaked spectrum (GPS) sources, high-frequency peaked (HFP) sources, compact steep-spectrum (CSS) sources, and megahertz-peaked spectrum (MPS) sources. These sub-classes differ in terms of linear size and turnover frequency, but are likely part of the continuum of radio galaxies, with HFP sources having the smallest and MPS sources having the largest linear sizes. A correlation between the projected linear size l in kpc, and the turnover frequency in GHz of these sources has been identified by O’Dea (1998) and is given by ν_rest, turn ∝ l^−0.65. Small PS sources thus have higher turnover frequencies.

HFP sources are sources with a spectral turnover above 5 GHz (O’Dea & Saikia 2021), and are the most compact PS sources with a linear size less than 500 pc (Dallacasa et al. 2001). GPS radio sources have spectral peaks between 500 MHz and 5 GHz (Gopal-Krishna 1983). They have typical linear sizes ≲1 kpc, and are powerful, with log P_1.4 GHz ≳ 25 W Hz⁻¹ (O’Dea 1998). CSS sources are just as powerful as GPS sources, with larger linear sizes (1–20 kpc; Fanti et al. 1990), and spectral peaks at < 500 MHz. MPS sources are sources that peak below 1 GHz in the observer’s frame. The MPS population is likely a combination of relatively nearby CSS and GPS sources, and compact high-frequency peaked sources at higher redshift, whose turnover frequency has been shifted to low frequencies due to the cosmological redshift (O’Dea & Saikia 2021).

There are two main hypotheses as to what produces the spectral turnovers and small linear sizes of PS sources. One of these hypotheses is that PS sources are the young precursors (e.g., Phillips & Mutel 1980; Wilkinson et al. 1994) of Fanaroff-Riley I and II (FRI and FRII; Fanaroff & Riley 1974) galaxies. PS sources often show evidence of compact double-lobed structures (O’Dea 1998). Furthermore, a relationship between radio power and linear size was found for PS sources (e.g., Kunert-Bajraszewska et al. 2010; An & Baan 2012). From the relation between linear size and turnover frequency, we find that if the youth model applies the majority of the time, we would expect HFP sources to be the youngest PS sources, evolving into GPS sources, then into MPS sources, and finally into FRI and FRII sources (e.g., Carvalho 1985; Kunert-Bajraszewska et al. 2010). Evidence from spectral break modeling (Callingham et al. 2015) and the motions of lobe hot spots (Owsianik & Conway 1998; Kaiser & Best 2007) supports this theory. In the youth hypothesis, synchrotron self-absorption (SSA; Snellen et al. 2000; de Vries et al. 2009) is the cause of the spectral turnover.

The other hypothesis is that the small linear sizes and spectral turnovers of PS sources are caused by their radio jets being contained within dense circumnuclear environments. Due to the dense environments, sources get “frustrated”, and they are not able to grow to larger spatial scales. The absorption mechanism associated with the frustration hypothesis is free-free absorption (FFA; Bicknell et al. 1997; Peck 1999; Callingham et al. 2015). Evidence for the frustration hypothesis is that the radio morphologies of CSS sources indicate strong interactions between the jets of a source and their environments (Wilkinson et al. 1984; van Breugel et al. 1984; Kunert-Bajraszewska et al. 2010). Furthermore, extended emission around PS sources was observed, indicating multiple epochs of activity (Baum et al. 1990; Stanghellini et al. 1990). In individual PS sources, unusually high densities were found, which also indicates FFA emission would be dominant (Peck 1999; Callingham et al. 2015; Sobolewska et al. 2019).

Research from the past decades has shown that the majority of PS sources are young rather than frustrated (e.g., Carvalho 1985; Kunert-Bajraszewska et al. 2010; O’Dea & Saikia 2021; Slob et al. 2022). It is always possible that a source is both young and frustrated, and individual PS sources could be associated with either or both mechanisms. In this work, we will assume that the youth hypothesis is dominant over the frustration hypothesis.

Most research that was done on PS sources in the last years (e.g., Snellen et al. 1998; Callingham et al. 2017; Keim et al. 2019; Slob et al. 2022) has focused on determining the cause of the spectral turnovers and the small linear sizes of PS sources. None of these studies compared the different types of PS sources against each other, which is vital in determining the evolution of sources from HFP sources to FRI and FRII sources. By studying the abundances of the different types of PS sources, we can determine the relative lifetimes of these phases and understand how the population of large-scale radio sources has evolved. Such an analysis requires statistically large samples of MPS, GPS, and HFP sources.

In the past decades, a revolution has taken place in wide-field radio astronomy. Surveys of unparalleled sensitivity, resolution, and sky area have been made. In this work, we will make use of all-sky surveys from the LOw Frequency ARray (LOFAR; van Haarlem et al. 2013) and the Karl G. Jansky Very Large Array (VLA). From LOFAR, we use the LOFAR Two-metre Sky Survey (LoTSS; Shimwell et al. 2022) and the LOFAR LBA Sky Survey (LoLSS; de Gasperin et al. 2023). Using these surveys, we can determine the low-frequency part of the spectral energy distribution (SED) of more sources than ever before. Constraining the low-frequency region of SEDs is vital in determining spectral turnovers and thus in identifying PS sources, specifically in identifying MPS sources, that have lower-frequency turnovers. The Very Large Array Sky Survey (VLASS; Lacy et al. 2020), which is a higher frequency survey, allows us to extend this research to GPS sources with turnovers at higher frequencies. Now that such surveys are available, we can construct larger samples of MPS and GPS sources than ever before, and we can compare these samples with each other. From these samples, we can find the relative abundances of these sources in the radio sky, and determine their relative lifetimes, with statistical robustness.

The surveys used in this work are further introduced in Sect. 2. They allow us to compile master samples of unresolved isolated radio sources, as will be described in Sect. 3. In Sect. 4, we use these master samples to identify samples of MPS and GPS sources. In Sect. 5 we construct Euclidean normalized source counts for our PS samples, which helps us understand the abundances of PS sources in the radio sky and allows us to draw conclusions about their relative lifetimes. In this paper we adopt H₀  =  70 km s⁻¹ Mpc⁻¹, Ω_M  =  0.27 and Ω_Λ  =  0.73 for a standard Lambda cold dark matter cosmological model (Hinshaw et al. 2013).

2. Surveys

To select our samples of PS sources, we used LoLSS, LoTSS, VLASS, and the NRAO VLA Sky Survey (NVSS; Condon et al. 1998). These surveys are among the most sensitive wide-field radio surveys to date, which allowed us to select faint PS sources. The sensitivity of common wide-field radio surveys is illustrated in Figure 1. We indicate the lowest reported flux density for many wide-field radio surveys from the last three decades. We also plot the faintest typical SSA PS sources that could be observed in the PS samples from Slob et al. (2022) and Callingham et al. (2017), as well as the faintest PS samples that could be observed in our GPS and MPS samples. We find that our MPS sample reaches almost four times lower flux densities than Slob et al. (2022), and our GPS sample reaches a 40 times lower flux density than Callingham et al. (2017). We can determine from Figure 1 that the sensitivity of LoLSS limits the faintest MPS sources we can identify, while the sensitivity of NVSS limits the faintest GPS sources selected.

Fig. 1. Observational limits of different radio surveys, as well as literary PS samples which are indicated at their central frequency. The lowest cataloged component flux density reported for each survey used in this research, and other relevant large radio surveys, is indicated in the plot. The GaLactic and Extragalactic All-Sky MWA Survey (GLEAM; Wayth et al. 2015) is represented by a line since it has variable limiting flux densities at different observing frequencies. The Callingham et al. (2017) sample is represented by an SSA model with a faintest peak flux density of 160 mJy at a turnover frequency of 190 MHz. The Slob et al. (2022) sample is represented by an SSA model with a faintest peak flux density 80 mJy at a turnover frequency of 100 MHz. We indicate the sensitivity of our samples with SSA models with a limiting flux density of 23 mJy at a turnover frequency of 100 MHz for the MPS sample, and with a limiting flux density of 3.9 mJy at 750 MHz for the GPS sample. In the figure, multiple wide-field radio surveys were indicated that have not been mentioned before; the Westerbork Northern Sky Survey (WENSS; Rengelink et al. 1997), Sydney University Molonglo Sky Survey (SUMSS; Mauch et al. 2003), Cambridge 7C (Hales et al. 2007), Australia Telescope 20 GHz (AT20G; Murphy et al. 2010), Very Large Array Low-Frequency Sky Survey Redux (VLSSr; Lane et al. 2014), TIFR GMRT Sky Survey Alternative Data Release 1 (TGSS; Intema et al. 2017), and the RAPID ASKAP Continuum Survey (RACS; Hale et al. 2021).

2.1. LOFAR surveys

LOFAR (van Haarlem et al. 2013) is a radio telescope that consists of 52 stations, centered in the Netherlands. The stations consist of high-band (HBA, 110–250 MHz) and low-band (LBA, 10–90 MHz) antennae.

LoTSS (Shimwell et al. 2022) is a deep wide-field radio survey, performed with the LOFAR HBA at 120–168 MHz. It aims to observe the whole northern sky. Data from the second data release (DR2) was used in this work, which consists of observations from 27% of the northern sky. The coverage is split into two regions centered at approximately 12h45m +44° 30′ and 1h00m +2° 00′, spanning 4178 deg² and 1457 deg², respectively. The DR2 catalog contains 4 396 228 radio sources, derived from the total intensity maps. It reaches a median rms sensitivity of 83 μJy beam⁻¹, a resolution of 6″, and is 90% complete at 0.8 mJy beam⁻¹.

In-band spectra are available for the sources in LoTSS, at 128, 144, and 160 MHz. The in-band spectra provide information about spectral properties along the 48 MHz wide band. However, these in-band spectra are not reliable for most sources, due to the narrow frequency range, combined with non-negligible uncertainties in the alignment of the flux density scale of the in-band images. Therefore, the in-band spectra were only used in this work as a visual guide to confirm spectral fits and were not used to select PS sources.

LoLSS (de Gasperin et al. 2023) is a low-band LOFAR wide-area survey at 42–66 MHz, centered at 54 MHz, which eventually will cover the whole sky above a declination of 24°. In this work, data from the first data release (DR1) was used, which has a typical rms of 1.55 mJy beam⁻¹, four times lower relative to the LoLSS preliminary data release (PDR; de Gasperin et al. 2021). So far, 650 deg² sky area has been released, centered around the Hobby-Eberly Telescope Dark Energy Experiment (HETDEX; Hill et al. 2008) Spring Field with 11 < RA < 16h and 45 < Dec < 62°. The LoLSS DR1 catalog consists of 42 463 sources, with an angular resolution of 15″, and 95% completeness at 11 mJy.

Separate LoLSS in-band spectra, at 44, 48, 52, 56, 60, and 64 MHz, are available. The in-band spectra were only used for visual confirmation that spectral fits seem reliable, and were not used for selecting PS sources.

2.2. VLA surveys

NVSS (Condon et al. 1998) is a 1.4 GHz sky survey performed with the VLA, covering the sky north of a declination of −40° (82% of the celestial sphere). Observations were done between September 1993 and October 1996. The NVSS catalog consists of 1 773 484 sources, has a resolution 45″, and is 99% complete at 3.4 mJy.

VLASS (Lacy et al. 2020) is a wide-field radio survey at 2–4 GHz, centered at 3 GHz, with an angular resolution of 2.5″, covering the whole sky above a declination of −40° (33 885 deg²). It has an rms goal of σ₁ = 70 μJy/beam in the coadded data, while the sensitivity for a single epoch is 120 μJy/beam. The B configuration of the VLA was used, with a maximum antenna separation of 10 km.

VLASS will observe the entire sky three times, which makes it possible to study variable and transient sources. By the time of writing, the first and second epochs were processed, and Quick look (QL) images for these epochs were released. The rapid CLEANing limits the quality of these images, but they are reliable for our science. In our research, the component catalog of epoch 2 was used. We only use the second epoch because the rapid CLEAN algorithm used was updated between epochs 1.1 and 1.2. There are several issues with epoch 1.1, which are related to a systematic under-measurement of flux density value in the QL images, issues with astrometry, and ghost images of bright sources. These issues were solved for the second epoch, so we disregard epoch 1 (Gordon et al. 2023). The VLASS QL epoch 2 catalog consists of 2 995 271 components.

3. Sample selection

In this work, we select two master samples, the low-frequency master sample (LF) and the high-frequency master sample (HF), defined by combining radio surveys and making cuts based on source confusion and resolution. In Sect. 4 we describe how we selected two sub-samples of PS sources from the master samples based on the shape of their SEDs. An overview of each step in the sample selection is provided in Table 1, which includes the number of sources in our samples after each selection step.

3.1. Selecting the LoTSS+NVSS sample

We started the selection of the master samples with the LoTSS DR2 catalog. To ensure that source confusion did not impact the derived spectra, we made a selection based on whether the sources were isolated. Of all surveys used in this research, NVSS has the lowest resolution of 45″. Therefore we only selected sources in LoTSS if there are no other sources within a radius of 45″.

Deconvolution errors around very bright (≳200 mJy) LoTSS sources can produce artifacts that can cause a source to appear as not isolated. To include these bright sources that are surrounded by complex noise structures, we included a source in the master sample if all neighboring sources had a total flux density less than 10% of the brightest source. The selection defined here contains 2 775 395 LoTSS sources.

We expect all PS sources to be unresolved in LoTSS since they generally are compact (≲1″ in linear size, O’Dea 1998; Chhetri et al. 2018), and LoTSS has a resolution of 6″. Therefore, all resolved sources in LoTSS were removed from the sample. We used the criteria for selecting resolved sources in LoTSS that were described by Shimwell et al. (2022). They defined an envelope that encompasses the 99.9 percentile of the R distribution, where $R=ln(\frac{S_I}{S_P})$ describes the ratio between the integrated flux density S_I and the peak flux density S_P. A source was defined as being resolved if R is larger than or equal to

$$\begin{array}{c}\hfill R_{99.9}=0.42+\left(\frac{1.08}{1+{\left(\frac{S/N}{96.57}\right)}^{2.49}}\right).\end{array}$$(1)

Here, S/N is the signal-to-noise ratio defined as $\frac{S_I}{{\sigma }_I}$, with σ_I the uncertainty on the integrated flux density. All sources with R > R_99.9 were deemed resolved, so we removed them from the sample. This step removed 6.8% of the sources in LoTSS, leaving us with 2 586 267 sources.

The above criterion has excluded the majority of the resolved sources, but it is based on the 99.9th percentile of the distribution. To remove resolved sources that might be left in the sample, we removed 10 sources with PyBDSF (Mohan & Rafferty 2015) S-code “C”, flagging complex sources, leaving 2 586 257 sources. The “C” sources correspond to multiple single-Gaussian sources in a single source-fitting island. By excluding these sources, only sources remain with S-code ‘S’, which are isolated sources fitted by a single Gaussian distribution, or S-code ‘M’, which are sources fit by multiple Gaussian distributions. The latter are kept in the sample because the fitting of multiple Gaussians might be caused by deconvolution errors.

The isolated, unresolved sources in LoTSS were then crossmatched to NVSS, with a crossmatching radius of 10″, which was found to be the best trade-off between completeness and reliability. The resulting LoTSS+NVSS sample contains 146 975 sources and forms the basis for both master samples defined in this work. Below, the LF and HF master samples are introduced separately.

3.2. Low-frequency sample

To ensure we only considered a single source in NVSS, we demanded in LoLSS DR1 there are no other sources within 45″, with the exception that sources ten times brighter than their neighbor are allowed in our sample. This criterion removed 1582 sources from the LoLSS sample, leaving us with 40 881 sources. We crossmatched the LoTSS+NVSS sample defined above to our selection of LoLSS, with a crossmatching radius of 5″, which produces our LF sample consisting of 12 962 sources. This crossmatching radius was chosen to be larger than the astrometric inaccuracy and smaller than the resolution of LoLSS (de Gasperin et al. 2023).

3.3. High-frequency master sample

We defined the HF sample by crossmatching the LoTSS+NVSS sample to the VLASS QL component catalog. The entire VLASS QL epoch 2 catalog contains 2 995 271 components and still contains duplicates due to the overlap of the QL images. As recommended by the user guide (Gordon et al. 2023), we selected only the sources without a duplicate, or sources with the best signal-to-noise ratio among multiple duplicate sources. We also select only sources with quality flags 0 or 4, according to their recommendations. This selection was designed to limit the contamination of spurious detections stemming from the limited quality of the QL images. After making these cuts to obtain reliable flux density measurements, we were left with a VLASS catalog containing 2 446 020 components.

We also selected only sources in VLASS that have no neighbors within a radius of 45″. However, if the total flux density of an unisolated source is at least ten times brighter than their neighbor, it would be included in the sample, since these neighbors are likely spuriously identified sources caused by sidelobe structures. By removing these clustered sources, we are left with 1 384 792 sources.

After removing the clustered sources, we also removed any source with a PyBDSF S-code “C”, which removed 13 559 sources from our sample, leaving us with only “S” or “M” S-code sources. Our final VLASS catalog contains 1 371 233 sources.

We then crossmatched this sample to the LoTSS+NVSS sample with a crossmatching radius of 2.5″, larger than the astrometric accuracy of VLASS (0.5–1″; Gordon et al. 2023). VLASS is not very sensitive to extended components due to the lack of short baselines of the VLA, while LoTSS is sensitive to extended lobe emission due to the abundance of short baselines of LOFAR. Furthermore, the core of a radio galaxy generally has a flatter spectrum than its lobes, causing the core to dominate over the lobes at higher frequencies. These effects might cause a shift in peak position larger than the astrometric inaccuracy, therefore a crossmatching radius larger than the astrometric accuracy of VLASS was chosen. We found 108 473 sources in this crossmatched sample, which we refer to as the HF master sample.

Because of the overlap in survey areas, it was expected that approximately all sources in the LoTSS+NVSS sample have a match in VLASS. However, we found that only 73% of the sources in the LoTSS+NVSS do. If the isolation criterion was not applied, this number would be 86%. Thus 13% of the sources without a VLASS match are likely isolated and unresolved at 6″ in LoTSS, and resolved in VLASS at 2.5″, which causes them to be removed from the sample.

We accounted for the other 14% of sources in LoTSS+NVSS that do not have a VLASS detection via several other factors. We found that 0.4% of the sources in the LoTSS+NVSS sample have missing VLASS detections because of missing VLASS coverage. 3.6% of the LoTSS+NVSS sources are too faint to be detected in VLASS, or VLASS has unusually high noise in that region of the sky.

For the remaining 10% of sources with a missing VLASS detection, we found some were resolved in LoTSS at 6″ while still passing our source size selection criterion. At the VLASS resolution of 2.5″, these sources are then strongly resolved. VLASS is more sensitive to compact structures, and radio lobes tend to have steeper spectra than the core. Therefore VLASS has more difficulty in detecting resolved lobes than LoTSS, and a large fraction of these sources are resolved out and go undetected in VLASS, which causes no VLASS flux density measurement to be available for these sources.

We investigated how these resolved sources got incorrectly labeled as unresolved in our selection process. Some look like standard double-lobed radio sources in LoTSS, where the lobes get fit as 2 separate components by PyBDSF. Since these components are further apart than 45″, and they pass the unresolved criteria from Eq. (1), they are not flagged as resolved sources. Others were only marginally resolved and therefore passed the selection criteria.

We have thus confirmed that the majority of sources with missing VLASS detections are resolved out in VLASS. This is not an issue for our science goal of identifying PS sources since the small angular sizes of almost all PS sources imply they must be unresolved in VLASS. We therefore assume no PS sources were removed from the LoTSS+NVSS sample due to missing VLASS detections.

3.4. Crossmatching to additional radio surveys

We crossmatched additional radio surveys to our LF and HF master samples, as well as the LoTSS and LoLSS inband spectra, using a crossmatching radius of 5″. The additional radio surveys are the TIFR GMRT Sky Survey Alternative Data Release 1 (TGSS; Intema et al. 2017) at 150 MHz, the Faint Images of the Radio Sky at Twenty-centimeters (FIRST; Becker et al. 1995 at 1.4 GHz, the Westerbork Northern Sky Survey (WENSS; Rengelink et al. 1997) at 326 MHz, and the Very Large Array Low-Frequency Sky Survey Redux (VLSSr; Lane et al. 2014) at 74 MHz. These surveys were not used for spectral index fitting, but were only used for the visual confirmation of the SEDs of the sources in our sample.

3.5. Uncertainties in flux density measurements

In this work, uncertainties are dominated by systematic uncertainties between surveys, rather than individual uncertainties on sources. Therefore, we combined in quadrature the uncertainty on flux density measurements with 10% of the flux density in LoTSS, NVSS, and VLASS, to account for the systematics between surveys. For sources in LoLSS, we used the reported flux density uncertainties (precision 6%; de Gasperin et al. 2023).

4. Peaked spectrum sources

From the LF and HF master samples, we identified PS sources from the shape of their SED. The SED between two flux density measurements can be described by a power law,

$$\begin{array}{c}\hfill S=a{\nu }^{\alpha },\end{array}$$(2)

where α is the spectral index, and a is the normalization parameter.

In the LF master sample, we defined α_low to be the spectral index between LoLSS and LoTSS, and α_high to be the spectral index between LoTSS and NVSS. In the HF master sample, we defined α_low to be the spectral index between LoTSS and NVSS, and α_high to be the spectral index between NVSS and VLASS.

Since we derived a and α from a fit to only two flux density measurements, the covariance matrix of the fit can not provide an uncertainty on the fit. To find the uncertainty σ of the spectral index, we used error-in-variables regression.

We defined a source to be a PS source if (α_low > σ_{α_low}) and (α_high < 0). These criteria ensured that α_low is always larger than zero with at least a 1σ certainty. We were less strict for α_high since any source with a positive α_low must have a turnover at some higher frequency, and if α_high < 0, we were sure that the turnover occurs below 1400 MHz for the MPS sample and below 3000 MHz for the GPS sample.

4.1. Megahertz-peaked spectrum sample

We present the color-color plot for the LF sample in Figure 2. The lower right quadrant contains PS sources, plotted in red. In the LF master sample of 12 962 sources, we find 506 PS sources (3.9%), which we will refer to as the megahertz-peaked spectrum (MPS) sample. The MPS sample is the largest sample of MPS sources identified to date, 1.4 times larger than the sample by Slob et al. (2022). The MPS sources have a spectral peak approximately at 144 MHz.

Fig. 2. Color-color plot for the 12 962 sources in our MPS sample, where αlow is the spectral index between LoLSS (54 MHz) and LoTSS (144 MHz), and αhigh is the spectral index between LoTSS and NVSS (1400 MHz). PS sources exist in the lower right quadrant and are indicated in red. The blue line represents the median of the error of αlow at 0.18. The contours represent 17, 60, 168, and 308 sources. It should be noted that there might be PS that are not indicated in red because they are hidden by the contours. The red dashed line represents the 1:1 ratio of the spectral indices. The median error bars are plotted instead of the individual errors. The normalized distributions of αlow and αhigh are plotted, with a median and standard deviation of −0.6 ± 0.4 respectively −0.7 ± 0.2. Mock SEDs are shown, to indicate the rough shape the SED of a source in each quadrant would have.

The method used so far for selecting the MPS sample is similar to that employed by Slob et al. (2022), where the main difference between our work and Slob et al. is that they used the LoLSS preliminary data release (PDR), while in this work we used Data Release 1 (DR1). For reference, the Slob et al. (2022) source counts values are indicated in Table B.2. DR1 is more sensitive and has a higher resolution (15″) than the PDR (47″) since direction-dependent effects have been corrected. We do not expect that the change in resolution has a significant effect on the classification of PS sources since the LoTSS+NVSS sample used in both of our works only includes isolated sources in 45″. Furthermore, the median uncertainty of α_low for the LF sample is 0.18, therefore our selection criteria to determine when a source is PS is stricter than that used by Slob et al., who defined a source to be PS if α_low > 0.1 and α_high < 0.

If we were to use the same criteria as Slob et al. to define when a source is a PS source, we would find 728 MPS sources in our LF sample. One would expect to find that most, if not all, of the 373 PS sources in the sample by Slob et al. are also present in this sample of 728 MPS sources. After accounting for sources that do not have a LoLSS observation in both the PDR and DR1, we find 148 sources were identified as PS in the work by Slob et al. and are not identified as MPS sources in our work. Conversely, there are 70 sources that we identify as PS sources in our MPS sample, that were not identified in Slob et al. as PS sources.

We find that this discrepancy is caused by a significant inconsistency in the reported flux densities for individual sources between the PDR and DR1. This inconsistency causes α_low to be different for sources in our sample and in the sample by Slob et al. As was outlined by de Gasperin et al. (2023), the inconsistency can be explained by systematic effects dominating over the noise in the PDR. These effects cause the flux density of individual sources to vary significantly between the PDR and DR1, which causes sources to be mislabeled in the Slob et al. sample. The PDR should thus be used with care. We assumed DR1 is more accurate and we will continue working with this survey in the rest of this work.

4.2. Gigahertz-peaked spectrum sample

In Figure 3 we plot the color-color plot for the HF master sample. We used the same selection criteria as for the MPS sample to define the gigahertz-peaked spectrum (GPS) sample. In the HF master sample of 108 473 sources, we identified 8032 GPS sources, which corresponds to 7.4%. These sources peak around 1400 GHz, and are indicated in red. The GPS sample is the largest sample of PS sources to date, over five times larger than the sample isolated by Callingham et al. (2017).

Fig. 3. Color-color plot for the 108 473 sources in our GPS sample, where αlow is the spectral index between LoTSS (144 MHz) and NVSS (1400 MHz), and αhigh is the spectral index between NVSS and VLASS (3 GHz). PS sources exist in the lower right quadrant and are indicated in red. The blue line represents the median of the error of αlow at 0.10, which corresponds to the median of the selection limit of PS sources. The contours represent 134, 416, 1349, and 2800 sources. The red dashed line represents the 1:1 ratio of the spectral indices. The median error bars are plotted instead of the individual errors, for readability. The normalized distributions of αlow and αhigh are plotted, with a median and standard deviation of −0.6 ± 0.4 respectively −0.7 ± 0.5. Mock SEDs are shown, to indicate the rough shape the SED of a source in each quadrant would have.

4.3. Comparison to literature PS samples

To confirm the validity of our selection process, we compared our MPS and GPS samples to previously identified CSS, GPS, CSO, and HFP samples (O’Dea 1998; Snellen et al. 1998, 2000; Peck & Taylor 2000a; Tinti et al. 2005; Labiano et al. 2007; Edwards & Tingay 2004; Randall et al. 2011). 45 of the literature sources were in our HF sample, 13 of which we identified as GPS sources. There are four reasons the remaining 32 were not identified as PS sources. Firstly, 13 of those sources are CSS sources, which do not have peaks in the frequency range of our surveys. Secondly, six sources had positive spectral indices over the whole frequency range of our surveys, indicating they have a turnover at frequencies ≳3 GHz. Thirdly, seven sources had an uncertainty on α_low that was too large, which indicates that the spectral peak is abnormally broad or occurs near the low end of the frequency window. Finally, six sources showed significant variability in their SEDs, indicating they are likely blazars. In light of this, we conclude that our selection process is valid for our GPS sample.

Of the literature PS sources, three were in our LF sample, which we did not identify as MPS sources. These three sources are also in the HF sample, therefore we will not discuss them separately here.

To highlight the diversity of PS sources we identified, we draw attention to the source B1315+415 in Figure 4. We fit the curve from Eq. (3) to the SED of this source, and find a peak frequency of 1.9 ± 0.5 GHz, and spectral indices α_low = 0.7 ± 0.1 and α_high = −0.58 ± 0.08. B1315+415 was identified as a PS source by Labiano et al. (2007), with a turnover at 2.3 GHz. The LoTSS image shows significant extended emission (extending to roughly 200 kpc) around a compact core with a peaked spectrum. We interpret the combination of extended structure around a PS core as evidence this source has restarted activity. This source, also indicated as J131739.21+411545.6, was also identified as a potential restarted PS source by Kukreti et al. (2023). They identify eight candidate sources that only show extended emission in their LoTSS images, and no extended emission at higher frequencies, which suggests that the extended emission is old. These candidates seem to have a broader [O III] profile than the non-peaked sources, suggesting that these candidate restarted AGN have more disturbed gas kinematics than evolved radio AGN.

Fig. 4. LoTSS image and SED of B1315+415, which is likely a restarted PS source. Top: The LoTSS view of B1315+415. The LoTSS beam size is indicated in the upper left corner, and the color bar indicates the flux density per beam. The lobes extend to roughly 200 kpc. Bottom: The SED of B1315+415, which contains flux density measurements from the B3 survey (Ficarra et al. 1985), the Very Long Baseline Array (VLBA) Calibrator Survey (VCS; Kovalev et al. 2007), the NRAO Green Bank telescope (Becker et al. 1991), the VLBA Imaging and Polarimetry Survey (VIPS; Helmboldt et al. 2007), the VLA (Patnaik et al. 1992), and the Effelsberg telescope (Vitale et al. 2015).

4.4. Sources with extreme spectra

In our LF sample, we find three PS sources with α_low > 2.5, and in the HF sample, we find one such source. These sources are interesting because such steep spectral indices imply that SSA can not be responsible for the turnover (O’Dea 1998). Therefore, these sources may be more consistent with the frustration hypothesis of PS sources. The SEDs of these extreme sources can be found in Figure 5. We indicate the exact spectral indices in the caption.

Fig. 5. SEDs of the four sources in our sample with extremely steep spectral indices. The curve from Eq. (3) was fitted to the radio flux density measurements available for these sources and is indicated as a black curve. Upper left: ILTJ144448.65+560615.6, which has a spectral index between LoLSS and LoTSS of αlow = 2.5 ± 0.4. Upper right: ILTJ151214.28+542710.1 with a spectral index between LoLSS and LoTSS of αlow = 2.9 ± 0.5. Lower left: ILTJ132201.76+564652.5 with a spectral index between LoLSS and LoTSS of αlow = 2.5 ± 0.3. Lower Right: 87GB B0903+6656 (ILTJ090723.44+664446.6), with a spectral index between LoTSS and NVSS of αlow = 2.7 ± 0.1. For this source, we additionally plot a 4.58 GHz observation taken by the NRAO Green Bank telescope (Becker et al. 1991), and an 8.4 GHz observation taken by the VLA (Healey et al. 2007).

To the flux density measurements available for these sources, we fit a generic curved model (Snellen et al. 1998):

$$\begin{array}{c}\hfill S(\nu )=\frac{S_{\mathrm{peak}}}{(1-e^{-1})}{\left(\frac{\nu }{{\nu }_{\mathrm{peak}}}\right)}^{{\alpha }_{\mathrm{thick}}}(1-exp{(-\nu /{\nu }_{\mathrm{peak}})}^{{\alpha }_{\mathrm{thin}}-{\alpha }_{\mathrm{thick}}}),\end{array}$$(3)

where S(ν) is the flux density at frequency ν in MHz, α_thin and α_thick indicate the spectral indices in the optically thin and optically thick regions of the SED, respectively, and S_peak is the flux density at the peak frequency ν_peak.

We want to particularly highlight 87GB B0903+6656 as its measurement of α_low > 2.5 is significant. Healey et al. (2007) previously identified this source as a flat-spectrum radio source. The extreme turnover is consistent with a non-detection of the source at 408 MHz (Marecki et al. 1999). The host of this source has been identified as a quasar by e.g. Souchay et al. (2009). Followup HI absorption or X-ray observations are warranted to confirm if the inferred dense circumnuclear medium surrounds the core-jet structure of the source.

4.5. Other sources with interesting spectra

We note that the upper right quadrant of the color-color plots contains sources with a positive spectral index across the entire frequency range probed by our surveys, therefore a spectral turnover must occur at some frequency higher than our observing window. There are 95 sources in the upper right quadrant of Figure 2 of the LF sample. These sources are likely either GPS or HFP sources, peaking above ∼1.4 GHz. Similarly, the upper right quadrant of Figure 3 of the HF sample contains 3015 sources, which are most likely HFP sources as they must peak above ∼2 GHz. Out of the 95 sources in the upper right quadrant of the LF sample, 54 are indeed classified as GPS sources in our sample.

Sources in the upper left quadrants of the color-color plots exhibit a convex spectrum, which is likely indicative of multiple epochs of AGN radio activity (Callingham et al. 2017). The SED of a source that shows such a convex spectrum is presented in Figure 6. Sources like this one are composite sources, with a steep-spectrum power-law component at low frequencies and an inverted component at high frequencies. In the LF sample, there are 134 convex spectrum sources, and in the HF sample, there are 5610. Follow-up observations at higher frequencies (≳5 GHz) are required to trace the higher frequency turnover and ensure the spectrum is not a product of intrinsic variability.

Fig. 6. SED of NGC 3894 (ILTJ114850.36+592456.2), which is an example of a source with a convex SED. Such an SED might be caused by multiple epochs of activity. Additional high-frequency measurements were added; a 5 Ghz NRAO 300-ft telescope measurement (Sramek 1975), a 10.45 GHz 100 m Effelsberg telescope measurement (Pasetto et al. 2016), and two VLBI measurements at 8 GHz and 415 GHz (Tremblay et al. 2016).

To the SED in Figure 6 we fit Eq. (3), with an additional power law term added to constrain the uptick. We find values of α_low = 0.18 ± 0.09 and α_high = −1.8 ± 0.9, and a spectral index of α = −2 ± 1 for the additional power law term.

This source, NGC 3894, was identified as a nearby low-power CSO by Peck & Taylor (2000b), Taylor et al. (1998). Tremblay et al. (2016) use Very Long Baseline Interferometry (VLBI) and found this source has a flat-spectrum core with diffuse lobes extending to the northwest and southeast. They find the core has a strongly inverted spectral index (α ∼ 1 − 1.5) between 5 and 8 GHz.

5. Euclidean normalized source counts

We studied the evolution of the PS sources in our samples via the differential source counts of the MPS sample at 144 MHz, and of the GPS sample at 1400 MHz. If PS sources are indeed young radio galaxies, there will be an offset between the source counts of a complete sample of radio-loud AGN and the source counts of our samples of PS sources, corresponding to the relative lifetimes of the PS sources.

5.1. MPS sample source counts

To compute the 144 MHz source counts for the MPS sample, we only included sources with LoTSS flux densities above S_144 MHz > 18 mJy. This limit corresponds to the 95% completeness limit of LoLSS (S_54 MHz = 11 mJy; de Gasperin et al. 2021), rescaled to 144 MHz assuming a power law with a spectral index of 0.52. The spectral index of 0.52 is the median of α_low in the MPS sample.

There are 484 MPS sources with a 144 MHz flux density brighter than 18 mJy. The resulting normalized source counts can be found in Figure 7 and in Table B.1. We removed the lowest flux density bin, as it was impacted by incompleteness.

Fig. 7. 144 MHz Euclidean normalized differential source counts. Top: The source counts of our MPS sample. Plotted are also the source counts from Callingham et al. (2017) and Slob et al. (2022), which are samples of PS sources. To compare them to a general population of AGN, the 150 MHz source counts from the LoTSS Deep Field (Mandal et al. 2021) are plotted. The 150 MHz source counts from the T-RECS simulation (Bonaldi et al. 2019) are plotted as well. The 150 MHz model of large-scale AGN from Massardi et al. (2010) is indicated by the solid gray line. This same model scaled down by a factor of 44 is indicated by the dashed gray line. Bottom: The residuals of our MPS source counts to the Massardi et al. (2010) 150 MHz model rescaled by a factor 44.

We compare our source counts with source counts from samples of PS sources from the literature. We plot the 144 MHz source counts of the MPS sample as obtained by Slob et al. (2022) and find they coincide with our measurements. We also plot the source counts from the PS sample of Callingham et al. (2017). Only sources with flux density S_143 MHz > 1 Jy were included, which corresponds to the reported 100% flux density completeness of the GLEAM survey. These source counts, though not extending to the same sensitivity, correspond well to the source counts we find for our sample in the higher signal-to-noise regime.

To compare the populations of PS sources with a general population of radio sources, we consider the source counts from the LoTSS Deep field (Mandal et al. 2021), which are the deepest 150 MHz source counts published to date. We find that the source counts of our MPS sample, as well as the other literary PS samples plotted, have a significant offset from the source counts of the AGN sample, while the shape of the distribution remains similar. This offset, with no shift in its peak, might indicate that MPS sources are indeed the young progenitors of FR I and FR II sources, with shorter lifetimes.

We plot the Massardi et al. (2010) model of the source counts of large-scale AGN, based on source counts compiled by De Zotti et al. (2009). In the models, sub-mJy sources were not taken into account as the contribution of star-forming galaxies becomes dominant over AGNs below roughly 1 mJy.

We further consider the source counts from the T-RECS simulation (Bonaldi et al. 2019) which simulates the radio sky in continuum, over the 150 MHz–20 GHz range. It includes polarization information for all radio sources, as well as realistic clustering properties by associating the sources to dark matter halos of a cosmological simulation. We can use the source counts of the total radio sky, as the AGN population dominates above 1 mJy (Massardi et al. 2010) which is the flux density range we are interested in. In Figure 7 we plot the 150 MHz T-RECS source counts.

We use orthogonal distance regression (ODR) to fit the 150 MHz Massardi et al. (2010) model, divided by a scaling parameter, to the source counts of our MPS sample. We find that the MPS sample has a best-fitting scaling parameter of 44 ± 2, which would indicate that MPS sources are 44 times less abundant than AGN at 144 Mhz. We plot the χ residuals, which are the residuals normalized by the uncertainties on the data, between our MPS source counts and the best fit for the scaled model in the bottom plot of Figure 7.

5.2. GPS sample source counts

We computed the 1400 MHz source counts for our GPS sample, which is limited by the completeness of NVSS (95% complete at S_1400 mHz = 3.24 mJy; Condon et al. 1998). We only included sources with NVSS flux density above this completeness limit. There are 6,924 sources in the GPS sample above the completeness limit, which we used to make the source counts. We removed the lowest two flux density bins, as they were impacted by incompleteness.

The Euclidean normalized source counts at 1400 MHz can be found in Figure 8 and in Table B.1. In this figure, we also plot the source counts for the sample of PS sources from Snellen et al. (1998). These source counts were evaluated at the individual peak frequencies for each PS source, corresponding to a median frequency of 2 GHz. We rescaled these source counts to 1400 MHz assuming the power law from Eq. (2) with a spectral index of –0.8.

Fig. 8. 1400 MHz Euclidean normalized differential source counts. Top: The source counts of our GPS sample. Plotted are also the source counts of the PS sample from Snellen et al. (1998). To compare our source counts to those of a general population of AGN, the 1400 MHz source counts from De Zotti et al. (2009) are plotted. The 1400 MHz source counts from the T-RECS simulation (Bonaldi et al. 2019) are indicated. The 1400 MHz model of large-scale AGN from Massardi et al. (2010) is indicated by the solid gray line. This same model scaled down by a factor of 28 is indicated by the dashed gray line. Bottom: residuals of our GPS source counts to the Massardi et al. (2010) 1400 MHz model rescaled by a factor 28.

To compare the source counts of our GPS sources with a general population of AGN, we plot the 1400 MHz model from Massardi et al. (2010), which was constructed using data presented by De Zotti et al. (2009). We further plot the source counts from the T-RECS simulation (Bonaldi et al. 2019).

We find our GPS samples have a significant offset from the AGN sample, while the shape of the distribution remains similar. The offset indicates the relative lifetimes of GPS sources. Through ODR fitting we find that the source counts of the GPS sample are scaled down by a factor 28 ± 1 compared to the AGN model from Massardi et al. (2010). From this offset we conclude that GPS sources are 28 times less abundant than radio AGN at 1400 MHz.

In Figure 8, a deviation at the higher flux density bins can be observed between the GPS sample and the Massardi et al. (2010) model. This is likely because of cosmic variance, caused by the limited sky area of our samples and the fact that sources with a higher flux density are more rare than sources with a lower flux density. Since LoTSS DR2 only observes 27% of the sky, we can expect this variance in the high flux-density bins.

5.3. Selection limits of the PS samples

We find that the source count offsets of our MPS and GPS samples (44 ± 2 and 28 ± 1 respectively) are significantly different from each other, which suggests that GPS sources are more abundant than MPS sources. However, we need to account for different selection effects between the two samples to ensure we are comparing equivalent distributions. We define a source to be a PS source if (α_low > σ_{α_low}) and (α_high < 0), therefore the distribution of σ_{α_low} influences how many sources will be identified as PS in our samples.

Sources in the LF sample generally have a larger σ_{α_low} (median 0.18) than sources in the HF sample (median 0.10). The distributions of σ_{α_low} are presented in Figure 9 for the HF and LF samples. As a consequence, a typical MPS source needs to have a steeper α_low than a typical GPS source to be reliably identified as a PS source. The effect will be less significant for high S/N sources since the uncertainty in the spectral index declines with increasing S/N.

Fig. 9. Distributions of σα_low for our LF and HF samples. The distribution of the LF sample is much broader than the distribution of the HF sample.

To guarantee the different σ_{α_low} distribution between the HF and LF samples is not the cause of the difference in the PS source count offsets, we needed to ensure both σ_{α_low} distributions are identical. We achieved that by randomly sampling new values of σ_{α_low} for the sources in our HF sample, from the distribution of σ_{α_low} of the LF sample. We accounted for the original σ_{α_low} dependence on S/N by binning the sources in S/N bins, such that high S/N sources would be assigned lower uncertainties than low S/N sources. With these new values of σ_{α_low} for the HF sample, we identified 5257 GPS sources, which is significantly lower than the 8032 GPS sources we identified with the original values of σ_{α_low}, consistent with the median uncertainty of σ_{α_low} of the HF sample roughly doubling.

We again constructed the source counts and find that the lowest flux density bins are affected more strongly than higher flux density bins, which is expected since the σ_{α_low} is larger for low S/N sources. After removing incomplete bins, we find the corrected GPS sample has an offset of 32 ± 2, which remains statistically different from the offset of 44 ± 2 that we found for the MPS sample.

To further test the robustness of the offset being different for our GPS and MPS samples, we also selected MPS and GPS sources using a hard limit of (α_low > 0.1) and (α_high < 0), as done in previous studies (e.g., Callingham et al. 2017; Slob et al. 2022). If we apply this same hard limit to our samples, we find an offset for the MPS sample of 39 ± 2 and an offset for the GPS sample of 29 ± 1. Using these selection criteria, the difference in the offset between the MPS and the GPS sample remains.

The reason we opted not to use this hard selection cut is that it does not account for the variable uncertainties in flux density measurements, which broadens the true distribution of the spectral indices. PS sources are located in the positive tail of the α_low distribution, therefore the broadening by uncertainty in α_low causes more sources to be shifted into the tail than out of it, which could cause an overestimation of the amount of PS sources. This is a fundamental problem with all samples of PS sources that are selected using a hard limit.

By using the S/N-dependent selection, the broadening of the spectral index distribution will have a smaller effect on how many PS sources are identified. When the broadening of the distribution of α_low occurs, the limit for a source to become a PS source also scales with σ_{α_low}. Therefore if the results from this work are compared with future work that has lower uncertainties in the flux density measurements, we would expect to find the same source counts offset, as long as the S/N-dependent selection limit is used.

A further robustness test for the difference in the offset of MPS and GPS sources would require simulating the noise properties of all the different surveys, and doing injection tests of sources with known spectra. Such a simulation would allow us to correct for any contamination of non-PS sources from the bulk of the spectral index distribution, which obviously would be worse for lower signal-to-noise sources.

However, several aspects make such an analysis challenging to conduct. Firstly, each of the four surveys used in this study has a different point-source response due to differing uv-coverage. Secondly, each survey would need to be searched using their corresponding survey source-finding algorithms, which would be particularly challenging for older surveys such as NVSS. Finally, “true” distributions for the parent LF and HF spectral indices would have to be assumed. There is no standardized spectral index distribution that is accepted as a true reflection of spectral indices from 54 MHz to 3 GHz, independent of signal-to-noise.

To conclude, we find that the offset of the MPS sample is larger than the offset of the GPS sample and is robust to different selection criteria. However, the exact offsets we find vary with different selection limits, which is a direct consequence of the tension between prioritizing either the completeness or the reliability of the isolated PS samples.

6. Discussion

We determine the abundance of MPS and GPS sources in the radio sky using the source count offsets we calculated in Sect. 5. We find that 2.22 ± 0.09% of the radio source population at 144 MHz are MPS sources, and 3.4 ± 0.1% of the radio source population at 1400 MHz are GPS sources. These percentages were calculated from the offset between the PS source counts and AGN population source counts for complete flux density bins.

Slob et al. (2022) computed the radio source counts of their MPS sample, and find they were scaled down by a factor of 40 compared to a general sample of radio-loud AGN. This number indicates that 2.4% of the radio sky at 144 MHz consists of MPS sources. Callingham et al. (2017) find that ∼4.5% of the radio source population is a PS source between 72 MHz and 1.4 GHz, which roughly corresponds to the combined frequency window of our MPS and GPS samples. Conversely, O’Dea (1998) indicates that ≈10% of the bright radio sky are GPS sources, which is 3 times more than the percentage we find. Snellen et al. (1998) report an offset of approximately 250 between the source counts of typical large-scale radio sources and the source counts of their GPS sample. This number would indicate that 0.4% of the sources in the radio sky at 2 GHz are GPS sources, 8 times lower than we report. Therefore, the abundances of GPS sources as reported by Snellen et al. (1998) and O’Dea (1998) vary by a factor of roughly 25 of each other.

The Snellen source counts offsets were calculated by rescaling the 325 MHz radio source counts from the WENSS mini-survey region (Rengelink et al. 1997) to 2 GHz using a spectral index of –0.85. Since Snellen et al. (1998) published, sensitive source counts at frequencies closer to 2 GHz have become available, therefore we calculate the offset of the Snellen et al. (1998) GPS source counts, rescaled with a spectral index of –0.8, to the 1400 Mhz model by Massardi et al. (2010). We find an offset of 51 ± 20, which is 5 times lower than what Snellen et al. (1998) originally reported, and corresponds to a GPS abundance of 1.9 ± 0.7% of the radio sky at 2 GHz. This number is closer to the value we report.

If we assume all PS sources evolve into FRI and FRII radio sources and there is no cosmological evolution, the source count offset captures the difference in lifetimes between PS sources and large-scale AGN. This assumption appears to be valid as Slob et al. (2022) did not find a significant cosmological evolution in the luminosity functions of peaked spectrum sources. We would then find that MPS sources have lifetimes roughly 44 times shorter than large-scale radio AGN, GPS sources have lifetimes roughly 28 times shorter than AGN, and MPS sources have lifetimes 1.6 times shorter than GPS sources.

To potentially understand the different relative lifetimes of MPS and GPS sources, we need to understand the jet formation and interaction with the gas in the host galaxy of young AGN. Sutherland & Bicknell (2007) produced three-dimensional simulations of the interaction of a light hypersonic jet with an inhomogeneous thermal and turbulently supported disk in an elliptical galaxy. The phase in which the jet breaks out of the energy-driven bubble occurs when the source has a linear size of roughly 1 kpc, corresponding approximately to when a GPS source transitions into an MPS source. The jets grow roughly twice as fast once they have escaped the optical host. This faster jet evolution suggests that for AGN, the evolution of the MPS phase is likely quicker than the evolution of the GPS phase. This difference in evolution rates could cause MPS sources to have shorter lifetimes, causing MPS sources to be less abundant than GPS sources in flux-density limited surveys.

It should be noted that the boundaries of the different classes of PS sources are not based on physical parameters, but rather on observational limits. Furthermore, we only compare the observer-frame turnover frequencies with each other, therefore our MPS and GPS samples might contain more populations than just MPS and GPS sources when considering the rest-frame turnover frequencies. However, if we assume that the redshift distribution of our MPS and GPS samples is similar, we can still conclude that PS sources with higher turnover frequencies are more abundant and have longer lifetimes than sources with lower turnover frequencies.

Note that our inference about the lifetime of PS sources from source counts ignores cosmological evolution. To account for this evolution, luminosity functions of our MPS and GPS samples are required, as it would allow us to calculate the relative abundances of MPS and GPS sources as a function of redshift. Slob et al. (2022) have computed luminosity functions for their sample of MPS sources, and found that there is no cosmological evolution between redshifts 0.1 and 0.3. From this result, we consider it plausible that the effect of cosmological evolution will be limited.

To further ensure a valid comparison between the MPS and GPS samples, we compare their redshift distribution. If there is a significant difference in the redshift distribution of our samples, the difference in the source counts would be impacted. In that case, luminosity functions would be required to determine the difference in abundance. However, we do not expect a difference between the redshift distributions between the MPS and GPS samples due to the similarity of their known optical host galaxies (Labiano et al. 2007; Slob et al. 2022). To ensure there was no significant difference between the redshift distributions of the GPS and MPS samples, we obtained the redshifts of the sources in our samples from the LoTSS DR1 and DR2 optical catalogs (DR1: Williams et al. 2019, DR2: Hardcastle et al. 2023). Redshifts were available for 237 of the MPS sources and 3685 of the GPS sources. We find no significant difference in the available redshift distribution of the two samples. Both distributions peak near z = 1. Since there is no significant difference between the redshift distribution of the MPS and GPS sample, we can conclude that the difference in source counts offset is caused by a difference in abundance.

7. Conclusions

The recent revolution in wide-field radio surveys, both in sensitivity, sky area, resolution, and frequency range, has made the study of statistical samples of PS sources possible for the first time. These samples allow us to study the abundance of GPS and MPS sources in the radio sky more accurately than ever before and allow us to draw conclusions about their relative lifetimes. In this work, we present a sample of 8032 GPS sources with spectral turnovers near 1400 MHz, and a sample of 506 MPS sources with turnovers near 144 MHz. These samples have been selected using LoTSS, LoLSS, NVSS, and VLASS. Our MPS sample is 1.4 times larger than the previously known largest sample of PS sources in the same frequency range as constructed by Slob et al. (2022). Our GPS sample is over five times larger than the previously known largest sample of PS sources (Callingham et al. 2017).

These samples of PS sources were defined from two master samples, consisting of unresolved and isolated radio sources with no assumption about their spectral shape. From our LF master sample, we classified sources as MPS sources if the power law spectral index between LoLSS and LoTSS is positive, and the spectral index between LoTSS and NVSS is negative. Similarly for the HF sample, a source is considered a PS source if the spectral index between LoTSS and NVSS is positive, and the spectral index between NVSS and VLASS is negative. This criteria selects sources with a concave SED. 3.9% of the sources in our LF sample were identified as MPS sources, and 7.4% of the sources in our HF sample were identified as GPS sources. We find the following:

• We calculated the 144 MHz respectively 1400 MHz Euclidean normalized source counts for our MPS and GPS samples and found that their shapes match those of the source counts of large-scale AGN, scaled down by a factor 44 ± 2 for the MPS sample and a factor 28 ± 1 for the GPS sample. From this offset, we find that 2.22 ± 0.09% of the radio source population at 144 MHz are MPS sources, and 3.4 ± 0.1% of the radio source population at 1400 MHz are GPS sources.

• If we interpret these offsets at face value, they imply that MPS sources have lifetimes approximately 44 times shorter than large-scale AGN, and GPS sources have lifetimes roughly 28 times shorter. They also imply that GPS sources live roughly 1.6 times longer than MPS sources.

• We were able to identify four PS sources with a lower spectral index of > 2.5, which is higher than the limit that can be associated with SSA. These sources must thus have spectral turnovers associated with FFA, indicating that these sources are likely associated with a dense circumnuclear environment.

• We identified B1315+415 as a possible restarted PS source since it has extended lobes and a compact core, and the SED of this source shows a spectral turnover.

It should be noted that we are unable to identify GPS sources with a spectral turnover at higher frequencies than 2 GHz, while the literature defines GPS sources as sources with a turnover up until approximately ∼5 Ghz. Therefore we are likely missing a large fraction of GPS sources. However, the definitions of the different categories of PS sources are based on arbitrary parameters, and our conclusions hold for the frequency windows we define.

The conclusions we draw from the source counts should be interpreted carefully, as they do not account for any redshift evolution. Using luminosity functions, we would be able to study the evolution of the abundance of PS sources as a function of redshift.

Once such luminosity functions are constructed, modeling them could help us understand the slope of the luminosity functions. We would then have a better understanding of how the number density of PS sources evolves with redshift, which allows us to understand the evolution between MPS and GPS sources.

As was discussed in Sect. 5.3, we have used an S/N-dependent selection limit in this work, to account for any potential overestimation of the amount of PS sources due to the broadening of the distribution of spectral indices. We have ensured that applying both a hard and an S/N-dependent selection limits results in a larger offset for the MPS sample than the GPS sample. Our conclusions would be further solidified by simulating the noise properties of the different surveys used in this work and doing injection tests of sources with known spectra. Such an analysis would rely on knowing the intrinsic distribution of the spectral indices, for both our LF and HF sample. It would also require we can perform the same source finding for all the surveys, which might pose a challenge, especially for the older surveys such as NVSS. Though we have shown that our conclusions about MPS sources being less abundant than GPS sources are robust, doing an analysis such as this would be useful to be sure about the true abundances of MPS and GPS sources, independent of any selection criteria.

Since we only used three flux density measurements to determine and classify the SEDs of our sources, we were limited in determining their exact turnover frequency. Spectral index fitting on more than three flux density measurements would allow us to determine the turnover frequency more precisely, which could help us determine the linear sizes of these sources according to the relation between linear size and turnover frequency by O’Dea (1998). We could then plot these sources in a diagram of, for example, linear size against radio power, for the sources with available redshift. Such a diagram provides a snapshot of radio source evolution, where individual sources will trace out trajectories on the plane. Comparing our MPS and GPS samples on this plane will be valuable to further study the evolution between these types of sources.

We have shown that the improvement in sensitivity of wide-field radio surveys in the past few decades allows us for the first time to make robust statements about the relative abundances of MPS and GPS sources. After completion of the radio surveys used in this work, we can expect an even further increase in PS sample size. The size of our GPS sample is limited by the survey area of LoTSS, which has currently covered 27% of the northern sky (Shimwell et al. 2022). Upon completion of LoTSS, we thus expect to find ∼32 ⋅ 10³ GPS sources. The MPS and LF samples are limited by the LoLSS sky area, which will also eventually cover the entire northern sky. Currently, 650 degree² has been observed (de Gasperin et al. 2023), which is only ∼3% of the northern sky. We would then expect to find ∼17 ⋅ 10³ MPS sources once both LoLSS and LoTSS have been completed.

Aside from the completion of radio surveys already being made, new radio telescopes that will be built in the coming years promise an even greater sensitivity and resolution. LOFAR is a pathfinder for the Square Kilometer Array (SKA), which is a radio telescope currently being built in Australia and South Africa. SKA’s mid-frequency array will be almost five times more sensitive than the VLA (SKAO 2022). SKA will provide opportunities for studying larger and more sensitive samples of PS sources than ever before.

Data availability

The catalogs containing the MPS and GPS sample, as described in Table A.1 are available at the CDS via anonymous ftp to cdsarc.cds.unistra.fr (130.79.128.5) or via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/689/A264

Acknowledgments

LOFAR data products were provided by the LOFAR Surveys Key Science project (LSKSP; https://lofar-surveys.org/) and were derived from observations with the International LOFAR Telescope (ILT). LOFAR van Haarlem et al. (2013) is the Low Frequency Array designed and constructed by ASTRON. It has observing, data processing, and data storage facilities in several countries, which are owned by various parties (each with their own funding sources), and which are collectively operated by the ILT foundation under a joint scientific policy. The efforts of the LSKSP have benefited from funding from the European Research Council, NOVA, NWO, CNRS-INSU, the SURF Co-operative, the UK Science and Technology Funding Council and the Jülich Supercomputing Centre. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. This research used the facilities of the Canadian Astronomy Data Centre operated by the National Research Council of Canada with the support of the Canadian Space Agency. This research has made use of "Aladin sky atlas" developed at CDS, Strasbourg Observatory, France, and of the NASA/IPAC Extragalactic Database (NED) operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. Furthermore, this research has made use of NASA’s Astrophysics Data System Bibliographic Services. This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France. This research also made use of TOPCAT, an interactive graphical viewer and editor for tabular data (Taylor 2005). This research made use of APLpy, an open-source plotting package for Python Robitaille & Bressert (2012), Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration 2018), and of matplotlib, a Python library for publication quality graphics (Hunter 2007). This research made use of NumPy (Harris et al. 2020), and of SciPy (Virtanen et al. 2020).

References

Appendix A: Column definitions

Appendix B: Source counts

All Tables

All Figures

Fig. 1. Observational limits of different radio surveys, as well as literary PS samples which are indicated at their central frequency. The lowest cataloged component flux density reported for each survey used in this research, and other relevant large radio surveys, is indicated in the plot. The GaLactic and Extragalactic All-Sky MWA Survey (GLEAM; Wayth et al. 2015) is represented by a line since it has variable limiting flux densities at different observing frequencies. The Callingham et al. (2017) sample is represented by an SSA model with a faintest peak flux density of 160 mJy at a turnover frequency of 190 MHz. The Slob et al. (2022) sample is represented by an SSA model with a faintest peak flux density 80 mJy at a turnover frequency of 100 MHz. We indicate the sensitivity of our samples with SSA models with a limiting flux density of 23 mJy at a turnover frequency of 100 MHz for the MPS sample, and with a limiting flux density of 3.9 mJy at 750 MHz for the GPS sample. In the figure, multiple wide-field radio surveys were indicated that have not been mentioned before; the Westerbork Northern Sky Survey (WENSS; Rengelink et al. 1997), Sydney University Molonglo Sky Survey (SUMSS; Mauch et al. 2003), Cambridge 7C (Hales et al. 2007), Australia Telescope 20 GHz (AT20G; Murphy et al. 2010), Very Large Array Low-Frequency Sky Survey Redux (VLSSr; Lane et al. 2014), TIFR GMRT Sky Survey Alternative Data Release 1 (TGSS; Intema et al. 2017), and the RAPID ASKAP Continuum Survey (RACS; Hale et al. 2021).

In the text

Fig. 2. Color-color plot for the 12 962 sources in our MPS sample, where αlow is the spectral index between LoLSS (54 MHz) and LoTSS (144 MHz), and αhigh is the spectral index between LoTSS and NVSS (1400 MHz). PS sources exist in the lower right quadrant and are indicated in red. The blue line represents the median of the error of αlow at 0.18. The contours represent 17, 60, 168, and 308 sources. It should be noted that there might be PS that are not indicated in red because they are hidden by the contours. The red dashed line represents the 1:1 ratio of the spectral indices. The median error bars are plotted instead of the individual errors. The normalized distributions of αlow and αhigh are plotted, with a median and standard deviation of −0.6 ± 0.4 respectively −0.7 ± 0.2. Mock SEDs are shown, to indicate the rough shape the SED of a source in each quadrant would have.

In the text

Fig. 3. Color-color plot for the 108 473 sources in our GPS sample, where αlow is the spectral index between LoTSS (144 MHz) and NVSS (1400 MHz), and αhigh is the spectral index between NVSS and VLASS (3 GHz). PS sources exist in the lower right quadrant and are indicated in red. The blue line represents the median of the error of αlow at 0.10, which corresponds to the median of the selection limit of PS sources. The contours represent 134, 416, 1349, and 2800 sources. The red dashed line represents the 1:1 ratio of the spectral indices. The median error bars are plotted instead of the individual errors, for readability. The normalized distributions of αlow and αhigh are plotted, with a median and standard deviation of −0.6 ± 0.4 respectively −0.7 ± 0.5. Mock SEDs are shown, to indicate the rough shape the SED of a source in each quadrant would have.

In the text

Fig. 4. LoTSS image and SED of B1315+415, which is likely a restarted PS source. Top: The LoTSS view of B1315+415. The LoTSS beam size is indicated in the upper left corner, and the color bar indicates the flux density per beam. The lobes extend to roughly 200 kpc. Bottom: The SED of B1315+415, which contains flux density measurements from the B3 survey (Ficarra et al. 1985), the Very Long Baseline Array (VLBA) Calibrator Survey (VCS; Kovalev et al. 2007), the NRAO Green Bank telescope (Becker et al. 1991), the VLBA Imaging and Polarimetry Survey (VIPS; Helmboldt et al. 2007), the VLA (Patnaik et al. 1992), and the Effelsberg telescope (Vitale et al. 2015).

In the text

Fig. 5. SEDs of the four sources in our sample with extremely steep spectral indices. The curve from Eq. (3) was fitted to the radio flux density measurements available for these sources and is indicated as a black curve. Upper left: ILTJ144448.65+560615.6, which has a spectral index between LoLSS and LoTSS of αlow = 2.5 ± 0.4. Upper right: ILTJ151214.28+542710.1 with a spectral index between LoLSS and LoTSS of αlow = 2.9 ± 0.5. Lower left: ILTJ132201.76+564652.5 with a spectral index between LoLSS and LoTSS of αlow = 2.5 ± 0.3. Lower Right: 87GB B0903+6656 (ILTJ090723.44+664446.6), with a spectral index between LoTSS and NVSS of αlow = 2.7 ± 0.1. For this source, we additionally plot a 4.58 GHz observation taken by the NRAO Green Bank telescope (Becker et al. 1991), and an 8.4 GHz observation taken by the VLA (Healey et al. 2007).

In the text

Fig. 6. SED of NGC 3894 (ILTJ114850.36+592456.2), which is an example of a source with a convex SED. Such an SED might be caused by multiple epochs of activity. Additional high-frequency measurements were added; a 5 Ghz NRAO 300-ft telescope measurement (Sramek 1975), a 10.45 GHz 100 m Effelsberg telescope measurement (Pasetto et al. 2016), and two VLBI measurements at 8 GHz and 415 GHz (Tremblay et al. 2016).

In the text

Fig. 7. 144 MHz Euclidean normalized differential source counts. Top: The source counts of our MPS sample. Plotted are also the source counts from Callingham et al. (2017) and Slob et al. (2022), which are samples of PS sources. To compare them to a general population of AGN, the 150 MHz source counts from the LoTSS Deep Field (Mandal et al. 2021) are plotted. The 150 MHz source counts from the T-RECS simulation (Bonaldi et al. 2019) are plotted as well. The 150 MHz model of large-scale AGN from Massardi et al. (2010) is indicated by the solid gray line. This same model scaled down by a factor of 44 is indicated by the dashed gray line. Bottom: The residuals of our MPS source counts to the Massardi et al. (2010) 150 MHz model rescaled by a factor 44.

In the text

Fig. 8. 1400 MHz Euclidean normalized differential source counts. Top: The source counts of our GPS sample. Plotted are also the source counts of the PS sample from Snellen et al. (1998). To compare our source counts to those of a general population of AGN, the 1400 MHz source counts from De Zotti et al. (2009) are plotted. The 1400 MHz source counts from the T-RECS simulation (Bonaldi et al. 2019) are indicated. The 1400 MHz model of large-scale AGN from Massardi et al. (2010) is indicated by the solid gray line. This same model scaled down by a factor of 28 is indicated by the dashed gray line. Bottom: residuals of our GPS source counts to the Massardi et al. (2010) 1400 MHz model rescaled by a factor 28.

In the text

	Fig. 9. Distributions of σα_low for our LF and HF samples. The distribution of the LF sample is much broader than the distribution of the HF sample.
In the text