© ESO 2022

1 Introduction

3C295 is one of the most famous jetted active galactic nuclei (AGN; Padovani 2017) in the sky, which was first detected in the 3C survey (Edge et al. 1959). It is associated with a nearby (z = 0.461, Ahn et al. 2013) massive cD galaxy (Mathieu & Spinrad 1981). Its electromagnetic spectrum has been well studied from the radio (Akujor et al. 1990; Perley & Taylor 1991; Napier et al. 1983; Perley et al. 2011), to optical (e.g. Thimm et al. 1994) to X-ray (e.g. Harris et al. 2000; Brunetti et al. 2001) frequencies. The radio integrated properties are so well known that 3C295 is a standard flux density calibrator from megahertz to gigahertz frequencies (Scaife & Heald 2012; Perley & Butler 2017), due to its extremely high flux density (90.87 Jy at 144 MHz, 19.42 Jy at 1.5 GHz) and compactness, with a largest angular size of 6′′. Although high frequency radio observations have resolved the sub-structure of 3C295 (e.g. Perley & Taylor 1991), the spatially resolved morphology at ~100 MHz frequencies remains unstudied.

3C295 is a scientifically interesting source for understanding the emission and absorption processes in radio-loud AGN. To understand these processes, we require spatially resolved information across a broad range of radio frequencies. While optically thin radio emission can often be adequately described using a simple power law model, sources such as 3C295 can be affected by processes such as spectral ageing at higher frequencies (≳500 MHz, depending on age) or absorption processes e.g. free-free absorption or synchrotron self-absorption at low frequencies (≲200 MHz). These emission and absorption processes can be modelled, but doing so requires measurements at multiple frequencies with enough spatial resolution to model distinct components.While the integrated radio spectrum of 3C295 is very well-constrained and shows a turn-over at ~80 MHz (Scaife & Heald 2012), it is impossible to tell without spatially resolved information at these frequencies which components of the source are driving this.

To answer this question, we use the International LOw Frequency ARray (LOFAR) Telescope (ILT, van Haarlem et al. 2013) to make sub-arcsecond resolution images of 3C295 at 132 MHz for the first time. The geographical spread of the ILT stations across Europe provide baselines of up to ~2000 km, and recent advances in calibration techniques have made it easier than ever to use the array to achieve tenths of arcsecond resolution at 150 MHz (Morabito et al. 2022). When combined with matched-resolution archival Very Large Array(VLA) and Multi-Element Radio-Linked Interferometer (MERLIN) observations, they allowed us to perform the first spatially resolved spectral modelling of 3C295 to investigate the physical conditions of its radio components.

In Sect. 2, we give an overview of the data and the resulting images. In Sect. 3 we analyse the spectral properties of these images. Finally, in Sect. 4 we close out with a discussion. Throughout this paper, we assume a standard ΛCDM (Cold Dark Matter) cosmologicalmodel, using Wright (2006) to calculate distances at 3C295’s redshift (Table 1). Throughout this paper, we define the spectral index α by S_ν ∝ ν^α, where S_ν is the flux density per unit frequency ν.

2 Observations and results

2.1 Data reduction

Each telescope has its own data calibration methods, and we describe our LOFAR, MERLIN, and VLA data reduction separately. Our target is one of the radio flux density calibrators described in Perley & Butler (2017). Although we sometimes used the flux scale described in Scaife & Heald (2012) during processing, which we note when relevant in the following section, this flux scale diverges at higher frequencies, and so we aligned all images to the Perley & Butler (2017) flux scale in the end. This results in an absolute flux scaling error of 5% for all frequencies.

2.1.1 LOFAR

The LOFAR HBA data were taken during a targeted observation of 3C 295 (PI: Sweijen; project code COM_010, L693723) on January 25, 2019. It was observed from 01:41:00 to 09:40:59 UT, for a total of 8 h on source and was book-ended by two 15 min calibrator observations. A typical observing setup was used where the data were recorded at 3.05 kHz frequency resolution and 1 s time resolution. The data were then archived to the Long Term Archive (LTA) after they were averaged in frequency to a resolution of 12.21 kHz and RFI was excised using AOFlagger Offringa (2010). The frequency coverage ranged from 120 to 168 MHz for a total of 48 MHz of bandwidth. Fifty one stations in total (24 core stations (CS), 14 remote stations (RS) and 13 international stations) participated in the observation, all of which provided usable data except for the remote station RS508.

