© The Authors 2022

1. Introduction

Gigahertz peaked-spectrum (GPS), compact steep-spectrum (CSS), and high-frequency peaked (HFP) sources are different classes of peaked-spectrum (PS) sources. PS sources are radio-loud active galactic nuclei (AGN), defined by a small linear size and a spectral peak in their broad-band radio spectra (O’Dea & Saikia 2021). The differentiation between the GPS, HFP, and CSS classes of PS sources is largely based on the frequency of their spectral peak and the maximum linear size of the source.

GPS sources are defined to have a spectral peak in the ∼0.4 to ∼5 GHz frequency range (Gopal-Krishna et al. 1983; O’Dea & Saikia 2021), while HFP sources have a spectral peak ≳5 GHz (Dallacasa et al. 2000). Both GPS and HFP sources have very small projected linear sizes of ≲-pagination1 kpc. On the other hand, CSS sources have larger linear sizes that exceed 1 kpc and have expected peaks below 400 MHz (Fanti et al. 1990). There is also a newly suggested class of PS sources with observed peak frequencies below 1 GHz; they are sometimes referred to as megahertz peaked spectrum (MPS) sources. These MPS sources have been placed in the continuum of PS sources (Falcke et al. 2004; Coppejans et al. 2015; Callingham et al. 2017).

There are two hypothesised scenarios for the small linear scales and the spectral properties of peaked spectrum sources. The first of these is the youth model, which argues that these PS sources are the precursors to massive radio-loud AGN. The youth model is supported by the morphology of PS sources, since many display a double-lobed structure on small scales. It has been shown that a relationship between radio power and linear size of PS sources exists (Kunert-Bajraszewska et al. 2010; An & Baan 2012), as well as a relation between turnover frequency and linear size (O’Dea & Baum 1997; Snellen et al. 2000). These relations suggest that double-lobed radio sources evolve from HFP sources into GPS sources, to CSS sources, and finally into Fanaroff-Riley I and II (FRI and FRII) sources (Carvalho 1985; Kunert-Bajraszewska et al. 2010). This evolutionary model has been further supported by age estimates from spectral break modelling and observations of the motion of hot spots from high-resolution imaging (Owsianik & Conway 1998; Kaiser & Best 2007).

A problem with the youth scenario is that it has been suggested that there is an overabundance of PS and CSS sources relative to large AGN (Kapahi 1981; Peacock & Wall 1982; O’Dea 1998; An & Baan 2012). This would imply the youth model cannot solely explain the existence of all PS sources. An alternative hypothesis is that PS sources are ‘frustrated’; in other words, these sources are not young galaxies, but are instead confined to small spatial scales because of extremely dense gas in their central environment. The frustration hypothesis is supported by observations of radio morphologies of CSS sources, which imply strong interactions between individual sources’ radio jets and their environments (Wilkinson et al. 1984; van Breugel et al. 1984; Kunert-Bajraszewska et al. 2010). Further evidence for the frustration hypothesis is given by the detection of extended emission around PS sources, implying multiple epochs of activity (Baum et al. 1990; Stanghellini et al. 1990). Studies of individual sources have found evidence of unusually high densities of the surrounding medium of these sources (e.g. Peck et al. 1999; Callingham et al. 2015; Sobolewska et al. 2019), further supporting the frustration hypothesis. For specific PS sources, both the youth and the frustration scenario may apply, since young sources with constant AGN activity could break through their dense surrounding medium (An & Baan 2012). Additionally, there is evidence that a fraction of PS sources do not grow into large AGN, but rather turn off and fade because their radio activity has stopped (e.g. Kunert-Bajraszewska et al. 2006; Orienti et al. 2010).

One way to differentiate between the two imposed hypotheses is by determining the absorption mechanism that causes the spectral peak. For most PS sources a synchrotron self-absorption (SSA) model describes the observed relation between peak frequency and linear size well (e.g. de Vries et al. 2009; Snellen et al. 2000), while for individual sources surrounded by a dense medium a free-free absorption (FFA) mechanism fits the observed turnover better (e.g. Bicknell et al. 1997; Peck et al. 1999; Callingham et al. 2015). In order to differentiate between the SSA and FFA mechanisms, accurate spectral data below the turnover frequency is required. For PS sources this implies that highly sensitive data at frequencies < 200 MHz is needed (Snellen et al. 2009). Additionally, an accurate low-frequency luminosity function of PS sources would be invaluable in testing how the luminosity of PS sources evolve relative to large-scale AGN, informing us whether the youth model is compatible with the observed luminosity evolution. Low-frequency wide-field surveys provide a way to homogeneously select samples of PS sources such that we can characterise incompleteness issues that have plagued previous attempts at computing luminosity functions (e.g. Snellen et al. 2000).

Recent developments in low radio frequency telescopes have led to more reliable characterisation of low-frequency spectra and better sensitivity. Wide-field surveys from these telescopes include the GaLactic and Extragalactic All-sky Murchison Widefield Array (GLEAM; Wayth et al. 2015) survey, the TIFR GMRT Sky Survey (TGSS; Intema et al. 2017), the LOFAR Two-Metre Sky Survey (LoTSS; Shimwell et al. 2022), and the LOFAR LBA Sky Survey (LoLLS; de Gasperin et al. 2021). In particular, Callingham et al. (2017) used the GLEAM survey to identify 1483 sources with spectral peaks between 72 MHz to 1.4 GHz, which doubled the number of known PS sources. Callingham et al. (2017) and Keim et al. (2019) both found that while SSA describes the turnover of a large subset of PS sources, a fraction of PS sources has to be described by a FFA model as the spectral slope below the peak violates the theoretical limit of SSA.

Recent LOw Frequency ARray (LOFAR; van Haarlem et al. 2013) radio surveys have started a revolution in high-sensitivity data at low frequencies. More sensitive observations enable us to potentially identify high-redshift PS sources, as well as low-luminosity PS sources. Previous PS samples have median 5 GHz radio luminosities of 10²⁶ − 10²⁷ WHz⁻¹ (e.g. O’Dea 1998; Snellen et al. 1998; Callingham et al. 2017). However, Shimwell et al. (2019) predicts that with LoTSS, PS sources with radio powers < 10²⁵ W Hz⁻¹ can be identified. It has been argued that these low-power compact sources can be the short-lived young radio sources that could explain the overabundance of PS sources compared to radio-loud sources (Kunert-Bajraszewska et al. 2010).

LOFAR has two wide-field low-frequency sky surveys in progress. The first of these surveys, LoTSS, has had a second data release (DR2) that covers 27% of the northern sky at 120–168 MHz, with a median sensitivity of 83 μJy beam⁻¹, and a resolution of 6″ (Shimwell et al. 2022). The second LOFAR survey used in this study is LoLSS, which has recently had a preliminary data release (PDR) that cover 3% of the northern sky, with a total of 25 247 identified sources at observing frequencies in the range 42–66 MHz, with a resolution of 47″ and a median sensitivity of 5 mJy beam⁻¹ (de Gasperin et al. 2021). The third radio survey that is important in this research is the NRAO VLA Sky Survey (NVSS; Condon et al. 1998). NVSS covers 82% of the northern sky at an observational frequency of 1.4 GHz, a resolution of 45″, and a median sensitivity of 0.45 mJy beam⁻¹.

The purpose of this paper is to combine observations from LoTSS, LoLSS, and NVSS to identify low-luminosity PS sources with spectral peaks at frequencies ∼150 MHz. We investigate whether the characteristics of these newly identified PS sources are consistent with the populations of PS sources identified by Callingham et al. (2017), and with PS source populations with spectral peaks at gigahertz frequencies. In particular, we aim to investigate the evolution of PS sources to ascertain their role in the evolution of radio-loud AGN.

The surveys and selection criteria used to select PS sources are outlined in Sects. 2 and 3, respectively. In Sect. 4 we cross-match our sample of PS sources to known GPS compact-symmetric objects (CSO) and HFP sources. The 5 GHz radio luminosities of our PS sources are presented in Sect. 5. In Sect. 6 the classification of PS sources according to their black hole accretion mechanism is outlined. Finally, we compute and analyse the radio source counts and luminosity functions in Sects. 7 and 8, respectively. Throughout this paper we adopt the standard Lambda cold dark matter cosmological model, with parameters Ω_M = 0.27, Ω_Λ = 0.73, and the Hubble constant H₀ = 70 km s⁻¹ Mpc⁻¹ (Hinshaw et al. 2013).

2. Surveys

In order to identify PS sources, flux density measurements from at least three different radio frequencies are needed. The sensitivity and observing frequency of the surveys used to identify PS sources dictate the peak frequencies and peak flux densities of the selected PS sources. Previous PS source studies were mostly limited to surveys above ∼500 MHz, and therefore identified sources with spectral peaks in the gigahertz frequency range (e.g. Snellen et al. 1998; O’Dea 1998). However, recently Callingham et al. (2017) have shown that it is now possible to identify large samples of PS sources with spectral peaks around ∼100 MHz.

The three surveys used to identify PS sources in this study were LoTSS, LoLSS, and NVSS. LoTSS and LoLSS are both low-frequency radio surveys conducted by LOFAR, with observing frequencies of 120–168 MHz and 42–66 MHz, respectively. The 1.4 GHz NVSS survey was used as the high-frequency survey. For further analysis, the PS sources were also cross-matched to other radio surveys with observational frequencies between 42 and 1400 MHz. From the combination of surveys in this study, our data is sensitive to detecting sources that have spectral peaks between ∼54 and ∼1400 MHz.

Of the three surveys used to identify PS sources in this research, LoLSS has the lowest sensitivity. Since LoLSS will provide the low-frequency data point for identifying PS sources, the completeness limit of our PS source sample will largely be set by the sensitivity of LoLSS. In addition, LoLSS PDR does not currently cover the full observation area of LoTSS DR2 and NVSS, which means not all of the available coverage of LoTSS DR2 and NVSS is used in this study. We show in Fig. 1 the PS spectra that this study and previous studies identify, based on the limiting flux densities of radio surveys.

