Multi-scale Radio and X-Ray Structure of the High-redshift Quasar PMN J0909+0354

The quasar J0909+0354 has been detected with ∼100 mJy level flux densities at various fRequencies in different sky surveys. Flux densities from the following observations are listed in Table 1. Karl G. Jansky Very Large Array (VLA) A-configuration imaging observations at 1.5, 4.9, and 8 GHz revealed that the radio emission of the quasar can be resolved into a double structure: a compact core and a secondary component at about 2'' angular separation in the north-northwestern direction (Gobeille et al. 2014). Archival Very Long Baseline Array (VLBA) observations (Petrov 2021, project code: BP171AB, PI: L. Petrov) of J0909+0354 at 4.3 and 7.6 GHz show compact, unresolved radio emission up to ∼10–20 pc, with 111 mJy and 76 mJy total flux densities, respectively. ⁹ Global VLBI observations at 5 GHz resolved the parsec-scale morphology, revealing a compact synchrotron self-absorbed core and a more diffuse jet structure visible up to ∼10 pc away from the core (Paragi et al. 1999).

X-ray emission of the quasar was first detected by ROSAT (ROentgen SATellite) in the 2–4 keV energy range with a flux of F = 9.7 ± 2.7 × 10⁻¹³ erg cm⁻² s⁻¹ (Brinkmann et al. 1997). Fluxes from observations by the BeppoSAX and Swift space telescopes are 1.9 × 10⁻¹² erg cm⁻² s⁻¹ (photon index γ = 1.16 ± 0.2, Donato et al. 2005) and 15.65 × 10⁻¹² erg cm⁻² s⁻¹ ($\gamma ={1.88}_{-2.28}^{+1.51}$, Oh et al. 2018), in the 0.1–50 keV and 14–195 keV energy ranges, respectively. The light curve of J0909+0354 observed in the framework of the 105 month Swift-BAT all-sky hard X-ray survey revealed some variability in the 14–195 keV energy range (Figure 1).

Zoom In Zoom Out Reset image size

The quasar J0909+0354 is also a source of electromagnetic radiation detected in various surveys in the ultraviolet, optical, and infrared wave bands, which are listed in Table 2.

The quasar was not detected in γ-rays with the INTErnational Gamma-Ray Astrophysics Laboratory (INTEGRAL, Winkler et al. 2003), the Large Area Telescope of Fermi Gamma-Ray Space Telescope (Fermi-LAT Atwood et al. 2009), the Compton Gamma Ray Observatory (CGRO, e.g., Kanbach et al. 1989; Meegan et al. 1998; Thompson et al. 1993; Goldstein et al. 2013), or the Astro-rivelatore Gamma aImmagini LEggero mission (AGILE, Tavani et al. 2009).

For a quantitative comparison of radio emission between parsec and kiloparsec scales, we utilized data obtained with the VLA at 1.5, 6.2, and 8.5 GHz.

The 8.5 GHz observations were carried out on 1998 March 14, in the framework of the Cosmic Lens All-Sky Survey (CLASS, Myers et al. 2003; project code: AM593, PI: S. Myers), in which 27 stations participated. The on-source time for J0909+0354 was 39 s. Two intermediate frequency channels (IFs) were used with one spectral channel in each, the total bandwidth was 50 MHz. The data were recorded in full polarization with a 3.3 s integration time. We calibrated the phases and amplitudes with the National Radio Astronomy Observatory (NRAO) Astronomical Image Processing System ¹⁰ (AIPS, e.g., Diamond 1995; Greisen 2003) package, following the steps of standard data reduction described in the cookbook ¹¹ for VLA continuum data, using 3C 48 as primary flux density calibrator.

The 6.2 GHz observations were conducted on 2012 November 18 (project code: 12B-230, PI: J. Wardle) with the participation of 26 telescopes. In the 3.47 hr long observing run, 39 sources were targeted including calibrators, from which the on-source time for J0909+0354 was 316 s. A total of 16 IFs were used with 64 spectral channels in each IF. The total bandwidth was 128 MHz. The data were recorded in full polarization, and were correlated with a 1 s averaging time.

The measurements at 1.5 GHz were conducted on 2016 October 25 (project code: 16B-015, PI: J. S. Farnes), with 26 antennae participating. From the total 10 hr of observation, the on-source time for J0909+0354 was 126 s. The 64 MHz total bandwidth was divided into 16 IFs, with 64 spectral channels in each. The data were recorded in full polarization, and were correlated with a 1 s averaging time.

We calibrated the phases and amplitudes of the 1.5 and 6.2 GHz VLA data sets with the Common Astronomy Software Applications ¹² (CASA, McMullin et al. 2007) software using 3C 48, J0738+1741, and J0831+0429 as calibrators, following the steps of standard data reduction for VLA continuum observations. ¹³ Then we exported the data to uvfits format.

The calibrated visibilities of all three VLA observations were then imported to the DIFMAP program ¹⁴ (Shepherd 1997), where we carried out hybrid mapping with cycles of phase and amplitude self-calibration and imaging (applying the clean deconvolution method by Högbom 1974). Finally, to quantitatively describe the brightness distribution of the source, we fitted circular Gaussian model components directly to the self-calibrated visibility data (Pearson 1995). Uncertainties for the fitted model parameters were calculated following Lee et al. (2008).