First the prefactor calibrator pipeline was run on one of the calibrator observations to obtain corrections for polarisation alignment, find the station bandpasses and correct for clock offsets (de Gasperin et al. 2019). After this all core stations were combined into a single super-station, ST001 to significantly reduce the field of view. Finally, the data were averaged to a time resolution of 4 s and a frequency resolution of 48.82 kHz. The data were then self-calibrated. We used PyBDSF (Mohan & Rafferty 2015) to generate a starting model from a VLA A-configuration image at 8.561 GHz shown in Fig. A.1 (Gilbert et al. 2004). Self-calibration then followed a strategy similar to that presented in van Weeren et al. (2021). First 5 cycles of phase-only calibration were performed at the resolution of the data. Subsequently the calibration switched to 5 cycles of phase-only calibration followed by amplitude calibration and a frequency interval of 195.28 kHz and a time interval of 15 min. Such short intervals were possible due to the exceptionally high S/N present on this class of source. Calibration was restricted to baselines 40 kλ or longer to help stabilise calibration. After each iteration a new image was made and > 7 σ_rms emission was used to update the model for the next iteration, where σ_rms is the local root-mean-square noise in the image. Imaging was done using multi-frequency synthesis, a robust − 1 weighting and multi-scale clean. WSClean (Offringa et al. 2014) was used for imaging and DP3 (van Diepen et al. 2018) was used for calibration. Finally, we note that international LOFAR data can be impacted by large dispersive delays, which are challenging to remove for low signal to noise sources where large bandwidths are necessary to find good solutions. Fortunately, 3C295 has a high enough flux density that we can perform self-calibration per channel, removing the need to solve for these dispersive delays separately. This data were further self-calibrated using the Wirtinger pack software (Smirnov & Tasse 2015; Tasse et al. 2018), which consists of the DDFacet imager and killMS calibration software.

We used small intervals in time and frequency (8 s and 4 channels, respectively) so as to maximise our tracking of underlying changes in the antenna gains, solving for full-Jones matrices in each interval. Intervals of increased ionospheric activity (corresponding to noisier measurements and thus worse S/N for our gain solutions) were corrected for by using the quality-based weighting scheme developed in Bonnassieux et al. (2018). This was key in maximising the relative contribution of the international baselines to the final image. Finally, drifts in the gain were corrected for by fixing the integrated flux density of the model to the Scaife & Heald (2012) flux density scale before calibration at every self-calibration pass.

The final LOFAR image was made with an inner uv-cut of 20 kλ (shortest baseline in the data set was ~950λ), a Briggs weighting with robust parameter of − 2 (Briggs 1995), a uv-taper with FWHM of 0.6′′, and a pixel size of 0.1′′. This ensures that the effective uv-coverage of all observations can be matched as much as possible. This image was made using wsclean2.10.0’s (Offringa et al. 2014) multi-scale CLEAN algorithm (Offringa & Smirnov 2017), as DDFacet lacked an implementation of uv-tapering at the time of reduction. We note that all other images used exactly the same parameters.

The final image sensitivity is 0.6 mJy beam⁻¹ at 132 MHz (which corresponds to a dynamic range – DR – of 13 000). We note that our noise levels are nevertheless still determined by DR limitations.

2.1.2 MERLIN

The archival MERLIN observation was taken on 26 May 1998, with a central observing frequency of 1.658 GHz and a total bandwidth of 26 MHz. The observation included 6 antennas. The data were extracted from the MERLIN archive, and flux-calibrated with 3C286. We then performed 8 passes of self-calibration using killMS and DDFacet to improve the calibration, followed by a final re-imaging using wsclean2.1.0, using exactly the parameters and initial model as for our LOFAR data.

The shortest baseline available in our MERLIN observation is 25 k λ, which corresponds to an angular scale of 8.25′′. It therefore does not lead to loss of flux, as 3C295 has a LAS of 6.1′′. Using the same imaging procedure as for LOFAR ensures homogeneous uv-coverage for the angular scales of interest.

The final noise level achieved was 25 mJy beam⁻¹, leading to a dynamic range of about 200 in the final science image. The hotspot positions were used to align the final image with the LOFAR images. This correction of ~0.4′′ was done in the image plane.

2.1.3 VLA