Fig. 1. Limiting flux densities and frequencies for major radio surveys. The limiting flux density for a survey is given by the faintest catalogued source flux density. The GLEAM survey is represented as a black line since it has variable limiting flux densities over its range of observing frequencies. The orange dashed curve represents the observational limit for the sample presented by O’Dea (1998). This is given by an SSA spectrum of a PS source with a peak at 750 MHz and an observed peak flux density of 300 mJy. The green dash-dotted curve represents the observational limit for the sample presented by Callingham et al. (2017). This is given by an SSA spectrum of a PS source with a peak at 190 MHz, with a faintest observed peak flux density of 160 mJy. The blue curve represents the theoretical observational limit of PS sources in this study. This is given by an SSA spectrum of a PS source with a peak at 100 MHz, with a faintest peak flux density of 80 mJy, based on the limiting flux densities of LoLSS-PDR, LoTSS, and NVSS. This figure also illustrates that in this study LoLSS is the survey that dictates the limiting flux density for identifying PS sources. The following plotted surveys were not previously mentioned: Cambridge 7C (Hales et al. 2007) survey, Westerbork Northern Sky Survey (WENSS; Rengelink et al. 1997), Texas Survey (TXS; Douglas et al. 1996), Molonglo Reference Catalogue (MRC; Large et al. 1991), Parkes (PKS; Wright & Otrupcek 1990) survey, Sydney University Molonglo Sky Survey (SUMSS; Mauch et al. 2003), MIT-Green Bank 5 GHz (87GB; Gregory & Condon 1991) survey, Parkes-MIT-NRAO (PMN; Wright et al. 1994) survey, Australia Telescope 20 GHz (AT20G; Murphy et al. 2010) survey, and the Rapid ASKAP Continuum Survey (RACS; Hale et al. 2021).

2.1. LOFAR surveys

Two new LOFAR surveys formed the basis of our study: LoTSS and LoLSS. LOFAR is a radio interferometer with 52 dipole-antenna stations across the Netherlands and Europe. Each LOFAR station consists of low-band and high-band antennae (LBA and HBA, respectively), which are used for 10–90 MHz and 110–250 MHz observations, respectively (van Haarlem et al. 2013).

The first LOFAR survey that was used in this study is LoTSS-DR2 (Shimwell et al. 2022). LoTSS-DR2 was formed from observations taken by LOFAR at 120–168 MHz between 2014 May 23 and 2020 February 05. This survey consists of 4 395 448 catalogued sources from two regions centred at 12h45m00s +44°30′00″ and 1h00m00s +28°00′00″ spanning 4178 and 1457 square degrees, respectively. In total this survey covers 27% of the northern sky. LoTSS has a 6″ resolution, and is 90% complete at 0.8 mJy (Shimwell et al. 2022). The full data reduction process for LoTSS is described in detail in Shimwell et al. (2022). For the sources detected in LoTSS, there is an optical catalogue available comprising two non-overlapping optical catalogues from LoTSS Data Release 1 (DR1; Williams et al. 2019; Duncan et al. 2019) and LoTSS Data Release 2 (DR2; Hardcastle et al., in prep.; Duncan 2022). This ancillary optical database provides us with optical source associations and multi-wavelength properties of the identified PS sources, including photometric and spectroscopic redshifts.

For sources in LoTSS-DR2, in-band spectra were also created. These in-band spectra consist of three flux density measurements with 16 MHz bandwidth and central frequencies of 128, 144, and 160 MHz. However, as shown by Shimwell et al. (2022), these in-band spectra are not reliable for most sources. We did therefore not use these in-band spectra for spectral modelling or identifying PS sources in this study. However, the in-band spectra are plotted in the spectral energy distributions to act as a visual guide in determining the reliability of a spectral fit.

The second LOFAR survey used in this research is the preliminary data release of LoLSS (de Gasperin et al. 2021), with observations at 42–66 MHz, formed from observations conducted between 2017 and 2019. This survey consists of 25 247 sources centred around the Hobby-Eberly Telescope Dark Energy Experiment (HETDEX) Spring Field with Right Ascensions from 11 to 16 h and Declinations from 45° to 62° (Hill et al. 2008), covering 740 square degrees. LoLSS has an angular resolution of 47″, and is 90% complete at 40 mJy. Since LoLSS is the survey with the smallest survey area, its sky footprint corresponds to the detection area of this study. In the final data release of LoLSS, a higher resolution of 15″ with better sensitivities of 1–2 mJy will be reached (de Gasperin et al., in prep.).

2.2. NVSS

NVSS is the high-frequency radio survey we used to identify PS sources in this study. NVSS is a continuum survey formed from observations conducted by the Very Large Array (VLA) at 1.4 GHz between September 1993 and October 1996 (Condon et al. 1998). The NVSS catalogue consists of 1 773 484 sources north of a declination of −40°, covering 82% of the northern sky. NVSS has an angular resolution of 45″ and is 99% complete at 3.4 mJy.

2.3. Additional radio surveys

In this study we used LoTSS, LoLSS, and NVSS to identify PS sources. However, after the initial identification of these PS sources, we also cross-matched these sources to other wide-field radio surveys. The additional radio surveys used in this research are the Very Large Array Low-Frequency Sky Survey Redux (VLSSr; Lane et al. 2014) at 74 MHz, the TIFR GMRT Sky Survey Alternative Data Release 1 (TGSS-ADR1; Intema et al. 2017) at 150 MHz, and the Faint Images of the Radio Sky at Twenty-centimeters (FIRST; Becker et al. 1995) survey at 1.4 GHz. If a PS source had counterparts in VLSSr and/or TGSS, the flux densities of these surveys were used for spectral modelling of the source. The source counterparts from the FIRST survey were only used as visual guides in the spectral energy distribution plots to confirm the accuracy of spectral fits to the other surveys.

3. PS source selection

Our PS source sample was formed by making cuts based on resolution, isolation, and whether a peak occurs in their spectra. The different selection criteria, and the number of sources left in our sample after each cut, are summarized in Table 1. The details and justifications of the five selection steps are provided below.

3.1. Source isolation, resolution, and cross-matching

1. To ensure that the derived radio spectra were not impacted by source confusion, any source that was not isolated was removed from our sample. A LoTSS source was deemed isolated if it had no other source within a 47″ radius. This isolation radius corresponds to the angular resolution of LoLSS, the lowest of all the surveys used in this study. A 47″ isolation radius ensures that sources in LoLSS are not composed of multiple independent LoTSS sources. However, some bright (≳200 mJy) sources have deconvolution errors introduced by the data reduction pipeline of LoTSS, which give rise to nearby sources that were incorrectly catalogued. To ensure we did not remove bright sources from the final sample due to deconvolution errors, the flux densities of nearby sources were also considered before flagging a source as not being isolated. If a neighbouring source had a flux density ≤-pagination10% of the central source, the source was deemed as isolated. Even if the faint source is a not a deconvolution artefact, a flux density ≲-pagination10% relative to the bright source means the impact on the final spectrum is negligible. This isolation selection criteria reduced the sample to roughly 60% of the total LoTSS catalogue, to 2 676 735 sources.

2. Since PS sources have small spatial scales (≲1″; e.g. O’Dea 1998; Chhetri et al. 2018), all sources that are resolved in LoTSS were removed from our sample as it implies the source has structure > 6″. Shimwell et al. (2022) have defined a criterion for identifying resolved sources in LoTSS. This criterion identifies a source as resolved when the natural logarithm of the ratio of the integrated flux density (S_I) to peak flux density (S_P), given by R = lnS_I/S_P, is greater than or equal to R_99.9, given by

   $$\begin{array}{c}\hfill R_{99.9}=0.42+\left(\frac{1.08}{1+{\left(\frac{\mathrm{SNR}}{96.57}\right)}^{2.49}}\right),\end{array}$$(1)

   where SNR is the signal-to-noise ratio, defined as $\frac{S_{\mathrm{I}}}{{\sigma }_{\mathrm{I}}}$, with σ_I the statistical error on the integrated flux density. All sources with R ≥ R_99.9 were classified as resolved, and were removed from our sample. This removed ≈6.8% sources from the previous cut, leaving us with 2 493 574 sources.

3. The above unresolved cut has removed most resolved sources. However, since the resolved criterion is based on the 99.9 percentile of a distribution (R_99.9), and the isolated criterion has a flux density cut, we need to include another selection step to ensure all unresolved, non-isolated sources are removed. This was done by removing all sources classified in the LoTSS catalogue as C type by PyBDSF (Mohan & Rafferty 2015). These are sources fit by a single Gaussian, but within an island of emission that contains other sources, such as radio relics around clusters. The remaining sample only contains sources classified by PyBDSF as S- and M-type sources. The S-type sources are isolated sources fitted with a single Gaussian, while the M-type sources are fitted with multiple Gaussians. These M-type sources were not removed from our sample since the extended emission of these sources is mostly caused by deconvolution errors around bright sources. The multiple Gaussians fitted to these sources thus generally do not represent actual extended emission. A total of nine C-type sources were removed from our sample.

4. The remaining sources from LoTSS were then cross-matched to LoLSS, NVSS, VLSSr, TGSS, and FIRST. This was done using the Tool for OPerations on Catalogues And Tables (TOPCAT; Taylor 2005) Starlink Tables Infrastructure Library Tool Set (STILTS) multi-table cross-matching tool. A cross-matching radius of 15″ was used for all surveys. As described in Sect. 2, cross-matches in LoLSS, LoTSS, and NVSS are required to identify whether a source is a PS source. We therefore remove all sources that do not have a cross-match in LoLSS or NVSS from our sample. Since LoLSS has a smaller sky coverage than LoTSS, and both LoLSS and NVSS have lower sensitivity than LoTSS, only ≈0.4% of the isolated unresolved LoTSS sources have counterparts in both surveys. After this cross-matching step we are left with 9768 sources in our sample, from which we can identify PS sources. We refer to this sample of sources as the Master sample.

3.2. Spectral classification

