© ESO 2022

1 Introduction

As the first quasar discovered, 3C 273 (J1229+0203) is one of the most well-studied objects of its type at almost every waveband with over a thousand peer-reviewed publications aiming to understand its properties. At its centre 3C 273 is observed to host a supermassive black hole (SMBH) of mass (890 ± 190) × 10⁶ M_⊙ (Peterson et al. 2004) with a complex one-sided jet observed on parsec scales (Zensus et al. 2020). The possibility that 3C 273 is intrinsically one-sided was considered during early investigations (Davis et al. 1985), but given active galactic nucleus (AGN) unification theory (Urry & Padovani 1995) and observations of other AGN, there is a general consensus in recent literature that 3C 273 is two-sided and intrinsically symmetrical observed at an angle of ≈ 5° (Meyer et al. 2016) with a jet-to-counterjet flux ratio greater than 5300 (Sambruna et al. 2001).

The core flux density is known to be highly variable at and above gigahertz (GHz) frequencies with the flux density at 37 GHz varying in the range of 10–57 Jy (Sukharev et al. 2016). Intraday variability has also been detected (Kalita et al. 2015), but no studies to date have been able to monitor time variability in the megahertz (MHz) regime due to the limitations currently imposed at these frequencies (e.g. resolution, artefact suppression, flux scale accuracy). The radio emission from the core is understood to be optically thick synchrotron radiation and, despite the observed variability, the radio spectrum is known to be flat between 0.1 and 100 GHz.

On larger scales, a faint inner jet extending 12 arcsec (33 kpc) joins the nucleus to an outer jet (Davis et al. 1985; Conway et al. 1993), which can be seen in the Karl G. Jansky Very Large Array (VLA) observations shown in Fig. 1. In addition to the radio waveband, this inner jet has been detected at optical wavelengths by the Hubble Space Telescope (HST; Martel et al. 2003). In contrast, the outer jet is bright and highly polarised, extending 27 kpc in the plane of the sky (Perley & Meisenheimer 2017). The jet has a knotty structure at radio, optical, and X-ray frequencies indicating modulation in the jet’s kinetic power. Georganopoulos & Kazanas (2004) also note a decrease in the ratio of the X-ray to radio brightness with distance from the core; in other words, as the radio jet increases in brightness with distance from the core, the X-ray jet becomes fainter.

Despite extensive research, there is no consensus regarding the origin of the X-ray emission associated with the relativistic jet (Georganopoulos et al. 2006) and there are varying interpretations of the multi-wavelength spectra for the knots in the jet of 3C 273. However, it is generally agreed that for each knot, a population of relativistic electrons gives rise to the radio to infrared emission via the synchrotron process. Radio and optical polarisation measurements are consistent in magnitude and orientation, which is further evidence that the emission across these wavebands is non-thermal in origin and stems from a single population of electrons (Roeser & Meisenheimer 1991). However, the mechanism behind the ultraviolet to X-ray emission, and how it fits into the overarching multi-wavelength picture, is less well understood.

While X-ray emission from Fanaroff and Riley class I radio galaxies (FR-Is; Fanaroff & Riley 1974) and BL Lacs is usually explained by extrapolating the radio to infrared synchrotron spectrum, in general the X-ray flux of Fanaroff and Riley class II radio galaxies (FR-IIs) and quasars are higher than expected. Observations of the knots in the jet show that 3C 273, a quasar believed to have an underlying FR-II morphology, is no exception with simple power law models underestimating the X-ray flux. To account for this shortfall, the favoured model for the X-ray emission has typically been inverse-Compton scattering of the cosmic microwave background (IC/CMB; Sambruna et al. 2001). However, several challenges to this interpretation have been made. Firstly, models seeking to reproduce the X-ray flux solely via IC/CMB require relativistic bulk motion on kiloparsec scales. Superluminal speeds (~ 10c) on parsec scales near the core as measured by Very Long Baseline Interferometry (VLBI) studies open the door to this possibility, but whether these speeds are maintained on larger scales is unclear. Arshakian & Longair (2004) suggest that jets should decelerate to mildly relativistic speeds on kiloparsec scales due to the entrainment of ambient gas, a finding later supported by the modelling of Mullin & Hardcastle (2009) who verify that jets decelerate, constraining the bulk Lorentz factor, Γ, for kiloparsec-scale radio jets of quasars to 1.18 < Γ < 1.49.

Ultimately, an IC/CMB origin of the X-ray emission was effectively ruled out by observational γ-ray results. Meyer et al. (2015) showed that IC/CMB models violate the upper limits on the γ-ray flux of 3C 273 and PKS 0637-752 at the 99.99% confidence level in more than one Fermi-LAT energy band. These deep upper limits at gigaelectronvolt (GeV) energies constrain the Doppler beaming factors, implying that the jet is not highly relativistic on kiloparsec scales, as is required for IC/CMB models.

One alternative to the upscattering of the CMB is the synchrotron self-Compton (SSC) processes that may also plausibly produce X-ray emission in the jet. However, Sambruna et al. (2001) find that for the X-ray component to be more luminous than the radio component, as is observed, weak magnetic fields (B ≲10⁻³ mG) orders of magnitude below the equipartition value are required, unless the jet is significantly debeamed (δ ≈ 0.3–0.5). Such significant debeaming is incongruous with the jet-to-counterjet flux ratio of 5300 (Sambruna et al. 2001), and so SSC is unlikely to be solely sufficient to produce the X-ray fluxes that are observed. A recent study by Wang et al. (2020) comprehensively demonstrated that the synchrotron, IC/CMB, and SSC components of a single population of electrons are unable to reproduce the expected X-ray emission. While IC/CMB and SSC are therefore likely to occur at some level in the jet of 3C 273, they are not the primary mechanisms behind the X-ray emission. Given that such ‘one-zone’ models struggle to explain high-energy emission not only in knots, but in the jets of FR-Is more generally (e.g. Harwood & Hardcastle 2012 and references therein), more complex models are likely required if one is to accurately describe the processes that drive the observed emission.