Our VLA data were acquired from the NRAO Science Data Archive. We restricted our search to observations from A-configuration and longer than two hours in duration to maximise our uv-coverage. We used two observations from the NRAO Science Data Archive, each in A-configuration and longer than two hours in duration to maximise our uv-coverage and each comprising two frequency bands, taken as part of project AP135 (PI: R. Perley) on 12 September 1987. The central frequencies of each observation were 4.76 and 14.975 GHz, with bandwidths of 50 MHz. Both use the full 27 VLA antennas. This configuration leads to shortest baselines of 9.3 k λ at 5 GHz and 35 k λ at 15.5 GHz, corresponding to LASs of ~22.2′′ and ~3′′, respectively. Because the latter is smaller than the LAS of 3C295, we added a calibrated C-configurationd dataset from the VLA science data archive in order to properly sample 3C295 at 15 GHz.

The archival data sets contained uncalibrated visibilities. We therefore used the Wirtinger pack for our self-calibration, with the exact same procedure and initial model as for LOFAR. We then imaged our final visibilities using wsclean to match the imaging parameters used at all other frequencies. This converged quickly for both frequencies, allowing us to acquire science images with rms noise levels of 15 mJy beam⁻¹ and 670 μJy beam⁻¹ at 5 and 15 GHz respectively (corresponding to dynamic ranges of 165 and 754, respectively).The hotspot positions were used to align the final image with the LOFAR images, as for the final eMERLIN image. This offset was ~ 0.2′′.

2.2 Radio images

In Fig. 1a, we show our main result: the first image ever produced of 3C295 at 132 MHz with a resolution of 0.4′′. Images of the source at higher frequencies, produced using the archival data described previously, are also shown in Fig. 1. The properties of these images (name, frequency, noise, dynamic range) are given in Table 2. They are all created using similar uv-tapers and uv-cuts and thus have matching resolutions and uv-coverage. They are also spatially aligned on the hotspots. They have all been tied to the flux density scale described in Perley & Butler (2017), since the Scaife &Heald (2012) flux scale used during self-calibration diverges at higher frequencies.The flux scale error is simply taken to be the reference 5% at all frequencies, as 3C295 is one of the calibrator sources in Perley & Butler (2017).

We note that, because all these images were made using a Briggs weighting with robust parameter − 2, which corresponds to uniform weighting, their noise levels are relatively high. The lower colour map thresholds start at the 5σ level for each image, usingthe σ values listedin Table 2 except for the LOFAR images, which instead use a minimum threshold of 100 mJy beam⁻¹. This is due to structure in the residuals near the source. The image is overlaid onto itself as a contour, starting at the 3σ level.

2.3 Morphology of 3C295

3C295 is a compact source with a largest angular scale of 6.1′′, which is fully covered by our shortest baselines at all frequencies. Due to its very high flux density (90.87 Jy at 144 MHz, 19.42 Jy at 1.5 GHz), its features are recovered at high S/N.

The morphology of 3C295 as observed in our new LOFAR image at 132 MHz agrees well with the morphology observed at higher frequencies: see Fig. 1 (e.g. Perley & Taylor 1991) We resolve two separate sources of emission, one in the north-east and one in the south-west, each about 10 kpc wide.

The Northern component can be divided into sub-components: a dominant compact hot-spot, dominant in flux density, anda fainter lobe surrounding it, which becomes more prominent at lower frequencies. The Northern hotspot is most compact at higher frequencies, and the North lobe emission closest to the host galaxy (situated between theNorthern and Southern components) becomes more apparent as frequency decreases.

In the Southern component, the hotspot has been reported in the literature as elongated along an East-West axis (e.g. Gupta et al. 2017). At the resolutions we achieve, we find a high-flux-density, compact source with an extension of more diffuse emission to the North-East. These are referred to, respectively, as the Southern hotspot and Southern lobe.

These four morphological components are identified with drawn regions in Fig. 2, and are all larger than the restoring beam (0.6′′∕3.6 kpc). Although we use only the regions identified in Fig. 2 in this work, the source itself is complex at 132 MHz.

Fig. 1 High resolution images of 3C295 at different frequencies. Flux density units are Jy bm−1, with colour scales starting at 5σ and tuned to show the source morphology. Each image includes an overlay of itself, with 10 levels starting at 3σ and increasing exponentially until reaching the maximum value in the image. The restoring beam is 0.6′′ in all cases, and the general properties are summarised in Table 2. Panel a: LOFAR, 132 MHz. Panel b: MERLIN, 1.658 GHz. Panel c: VLA, 4.760 GHz. Panel d: VLA, 14.975 GHz.