From the Master sample, PS sources were identified by their spectra. A defining feature of PS sources are the power-law components above and below the spectral peak. Firstly, we fit a generic non-thermal power-law model of the form

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

where S_ν is flux density, a is the amplitude of the synchrotron emission in Jy, ν is the frequency in MHz, and α is the spectral index. This model was fit to the flux densities from LoLSS and LoTSS to obtain the low-frequency spectral index α_low. The power-law model was also fit to the flux densities of LoTSS and NVSS to obtain the high-frequency spectral index α_high. Since our fitting routine uses too few data points to calculate the uncertainties on α_low and α_high, we estimate these uncertainties by fitting a power law to the 1σ upper and lower uncertainties of the flux densities from one survey to the lower and upper uncertainties of the second survey, respectively. The corresponding limiting values of α_low and α_high for these fits then represent the 1σ limits of each respective α_low and α_high.

All sources can now be situated in radio colour-colour phase space, as given by the α_low and α_high of the sources. In this radio colour-colour phase space, spectral peaks can be identified (e.g. Sadler et al. 2006; Callingham et al. 2017). The radio colour-colour phase space for the 9768 sources in our Master sample, as obtained from the flux density points of LoLSS, LoTSS, and NVSS, is presented in Fig. 2. As expected, most sources cluster around a median of (α_low, α_high) =(−0.6 ± 0.2, −0.8 ± 0.1) in the third quadrant of Fig. 2. This is consistent with previous spectral index studies at similar observing frequencies, although our sample has a higher median value and a larger semi-interquartile range for α_low than previous studies (e.g. Tasse et al. 2007; Lane et al. 2014; Callingham et al. 2017). This larger standard deviation in α_low is likely due to larger uncertainties in calibrating the flux density scales for LoLSS compared to the higher frequency surveys that were used in previous studies (de Gasperin et al. 2021).

Fig. 2. Radio colour-colour diagram for the 9768 LoTSS sources left in our Master sample. αlow represents the spectral index between the LoLSS (54 MHz) and LoTSS (144 MHz) flux density points. αhigh represents the spectral index between the LoTSS and NVSS flux density points. To better illustrate the large number of sources around the median of αlow and αhigh, given by (αlow, αhigh) =(−0.6 ± 0.2, −0.8 ± 0.1), contours and a density map are plotted in this region. The contour levels represent 15, 75, 145, and 265 sources. The colours in the underlying density map show the number of sources for each pixel. The number of sources corresponding to each corresponding shading are illustrated by the colour bar at the top left of the plot. The grey lines show spectral indices of zero to illustrate the four quadrants of the plot. The dashed red line represents a one-to-one relation between αlow and αhigh. The blue line at αlow = 0.1 illustrates the selection limit used to identify PS sources in this study. In the corner of each quadrant the shape of a typical spectrum in that quadrant is shown in grey. To avoid confusion, individual error bars are not plotted, but in the top left of the plot the median error bar size is shown. The histograms at the top and right of the diagram illustrate the distributions of αlow and αhigh, respectively. These distributions have been normalised to the maximum value in the distribution. In these histograms Gaussian fits to these distributions are overplotted in red, and the dashed black lines show the median values for these distributions.

Sources in the third quadrant of Fig. 2 have spectra that are described by an optically thin synchrotron power law. Sources in the first quadrant of Fig. 2 follow a positive power law from 54 MHz to 1.4 GHz. These sources are expected to have a spectral peak in the gigahertz range consistent with archetypal GPS sources. Sources in the fourth quadrant have convex spectra. These sources could have another spectral turnover at a frequency ≳1 GHz, which could indicate multiple epochs of AGN activity. Sources with a spectral turnover between 54 MHz and 1400 MHz are located in the second quadrant of Fig. 2. The PS sources we are interested in are therefore selected from this section of the radio colour-colour phase space.

1. Sources in the second quadrant have a peak in their spectrum around 144 MHz. Not all sources in this second quadrant have been classified as PS sources, but instead the cut for a source to be a PS source was made at α_low ≥0.1. This cut minimises the contamination of flat spectrum sources in the selected PS sample. In the literature, most previous studies have also made a cut for PS sources at α_high ≤ − 0.5 (O’Dea 1998). However, from the continuous distribution of α_high around α_high = − 0.5 in Fig. 2, this limit appears to be arbitrary, as also concluded by Callingham et al. (2017). In order to compare the results from this study to previous studies, we make a distinction between a hard PS sample, containing PS sources with α_high ≤ − 0.5, and a soft PS sample, containing PS sources with 0.0≥ α_high ≥ − 0.5. The hard PS sample contains a total of 212 sources, and the soft PS sample contains 161 sources. The full PS sample is obtained from the combination of the soft and hard PS samples, and contains 373 sources. The spectra of a source from the soft and hard samples are shown in Fig. 3. The table with the optical and radio characteristics for the PS sample is available at the CDS in the style of the table presented in Appendix A.

   Fig. 3. Example spectra of a soft sample PS source (top) and a hard sample PS source (bottom). The flux densities of LoLSS (blue square), VLSSr (red circle), LoTSS-inband spectra (purple diamonds), LoTSS (green left-pointing triangle), TGSS (brown right-pointing triangle), NVSS (indigo downward-pointing triangle), and FIRST (pink upward-pointing triangle) are plotted where available. The black curve represents the generic curved model from Eq. (3), which was fitted to the LoLSS, LoTSS, VLSSr, TGSS, and NVSS data points. The power-law fits used to determine αlow and αhigh are shown in orange.

After applying these selection steps, we identified a sample of 373 candidate PS sources. To verify the peak in the spectra of these sources, and to better sample their spectra, we cross-matched these sources to VLSSr and TGSS. With these additional spectral points, a curved model was fit to the spectra of a subsample of our PS sources using the least-squares method. The generic curved model used for this is of the form

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

where S_ν is the flux density at frequency ν, in MHz. S_p is the flux density at the peak frequency ν_p. α_thin and α_thick are the spectral indices in the optically thin and optically thick parts of the spectrum, respectively (Snellen et al. 1998).

Since this model depends on four parameters, a cross-match to at least two more surveys in addition to LoLSS, LoTSS, and NVSS was needed to obtain a fit for this model. In total, 36 out of 373 PS sources had enough cross-matches to accurately fit a curved model. We note that we do not use the results of the curved spectral model for PS sources for further statistical analysis in this paper, but the curved spectral model can provide an extra confirmation for the identification of the PS nature for an individual source.

3.3. Flux density completeness of PS sample

For further analysis of the population, the flux density completeness of the PS sample needs to be known. Since our radio detection limit is dominated by the LoLSS flux density limit, we use this limit to compute the flux density limit at 144 MHz for our sample. We extrapolate the LoLSS 90% completeness limit at 54 MHz of S_{54 MHz, 90%} = 40 mJy (de Gasperin et al. 2021) to 144 MHz using a simple power law to find the flux limit at this selection frequency. For the PS sample, we use a spectral index for extrapolation of α_low = 0.1, corresponding to the selection limit of PS sources. This results in an estimated limiting flux density for our PS sample of S_{144 MHz, PSlim} = 44 mJy. This means a LoTSS source can only be identified as peaked-spectrum in our analysis if it has a LoTSS flux density exceeding 44 mJy.

3.4. Redshift information

For the Master sample, the redshift information from the LoTSS optical catalogues was obtained. There were two separate non-overlapping optical catalogues available for Data Release 1 (DR1; 2019 Williams et al. 2019) and DR2 (Hardcastle et al., in prep.) of LoTSS. The photometric redshifts in these optical catalogues were obtained using the methods outlined by Duncan et al. (2019) for DR1 and Duncan (2022) for DR2. For the sources in DR1 the spectroscopic redshift and the median photometric redshift were used. For the sources in DR2 the spectroscopic redshift from the Sloan Digital Sky Survey (SDSS) and the estimated photometric redshift were used if they were flagged as good quality. For both DR1 and DR2, the spectroscopic redshift was preferred when available. Optical counterparts were identified for 5026 sources in our Master sample. This resulted in available redshift data for 3303 out of 9768 sources in the Master sample. For the sources classified as PS sources in selection step 5, redshift information is available for 138 out of 373 sources. Spectroscopic redshift information was available for 54 sources in our PS sample. The corresponding redshift distribution for the PS sources is shown in Fig. 4. This distribution has a median redshift of 0.80, and a highest redshift of 5.01.

Fig. 4. Redshift distribution of the 138 PS sources with available redshift data. The red histogram is the distribution of the spectroscopic redshifts available for 54 sources. The black histogram is the distribution of the best redshifts, from the combination of spectroscopic and photometric redshifts.

4. Cross-matching to known PS samples

We tested the reliability of our PS source selection criteria by checking if our sample contains previously identified GPS, CSO, and HFP sources. This sample of previously known GPS, CSO, and HFP sources was obtained by collating the known GPS, CSS, and HFP source samples described by Callingham et al. (2017), and removing the CSS sources from this sample. Our final known GPS, CSO, and HFP sample consists of the samples isolated by O’Dea (1998), Snellen et al. (1998, 2002), Peck & Taylor (2000), Tinti et al. (2005), Labiano et al. (2007), Edwards & Tingay (2004), and Randall et al. (2011). In this sample of known PS sources, seven sources are located in the survey area of LoLSS used in this study. Of these seven sources, two have counterparts in our Master sample. The remaining five sources are too faint to be detected in LoLSS. Neither of the two known PS sources was identified via our PS sample criteria.

To understand why these two sources were not selected by our PS sample selection, we inspected the spectral properties by characterising their radio spectra between 54 MHz and 20 GHz. The spectra of these two sources are shown in Fig. 5.

Fig. 5. Spectra of the two known GPS, CSO, and HFP sources in the Master sample. The symbols represent data from the same surveys as in Fig. 3. The grey crosses show the data presented by Marvil et al. (2015), the dark green circles the data from the TXS survey (Douglas et al. 1996), and the teal hexagons the data presented by Tremblay et al. (2016). The black curve in the left plot shows the generic curved model from Eq. (3).