Many competing hypotheses have been proposed to resolve these inconsistencies. Cara et al. (2013) favour a second population of particles, whether leptonic or hadronic, to explain the X-ray emission for the knots. Wang et al. (2020) also found that the multi-wavelength spectra can be accurately reproduced by including synchrotron processes associated with a second population of particles. It is unclear, however, whether this second population consists of electrons or protons as the SED predictions in both cases satisfy the observations they considered. Aharonian (2002) proposed a broken power law spectrum of accelerated protons (e.g. from the jet shear boundary) that produce synchrotron radiation as the origin of the X-ray and γ-ray flux from the knots. However, super-Eddington jet powers or strong (~10 mG) magnetic fields are required in the hadronic case as more energy is needed to accelerate the protons sufficiently. A secondary leptonic population would overcome this issue, although the origin of such a second population is unknown.

It has also been suggested that two populations of electrons could have different acceleration mechanisms. Liu et al. (2015) proposed separate shock regions in the knots for the radio and X-ray emission with distinct populations formed upstream and downstream because of the shock compression effect, where the maximum energy of the injection spectrum into the downstream emission region is lower than the maximum energy of the upstream region. Assuming that particle acceleration in the knots occurs in situ (i.e. that the acceleration and emission zones are co-spatial), testing whether the radio and X-ray knots are spatially coincident can provide evidence of multiple acceleration sites, and hence multiple populations. Marchenko et al. (2017) found that the X-ray knots are upstream of the radio knots by 0.2–1.0 kpc. Alternatively, Liu et al. (2017) have proposed shear acceleration to produce high-energy electrons, but a rigorous comparison between the predictions of such models has yet to be performed.

As highlighted by Sambruna et al. (2001), an underlying issue exists with the ad hoc nature of invoking a second component in that it produces a scenario where there are enough free parameters that the physical state of the jet is essentially unconstrained without comprehensive observations, as complex models invariably provide better fits to data due to overfitting. Therefore, despite extensive research, there is no consensus regarding the origin of the X-ray emission associated with the relativistic jet.

One area of investigation that has yet to be fully explored is that of low-frequency emission. Historically, observations in the MHz regime have lacked either the resolving power and/or UV coverage to investigate the knot and diffuse jet emission at long wavelengths. However, the international baselines provided by the LOw-Frequency ARray (LOFAR; van Haarlem et al. 2013) provide the opportunity to investigate this previously unexplored area for the first time, and provide a baseline reference from which to measure the proper motion of the knots at low frequencies. In this paper, we present high-resolution images at 150 MHz, the first sub-arcsecond images of 3C 273 below 1 GHz, to investigate five key unresolved questions: What compact structures are observable in the jet of 3C 273 at low frequencies? What processes result in the observed low-frequency emission of the jet knots? How does the jet’s kinetic power determined at low frequencies compare to previous estimates? Is the previously unseen counterjet detected at low frequencies? What constraints can be placed on the bulk speed of the jets and their inclination angle?

In identifying the features of the jet, we adopt the naming convention of Jester et al. (2005), with the addition of features fs, Ex1, and H1 (Fig. 2). Throughout this paper we define the spectral index such that S ∝ ν^−α and use a concordance model in which H₀ = 71 km s⁻¹ Mpc⁻¹, Ω_m = 0.27, and Ω_Λ = 0.73 (Spergel et al. 2003).

Fig. 1 VLA C and P band images of 3C 273. Left: VLA C band at 4885 MHz at a resolution is 0.5 arcsec. Contours are such that Scont = 2 × 2n mJy beam−1 with a peak brightness of 35.99 Jy beam−1. The inner jetis marked for clarity. Right: VLA P band at 327 MHz at a resolution of 0.5 arcsec. Contours are such that Scont = 80 × 2n mJy beam−1 with a peak brightness of 36.37 Jy beam−1. These data are reproduced from Perley & Meisenheimer (2017).

Fig. 2 HST image of the 3C 273 outer jet masked at 0.7 e−∕ s to only show the jet. The direction of motion is left-to-right, with the core ~ 11 arcsec to the left of the image. In1, In2, and Ex1 are galaxies unrelated to the jet; fs is a foreground star. The naming convention is from Jester et al. (2005), with the addition of features fs, Ex1, and H1.

2 Observations and data reduction

2.1 LOFAR observations

The quasar 3C 273 was observed during LOFAR cycle 8 using 24 core stations, 14 remote stations, and 12 international stations (sixstations in Germany, three in Poland, and one in each in France, the UK, and Sweden). The observations were made using a four-hour track when 3C 273 was highest in the sky. This is particularly desirable when using LOFAR due to the source’s low declination (+ 02°), where increased ionospheric decorrelation of the signal increases the difficulty in achieving a robust calibration of the source. Since 3C 273 is radio bright, these shorter observations will provide a sufficient signal-to-noise ratio to image the source with the added benefit of increasing the likelihood of observing during optimal ionospheric conditions.

3C 295 and 3C 280 were chosen as calibrators. While 3C 295 and 3C 280 are 50° and 46° from 3C 273, respectively, as they were only used to calibrate the direction-independent effects it does not significantly impact the overall robustness of the calibration. 3C 295 is a bright standard flux density calibrator at low frequencies (Scaife & Heald 2012) and situated at a high declination. It is therefore an ideal calibrator source. A summary of the observations, including the observationsof the two calibrators, is given in Table 1.

2.2 Preprocessing