3 Resolved spectral analysis

3.1 Spectral index maps and radio colour-colour analysis

To study the resolved spectral trends within 3C295, we produce two spectral index maps using the Broadband Radio Astronomy Tools (BRATs; Harwood et al. 2013, 2015), shown in Fig. 3: one between 132 MHz and 1.658 GHz, and the other between 4.76 GHz and 14.975 GHz. These are made using only pixels with a flux density value above 15σ.

It is immediately apparent that the overall 0.132–1.658 GHz spectral index values are significantly flatter than the 4.76–14.975 GHz spectral index values. The observed spatial distribution trends of the spectral indices are consistent with those typically observed in classical FRII sources (Orrù et al. 2010; McKean et al. 2016; Harwood et al. 2017): the flatter indices correspond to the radio hotspots, where particle acceleration takes place, and the overall steepening of the source is in the direction of the host galaxy. In general, we find an overall flattening of the spectral indices towards the lower frequencies, where the spectrum gets closer to the injection value and absorption processes start playing a role.

We also use these spectral index maps to perform a radio colour-colour analysis (Katz-Stone et al. 1993). These plots show the distribution of spectral indices between one pair of low frequencies (α_low) versus the spectral indices of another pair of higher frequencies (α_high). This provides crucial insights into the spectral shape of radio sources. Furthermore, the presence of emission due to mixing of different electron populations and regions with different magnetic field strengths or radiation losses lead to different curves in the colour–colour plot (e.g. Rudnick 2001; Rajpurohit et al. 2020, 2021). Radio colour-colour plots are similar to spectral curvature maps - the further a data point lies from the unit line (where α_low= α_high), the more spectral curvature is measured for that point. Similarly, the stronger the curvature, the older the age along a given spectral shape.

We first create a colour-colour plot extracting the spectral index from square shaped boxes of width 0.6′′, which correspond to 3.6 kpc in physical units. This size was chosen because it is equivalent to the beam size of our observations. These boxes are shown in Fig. 4b. The low frequency spectral index values are extracted between 132 MHz and 1.658 GHz and the high frequency one between 4.760 and 14.975 GHz. The curvature is negative for a convex spectrum. The resulting plot is shown in Fig. 4a, and the relevant regions are shown in Fig. 4b. All data points lie below the unit line (solid black line, where ), which indicates a clear overall negative curvature for 3C295. All regions seem to follow a single, continuous trend in the colour-colour plot, similar to trends observed for other radio galaxies in the literature (e.g. Figs. 3, 10 and 5 in Katz-Stone et al. 1993; Brienza et al. 2020; Shulevski et al. 2015, respectively). The two data points furthest from the power-law line – which are therefore associated with the oldest emission in the source – correspond to the South-Westernmost boxes in the North lobe: in other words, with the regions nearest to the host galaxy.

Because we only have access to one box per hotspot, we create a more illustrative pixel-based colour-colour plot, which follows the same trends, albeit more noisily. This allows us to more finely follow the spectral properties of the source. The resulting plots are shown in Figs. 4c,d. It is evident that the spectral shape of the Northern and Southern hotspots are different than the Northern and Southern lobes.

To better understand the curvature distribution within the individual components of 3C295, we superimposed the observed data with standard spectral ageing models. Indeed, radio colour-colour plots have the advantage that the standard spectral models depend only on the injection index, and are independent of magnetic field strength; adiabatic expansion, compression; or radiative losses (Katz-Stone et al. 1993; van Weeren et al. 2012). They represent an easy way to visualise the ageing models and to trace back the data to injection properties. The extrapolation of the trajectory in the colour-colour plot to the power-law line allows us to estimate the injection index of radio source.

We consider two standard spectral models. In the Jaffe-Perola (JP) spectral model (Jaffe & Perola 1973), the electron population of a given region is accelerated once, resulting in a power-law energy distribution and therefore a power-law flux density profile. This model also assumes that electrons undergo a series of scattering events which randomise their pitch angle relative to the magnetic field lines. The continuous injection (CI) model (Pacholczyk 1970) assumes instead, as the name implies, ongoing injection of fresh particles into a given region.

As a reference, we overlay four curves to the data, which represent JP models with injection indices equal to − 0.50, − 0.70, and − 1.0, and a CI model with injection index equal to − 0.50. The − 0.50 is the minimum injection index commonly used as a lower bound and the other two values demonstrate how the curve wouldshift in the plot when increasing the injection index. We note that most of the points related to the source’s lobe in Figs. 4a,c can be described with a JP model with injection index of approximately − 0.70.