The source ILT J111036.37+481752.4 has a counterpart in the CSO sample presented by Peck & Taylor (2000), and has a peaked spectrum near ≈250 MHz. However, using our spectral fitting routine described in Sect. 3, we computed the spectral indices of this source to be α_low = 0.48 ± 0.14 and α_high = 0.07 ± 0.05, thus not classifying it as a PS source. This is because the spectrum appears flat between 144 MHz and 1.4 GHz. The source has an abnormally large spectral width of ≈-pagination1.8 GHz, which is substantially larger than the median FWHM of 750 MHz of the PS sample isolated by Callingham et al. (2017).

The source ILT J114850.36+592456.2 also has a counterpart in the CSO sample presented by Peck & Taylor (2000). In Fig. 5 we see that this source potentially has a convex spectrum in the range 54–1400 MHz, with a possible spectral turnover at ∼4 GHz. However, since spectral variability can be significant at high frequencies, more high-frequency data points after the turnover are needed in order to confirm this spectral turnover. A convex source with a spectral turnover above 1 GHz is suggestive of multiple epochs of AGN activity, where the low-frequency section of the spectrum is dominated by emission from aged electrons, while the higher frequency peak comes from recent core activity (Callingham et al. 2017, and references therein).

In summary, for the two known PS sources that were in our Master sample, one had a convex spectrum at our selection frequencies and the other had a relatively flat peak between our LoTSS and NVSS data points. The latter demonstrates that our PS sample will thus be a slight underestimation of the total number of PS sources in the Master sample, especially to those with wide spectral peaks. We can also conclude that the 373 sources in our PS sample are all newly identified PS candidates, increasing the number of known PS sources in our detection area by a factor of 50.

5. 5 GHz radio powers

To investigate how the radio luminosity distribution of our PS sample compares to literature PS samples, we computed the 5 GHz luminosity of the 3303 sources in the Master sample that have redshift information available. Of these 3303 sources, 138 are PS sources. The frequency of 5 GHz was chosen in order to compare the luminosities of our sample with the samples of O’Dea (1998), Snellen et al. (1998), and Callingham et al. (2017), who all evaluated their radio luminosities at 5 GHz. The 5 GHz radio luminosity, P_5 GHz, was computed using

$$\begin{array}{c}\hfill {\mathrm{P}}_{5\mathrm{GHz}}=\frac{4\pi D_L^2{S}_{5\mathrm{GHz}}}{{(1+z)}^{1+{\alpha }_{\mathrm{high}}}},\end{array}$$(4)

where D_L is the luminosity distance, S_5 GHz is the flux density of a source at 5 GHz, and the factor 1/(1 + z)^{1 + α_high} is the k-correction. We computed S_5 GHz by extrapolating the power law from Eq. (2) up to 5 GHz with spectral index α_high. This assumes the spectrum of a source follows the same power law fitted between 140 MHz and 1.4 GHz, up to frequencies of 5 GHz, without strong deviations. We expect this assumption to be valid, since deviations due to spectral curvature are not significant for our sources at frequencies above 1 GHz but below 10 GHz (Chhetri et al. 2012; Callingham et al. 2017).

The distribution of the 5 GHz radio power for the PS sample and the Master sample is provided in Fig. 6. We find that our PS sample has a median P_5 GHz value that is 10^0.7 W Hz⁻¹ higher than the median P_5 GHz value of the Master sample. The same is found for the 90% complete subsamples of the Master sample and the PS sample, which are not limited by incompleteness of LoLSS. This difference in luminosities between the two samples is likely due to the fact that for a given redshift and flux density in LoLSS, a PS source will have a higher flux density in LoTSS than a simple power-law spectrum source, due to its spectral shape. We are thus selecting relatively high flux-density sources in the PS sample, corresponding to the higher average 5 GHz radio power for PS sources.

Fig. 6. Distribution of the 5 GHz radio power of the Master sample and PS sample. The upper histogram shows the distribution of the 5 GHz radio power for the 3303 sources from the Master sample that have redshift information available in red. The black histogram in the bottom plot shows the distribution of the 5 GHz radio power of the 138 sources from the PS sample for which redshift data is available. The median redshifts of these distributions are plotted as dashed lines in their respective colours. For the total sample, the median P5 GHz value and the range of the 16th and 84th percentiles of the distribution is given by $25.0_{-0.7}^{+1.0}$ log10 W Hz−1, and for the PS source sample this is given by $25.7_{-1.1}^{+0.9}$ log10 W Hz−1.

In Fig. 7 we compare the 5 GHz radio power of our PS sample to the literature PS samples presented by Snellen et al. (1998), O’Dea (1998), and Callingham et al. (2017). The variations in 5 GHz radio power with redshift for these samples are shown in Fig. 8. We note that the 5 GHz radio powers we find for the sample presented by Callingham et al. (2017) are higher than those in the original paper, due to a missing factor of four in their radio luminosity calculation. Our PS sample contains the lowest luminosity PS source identified to date with P_5 GHz = 2.0 × 10²² W Hz⁻¹.

Fig. 7. Distribution of the 5 GHz radio power for different samples of PS samples. The distribution of the 5 GHz radio power for the 138 sources from the PS sample for which redshift data is available is shown in black (upper histogram). The green, blue, and orange histograms in the lower plot represent the 5 GHz radio power of the Snellen et al. (1998) PS sample, the O’Dea (1998) PS sample, and the Callingham et al. (2017) PS sample, respectively. The median redshifts of all distributions are plotted as dashed lines in their respective colours, in both the upper and lower plots. The median P5 GHz value and the range of the 16th and 84th percentiles of the distributions of the samples in this study, Snellen et al. (1998), O’Dea (1998), and Callingham et al. (2017), are given by $25.7_{-1.1}^{+0.9}$, $26.2_{-0.4}^{+0.5}$, $27.6_{-1.1}^{+0.7}$, and $27.2_{-1.2}^{+0.6}$ (log10 W Hz−1), respectively.

Fig. 8. 5 GHz radio power against redshift for the 138 sources from the PS sample for which redshift data is available (black circles). The green triangles, blue squares, and orange diamonds represent the PS samples of O’Dea (1998), Snellen et al. (1998), and Callingham et al. (2017), respectively. The dashed grey line corresponds to the 5 GHz luminosity limit for a source that has a peak flux of 44 mJy at 144 MHz, corresponding to the 90% completeness limit of the PS source selection. To compute this 5 GHz luminosity limit, a median spectral index of −0.87, corresponding to the 90% limit of αhigh, was used.

In Fig. 8 we also plot the curve for the estimated 90% limit of luminosity against redshift. To compute this curve, we used the estimated 90% flux density completeness at 144 MHz of 44 mJy described in Sect. 3.3 and extrapolate this to 5 GHz using a simple power law. The spectral index used for this extrapolation (α_high = −0.87) corresponds to the 90% limit of the α_high distribution for PS sources. We then use this 5 GHz flux density limit of 2 mJy to compute the corresponding estimated 90% limit of the luminosity as a function of redshift.

Compared to the PS samples from Snellen et al. (1998), O’Dea (1998), and Callingham et al. (2017), the PS sources presented in this study have 5 GHz radio powers that are, on average, roughly an order of magnitude smaller. Compared to the O’Dea (1998) and Callingham et al. (2017) samples, our PS sample does not identify any high-luminosity (P_5 GHz >  1 × 10²⁸ W Hz⁻¹) PS sources. However, as can be seen in Fig. 6, our Master sample also does not contain these high-luminosity sources. Therefore, the lack of high-luminosity PS sources is likely due to cosmic variance since the LOFAR surveys we use to select PS sources have not yet surveyed the whole sky.

6. High or low excitation classification of the PS sample

To investigate the dominant accretion mode for PS sources, we classify sources in our PS sample as high-excitation or low-excitation radio sources (HERGs or LERGs). HERGs are sources that have a radiatively efficient accretion mode, and radiate strongly across their electromagnetic spectrum (e.g. Best & Heckman 2012, and references therein). LERGs, on the other hand, have an accretion mode that leads to less strong radiative emission throughout the electromagnetic spectrum (e.g. Hardcastle et al. 2007). HERGs are shown to be a strongly evolving population, while LERGs show little cosmic evolution (Best & Heckman 2012; Pracy et al. 2016). Therefore, the classification of PS sources into HERGs and LERGs provides insight into the evolution of the PS source population.

We classified our PS sources into HERGs and LERGs using available optical spectra of our PS sources. We obtained the optical spectra for 54 of our PS sources from the 17th data release of the Sloan Digital Sky Survey (SDSS DR17; Abdurro’uf 2022).

Previous studies classified sources as HERGs and LERGs based on the relative strength and equivalent width of the 5007 Å [O III] line, with LERGs having significantly less [O III] emission than HERGs (e.g. Laing et al. 1994; Tadhunter et al. 1998). In more recent studies, sources have been classified as HERGs and LERGs based on multiple optical emission lines from separation between Seyfert and LINER galaxies proposed by Kewley et al. (2006; e.g. Baldi & Capetti 2010; Buttiglione et al. 2010). However, Buttiglione et al. (2010) noted that for a small sample of sources, the use of the Seyfert and LINER diagnostic diagrams can give ambiguous results. We therefore adapt the excitation index (EI) classification scheme defined by Buttiglione et al. (2010).

The EI combines the emission-line ratios of four emission lines, and is defined as EI = log₁₀([O III] / Hβ) – $\frac{1}{3}$[log₁₀([N II] / Hα) + log₁₀([S II] / Hα) + log₁₀([O I] / Hα)]. The division between HERGs and LERGs is made at a value of EI = 0.95 (Buttiglione et al. 2010), where HERGs have an EI above this limit and LERGs have EIs below this limit. However, the EI classification method requires the presence of all four emission lines. Unfortunately, these four emission lines were not all available for many of the SDSS spectra of the PS sources. Therefore, if a classification based on the EI was not possible, we used the [O III] emission line only. Here we identified a source as a HERG when the equivalent width of its [O III] emission line was above 5 Å, in accordance with selection step (iii) outlined by Best & Heckman (2012).