The raw data were initially processed by the observatory in the standard manner. The preprocessing pipelines automatically flagged RFI, averaged the data in time and frequency, and performed demixing of the A-team sources (Cassiopeia A, Cygnus A, Hydra A, Taurus, Virgo A) where required. Due to the large size of LOFAR observations (4–20 terabytes per pointing), the data were averaged from the raw time-resolution of one second and 64 channels per sub-band (ch/sb) to 2 s and 16 ch/sb providing a bandwidth of 12.2 kHz per channel. This reduced the computational requirements required for the data reduction and analysis, and had negligible impact on the science cases presented in this paper.

2.3 Calibrator processing for direction independent effects

Upon completion of the preprocessing, the data were downloaded and initial calibration performed using the LOFAR pre-facet calibration pipeline (PREFACTOR¹; de Gasperin et al. 2019) on the core and remote stations. As only one calibrator can be used by the pipeline, the PREFACTOR calibrator pipeline was first used to process the 3C 295 and 3C 280 data.

The solutions from the pipeline runs were then inspected. The LOFAR station bandpass solutions associated with the 3C 280 run had an incorrect shape compared to the theoretical bandpass (likely due to calibrator model inaccuracies), and so 3C 295 was used for the remainder of this analysis. A two-component model of 3C 295 is packaged with PREFACTORwhich, while sufficient for studies that only use the core and remote stations, is not suitable for long-baseline analysis. A high-resolution model provided by the LOFAR VLBI working group, consisting of CLEAN components (Högbom 1974) was therefore used.

The observations were put through the PREFACTOR pipeline. The calibrator phase solutions are shown in Fig. 4 where CS001HBA0 is the reference station. All other core stations have small phase offsets with respectto this reference station. The remote stations show more variation than the core stations with the international stations showing the fastest variations with frequency. The stations in this figure are sorted by baseline length from the reference station. An interactive map showing the station locations can be found online² The data were then compressed using DYSCO³ (Offringa 2016), allowing LOFAR data to be compressed by approximately a factor of six with only a one per cent increase in system noise for typical correlator time and frequency resolutions.

Before calibrating the international stations, the data were processed through the PREFACTOR target pipeline, which flags the data, removes contributions from bright off-axis sources, and calibrates the phases on the core and remote stations. Following this, the data areready for international station calibration.

Fig. 3 HST image of 3C 273 at 850–1700 nm with a peak brightness of 28606 e− s−1. The image combines four exposures providing a 37-s integration time.

2.4 International station calibration

The LOFAR long-baseline pipeline (Jackson et al. 2016, 2022; Morabito et al. 2022) was run to calibrate the international stations. The pipeline takes the 3C 273 data and the PREFACTOR solutions as input, specifically the gain solutions from the calibrator pipeline and the phase solutions from the target pipeline.

Based on a 6 arcsec resolution image, the hotspot of 3C 273 (i.e. the jet terminus) was determined to have the largest integrated flux density. The pipeline automatically determined this to be the best in-field calibrator, and after applying the PREFACTORcalibration solutions the data were phase-shifted to the position of this component. The core LOFAR station beams were then added coherently to form station ST001. A phase-only self-calibration of the in-field calibrator was then performed against an initial model. The pipeline, which is optimised to self-calibrate compact calibrators with simple structure from the LBCS (Jackson et al. 2022), attempts to initialise self-calibration with a point source. In this case, a model using a single Gaussian with major and minor full width at half maximum (FWHM) values of 0.1 arcsec was initially used, but the results were poor, given that a point source model neglects the large-scale jet. An initial model of the whole source was therefore built using PYBDSF⁴ from a lower resolution image (~1 arcsec) produced from four hours of LBA data (Groeneveld et al. 2022). All of the calibration solutions were then applied to the data.

Fig. 4 Calibrator phase solutions for 3C295 for the XX polarisation. The reference station is CS001HBA0, and the plots increase by baseline length from the reference station reading from left to right. All plots are on the same time and frequency axes. The first remote station is labelled, as are individual international stations.

Fig. 5 Components to the 3C 273 spectral energy distribution. Data are from Perley & Meisenheimer (2017), except for the 151 MHz flux density of the halo, which is from Punsly & Kharb (2016). The flux density of the jet includes all knots. The halo refers to the diffuse cocoon of radio emission encompassing the jet (A, B, and R of Fig. 9). The radio spectral indices of the core, jet, and halo are 0.30 ± 0.07, − 0.81 ± 0.03, and − 0.79 ± 0.11, respectively. The spectral indices of the jet and the halo are consistent with optically thin synchrotron emission, as expected. The LOFAR HBA bandwidth is shown by the shaded region.

2.5 Self-calibration and imaging

From the initial model for the brightness distribution of 3C 273 used in Sect. 2.4, the relative gains (amplitude and phase) were calibrated. The calibration solutions were applied to the data and were then imaged, which provided an improved model for the observed sky surface brightness distribution. WSCLEAN was used for imaging stage (Offringa et al. 2014), which uses both the multi-scale clean algorithm (MSCLEAN; Cornwell 2008; Rau & Cornwell 2011) to improve the diffuse emission fidelity, and multi-frequency synthesis (MFS; Rau & Cornwell 2011) to account for spectral variations across the bandwidth. These steps were iterated in the standard manner until the image produced did not result in an improvement in quality⁵.

Fig. 6 LOFAR images of 3C 273 at 151 MHz with a beam size of 0.31 × 0.21 arcsec. Contours levels are set such that Scont = 40 × 2n mJy beam−1. Shown are the full source (left) and the outer jet with the components labelled (right).

2.6 Astrometry and flux-scale bootstrapping

Attaining absolute astrometric positions with LOFAR is inherently difficult due to ionospheric refraction causing sources to shift by several arcseconds in the imageplane at 150 MHz. Accurate total electron content (TEC) solutions should mostly correct for this effect, but, combined with other factors such as positional shifts during self-calibration, additional alignment was also performed to ensure that robust astrometric positions were achieved.