Interestingly, the data points associated with 3C295’s hotspots have a significantly different distribution from the lobes (Fig. 4d) with values varying mainly in the range of − 0.3 < α < −0.5. This distribution is inconsistent with the JP and CI models. Such flat spectral indices could be indicative of synchrotron self-absorption, thermal absorption, or a low-energy cut-off (e.g. Katz-Stone & Rudnick 1997); in this latter case, given the magnetic field of the hotspots constrained by equipartition & SSC models (Harris et al. 2000), a low-energy cutoff or flattening in the hotspot electron spectrum would be present at γ ~ 300. The hotspot data distribution in the colour-colour plot can thus best be understood as evidence of one or more of these low-frequency absorption processes taking place within the hotspots.

It is worth mentioning that some points related to the lobes shown in Fig. 4c appear to be located in the same area of the plot where the hotspot points lie (on the left side of the JP model with injection index −0.5). These points are likely coming from regions where there is an overlap between the lobe and hotspot flux density contribution.

In conclusion, the spectral trend in the lobes is consistent with a JP model, indicative of particle radiative cooling. To the contrary, hotspot data provides evidence for the presence of absorption processes, as we further investigate in the next section.

Fig. 2 Regions associated with each component: inner ellipses are the “hotspots”, outer circles are the “lobes”. These are used in combination with a 5σ flux density threshold, indicated by the overlays in Fig. 1.

Fig. 3 Spectral index maps with associated error maps. Panel a: 0.132–1.658 GHz spectral indices. Panel b: 0.132–1.658 GHz spectral index error map. Panel c: 4.76–14.975 GHz spectral indices. Panel d: 4.76–14.975 GHz spectral index error maps.

3.2 Component spectral analysis

To further investigate the signs of absorption in the low-frequency spectra of the hotspots identified in the previous section, we perform a spectral fit to the radio spectral energy distributions of the hotspots. Spectral modelling is often done using BRATs, but this software package has not yet implemented low-frequency absorption models. We thus fit our hotspot radio spectral energy distributions using two classical absorption models: synchrotron self-absorption and free-free absorption.

Using the morphological components defined in Fig. 2, we integrate the flux density in our four morphological regions, subtracting the hotspot region flux density from the total lobe region flux density for both North and South. We assume that the underlying spectrum before absorption does not change over time or within a given region, which is necessary to reduce the number of degrees of freedom to fewer than the number of data points for all our absorption models – synchrotron self-absorption (SSA), free-free absorption (FFA), and power law (PL) spectra. The latter is only used as a point of comparison.

Here we briefly outline the spectral models used, given in Eqs. (1)–(3) below. Further details on these models can be found in Callingham et al. (2015), for example.

For the power-law spectral model, given by (1)

there are two free parameters: a flux density normalisation S₀ (corresponding to the flux density of the source at 150 MHz in our case) and the spectral index (α_PL). It is therefore well-constrained by our data. From Fig. 4, we expect this model to be a poor fit for all components – no part of our source lies on the unit line of our colour-colour plot. It is therefore used as a point of comparison and contrast with the other models.

The free-free absorption model, given by: (2)

has 3 free parameters: flux density normalisation S₀, spectral index α_FFA, and the absorption frequency ν_c, which is the frequency atwhich the optical depth is unity. This model is also well-constrained by our available data.

Finally, for the synchrotron self-absorption model we have (3)

Equation (3) describes synchrotron self-absorption with 3 free parameters: S₀, the flux density at 150 MHz, ν_c, the absorption frequency, which is the frequency at which the optical depth is unity, and β, the power law index of the energy distribution of the emitting electron population. We note that β is related to the synchrotron spectral index as .

We perform fits using scipy’s curve_fit function, which performs a Levenberg-Marquardt fit between our function and our data points. These fits are then used as the initial parameters of our Markov-chain Monte-Carlo (MCMC) fitting method. We employed an ‘affine invariant’ MCMC ensemble sampler (Goodman & Weare 2010) as implemented by the EMCEE package (Foreman-Mackey et al. 2013), creating 500 ‘walkers’ to fit the data, with the 1σ confidence interval for each of our fits used as a final error estimate. The points and the best-fit curves are all shown in Fig. 5. The best-fit parameter values and their errors are summarised in Table 3.

