Discovery of a radio jet in the Cloverleaf quasar at z = 2.56

ABSTRACT

The fast growth of supermassive black holes and their feedback to the host galaxies play an important role in regulating the evolution of galaxies, especially in the early Universe. However, due to cosmological dimming and the limited angular resolution of most observations, it is difficult to resolve the feedback from the active galactic nuclei (AGNs) to their host galaxies. Gravitational lensing, for its magnification, provides a powerful tool to spatially differentiate emission originating from AGN and host galaxy at high-redshifts. Here, we report a discovery of a jet-like radio structure in a strongly lensed starburst quasar, H1413+117 or Cloverleaf at redshift z = 2.56, based on observational data at optical, submillimetre, and radio wavelengths. With both parametric and non-parametric lens models and with reconstructed images in the source plane, we find a well-separated, kpc-scaled, single-sided radio jet located at projected |${\sim }1.2\, \mathrm{kpc}$| to the north-west of the host galaxy in the source plane. This could indicate the co-existence of feedback from the AGN by both wind and jet in the Cloverleaf quasar.

gravitational lensing: strong, galaxies: jets, submillimetre: galaxies

1 INTRODUCTION

The correlated growth of the bulge mass of a galaxy and its supermassive black hole (SMBH) mass throughout cosmic time has been described as their co-evolution (Ferrarese & Ford 2005; Kormendy & Ho 2013). When the SMBH actively accretes, the active galactic nucleus (AGN) turns on, which can drive powerful feedback to the surrounding interstellar medium (ISM) of the host galaxy. The radiation pressure and the energetic particles in the AGN wind and jet behave as the main forms of feedback. The feedback process can sweep and heat up the cold gas along its radial paths and, as a result, quench the star formation (Springel, Di Matteo & Hernquist 2005; McNamara & Nulsen 2007; Cattaneo et al. 2009). For the high-luminosity quasars, the host galaxy and surrounding environment receive powerful feedback from its AGN through electromagnetic radiation (known as quasar mode), kinetic jets (known as radio mode), or both (Churazov et al. 2005; Fabian 2012; Kormendy & Ho 2013).

Similar to the cosmic star formation, the evolution of quasi-stellar object (QSO) activity also shows a strong cosmic evolution, reaching its peak at z ∼ 2 (e.g. Schmidt & Green 1983; Wolf et al. 2003). The feedback of QSOs in this era is therefore considered as the key mechanism to quench star formation. Galactic outflows driven by AGNs (e.g. Nesvadba et al. 2008, 2011; Prochaska & Hennawi 2009; Alexander et al. 2010; Cano-Díaz et al. 2012; Farrah et al. 2012) have shown evidence of quenching. High-velocity galactic winds, which are commonly detected with the broad absorption lines against a strong quasar continuum (e.g. Moe et al. 2009; Saez, Chartas & Brandt 2009; Dunn et al. 2010) and also evidenced with highly ionized Fe lines in low-luminosity AGNs (Shi et al. 2021), can efficiently remove the ISM from the host galaxies. On the other hand, radio lobes/jets (e.g. Gopal-Krishna & Wiita 2001; McNamara & Nulsen 2007; Cattaneo & Best 2009; Cattaneo et al. 2009; Fabian 2012), which are often discovered on scales from several tens of kpc to Mpc, can also directly impact star formation.

However, the majority of high-redshift QSOs have very compact sizes (Silverman et al. 2019; Li et al. 2021; Ding, Silverman & Onoue 2022), which limits the separation between their host galaxy and AGN-related components. Multiple wavelength observations with spatial resolutions on scales of sub-kpc are therefore needed, which are still extremely difficult, given the capability of current instruments and cosmological dimming effects (Lanzetta et al. 2002).

Gravitational lensing serves as a powerful tool to improve both sensitivity and spatial resolution. With proper demagnification modelling, the lensing system allows us to resolve the detailed galactic structures of high-z QSOs for case studies. Cloverleaf (also known as H1413+117) is one such target. Discovered by Magain et al. (1988), it is a strongly lensed starburst galaxy with a bright AGN at z = 2.558. It contains four images at both optical and submillimetre wavelengths (Kneib et al. 1998). With the detections of CO, HCN, and far-infrared (FIR) emission, Cloverleaf is confirmed to contain an extensive molecular disc of ∼1.6 kpc in diameter, a molecular gas mass of |$\sim 10^{10}\, \rm {\rm M}_{\odot }$|⁠, and a star formation rate of |$10^3\, \rm {\rm M}_{\odot }\, yr ^{-1}$| (Barvainis et al. 1997; Solomon et al. 2003; Venturini & Solomon 2003; Weiß et al. 2003). On the other hand, information about this system in the lens plane is limited due to the lens galaxy’s faint emission (Angonin et al. 1990; Lawrence 1996; Turnshek et al. 1997), while Chantry & Magain (2007) suggested a deconvolved lens image from Hubble Space Telescope (HST) F160W and F160W data. The redshift of the lens galaxy is derived from the H i narrow absorption systems (Magain et al. 1988; Turnshek et al. 1988) at z ∼ 1.7.

With the latest high-sensitivity, high-resolution, and high-fidelity data of Cloverleaf system from ALMA, Very Large Array (VLA), and HST, we are able to not only compare morphologies from different gas phases but also constrain the lens galaxy mass models more accurately than previous models. The traditional method (e.g. Kneib et al. 1998), which is based on parametric modelling, has limitations in the cases of irregular sources. The newly developed non-parametric modelling, which is required for high angular resolution and high-sensitivity analysis, greatly improves the reconstructed source plane (Nightingale & Dye 2015; Nightingale, Dye & Massey 2018) and thus, the structures of the Cloverleaf system can be well distinguished.