The number of HERGs and LERGs classified using the described classification methods are presented in Table 2. In total, we can classify roughly half of the PS sources with spectral information available as HERGs or LERGs using emission line criteria. For the unclassified sources, not enough emission lines were detected in order to class them. This could be due to a genuine lack of emission lines, suggesting that these sources are LERGs. However, the lack of detected emission lines can be due to low signal-to-noise spectra.

We therefore investigate the classification mechanism based on the [O III] line luminosity versus radio luminosity proposed by Best & Heckman (2012) for the unclassified sources. Figure 9 shows the distribution of PS sources in the [O III] line luminosity versus radio luminosity plane. From Fig. 9 we find that three unclassified PS sources can be placed in the [O III] line luminosity versus radio luminosity plane. However, all three unclassified sources lie above the line defined by Best & Heckman (2012) as the lower limit to the distribution of HERGs. Sources below this line can be classified as LERGs, while sources above this line can be either HERGs or LERGs. This classification mechanism can therefore not definitively classify these three sources as HERGs or LERGs. We note that none of the LERGs in our PS sample lie below the HERG–LERG division line, indicating that the detected PS LERGs have relatively strong [O III] line emission.

Fig. 9. [O III] line luminosity vs 1.4 GHz radio luminosity of sources classified as HERGs and LERGs, and unclassified PS sources. The red diamonds, orange squares, and blue triangles indicate PS HERGs, PS LERGs and unclassified PS sources, respectively. The [O III] line luminosities vs 1.4 GHz radio luminosities of the sample of HERG and LERG sources presented by Best & Heckman (2012) are plotted as red and black points, respectively. In order to avoid overcrowding of the figure, only 20% of the LERGs in the Best & Heckman (2012) sample are plotted. The solid black line indicates the lower limit to the distribution of HERGs proposed by Best & Heckman (2012).

In total, 20 out of the 25 classified PS sources are HERGs, suggesting that our PS sources are more likely to be HERGs than LERGs. This corresponds well to the hypothesised youth model for PS sources, as HERGs are shown to be strongly evolving sources at z ≲ 1 (Best & Heckman 2012; Pracy et al. 2016). Additionally, roughly one-third of the PS sources that we were able to classify as HERGs or LERGs have L_1.4 GHz ≳ 10²⁶ W Hz⁻¹, all of which have been classified as HERGs. This supports previous observations that HERGs dominate the population at high luminosities (e.g. Best & Heckman 2012; Butler et al. 2018).

In our sample of PS sources we have an overabundance of HERG sources compared to a general sample of AGN, which is dominated by LERGs (Best & Heckman 2012; Pracy et al. 2016). However, in our characterisation of PS sources we did not take redshift into account. At high redshifts we expect to find more HERGs since these will have more active spectra and stronger optical emission. This will cause them to have a higher likelihood of being detected by optical surveys.

To compare the classification of our PS sample without this redshift bias, we apply the redshift restriction z ≤ 0.3 from the HERG and LERG classification of Best & Heckman (2012). In this redshift range, we classify five PS sources as EI-LERGs, and six PS sources as HERGs (five EI-HERGs, one [O III] HERG); one source cannot be classified. Furthermore, the PS sources that were identified as HERGs in this local sample all have a 1.4 GHz luminosity < 10²⁵ W Hz⁻¹. This shows that without a redshift and luminosity bias we still find a significant overabundance of HERGs in our PS sample, as Best & Heckman (2012) and Pracy et al. (2016) identify HERGs and LERGs in a roughly one-to-ten ratio in this redshift and luminosity range, instead of the one-to-one ratio for our PS sample. Therefore, the fact that most of our PS sources are classed as HERGs indicates the population is a quickly evolving one, as is consistent with the youth model of the radio emission. However, we note that more PS sources with optical spectral information are needed to confirm if this overabundance is observed for our entire PS sample.

7. Euclidean normalised source counts

In order to investigate the evolution of PS sources, we determined the radio source counts for different samples of PS sources and compared these with the counts predicted by evolutionary models of radio populations.

The differential source counts, dN/dS, were computed for our sample of PS sources with 144 MHz flux densities ≥44 mJy, corresponding to the flux density completeness limit of our sample. The value of dN/dS was computed by summing the observed number of PS sources in six 144 MHz flux bins and dividing these by the detected area of the sky, A. In the case of our PS sample, we used the detection area of LoLSS (A = 740 deg²; de Gasperin et al. 2021) as this is the survey that limits the detection area for our PS sources. We then normalised the source counts with a factor S^2.5, which corresponds to normalising to a uniformly distributed Euclidean space. The resulting normalised differential source counts are shown in Fig. 10 and are listed in Table 3.

Fig. 10. Euclidean normalised differential source counts for different samples of PS sources. The black squares show our sample of PS sources. The blue circles show the Snellen et al. (1998) PS sample converted to 144 MHz. The orange diamonds show the Callingham et al. (2017) PS sample. For reference, the source counts for the LoTSS Deep Field derived by Mandal et al. (2021) are plotted as red triangles. The solid grey line shows the model for 150 MHz AGN source counts proposed by Massardi et al. (2010). The dashed grey line shows the Massardi et al. (2010) model scaled down by a factor of 40, which is consistent with the PS source data.

We can use this same method to construct the normalised radio source counts for the PS sample presented by Callingham et al. (2017). To construct the source counts at 144 MHz for this sample, we used the reported GLEAM flux densities at 143 MHz, and assumed the flux density does not change significantly over 1 MHz. Only sources with 143 MHz flux densities ≥1 Jy, corresponding to the estimated 100% completeness limit of this PS sample, were considered when evaluating these source counts. The computed source counts for the Callingham et al. (2017) PS sample are shown in Fig. 10 and are listed in Table 3.

We also compared the radio source counts for our sample of PS sources to the radio source counts derived by Snellen et al. (1998) for their sample of PS sources. These source counts were evaluated at the peak frequency of the individual PS sources, which was assumed to correspond to a median frequency of 2 GHz. We therefore determined the radio source counts for this sample at 144 MHz by shifting the counts from 2 GHz to 144 MHz assuming the power-law model from Eq. (2) with a spectral index of −0.80. The resulting source counts are shown in Fig. 10. We note that we excluded the highest flux density bin since this was not derived directly from the sample presented by Snellen et al. (1998) and likely contains uncharacterisable systematics.

In order to compare the radio source counts for the different PS samples to a general sample of AGN, the 150 MHz LoTSS Deep Field radio source counts determined by Mandal et al. (2021) are also plotted in Fig. 10, along with the 150 MHz model for AGN source counts presented by Massardi et al. (2010). Moving this model down by a factor of 40 roughly agrees with the observed curve for the different observed PS source counts.

The interpretation of this factor of 40 between the total and PS source counts needs to be treated with caution. While it is tempting to assume that if PS sources evolve into the larger scale sources sampled by Mandal et al. (2021), they should undergo similar cosmological evolution since their lifetime is significantly less than Hubble time. Therefore, this factor of 40 would encode the ratio of the lifetime of the two source classes and the luminosity function of PS sources. If taken at face value, this factor implies that the lifetime of the PS phase is ≈40 times shorter than the lifetime of a large radio galaxy at low frequencies.

Snellen et al. (1998) found that the source counts at 2 GHz for their PS sample are scaled down by a factor of 250 compared to the source counts of large-scale radio galaxies. If we assume this factor also encodes the lifetime of this sample of PS sources, that would mean that PS sources selected at 2 GHz have much shorter lifetimes, and thus stronger evolution than those selected at 144 MHz. This stronger evolution could be related to the jet power of the sources since at 2 GHz the detected radiation is much closer to the core than at 144 MHz, causing the detected jet power to be stronger too. Closer to the core, evolution of the galaxy would then take place more quickly than in the outer regions. However, this interpretation of the scaling factor does not take into account any redshift evolution, and is probably too simplistic as the redshift evolution of AGN and PS sources are not expected to be identical (Labiano et al. 2007; Kunert-Bajraszewska et al. 2010). To decouple the impact of any potential luminosity evolution from the source counts, we investigate the luminosity function in the following section.

8. Peaked-spectrum radio luminosity function

The second method we used to characterise the evolution of our PS sources is the luminosity function. We compute the luminosity function for both our Master sample and our PS sample at 144 MHz using the standard 1/V_max method (Schmidt 1968), with the luminosity function in a given luminosity bin centred at L given by

$$\begin{array}{cc}\hfill \frac{\mathrm{d}N(L_j)}{\mathrm{d}logL}& =\frac{1}{\mathrm{\Delta }log{L}_j}{\displaystyle \sum _{i=1}^N}\frac{1}{V_i},\hfill \end{array}$$(5)

where the sum is over all N sources in the luminosity bin, and where V_i corresponds to the volume over which a galaxy can be detected given the optical and radio selection criteria for the sample. This is calculated as V_i = V_max − V_min, where V_max and V_min are the volumes corresponding to the upper and lower redshift limits, respectively, for which a source could be detected in our sample. Below, we outline our methods for estimating our selection criteria that determine the optical and radio completeness, and V_max.

8.1. Estimating selection effects

In order to calculate the maximum volume over which a galaxy can be detected, both the radio and optical detection limits need to be characterised. For the radio detection limit of PS sources we can use the flux density limit of S_{144 MHz, PSlim} = 44 mJy described in Sect. 3.3. However, since the spectral shape of a source determines the extrapolated flux density at 144 MHz, we have to use a different spectral index for the extrapolation for the Master sample and the PS sample. For the flux extrapolation of the Master sample we used a spectral index α_low = −1.1, corresponding to the lower 95% limit to the distribution of α_low of the Master sample. This resulted in a limiting flux density for the Master sample of S_{144 MHz, lim} = 13 mJy.