Figure 5 shows the spectral behaviour of our 4 components along with the spectral model best-fit curves. We see that the SSA and FFA models both seem give good fits to the data. However, only the LOFAR data point lies below the turnover frequency, which means that we lack the low-frequency constraints to discriminate between FFA and SSA.

We also see that the power law fit fails to describe the low-frequency properties of the hotspots, or the high-frequency properties of the lobes. This is also apparent in the power-law fit values, which are systematically higher than 1 – the model is not complex enough to capture the spectral behaviour of the data.

For all components, the FFA and SSA spectra are in very close agreement at high frequencies. The flux normalisation at 150 MHzvaries between the fits for all components, but this seems to be strongly tied to the un-physically sharp decrease in the best-fit SSA spectra at lower frequencies (see Fig. 5). An additional ILT measurement at 55 MHz will allow further resolved spectral analysis of 3C295. This will be more challenging due to the lower resolution that comes with lower frequencies, but it will help that the integrated spectrum of 3C295 peaks at about 50 MHz (Scaife & Heald 2012). The lack of a low-frequency constraint also appears when we look at the integrated flux and best-fit models, shown in Fig. 6. There, we show the total flux of 3C295 at each frequency, and plot the sum of the best-fit curves for both the lobes (power law) and the hotspots (one line per model, FFA or SSA, specified in the figure). We see that, overall, the flux drop-off of the SSA model is mitigated, and that both models appear to fit the data equivalently well: clearly, a further low-frequency constraint is necessary.

Fig. 4 Radio colour-colour plots of different regions of 3C295. The spectral index values of Fig. 3a are plotted in the x-axis and the values of Fig. 3c in the y-axis. Panel a: radio colour-colour plot of 3C295 overlaid with the JP and CI spectral aging models with different injection indices. The regions used for extracting the spectral indices are shown in Fig. 4b. The two North lobe outliers correspond to regions near the host galaxy. Panel b: per-pixel colour-colour map to the lobes of Fig. 3. We see that the emission is largely consistent with the JP model. Outliers are pixels which are near the hotspots, and therefore likely to be contaminated by hotspot regions. Panel d: Per-pixel colour-colour plot of the hotspot data of Fig. 3. We see that the emission in hotspots is inconsistent with all three ageing models.

Fig. 5 Per-component spectral fitting, using synchrotron self-absorption, free-free absorption and simple power law models. The model colours are the same in all plots. The dashed line in the lobe plots represents a power-law fit done using only the two lowest-frequency data points.

4 Conclusion and future work

In this paper, we have used the sub-arcsecond capabilities of the ILT to create a spatially resolved map of 3C295 at 132 MHz. Combined with archival data from the VLA and MERLIN at GHz frequencies, we analyse four components: the Northern and Southern hotspots and two lobes. While a low-frequency turnover had previously been identified in the integrated spectra for 3C295, here we clearly show that this is due to low-frequency absorption processes in the hotspots. The lobes remain consistent with un-absorbed, optically thin synchrotron emission, which can be described by a JP particle ageing model.

While the hotspot spectral measurements are consistent with the presence of the low-frequency absorption, we are unable to distinguish between synchrotron self-absorption and free-free absorption. Doing so will require independent observations at lower frequencies. ILT Low Band Antenna (LBA) observations at 55 MHz would provide a strong constraint on the low frequency spectral behaviour, but currently the degradation in resolution by a factor of 2 when compared to the HBA means we cannot differentiate between the lobe and hotspot emission. The ILT is expected to expand in the coming years, notably with the recent addition of a Latvian station to the array, and the immenent addition of an Italian station. These will further improve uv-coverage and resolution of the ILT.

Finally, in order to maximise our scientific returns, the following advances, in our estimation, be made; they are listed in order of expected difficulty and significance:

• Formalise the limitations imposed by the pixel statistics in interferometric images for spectral analysis,

• Formalise the limitations imposed by joint deconvolution processes in very-large-bandwidth analysis,

• Formalise the limitations imposed by the conformity of uv-coverage across multi-frequency observations.

Sub-arcsecond imaging with the ILT is paving the way for spatially resolved spectral modelling of compact radio-loud AGN such as 3C295. Future implementation of low-frequency absorption models to BRATs and elsewhere will allow physically motivated modelling that already takes into account the above list of limitations. By continuing to work on LoTSS post-processing at sub-arcsecond resolution (Morabito et al. 2022), we will be able to perform this kind of spectral modelling for large numbers of sources in the near future.