This work focuses on observational data at three wavelengths: (i) optical, which is dominated by the QSO emission and is seen as almost four point sources from HST images (Kneib et al. 1998); (ii) submillimetre, which traces the gaseous dusty disc of the host galaxy (Granato, Danese & Franceschini 1996; Solomon et al. 2003), with current ALMA observations having enough resolution to resolve the four images and; (iii) radio, where the radio emission traces relativistic charged particles accelerated by the radio jet (Kayser et al. 1990) launched by the central AGN or starburst. With these data, we can separate the AGN contribution from the host galaxy.

This paper is organized as follows: Observations and data reduction are presented in Section 2. The modelling processes, the lens models of the host galaxy, and the newly discovered radio jet, which are built up with both parametric and non-parametric models, are shown in Section 3. In Section 4, we analyse the spatial distribution of multiwavelength images and the spectral energy distribution (SED) of radio emission. In Section 5, we discuss the confirmation, emission components, and feedback of this radio jet. Last, we conclude the summary in Section 6.

Throughout this paper, we assume a cosmology with Ω_m = 0.310, Ω_Λ = 0.689, and |$H_0 = 67.7\, \mathrm{km\, s^{-1}\, Mpc^{-1}}$| (Planck Collaboration 2020). In this case, 1 arcsec corresponds to ∼8.7 kpc in the lens plane at z = 1.7 and ∼8.22 kpc in the source plane at z = 2.56.

2 MULTIBAND OBSERVATIONAL DATA AND DATA REDUCTION

We obtain multiwavelength, high-angular resolution data from the archival systems of HST (optical), ALMA (submillimetre), VLA (radio), and e-Merlin (radio).

2.1 HST F814W Image

The optical image is obtained from the HST archive. This three-orbit observation was performed on 2000 January 26, with the Advanced Camera for Surveys (ACS) and the F814W filter. We choose the ACS images for their high spatial sampling rate of 0.05 arcsec pix⁻¹ (Lucas & Rose 2021). We stacked all ACS images with swarp (Bertin et al. 2002). We use imfit in Common Astronomy Software Applications version 6.1.2.7 (casa; McMullin et al. 2007) to fit a star located at [14: 15: 47.39856, +11: 29: 50.7480; 2MASS id: J14154739+1129507 (Cutri et al. 2003)] in the field and obtain its full width at half-maximum (FWHM) to be ∼0.08 arcsec, which is adopted as the image resolution. The rms noise is evaluated to 0.2 |$\mathrm{erg\, s^{-1}\, cm^{-2}\, {\mathring{\rm A}}^{-1}}$| (below referred as counts).

2.2 ALMA 300 GHz continuum

We adopt the ALMA archival data from projects 2012.1.00175.S (PI: van der Werf, Paul) and 2017.1.00963.S (PI: Sharon, Chelsea), both observed at Band 7 (see Table 1 for observational details).

Table 1.

Details of the observations on Cloverleaf from the ALMA data archive.

Project ID	Date	Calibrators		Frequency coverage		Baseline	On-source time
		Gain	Flux/BP	LSB (GHz)	USB (GHz)	kλ	min
2012.1.00175.S	2015 June 25	Titan	J1415+1320	325.86–329.28	338.43–341.46	28–1433	68
2012.1.00175.S	2015 June 30	J1550+054	J1415+1320	276.68–279.62	289.50–292.43	36–1530	9
2012.1.00175.S	2015 September 27	J1550+054	J1415+1320	311.93–314.43	325.86–327.85	40–2414	17
2017.1.00963.S	2018 May 22	J1337+1257	J1347+1217	341.42–345.04	353.42–357.04	14–372	5

We calibrate the raw data with casa (version 6.1.2.7, McMullin et al. 2007) with the associated standard pipeline. Then, we combine all the calibrated on-source data and omit the frequency range of CO J = 11–10. We invert the combined visibility data by task tclean with the ‘mfs’ mode, a gridder of mosaic, a cell size of 0.02 arcsec, and a Briggs weighting with a robust of 1.5. The synthesized image (see Fig. 1 top right) has a beam size of 0.26 arcsec × 0.20 arcsec and a position angle of −31°. The noise level is |$35\, \rm {\mu Jy\, beam^{-1}}$|⁠, estimated from the emission-free area in the field.

Figure 1.

Images of Cloverleaf at optical, submillimetre, and radio bands. Top left: HST optical data where components A, B, C, and D are four images of this lensing system; the black pluses and the white pluses mark the peak positions of the model images with the parametric and the non-parametric lens model respectively. Top right: dust continuum data around 300 GHz where the black crosses mark the peak positions of four components in the HST data. Bottom: Radio images of Cloverleaf at 1.5 (left), 8.4 (middle), and 33 GHz (right), respectively. Black crosses show the peak positions of four components in the HST data. White lines are the contours of the ALMA dust continuum, starting from 25σ, with steps of 30σ. The black ellipse presents the region we extract fluxes of the radio jet. We use the same restoring beams of 0.25 arcsec in all radio images to generate the same source planes for spectral index analysis.

2.3 VLA archival data

2.3.1 VLA 8.4 GHz continuum