To study the parsec-scale radio structure of the quasar J0909+0354, we used archival data as well as new observations made by various VLBI arrays. The latest and most sensitive data set was acquired at 5 GHz with the European VLBI Network (EVN) on 2019 March 1 (project code: EP115, PI: K. Perger). The observation lasted for a total of 6 hr and involved 15 radio telescopes: Jodrell Bank Mk2 (United Kingdom), Westerbork (The Netherlands), Effelsberg (Germany), Medicina, Noto (Italy), Onsala (Sweden), Tianma, Nanshan (China), Toruń (Poland), Yebes (Spain), Svetloe, Zelenchukskaya, Badary (Russia), Hartebeesthoek (South Africa), and Irbene (Latvia). The on-source integration time was 5.24 hr. The data were recorded at a rate of 1024 Mbit s⁻¹ in left and right circular polarizations. The total bandwidth was 16 MHz per polarization in 32 spectral channels per IF, and a total of 8 IFs were used. The data were correlated with a 2 s averaging time at the EVN Data Processor at the Joint Institute for VLBI European Research Infrastructure Consortium (Dwingeloo, the Netherlands).

We calibrated the phases and amplitudes of the visibilities in the AIPS package. We applied a priori amplitude calibration based on radio telescope gain curves and measured system temperatures, then removed interchannel delay and phase offsets using a 1 minute data segment of a bright calibrator source (J0909+0121). Visual inspection and flagging of the visibilities were followed by global fringe fitting (Schwab & Cotton 1983) on the target source J0909+0354. The calibrated data were exported into uvfits format for further analysis.

For a comparison, we also recovered and analyzed the archival data observed by a global VLBI network on 1992 September 27–28 (Paragi et al. 1999). Nine telescopes, Effelsberg, Medicina, Onsala, the phased array of the Westerbork Synthesis Radio Telescope, as well as Green Bank, Haystack, Owens Valley, and the phased VLA (the latter four in the USA) participated in the observations, which were part of a 24 hr long experiment. The on-source time for J0909+0354 was 3 hr. The data were recorded in left circular polarization with a total bandwidth of 28 MHz in 7 IFs, and were correlated in the Max Planck Institute for Radio Astronomy (Bonn, Germany). For further analysis, we used the visibility data calibrated by Paragi et al. (1999).

To supplement our long-track VLBI observations for studying possible changes in the parsec-scale radio structure of J0909+0354, we also analyzed archival "snapshot" data obtained with the VLBA. These observations were conducted in the framework of the 7th VLBA Calibrator Survey ¹⁵ (Petrov 2021, project code: BP171AB, PI: L. Petrov) on 2013 April 28. All ten telescopes of the array participated in the dual-frequency (4.3 and 7.6 GHz) observation that was carried out in right circular polarization, with an on-source time of 1 minute. The total bandwidth was 32 MHz in 8 IFs. The a priori calibrated visibility data sets were produced by the PIMA v2.03 software (Petrov et al. 2011) and were obtained from the Astrogeo VLBI image database. ¹⁶

All four calibrated VLBI data sets were then imported to the DIFMAP program for phase and amplitude self-calibration, imaging, and model fitting, similar to the VLA data treatment described above. Errors for the fitted model parameters were calculated following Lee et al. (2008), considering an additional 5% calibration uncertainty for flux densities.

The X-ray emission associated with the quasar J0909+0354 was observed with the Chandra Advanced CCD Imaging Spectrometer (ACIS) as part of a survey of radio-loud quasars at z > 3 (ObsID 20404, PI: D. Schwartz). We then used 77.5 ks follow-up observations (ObsIDs 22568, 23161, and 23162, PI: D. Schwartz) to reveal the extended X-ray features. The latter observations took place on 2020 February 17, 18, and 20. The data were reduced with CIAO (Chandra Interactive Analysis of Observations) version 4.12 (Fruscione et al. 2006), using SHERPA version 2 (Doe et al. 2007). Imaging used SAOImage ds9 version 8.2b1 (Joye & Mandel 2005). Background was determined to be 0.0591 ± 0.0013 counts arcsec⁻² from two rectangular regions totaling 34,131 arcsec² away from the quasar and not including sources.

Model fitting to the visibilities of the 1.5, 6.2, and 8.5 GHz VLA observations resulted in two components for all three data sets: a compact core and a slightly more diffuse feature (discussed in detail in Section 5.5) at ∼233 to the north-northwestern (NNW) direction (position angle ∼−165 as measured from north through east, Figure 2). The sums of the flux densities of the circular Gaussian model components describing the kiloparsec-scale radio structure is 153.1 ± 14.4 mJy, 198.9 ± 12.6 mJy, and 211.1 ± 13.3 mJy, for the 1.5, 6.2, and 8.5 GHz data, respectively. ¹⁷ The properties of the two fitted components for the VLA data sets are listed in Table 3. As it was previously noted by Gobeille et al. (2014), we examined the nature of NNW, and found that its three-point spectral index ¹⁸ between 1.5 and 8.5 GHz is α_NNW = −1.08 ± 0.17. The three-point spectral index for the core between 1.5 and 8.5 GHz is ${\alpha }_{\mathrm{core}}=0.19\pm 0.01$. To avoid the possible effect of time variability, we also determined the spectral index for this component by separately processing the first and last two IFs from the same 6.2 GHz observation. Repeating the hybrid mapping and Gaussian model fitting on these two subsets, we found the values ${\alpha }_{7.77,\mathrm{core}}^{4.62}=-0.51\pm 0.05$ and ${\alpha }_{7.77,\mathrm{NNW}}^{4.62}=-1.03\pm 0.10$ for the core and NNW components, respectively. The latter value is consistent with the identification of the NNW component as a steep-spectrum jet hotspot—see further discussion on this in Section 5.5 below.