In addition to the radio survey selection limits, we also need to take the optical selection limits into account. However, the LoTSS optical catalogues are a combination of different optical surveys, and thus have no defined detection limit. In order to define a consistent optical detection limit, we only used sources from our Master sample that have counterparts in SDSS-DR7 (Abazajian et al. 2009). These counterparts were obtained by cross-matching the position of SDSS sources with the positions of LoTSS sources using a cross-matching radius of 2″s. To compute the optical V_max, we used the SDSS g- and i-band photometry, implying that we only consider sources for which both the SDSS g- and i-band magnitudes are measured. In total, this Master subsample with SDSS counterparts consists of 1909 Master sample sources, of which 105 are PS sources. We then set the optical detection limit as the 95% completeness-limit of SDSS in the i-band (m_i, lim = 21.3 mag). We refer to this sample as the SDSS-selected sample.

To test the influence of the optical incompleteness on the resulting luminosity function, we also computed the luminosity functions using the Pan-STARRS g- and i-band photometry included in the LoTSS DR1 optical catalogue (Williams et al. 2019). To ensure uniform selection effects for this sample, we use only the photometric redshifts. We refer to this sample as the photo-z selected sample. In total, this results in a Master subsample of 2156 sources, of which 88 sources are PS sources. We set the detection limit of this sample as a conservative estimation for where the photometric redshifts are still well calibrated (m_i, lim = 22.5 mag).

8.2. Estimating V_max

Using the detection limits described in Sect. 8.1, we were able to calculate the maximum volume in which a source can be detected by our optical and radio surveys, V_max, opt and V_max, radio, respectively. In order to calculate V_max, opt, we first calculated the absolute i-band magnitude, M_i, for each source as

$$\begin{array}{c}\hfill M_i=m_i-\mathrm{DM}-K_i(z),\end{array}$$(6)

where m_i is the apparent magnitude, DM is the distance modulus, and K_i(z) is the k-correction. The k-corrections are calculated with the K-corrections calculator (Chilingarian et al. 2010; Chilingarian & Zolotukhin 2012) for the g-i colour. To determine V_max, opt, we evaluated Eq. (6) at m_i = m_i, lim, at a series of redshifts (Δz = 0.0001) to find the greatest redshift, z_{max, optical}, at which M_i, lim ≤ M_i, source. Above this redshift, the source is too faint to be detected by the optical survey.

Analogously, the maximum radio redshift, z_max, radio, is found by first computing the 144 MHz radio luminosity, L_144 MHz, of each source using Eq. (4), evaluated at 144 MHz instead of 5 GHz. We then compute Eq. (4) at S_144 MHz = S_{144 MHz, (PS)lim} at the same series of redshifts as the optical method. From this we can identify z_max, radio as the smallest redshift where L_{144 MHz, (PS)lim} ≥ L_{144 MHz, source}. Above this redshift, the source is too faint to be detected in our sample.

The combined maximum optical and radio redshift limits for each source z_max is now the minimum of z_max, radio and z_max, opt. We then calculated V_max from the integrated comoving volume corresponding to z_max and multiplying this volume by the fraction of the sky covered by our sample. This area of detection corresponds to the fractional sky coverage of LoLSS (740 deg²), which is the survey that limits our detection area.

8.3. Computing the radio luminosity function and implications

We computed the luminosity function using Eq. (5) for sources in five different redshift ranges to investigate the evolution of our PS sample. The edges of these redshift ranges were chosen as 0.1, 0.5, 1.0, 1.5, and 3.0, to allow for a sufficient number of sources in each redshift range, while being small enough to probe any evolution in the luminosity functions between redshift ranges. These redshift limits correspond to the lower redshifts for each respective redshift bin used to calculate V_min. We did not evaluate the luminosity function for luminosities below 10²³ W Hz⁻¹, since star-forming galaxies (SFGs) dominate the luminosity function at these low luminosities (e.g. Sabater et al. 2019; Franzen et al. 2021), while for this study we are only interested in the AGN population.

We computed the luminosity function for each redshift range, and for the Master sample and the PS sample. The size of the luminosity bins, Δlog L, is defined separately for the Master and PS sample in each redshift bin (see Table B.1) in order to have the most robust and uniform number of sources as possible in all luminosity bins. For easy comparison, the same luminosity bins are used for the SDSS and photo-z selected samples. We note that we are limited in our analysis by the small number of PS sources in all redshift ranges. This is most notable in the ‘local’ redshift bin (z <  0.1) where the PS source luminosity function could only be evaluated in one luminosity bin.

We estimated the Poissonian counting error on dN(L_j)/dlog L as

$$\begin{array}{cc}\hfill {\sigma }_j& ={\left[\frac{1}{{(\mathrm{\Delta }log{L}_j)}^2}{\displaystyle \sum _{i=1}^N}\frac{1}{V_i^2}\right]}^{1/2}.\hfill \end{array}$$(7)

In luminosity bins with a small number of sources (N <  5), we used the 84% upper and lower confidence limits estimated from Poisson statistics, using Eqs. (9) and (12) of Gehrels (1986). Our resulting luminosity functions are shown in Fig. 11, and tabulated in Appendix B.

Fig. 11. 144 MHz luminosity function, Φ, for the Master sample (black circles) and PS sample (red squares) at different redshifts. The filled symbols show the SDSS selected sample, and the open symbols show the photo-z selected sample. The number of sources in the Master sample and PS sample for each redshift range are shown at the bottom left of their respective plots. The grey curve corresponds to the double power-law model for the AGN local (z ≃ 0.1) luminosity function at 1.4 GHz parametrised by Heckman & Best (2014), converted to 144 MHz using a spectral index −0.7. The 150 MHz luminosity functions for the SF+AGN sample presented by Williams et al. (2018), and the AGN sample presented by Kondapally et al. (2022) are shown as blue squares and green diamonds, respectively. For the z >  1 bins, both the PS and Master sample luminosity functions are not representative of a complete sample due to their optical incompleteness. However, since the optical incompleteness is identical between the PS and Master samples, the offset is a true reflection of evolution between the samples.

The luminosity functions for the SDSS and photo-z selected samples in Fig. 11 agree within ∼1σ throughout redshift and luminosity space, implying that it is the radio incompleteness that dominates the calculation for V_max. However, both optical samples only include roughly one-third of our PS sample, with the rest of the sample likely being too faint for optical detection. Optical incompleteness will impact the luminosity functions at high redshifts, ensuring that the luminosity functions at z >  1 do not represent a complete sample. However, since the optical characteristics of our PS sample and a general sample of radio-loud AGN are similar (Labiano et al. 2007; Nascimento et al. 2022) we can assume that this optical incompleteness will be similar for the PS and Master sample. We conclude that the observed differences between these two samples will not be dominated by selection effects.

In Fig. 11 we also plot the double power-law model for the local (z ≃ 0.1) AGN luminosity function derived at 1.4 GHz by Heckman & Best (2014), but shifted to 144 MHz by using a spectral index of −0.7. We plot the Heckman & Best (2014) model for all redshift ranges to help guide the readers eye. For a complete sample of AGN, the measured luminosity functions in our local redshift range are expected to line up with this model, as was shown by Sabater et al. (2019) for a sample of AGN in LoTSS-DR1. However, as can be seen in the two upper plots of Fig. 11, the local luminosity functions of both the Master sample and PS sample do not line up with this model, but are shifted down by a factor of ∼4.

For the z >  0.5 redshift bins, we plot the Williams et al. (2018) SF+AGN, and the Kondapally et al. (2022) AGN luminosity functions at 150 MHz. For the Kondapally et al. (2022) sample, the redshift bins do not align with the redshift bins in this work, so we used the redshift ranges 0.7 <  z <  1.0, 1.3 <  z <  1.7, and 1.7 <  z <  2.1. For the Williams et al. (2018) sample the edges of the redshift bins are 0.50, 1.00, 1.50, and 2.00.

In the 0.5 <  z <  1.0 redshift bin we also find an offset of a factor ∼4 between our Master sample luminosity function and the luminosity functions presented by Williams et al. (2018) and Kondapally et al. (2022). This offset is present because our first two selection criteria for the Master sample (the source has to be isolated in LoTSS, and the source has to be unresolved in LoTSS), have not been corrected for when computing the luminosity function. We therefore only use a subset of the total AGN population in our detection area to compute the number densities, causing the total luminosity function to be shifted down. However, we note that in the first redshift bin (z <  0.1) our ability to compare the Master sample luminosity function to the model is limited by the small number of sources available (N = 47 for the Master sample, N = 6 for the PS sample).

However, in the two highest redshift bins (z >  1.0) the offset between the Master sample and the literature sample luminosity functions is not present. This can be explained from the fact that at high redshifts most sources will be unresolved. This means that at high redshifts, the Master sample will be close to a complete sample, causing the luminosity function to align with the literature samples.

In the four highest (z ≥ 0.1) redshift ranges, the Master sample luminosity functions flatten and even display a turnover at low luminosities for the redshift ranges > 0.5. This feature is not observed in previous studies of the 150 MHz luminosity function for complete AGN samples at high redshifts (Bonato et al. 2021; Kondapally et al. 2022). We conclude that the turnover at low luminosities in our sample is caused by incompleteness, mostly because of our optical surveys.

However, as this incompleteness is due to survey and selection effects, it is similar for both the Master and PS samples. Therefore, we can use the relative offset between the luminosity functions of the Master sample and the PS sample in each redshift range to estimate the cosmological evolution of our PS sample. These offsets, Δlog(Φ), were found by interpolating the Master sample luminosity functions to the centres of the PS sample luminosity bins, and computing Δlog(Φ)=log(Φ_Master)−log(Φ_PS). Here Φ_Master and Φ_PS are the luminosity functions at the centres of the PS sample luminosity bins for the Master sample and PS sample. The errors on Δlog(Φ) were computed as the quadratic sum of the Poissonian counting errors of the two luminosity functions. The resulting offsets between the Master and PS samples are shown for the SDSS selected sample in Fig. 12, and tabulated in Table B.1. The photo-z selected offsets are not shown in Fig. 12 since the difference in the relative offsets for the two different samples is negligible.

Fig. 12. Difference between the Master sample luminosity function and the PS sample luminosity function, at each PS sample luminosity bin for different redshift ranges.