The VLA 8.4 GHz data were obtained from the National Radio Astronomy Observatory (NRAO) archive.¹ This data have been originally presented by Kayser et al. (1990). The project IDs are AS357 and AC243, which were performed on 1989 January 13 and 1989 February 4, respectively. Both observations were done with the A-configuration of the VLA. The total observing time was ∼5 and 4 h, with ∼1.5 and 1 h effective on-source time, respectively. Each on-source scan took 15 and 10 min with a sampling time of 1 s. In both observations, 3C 286 (1328+307) was adopted as both flux and bandpass calibrators, while the phase calibrator was 1413+135 in B1950 (or J1415+1320 in J2000). The receivers were configured with two spectral line windows (SPW), centring at 8.415 and 8.464 GHz, respectively. Each SPW covers a 50 MHz bandwidth, offering a total bandwidth of 100 MHz, with full stock polarizations. The baselines of both observations range from 15 to 950 kλ in the ultraviolet (UV) space.

The data was first imported to Astronomical Image Processing System (AIPS) and converted to the uvfits format. We use casa (ver. 6.1.0) to calibrate the data manually following standard procedures. We then image the data with task tclean in casa, with the ‘mfs’ mode. We use a Briggs weighting with a robust of 0.5. The final synthesized beam is 0.25 arcsec × 0.23 arcsec with a position angle of −47°. We then convert the equatorial coordinate from the B1950 epoch to the J2000 epoch with task imregrid. The rms noise level is |${\sim }14\, \rm {\mu Jy\, beam^{-1}}$|⁠.

2.3.2 VLA 33 GHz continuum

We obtained the Ka band Jansky Very Large Array (JVLA) data from the National Radio Astronomy Observatory (NRAO) archival system. The project number is 13B-051 (PI: Sharon, Chelsea). This project consists of a total of 18 execution blocks (EBs) targeting Cloverleaf. However, three of them, taken before 2013 October 22, were corrupted and not usable due to incomplete execution. The remaining 15 EBs were taken from 2013 October 27 to 2014 January 5, with the B-array configuration. The pointing observations were performed at the X band, and the observations were performed at the Ka band. These observations adopted 27 antennas and took a total of ∼6.8 h on-source time.

The receiver covers a frequency range of 31.95–34.00 GHz, with 16 spectral line windows (128 MHz each) and two polarizations. The total bandwidth is ∼2 GHz. The bandpass/flux calibrator and the gain calibrator are 3C 286 and J1415+1320, respectively. The baselines range from 16 to 1192 kλ.

We calibrate the raw data using casa (ver. 5.6.2–2) with the standard VLA data reduction pipeline. Then we combine all calibrated data and abandon the spectral windows that contain the CO J = 1→0 line emission. To save computational resources, we bin all channels for each spectral line window. Considering a maximum offset of ∼2 arcsec from the phase centre, the bandwidth smearing effect is negligible ((Δν/ν₀) × (θ_offset/θ_HPBW) ≤ 0.026) for Cloverleaf. Then we image the visibility data with casa (ver. 6.1.2.7), using task tclean. We adopted the ‘mfs’ mode, a gridder of standard, a cell size of 0.05 arcsec, and a Briggs weighting with a robust of 1.5. The cleaned image has a beam size of 0.26 arcsec × 0.23 arcsec with a position angle of −37°. The rms noise level is |${\sim }4.6\, \rm {\mu Jy\, beam^{-1}}$|⁠.

2.4 e-Merlin 1.5 GHz continuum data

The e-Merlin 1.5 GHz data was obtained from the e-Merlin archival system. The project number is CY4215. The observation was performed on 2017 February 9, with a total observing time of 13 h. The receiver covers a frequency range of 1.25–1.77 GHz, with eight spectral windows, of which each has 64 MHz sampled with 128 channels. The total bandwidth is 512 MHz. The average sampling time is 2.0 s.

In total, seven antennas were adopted during the observations. The bandpass and flux calibrator was 1331+305 (3C286), and the gain calibrator was J1415+1320. Every 10 min, the telescope switch between the target (7 min) and the gain calibrator (3 min). The observations spanned ∼7.7 h on source in total. The baselines range from 24 to 1276 kλ in the UV space.

The archival data has been correlated but not calibrated. We first convert the common uvfits format to the format of ‘measurement set’ and use casa (ver. 6.1.0) to flag the data and calibrate manually. For imaging, we adopted the ‘mfs’ mode, a gridder of standard, a cell size of 0.05 arcsec, and a Briggs weighting scheme with a robust of 1.5. The cleaned image has a default beam size of 0.26 arcsec × 0.18 arcsec with a position angle of 25⁰. The rms noise level is |${\sim }47\, \rm {\mu Jy\, beam^{-1}}$|⁠.

We also clean the e-Merlin data with a nterm of 2 to get the spectral index map in the lens plane (see Fig. 4).

In the spectral slope index analysis, we use a uniform restoring beam of 0.25 arcsec × 0.25 arcsec for all radio images, to analyse these data with the same reconstructed source plane.

2.5 Astrometric correction

We notice a systematic astrometric position offset between the HST image and those at radio/sub-mm bands. We fit the peak positions of the four images with SAOIMAGEDS9 (Joye & Mandel 2003) for both the HST data and the ALMA continuum data. From these, we get an average offset of 0.126 arcsec in the right ascension (RA) and 0.219 arcsec in the declination (Dec.). After we corrected the offset, the peak positions of the images obtained from HST and ALMA (see Fig. 1) are still not yet perfectly aligned (though the current offsets are within one pixel of the HST image), which may originate from different emission areas and/or imperfections of the observations.