Zoom In Zoom Out Reset image size

Table 3. Model Parameters of J0909+0365 from the VLA Observations

		R.A.	Decl.	ϑ	S	R	ϕ
					(mJy)	('')	(°)
1.5 GHz	core	9h9m1591 ± 0001	$3^\circ 54^{\prime} 43\buildrel{\prime\prime}\over{.} 3\pm 0\buildrel{\prime\prime}\over{.} 01$	<0013 ± 001	133.1 ± 13.75	0	0
	NNW	9h9m1587 ± 002	$3^\circ 54^{\prime} 45\buildrel{\prime\prime}\over{.} 5\pm 0\buildrel{\prime\prime}\over{.} 11$	082 ± 022	20.0 ± 4.4	2.27 ± 0.11	−13.5 ± 2.7
6.2 GHz	core	9h9m1591 ± 0001	$3^\circ 54^{\prime} 43\buildrel{\prime\prime}\over{.} 0\pm 0\buildrel{\prime\prime}\over{.} 01$	0048 ± 0002	193.6 ± 12.5	0	0
	NNW	9h9m1587 ± 004	$3^\circ 54^{\prime} 45\buildrel{\prime\prime}\over{.} 2\pm 0\buildrel{\prime\prime}\over{.} 07$	045 ± 013	5.3 ± 1.2	2.33 ± 0.04	−16.5 ± 1.0
8.5 GHz	core	9h9m1591 ± 0001	$3^\circ 54^{\prime} 43\buildrel{\prime\prime}\over{.} 2\pm 0\buildrel{\prime\prime}\over{.} 01$	<0010 ± 001	208.3 ± 13.25	0	0
	NNW	9h9m1586 ± 001	$3^\circ 54^{\prime} 45\buildrel{\prime\prime}\over{.} 5\pm 0\buildrel{\prime\prime}\over{.} 02$	<009 ± 004	2.8 ± 1.0	2.37 ± 0.02	−17.0 ± 0.4

