# Characterising the Extended Morphologies of BL Lacertae Objects at 144 MHz with LOFAR

Seán Mooney, Francesco Massaro, John Quinn, Alessandro Capetti, Ranieri D. Baldi, Gülay Gürkan, Martin J. Hardcastle, Cathy Horellou, Beatriz Mingo, Raffaella Morganti, et al.

► To cite this version:

Seán Mooney, Francesco Massaro, John Quinn, Alessandro Capetti, Ranieri D. Baldi, et al. Char-acterising the Extended Morphologies of BL Lacertae Objects at 144 MHz with LOFAR. The Astro-physical Journal Supplement, 2021, 257 (2), pp.30. 10.3847/1538-4365/ac1c0b . hal-03354023

HAL Id: hal-03354023

https://hal.science/hal-03354023v1

Submitted on 15 Jul 2025

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Characterising the extended morphologies of BL Lacs at 144 MHz with LOFAR

SEÁN MOONEY,¹ FRANCESCO MASSARO,² JOHN QUINN,¹ ALESSANDRO CAPETTI,³ RANIERI D. BALDI,⁴,⁵ GÜLAY GÜRKAN,⁶

MARTIN J. HARDCASTLE,⁷ CATHY HORELLOU,⁸ BEATRIZ MINGO,⁹ RAFFAELLA MORGANTI,¹⁰,¹¹ SHANE O'SULLIVAN,¹²

URSZULA PAJDOSZ-ŚMIERCIAK,¹³ MAMTA PANDEY-POMMIER,¹⁴ AND HUUB RÖTTGERING¹⁵

¹School of Physics, University College Dublin, Belfield, Dublin 4, Republic of Ireland

²Dipartimento di Fisica, Università degli Studi di Torino, via Pietro Giuria 1, 10125 Torino, Italy

³INAF - Osservatorio Astrofisico di Torino, via Osservatorio 20, Pino Torinese, 10025 Italy

⁴Istituto di Radioastronomia - INAF, Via P. Gobetti 101, I-40129 Bologna, Italy

⁵School of Physics and Astronomy, University of Southampton, Southampton, SO17 1BJ, UK

⁶Thüringer Landessternwarte, Sternwarte 5, D-07778 Tautenburg, Germany

⁷Centre for Astrophysics Research, School of Physics, Astronomy & Mathematics, University of Hertfordshire, College Lane, Hatfield AL10 9AB, United Kingdom

⁸Chalmers University of Technology, Department of Space, Earth and Environment, Onsala Space Observatory, SE-43992 Onsala, Sweden

⁹School of Physical Sciences, The Open University, Walton Hall, Milton Keynes, MK7 6AA, United Kingdom

¹⁰ASTRON, the Netherlands Institute for Radio Astronomy, Postbus 2, 7990 AA Dwingeloo, the Netherlands

¹¹Kapteyn Astronomical Institute, University of Groningen, PO Box 800, 9700 AV Groningen, the Netherlands

¹²School of Physical Sciences & CfAR, Dublin City University, Glasnevin, Dublin 9, Republic of Ireland

¹³Astronomical Observatory of the Jagiellonian University, Orla 171, 30-244 Cracow, Poland

¹⁴USN, Station de Radioastronomie de Nançay Observatoire de Paris route de Souesmes 18330 Nançay France & Univ Lyon, Univ Lyon1, Ens de Lyon, CNRS, Centre de Recherche Astrophysique de Lyon UMR5574, 9 av Charles André F- 69230, Saint-Genis-Laval, France

¹⁵Leiden Observatory, Leiden University, PO Box 9513, NL-2300 RA Leiden, the Netherlands

(Received December 22, 2020; Revised July 11, 2021; Accepted August 1, 2021)

Submitted to The Astrophysical Journal Supplement Series

ABSTRACT

We present a morphological and spectral study of a sample of 99 BL Lacs using the LOFAR Two-Metre Sky Survey Second Data Release (LDR2). Extended emission has been identified at GHz frequencies around BL Lacs, but with LDR2 it is now possible to systematically study their morphologies at 144 MHz, where more diffuse emission is expected. LDR2 reveals the presence of extended radio structures around 66/99 of the BL Lac nuclei, with angular extents ranging up to 115'', corresponding to spatial extents of 410 kpc. The extended emission is likely to be both unbeamed diffuse emission and beamed emission associated with relativistic bulk motion in jets. The spatial extents and luminosities of the extended emission are consistent with the AGN unification scheme where BL Lacs correspond to low-excitation radio galaxies with the jet axis aligned along the line-of-sight. While extended emission is detected around the majority of BL Lacs, the median 144-1400 MHz spectral index and core dominance at 144 MHz indicate that the core component contributes ~42% on average to the total low-frequency flux density. A stronger correlation was found between the 144 MHz core flux density and the γ-ray photon flux (r = 0.69) compared to the 144 MHz extended flux density and the γ-ray photon flux (r = 0.42). This suggests that the radio-to-γ-ray connection weakens at low radio frequencies because the population of particles that give rise to the γ-ray flux are distinct from the electrons producing the diffuse synchrotron emission associated with spatially-extended features.

Keywords: BL Lacertae objects: general — catalogs — surveys

Corresponding author: Seán Mooney

sean.mooney@ucdconnect.ie

1. INTRODUCTION
Blazars, interpreted to be radio-loud active galactic nuclei (AGN) with relativistic jets oriented at small angles along the line of sight, are the largest known population of sources in the extragalactic γ-ray sky. The small inclination angles result in strong Doppler beaming effects that significantly boost the observed flux of the jet moving towards the Earth and deboost the flux of the counter jet. The jets eventually dissipate by either decollimation or by depositing their kinetic energy in a terminal shock. These processes produce spatially-extended diffuse unbeamed emission which is potentially observable at low frequencies in addition to the beamed emission (as indicated by simulations in, e.g., Hardcastle & Krause 2014; Massaglia et al. 2016).

Blazars are classified as either BL Lacs, which typically exhibit narrow emission or absorption lines in their optical spectra, or flat-spectrum radio quasars, which have broad emission lines typical of quasars (Landt et al. 2001). Hereinafter we adopt the nomenclature of the blazar catalogue Roma-BZCAT (Massaro et al. 2015), labelling the former class of BL Lacs as BZBs and the latter as BZQs. BZBs are the focus of this study their emission is more dominated by non-thermal components when compared to BZQs, which can have features such as a dusty torus and a big blue bump in their broad-band spectral energy distribution (Wilkes 2004).

According to the unification scheme of radio-loud AGN (Urry & Padovani 1995), source orientation with respect to the line-of-sight is key to explaining observational differences between blazars and radio galaxies, with low-excitation and high-excitation radio galaxies (LERGs and HERGs; Hine & Longair 1979) believed to be the parent populations of BZBs and BZQs, respectively. A strong prediction of the unification theory is that at low frequencies the morphologies of BZBs should conform with LERGs viewed at small angles to the jet (see the review by Best & Heckman 2012). For highly aligned sources the low-frequency radio emission is expected to consist of an unresolved beamed core surrounded by extended diffuse emission related to the jet terminus, while for jets with non-zero inclination angles, in addition to the core, there may be emission associated with the large-scale jet and potentially the counterjet including both termination regions.

While there have been GHz studies that detected extended emission (e.g Ulvestad et al. 1983; Antonucci & Ulvestad 1985; Laurent-Muehleisen et al. 1993), to date, the BZBs morphologies at $\le 144$ MHz have not been systematically explored because no wide-field survey has had the spatial resolution, dynamic range, and sensitivity required to resolve all the components of the sources. Information about the beamed and extended