The radio band observations use the same phase calibrator J1415+1320. Their coordinates are almost identical to each other. The largest offset is less than 13 mas, which happens between the VLA 8.4 and 33 GHz data. Since the offset is smaller than the pixel scale, we neglect it in the lens modelling below.

3 LENS MODELLING

To analyse the detailed morphology of multiwavelength images in the source plane, we adopt the lens model software PyAutoLens (AutoLens²; Nightingale & Dye 2015; Nightingale et al. 2018, 2021b). We model Cloverleaf with both parametric and non-parametric models for the lensed source galaxy, following the visual step-by-step guide provided in PyAutoLens³

We take the assumption from Kneib et al. (1998) that the lens system is contributed by a single galaxy (the narrow absorption systems) and a distant cluster both at z ∼ 1.7.

3.1 Overview of lens modelling strategy

We perform lens modelling with the following series of linked heuristic steps:

• Fit the HST data using a simple lens model which assumes point source emission for the source galaxy. This provides an initial estimate for the parametric lens mass model.

• Take the highest likelihood results from the previous step as the initial input to fit the ALMA 300 GHz continuum data with another parametric mass model, using an extended source model which follows a smoothly parametrized form of elliptical Sérsic profile.

• Using the mass model inferred in the previous step to initialize the fit, fit different waveband data sets with a non-parametric source reconstruction that uses a Voronoi mesh.

In our lens model, the redshift of the lens plane is set as 1.7 (Kneib et al. 1998).

3.2 Paramtric fitting – optical and sub-mm images

3.2.1 The HST optical data

The HST optical (814 nm) image is shown in the top left-hand panel of Fig. 1. We use imfit to measure the sizes (FWHM) of the four components (see Table 2), by assuming that they follow Gaussian distributions. The fitted FWHMs of all four components are almost identical to the FWHM of a point source (∼0.08 arcsec), so all these images can be treated as compact point sources without extended structures, which is also consistent with the result that the optical emission is dominated by a point like AGN in the literature (Magain et al. 1988).

Table 2.

Apparent sizes of four images of Cloverleaf at different wavelengths, fitted with Gaussian profiles. The sizes are described as (the FWHM of the major axis) × (the FWHM of the minor axis).

Band	ID	Size (⁠|$\rm mas\, \times \, mas$|⁠)	Peak S/N
814 nm	A	89.4 × 85.6	173
814 nm	B	84.3 × 77.5	151
814 nm	C	91.3 × 86.7	106
814 nm	D	82.4 × 76.6	106
300 GHz	A	540 × 183	116
300 GHz	B	490 × 254	122
300 GHz	C	478 × 216	112
300 GHz	D	426 × 180	75.8
1.5 GHz	A	540 × 257	17.4
1.5 GHz	B&D	443 × 121	28.4
1.5 GHz	C	262 × 161	10.8

First, to model the HST data, we adopt an isothermal ellipsoid with an external shear for the lens galaxy’s mass and use a ‘point’ to represent the source galaxy. With this optical-band lens modelling, we can preliminary estimate the lens mass model, as well as the source’s flux and position.

The four optical images of Cloverleaf point sources may suffer from microlensing (Kayser et al. 1990), which may lead to a biased result since the microlensing effect is not included in the strong lens modelling. Additionally, the limited number of data (i.e. flux and position) provided by the four lensed images only barely constrain the lens model, resulting in considerable uncertainties. These two limitations beg for a refined analysis at the submillimetre band, on which we perform the lens modelling with an extended lensed arc and use the modelling result of the optical band as a prior. The extended lensed arc at the submillimetre wavelength is not subject to the microlensing effect and its entire brightness distribution can be used to improve the precision of the lens model.

3.2.2 The ALMA 300 GHz continuum

We next fit the ALMA 300 GHz submillimetre data with a more complex parametric lens model. We adopt a power-law ellipsoid (Tessore & Metcalf 2015) with an external shear as the lens galaxy mass profile and model the source with an elliptical Sérsic profile.

We fix the centre of the lens mass model (from the result of the HST data modelling) to reduce the dimensions of parametric space. All remaining parameters from the previous model are then used to create the initial input for this modelling. The inferred parameters of the best-fitting model at 1σ confidence are shown in Table 3, including the centre positions (x and y), the elliptical components (e_x and e_y), the Einstein radius (r), the slope of the power-law profile (α), and the elliptical components of the external shear (γ_x and γ_y).

Table 3.

Best-fitting parameters with 1σ confidence of the lens object models derived from the ALMA 300 GHz continuum data.

Parameters	Definition	Parametric model	Non-parametric model
x	RA (hh:mm:ss) at mass centre	14:15:46.2384 |$^{+\ \lt 0.0001}_{-\ \lt 0.0001}$|	14:15:46.2325 |$^{+\ \lt 0.0001}_{-0.0002}$|
y	Dec. (° ′ ″) at mass centre	11:29:43.705|$^{+\ \lt 0.001}_{-\ \lt 0.001}$|	11:29:43.694|$^{+\ \lt 0.001}_{-0.002}$|
ex	Elliptical component x	|$-0.136^{+0.002}_{-0.002}$|	|$-0.126^{+\ \lt 0.001}_{-0.001}$|
ey	Elliptical component y	|$-0.263^{+0.002}_{-0.002}$|	|$-0.270^{+0.003}_{-\ \lt 0.001}$|
r	Einstein radius (arcsec)	0.645|$^{+0.001}_{-0.001}$|	0.636|$^{+\ \lt 0.001}_{-0.001}$|
α	Slope	1.897|$^{+0.003}_{-0.003}$|	1.884|$^{+\ \lt 0.001}_{-0.014}$|
γx	Elliptical components x of shear	0.034|$^{+0.001}_{-0.001}$|	0.043|$^{+\ \lt 0.001}_{-0.004}$|
γy	Elliptical components y of shear	|$-0.054^{+0.001}_{-0.001}$|	|$-0.057^{+\ \lt 0.001}_{-0.003}$|