Note. Column 1—model component name, Columns 2, 3—coordinates (R.A., and decl.), Column 4—circular Gaussian model component size (FWHM); due to its small size, the uncertainty of the core component is given as the relative astrometric precision of the VLA (calculated following the error estimation of Brogan et al. 2018). We note that these uncertainties are in the same order of magnitude as the values calculated following the VLA observation guide (http://science.nrao.edu/facilities/vla/docs/manuals/oss/performance/positional-accuracy), i.e., 10% of the FWHM of the restoring beam, Column 5—flux density, Column 6—radial distance from the core component, Column 7—model position angle with respect to the core measured from north through east.

Download table as:  ASCIITypeset image

The Chandra 0.5−7 keV X-ray image of J0909+0354 is shown in Figure 3, and reveals three distinct features. There is bright emission from the quasar, from a point-like source at a 648 angular distance at position angle 55°, and emission extending from the quasar 23 toward the NNW component at position angle −174. We interpret the latter feature as a kiloparsec-scale jet, with enhanced emission at its end coincident with NNW.

Zoom In Zoom Out Reset image size

To establish the reality of the jet, we performed a high fidelity simulation of the quasar, using saotrace-2.0.4_03 ¹⁹ to generate rays that are passed to marx-5.5.0 ²⁰ (Davis et al. 2012) to simulate an ACIS-S image. We measure 9539 counts in a 14 circle around the quasar, 54 in a 24 × 06 box region defined between the core and NNW, and 28 in the circle around the NNW radio emission at the end of the jet. Normalizing the simulated counts to the number in the quasar (i.e., the core component), we predict that 36.5 counts from the quasar scatter into the jet box, and that 6.4 scatter into the NNW region, giving a net 17.2 ± 9.5 counts for the box area, and 21.5 ± 5.9 in the NNW region. Interpreting all of the emission to the NNW as a single jet would result in 39.0 ± 10.4 counts. With the average background of 0.03 counts, detections of X-ray photons are significant (≥3σ) in the box and NNW regions, either considering the emission as an ensemble continuous jet or as two separate features. Taking the dimensions of the box area and the 04 extraction circle around NNW into account, there are 11.9 counts arcsec⁻² and 42.8 counts arcsec⁻² in the two features, respectively, implying an enhancement at the latter component.

The quasar fits a power-law photon spectrum ∝ E^{−1.24±0.03}, using the galactic absorption n_H = 3.5 × 10²⁰ cm⁻² and also fitting an intrinsic absorption of $({2.7}_{-1.7}^{+1.8})\times {10}^{22}$ cm⁻² at the source. The measured 0.5–7 keV flux of 1.76 × 10⁻¹² erg cm⁻² s⁻¹ corresponds to a luminosity of 1.78 × 10⁴⁷ erg s⁻¹ in the 2.1–30 keV rest-frame band. This value is consistent with previous X-ray measurements in similar bandwidths. Attributing all of the emission to the NNW to an X-ray jet, it would have a photon spectrum ∝ E^{−1.71±0.30} and a flux of (8.5 ± 2.5) × 10⁻¹⁵ erg cm⁻² s⁻¹. Such a jet flux would be 0.5% of the quasar flux. This is consistent with but somewhat lower than the median 2% jet-to-core ratio from a sample of quasars at lower redshifts (Marshall et al. 2018). We cannot exclude a significant contribution to the jet flux from a point-like X-ray source at the NNW component. Taking the extraction circle of 04 radius about the NNW region, we fit a photon spectrum ∝ E^{−1.42±0.44}. Considered as an upper limit, such a source could be contributing 4.5 × 10⁻¹⁵ erg cm⁻² s⁻¹ of the flux extending to the NNW.

The source to the NE fits a power law with a photon index of 1.81 ± 0.33, giving a flux of 2.03 × 10⁻¹⁴ erg cm⁻² s⁻¹.

In the original publication of the global VLBI observations, Paragi et al. (1999) identified multiple jet components in the parsec-scale radio structure of the quasar J0909+0354. However, a conclusive analysis of the jet structure extended to the north was hampered by the lack of the long north–south baselines in the global VLBI array. This resulted in a relatively poorer angular resolution in that position angle. After our re-analysis of the 1992 data, the best Gaussian fit to the visibility data provided three model components (Figure 4(c), Table 4). The integral flux density based on the fitted circular Gaussian model components is 52.5 ± 5.2 mJy. Since the major axis of the elongated restoring beam coincided with the jet direction, additional observations were needed for the reliable characterization of the inner jet structure in J0909+0354.

Zoom In Zoom Out Reset image size

The longer north–south interferometric baselines of the 2019 EVN observation resulted in a nearly circular restoring beam (Figure 4(b)). Three circular Gaussian components were found as the best-fit model of the brightness distribution (Table 4). The higher north–south resolution of the EVN observations allows us to conduct a detailed analysis of the parsec-scale morphology of the quasar. The jet propagates toward the north, and is clearly detectable up to ∼20 mas (∼150 pc) away from the core. It becomes more diffuse with increasing distance with respect to the core, apparently splitting into two branches like a fork at ∼10 mas. The peak intensity in the image is 64.1 ± 1.5 mJy beam⁻¹, the integral flux density of the fitted model components is 70.2 ± 3.6 mJy (within a ∼65 pc projected distance; Figure 4(b), Table 4).

For a qualitative comparison of the 1992 and 2019 images obtained at the same observing frequency, we restored the clean map of the most recent 2019 observations with the same elongated beam of the 1992 global VLBI observations (Figure 4(d)). The overall structure of the source is clearly similar at both epochs, with more diffuse emission detected in the outer regions of the jet in 2019 (Figures 4(c) and (d)). This can be attributed to the higher sensitivity of the new EVN observations compared to that of the old global VLBI observations.

The structure of the quasar appears less resolved in the 4.3 (Figure 4(a)) and 7.6 GHz VLBA (Figure 5) images compared to that in the global (Figure 4(c)) and EVN (Figure 4(b)) measurements, due to the larger restoring beam that is more elongated to the north–south direction. However, the jet direction toward the north is still discernible. Model fitting to the 4.3 GHz VLBA data set resulted in 3 components (Table 4); the peak intensity is 101.7 ± 4.3 mJy beam⁻¹, while the integral flux density of the model components is 108.9 ± 7.2 mJy (within a ∼50 pc projected distance). Our best-fit model for the 7.6 GHz VLBA data consists of 2 components (Table 4). The integral flux density of the fitted Gaussian model components is 75.1 ± 5.5 mJy (within a ∼10 pc projected distance), while the peak intensity of the clean map is 70.0 ± 3.5 mJy beam⁻¹.

Zoom In Zoom Out Reset image size

Using the fitted model parameters (Table 4), we calculated the apparent brightness temperatures (T_b) for the core components at each VLBI epoch, applying the formula

$$\begin{eqnarray}&&{T}_{{\rm{b}}}=1.22\times {10}^{12}\,(1+z)\,\displaystyle \frac{S}{{\vartheta }^{2}{\nu }^{2}}\,{\rm{K}}\end{eqnarray} \tag{ 1 }$$

(e.g., Condon et al. 1982). Here, S is the flux density in Jy, ϑ is the angular size of the fitted circular Gaussian model component (FWHM) in mas, and ν is the observing frequency in GHz. Values for the global, VLBA, and EVN observations are T_B,G > 0.13 × 10¹¹ K, T_B,VLBA,4.3 = 3.2 × 10¹¹ K, T_B,VLBA,7.6 = 2.1 × 10¹¹ K, and T_B,EVN = 2.9 × 10¹¹ K, respectively.

We also calculated Doppler factors according to the equation

$$\begin{eqnarray}&&\delta =\displaystyle \frac{{T}_{{\rm{b}}}}{{T}_{{\rm{b}},\mathrm{int}}},\end{eqnarray} \tag{ 2 }$$

where T_b,int is the intrinsic brightness temperature of the source. Lower and upper limits of the Doppler factor were given by applying T_b,int ≈ 5 × 10¹⁰ K (Readhead 1994) and T_b,int ≈ 3 × 10¹⁰ K (Homan et al. 2006). The former corresponds to the equipartition state between the energy densities of the emitting plasma and the magnetic field, while the latter is a characteristic value determined for parsec-scale AGN jets based on VLBI observations. The resulting values of the Doppler factor for the VLBA and EVN data sets are between δ_eq = 4–6 and δ_char = 7–10 (Table 4, Cols. 8 and 9). Constrained by the upper limit on the model component FWHM of the global VLBI measurement, only lower limits were found for the values of brightness temperatures and Doppler factors in 1992.

Applying the following equation (Ghisellini et al. 1993),

$$\begin{eqnarray}&&{\delta }_{\mathrm{IC}}=f(\alpha ){S}_{{\rm{r}}}{\left(\displaystyle \frac{\mathrm{ln}\tfrac{{\nu }_{{\rm{b}}}}{{\nu }_{{\rm{r}}}}}{{S}_{{\rm{X}}}{\vartheta }^{6-4\alpha }{\nu }_{{\rm{X}}}^{-\alpha }{\nu }_{{\rm{r}}}^{5-3\alpha }}\right)}^{\displaystyle \frac{1}{4-2\alpha }}(1+z),\end{eqnarray} \tag{ 3 }$$

a lower limit to the Doppler factor can be calculated, assuming that the X-ray emission originates from the synchrotron self-Compton process of the quasar jet, ²¹ where S_r, ν_r, and ϑ are the flux density (Jy), frequency (GHz) and FWHM diameter (mas) of the core model component from the given VLBI observation, respectively, S_X and ν_X are the flux density (Jy) and energy (keV) of the X-ray emission, and α = −0.75 is the spectral index (assumed value, e.g., Ghisellini et al. 1993), ν_b = 10⁵ GHz is the cut-off frequency of the high-energy radiation, and f(α) = −0.08α + 0.14. Values of the inverse-Compton Doppler factor are δ_IC > 2 for both the EVN and VLBA observations, while it is δ_IC > 0.1 for the global VLBI data set, making the two independent estimates of the Doppler factor consistent with each other.

Both the 1992 global VLBI and 2019 EVN observations show a complex parsec-scale morphology with a compact core and a jet extending northwards, which is also hinted at by the morphology in the 4.3 GHz VLBA image (Figures 4(a)–(c), respectively). However, as the major axis position angle of the elongated restoring beams in both the global and the VLBA snapshot observations almost perfectly coincides with the jet direction, in the following we will only discuss the jet structure in detail based on the new EVN observations. The milliarcsecond-scale radio image at 5 GHz (Figure 4(b)) reveals a jet morphology extending up to ∼20 mas (150 pc) toward the north, with respect to the core component (i.e., the synchrotron self-absorbed base of the radio jet).

The jet becomes more diffuse at ∼10 mas (∼75 pc), and apparently splits into two branches. The observed division can be explained by the radio jet interacting with a denser region of the surrounding interstellar medium (e.g., Attridge et al. 1999; Dallacasa et al. 2013). Alternatively, the fork-like morphology can be attributed to the spine–sheath structure of the jet. The model, supported by numerical simulations as well (e.g., McKinney 2006; Komissarov et al. 2007), states that the inner region (spine) of the AGN jet propagates with relativistic speeds, while it is surrounded by a slower sub-relativistic sheath (e.g., Komissarov 1990). A similar structure can be found in the parsec-scale jets of, e.g., Mrk 501 (Giroletti et al. 2004), 3C 66A (0219+428) and 3C 380 (1828+487; Figures 2 and 4 in Lister et al. 2013), 4C 76.03 (0404+768; Figure 5 in Dallacasa et al. 2013), 3C 84 (Nagai et al. 2014), 1308+326 (Britzen et al. 2017), and S5 0836+710 (Vega-García et al. 2020). The spine–sheath structure has recently been reported in one of the best studied cases of the core–jet morphology in the quasar 3C 273 (Bruni et al. 2021). Seven out of these eight quasars have complex morphology on kiloparsec scales, i.e., hotspots, multiple components, extended jets, or radio lobes (Kellermann et al. 1971; Pedlar et al. 1990; Wilkinson et al. 1991; Hummel et al. 1992; Murphy et al. 1993; Price et al. 1993; Akujor & Garrington 1995; Xu et al. 1995; Taylor et al. 1996; Cassaro et al. 1999; Perucho et al. 2012; Perley & Meisenheimer 2017), while 4C 76.03 has a compact structure, unresolved with the VLA (Xu et al. 1995). Similar to our target source, J0909+0354, a single hotspot was identified in the kiloparsec-scale jet of 3C 273 (e.g., Meisenheimer & Heavens 1986; Perley & Meisenheimer 2017).

The structured jet has a pronounced footprint in the linearly polarized emission (e.g., Pushkarev et al. 2005; Murphy et al. 2013), and it is also visible in the full polarized intensity in the manner of the relative brightening of the outer regions further away from the jet axis (Giroletti et al. 2004; Ghisellini et al. 2005; Nagai et al. 2014; Giovannini et al. 2018; Ros et al. 2020). The influence of external processes (such as the effect of the surrounding medium) on the observed properties of the jet are negligible in the spine–sheath model (e.g., McKinney 2006).

A similar approach was discussed in the framework of the two-fluid jet model by Pelletier & Roland (1989), in which both the superluminal motion at parsec scales and hotspots at kiloparsec scales are explained by an outer thermal electron–proton flow (propagating at nonrelativistic speeds, called jet) and an inner relativistic electron–positron plasma (called beam). The apparent split in the parsec-scale radio jet (also referred to as limb-brightening) occurs at a distance where the magnetic field becomes weaker than a critical value, hence allowing the relativistic and thermal components to interflow. The two-fluid model accounts the observed one-sidedness to the different fraction of the relativistic flow components in the jet and counter-jet (at parsec scales), and thus the asymmetric re-acceleration of the thermal flow (resulting in hotspots at kiloparsec scales), rather than the effect of Doppler-beaming/deboosting. The model was recently applied to, e.g., 3C 273 (Bruni et al. 2021), was proposed to explain (Kharb et al. 2015) the structure of multiple sources showing both blazar and Fanaroff–Riley type II characteristics (e.g., Landt et al. 2006; Kharb et al. 2010, 2015), and was addressed in the discussion of the blazar PKS 0735+178, which is reported to show a Fanaroff–Riley type II kiloparsec structure (Goyal et al. 2009).

To further emphasize the diffuse emission of the parsec-scale jet, we applied a Gaussian (u, v) taper (i.e., a scheme where the weights of the visibilities decrease as a function of (u, v) radius) with a value of 0.2 at 10 million wavelength radius to the EVN data set, and repeated the hybrid mapping and Gaussian model fitting procedure. The tapered image (Figure 6) confirms the bending of the jet, which is indicated by the weak (0.35 ± 0.08 mJy flux density) model component found at 34.5 ± 2.4 mas from the core (at a 2° ± 4° position angle).

Zoom In Zoom Out Reset image size

The relative positions of VLBI model components (Table 4) and the overall shape of the jet imply that the emission continues toward the NNW component identified in the VLA radio map (as well as in the Chandra image), but is resolved out by the EVN between ∼100 pc and kiloparsec scales. Moreover, model component positions indicate a bending trajectory: the jet shows a slight turn toward the northwestern direction by 30° between 65 pc (C1) and 20 kpc (NNW). Such an apparent bending has been detected in numerous blazar jets (e.g., Conway & Murphy 1993; Hong et al. 2004; Kharb et al. 2010; Zhao et al. 2011; Singal 2016; Perger et al. 2018). In our case, such a morphology might appear interrupted between the components C1 and NNW due to an insufficient brightness sensitivity of the observing system.

We note that the positional discrepancy between the two VLBA data sets, the absence of the outer jet component at 7.6 GHz, and the large difference between the core flux densities at the two observing frequency bands can be explained as some of the inner jet components at 7.6 GHz are blended into the core at 4.3 GHz, and that the outer components (i.e., the counterparts of B1 and B2) are too faint to be detected at 7.6 GHz due to the steepening of the spectrum of the jet further away from the core.

Brightness temperatures determined from VLBA and EVN measurements well exceed both the theoretical (T_b,int ≈ 5 × 10¹⁰ K, Readhead 1994) and the somewhat lower empirical (T_b,int ≈ 3 × 10¹⁰ K, Homan et al. 2006) limits. The relativistic enhancement is thus clearly indicated by the high values of the Doppler factor (δ_eq = 4–6, δ_char = 7–10). Flux densities at different-epoch VLBI observations reveal the variability of the source at parsec scales. We note that due to the improper resolution of the interferometer (i.e., the upper limit on the FWHM diameter of the core), only lower limits could be determined for the brightness temperature and Doppler factors of the global VLBI observation; therefore, we excluded these data from further analysis of the jet parameters.

Using parameters derived from the models fitted to the visibility data of the EVN observation, we estimated the inclination angle of the jet with respect to the line of sight of the observer. We chose the value of the bulk Lorentz factor to be between Γ = 5 and Γ = 15 (typical values found for high-redshift AGN jets, e.g., Volonteri et al. 2011). The estimated inclination angle of the jet is in the range 8° ≤ θ_eq ≤ 14° and 0° ≤ θ_char ≤ 7° assuming the values of the equipartition and empirical Doppler factors, respectively. Assuming the empirical upper limit on the bulk Lorentz factor (Γ = 25, determined in a parsec-scale proper-motion study of a large sample of AGN jets, Kellermann et al. 2004) for our calculation, the viewing angle of the jet is constrained to the ranges of 6° ≤ θ_eq ≤ 8° and 5° ≤ θ_char ≤ 6°, for the equipartition and empirical Doppler factors, respectively. The possible ranges of jet parameters are illustrated in Figure 7. We note that using the inverse-Compton Doppler factor (δ_IC > 2), upper limits on the jet inclination angle of θ_IC ≤ 23°, θ_IC ≤ 14°, and θ_IC ≤ 11° can be derived by applying bulk Lorentz factors of Γ = 5, Γ = 15, and Γ = 25, respectively.

Zoom In Zoom Out Reset image size

Based on the circular Gaussian model components fitted to the 2013 VLBA and 2019 EVN visibility data, we estimated the apparent proper motion in the parsec-scale jet of J0909+0354. Although including the 1992 global VLBI observations could make the estimation more accurate (due to the longer time span), we do not consider this model as a starting point, because of the unfortunate network geometry.

Since both the 2013 and 2019 data are characterized by a core and two additional jet components, we assumed that the B1 and B2 components in 2013 (VLBA data, 4.3 GHz) correspond to C1 and C2 in 2019 (EVN data, 5 GHz), respectively. This is further supported by the fact that the component position angles are equal within the uncertainties (Table 4, Col. 5). Over the 5.84 yr time span, the calculated proper-motion values for the components B1–C1 and B2–C2 are μ₁ = 0.30 ± 0.03 mas yr⁻¹ and μ₂ = 0.37 ± 0.02 mas yr⁻¹, respectively. These correspond to β₁ = (31 ± 3)c and β₂ = (39 ± 3)c apparent superluminal speeds, considering the cosmological time dilatation. Although apparent proper-motion values as high as our estimates (β₁ and β₂) are presented in the literature (e.g., 13 AGN with β_app > 20, 3 AGN with β_app > 30, Kellermann et al. 2004), both the lowest and highest values of the bulk Lorentz factor determined for either of the component transverse speeds (${{\rm{\Gamma }}}_{\min }=53$ and ${{\rm{\Gamma }}}_{\max }=192$) are much higher than the range of 5 ≲ Γ ≲ 15 determined for z > 3 AGN jets (Volonteri et al. 2011), and even the lowest value well exceeds the empirical maximum of Γ = 25 (e.g., Kellermann et al. 2004; Lister et al. 2016; Pushkarev et al. 2017), thus making the determined proper motions rather questionable. We note that due to the slightly different frequency and the different resolutions, the apparent core position (with the core–jet emission blended together) may be different in the VLBA and EVN radio maps, resulting in an apparently different core–jet separation. In any case, new sensitive follow-up VLBI observations to be conducted at 5 GHz in the next 5–10 yr could settle the issue.

The broadband and radio spectra of the quasar are shown in the upper and lower panels of Figure 8, respectively. Data were obtained from the NASA/IPAC Extragalactic Database, ²² the photometry tool of the VizieR service, ²³ and from Chandra, VLA, global VLBI, VLBA, and EVN observations discussed in this paper.

Zoom In Zoom Out Reset image size

The differences in total flux densities at similar frequencies but different epochs can be partly attributed to the different angular resolutions of the interferometer arrays (FIRST, NVSS, and CLASS, Condon et al. 1998; Myers et al. 2003; Helfand et al. 2015) and single-dish (PMN, GBT, Becker et al. 1991; White & Becker 1992; Griffith et al. 1995; Gregory et al. 1996) observations. We fitted the logarithm of the flux density values of the entire radio wave band (from 150 MHz to 8.4 GHz) with a linear curve, and found a spectral index of α_kpc = −0.13 ± 0.06 for the kiloparsec-scale structure.

The flux density variability of the parsec-scale emission is 35–50% between the epochs of the global VLBI, VLBA, and EVN observations at ∼5 GHz. We calculated radio spectral indices for the parsec-scale structure using the flux densities determined from the simultaneous dual-frequency VLBA observations and the new EVN data. The power-law spectral index in the 4.34–7.62 GHz frequency range is ${\alpha }_{4.34}^{7.62}\,=-0.66\pm 0.17$ and α_pc = − 0.46 ± 0.13, for the VLBA and VLBA–EVN total flux density data, respectively. Similarly, spectral indices for the core components are ${\alpha }_{4.34,\mathrm{core}}^{7.62}\,=-0.74\pm 0.22$ and ${\alpha }_{\mathrm{pc},\mathrm{core}}=-0.55\pm 0.15$. On the one hand, as the quasar shows significant flux density variability on parsec scales (Table 4), these values should be treated with caution. Although parallel observations with the VLBA are not affected by variability, the steep spectrum can be the result of the different angular resolutions at the two frequencies, thus the comparisons of the components and flux densities remain uncertain. On the other hand, steep spectra with $-0.52\,\leqslant {\alpha }_{{\rm{c}}{\rm{o}}{\rm{r}}{\rm{e}}}\leqslant -1$ were previously found for high-redshift (z ≥ 4) quasars, albeit not those with Doppler-boosted radio emission (e.g., Frey et al. 2003, 2005, 2008, 2010, 2011; Coppejans et al. 2016; Cao et al. 2017), the blazar PSO 0309+27 (at z = 6.1, Spingola et al. 2020) between 1.5 GHz and 5 GHz frequencies, as well as J0906+6930 (z = 5.47) in the 15 GHz ≤ ν ≤ 43 GHz frequency range (although its core spectral index flattens to ${\alpha }_{\mathrm{core}}=0.2$ below 8.4 GHz, Zhang et al. 2017).

As noted in Section 4.1, the radio spectrum of the NNW component between 1.5 and 8.5 GHz is consistent with its identification as a jet hotspot. Not that many high-redshift sources demonstrate well-pronounced kiloparsec-scale jets and/or extended morphological features that might be suspected as being jets. An increasing rareness of detectable jets at both parsec and kiloparsec scales in AGNs with increasing redshift is rather natural (e.g., Gurvits 1999; Gurvits et al. 2015). Nevertheless, such examples are known: J1405+0415 at z = 3.215 (Gurvits et al. 1992; Yang et al. 2008), J2217+0220 at z = 3.572 (Lobanov et al. 2001), and J2134−0419 at z=4.33 (Perger et al. 2018), to list just three of them.

Two-point spectral indices (${\alpha }_{6.2\mathrm{NNW}}^{1.5}=-0.94\pm 0.19$, ${\alpha }_{7.77,\mathrm{NNW}}^{4.62}=-1.03\pm 0.10$, and ${\alpha }_{8.5,\mathrm{NNW}}^{6.2}=-2.05\pm 0.25$) calculated from the VLA observations for the NNW component show a spectral steepening with increasing frequency, which is also indicated in the 1.5–6.2 GHz and 6.2–8.5 GHz spectral index maps in Figure 9. Interpreting NNW as the approaching hotspot of the quasar, we expect it to have a flatter spectrum than the receding one (e.g., Dennett-Thorpe et al. 1997; Ishwara-Chandra & Saikia 2000). Thus a nondetection of the other (receding) hotspot can be explained by the fact that the spectral index difference is enhanced resulting from the small inclination angle of J0909+0354 (e.g., Dennett-Thorpe et al. 1997; Ishwara-Chandra & Saikia 2000). We can apply the formula

$$\begin{eqnarray}&&K={\left(\displaystyle \frac{1+\beta \cos \theta }{1-\beta \cos \theta }\right)}^{2-\alpha },\end{eqnarray} \tag{ 4 }$$

where K is the ratio of the flux densities of the approaching and receding sides of the jet, β = vc⁻¹ is the speed of the jet, θ is the jet inclination angle with respect to the line of sight of the observer, and α is the spectral index of the hotspot. With β = 0.3 (e.g., from Dennett-Thorpe et al. 1997), θ = 23° (upper limit calculated for J0909+0354), and single-epoch NNW spectral index ${\alpha }_{7.77,\mathrm{NNW}}^{4.62}=-1.03$, the flux density ratio for the approaching (NNW) and receding hotspots of the kiloparsec-scale jet is K = 5.42. Flux densities for the hotspot on the receding side are then expected to be approximately S_1.5,r = 4 mJy, S_6.2,r = 1 mJy, and S_8.5,r = 0.5 mJy, for the 1.5, 6.2, and 8.5 GHz frequency bands, respectively. These values are in the same order of magnitude as the 3σ rms noise of the VLA clean maps (Figure 2), which can naturally explain the nondetection of the hotspot on the receding side of the kiloparsec-scale jet by the array. Considering higher jet speeds or lower jet viewing angles results in even lower flux density estimates. The spectral steepening of NNW with increasing frequency implies radiation losses due to spectral aging (e.g., Krolik & Chen 1991; Blundell et al. 1999; Ishwara-Chandra & Saikia 2000; Vaddi et al. 2019).

Zoom In Zoom Out Reset image size

There is an apparent contradiction between identifying NNW as a hotspot and the proposed spine–sheath feature observed on parsec scales (Section 5.1). This can be resolved by taking into account the different timescales of the parsec- and kiloparsec-scale structures. Although the injection to hotspots is expected to originate from a continuous supply of particles, there is a large physical distance between the inner jet (∼10 pc) and the kiloparsec-scale NNW component (∼17 kpc). The 17 kpc projected distance of NNW with respect to the core component translates to a 43.5 kpc length considering the upper limit of ${\theta }_{\max }=23^\circ$ on the jet inclination angle. The presently observed hotspot-like characteristics are not expected to be affected by the more recent state of the parsec-scale structure. Assuming the speed of the jet fueling NNW is 0.3c, it takes ∼0.5 Myr in the rest frame of the quasar for the newly developed changes in the parsec-scale jet to propagate to the kiloparsec-scale structure. The hotspot scenario is further supported by the X-ray detection of enhanced emission at the NNW region, and a jet-like feature connecting it to the core (Section 4.1). Concluding the discussion above, identification of NNW as a hotspot cannot be excluded.

From the survey of Civano et al. (2016), the post facto probability of such an unrelated strong source being within 7'' of our target is 0.03%. However, with more than 300,000 distinct X-ray sources in the second Chandra source catalog (Evans et al. 2020), such a probability is not evidence of association with the quasar J0909+0354. Indeed, there is a faint source in the Pan-STARRS images coincident with the NE X-ray source. We estimate g, r, i, z, and y magnitudes of 23.89, 22.74, 22.73, 22.63, and 21.50, respectively, from the prescription of Waters et al. (2020). An optical object is also detected at the position of the NE feature in the g, r, and z bands in the Dark Energy Camera Legacy Survey ²⁴ (DECaLS, Dey et al. 2019). We note that, contrary to NE, the NNW component has no optical counterpart. With no radio emission at or in the direction of the NE feature, there is no evidence to associate it with J0909+0354, and we will assume here that it is a foreground or background object.

Mid-infrared Wide-field Infrared Survey Explorer (WISE) maps (from the AllWISE data release, Wright et al. 2010; Mainzer et al. 2011) centered at the position of J0909+0354 at 3.4 μm (W1) and 4.6 μm (W2) show an extended emission surrounding the core component and NE (Figure 10; DN-to-magnitude conversion was carried out as described by Cutri et al. 2011). The emission can be traced up to 15–20''(∼130–160 kpc projected linear distance) with respect to the position of the core. We note that the extended emission can be a blend of individual sources located at different cosmological distances but seen in projection, considering the angular resolution of WISE (61 and 64 for W1 and W2, respectively). But a physical connection between the two objects (J0909+0354 and NE) cannot be ruled out based on the data available.

Zoom In Zoom Out Reset image size