emission has been inferred from spectral studies, where the beamed emission is expected to follow a power law with a flat spectral index (i.e. $-0.5 \le \alpha \le 0.5$, where $S_\nu \propto \nu^\alpha$ throughout) and the diffuse emission is believed to have a power law spectral index more typical of optically-thin synchrotron emission ($\alpha \approx -0.8$). The flat-spectrum beamed component dominates the radio spectrum above $\sim 1$ GHz (e.g. Healey et al. 2007) and in the last decade it has been shown that the spectral indices of blazars are generally flat below $\sim 1$ GHz as well (Massaro et al. 2013b). Indeed, this characteristic flat spectrum at low frequencies has been successfully leveraged in identifying γ-ray sources (Massaro et al. 2014).

Blazars tend to be core-dominated in the GHz regime, where the core dominance is a proxy for the beaming factor and hence the inclination angle (Antonucci & Ulvestad 1985; Perlman & Stocke 1993). In the $\sim 100$ MHz regime, it has been inferred that blazars are core-dominated, on the basis of the flat spectral indices (Massaro et al. 2014), but a direct measurement of the core dominances for a sample of blazars has yet to be ascertained. Giroletti et al. (2016) and d'Antonio et al. (2019) estimated the low frequency core dominances of blazars using a spectral decomposition technique, relying on assumptions about the spectral indices of the core and extended components respectively.

Given the broadband nature of blazar spectral energy distributions, a comprehensive understanding of their low frequency spectra is paramount to obtain a complete view of blazars. This is highlighted by the radio-to-γ-ray connection, where there is a well-established link between the $\gtrsim 1$ GHz flux and the γ-ray flux, spanning $\sim 17$ decades of energy (e.g. Ackermann et al. 2011). From this, it is inferred that the radio and γ-ray emission are produced by the same population of relativistic electrons. The relationship between the flux density at hundreds of MHz and the γ-ray flux is less clear however (Giroletti et al. 2016; Mooney et al. 2019), possibly because the γ-ray-emitting particles are distinct from the electrons that give rise to the large-scale diffuse emission that is expected to be prevalent at low frequencies.

LOFAR is conducting a high angular resolution, highly-sensitive survey of the northern hemisphere sky at 144 MHz (LoTSS; the LOFAR Two-Metre Sky Survey; Shimwell et al. 2017), with more than 21% of the sky observed to date. The LoTSS Second Data Release (LDR2) presents a unique opportunity to make improved low-frequency radio measurements of blazars and, for the first time, investigate the morphology of blazars at 144 MHz. We present LDR2 data for a sample of BZBs and measure the spatial extents, core dominances, and spectral indices to characterise the diffuse
emission. The radio-to-$\gamma$-ray connection is investigated
for the low frequency core and extended components
separately in order to understand why this trend, which
is clear at GHz frequencies, tends not to be significant
at 144 MHz.

This paper is organised as follows: In § 2 the sample
selection is outlined, in § 3 the analysis is detailed, in
§ 4 the results are presented, in § 5 the findings are
discussed, and in § 6 a summary is provided. A ΛCDM
cosmological model is used in this paper with h = 0.70
(The LIGO Scientific Collaboration et al. 2017), Ω_m =
0.26, and Ω_Λ = 0.74, where H₀ = 100h km s⁻<fnref n="1" /> Mpc⁻<fnref n="1" /> is
the Hubble constant.

2. DATASETS AND SAMPLE SELECTION

2.1. *LOFAR dataset*

The Low Frequency Array (LOFAR; van Haarlem
et al. 2013) is a radio interferometer with stations lo-
cated throughout Europe. One goal of the LOFAR Sur-
veys Key Science Project (SKSP) is to map the north-
ern hemisphere sky at 120–168 MHz; this is known as
the LOFAR Two Metre Sky Survey (LoTSS; Shimwell
et al. 2017). In this paper we use data from a subset of
LoTSS Second Data Release (LDR2).

LDR2 data were processed with the SKSP pipeline<fnref n="1" />
version 2.2 as described by Shimwell et al. (2019). In the
pipeline, direction-dependent calibration was carried out
using KILLMS (Tasse 2014; Smirnov & Tasse 2015) and
imaging was done using DDFACET (Tasse et al. 2018),
where a novel self-calibration strategy was employed (§ 5
of Shimwell et al. 2019; Tasse et al. 2020). With this
strategy, extended emission that is undetected in early
cycles of self-calibration is less likely to modelled out,
and more attention is paid to properly deconvolving that
emission. This more complex self-calibration process im-
proved the sensitivity to extended emission with respect
to LoTSS DR1 (Tasse et al. in prep.). Source catalogues
were extracted by the pipeline using PYBDSF (Mohan
& Rafferty 2015).

LDR2 observations are in progress. We use the LDR2
data that were available as of 2020 February 01, amount-
ing to 4240 deg<fnref n="2" /> of sky coverage, which is 21% of the
northern hemisphere sky (Fig. 1). This LDR2 subset
encompasses the publicly-available LDR1 and the cata-
logue contains 3.6 M sources. The resolution is 6" with
a pixel size of 1.5" and the median RMS noise level is
~70 µJy beam⁻<fnref n="1" />.

2.2. *Ancillary datasets*

<figure><img src="image_2.png" /></figure>

**Figure 1.** LoTSS DR2 multi-order coverage (blue) overlaid on the Fermi-LAT 0.3–1 GeV HEALPix survey (greyscale; Atwood et al. 2009). Aitoff projection is used. LoTSS DR2 covers 27% of the northern hemisphere sky as of 2020 February 01. The goal of LoTSS is to survey the entire Northern Hemisphere sky.