3.3 Non-parametric fitting

Non-parametric source reconstructions with PyAutoLens are performed using a Voronoi mesh, where the Voronoi cells adapt to the source morphology. In this way, a relatively higher resolution is dedicated to its bright central regions than that in the outer region with weaker signals. The source reconstruction is performed after an assumed parametric mass model maps coordinates from the lens plane to the source plane with ray-tracing.

3.3.1 The ALMA 300 GHz continuum

When modelling the 300 GHz submillimetre data, we adopt iterative steps. To fit the non-parametric source, we first fix the parameters of the lens object according to our parametric mass model (see Section 3.2.2), and then use a pixelized Voronoi grid that adapts to the mass model’s magnification pattern, giving higher resolutions to areas with higher magnification.

Then, we perform a fit where the source plane is reconstructed by a brightness-weighted Voronoi (Fig. 2 left) mesh. We fit a new mass model, a power-law ellipsoid (Tessore & Metcalf 2015) with an external shear (shown in Table 3). With the Bayesian inference (dynesty) integrated within PyAutoLens, we adopt the best-fitting results and their associated errors as the lens profile (Table 3), which is further used for the non-parametric modelling of the radio data.

Figure 2.

Reconstructed source planes of Cloverleaf at multiple bands: Top left: submillimetre continuum at 300 GHz; top right: radio continuum at 1.5 GHz; bottom left: radio continuum at 8.4 GHz; bottom right: radio continuum at 33 GHz. The solid black line is the inner caustics. The centre is on (14:15:46.2384^h, 11:29:43.705) in the J2000.0 equatorial system of coordinates. The white solid ellipse and the red dashed ellipse are the areas to draw fluxes for SEDs of the radio jet and the host galaxy, respectively.

3.3.2 The 1.5, 8.4, and 33 GHz continuum

The ALMA 300 GHz continuum image provides strong constraints on our lensing galaxy mass model, therefore we fix the mass model parameters to values derived from this data in order to perform non-parametric fitting at radio bands. Using the same mass model across all wavelengths also ensures source reconstructions can be compared in a uniform way.

For the data from VLA and e-Merlin, the source plane is reconstructed with a magnification-weighted Voronoi mesh instead of a brightness-weighted one. This is because a magnification based mesh produces an identical grid of Voronoi cells across all three frequency channels, ensuring consistent sampling of the source when calculating the spectral indices pixel-by-pixel to obtain a spectral indices map. The residual maps at multiple wavelengths are presented in Appendix  B. We measure the radio fluxes within two regions in the source plane, shown as white and red ellipses in Fig. 2. The error compromises a 15 per cent absolute flux uncertainty, a |${\sim }\, 5~{{\ \rm per\ cent}}$| error from the uncertainty of the model by running multiple models with all parameters randomly sampled with their probability density distributions, and a statistical error measured in emission-free regions. The measured quantities of different components in the reconstructed source plane are listed in Table 4.

Table 4.

Measured quantities of the Cloverleaf components in the source plane at different wavelengths, where the fluxes are obtained within the ellipses in Fig. 2, and the intrinsic sizes are described in (the FWHM of the major axis) × (the FWHM of the minor axis). The error (1σ) consists of contributions from statistical error and model uncertainty and from the absolute calibration error (15 per cent, shown in the squared bracket).

Component	Frequency	Flux	Intrinsic size	Peak position
	(GHz)	(µJy)	(arcsec × arcsec)	(RA, Dec.)
Host galaxy	300	2.0 ± 0.3 [± 0.3] × 103	0.12 × 0.11	(14:15:46.2378, 11:29:43.659)
Host galaxy	33	4.8 ± 1.3 [± 0.7]	0.14 × 0.11	(14:15:46.2364, 11:29:43.659)
Host galaxy	8.4	7 ± 2 [± 1]	0.14 × 0.095	(14:15:46.2383, 11:29:43.698)
Host galaxy	1.5	26 ± 10 [± 4]	–	(14:15:46.2388, 11:29:43.689)
Radio jet	33	4.7 ± 1.0 [± 0.7]	0.15 × 0.086	(14:15:46.2315, 11:29:43.787)
Radio jet	8.4	14 ± 3 [± 2]	0.11 × 0.049	(14:15:46.2315, 11:29:43.781)
Radio jet	1.5	1.1 ± 0.2 [± 0.2] × 102	0.12 × 0.045	(14:15:46.2299, 11:29:43.781)

To verify the consistency between the optical data and our parametric/non-parametric models, we invert a simple point source to the lens plane, with lens mass models derived from parametric and non-parametric fitting. Both positions and fluxes of the point source are adopted from the result first fitted with the HST data. For both parametric and non-parametric models, the positions of the simulated point source (Fig. 1) match well with all four images in the HST data, within 0.5 pixel size. This indicates that our models are reliable.

4 RESULTS

4.1 Spatial distributions of multiwavelength images