In order to achieve robust astrometric positions, the core of 3C 273 was therefore aligned for the LOFAR, VLA, and HST observations to the position reported by Johnston et al. (1995), which is accurate to better than 3 mas. For the radio data the core position was taken to be the pixel of the peak brightness; for the near-infrared data the intersection of the diffraction spikes was used to define the core position. For AGN the position of the core can drift across frequencies; however, this effect is seen on milliarcsecond scales, and so is unlikely to have a significant impact on our analysis. Given the resolution of the LOFAR, VLA, and HST data, we conclude the accuracy of the astrometry in the composite images is < 0.05 arcsec.

Initially, the flux calibrator 3C 295 and the Scaife & Heald (2012) flux scale were used to set the flux scale in the LOFAR observations. However, imperfect beam models mean that the relative gains between 3C 295 and 3C 273 are not well determined, and introduce an error on the absolute flux density scale. To determine the scaling factor required to correct for these errors, the flux measurements of Perley & Meisenheimer (2017) and Punsly & Kharb (2016) were used to perform a linear least squares regression and ascertain the spectral indices (Fig. 5). PYBDSF was then run on the LOFAR image to model the source and to find the total flux density of the core component. The fit of the core at 150 MHz was then used to bootstrap the flux scale in the LOFAR data.

Due to the relatively high surface brightness of 3C 273, maximising image resolution rather than image sensitivity was the primary driver in determining the imaging parameters. Uniform weighting (robust − 2) was therefore chosen to produce final images at the highest resolution, suitable for determining the properties of the small-scale structure (e.g. knots). This provided a beam size of 0.31 × 0.21 arcsec with a position angle of 164°. The resulting radio maps are shown in Fig. 6.

2.7 Ancillary data

Higher frequency radio and near-infrared data were used to aid the LOFAR analysis and interpretation. The fully calibrated VLA observations of Perley & Meisenheimer (2017) were retrieved, which include images from 0.3 to 43 GHz at a range of resolutions. Near-infrared data taken with the Hubble widefield camera 3 (WFC3) were also obtained from the Hubble Legacy Archive⁶ (Jenkner et al. 2006) and are presented in Fig. 3.

3 Results and discussion

3.1 Morphology and brightness distribution

The overall large-scale morphology of 3C 273 shown in Fig. 6 is in good agreement with previous studies performed with the VLA (e.g. Perley & Meisenheimer 2017). While this suggests that an overall good level of image fidelity has been achieved, some unusual features on smaller scales are present in the vicinity of the core region, particularly in the north–south direction. Due to the low declination of 3C 273, the uv tracks for 3C 273 are almost perfectly horizontal, hence any slowly varying errors of individual baselines will result in artefacts in the north–south direction, as is observed in Fig. 6. These features are therefore excluded from our analysis.

Determining the RMS noise from a blank region of sky well away from the source, we find that σ_rms= 4.8 mJy beam⁻¹. The peak brightness is located in the core region with a peak brightness of 7.8 Jy beam⁻¹ giving a dynamic range of 1614. This is lower than that obtained by de Gasperin et al. (2012) who achieved a dynamic range of ~ 5000 for their (non-international baseline) LOFAR observations of M87. While the sources are relatively close on thesky (with both located in the Virgo constellation), the declination of M87 is 10° higher than 3C 273, with the M87 study also using a longer integration time and larger bandwidth (116–162 MHz) than for the observationsused in this paper. The dynamic range is therefore at the expected level given the increased resolution provided by the inclusion of the international baselines, which is also impacted by the super-station, resulting in an increased noise level in the image due to the loss of the core-to-core station baselines. To estimate the theoretical noise, we account for the loss of core-to-core baselines, uniform weighting for imaging, the amount of data flagging, and the low declination of the source to arrive at a value of ~ 0.23mJy beam⁻¹. However, this value must be taken with a grain of salt as there are several important factors, which we have not taken into account, that could bridge the gap to the final σ_rms we measure. First, and most critical, is the fact that M87 is ~ 10 degrees away with a flux density in excess of 10³ Jy at 151 MHz.Although demixing was performed to remove the contributions of this off-axis source, the results are only as good as the calibration in the direction of this bright off-axis source, and we do not yet have a high-resolution model of M87 at these frequencies. Even optimistically, one per cent residuals would still be > 10 Jy. It is crucial to note that the core stations are combined into the super-station after demixing but before self-calibration; any effects from improper subtraction of M87 during the demixing process will be ‘baked in’ to the super-station, and these effects cannot be removed by subsequent self-calibration. Better demixing of M87 is not possible until we have a high-resolution model, and even then will be computationally challenging as the quality of demixing scales with the number of components in the model source. Second, the combination of core stations is expected to increase the theoretical noise slightly due to the combined station beam, which is not yet well mapped. This effect is at the level of approximately two per cent for a high-declination source, but we have not assessed this yet for low-declination sources. Finally, the elevation correction we implemented only accounts for the projection of the stations and does not include complications such as the increased ionospheric path length at low declination.

Shimwell et al. (2019) show the LOFAR flux calibration error for HBA observations to be ≈ 15 per cent. Due to the additional core flux scaling performed, the observations presented within this paper are likely to improve on this value, but we opt to use 15% to provide a conservative estimate. We therefore find a total flux density for the source of 91 ± 14 Jy. This agrees well with the previous VLBI study at 151 MHz by Arshakian et al. (2010) who find a total flux density of 97.95 ± 4.90 Jy. We are therefore confident that our flux calibration is robust.