The Faint Images of the Radio Sky at Twenty-cm
(FIRST; Becker et al. 1995)<fnref n="2" /> survey was used, where
available, to visually check for the presence of extended
emission at 1.4 GHz. FIRST data were only unavailable
for BL Lacs lying outside its footprint. These data have
a comparable resolution to LDR2 (~5").

The 144–1400 MHz spectral indices are computed us-
ing the NRAO VLA Sky Survey (NVSS; Condon et al.
1998)<fnref n="3" /> with LDR2. NVSS is preferred to FIRST because
for extended sources there can be a deficit in the FIRST
total flux densities compared to NVSS, particularly for
sources with total fluxes in the 2–20 mJy range. This is
due to extended flux being undetectable given the lack
of short baselines in the FIRST UV coverage (Helfand
et al. 2015). However, the FIRST resolution is more
favourable than the 45" resolution offered by NVSS, so
FIRST was used over NVSS to look for large-scale jets
at 1.4 GHz.

The difference between the spatial resolutions of
NVSS and LoTSS are not likely to be problematic (de
Gasperin et al. 2018). All BL Lacs are mainly point-like
in FIRST and NVSS (i.e. consistent with the beam size)
and flux densities were computed by integrating over the
beam. We also checked the flux densities against those
reported in Giroletti et al. (2016) and Massaro et al.
(2013a) to verify that no flux was missing (where NVSS
had similar beam sizes to these surveys).

We used *i*-band data from the Panoramic Survey Tele-
scope and Rapid Response System Data Release 1 (Pan-
STARRS1; Chambers et al. 2016)<fnref n="4" /> to check that, for
each BZB, the emission at detected MHz frequencies is
likely associated with its optical counterpart correspond-

<footnote n="1">¹ https://github.com/mhardcastle/ddf-pipeline</footnote>

<footnote n="2">² http://sundog.stsci.edu/</footnote>

<footnote n="3">³ https://www.cvnrao.edu/nvss</footnote>

<footnote n="4">⁴ https://panstarrs.stsci.edu</footnote>
ing to the coordinates reported in Roma-BZCAT. Both FIRST and Pan-STARRS1 images for the sample can be accessed through the online appendices.

## 2.3. BZB sample selection

We started with the Roma-BZCAT 5th edition (Massaro et al. 2015) blazar catalogue and took the following steps to identify the BZBs that are in LDR2.

* We first considered all 3561 blazars in Roma-BZCAT. We added to this 331 newly identified blazars from a recent optical spectroscopic campaign that will feature in the next release of Roma-BZCAT (Landoni et al. 2015; Ricci et al. 2015; Álvarez Crespo et al. 2016; Massaro et al. 2016).

* We selected only the BZBs from this enhanced version of Roma-BZCAT that have a counterpart in LDR2 within a few arcseconds, where the criterion depended on the spatial extent of the source. The Roma-BZCAT position, which refers to the central engine, and the LDR2 position, which relates to the centroid of the radio emission, are not necessarily aligned. The majority (55%) of BZBs had a crossmatch separation of < 1" and 85% of sources had a separation of < 6", which is the LDR2 beam width. The largest crossmatch separations pertain to highly extended sources. A counterpart was identified for all known BZBs within the LDR2 footprint and all crossmatches were visually confirmed. We then excluded BZBs that have uncertain redshift estimates, and 99 BZBs remained.

* We also identified the BZBs in the sample that have a γ-ray counterpart by crossmatching the Roma-BZCAT positions with the Fourth Catalog of AGN detected by the Fermi Large Area Telescope (4LAC; Abdollahi et al. 2020; Ajello et al. 2020) within 5", where we used the position of the associated counterpart in 4LAC. A total of 53 BZBs are γ-ray detected in 4LAC.

* Using the Roma-BZCAT positions, the BZB sample was cross-matched with FIRST (Becker et al. 1995) within 5". There were FIRST counterparts to 83/99 BZBs.

* Again using the Roma-BZCAT positions, the BZB sample was cross-matched with Pan-STARRS1 (Chambers et al. 2016) within 5", and a counterpart was identified in all cases.

The final sample consists of 99 BZBs; 91 are from Roma-BZCAT (5th edition) and 8 are from recent opti-

<table><caption>**Table 1.** Overview of sources in our sample.</caption><thead><tr><th></th><th>Total</th><th>Selected*</th></tr></thead><tbody><tr><td>Roma-BZCAT v5.0</td><td>3561</td><td>91</td></tr><tr><td>Optical campaigns of unassociated γ-ray sources</td><td>331</td><td>8</td></tr><tr><td>Total</td><td></td><td>99**</td></tr></tbody></table>

NOTE— *BZBs in the LDR2 footprint with a reliable redshift estimate.

**All 99 have an optical counterpart in Pan-STARRS, 83 have a radio counterpart in FIRST, and 53 have a γ-ray counterpart in 4LAC.

cal spectroscopic campaigns (Paiano et al. 2019; Peña-Herazo et al. 2019; de Menezes et al. 2019; Massaro et al. 2016). An overview of the selection is given in Table 1.

For the analysis presented in § 3, we used the data in the LDR2 catalogue. However, for J1231+3711, the source extraction software (PyBDSF) failed to accurately characterise the widespread diffuse emission that surrounds the compact core. Therefore, instead of using the LDR2 catalogue values for J1231+3711, we reran PyBDSF with a reduced threshold for the flux detection, and used these recalculated values for our analysis.

## 3. DATA ANALYSIS

The images and the catalogue we constructed based on LDR2 were the starting point of this analysis. The extents, core dominances, and spectral indices of the BZBs were calculated as follows.

### 3.1. Angular and spatial extents at 144 MHz

We used the empirically-derived equation in § 3.1 of Shimwell et al. (2019) to determine which BZBs in the sample are resolved (rBZBs) and which BZBs are unresolved (uBZBs). This equation defines an envelope separating resolved and unresolved sources,

$$ \mathcal{R} = \frac{S_{\text{total}}}{S_{\text{peak}}} - 1.25 - 3.1 \left( \frac{S_{\text{peak}}}{\sigma} \right)^{-0.53} \quad (1) $$

where $S_{\text{total}}$ is the total flux density of a source, $S_{\text{peak}}$ is the peak flux, and $\sigma$ is the RMS noise as per the LDR2 catalogue. If $\mathcal{R} > 0$, then the source is classified as an rBZB, otherwise the source is a uBZB. In total, there are 66/99 rBZBs and 33/99 uBZBs in the sample.

The intrinsic angular extent, $\Phi$, of each BZB was taken to be the FWHM of the major axis of the source after Gaussian deconvolution of the point spread function (PSF). For BZBs that consist of multiple Gaussians
grouped, the FWHM in the catalogue refers to the com-
bination of the major axes of the individual components
combined using moment analysis. The median intrinsic
angular extent for the rBZBs is (13 ± 2)" and the me-
dian upper limit on the intrinsic angular extents for the
uBZBs is (4 ± 1)".

We computed the spatial extent, *D*, of each BZB using
the intrinsic angular extent and the redshifts, assuming
the flat cosmology described in § 1 (Wright 2006). For
uBZBs, the angular and spatial extents are reported as
upper limits.

3.2. *Core dominances at 144 MHz*

The core dominance is the ratio of fluxes from the core and extended regions (Giovannini et al. 1990), and the distribution of relative intensities of the core to the extended emission should reflect the distribution of Doppler factors (Antonucci & Ulvestad 1985; Perlman & Stocke 1993). Following on from studies such as that of Morganti et al. (1997), we define the core dominance at 144 MHz, $\rho_{144}$, as

$$
\rho_{144} = \frac{S_{\text{core}}}{S_{\text{ext}}} \qquad (2)
$$

where $S_{\text{core}}$ and $S_{\text{ext}}$ are the core and extended flux densities respectively and the core is defined to be coincident with the AGN central engine. The core dominances are calculated for the rBZBs only.

To calculate $\rho_{144}$ for each rBZB, first the core and extended flux components were separated. We simulated a point source that, when convolved with the LDR2 PSF, had the same peak flux as the core in the LDR2 image. The integration of this Gaussian component represents the core flux density, $S_{\text{core}}$. We then computed the extended flux density as $S_{\text{ext}} = S_{\text{total}} - S_{\text{core}}$. Finally, $S_{\text{core}}$ and $S_{\text{ext}}$ were used with Eq. 2 to yield $\rho_{144}$.

3.3. *Spectral properties at 144 MHz*

As mentioned in § 2.2, NVSS fluxes were used with LDR2 to compute the 144–1400 MHz spatially integrated spectral indices of the BZBs. These spectral indices were then used to calculate the total 144 MHz specific luminosities of the sources in units of W Hz⁻¹.

4. RESULTS

4.1. *Overview*

The median properties of the sample are provided in Table 2, with uncertainties derived from bootstrapping. Values for each BZB are given in Table 5, with the full

<table><caption>**Table 2.** Median properties of the resolved and unresolved BL Lacs.</caption>
  <thead>
    <tr>
      <th></th>
      <th>Resolved</th>
      <th>Unresolved</th>
      <th>Total</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <td>N</td>
      <td>68</td>
      <td>31</td>
      <td>99</td>
    </tr>
    <tr>
      <td>z</td>
      <td>0.38 ± 0.02</td>
      <td>0.37 ± 0.04</td>
      <td>0.38 ± 0.02</td>
    </tr>
    <tr>
      <td>log<sub>10</sub> L<sub>144</sub></td>
      <td>25.6 ± 0.1</td>
      <td>24.8 ± 0.2</td>
      <td>25.3 ± 0.2</td>
    </tr>
    <tr>
      <td>α<sup>1400</sup><sub>144</sub></td>
      <td>-0.37 ± 0.06</td>
      <td>-0.12 ± 0.05</td>
      <td>-0.30 ± 0.03</td>
    </tr>
    <tr>
      <td>ρ<sub>144</sub></td>
      <td>0.74 ± 0.06</td>
      <td>...</td>
      <td>...</td>
    </tr>
    <tr>
      <td>D (kpc)</td>
      <td>69 ± 4</td>
      <td>...</td>
      <td>...</td>
    </tr>
    <tr>
      <td>Φ (")</td>
      <td>13 ± 2</td>
      <td>...</td>
      <td>...</td>
    </tr>
  </tbody>
</table>

NOTE—Number of sources is *N*, redshift is *z*, 144 MHz luminosity is *L*<sub>144</sub>, 144–1400 MHz spectral index is *α*<sup>1400</sup><sub>144</sub>, 144 MHz core dominance is *ρ*<sub>144</sub>, spatial extent is *D*, and angular extent is *Φ*.

Table accessible with the supplementary material on-line<fnref n="5" />.

4.2. Visual inspection of the BZB images

LDR2, FIRST and Pan-STARRS1 *i*-band images of all 99 BZBs are included with the supplementary material online². Broadly, we identify three classes of morphology for the BZBs. However, several sources have highly com- plex, asymmetrical morphologies, and in-depth studies using spatially-resolved spectral information would be necessary to fully interpret the source structures. While imaging artefacts cannot be categorically ruled out for all sources, they do not appear to be a significant factor in any case.

Regarding morphology, firstly there are the 33 BZBs that do not appear to be extended in LDR2; these sources are the uBZBs by definition.

The second group of sources are the rBZBs where the extended morphology seems to consist exclusively of dif- fuse emission; there are 52/99 such BZBs. This un- beamed emission is typically of low surface brightness and is believed to be emission from the extended lobes, where the jet interacts with the intergalactic medium. These BZBs typically appear as point sources in FIRST because the diffuse emission is not detected. When the inclination angle between the jet and our line-of-sight is small, the diffuse emission surrounds the beamed core component (e.g. J0911+3349 in Fig. 2). Alternatively, the inclination angle could be nonzero (but small) or there could be some level of jet bending, in which case there can be a distinct region of diffuse emission separate

<footnote n="5">⁵ Accessible at http://tiny.cc/bzb-images-mooney-2020 for review.</footnote>
from the core. For example, for J1000+5746, termination regions associated with the jet and counterjet are potentially observed (Fig. 2).

The remaining 14 BZBs are also extended, but some fraction of the extended emission is likely to be beamed. These sources are either extended in FIRST or have compact regions of high surface brightness in LDR2, with fluxes comparable to the core (e.g. J1340+4410 in Fig 2). This beamed emission is believed to be associated with relativistic bulk motion in the jet, in contrast to the unbeamed emission that is linked with the terminus of the jet. These sources have one-sided jets in FIRST but emission associated with a counterjet is detected in LDR2 in some cases.

### 4.3. Angular and spatial extents at 144 MHz

The angular extent distribution is shown in Fig. 3 (bottom) and redshift as a function of angular extent is shown in Fig. 3 (top). The BZBs with the most extended emission tend to be at lower redshifts, but generally there is no trend. The median angular extent of the rBZBs is $(13 \pm 2)''$. Note that we classify the sources as rBZBs or uBZBs only to aid the analysis. That is, the uBZBs and rBZBs likely form a continuous distribution rather than representing distinct BZB subclasses.

Fig. 4 shows the distribution of the spatial extents for the rBZBs, where the median rBZB spatial extent is $(69 \pm 4)$ kpc. There is a large spread in the extents, ranging up to 410 kpc, and two sources, J1231+3711 ($z = 0.219$), and J1340+4410 ($z = 0.546$) are more than 300 kpc. While this implies large deprojected (although not necessarily unphysical) sizes, it is possible that these jets are bent by some degree (as is typical for astrophysical jets), which could have the effect of reducing the deprojected sizes.

### 4.4. Core dominances at 144 MHz

We find 19/66 of the rBZBs are core-dominated (i.e. $\rho_{144} \ge 1$). The median core dominance of the rBZBs is $0.74 \pm 0.06$. Fig. 5 (bottom) shows the distribution of $\rho_{144}$ for the rBZBs, which spans two orders of magnitude, reflecting the sensitivity of the beaming factor with respect to the inclination angle and the Lorentz factor. There is a trend$^6$ ($r = 0.63, N = 64$, and $p = 3 \times 10^{-8}$) between $\rho_{144}$ and $\alpha_{144}^{1400}$ (Fig. 5; top), where the least core-dominated BZBs tend to have steeper $\alpha_{144}^{1400}$.

### 4.5. Spectral properties at 144 MHz

Of the 99 sources, 22 BZBs have steep spectra ($\alpha_{144}^{1400} < -0.5$), 74 are flat ($-0.5 \le \alpha_{144}^{1400} \le 0.5$), and 3 are inverted ($\alpha_{144}^{1400} > 0.5$). The median 144–1400 MHz spectral index for the sample is also flat ($-0.30 \pm 0.03$). A permutation test (Good 2013) was used to show that the median spectral index of the well-resolved sources ($\alpha_{144}^{1400} = -0.37 \pm 0.06$) is steeper than that of the unresolved sources ($\alpha_{144}^{1400} = -0.12 \pm 0.05$) at a significance level of $p < 3 \times 10^{-4}$. The median spectral indices can be classed as flat for both subsamples. Fig. 6 shows the distribution of $\alpha_{144}^{1400}$ for the uBZBs (top) and rBZBs (bottom).

### 4.6. Radio luminosities and redshifts

The integrated specific luminosities (W Hz$^{-1}$) of the sources were taking into account using $\alpha_{144}^{1400}$. Fig. 7 (top) shows the luminosity at 144 MHz, $L_{144}$, as a function of redshift, $z$. While $L_{144}$ loosely increases with $z$, this could be attributed to Malmquist bias (Malmquist 1925). The redshift distributions for the uBZBs and rBZBs are shown in Fig. 7 (middle; bottom). All BZBs in the sample are at $z \le 0.761$, and the median redshifts of the uBZBs and rBZBs are $0.38 \pm 0.02$ and $0.37 \pm 0.04$ respectively.

The distribution of luminosities for the uBZBs and the rBZBs is shown in Fig. 8 (top two panels respectively). There is a clear distinction between these distributions, where rBZBs tend to be more luminous. The rBZBs have higher total luminosities due to the presence of more extended radio emission. The luminosities of the core and extended components were then calculated for the rBZBs, where the core and extended spectral indices were assumed to be 0 and $-0.8$ respectively. This assumption is relatively uncontroversial because the cores of blazars are known to have flat spectral indices, whereas optically-thin synchrotron emission has a spectral index close to $-0.8$. The distributions for the core and extended luminosities are shown in Fig. 8 (bottom two panels respectively). While the rBZBs tend to be more luminous than the uBZBs, the distributions of the core luminosities of the rBZBs and the total luminosities of the uBZBs are similar. The total and core luminosities are shown as a function of redshift in Table 3.

### 4.7. The radio-to-$\gamma$-ray connection

In total, 53/99 BZBs had $\gamma$-ray detections; properties of the $\gamma$-ray BZBs are given in Table 4. The medians of the logarithms of the integrated photon fluxes at 1–100 GeV for the rBZBs ($-9.61 \pm 0.08$) and uBZBs ($-9.61 \pm 0.13$) are the same within uncertainty. The $\gamma$-ray-detected and non-$\gamma$-ray-detected BZBs are similar

<footnote>$^6$ The Pearson correlation coefficient, $r$, is reported with the associated probability value, $p$.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Multiwavelength images of J0911+3349 (i–iii), J1000+5746 (iv–vi), and J1340+4410 (vii–ix). The complete figure set (281 images) is available in the online journal. (i, iv, vii) LDR2 images with LDR2 contours. Contours mark 5σ, 10σ, 20σ, and 40σ levels at 144 MHz, where σ is the local RMS noise. Dashed grey contours mark the −3σ level. Colour bars range from 0 mJy beam⁻¹ to the peak flux in each image. The dashed blue circle is the FWHM of the LDR2 PSF. (ii, v, viii) FIRST maps with LDR2 contours. (iii, vi, ix) Pan-STARRS i-band images with LDR2 contours.</figcaption></figure>

in terms of spatial extents, angular extents, luminosities and core dominances. However, the γ-ray BZBs tend to have smaller redshifts. The median radio spectral index of the γ-ray BZBs is also slightly flatter than that of the BZBs without a detection.

Fig. 9 shows the total (top), core (middle), and extended (bottom) 144 MHz flux density against $S_{\text{GeV}}$ for the rBZBs. The Pearson correlation coefficient between $S_{\text{total}}$ and $S_{\text{GeV}}$ is 0.52 ($N = 36$ and $p = 1 \times 10^{-3}$). For $S_{\text{core}}$ and $S_{\text{GeV}}$, $r$ increases to 0.69 ($N = 36$ and $p = 3 \times 10^{-6}$), while for $S_{\text{ext}}$ and $S_{\text{GeV}}$, $r = 0.42$ ($N = 35$ and $p = 1 \times 10^{-2}$).

5. DISCUSSION

Previous studies identified extended emission at GHz frequencies around BZBs (e.g. Ulvestad et al. 1983), but the preliminary LoTSS Second Data Release has made it possible to systematically study the morphology of BZBs at 144 MHz for the first time, where we expect to find a larger fraction of diffuse emission due to the steeper spectra of the extended structures. We identified extended emission around 66/99 BZBs at 144 MHz. We compare BZBs and LERGs because radio-loud AGN unification theories predict that these
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Redshift, *z*, versus intrinsic angular extent, $\Phi$, (top) and the distribution of angular extents (bottom) of the resolved and unresolved BL Lacs.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Distribution of spatial extents, *D*, of the resolved BL Lacs.</figcaption></figure>

observational source classifications are the same source type intrinsically (Massaro et al. 2020a,b; Capetti et al. 2020). The spatial extents spanned a wide range (up to 410 kpc), and this is in line with predictions of the AGN unification scenario, where BZBs are LERGs viewed end-on, because LERGs range from several kpc to several Mpc in the plane of the sky. For the BZBs that were not found to be extended, it is possible that the resolution is insufficient to spatially resolve these sources, or that extended low-surface brightness radio emission exists below the LDR2 sensitivity. Particularly for the BZBs with $S_{\text{total}} \lesssim 40$ mJy, the extended emission may be below the detection threshold of LDR2.

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Spectral index, $\alpha_{144}^{1400}$, versus core dominance, $\rho_{144}$, (top) and the distribution of core dominances (bottom) for the resolved BL Lacs.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Distribution of 144–1400 MHz spectral indices, $\alpha_{144}^{1400}$, of the unresolved (top) and resolved (bottom) BL Lacs.</figcaption></figure>

The core luminosities per redshift bin of the rBZBs are comparable to the total luminosities of the uBZBs, suggesting no major difference in intrinsic core power between the two sources. That is, at a given core power we find large ranges of *D* and $L_{\text{ext}}$. A similar effect has also been seen in Faranoff-Riley class 0 sources (FR-0s; Baldi et al. 2015) and FR-Is, where Capetti et al. (2020) found that for a given core power, FR-0s tend to have fainter diffuse emission than FR-Is. The key vari-
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Total luminosity at 144 MHz, $L_{144}$, versus redshift, $z$, (top) and the distribution of redshifts for the unresolved (middle) and resolved (bottom) BL Lacs.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Distribution of (from top to bottom) the 144 MHz luminosities for the unresolved BL Lacs, the luminosities for the resolved BL Lacs, the core luminosities for the resolved BL Lacs, and the extended luminosities for the resolved BL Lacs.</figcaption></figure>

<table><caption>**Table 3.** Median total and core luminosities of resolved and unresolved BL Lacs.</caption><thead><tr><th rowspan="2">z</th><th colspan="3">Resolved</th><th colspan="2">Unresolved</th></tr><tr><th>N</th><th>L<sub>total</sub></th><th>L<sub>core</sub></th><th>N</th><th>L<sub>total</sub></th></tr></thead><tbody><tr><td>0.0–0.2</td><td>8</td><td>3.7 × 10<sup>24</sup></td><td>1.4 × 10<sup>24</sup></td><td>5</td><td>2.0 × 10<sup>24</sup></td></tr><tr><td>0.2–0.4</td><td>29</td><td>3.2 × 10<sup>25</sup></td><td>9.2 × 10<sup>24</sup></td><td>14</td><td>4.9 × 10<sup>24</sup></td></tr><tr><td>0.4–0.6</td><td>18</td><td>8.4 × 10<sup>25</sup></td><td>1.9 × 10<sup>25</sup></td><td>12</td><td>1.0 × 10<sup>25</sup></td></tr><tr><td>0.6–0.8</td><td>11</td><td>5.5 × 10<sup>25</sup></td><td>2.3 × 10<sup>25</sup></td><td>2</td><td>1.5 × 10<sup>26</sup></td></tr><tr><td>All</td><td>66</td><td>3.5 × 10<sup>25</sup></td><td>1.2 × 10<sup>25</sup></td><td>33</td><td>6.2 × 10<sup>24</sup></td></tr></tbody></table>

NOTE—Redshift bin is $z$, number of sources is $N$, total luminosity is $L_{total}$, and core luminosity is $L_{core}$ (both in W Hz$^{-1}$).

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Total flux density, $S_{total}$, (top), core flux density, $S_{core}$ (middle), and extended flux density, $S_{ext}$ (bottom) at 144 MHz versus γ-ray photon flux, $S_{GeV}$, for the γ-ray-detected resolved BL Lacs. The correlation coefficient and statistical significance are $r$ and $p$ respectively.</figcaption></figure>

ables controlling whether a BZB is extended or not is unclear, and could relate to the initial conditions or the environment. It is possible that the uBZBs have either only switched on recently, where the AGN activity is too recent to have produced plumes, or perhaps they have been switched off for $\gtrsim 10^8$ years, as suggested by Pun- sly et al. (2015), where this timescale is an upper limit based on the fact that the synchrotron break frequency
<table><caption>**Table 4.** Median properties of the γ-ray and non-γ-ray detected BL Lacs. Reported uncertainties are the standard error.</caption><thead><tr><th></th><th>γ-ray</th><th>Non-γ-ray</th></tr></thead><tbody><tr><td>N</td><td>53</td><td>46</td></tr><tr><td>z</td><td>0.29 ± 0.04</td><td>0.46 ± 0.02</td></tr><tr><td>log<sub>10</sub> L<sub>144</sub></td><td>25.3 ± 0.3</td><td>25.3 ± 0.1</td></tr><tr><td>α<sub>144</sub><sup>1400</sup></td><td>-0.26 ± 0.05</td><td>-0.38 ± 0.06</td></tr><tr><td>log<sub>10</sub> S<sub>GeV</sub></td><td>-9.6 ± 0.1</td><td>...</td></tr><tr><td>N<sub>ext</sub></td><td>36</td><td>30</td></tr><tr><td>ρ<sub>144</sub></td><td>0.83 ± 0.17</td><td>0.72 ± 0.05</td></tr><tr><td>D (kpc)</td><td>70 ± 15</td><td>69 ± 6</td></tr><tr><td>Φ (′′)</td><td>15 ± 3</td><td>12 ± 1</td></tr></tbody></table>

NOTE—Number of sources is *N*, redshift is *z*, 144 MHz luminosity is $L_{144}$, 144–1400 MHz spectral index is $\alpha_{144}^{1400}$, γ-ray flux is $S_{GeV}$, number of extended BL Lacs is $N_{ext}$, 144 MHz core dominance of extended BL Lacs is $\rho_{144}$, spatial extent of extended BL Lacs is $D$, and angular extent of extended BL Lacs is $\Phi$.

drops to $\lesssim 100$ MHz, assuming $\alpha_{144}^{1400} \approx -0.8$, $a \sim 1 \mu G$ magnetic field, and no new injections of energetic particles (Cordey 1986; Jamrozy et al. 2007).

We directly measured the core dominances of the extended BZBs by separating the core and extended flux components. Our results are broadly in line with Giroletti et al. (2016) and Fan & Wu (2018), where the core dominances of blazars were calculated by decomposing $S_{total}$ using estimates of the core and extended spectral indices. While the core component contributes more than half of the total flux for 19/66 of the rBZBs, the majority of these extended BZBs were not core-dominated. This finding is in agreement with d’Antonio et al. (2019) where they state that deeper low frequency surveys should result in lower core dominances because more radio lobe emission can be detected. Thus, the core dominance may not be a reliable parameter for selecting new blazar candidates in future low frequency surveys.

The median radio spectral index for the sample is flat, which confirms several recent studies that have found that BZBs, and blazars generally, have flat spectral indices down to $\sim 100$ MHz at least (Giroletti et al. 2016; Mooney et al. 2019; d’Antonio et al. 2019). However, the most extended BZBs tended to have spectral indices closer to $-0.8$, similar to classic radio galaxies, likely due to the presence of diffuse emission.

The population of relativistic electrons that give rise to the beamed radio emission are thought to be responsible for producing the γ-ray emission, by upscattering seed photons from external radiation fields to γ-ray energies. This is believed to be also the mechanism behind the radio-to-γ-ray connection at GHz frequencies. Previous studies have shown that this correlation weakens with decreasing radio frequency (e.g. Giroletti et al. 2016; Mooney et al. 2019). By calculating the radio-to-γ-ray correlation using the core and extended flux densities at 144 MHz separately, we have shown that this connection weakens at low radio frequencies because there is an increase in diffuse emission. The diffuse MHz emission is not expected to correlate strongly with the γ-ray emission because the diffuse MHz emission is unbeamed, and the electron population producing the diffuse MHz emission is likely to be distinct from the electrons that give rise to the γ-rays, reflecting the time-integrated history of jet activity rather than the instantaneous view given by the core emission. In contrast, the beamed core emission at low frequencies still correlates reasonably well with the γ-ray photon flux.

## 6. CONCLUSIONS AND FUTURE PROSPECTS

We have presented the first morphological study of a sample of BZBs at 144 MHz. Our findings are as follows:

* Extended emission was revealed around 66/99 BZBs. The distributions of spatial extents and the extended component luminosities at 144 MHz are consistent with expectations based on the AGN unification paradigm, where BZBs are the aligned counterparts of LERGs. For a given core luminosity and redshift, a range of spatial extents were found.

* The median integrated 144–1400 MHz spectral index for the BZBs was $-0.30 \pm 0.03$, confirming that the low frequency spectra of BZBs are still dominated by the beamed core component, but this depends on the spatial extent, with the spectral index being steeper for more extended BZBs.

* Of the 66 rBZBs, 19 were core-dominated at 144 MHz, and the most core-dominated sources tended to have smaller spatial extents and flatter spectral indices.

* Positive correlations were identified between the 1–100 GeV flux and both the 144 MHz total ($r = 0.52$) and core ($r = 0.69$) flux densities for the rBZBs. The correlation between the γ-ray flux and the extended emission was weaker ($r = 0.42$). This suggests that the 144 MHz emission from the
core is more likely to be directly related to the
present nuclear activity than the large-scale ex-
tended emission.

Characterising the diffuse emission around BZBs can
help to explain some of the outstanding questions re-
garding blazars. For example, measurements of the dif-
fuse flux allow for estimates of the jet power via scal-
ing relations (Meyer et al. 2011). Calculating the jet
power can shed light on the relationship between the
jet power and the AGN accretion rate, while estimates
of the diffuse flux and the spatial extent can also be
used to characterise the environments. We focused on
BZBs using LoTSS, the most sensitive low frequency
survey in existence. In future, LoTSS in-band spec-
tral index maps will make it possible to distinguish be-
tween the beamed jet flux and the unbeamed diffuse flux

that contribute to the extended emission. Other forth-
coming surveys will help further constrain the 144 MHz
morphology of BZBs. The LOFAR LBA Sky Survey
(LoLSS; de Gasperin et al., in prep) will eventually
cover the northern hemisphere sky at 42–66 MHz at a
sensitivity of 1 mJy beam⁻¹ and 15" resolution. With
LoLSS, we will be able to identify possible ultra-steep
spectrum radio halos around the unresolved BZBs that
are undetectable in LDR2. The future sky survey with
the LOFAR international stations (Morabito et al. in
prep) will image the 122–168 MHz sky at 0.3" resolu-
tion, potentially spatially resolving some of the unre-
solved BZBs. In addition, optical spectra provided by
the WEAVE-LOFAR survey (Smith et al. 2016) will be
beneficial in identifying LDR2 counterparts to uniden-
tified γ-ray sources. In the long term, surveys with the
Square Kilometre Array (SKA) will also be key in study-
ing this population (Kharb et al. 2016).
<table><caption>**Table 5.** Details of the BZB sample.</caption><thead><tr><th>Name</th><th>RA</th><th>Dec</th><th>z</th><th>Resolved?</th><th>S<sub>total</sub></th><th>S<sub>GeV</sub></th><th>α<sub>144</sub><sup>1400</sup></th><th>ρ<sub>144</sub></th><th>D</th><th>Φ</th><th>L<sub>144</sub></th></tr><tr><td></td><td>°</td><td>°</td><td></td><td></td><td>mJy</td><td>ph cm<sup>-2</sup>s<sup>-1</sup></td><td></td><td></td><td>kpc</td><td>"</td><td>WHz<sup>-1</sup></td></tr><tr><td></td><td>(1)</td><td>(2)</td><td>(3)</td><td>(4)</td><td>(5)</td><td>(6)</td><td>(7)</td><td>(8)</td><td>(9)</td><td>(10)</td><td>(11)</td><td>(12)</td></tr></thead><tbody><tr><td>J0047+3948</td><td>11.9801</td><td>39.8160</td><td>0.252</td><td>T</td><td>264 ± 53</td><td>(7.4 ± 0.5) × 10<sup>-10</sup></td><td>-0.46 ± 0.09</td><td>0.36 ± 0.07</td><td>68.1 ± 0.4</td><td>17.1 ± 0.1</td><td>(4.6 ± 0.9) × 10<sup>25</sup></td></tr><tr><td>J0112+2244</td><td>18.0243</td><td>22.7441</td><td>0.265</td><td>T</td><td>233 ± 47</td><td>(8.0 ± 0.2) × 10<sup>-9</sup></td><td>0.22 ± 0.04</td><td>1.21 ± 0.25</td><td>38.3 ± 0.1</td><td>9.3 ± 0.1</td><td>(3.9 ± 0.8) × 10<sup>25</sup></td></tr><tr><td>J0123+3420</td><td>20.7860</td><td>34.3468</td><td>0.272</td><td>T</td><td>169 ± 34</td><td>(2.4 ± 0.4) × 10<sup>-10</sup></td><td>-0.60 ± 0.12</td><td>0.20 ± 0.04</td><td>95.9 ± 0.4</td><td>22.9 ± 0.1</td><td>(3.6 ± 0.7) × 10<sup>25</sup></td></tr><tr><td>J0203+3042</td><td>30.9345</td><td>30.7105</td><td>0.761</td><td>T</td><td>345 ± 69</td><td>(2.24 ± 0.09) × 10<sup>-9</sup></td><td>-0.30 ± 0.06</td><td>0.84 ± 0.17</td><td>130.7 ± 0.3</td><td>17.4 ± 0.1</td><td>(6.4 ± 1.3) × 10<sup>26</sup></td></tr><tr><td>J1808+3520†</td><td>272.2071</td><td>35.3452</td><td>0.141</td><td>T</td><td>94 ± 19</td><td>(2.2 ± 0.4) × 10<sup>-10</sup></td><td>-0.50 ± 0.10</td><td>0.19 ± 0.04</td><td>72.7 ± 0.8</td><td>29.1 ± 0.3</td><td>(4.7 ± 0.9) × 10<sup>24</sup></td></tr><tr><td>J2229+2255†</td><td>337.2966</td><td>22.9166</td><td>0.440</td><td>T</td><td>18 ± 4</td><td>...</td><td>0.35 ± 0.07</td><td>0.31 ± 0.06</td><td>58.3 ± 1.7</td><td>10.1 ± 0.3</td><td>(8.9 ± 1.8) × 10<sup>24</sup></td></tr><tr><td>J2237+1840†</td><td>339.2701</td><td>18.6822</td><td>0.722</td><td>F</td><td>7 ± 1</td><td>...</td><td>-0.08 ± 0.02</td><td>...</td><td>≤ 28.1 ± 0.7</td><td>≤ 3.8 ± 0.1</td><td>(9.9 ± 2.0) × 10<sup>24</sup></td></tr><tr><td>J2343+3439</td><td>355.8899</td><td>34.6641</td><td>0.366</td><td>T</td><td>203 ± 41</td><td>(3.9 ± 0.5) × 10<sup>-10</sup></td><td>-0.79 ± 0.16</td><td>0.11 ± 0.02</td><td>151.9 ± 1.7</td><td>29.6 ± 0.3</td><td>(8.9 ± 1.8) × 10<sup>25</sup></td></tr></tbody></table>

NOTE—(1) Source name. (2-3) Coordinates in the J2000 frame. (4) Redshift. (5) A flag stating whether the source is resolved ('T') or not ('F'). (6) Integrated flux density at 144 MHz from LDR2. (7) γ-ray flux from 4LAC. (8) 144-1400 MHz spectral index with LDR2 and NVSS. (9) Core dominance. (10) Spatial extent. (11) Angular extent. (12) Radio luminosity.

This table is available in its entirety in machine-readable form.

† BZBs not in Roma-BZCAT v5.0.
ACKNOWLEDGMENTS

SM acknowledges support from the Irish Research Council Postgraduate Scholarship and UCD U21 funding. BM acknowledges support from the UK Science and Technology Facilities Council (STFC) under grants ST/R00109X/1 and ST/R000794/1. MJH acknowledges support from STFC [ST/R000905/1]. FM acknowledges financial contribution from the agreement ASI-INAF n.2017-14-H.0. This work is supported by the “Departments of Excellence 2018–2022” grant awarded by the Italian Ministry of Education, University and Research (MIUR) (L. 232/2016). This research has made use of resources provided by the Compagnia di San Paolo for the grant awarded on the BLENV project (S1618_L1_MASF_01) and by the Ministry of Education, Universities and Research for the grant MASF FFABR 17.01. This investigation is supported by the National Aeronautics and Space Administration (NASA) grants GO4-15096X, AR6-17012X, GO6-17081X, GO9-20083X, and GO0-21110X.

LOFAR 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 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 LO-
FAR 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 Insti-
tute for Computing (NIC).

This research made use of the University of Hert-
fordshire high-performance computing facility and the
LOFAR-UK computing facility located at the Uni-
versity of Hertfordshire and supported by STFC
[ST/P000096/1], and of the Italian LOFAR IT com-
puting infrastructure supported and operated by INAF,
and by the Physics Department of Turin university (un-
der an agreement with Consorzio Interuniversitario per
la Fisica Spaziale) at the C3S Supercomputing Centre,
Italy.

The Pan-STARRS1 Surveys (PS1) and the PS1 public science archive have been made possible through contributions by the Institute for Astronomy, the University of Hawaii, the Pan-STARRS Project Office, the Max-Planck Society and its participating institutes, the Max Planck Institute for Astronomy, Heidelberg and the Max Planck Institute for Extraterrestrial Physics, Garching, The Johns Hopkins University, Durham University, the University of Edinburgh, the Queen's University Belfast, the Harvard-Smithsonian Center for Astrophysics, the Las Cumbres Observatory Global Telescope Network Incorporated, the National Central University of Taiwan, the Space Telescope Science Institute, the National Aeronautics and Space Administration under Grant No. NNX08AR22G issued through the Planetary Science Division of the NASA Science Mission Directorate, the National Science Foundation Grant No. AST-1238877, the University of Maryland, Eotvos Lorand University (ELTE), the Los Alamos National Laboratory, and the Gordon and Betty Moore Foundation.

The National Radio Astronomy Observatory is a facil-
ity of the National Science Foundation operated under
cooperative agreement by Associated Universities, Inc.

*Facilities:* LOFAR, Haleakala Observatory (Pan-STARRS), and VLA.

*Software:* DDFacet (Tasse et al. 2018), KillMS (Tasse 2014; Smirnov & Tasse 2015), Python (van Rossum 1995), NDPPP (van Diepen et al. 2018), NumPy (van der Walt et al. 2011), Astropy (Astropy Collaboration et al. 2013), TOPCAT (Taylor 2013), PyBDSF (Mohan & Rafferty 2015), and Matplotlib (Hunter 2007).

REFERENCES

Abdollahi, S., Acero, F., Ackermann, M., et al. 2020, ApJS,

247, 33, doi: 10.3847/1538-4365/ab6bcb

Ackermann, M., Ajello, M., Allafort, A., et al. 2011, ApJ,

741, 30, doi: 10.1088/0004-637X/741/1/30
Ajello, M., Angioni, R., Axelsson, M., et al. 2020, ApJ, 892, 105, doi: 10.3847/1538-4357/ab791e

Álvarez Crespo, N., Massaro, F., Milisavljevic, D., et al. 2016, AJ, 151, 95, doi: 10.3847/0004-6256/151/4/95

Antonucci, R. R. J., & Ulvestad, J. S. 1985, ApJ, 294, 158, doi: 10.1086/163284

Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068

Atwood, W. B., Abdo, A. A., Ackermann, M., et al. 2009, ApJ, 697, 1071, doi: 10.1088/0004-637X/697/2/1071

Baldi, R. D., Capetti, A., & Giovannini, G. 2015, A&A, 576, A38, doi: 10.1051/0004-6361/201425426

Becker, R. H., White, R. L., & Helfand, D. J. 1995, ApJ, 450, 559, doi: 10.1086/176166

Best, P. N., & Heckman, T. M. 2012, MNRAS, 421, 1569, doi: 10.1111/j.1365-2966.2012.20414.x

Capetti, A., Massaro, F., & Baldi, R. D. 2020, A&A, 633, A161, doi: 10.1051/0004-6361/201935962

Chambers, K. C., Magnier, E. A., Metcalfe, N., et al. 2016, arXiv e-prints, arXiv:1612.05560. https://arxiv.org/abs/1612.05560

Condon, J. J., Cotton, W. D., Greisen, E. W., et al. 1998, AJ, 115, 1693, doi: 10.1086/300337

Cordey, R. A. 1986, MNRAS, 219, 575, doi: 10.1093/mnras/219.3.575

d'Antonio, D., Giroletti, M., Giovannini, G., & Maini, A. 2019, MNRAS, 490, 5798, doi: 10.1093/mnras/stz2967

de Gasperin, F., Intema, H. T., & Frail, D. A. 2018, MNRAS, 474, 5008, doi: 10.1093/mnras/stx3125

de Menezes, R., Peña-Herazo, H. A., Marchesini, E. J., et al. 2019, A&A, 630, A55, doi: 10.1051/0004-6361/201936195

Fan, X.-L., & Wu, Q. 2018, ApJ, 869, 133, doi: 10.3847/1538-4357/aaecee

Giovannini, G., Feretti, L., & Comoretto, G. 1990, ApJ, 358, 159, doi: 10.1086/168970

Giroletti, M., Massaro, F., D'Abrusco, R., et al. 2016, A&A, 588, A141, doi: 10.1051/0004-6361/201527817

Good, P. 2013, Permutation Tests: A Practical Guide to Resampling Methods for Testing Hypotheses, Springer Series in Statistics (Springer New York). https://books.google.ie/books?id=pK3hBwAAQBAJ

Hardcastle, M. J., & Krause, M. G. H. 2014, MNRAS, 443, 1482, doi: 10.1093/mnras/stu1229

Healey, S. E., Romani, R. W., Taylor, G. B., et al. 2007, ApJS, 171, 61, doi: 10.1086/513742

Helfand, D. J., White, R. L., & Becker, R. H. 2015, ApJ, 801, 26, doi: 10.1088/0004-637X/801/1/26

Hine, R. G., & Longair, M. S. 1979, MNRAS, 188, 111, doi: 10.1093/mnras/188.1.111

Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90, doi: 10.1109/MCSE.2007.55

Jamrozy, M., Konar, C., Saikia, D. J., et al. 2007, MNRAS, 378, 581, doi: 10.1111/j.1365-2966.2007.11782.x

Kharb, P., Lal, D. V., Singh, V., et al. 2016, Journal of Astrophysics and Astronomy, 37, 34, doi: 10.1007/s12036-016-9411-z

Landoni, M., Massaro, F., Paggi, A., et al. 2015, AJ, 149, 163, doi: 10.1088/0004-6256/149/5/163

Landt, H., Padovani, P., Perlman, E. S., et al. 2001, MNRAS, 323, 757, doi: 10.1046/j.1365-8711.2001.04269.x

Laurent-Muehleisen, S. A., Kollgaard, R. I., Moellenbrock, G. A., & Feigelson, E. D. 1993, AJ, 106, 875, doi: 10.1086/116691

Malmquist, K. G. 1925, Meddelanden fran Lunds Astronomiska Observatorium Serie I, 106, 1

Massaglia, S., Bodo, G., Rossi, P., Capetti, S., & Mignone, A. 2016, A&A, 596, A12, doi: 10.1051/0004-6361/201629375

Massaro, E., Maselli, A., Leto, C., et al. 2015, Ap&SS, 357, 75, doi: 10.1007/s10509-015-2254-2

Massaro, F., Capetti, A., Paggi, A., et al. 2020a, ApJL, 900, L34, doi: 10.3847/2041-8213/abac56

Massaro, F., D'Abrusco, R., Giroletti, M., et al. 2013a, ApJS, 207, 4, doi: 10.1088/0067-0049/207/1/4

Massaro, F., Giroletti, M., D'Abrusco, R., et al. 2014, ApJS, 213, 3, doi: 10.1088/0067-0049/213/1/3

Massaro, F., Giroletti, M., Paggi, A., et al. 2013b, ApJS, 208, 15, doi: 10.1088/0067-0049/208/2/15

Massaro, F., Álvarez Crespo, N., D'Abrusco, R., et al. 2016, Ap&SS, 361, 337, doi: 10.1007/s10509-016-2926-6

Massaro, F., Capetti, A., Paggi, A., et al. 2020b, ApJS, 247, 71, doi: 10.3847/1538-4365/ab799e

Meyer, E. T., Fossati, G., Georganopoulos, M., & Lister, M. L. 2011, ApJ, 740, 98, doi: 10.1088/0004-637X/740/2/98

Mohan, N., & Rafferty, D. 2015, PyBDSF: Python Blob Detection and Source Finder. http://ascl.net/1502.007

Mooney, S., Quinn, J., Callingham, J. R., et al. 2019, A&A, 622, A14, doi: 10.1051/0004-6361/201833937

Morganti, R., Oosterloo, T. A., Reynolds, J. E., Tadhunter, C. N., & Migenes, V. 1997, MNRAS, 284, 541, doi: 10.1093/mnras/284.3.541

Paiano, S., Falomo, R., Treves, A., Franceschini, A., & Scarpa, R. 2019, ApJ, 871, 162, doi: 10.3847/1538-4357/aaf6e4

Peña-Herazo, H. A., Massaro, F., Chavushyan, V., et al. 2019, Ap&SS, 364, 85, doi: 10.1007/s10509-019-3574-4
Perlman, E. S., & Stocke, J. T. 1993, ApJ, 406, 430,
doi: 10.1086/172454

Punsly, B., Marziani, P., Kharb, P., O'Dea, C. P., &
Vestergaard, M. 2015, ApJ, 812, 79,
doi: 10.1088/0004-637X/812/1/79

Ricci, F., Massaro, F., Landoni, M., et al. 2015, AJ, 149, 160, doi: 10.1088/0004-6256/149/5/160

Shimwell, T. W., Röttgering, H. J. A., Best, P. N., et al.
2017, A&A, 598, A104,
doi: 10.1051/0004-6361/201629313

Shimwell, T. W., Tasse, C., Hardcastle, M. J., et al. 2019, A&A, 622, A1, doi: 10.1051/0004-6361/201833559

Smirnov, O. M., & Tasse, C. 2015, MNRAS, 449, 2668,
doi: 10.1093/mnras/stv418

Smith, D. J. B., Best, P. N., Duncan, K. J., et al. 2016, in
SF2A-2016: Proceedings of the Annual meeting of the
French Society of Astronomy and Astrophysics, 271-280.
https://arxiv.org/abs/1611.02706

Tasse, C. 2014, A&A, 566, A127,
doi: 10.1051/0004-6361/201423503

Tasse, C., Hugo, B., Mirmont, M., et al. 2018, A&A, 611, A87, doi: 10.1051/0004-6361/201731474

Tasse, C., Shimwell, T., Hardcastle, M. J., et al. 2020, arXiv e-prints, arXiv:2011.08328.
https://arxiv.org/abs/2011.08328

Taylor, M. 2013, Starlink User Note, 253

The LIGO Scientific Collaboration, The VIRGO Collaboration, The 1M2H Collaboration, et al. 2017, Nature, 551, 85, doi: 10.1038/nature24471

Ulvestad, J. S., Johnston, K. J., & Weiler, K. W. 1983, ApJ, 266, 18, doi: 10.1086/160755

Urry, C. M., & Padovani, P. 1995, PASP, 107, 803,
doi: 10.1086/133630

van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science and Engineering, 13, 22, doi: 10.1109/MCSE.2011.37

van Diepen, G., Dijkema, T. J., & Offringa, A. 2018, DPPP: Default Pre-Processing Pipeline. http://ascl.net/1804.003

van Haarlem, M. P., Wise, M. W., Gunst, A. W., et al. 2013, A&A, 556, A2, doi: 10.1051/0004-6361/201220873

van Rossum, G. 1995, Python tutorial, Tech. Rep. CS-R9526, Centrum voor Wiskunde en Informatica (CWI), Amsterdam

Wilkes, B. 2004, in Astronomical Society of the Pacific Conference Series, Vol. 311, AGN Physics with the Sloan Digital Sky Survey, ed. G. T. Richards & P. B. Hall, 37.
https://arxiv.org/abs/astro-ph/0310905

Wright, E. L. 2006, PASP, 118, 1711, doi: 10.1086/510102