As shown in Section 3, the HST optical images are point sources dominated by the central AGN. On the other hand, the 300 GHz image (see Fig. 1 top right) shows that all four components are much larger than the beam size (see Table 2), indicating that the sub-mm continuum emission has extended structures.

The radio continuum images (shown in Fig. 1) also show resolved extended structures. The B and D components tend to merge when moving to the lower frequency, indicating that the source is lying on the caustics. The other two, components A and C, both also have spatial deviations, not only between different radio bands but also between the images from HST and ALMA.

At 1.5 GHz band, the radio-to-sub-mm spatial offsets are (0.24, −0.01 arcsec) and (−0.03, 0.24 arcsec) in RA and Dec., for components A and C, respectively. All components in the images of 8.4 and 33 GHz are irregular with respect to the synthesized beam shapes at the corresponding wavelengths; thus, we do not fit their locations with imfit.

Among the three radio bands (1.5, 8.4, and 33 GHz), the image morphologies are not identical in the source plane. The 1.5 GHz image shows almost no emission from the host galaxy, while both 8.4 and 33 GHz radio images show detections at the same position. This suggests the multiple energy sources among different wavelengths.

4.2 The discovery of a radio jet

As shown in the reconstructed image of the source plane at 300 GHz of Cloverleaf (Fig. 2), the cold dust emission shows an extended disc-like structure of 0.115 arcsec × 0.109 arcsec in FWHM (⁠|$0.946\, \mathrm{kpc}\times 0.897\, \mathrm{kpc}$|⁠) rather than a point source.

The peak of this dust continuum has a signal-to-noise ratio (S/N) ∼27, with substructures around it. The majority of the host galaxy emission is distributed in the central region, inside the inner caustics, while a small clump is located in the south-west to the centre. Due to the uncertainty of the lens model, several emission clumps outside the caustics have less fidelity.

All reconstructed radio images show extended emission to the north-west of the dust continuum core, at a projected distance of ∼1.15 kpc (0.14 arcsec) away from the centre of the host galaxy. The peaks of this radio component are roughly consistent among different radio wavebands. This indicates a single-sided radio jet launched from the host galaxy, with a comparable size to the host galaxy at 300 GHz. The FWHMs of this radio clump, fitted with Gaussian profiles, are |$0.99\, \mathrm{kpc}\times 0.37\, \mathrm{kpc}$|⁠, |$0.90\, \mathrm{kpc}\times 0.40\, \mathrm{kpc}$|⁠, and |$1.23\, \mathrm{kpc}\times 0.71\, \mathrm{kpc}$|⁠, at 1.5, 8.4, and 33 GHz, respectively.

On the other hand, in the source plane, the 1.5 GHz radio emission from the host galaxy is almost negligible, while the 8.4 and 33 GHz images show radio emission at the location of the host galaxy, with their S/Ns of 3.5 and 3.7, respectively.

4.3 Spectral indices of the radio continuum

The multiple-band radio maps allow us to probe the SEDs of both the host galaxy and the radio jet. To evaluate the impact of the differential magnification effect and to probe the power sources of the radio emission, we construct SEDs derived from both the lens plane and the source plane (see Fig. 3).

Figure 3.

Spectral energy distributions (SEDs) of the radio jet and the host galaxy of Cloverleaf. Top left: SED of the radio jet in the lens plane. Top right: SED of the radio jet in the source plane. Bottom left: SED of the host galaxy in the lens plane. Bottom right: SED of the host galaxy in the source plane. Radio fluxes are measured from the regions shown in the lens plane (Fig. 1) and the source plane (Fig. 2). Orange lines show the best linear fitting result for all three frequencies. The black dashed line in the top left-hand panel shows the spectral index of the radio jet in the lens plane derived from the 1.5 GHz e-Merlin data only. All derived spectral indices are listed in Table 5. The error bars are at 1σ confidence consisting of statistical error, lens model uncertainty, and absolute calibration error (15 per cent).

We first measure radio fluxes from regions shown in the lens plane (Fig. 1) and the source plane (Fig. 2) for 1.5, 8.4, and 33 GHz, respectively.

The regions adopted to measure fluxes in the source plane are shown in Fig. 2. In the lens plane, however, the emission from the host galaxy and the radio jet do not show distinct spatial offsets between each other. The radio emission of components B and D tend to merge together, indicating that this emission component is lying on the caustics (Appendix  A), where the radio jet populates while the host galaxy does not. Therefore, to eliminate contamination from the host galaxy, we choose this region (shown as black ellipses in Fig. 1) in the lens plane and construct the SED of the lensed radio structure. On the other hand, we adopt beam-sized regions at the peaks of the four image components to measure fluxes from the host galaxy. We then sum up fluxes from all four components to increase the S/N.

Then we fit spectral indices with fluxes from all three radio bands, which are shown in orange lines of Fig. 3. All SEDs can be fitted with one single power-law profile.

We also use the e-Merlin data only, which covers 1.25–1.77 GHz bandwidth, to fit a spectral index map (Fig. 4). We use a 5σ mask from the e-Merlin continuum to mask the spectral index map and the error map. The average spectral index of the merging components B and D, which are labelled with the black dashed ellipse shown in Fig. 4, same as the black ellipse in Fig. 1, to be |$\alpha ^{1.25}_{1.77}=1.35 \pm 0.18$|⁠. It represents the average spectral index of the lensed radio jet at 1.5 GHz in the lens plane. All derived spectral indices for both the host galaxy and the radio jet, are listed in Table 5.