The core has the highest peak brightness in the image (7.8 Jy beam⁻¹), which is 1.8 times that of the next brightest feature, the radio hotspot where the jet terminates (4.4 Jy beam⁻¹). The relative brightness of the features can be seen in Fig. 7. While the inner jet is not visible, a structure in the outer jet is apparent with all but one of the jet structures observed at optical wavelengths (Fig. 2) identified in the LOFAR image (Fig. 6, bottom). Knot B2 is the exception, for which we set an upper limit on its brightness of 40 mJy beam⁻¹. A summary of the measured flux densities for each of the knots is shown in Table 2.

Figure 8 shows the LOFAR images of the jet overlaid on HST data. The radio jet extends ~ 2 arcsec beyond the termination point of the jet with the shape of the radio jet closely following that of the optical jet. It is interesting to note that the galaxies In1 and In2 (Fig. 2), while detected with HST, are not present in the LOFAR observation, which suggests that these sources are radio quiet AGN.

Figure 9 shows the LOFAR images overlaid on VLA images at 325 MHz, 1.4 GHz, and 15 GHz. The 325 MHz P-band image is at 4 arcsec resolution and illustrates the resolution that has been attainable prior to LOFAR at similar frequencies. At the VLA resolution, low-frequency radio lobes are detectable on either side of the jet close to where the inner and outer jets meet (A and B in Fig. 9), and the lobe on the northern side of the jet (A) is particularly extended. For the VLA 1.4 GHz L-band images at 1.4 arcsec resolution, the inner jet is apparent in the VLA data (denoted P) and, while the radio lobe on the northern side of the jet is not detected (denoted Q), extended emission on the other side of the jet is present (denoted R). These inner jet features are not detected in the LOFAR images; however, at VLA U-band frequencies (15 GHz) which provide comparable resolution, the morphology of the jet is similar to that observed in the LOFAR data with knot B1 (denoted X) undetected in either image.

In the LOFAR data, we observe a minor extension to the south-west of the hotspot (denoted Z) with a peak brightness of 102 mJy beam⁻¹ that is not detected in the VLA observation. The resolution between LOFAR and VLA at U-band frequencies is similar, but there are two orders of magnitude between the image frequencies. Measuring the RMS noise of the VLA image to be 0.3 mJy beam⁻¹, the emission would require a spectral index steeper than 1.2 to be detectable with LOFAR, butnot in the VLA U-band image. From our observations it is not possible to deproject 3C 273, while the observed one-sided jet and a hotspot-like jet termination suggest FR-II as a plausible morphology. In such cases it has been shown that the injection index, the initial electron energy distribution at the point of acceleration, is commonly observed in the spectrum of FR-IIs to range between α_inj ≈ 0.6 and 1.0 (Harwood et al. 2013, 2015, 2016, 2017). Given that at least some spectral ageing (the preferential cooling of high-energy electrons resulting in spectral curvature) can be assumed, and the wide frequency range considered, a spectral index greater than 1.2 is certainly achievable. While spectral steepening will be less significant at frequencies below 15 GHz, causing this structure to be brighter in the P-band and L-band images, the resolution of the archival data is insufficient to discern structures on these scales. While it is therefore plausible that the extension is real, further investigation is needed to ascertain whether it is intrinsic to 3C 273, a background galaxy, or an imaging artefact.

Fig. 7 3C 273 ridgeline histogram with the brightness measured relative to the core. The image has been rotated counterclockwise by 42° such that the top of the image points north-west. The features around the core are likely to be imaging artefacts. The green dashed line shows the results of Davis et al. (1985), and the dotted black line the 3σ RMS noise of the LOFAR image.

Fig. 8 LOFAR data overlaid on HST near-infrared data. Black contours indicate the HST emission, and are set such that Scont = 0.7 × 2n e− s−1. The filled contours show the LOFAR dataset such that Scont = 40 × 2n mJy beam−1.

3.2 Spectral energy distributions of the knots

The resolution of our observations provides the first opportunity to investigate the spectrum of the jet knots at low frequencies. Regions coincident with the optical features observed in Fig. 2 were defined for the radio images (Fig. 6) and, combined with the VLA measurements of Jester et al. (2007), the spectral index for each region was calculated between 150 MHz and 8.3 GHz (Table 2). While the outer most regions (H1/H2) are in line with the expected spectrum of the hotspot and surrounding emission at low frequencies (Harwood et al. 2016), the origin of the spectrum of regions closer to the core are less clear. Three of the regions (A, B1, and C2) fall below the α = 0.5 physical minimum for first-order Fermi acceleration, which is currently the favoured mechanism by which particles responsible for the radio emission are primarily accelerated. Two possible causes for the observed emission are free-free absorption (FFA) and synchrotron self-absorption (SSA), which can lead to a turnover in the spectrum at low frequencies. In both cases attenuation of the emission is a function of frequency either through thermal absorption as a result of passing through an ionised screen (FFA) or because the brightness temperature cannot exceed the temperature of the non-thermal electrons (SSA). These models are discussed in detail in the context of radio-loud AGN elsewhere (e.g. Callingham et al. 2015; McKean et al. 2016), and so we do not repeat that process here, but for FFA the expected emission is given by (e.g. Kassim 1989) $$S_{\text{FFA},\nu }=S_{{\nu }_0}{\left(\frac{\nu }{{\nu }_0}\right)}^{-\alpha }\exp \left(-{\tau }_{{\nu }_0}{(\nu /{\nu }_0)}^{-2.1}\right),$$(1)

where $S_{{\nu }_0}$ and ${\tau }_{{\nu }_0}$ is respectively the flux density and optical depth at a reference frequency of 150 MHz and α is the optically thin spectral index. Similarly, assuming the emission region is homogeneous, the expected emission for the SSA model is given by (Kellermann 1966; Tingay & de Kool 2003) $$S_{\text{SSA},\nu }=a{\left(\frac{\nu }{{\nu }_{\text{p}}}\right)}^{-\alpha }\left(\frac{1-{\text{e}}^{-\tau }}{\tau }\right),$$(2)