Figures 11 and 12 show that for the lowest redshift range (z <  0.1) the offset of the PS sample luminosity function to the Master sample luminosity function is smaller than for the four highest redshift ranges. In this lowest redshift bin the offset between the two luminosity functions is close to zero, whereas at higher redshifts the offset is about one order of magnitude. This result is somewhat surprising as it has been generally assumed that the local Universe is underdense in PS sources (e.g. Snellen et al. 2000; Labiano et al. 2007; Kunert-Bajraszewska et al. 2010). However, since we only identified six PS sources with z <  0.1, caution needs to be taken to avoid overinterpreting the data.

More importantly, the offset for all four redshift ranges with z >  0.1 is similar, with the PS luminosity function shifted down by a factor ∼10 compared to the Master sample. As all these redshift ranges only probe sources with L_144 MHz, >  10²⁵ W Hz⁻¹, this informs us that there are approximately ten Master sample sources for every one PS source in a given volume at L_144 MHz, >  10²⁵ W Hz⁻¹.

From Fig. 11 we conclude that the Master sample luminosity function is a factor of ∼4 smaller than that of a complete sample (Williams et al. 2018; Kondapally et al. 2022) in the redshift ranges z <  1.0. From this we can conclude that there are ∼40 radio-loud AGN sources for every PS source for L_144 MHz >  10²⁵ W Hz⁻¹ and z <  1.0. For z >  1.0 this offset reduces to a factor ∼10 between a complete AGN sample and the PS sample.

If we assume the youth scenario for PS sources is correct, PS sources will evolve into Master sample sources on smaller timescales than the redshift ranges probed. Because we observe similar offsets between the Master and PS luminosity functions for L_144 MHz >  10²⁵ W Hz⁻¹ at all redshifts ranges, we can conclude that the birth and death rates of PS sources stay the same between redshifts of 0.1 to 3. This confirms that PS source lifetimes are short compared to cosmological timescales.

The shape of the PS luminosity function also remains approximately constant relative to the Master sample and to the complete AGN populations presented by Williams et al. (2018) and Kondapally et al. (2022). A similar luminosity function for PS sources and large-scale radio-loud AGN implies that PS sources evolve into large-scale radio sources as their cosmological evolution is the same. This result lends strong support to the youth model for PS sources. A constant birth rate for PS sources across z >  0.1 also implies that a fraction of radio-loud AGN fade to the point that they will no longer be identified as high luminosity (L_144 MHz, >  10²⁵ W Hz⁻¹). If this were not the case, there would be an overabundance of radio-loud AGN in the local Universe relative to the distant Universe.

However, we note that the offset between the Master and PS sample luminosity functions is less significant for sources with L_144 MHz <  10²⁵ W Hz⁻¹. This suggests that there are differences between PS sources with low luminosities and those with high luminosities. It is possible that at these low luminosities we are potentially also detecting frustrated PS sources, which have such a low jet power that the jet may not ever penetrate the interstellar medium of its host galaxy. In this case the birth and death rates of the young population of PS sources would be the same as in the more distant universe, but the local Universe has an additional number of low-luminosity frustrated PS sources that we are not sensitive to at high redshifts in our flux density-limited surveys. However, additional detections over a larger detection area of local low-luminosity (L_144 MHz <  10²⁵ W Hz⁻¹) PS sources are necessary to confirm this hypothesis.

9. Conclusions

In this study we have identified 373 PS sources with spectral peaks around 144 MHz using the LOFAR surveys LoTSS and LoLSS, as well as NVSS. A source was identified as a PS source when the power-law fit between its LoLSS and LoTSS flux density measurements was positive, while the power law between the LoTSS and NVSS flux densities was negative. For comparison purposes, we also defined a Master sample of unresolved radio sources that made no assumption about the spectral shape of the source, and is the sample from which the PS sample was drawn.

Our PS sample has increased the number of known PS sources around the HETDEX area by a factor of 50. Using the LoTSS ancillary optical catalogue, we are able to identify the redshift for 138 out of 373 PS sources, of which 54 were spectroscopic redshifts. The 5 GHz radio luminosity distribution computed for our PS sample shows that we have identified the lowest average radio power PS sample to date, by roughly an order of magnitude.

Using this sample, we investigated the evolution of our PS sample using HERG and LERG classification, source counts, and luminosity functions. We found the following:

• HERG–LERG classification: Using optical line emission criteria we were able to identify 25 PS sources as HERGs (20 sources) or LERGs (5 sources). When we take redshift bias into account, we find a one-to-one ratio of HERGs and LERGs in our PS sample, showing our PS sources are predominantly HERGs compared to a general AGN population. Since HERGs are a quickly evolving population, this suggests our PS sample will be dominated by quickly evolving sources.

• Source counts: We evaluated the Euclidean normalised source counts at 144 MHz for our PS sample, as well as two known PS samples (Snellen et al. 1998; Callingham et al. 2017). When compared to the source counts for a general sample of AGN, the PS source counts are scaled down by a factor of 40. From this we conclude that the lifetime of PS sources is about 40 times shorter than that of low-frequency radio-loud AGN, assuming that the luminosity function of both populations are the same.

• Luminosity function: With the redshift information available, we computed the radio luminosity function for our PS sample. This is the first time a PS luminosity function has been produced for a homogeneously selected PS sample. We demonstrate that for redshifts > 0.1 and for sources with L_144 MHz ≳ 10²⁵ W Hz⁻¹, the offset between the PS luminosity function and that of the Master sample remains constant. We interpret this as strong evidence that these high-luminosity PS sources evolve into large-scale radio-loud AGN. This conclusion also implies that there is one PS source for about every ten unresolved high-luminosity radio-loud AGN, and the rate at which PS sources enter the later population is roughly consistent with the rate at which the large AGN fade to lower luminosities. If this were not the case, there would be an overabundance of PS sources relative to their evolved counterparts. However, in the local (z <  0.1) Universe, we note that the offset between the PS and Master sample luminosities is smaller at luminosities less than 10²⁵ W Hz⁻¹. We interpret this as the potential addition of frustrated sources in the population, which do not have the power to evolve into large-scale radio-loud AGN. We do not see this population at higher redshifts as our flux-limited surveys are not sensitive to these low-luminosity sources at those distances.

We conclude that our HERG–LERG classification, source counts, and luminosity function analyses all indicate that our population of PS sources is a quickly evolving radio-loud AGN population. We suggest that this provides strong support that the youth scenario applies to the majority of our PS sources with low-frequency luminosities ≳10²⁵ W Hz⁻¹. In addition, the relative lifetimes and abundances of PS sources found using the source counts and luminosity functions are in agreement with each other, further supporting this hypothesis. However, we note that the surplus of PS sources at luminosities lower than this still allows the frustration hypothesis of PS sources to apply.

To test if these sources are in fact frustrated, precise broad-band spectral modelling is required. Furthermore, a larger sample of PS samples, which will be derivable once the LoLSS and LoTSS all-sky surveys are completed, is needed to test the robustness of these results. In particular, accurate and complete spectroscopic redshift information from WEAVE-LOFAR (Smith et al. 2016) will be especially important to reduce any optical completeness issues in our analysis. Our conclusions also imply that approximately 1 in 40 sources in a flux-density limited survey will be a PS source.

Future very long baseline interferometry (VLBI) observations with the LOFAR international baselines are needed to identify specific high-resolution structures and morphologies for the PS sources in our sample. With this, the environments and symmetries of our sample of PS sources can be investigated to confirm our hypothesis that the majority of these sources are not frustrated, but will instead grow to be large radio galaxies. Furthermore, these morphological studies could identify transient short-lived PS sources by searching for multi-epoch activity.

Detailed studies of the optical hosts of our sample of PS sources are needed to investigate the optical nature of these sources. This can give insight into whether PS sources with quasar and galaxy hosts have different characteristics, which could be suggestive of a different spectral turnover mechanism.

To improve the selection of PS sources in this detection area, a better sensitivity in our low frequency survey is needed as this is the most limiting factor in identifying PS sources in this study. The final LoLSS data release will enable us to improve this sensitivity, as it will be approximately five times deeper, and have a resolution that is about three times higher than that in the preliminary data release used in this study (de Gasperin et al., in prep.). The methods of selecting PS sources in LOFAR surveys outlined in this study can then be easily applied to future data releases with these improvements as well as a larger sky coverage of LoTSS and LoLSS to identify a large number of new PS sources. Finally, our empirically derived luminosity functions need to be tested against evolutionary models for radio-loud AGN (Bicknell et al. 2018) to ensure that the models are consistent with the evolution we are suggesting.

Acknowledgments

We thank E. M. Sadler (CSIRO/University of Sydney) and R. D. Ekers (CSIRO) for their helpful suggestions on this work. JRC thanks the Nederlandse Organisatie voor Wetenschappelijk Onderzoek (NWO) for support via the Talent Programme Veni grant. KJD acknowledges funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 892117 (HIZRAD). MJH acknowledges support from the UK STFC [ST/V000624/1]. 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 that are collectively operated by the ILT foundation under a joint scientific policy. The ILT resources have benefited from the following recent major funding sources: CNRS-INSU, Observatoire de Paris and Université d’Orléans, France; BMBF, MIWF-NRW, MPG, Germany; Science Foundation Ireland (SFI), Department of Business, Enterprise and Innovation (DBEI), Ireland; NWO, The Netherlands; The Science and Technology Facilities Council, UK; Ministry of Science and Higher Education, Poland; The Istituto Nazionale di Astrofisica (INAF), Italy. This research made use of the Dutch national e-infrastructure with support of the SURF Cooperative (e-infra 180169) and the LOFAR e-infra group. The Jülich LOFAR Long Term Archive and the German LOFAR network are both coordinated and operated by the Jülich Supercomputing Centre (JSC), and computing resources on the supercomputer JUWELS at JSC were provided by the Gauss Centre for Supercomputing e.V. (grant CHTB00) through the John von Neumann Institute for Computing (NIC). This research made use of the University of Hertfordshire high-performance computing facility and the LOFAR-UK computing facility located at the University of Hertfordshire and supported by STFC [ST/P000096/1], and of the Italian LOFAR IT computing infrastructure supported and operated by INAF, and by the Physics Department of Turin university (under an agreement with Consorzio Interuniversitario per la Fisica Spaziale) at the C3S Supercomputing Centre, Italy. This research has made use of the NASA/IPAC Extragalactic Database (NED), which is funded by the National Aeronautics and Space Administration and operated by the California Institute of Technology, the VizieR catalogue access tool, CDS, Strasbourg, France, and the “K-corrections calculator” service available at http://kcor.sai.msu.ru/. This research has also made use of “Aladin sky atlas” developed at CDS, Strasbourg Observatory, France, TOPCAT (Taylor 2005), NASA’s Astrophysics Data System bibliographic services, and Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration 2018). This work has also made use of the Python packages Numpy (Harris et al. 2020), SciPy (Virtanen et al. 2020), and Matplotlib (Hunter 2007).