Figure 4.

Spectral index map (left) and its error map (right) in the lens plane at 1.5 GHz. The white and red lines are the contours of the images at 1.5 GHz (8, 14, 20σ) and 300 GHz (50σ), respectively. The red circles represent the beam sizes. The mask we use is a 5σ cutoff from the e-Merlin continuum. The black dashed line is the area we derive the spectral index at 1.5 GHz.

Table 5.

Spectral indices of different components in the lens plane and the source plane within 1σ: α_fitted is the spectral index of the best linear fit, considering fluxes in all three radio bands.

Component	Plane	|$\alpha ^{1.5}_{8.4}$|	|$\alpha ^{8.4}_{33}$|	αfitted	|$\alpha ^{1.25}_{1.77}$|
Radio jet	lens	1.10 ± 0.12	1.12 ± 0.16	1.11 ± 0.01	1.35 ± 0.18
Radio jet	source	1.17 ± 0.16	1.04 ± 0.19	1.11 ± 0.04	–
Host galaxy	lens	0.84 ± 0.12	0.76 ± 0.16	0.80 ± 0.02	–
Host galaxy	source	0.79 ± 0.15	0.58 ± 0.19	0.70 ± 0.06	–

Although the average slope of the e-Merlin spectral index map is consistent with that obtained with |$\alpha ^{1.5}_{8.4}$| within 1σ in the whole region, it also shows a systematic gradient of the radio slope from the east to the west direction. This indicates that the radio jet may have complex structures or physical processes that need higher-resolution data to study in detail.

From the flux maps shown in Fig. 2, we fit a spectral index map (Fig. 5, by the opposite sign convention of α) using delensed flux maps in the source plane. As shown in Fig. 5, we fit power-law profiles for all pixels at maps of 8.4, 1.5, and 33 GHz, with a mask of S/N > 3.5 at the 33 GHz source plane image. The spectral index map shows flatter SED at the host galaxy position, while steeper SED at the radio jet.

Figure 5.

Spectral index map (derived from the reconstructed images at 1.5, 8.4, and 33 GHz) in the source plane: The black solid line is the inner caustics. The centre is on (14:15:46.2384^h, 11:29:43.705) in the J2000 equatorial system of coordinates. The white solid ellipse and the red dashed ellipse mark the positions of the radio jet and the host galaxy, respectively.

5 DISCUSSION

5.1 Position offsets between submillimetre and radio emission

Compared with the lensed image of Cloverleaf at ∼300 GHz, those at radio bands have different spatial distributions, which indicates that the emission distributions at different bands are different from one another in the source plane, similar to the case of strongly lensed quasar MG J0751+2716 with a radio jet (Spingola et al. 2018; Powell et al. 2022). This can be inferred directly from the reconstructed source planes (see Fig. 2). All reconstructed radio images contain radio structures on scales of |${\sim }0.9\times 0.4\, \mathrm{kpc}$|⁠, on the north-west part of the inner caustics, which are absent at the submillimetre wavelength.

5.2 Possibility to be a jet in the foreground?

One may also ask if the radio feature comes from the foreground lens galaxy, instead of the background source galaxy. First, the multiple components shown in all radio bands indicate that it is unlikely to be multiple jets from the centre. Second, the spatial distributions of all radio components are highly consistent with the scenario that the radio structure locates on the caustics, i.e. merging between lensed components (Appendix  A). Therefore, we suggest that this radio structure should reside in the source plane and originate from the position on the north-west part of caustics.

5.3 Why one-sided?

On the other hand, the radio source only presents on one side of the host galaxy, while the other side shows almost no signal. Depending on the geometry, the other side of the radio jet might be too faint to be detected due to being away from the caustics, or simply does not exist. In any case, if the other side of the radio jet exists, it should lead to emission populating the area between components A and C (see Appendix  A). Future deeper observations would be the key to verifying this scenario.

5.4 The power-law profiles of the SEDs

As shown in Fig. 3 and Table 5, the SEDs of the radio jet in both lens plane and source plane exhibit power-law profiles, and have similar spectral indices, which indicates that the emission is dominated by the synchrotron from electrons in the jet material.

On the position of the host galaxy, the SED can be also fitted with one single power-law profile with a flatter spectral index (see Fig. 5 and Table 5) compared with the one at the position of the radio jet. The emission (rest frequency, f_rest ∼5–120 GHz) at this position may have a contribution from free–free and/or dust continuum (Condon 1992), together with a contribution from the synchrotron emitted from the relativistic electron accelerated by supernovae remnants.

5.5 AGN feedback of Cloverleaf

Based on the broad absorption line profile (Hazard et al. 1984; Magain et al. 1988) and the HST optical images (Fig. 1) of this quasar, the wind/radiation feedback could take place in this galaxy.

On the other hand, this radio jet/structure is located at a projected distance of ∼1.15 kpc (0.14 arcsec) from the galaxy centre, with a size similar to that of the host galaxy at 33 GHz. These are consistent with radio lobe features found in massive radio galaxies (Fabian 2012), however, Cloverleaf seems to have a lower galaxy mass than those massive galaxies (Solomon et al. 2003). The kinetic energy residing in the jet could severely heat up the gas and drive strong feedback to the star formation in Cloverleaf.