where $$\tau ={\left(\frac{\nu }{{\nu }_{\text{p}}}\right)}^{(\delta -4)/2}$$(3)

with ν_p representing the frequency at which the region becomes optically thick and δ the power law index of the initial electron energy distribution such that α = (δ − 1)∕2. As measurements for the knots D2 and H3 are combined by Jester et al. (2007), as are B2 and B3 in this paper, we are only able to provide upper limits for these regions in addition to knot B1 which is below our detection limit. Regardless of these limitations, we find that fitting the FFA and SSA models is able provide a consistent interpretation of the observed spectra, the results of which are shown in Fig. 10 and Table 3.

The H1–H2 region is in line with expectations for being the primary acceleration region (i.e. the hotspot), being well described within the errors by a power law and optically thin at 150 MHz in the FFA (${\tau }_{{\nu }_0}=-0.13\pm 0.20$) and SSA (ν_p = 97.86 ± 2.47 MHz) models. Theoptically thin spectral index component is relatively consistent for both models and across all knots, being dominated by the GHz observations, which are well above the frequency at which the regions becomes optically thick. For regions closer to the core where the observed brightness is significantly lower than would be expected at 150 MHz, the models are able to provide a good description of the observed spectrum in all cases. For knot A, where we are not constrained by upper limits, we find our LOFAR observations are likely just below the frequency at which the region becomes optically thick (ν_p = 311.3 ± 17.2 MHz) with an optical depth of ${\tau }_{{\nu }_0}=2.45\pm 0.20$.

The wide variation in ${\tau }_{{\nu }_0}$ across the knot regions suggest that, if FFA is the mechanism by which absorption is occurring, the ionised screen cannot be homogeneous. While an inhomogeneous screen external to 3C 273 is possible, the scales on which the jet and knots are observed suggests theabsorption is likely occurring in situ. Such a local ionised medium is perhaps not surprising given the prominent nature of the observed knots. The jets of FR-IIs typically remain uninterrupted out to large distances and interaction, for example between the jet sheath and an inhomogeneous ionised medium, may explain the prominent knotty structure of the jet at radio frequencies.

With only one data point at low frequencies it is currently not possible to differentiate between the two models; however, both the FFA and SSA models are able to provide a good description for all the knots. We therefore suggest that, while the precise mechanism is not clear, absorption is present within the knot regions and should be carefully considered when undertaking future investigations. The planned upgrade of the LOFAR Low-Band Antennas (LBA) will provide observations at ~ 50 MHz that will better constrain these models and may also be able to differentiate between the two models for the most heavily absorbed knots. Ultimately, the Square Kilometre Array (SKA) – GMRT VLBI will provide the observations capable of constrainingthe spectrum at a few hundred MHz where the difference between SSA and FFA is most significant. Upcoming and future instruments, when combined with LOFAR results, will be a key to continuing this form of investigation.

3.3 Estimate of the diffuse emission

To determine the flux density of the diffuse emission surrounding the knots, we blanked the entire image at 150 MHz below a threshold of 5σ_rms and summed over the pixels containing the core and jet emission (thus avoiding bright artefacts). Subtracting this from the total flux density of 91 Jy, we arrive at an estimate of 2.6 ± 0.4 Jy for the diffuse flux density.

To estimate the diffuse emission the flux density of the entire jet was measured to a threshold of 10σ_rms and the sum of the knots flux densities subtracted, providing an estimate for the diffuse flux density of 2.8 ± 0.4 Jy. Some caution should be exercised with respect to the robustness of this value as it is possible that the diffuse emission is overestimated due to non-Gaussian shoulders present in the PSF. Conversely, due to observing the jet at a small angle to the line of sight, a significant fraction of the diffuse emission is likely to be in superposition with knots. This will act to both reduce the measured diffuse emission and increase the knot emission. Due to the relative brightness of the compact knots compared to the surrounding diffuse emission, this is unlikely to have a significant impact on the knot values discussed in Sect. 3.1, but acts here to reduce the amount of intrinsic diffuse emission we are able to measure. Determining the relative impact of each of these factors would require detailed modelling of the source and so is beyond the scope of this paper; however, we can make comparison to previous studies where such an estimate has been made.

Punsly & Kharb (2016) measured the diffuse flux density at 327 MHz to be 1.33 ± 0.13 Jy and extrapolated this to 151 MHz to arrive at an estimate of 2.7 Jy. Their image at 327 MHz had a resolution of 7 × 6 arcsec so the core, jet, and diffuse emission were not well-resolved. Punsly & Kharb (2016) determined that the diffuse flux density was underestimated by 50–100% and increased the value accordingly. In the LOFAR data, however, the core, jet, and diffuse emission are well resolved so we make no adjustments to the calculated flux density. Nevertheless, the range of jet powers still overlaps with the final estimates of Punsly & Kharb (2016) and so the more naive methods used in both papers appear to provide a reasonable estimate of the observed emission.

To convert optically thin 150 MHz flux densities from the lobes into estimates of the long term time-averaged jet power, Q, we apply the relations of Willott et al. (1999) and Punsly (2005) who find $$Q\approx {\left(\frac{f}{15}\right)}^{1.5}(1.1\times {10}^{45}){(X^{1+\alpha }{Z}^2{S}_{150})}^{0.857}$$(4) $$Z\equiv 3.31-3.65{(X^4-0.203X^3+0.749X^2+0.444X+0.205)}^{-0.125}$$ $$X\equiv 1+z,$$