References

Appendix A: List of column headings for table of PS sample

The numbers, names, and units for the columns of the table presenting our PS sample are outlined below.

Appendix B: Peaked-spectrum radio luminosity function

All Tables

All Figures

Fig. 1. Limiting flux densities and frequencies for major radio surveys. The limiting flux density for a survey is given by the faintest catalogued source flux density. The GLEAM survey is represented as a black line since it has variable limiting flux densities over its range of observing frequencies. The orange dashed curve represents the observational limit for the sample presented by O’Dea (1998). This is given by an SSA spectrum of a PS source with a peak at 750 MHz and an observed peak flux density of 300 mJy. The green dash-dotted curve represents the observational limit for the sample presented by Callingham et al. (2017). This is given by an SSA spectrum of a PS source with a peak at 190 MHz, with a faintest observed peak flux density of 160 mJy. The blue curve represents the theoretical observational limit of PS sources in this study. This is given by an SSA spectrum of a PS source with a peak at 100 MHz, with a faintest peak flux density of 80 mJy, based on the limiting flux densities of LoLSS-PDR, LoTSS, and NVSS. This figure also illustrates that in this study LoLSS is the survey that dictates the limiting flux density for identifying PS sources. The following plotted surveys were not previously mentioned: Cambridge 7C (Hales et al. 2007) survey, Westerbork Northern Sky Survey (WENSS; Rengelink et al. 1997), Texas Survey (TXS; Douglas et al. 1996), Molonglo Reference Catalogue (MRC; Large et al. 1991), Parkes (PKS; Wright & Otrupcek 1990) survey, Sydney University Molonglo Sky Survey (SUMSS; Mauch et al. 2003), MIT-Green Bank 5 GHz (87GB; Gregory & Condon 1991) survey, Parkes-MIT-NRAO (PMN; Wright et al. 1994) survey, Australia Telescope 20 GHz (AT20G; Murphy et al. 2010) survey, and the Rapid ASKAP Continuum Survey (RACS; Hale et al. 2021).

In the text

Fig. 2. Radio colour-colour diagram for the 9768 LoTSS sources left in our Master sample. αlow represents the spectral index between the LoLSS (54 MHz) and LoTSS (144 MHz) flux density points. αhigh represents the spectral index between the LoTSS and NVSS flux density points. To better illustrate the large number of sources around the median of αlow and αhigh, given by (αlow, αhigh) =(−0.6 ± 0.2, −0.8 ± 0.1), contours and a density map are plotted in this region. The contour levels represent 15, 75, 145, and 265 sources. The colours in the underlying density map show the number of sources for each pixel. The number of sources corresponding to each corresponding shading are illustrated by the colour bar at the top left of the plot. The grey lines show spectral indices of zero to illustrate the four quadrants of the plot. The dashed red line represents a one-to-one relation between αlow and αhigh. The blue line at αlow = 0.1 illustrates the selection limit used to identify PS sources in this study. In the corner of each quadrant the shape of a typical spectrum in that quadrant is shown in grey. To avoid confusion, individual error bars are not plotted, but in the top left of the plot the median error bar size is shown. The histograms at the top and right of the diagram illustrate the distributions of αlow and αhigh, respectively. These distributions have been normalised to the maximum value in the distribution. In these histograms Gaussian fits to these distributions are overplotted in red, and the dashed black lines show the median values for these distributions.

In the text

Fig. 3. Example spectra of a soft sample PS source (top) and a hard sample PS source (bottom). The flux densities of LoLSS (blue square), VLSSr (red circle), LoTSS-inband spectra (purple diamonds), LoTSS (green left-pointing triangle), TGSS (brown right-pointing triangle), NVSS (indigo downward-pointing triangle), and FIRST (pink upward-pointing triangle) are plotted where available. The black curve represents the generic curved model from Eq. (3), which was fitted to the LoLSS, LoTSS, VLSSr, TGSS, and NVSS data points. The power-law fits used to determine αlow and αhigh are shown in orange.

In the text

	Fig. 4. Redshift distribution of the 138 PS sources with available redshift data. The red histogram is the distribution of the spectroscopic redshifts available for 54 sources. The black histogram is the distribution of the best redshifts, from the combination of spectroscopic and photometric redshifts.
In the text

Fig. 5. Spectra of the two known GPS, CSO, and HFP sources in the Master sample. The symbols represent data from the same surveys as in Fig. 3. The grey crosses show the data presented by Marvil et al. (2015), the dark green circles the data from the TXS survey (Douglas et al. 1996), and the teal hexagons the data presented by Tremblay et al. (2016). The black curve in the left plot shows the generic curved model from Eq. (3).

In the text

Fig. 6. Distribution of the 5 GHz radio power of the Master sample and PS sample. The upper histogram shows the distribution of the 5 GHz radio power for the 3303 sources from the Master sample that have redshift information available in red. The black histogram in the bottom plot shows the distribution of the 5 GHz radio power of the 138 sources from the PS sample for which redshift data is available. The median redshifts of these distributions are plotted as dashed lines in their respective colours. For the total sample, the median P5 GHz value and the range of the 16th and 84th percentiles of the distribution is given by $25.0_{-0.7}^{+1.0}$ log10 W Hz−1, and for the PS source sample this is given by $25.7_{-1.1}^{+0.9}$ log10 W Hz−1.

In the text

Fig. 7. Distribution of the 5 GHz radio power for different samples of PS samples. The distribution of the 5 GHz radio power for the 138 sources from the PS sample for which redshift data is available is shown in black (upper histogram). The green, blue, and orange histograms in the lower plot represent the 5 GHz radio power of the Snellen et al. (1998) PS sample, the O’Dea (1998) PS sample, and the Callingham et al. (2017) PS sample, respectively. The median redshifts of all distributions are plotted as dashed lines in their respective colours, in both the upper and lower plots. The median P5 GHz value and the range of the 16th and 84th percentiles of the distributions of the samples in this study, Snellen et al. (1998), O’Dea (1998), and Callingham et al. (2017), are given by $25.7_{-1.1}^{+0.9}$, $26.2_{-0.4}^{+0.5}$, $27.6_{-1.1}^{+0.7}$, and $27.2_{-1.2}^{+0.6}$ (log10 W Hz−1), respectively.

In the text

Fig. 8. 5 GHz radio power against redshift for the 138 sources from the PS sample for which redshift data is available (black circles). The green triangles, blue squares, and orange diamonds represent the PS samples of O’Dea (1998), Snellen et al. (1998), and Callingham et al. (2017), respectively. The dashed grey line corresponds to the 5 GHz luminosity limit for a source that has a peak flux of 44 mJy at 144 MHz, corresponding to the 90% completeness limit of the PS source selection. To compute this 5 GHz luminosity limit, a median spectral index of −0.87, corresponding to the 90% limit of αhigh, was used.

In the text

Fig. 9. [O III] line luminosity vs 1.4 GHz radio luminosity of sources classified as HERGs and LERGs, and unclassified PS sources. The red diamonds, orange squares, and blue triangles indicate PS HERGs, PS LERGs and unclassified PS sources, respectively. The [O III] line luminosities vs 1.4 GHz radio luminosities of the sample of HERG and LERG sources presented by Best & Heckman (2012) are plotted as red and black points, respectively. In order to avoid overcrowding of the figure, only 20% of the LERGs in the Best & Heckman (2012) sample are plotted. The solid black line indicates the lower limit to the distribution of HERGs proposed by Best & Heckman (2012).

In the text

Fig. 10. Euclidean normalised differential source counts for different samples of PS sources. The black squares show our sample of PS sources. The blue circles show the Snellen et al. (1998) PS sample converted to 144 MHz. The orange diamonds show the Callingham et al. (2017) PS sample. For reference, the source counts for the LoTSS Deep Field derived by Mandal et al. (2021) are plotted as red triangles. The solid grey line shows the model for 150 MHz AGN source counts proposed by Massardi et al. (2010). The dashed grey line shows the Massardi et al. (2010) model scaled down by a factor of 40, which is consistent with the PS source data.

In the text

Fig. 11. 144 MHz luminosity function, Φ, for the Master sample (black circles) and PS sample (red squares) at different redshifts. The filled symbols show the SDSS selected sample, and the open symbols show the photo-z selected sample. The number of sources in the Master sample and PS sample for each redshift range are shown at the bottom left of their respective plots. The grey curve corresponds to the double power-law model for the AGN local (z ≃ 0.1) luminosity function at 1.4 GHz parametrised by Heckman & Best (2014), converted to 144 MHz using a spectral index −0.7. The 150 MHz luminosity functions for the SF+AGN sample presented by Williams et al. (2018), and the AGN sample presented by Kondapally et al. (2022) are shown as blue squares and green diamonds, respectively. For the z >  1 bins, both the PS and Master sample luminosity functions are not representative of a complete sample due to their optical incompleteness. However, since the optical incompleteness is identical between the PS and Master samples, the offset is a true reflection of evolution between the samples.

In the text

	Fig. 12. Difference between the Master sample luminosity function and the PS sample luminosity function, at each PS sample luminosity bin for different redshift ranges.
In the text