Fig. 6 Integrated spectral energy distribution, created by taking the sum of the component fluxes at each frequency, and the sum of the power-law best-fit curves of the lobes with the FFA and SSA best-fit curves, respectively, to show the overall fit tothe source in both cases.

Acknowledgements

The authors thank F. Massaro for fruitful discussions over the course of this work, which helped improve the form and contents of this paper. E.B., A.B., M.B., C.J.R. acknowledge support from the ERC-Stg grant DRANOEL, n.714245. K.R. and M.B. acknowledge financial support from the ERC Starting Grant “MAGCOW”, no. 714196. L.K.M. is grateful for support from the UKRI Future Leaders Fellowship (grant MR/T042842/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. The VLA is operated by the US National Radio Astronomy Observatory. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. Based on observations made with MERLIN, a National Facility operated by the University of Manchester at Jodrell Bank Observatory on behalf of STFC. We thank Anita Richards and the e-MERLIN team for the maintenance and automatic processing associated with the MERLIN archive. R.Jv.W. acknowledges support from the ERC Starting Grant ClusterWeb 804208.

Appendix A Initial calibration model of 3C295

Fig. A.1 VLA map of 3C295 at 8.561 GHz, courtesy of M. Hardcastle. Pixel size is 0.05′′, beam size is 0.175′′, and flux is in Jy. The image is overlaid on itself in yellow contours starting from 0.5σ and going to 0.9 of the maximum flux present. 3C295 is clearly well-resolved and shows complex structure with both compact and diffuse emission present.

References

All Tables

All Figures

Fig. 1 High resolution images of 3C295 at different frequencies. Flux density units are Jy bm−1, with colour scales starting at 5σ and tuned to show the source morphology. Each image includes an overlay of itself, with 10 levels starting at 3σ and increasing exponentially until reaching the maximum value in the image. The restoring beam is 0.6′′ in all cases, and the general properties are summarised in Table 2. Panel a: LOFAR, 132 MHz. Panel b: MERLIN, 1.658 GHz. Panel c: VLA, 4.760 GHz. Panel d: VLA, 14.975 GHz.

In the text

	Fig. 2 Regions associated with each component: inner ellipses are the “hotspots”, outer circles are the “lobes”. These are used in combination with a 5σ flux density threshold, indicated by the overlays in Fig. 1.
In the text

	Fig. 3 Spectral index maps with associated error maps. Panel a: 0.132–1.658 GHz spectral indices. Panel b: 0.132–1.658 GHz spectral index error map. Panel c: 4.76–14.975 GHz spectral indices. Panel d: 4.76–14.975 GHz spectral index error maps.
In the text

Fig. 4 Radio colour-colour plots of different regions of 3C295. The spectral index values of Fig. 3a are plotted in the x-axis and the values of Fig. 3c in the y-axis. Panel a: radio colour-colour plot of 3C295 overlaid with the JP and CI spectral aging models with different injection indices. The regions used for extracting the spectral indices are shown in Fig. 4b. The two North lobe outliers correspond to regions near the host galaxy. Panel b: per-pixel colour-colour map to the lobes of Fig. 3. We see that the emission is largely consistent with the JP model. Outliers are pixels which are near the hotspots, and therefore likely to be contaminated by hotspot regions. Panel d: Per-pixel colour-colour plot of the hotspot data of Fig. 3. We see that the emission in hotspots is inconsistent with all three ageing models.

In the text

	Fig. 5 Per-component spectral fitting, using synchrotron self-absorption, free-free absorption and simple power law models. The model colours are the same in all plots. The dashed line in the lobe plots represents a power-law fit done using only the two lowest-frequency data points.
In the text

	Fig. 6 Integrated spectral energy distribution, created by taking the sum of the component fluxes at each frequency, and the sum of the power-law best-fit curves of the lobes with the FFA and SSA best-fit curves, respectively, to show the overall fit tothe source in both cases.
In the text

	Fig. A.1 VLA map of 3C295 at 8.561 GHz, courtesy of M. Hardcastle. Pixel size is 0.05′′, beam size is 0.175′′, and flux is in Jy. The image is overlaid on itself in yellow contours starting from 0.5σ and going to 0.9 of the maximum flux present. 3C295 is clearly well-resolved and shows complex structure with both compact and diffuse emission present.
In the text