where z is the redshift and f is an empirical multiplicative factor incorporating uncertainties associated with departures from minimum energy and variations in geometric effects, filling factors, hadronic contributions, and the low-frequency cutoff. Punsly & Kharb (2016) state that f is likely to be 10–20, and so we therefore setf = 20, treating thisas an upper limit. Punsly (2005) independently derive an estimator for the jet power, where the lobe energy is primarily inertial (i.e. thermal, turbulent, and kinetic energy). This assumption was validated by X-ray data that indicate the energy in radio lobes is dominated by inertial energy, as opposed to magnetic field energy, unlike the hotspots, which are usually near equipartition. Punsly (2005) show that $$Q\approx (5.7\times {10}^{44})\times (X^{1+\alpha }{Z}^2{S}_{150})$$(5)

represents a lower estimate on the jet power. Following the methodology of Punsly & Kharb (2016), we find that between the limits given by Eqs. (4) and (5) where f = 20, the jet power of 3C 273 is in the range 3.5 × 10⁴³–1.5 × 10⁴⁴ erg s⁻¹. This is in agreement with Punsly & Kharb (2016), who derive the jet power of 0.7–3.7 × 10⁴⁴ erg s⁻¹, which falls within the typical range for radio-loud quasars.

Fig. 9 LOFAR data overlaid on VLA data taken from Perley & Meisenheimer (2017) at P, L, and U band (top left, top right, and bottom, respectively). The beams of the VLA images are 4 arcsec (P band), 1.4 arcsec (L band), and 0.125 arcsec (U band), and are shown as dotted circles. Black contours indicate the VLA emission and are set such that Scont = σrms × 2n. The filled contours show the LOFAR dataset such that Scont = 40 × 2n mJy beam−1.

3.4 Presence of the counterjet

Despite the increased sensitivity and excellent uv coverage including short baselines provided by LOFAR observations, the counterjet remains undetected. Assuming intrinsically symmetrical jets where the observed asymmetry is due to relativistic beaming, from the 3σ_rms local to theexpected counterjet region, we determine an upper limit of S_cj = 40 mJy beam⁻¹.

Using the counterjet peak brightness limit along with the measured peak brightness of the jet discussed in Sect. 3.1 we are also able to provide constraints on the bulk speed and the angle to our line of sight. For a a pair of relativistic jets, where one jet is moving towards the observer (j) and one is moving away (cj), the ratio of the observed flux densities is given by $$\frac{S_{\text{j}}}{S_{\text{cj}}}={\left(\frac{1+\beta \cos \theta }{1-\beta \cos \theta }\right)}^{m-\alpha },$$(6)

where β ∈ [0, 1], θ ∈ [0, π∕2], m = 3 for a knotty jet, and α =0.8 for optically thin synchrotron emission. Substituting in our values of 4400 and 40 mJy beam⁻¹ for the jet and counterjet respectively, we find for 3C 273 that β cosθ ≥ 0.55.

Assuming an inclination angle of 5.5°, the midpoint of the range determined by Meyer et al. (2016; 3.8–7.2°), we are able to place a lower limit on the bulk speed of β ≳ 0.55. This lower limit is in line with studies by Meyer et al. (2016, 2017), where it was asserted that the bulk motion of the kiloparsec-scale jet is mildly relativistic at most. Using the standard equation for the bulk Lorentz factor $$\Gamma =\frac{1}{\sqrt{1-{\beta }^2}}$$(7)

for β ≥ 0.55 for 3C 273, we find that Γ ≥ 1.2. As discussed in Sect. 1, The IC/CMB mode requires significant bulk velocities on kiloparsec scales (Γ ≈ 10) that cannot be excluded by this lower limit, although a smaller angle to the line of sight is likely still required if such an emission mechanism were to dominate.

While we are able to provide some constraints on the jet properties using the upper limits given above, the absence of any counterjet emission raises some interesting questions. Given the steep spectrum of lobe and hotspot emission, LOFAR is arguably the instrument best placed to make such detection given its UV coverage and resolving power at low frequencies. The brightest regions of the jet (the H1–H2 lobe) should have an analogous counterpart if jet symmetry is assumed that should be the least impacted by beaming effects in the counterjet from the observers frame of reference. Itmay be the case that the counterjet emission remains below the detection limit of our observations, although with such a consistent lack of detection across all wavelengths, alternative interpretations should be considered.

One possible cause is that while the jets may be intrinsically similar, light travel effects may cause an observational asymmetry. Assuming the inclinations determined by Meyer et al. (2016) described above, the deprojected distance between the furthest extent of the two jets is between 479 and 905 kpc, which translates tolight travel time between ≈1.6 and 3.0 Myr. It is therefore possible that the counterjet is not observationally symmetric due to being observed at an earlier time and the location and/or morphology of the peak emission may differ from that of the main jet.

Another possibility that cannot be ruled out is that the counterjet is located behind the core, relative to our line of sight. This may be due to either a slight asymmetry in the jet launch angle or bending of the jet on kiloparsec scales. The narrow angle to the line of sight means that only a relatively small deviation would be required to cause the brightest outer lobe regions to be obscured by the core, and this provides a plausible explanation for the non-detection of the counterjet.

As is often the case, the true cause is likely some combination of beaming and light travel time effects, with possibly some jet asymmetry that is the cause of the persistent non-detection of the counterjet. It is beyond the scope of this paper to determine which, if any, of these scenarios is the dominant cause, but the lack of detection at the low frequencies suggests that interpretations other than a simple lack of sensitivity in the observations should be considered. Deeper observations using the LOFAR HBA (150 MHz) and at lower frequencies using the LOFAR LBA (~ 50 MHz) would proveinteresting in this regard by either detecting or further constraining the counterjet’s peak brightness.

Fig. 10 Spectral energy distributions of the knots in 3C273 with absorption models overlaid. The solid dark blue line indicates the FFA model and the dashed magenta line the SSA model. GHz data points taken from Jester et al. (2007).

4 Conclusions