Chartas et al. (2004) constrained the mass of the SMBH in the Cloverleaf quasar by the microlensed image A (the one located at the south in the quadruple images shown in Fig. 1) at X-ray band. They also adopted the luminosity from Granato et al. (1996)’s dust-enshrouded AGN models. Then the Eddington ratio of this AGN is found to be ∼0.1, assuming the magnification of the AGN to be 11 (Venturini & Solomon 2003). The Eddington ratio of 0.1 indicates that the AGN is not radiating close to the Eddington limit, yet is still a powerful quasar.

Therefore, given its strong ongoing starburst, likely driven by the interaction with a neighbouring gas rich galaxy (Stacey & Arrigoni Battaia 2022), together with the broad absorption line and the newly discovered radio jet, we suggest that the AGN feedback in Cloverleaf consists of coexistent contributions from both radiation-driven wind and AGN-launched jet. High angular resolution observations of both radio continuum and molecular gas tracers in the future, which help resolve the jet structure and gas kinematics, could throw light on the co-existence of both mechanisms of AGN feedback.

6 CONCLUSION

With high angular resolution photometric data at optical, sub-mm, and radio wavelengths, we reconstructed multiband source images of the gravitational lensed Cloverleaf QSO, with both parametric and non-parametric modelling. From different emission locations, we discovered a single-sided radio jet, with a size of ∼1 kpc, located at projected ∼1.2 kpc to the north-west of Cloverleaf. The co-existence of a radio jet and broad absorption line indicates that the Cloverleaf quasar likely feedbacks its host galaxy with both AGN wind and radio jet. This provides a unique case for the galaxy-SMBH co-evolution studies at high-redshifts.

ACKNOWLEDGEMENTS

Based on observations made with the NASA/ESA Hubble Space Telescope, and obtained from the Hubble Legacy Archive, which is a collaboration between the Space Telescope Science Institute (STScI/NASA), the Space Telescope European Coordinating Facility (ST-ECF/ESAC/ESA) and the Canadian Astronomy Data Centre (CADC/NRC/CSA). LZ acknowledges the referee for pointing out the mistake in data reduction and for kind and helpful suggestions, which heavily polish this study. LZ appreciates Yunwei Deng and Chengjiang Yin for their technical helps. ZYZ and LZ acknowledge support of the National Natural Science Foundation of China (NSFC) under grants numbers 12173016, 12041305. ZYZ and LZ acknowledge the science research grants from the China Manned Space Project (number CMS-CSST-2021-A08). ZYZ and LZ acknowledge the Program for Innovative Talents, Entrepreneur in Jiangsu. CY acknowledges support from ERC advanced grant number 789410. DDX acknowledges the NSFC grant number12073013. This work is supported by the China Manned Space Project (number CMS-CSST-2021-A07). RL acknowledges the support of National Nature Science Foundation of China (numbers 11988101,11773032,12022306), the support from the Ministry of Science and Technology of China (number 2020SKA0110100), the science research grants from the China Manned Space Project (numbers CMS-CSST-2021-B01, CMS-CSST-2021-A01), CAS Project for Young Scientists in Basic Research (number YSBR-062). Z-CZ was supported by the National Natural Science Foundation of China (grant number 12233002).

This work used the following software: Autolens (Nightingale & Dye 2015; Nightingale et al. 2018, 2021b), casa (McMullin et al. 2007), SAOimageDS9 (Joye & Mandel 2003), astropy (Astropy Collaboration 2013; Price-Whelan et al. 2018), dynesty (Speagle 2020), corner (Foreman-Mackey 2016), emcee (Foreman-Mackey et al. 2013), matplotlib (Hunter 2007), numpy (van der Walt, Colbert & Varoquaux 2011), pyautofit (Nightingale, Hayes & Griffiths 2021a), pylops (Ravasi & Vasconcelos 2020), python (Van Rossum & Drake 2009), and scipy (Virtanen et al. 2020).

DATA AVAILABILITY

All data adopted in this work are publicly available in corresponding data archival systems of HST (https://hla.stsci.edu/hlaview.html), ALMA (https://almascience.eso.org/aq/), VLA (https://data.nrao.edu/portal/), and e-Merlin (https://www.e-merlin.ac.uk/archive/).

Footnotes

1

https://data.nrao.edu/portal

2

https://github.com/Jammy2211/PyAutoLens

3

Both parametric and non-parametric fitting methods are provided as Jupyter notebooks at the following link:https://github.com/Jammy2211/autolens_likelihood_function

APPENDIX A: MOCK DATA

In Fig. A1, we present the mock data of the lensing system of Cloverleaf, to simulate the response of a background source with a Sérsic profile. The lens model is the same as the one in our non-parametric model. The source is a Sérsic ellipsoid moving from the north-west part to the south-east part of the inner caustics in the source plane, and its centre is marked with red crosses as moving.

Figure A1.

Mock data of this lensing system: The lens model is the same as the one in our non-parametric model. The source is a Sérsic ellipsoid moving from the north-west part to the south-east part of the inner caustics in the source plane, and its centre is marked with red crosses as moving. The solid black line and the red solid line are the outer and the inner critical curve (which is incorrect due to the numerical issue and should be a point) in the lens plane, respectively. The solid black line (diamond-shaped) in the source plane is the inner caustics.

APPENDIX B: RESIDUAL

In Fig. B1, we present the residual maps of our lensing models of the 300, 1.5, 8.4, and 33 GHz data.

Figure B1.

Residual maps of non-parametric models at 300, 1.5, 8.4, and 33 GHz: The red dot circle is the mask. The solid black line and the red solid line are the outer critical curve and the inner critical curve (which is incorrect due to the numerical issue and should be a point), respectively.