In this paper, we have presented the first high-resolution observations of 3C 273 at MHz frequencies, recovering both compact and diffuse emission on sub-arcsecond scales. We have shown that robust high-fidelity imaging of low-declination complex sources is possible with the LOFAR international baselines, and that these observations are key to understanding the physics that underlie the small-scale structures of radio-loud AGN. By analysing the kiloparsec-scale morphology and spectrum of 3C 273 at MHz frequencies to answer the questions posed in Sect. 1 we find the following:

1. The small-scale structure of 3C 273 matches that at higher frequencies;

2. The low-frequency spectrum of the jet knots can be well described by free-free absorption (FFA) and synchrotron self-absorption (SSA) models;

3. We derive a kinetic power for the jet in the range 3.5 × 10⁴³–1.5 × 10⁴⁴ erg s⁻¹ which agrees with previous estimates derived from higher frequency observations;

4. The counterjet remains undetected at 150 MHz, placing limits on the peak brightness of S_{cj_150} < 40 mJy beam⁻¹;

5. We derive lower limits of β ≳ 0.55 and Γ ≥ 1.2 for the for the bulk speed and Lorentz factor respectively.

Further low-frequency investigations (e.g. using LOFAR LBA observations) are required to determine many of the details of 3C 273, such as counterjet detection and to better differentiate between absorption models. However, our investigation provides the first constraints on the characteristics of 3C 273 at low frequencies, and is the first step towards understanding the detailed mechanics of 3C 273 at these wavelengths, a step that will be vital moving forward with both the LOFAR international baselines and towards the SKA era.

Acknowledgements

We wish to thank the referee, Rick Perley, for his constructive comments and suggestions which have helped improve this paper. S.M. acknowledges support from the Irish Research Council’s Government of Ireland Postgraduate Scholarship scheme. LKM is grateful for support from the UKRI Future Leaders Fellowship (grant MR/T042842/1). C.G. acknowledges support from the ERC Starting Grant ClusterWeb 804208. J.M. acknowledges financial support from the State Agency for Research of the Spanish MCIU through the “Center of Excellence Severo Ochoa” award to the Instituto de Astrofísicade Andalucía (SEV-2017-0709) and from the grant RTI2018-096228-B-C31 (MICIU/FEDER, EU).

References

All Tables

All Figures

Fig. 1 VLA C and P band images of 3C 273. Left: VLA C band at 4885 MHz at a resolution is 0.5 arcsec. Contours are such that Scont = 2 × 2n mJy beam−1 with a peak brightness of 35.99 Jy beam−1. The inner jetis marked for clarity. Right: VLA P band at 327 MHz at a resolution of 0.5 arcsec. Contours are such that Scont = 80 × 2n mJy beam−1 with a peak brightness of 36.37 Jy beam−1. These data are reproduced from Perley & Meisenheimer (2017).

In the text

	Fig. 2 HST image of the 3C 273 outer jet masked at 0.7 e−∕ s to only show the jet. The direction of motion is left-to-right, with the core ~ 11 arcsec to the left of the image. In1, In2, and Ex1 are galaxies unrelated to the jet; fs is a foreground star. The naming convention is from Jester et al. (2005), with the addition of features fs, Ex1, and H1.
In the text

	Fig. 3 HST image of 3C 273 at 850–1700 nm with a peak brightness of 28606 e− s−1. The image combines four exposures providing a 37-s integration time.
In the text

	Fig. 4 Calibrator phase solutions for 3C295 for the XX polarisation. The reference station is CS001HBA0, and the plots increase by baseline length from the reference station reading from left to right. All plots are on the same time and frequency axes. The first remote station is labelled, as are individual international stations.
In the text

Fig. 5 Components to the 3C 273 spectral energy distribution. Data are from Perley & Meisenheimer (2017), except for the 151 MHz flux density of the halo, which is from Punsly & Kharb (2016). The flux density of the jet includes all knots. The halo refers to the diffuse cocoon of radio emission encompassing the jet (A, B, and R of Fig. 9). The radio spectral indices of the core, jet, and halo are 0.30 ± 0.07, − 0.81 ± 0.03, and − 0.79 ± 0.11, respectively. The spectral indices of the jet and the halo are consistent with optically thin synchrotron emission, as expected. The LOFAR HBA bandwidth is shown by the shaded region.

In the text

	Fig. 6 LOFAR images of 3C 273 at 151 MHz with a beam size of 0.31 × 0.21 arcsec. Contours levels are set such that Scont = 40 × 2n mJy beam−1. Shown are the full source (left) and the outer jet with the components labelled (right).
In the text

	Fig. 7 3C 273 ridgeline histogram with the brightness measured relative to the core. The image has been rotated counterclockwise by 42° such that the top of the image points north-west. The features around the core are likely to be imaging artefacts. The green dashed line shows the results of Davis et al. (1985), and the dotted black line the 3σ RMS noise of the LOFAR image.
In the text

	Fig. 8 LOFAR data overlaid on HST near-infrared data. Black contours indicate the HST emission, and are set such that Scont = 0.7 × 2n e− s−1. The filled contours show the LOFAR dataset such that Scont = 40 × 2n mJy beam−1.
In the text

Fig. 9 LOFAR data overlaid on VLA data taken from Perley & Meisenheimer (2017) at P, L, and U band (top left, top right, and bottom, respectively). The beams of the VLA images are 4 arcsec (P band), 1.4 arcsec (L band), and 0.125 arcsec (U band), and are shown as dotted circles. Black contours indicate the VLA emission and are set such that Scont = σrms × 2n. The filled contours show the LOFAR dataset such that Scont = 40 × 2n mJy beam−1.

In the text

	Fig. 10 Spectral energy distributions of the knots in 3C273 with absorption models overlaid. The solid dark blue line indicates the FFA model and the dashed magenta line the SSA model. GHz data points taken from Jester et al. (2007).
In the text