Parsec-scale evolution of the gigahertz-peaked spectrum quasar PKS 0858 − 279

ABSTRACT

We conducted multi-epoch, multifrequency parsec-scale studies on the gigahertz-peaked spectrum quasar PKS 0858 − 279 with the Very Long Baseline Array (VLBA). Our observations on 2005 November 26 elucidated a weak core, characterized by an inverted spectrum, and a distinctly bent jet that exhibited a notable bright feature in its Stokes I emission. Through comprehensive analysis of polarization and spectral data, we inferred the formation of a shock wave within this feature, stemming from interactions with a dense cloud in the ambient medium. In this paper, Very Long Baseline Interferometry-Gaia astrometry further reinforces the core identification. With a deep analysis of six additional VLBA epochs spanning from 2007 to 2018, we observed that while the quasar’s parsec-scale structure remained largely consistent, there were discernible flux density changes. These variations strongly imply the recurrent ejection of plasma into the jet. Complementing our VLBA data, RATAN-600 observations of the integrated spectra suggested an interaction between standing and travelling shock waves in 2005. Moreover, our multi-epoch polarization analysis revealed a drastic drop in rotation measure values from 6000 to 1000 rad m⁻² within a single year, attributable to diminishing magnetic fields and particle density in an external cloud. This change is likely instigated by a shock in the cloud, triggered by the cloud’s interaction with the jet, subsequently prompting its expansion. Notably, we also observed a significant change in the magnetic field direction of the jet, from being perpendicular post its observed bend to being perpendicular prior to the bend – an alteration possibly induced by the dynamics of shock waves.

1 INTRODUCTION

The evolution of active galactic nuclei (AGNs) is pivotal to understanding the formation of galaxies and studying the environment in which they evolve. Gigahertz-peaked sources (GPS) are key candidates believed to be in early developmental phases. Probing these sources can deepen our comprehension of galaxies during their initial stages of growth. While the general framework of AGN evolution encompasses vast cosmic time-scales, our specific focus in this article is on a time frame of roughly 20 yr. This allows us to highlight short-term dynamics that may offer insights into broader evolutionary patterns. Several studies have cast light upon GPS sources, unravelling their temporal properties (e.g. Carvalho 1998; O’Dea 1998; Lister 2003; Stawarz et al. 2008; O’Dea & Saikia 2021). These sources manifest an inverted spectrum below the turnover frequency, attributed to synchrotron self-absorption or free–free absorption. By scrutinizing the short-term variations in their properties and spectra, we aim to address many challenges surrounding GPS sources, such as classification issues, their interaction dynamics with the surrounding medium, and the inherent nature of absorption mechanisms.

Very Long Baseline Interferometry (VLBI) offers an unparallelled perspective, allowing us to delve deep into the processes occurring near the central engines of quasars. VLBI emphasizes explosive events and the intricate birth and formation dynamics of jets. These jets can interact dynamically with their surrounding environment. The capacity of VLBI to distinguish moving jet components offers a detailed glimpse into the kinematics of relativistically moving plasma (e.g. Lister et al. 2009). Typically bright knots, whose origin remains an open question, are commonly observed. In some cases, they have superluminal apparent speeds far from the core, but they can also be quasi-stationary (e.g. Lobanov 2011; Hervet, Boisson & Sol 2016). Several models have been proposed to explain the properties of the moving plasma (e.g. Blandford & Königl 1979; Marscher & Gear 1985). Furthermore, numerical simulations are used to investigate how jets evolve and to exhibit these bright features (e.g. Gómez et al. 1997; Komissarov 1999; Agudo et al. 2001; Fromm et al. 2016; Komissarov & Porth 2021).

In a previous article (Kosogorov et al. 2022), which we will refer to as Paper I, we investigated the properties of quasar PKS 0858 − 279 in detail. The source captured our interest due to its high amplitude variability over time-scales of several months from the RATAN-600 data, while displaying extended structure in the early VLBI images. We resolved both the core and a prominent bright jet feature. Analysing the spectral and polarization properties, we found out that a shock wave was formed in that feature. This shock formation may explain the notably stronger magnetic field in the jet feature, estimated to be ∼1 G, compared to the core’s magnetic field of about 0.1 G. Furthermore, from the multifrequency light curves, we discerned a variability Doppler factor of around 2. This estimation method leverages the Doppler boosting of flux density variations in the quasar’s jet, representing one approach to find the Doppler boosting factor (see eq. 3 in Paper I). In this work, our aim is to understand the changes in the properties of the source over time, pinpoint the source of variability, ascertain the kinematics of the jet plasma, and track the evolution of the brightest jet feature. To achieve these goals, we examine the source evolution using multi-epoch polarization VLBI observations.

This paper is structured as follows: in Section 2, we provide the observation details and the data reduction process. Section 3 showcases the joint analysis of optical Gaia positional observations compared with the VLBI ones. In Section 4, we present the results of the VLBA and RATAN-600 data analysis, including Stokes I images, modelling of the source structure, kinematics, spectral properties, and variability and polarization properties. In Section 5, we discuss our results in light of possible scenarios of the changes which took place during the time-span of our observations. Finally, Section 6 sums up all the findings. Through the paper, the model of flat lambda-cold dark matter cosmology is used with |$H_0 = 68\, \mathrm{km\, s}^{-1}$|⁠, Ω_M  = 0.3, and |$\Omega _\Lambda =0.7$| (Planck Collaboration XIII 2016). The position angles are measured from north through east.

2 OBSERVATIONS AND DATA REDUCTION

2.1 VLBA observations

As part of a dedicated campaign to study the quasar, it was observed with the Very Long Baseline Array (VLBA) of the National Radio Astronomy Observatory (NRAO) on 2005 November 26, and 2007 February 09, May 23, and August 03. The project codes are BK128, BK133a, BK133b, and BK133c. Details of the observation and the reduction of the BK128 project data can be found in Paper I. Here, we describe the BK133 project data. The observations were made in a dual-polarization mode. The observing bands were centred at 1.5, 2.3, 4.8, 8.3, 15.4, and 22.2 GHz. Each band comprised four 8 MHz-wide frequency channels (IFs) per polarization. We combined all IFs at all frequencies in the consequent analysis, except for the analysis presented in Section 4.6. In the polarization analysis, the four lowest frequency bands were split into two sub-bands, two IFs per sub-band, centred at 1.41, 1.66, 2.28, 2.39, 4.64, 5.04, 8.15, and 8.50 GHz. The data at each VLBA station were recorded with two-bit sampling and a total rate of 256 Mbps in right and left circular polarization.

The observed calibrators consisted of the fringe check, cross-polarization, and bandpass calibrator 1226 + 023 (3C 273) and the leakage-term calibrator 0919 − 260. Apart from that, the electric vector position angle (EVPA) calibrators included 1328 + 307 (3C 286) during all three BK133 epochs, 0851 + 202 (OJ 287) in February and 0923 + 392 in May and August.

The initial data processing steps were performed in the Astronomical Image Processing System (aips) package (Greisen 2003) for each sub-band independently, paying special attention to accurate amplitude calibration. They comprised typical VLBI data processing procedures which can be found in the aips cookbook. For the first phase self-calibration step, the CALIB task in aips was utilized. We used the difmap software package (Shepherd 1997) to iteratively do phase and amplitude self-calibration as well as obtain clean images. The resulting visibility amplitude calibration accuracy after a priori and self-calibration steps was estimated to be approximately 10 per cent for 22 GHz and around 5 per cent for frequencies in the 1.4–15 GHz range. These estimations agree well with the calibration accuracy described in the VLBA observational status summary.¹ Typical processing of the 8 GHz astrometric VLBI sessions (Table 1) might result in somewhat lower accuracy of amplitude calibration due to the IFs being widespread across the band and shorter on source observing time.

Apart from using the traditional clean algorithm (Högbom 1974), we used regularized maximum-likelihood methods, which are implemented in the eht-imaging package (e.g. Chael et al. 2018), to check the fidelity of our images. In Fig. 1, we present the images at the highest frequency of our observations of 22 GHz in 2005 reconstructed with these two approaches. When reconstructed with the same beam parameters, the results were consistent, with only minor differences noted in the southeast emission.

We calibrated the EVPA of 0858 − 279 using 1328 + 307 for 1.5 and 2.2 GHz bands and 0851 + 202 and 0923 + 392 for all other frequencies. The calibration procedure was performed in the same manner as described in the previous analysis. We made use of the calibrators absolute EVPA information from the VLA (Very Large Array) Monitoring Program² (Taylor & Myers 2000) and the MOJAVE³ (Monitoring Of Jets in Active galactic nuclei with VLBA experiments) observations.

Additionally, we accessed VLBI images of the source at 8.3 GHz from three other epochs: 2017 May 27, and 2018 April 08 and July 21, which are publicly available from the Astrogeo VLBI fits image database.⁴ The summary of all VLBI observations used in this paper can be found in Table 1.

2.2 RATAN-600 observations

Furthermore, we observed the source at 5–7 frequencies from 1 to 22 GHz with the RATAN-600 radio telescope of the Special Astrophysical Observatory (e.g. Parijskij 1993) since 1997, conducting 2–4 observing sessions annually. Within each session, the continuum spectra for the source were typically measured 2–4 times and averaged. We utilized the data from 1997 to the present date, extending the previously used time period of the first 13 yr. The details on the observations, calibration, and presented errors can be found in the corresponding part of Paper I.

We compared the integrated flux densities from the VLBI images with RATAN-600 flux densities for the nearest dates to their VLBI observation. They turned out to be in good agreement with all epochs but 2017 May 27. This specific epoch demonstrated that the flux density is 30 per cent lower than the corresponding one from RATAN-600. To account for that, we applied a correction in difmap to make the integrated flux density consistent. We note that, for some data sets within the Astrogeo database, amplitude calibration might not be the primary focus, especially for astrometric/geodetic-type observations with wide-spread IFs. At the same time, every observing RATAN-600 session utilizes observations of several calibrators. Given the meticulous calibration procedures associated with RATAN-600 observations, we suggest that it is the VLBI flux density from this specific data set in the Astrogeo database that required correction.

2.3 Gaia and VLBI astrometry observations

This work used the Gaia Data Release 3 (Gaia DR3, Gaia Collaboration 2016, 2023). This catalogue provides extremely accurate measurements of optical positions of various sources, including AGNs. Similar submilliarcsecond accuracy can be obtained with the VLBI observations. In further analysis, we used the VLBI-based Radio Fundamental Catalog (RFC)⁵ to cross-identify and study 0858 − 279.

3 VLBI-GAIA ASTROMETRY ANALYSIS

In most cases, we expect Gaia positions to be close to the jet base, within the VLBI core. At the same time, the previous analyses of the VLBI-to-Gaia offsets showed that there are about 10 per cent of AGN with significant differences between their optical and VLBI coordinates (e.g. Kovalev, Petrov & Plavin 2017; Petrov & Kovalev 2017; Petrov, Kovalev & Plavin 2019; Plavin, Kovalev & Petrov 2019b). These offsets can vary and be as large as tens of mas. Among AGN with significant offsets, the majority of sources display Gaia positions further downstream along the jet direction compared to the VLBI coordinates. Such behaviour is expected when strong and extended parsec-scale optical features are found in jets. On the other hand, some sources exhibit offsets in the opposite direction. This can be caused by the core-shift effect or the dominant impact of the accretion disc.

Since 0858 − 279 displays an unusual VLBI parsec-scale radio structure characterized by a weak core and a pronounced jet feature, there was significant interest in comparing its optical coordinates with the VLBI ones. Under the assumption that the RFC VLBI coordinates align with the dominant jet feature J₀, the RFC coordinates of the source correspond to this position. This presumption is based on the consistent observation of the feature for over a decade (further elaborated in Section 4.1). The Gaia and VLBI coordinates can be found in Table 2. The Gaia position does not match the dominant jet feature. The offset between Gaia and VLBI is 2.4 ± 0.2 mas. By understanding this offset and comparing it with the RFC catalogue’s indication for J₀, we positioned the Gaia coordinates on our map. As a result, the Gaia position appears close to the core region, suggesting its association with the accretion disc. The Gaia and the VLBI core positions matched with good accuracy, within the error limits (refer to Fig. 2). Such a close alignment emphasizes that radiation from the accretion disc and the beginning of the jet dominates at optical wavelengths.

4 VLBA AND RATAN-600 DATA ANALYSIS RESULTS

4.1 Stokes I maps and modelling of the quasar structure

The general structure of the source in 2007, 2017, and 2018 remained consistent with that of 2005, with no new features appearing in its Stokes I emission. One can find the images at all frequencies from the epoch 2005 November 26 in fig. 4 in Paper I. In Fig. 3, we display the stacked Stokes I images of the quasar from all available epochs at 8 and 15 GHz. The resolution of these images was set to an average resolution of all maps used for the reconstruction. We do not show stacked images at 22.2 GHz, since there are only three epochs available for stacking, and the resolution is comparable to the resolution of the 8 GHz images. Such low resolution arises because data from some antennas were corrupted. The problems for different time periods included the absence of data, communication problems, bad weather conditions, and power glitches (see Table 1 for more detail). Additionally, in Table 3, we present the main parameters of all maps at all frequencies.

The stacked image at 15.4 GHz shows an emission in the northeast region relative to the central brightest jet feature. To assess the significance of this radiation, we merged the UV data sets using the DBCON task in aips and subsequently reconstructed the source image. This image did not present the signs of the corresponding emission. Thus, we concluded that the radiation does not match the real physical plasma emission from this region.

In our earlier analysis, we determined that the southwest region was associated with the opaque core, and the prominent central feature was actually a part of the jet. Lower frequencies are not discussed in this analysis due to their limited resolution, which made detailed structural study challenging. To study the structure and its evolution in more detail, we fitted circular Gaussian components to the visibility data in difmap with the model fit independently for different epochs. We continued to add new Gaussian components only when they led to a significant reduction in χ². Thus, this approach enabled us to find the optimal minimum number of components. In Table 4, we give their main parameters at 8.3, 15.4, and 22.2 GHz for all epochs. Typical errors turned out to be 5 per cent–10 per cent. In addition, Fig. 2 shows the results of our modelling at 15.4 GHz. We note that we were able to model the structure with a minimum of four components at 15.4 GHz. Meanwhile, for the 2007, 2017, and 2018 epochs, three components best represented the data at 8.3 GHz. We tried to distinguish between the core component C and the next one along the jet, S, further downstream, by fitting additional delta components (also referred to as clean components) in this region, alongside the Gaussian components. However, the fit did not converge due to low resolution.

A notable observation is that the number of Gaussian components best describing the brightest jet feature transitioned from two (J₀ and J₁ components) to just one (J component) over a year-long period. On the one hand, this might be directly related to the evolution of a shock wave across that specific part of the jet (see Section 5 for more detail). On the other hand, it can be an artificial effect since the fitted components do not always represent real physical features in the jet (e.g. Lister et al. 2009).

4.2 Kinematics

We analysed the variation in distance between the core component C and the jet component J over time, as illustrated in Fig. 4, utilizing the 8.3 GHz data which encompasses all seven epochs. The results show that the distance increased with the jet feature apparent speed of (1.9 ± 0.5)c assuming that the core component is stationary. Such behaviour might hint at how the component slowly carves its way through the external medium. If we assume that the jet is observed at critical angle (e.g. Homan et al. 2021), the relativistic Doppler factor can be estimated as |$\delta = \sqrt{1+\beta ^2} \sim \beta \sim 2$|⁠. This is consistent with the results obtained in Paper I, where the variability Doppler factor was estimated. This result, however, does not conclusively rule out the possibility of a stationary feature where the apparent motion is caused by changing opacity effects in the source.

4.3 Variability of the parsec-scale jet structure

To study the variability properties of the quasar, we need to align the images from different epochs. That was not a trivial process since the absolute positions of sources in the sky are lost during the self-calibration process. Apart from that, different components can evolve between multi-epoch observations, and new ones can emerge due to the ejection of plasma material. Finally, variable synchrotron self-absorption can have an impact on visible positions. As a consequence, this ambiguity in positions and the number of components complicates the process of component identification across epochs. However, modelling of the quasar structure demonstrated that the core component C was consistently visible across all epochs (see also spectral index maps in Section 4.5).

We tried several approaches to align the images at 15.4 GHz first between 2005 November 26 and 2007 February 09. They included aligning based on the position of the core component as well as by the positions of the J and J₀ or J and J₁ components. After visually inspecting the differences between the Stokes I images, we determined that aligning based on the core component yielded the best results. In Fig. 5, we present the difference between the Stokes I images of the quasar at 15.4 GHz, as reconstructed with the same pixel sizes, uv-coverage, and an epoch averaged beam. At the same time, other alignment methods produced results with minor artefacts in the image, but these did not significantly deviate from the results obtained by aligning with the core component.

The most drastic changes in the 2005 and 2007 epochs occurred in the brightest jet feature. The changes in the core’s flux density and in the southeast component were almost negligible. Additionally, the changes in the emission between three epochs in 2007 turned out to be minor, showing that the flux density was basically unaltered both in the core region and in the jet. Therefore, we do not show these images.

We investigated the variation in radiation between 2007, 2017, and 2018 by aligning the Stokes I images at 8.3 GHz. The best results for the intensity difference maps were obtained by the jet feature alignment. We present them in Fig. 5. The core alignment yielded similar findings, meaning that in both cases, the most drastic changes were observed mainly in the jet feature. However, the core alignment resulted in a difference image where the region of flux change was not as symmetric relative to the centre of the component as it was in the jet alignment. We note that the flux density dropped in the jet feature J between 2007 and 2017. There was a small flare in 2018 with flux density rising in both core and jet compared to 2017. In accordance with the causality principle, the simultaneous growth in these two regions cannot be linked. The projected distance between the core and the component J is approximately 2 mas, which is about 17 pc. Such a distance makes it improbable that the flare had time to reach the component J. Finally, the core flux started to drop in the same year, while the jet flux density increased slightly. Due to the minimal impact of these changes, we have chosen not to display them. Table 4 shows more detail on the amount of flux density changes in the core and the jet features.

4.4 RATAN-600 single-dish spectrum variability

We conducted a long-term variability analysis of the quasar using RATAN-600 data. Light curves for the three highest frequencies are shown in Fig. 6. We modelled the integrated spectrum of the quasar at each epoch using the formula that describes the spectrum of a synchrotron homogeneous source (e.g. van der Laan 1966; Pacholczyk 1970):

where I_m is the maximum intensity in the spectrum reached at the corresponding frequency ν_m, and τ_m, the optical depth at the turnover frequency ν_m, is approximated as:

In Fig. 7, we illustrate the evolution of these parameters on the (ν_m, S_m)-plane. The errors were derived using the variance–covariance matrix of the fit. We notice that at some epochs the fit did not converge. Therefore, the number of points in Fig. 7 is lower than the total number of epochs in the RATAN-600 observations. We acknowledge that the assumption of homogeneity for the integrated spectrum is somewhat simplistic, given that the spectral distribution includes contributions from various jet regions. However, we utilized this conservative approach in further analysis instead of the generalized one with a varying index of optically thick radiation (see eq. 7 in Paper I) to keep the number of free parameters to a minimum. Apart from that, the generalized approach showed the same evolutionary trends for the evolution of I_m and ν_m.

4.5 Spectral index maps

Throughout this paper, we define the spectral index α as S ∝ ν^+α, where S is the flux density observed at frequency ν. This convention is used for all spectral index calculations and discussions. Fig. 8 displays the reconstructed spectral index and the corresponding error maps for the 15–22 GHz range during two epochs: 2005 November 26 and 2007 February 09, both with identical angular resolution. We aligned the images at different frequencies by using the optically thin jet components (e.g. Kovalev et al. 2008; Sokolovsky et al. 2011; Fromm et al. 2013) as a reference similar to the 2005 epoch analysis, because their absolute positions do not depend on frequency. The derived spectral index maps exhibit a similar structure with a clearly visible opaque core in the southwest region. Additionally, the values increase in the core region compared with the 2005 November 26 epoch. This could be partially caused by boundary effects, which refer to the increased uncertainty in areas where the signal-to-noise ratio is lower. To test if this increase is significant, we calculate spectral indices for the corresponding core components taken from Table 4. The core spectral index value in 2005 turns out to be 0.4 ± 0.4, while in 2007, it is 1.1 ± 0.4. We conclude, therefore, that the increase is significant. Apart from that, the spectral indices also increase in the brightest jet feature. This can be related to the opacity properties, which are critical in the subsequent analysis. This is considered carefully when studying the polarization properties of the source in Section 4.6. The spectral index error is calculated in the same manner as in Paper I. Owing to the low resolution at 22 GHz and the varying opacity effects at 8 GHz for the 2007 epochs (Section 2.1), we did not measure the core shift (Blandford & Königl 1979; Konigl 1981; Lobanov 1998; Pushkarev et al. 2012; Plavin et al. 2019a), and as a result, magnetic field estimates for the core region during these epochs were not determined.

4.6 Polarization properties

We performed polarization analysis for the three epochs in 2007. First, we reconstructed fractional polarization images. They do not demonstrate any significant changes in the structure compared to the 2005 epoch, but, over time, the radiation becomes systematically less polarized at the highest frequencies. Table 5 shows the degree of linear polarization for the dominant component J. Second, some of the EVPA images showed clear presence of opacity effects. For example, in Fig. 9, a distinct region in the centre of the 8.5 GHz map has EVPAs perpendicular to those in the surrounding areas. This can be explained in terms of an emission regime which becomes optically thick at around 8 GHz in the centre region, while being optically thin in the outer parts. The evidence for opacity in this central part is a low degree of polarization and high spectral index. This effect is stable across all epochs in 2007 and the 2005 epoch, with EVPAs in the upper region of the map being perpendicular to EVPAs in the lower part of the images. Apart from that, when considering this effect and applying a 90° correction to this pronounced region (see the right part of fig. 9), the reconstructed magnetic field direction and the RM maps turned out to be consistent between different frequencies and epochs. Additionally, a large region in the centre of the 8.2 GHz EVPA map showed no signal above the 3σ_P level. This is likely caused by changing opacity and a potentially strong increase in absorption at this specific frequency.

Furthermore, we reconstructed the RM maps for frequency ranges 1.4–2.4, 2.3–5.0, 8.2–15.4, and 15.4–22.2 GHz for two epochs 2007 February 09 and May 23. For the 2007 August 03 epoch, instead of taking the 8.2–15.4 and 15.4–22.2 GHz ranges, we used the range 8.5–15.4 GHz, since the 22.2 GHz data were not optimal for analysis (Table 1), and polarization signal in the 8.2 GHz EVPA map is small. The RM values were calculated in the same manner as in the analysis of the 2005 November 26 epoch. We then derived the magnetic field direction images, considering the radiation regime and applying a 90° correction when necessary. In Fig. 10, we show the rotation measure (RM) maps and the magnetic field directions for the highest frequency intervals from epochs in 2007. Apart from that, in Table 6, we present the corresponding RM and magnetic field direction values. The results at other lower frequencies turned out to be the same within the error limits. However, we do not show them since the uncertainties in the direction of the EVPAs at low frequencies were drastically larger. In addition, the results are in good agreement between different epochs in 2007.

We used several criteria to analyse the RM gradients in the reconstructed maps. They are described by Hovatta et al. (2012) and include the jet width, which should be preferably more than two beams wide based on the simulations performed to investigate how large spurious gradients can arise due to noise and finite size of the restoring beam. In our study, the size of the jet width is approximately two beams, indicating we have enough resolution. Apart from that, gradients should exceed the 3σ level, where σ is the largest RM error. Similar criteria can be found in, for example, Taylor & Zavala (2010) or Kravchenko, Kovalev & Sokolovsky (2017). For our data, the gradients turned out to be insignificant for all epochs and frequency ranges.

5 DISCUSSION

5.1 Kinematics

The structure of the quasar revealed no significant changes over a period of more than 10 yr. The brightest jet feature was observed during all epochs from 2005 to 2018. Fig. 8 shows that the spectral index increases in this region, which must correspond to an energy gain in the emitting particles at the observed interval. The energy which goes into particle acceleration is probably being generated at the expense of kinetic energy and, taking into account that the feature is steady, we could interpret it as a strong interaction between the jet flow and a dense interstellar medium (ISM) region. The impact from a cloud could trigger a shock, which would fit in the picture of the apparent observed bend in this region where the J component is located. In addition, we observe the component S between the core C and the brightest jet feature, J₁ in 2005 and J in 2007, at the highest frequencies (see Fig. 2). The low resolution during other epochs did not allow us to detect it in 2017 and 2018. This component was located at a distance of ∼0.7 mas from the core in 2005 and 2007 May and August. According to its properties, we suggest that this component could be a stationary recollimation shock wave. Its position shifted a little further downstream in 2007 February with the distance from the core becoming ∼1 mas. Such behaviour could be a consequence of the standing shock being dragged downstream and then returning to its initial position upstream (e.g. Gómez et al. 1997; Fromm et al. 2016). However, changes in the opacity of the core region could also affect the relative positions of the components.

The two components J₀ and J₁ observed in 2005 (see Fig. 2) could hint that there were two shock waves in 2005, one of which was stationary and the other was travelling. In this scenario, J₁ in 2005 could be cross-identified with the J component in 2007 by the same distance to the core – assuming that the opacity effects did not affect the core positions. This would suggest that J₀ was the travelling component. However, it is not clear why the J₁ flux density in 2005 was not consistent with the flux densities of J in 2007 (see Table 4). One feasible explanation could be related to non-linear effects of interacting shock waves, which is the scenario we favour (see also further discussion of the spectral peak values evolution).

If the J and J₀ components were the same ones, the apparent velocity would be (8.9 ± 1.9)c, where c is the speed of light, which again made use of an assumption that the core position was stationary. In principle, the core position can change due to significant nuclear flares (e.g. Plavin et al. 2019a). The time-scale of such changes may align with our observed time-scales, on the order of a year. However, considering the relatively small variation in flux density in our core component, we propose that the assumption of a stationary core position in this context is justified.

Although superluminal motion is often found in the jets (e.g. Kellermann et al. 2003; Lister et al. 2016), in this case, the distance between the core and the component J would have decreased, implying a backward motion of the brightest jet feature. Such behaviour would be exceptionally peculiar at these velocities, as slower motions are typically observed in sources exhibiting inward movements (e.g. Lister et al. 2013, 2016, 2021). At the same time, it would not be consistent with the results obtained in Section 4.2. Nevertheless, significant superluminal component motions towards the core can sometimes occur (e.g. Sawada-Satoh et al. 2015). They are mostly explained by virtue of highly curved jets geometry or emerging unresolved components in the core. In the context of a highly curved geometry, Sawada-Satoh et al. suggest that if the apparent inward motion is due to a highly curved trajectory, the flux density of the corresponding moving component would exhibit time variability due to a variation in the geometrical-origin Doppler factor. However, the radio light curve in their study reveals a constant flux density during the inward motion, rendering this scenario less likely. In our observation, significant changes in flux density occur when J is associated with J₀. Simultaneously, in this situation, the jet should curve back towards the line of sight (e.g. Kellermann et al. 2003), making emission in the southeast direction unlikely.

Apart from that, clear helical patterns are sometimes observed in AGN jets (e.g. Hong et al. 1999, 2004; Zhao et al. 2011; Perucho et al. 2012; Bruni et al. 2021; Gómez et al. 2022; Fuentes et al. 2023) where a curved jet geometry is noticeable and could potentially explain the inward motion in our scenario. Typically, these instances exhibit drastic changes in the position angle of jet components (e.g. Lister et al. 2021). Contrarily, in our case, the position angle remains relatively stable over extended periods, leading us to favour alternative scenarios. Moreover, the field is usually perpendicular to the direction of propagation, whereas in our case it is aligned with the hypothetical interaction region. These observations suggest that this scenario may not fully apply for our data.

Finally, the last option is that the two components in 2005 could not represent any real physical features, as discussed in Section 4.1. The numerical simulations of the interaction between a stationary shock and a travelling one (e.g. Fromm et al. 2016) and further high-frequency VLBI observations could potentially help us distinguish between different scenarios.

5.2 Spectral changes and flux density variability in the core and jet feature

The observed changes in the flux density in the core and the jet feature suggested that occasionally new material was ejected in the jet. Apart from the visible effects on the VLBI variability images, this could also be a natural explanation for the variation of flux densities in the integrated spectrum (see Fig. 6).

In principle, scattering effects could impact the flux variability of the quasar. As elucidated in Paper I, an analysis of the size–frequency dependence for the quasar was conducted, primarily focusing on the core at high frequencies (15 and 22 GHz). The findings revealed a good agreement with the conical synchrotron jet model, where core size depends on frequency as ∝ ν⁻¹ (Blandford & Königl 1979), suggesting that scattering effects were not dominant at these frequencies. In this paper, upon extending the analysis to other epochs, a deviation from the conical synchrotron jet model was observed, hinting at the potential presence and influence of scattering. While extreme scattering events can occasionally modify the flux by a factor of a few, exhibiting a distinct U-shaped pattern (e.g. Vedantham 2018), they generally do not contribute significantly to the overall observed variability. Furthermore, beyond these extreme events, scattering does not notably alter the flux, as per typical scattering behaviour. Apart from these considerations, strong scattering is deemed not very likely given the sky coordinates of the quasar (e.g. Pushkarev & Kovalev 2015; Koryukova et al. 2022), leading us to maintain the assumption that scattering does not significantly impact our main results.

The (ν_m − S_m)-plane shows two distinct stages (see Fig. 7). The first stage, lasting from 1997 to ∼2010, displays a general trend of both the flux density and the peak frequency dropping, indicative of an adiabatic stage in the classical shock-in-jet model (Marscher & Gear 1985). The subsequent stage presents a scenario where the flux density remains roughly constant, but the peak frequency increases. This latter behaviour is not typical for the classical model, given the recognized loss mechanisms such as adiabatic, synchrotron, or inverse Compton losses. Such observations suggest a more intricate picture, hinting that the simple shock-in-jet model may not fully encapsulate the dynamics at play in the jet. Moreover, the high errors of ν_m and S_m do not allow us to distinguish between other models of shock–shock interactions and extract any parameters. In Fromm et al. (2016), the authors show that the evolution of the (ν_m − S_m)-plane depends on whether the spectral peak is produced at the shock-perturbed region in the case of a conical expanding jet or in the case of shock–shock interactions. According to these results, the evolution in the aforementioned plane is clockwise in the latter case and counterclockwise in the former (see fig. 11 in that paper). Although the numerical simulations simplify the scenarios, and the spectra that we observe correspond to the whole source rather than the J components, we can derive interesting conclusions from a direct comparison. Indeed, the global evolution in the (ν_m − S_m)-plane (see Fig. 7) is in the counterclockwise direction, which means that the spectral peak is probably produced by the evolution of a perturbation along the conically expanding region. However, in 2004, the peak coordinates displace to the right in the plane (see points in black squares in Fig. 7) and come back to follow the original track along the C-shaped evolution with oscillations to lower peak frequencies. The sudden rightward motion is perhaps due to the interaction between the components J₀ and J₁. It is possible that precisely in 2005, we started to see these shocks separating. The rest of the time-span between 1997 and 2020, the evolution of the spectral peak would be dominated by the expansion of the travelling perturbation along the jet, until the component fades out. At this point, the spectral peak could be produced either by the core or by the bright stationary J component.

5.3 Rotation measure results and implications

We note that the RM values dropped significantly from above 6000 rad m⁻² in 2005 to around 1000 rad m⁻² in 2007. This behaviour is not common for jets but possible for various reasons, and it has been observed before (e.g. Zavala & Taylor 2001). One possibility to explain this could be that properties of the external medium in the screen changed significantly, that is, its density, magnetic field, or the direction of the magnetic field, causing a change in the projected value along the line of sight. In this case, the hypothetical colliding external cloud should be rather small to cause such changes over year-long periods in the observed RM. When we tried to calculate its size, taking into account causality, it turned out to be around 0.3 pc. If projection effects were added, the cloud size would have parsec scales. Nevertheless, both the magnetic field and the density of the interacting regions might have dropped within the span of just over one year separating the observing epochs, causing this change. Given the causality argument, this implies that the region undergoing these changes must be smaller than just over one light-year in size. It is feasible that the hypothetical interaction between the jet and a dense cloud affected the cloud properties and the electron density decreased. The numerical simulations by Bosch-Ramon, Perucho & Barkov (2012) showed that the interaction between a jet and a cloud triggers a shock in the latter, heating it and causing its expansion. The cloud thus becomes less dense. Apart from that, under this interpretation, we can predict that the RM values should keep falling and, at some point, the component J will fade out. Within this scenario, quasar 0858 − 279 represents a good candidate to show shock-ionization lines as those observed in other young sources (e.g. Labiano et al. 2005, 2006; Holt, Tadhunter & Morganti 2009; Mahony et al. 2016; Morganti et al. 2021).

Alternatively, a different part of the jet might have been illuminated, leading us to measure the RM values in that distinct region. However, the distance between the most illuminated parts of the jet should not exceed the distance between the J₀ and J₁ components, which is ∼0.2 mas (or 1.5 pc), quite short for such drastic changes. If we assume that there is one cloud with constant properties and a size of more than 1 pc, this option is highly unlikely.

The direction of the magnetic field is also considerably altered: it became perpendicular to the propagation direction of the jet from the central engine to the observed bend (Fig. 10). This could fit well into the picture of two interacting shock waves, in the sense that the polarized emission could be concentrated at the head of the moving shock. However, the interaction between the jet and the Faraday screen could also be the reason for a change in polarization properties.

6 SUMMARY

The GPS quasar 0858 − 279 initially exhibited several bizarre peculiarities, including strong variability, which could not be explained by its seemingly stable structure in the early VLBI images. In the analysis of the 2005 November 26 epoch, we were able to resolve a typical core-jet structure. We also came to the conclusion that a shock wave is formed in the jet via a strong interaction with the ambient medium. That made the emission from this part brighter than from other parts of the source.

In this paper, we have presented multi-epoch analysis of the quasar with additional six epochs of observations. The brightest jet feature was detected during all epochs covering a range of almost 15 yr. VLBI-Gaia astrometry confirmed the interpretation of the source structure and the identification of the core. The variability of emission demonstrated that most of the changes occurred in the brightest jet feature, yet for some epochs, flux densities varied slightly in the core region. Apart from that, the spectral properties of the source changed within a year-long period. This included flattening of the spectrum in the jet feature – components J – and spectrum becoming more inverted in the core region. Furthermore, we studied the polarization properties for three epochs in 2007. The most noticeable changes were the RM values decreasing from 6000 to 1000 rad m⁻² in one year.

We studied the spectral evolution of the source along the epochs and found that it fits well into the picture of evolving perturbation along the conically expanding jet, presumably between the components C and J. However, we observe rapid changes in peak frequency and flux density around 2005. This is interpreted as a hint of interaction between the evolving perturbation and the bright component J, that is, a possible shock–shock collision, as supported by numerical simulations (Fromm et al. 2016).

The changes in the RM values can be interpreted as being produced by the changes in the screening medium. In particular, this could be induced in the hypothetical ISM cloud by the shock. In this frame, as the cloud is shocked and heated, it expands and becomes more dilute, changing the electron densities and therefore, the measured Faraday rotation. Interestingly, this is also in good agreement with the interpretation of a strong collision between the jet and a dense cloud, which would act as a Faraday screen. In this scenario, we expect the RM to keep falling in time and, at some point, the component J to become dimmer.

In one year period, we also observe a change in the direction of the jet magnetic field, becoming perpendicular to the propagation of the jet before its observed bend. If we take into account that the jet kinematics and spectral evolution could also be interpreted in terms of travelling perturbation interacting with recollimation, standing shock at the location of the component J, the changes in the magnetic field direction could be associated to the reinforcement of the perpendicular components by the travelling shock, after the passage.

Further VLBI observations with improved uv-coverage with polarization analysis, search for hints of shock ionization of the ISM, as indicative of jet–medium interaction, and numerical simulations could shed light on the evolution of this peculiar quasar.

ACKNOWLEDGEMENTS

We thank Nicholas MacDonald for reading the manuscript and useful comments. We thank E. Bazanova for the English language editing of this manuscript.

YYK was supported by the M2FINDERS project which has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 Research and Innovation Programme (grant agreement no. 101018682). MP acknowledges support by the Spanish Ministry of Science through grants PID2019-105510GB-C31, PID2019-107427GB-C33, and from the Generalitat Valenciana through grant PROMETEU/2019/071.

The VLBA is an instrument of the National Radio Astronomy Observatory, a facility of the National Science Foundation operated under the cooperative agreement by Associated Universities, Inc. This work made use of the Swinburne University of Technology software correlator, developed as part of the Australian Major National Research Facilities Programme and operated under license (Deller et al. 2011). The observations with the SAO RAS telescopes (RATAN-600) are supported by the Ministry of Science and Higher Education of the Russian Federation. This research made use of NASA’s Astrophysics Data System Bibliographic Services. This research made use of the data from the MOJAVE database, which is maintained by the MOJAVE team (Lister et al. 2018) and the University of Michigan Radio Astronomy Observatory. This work has made use of the data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

DATA AVAILABILITY

The correlated data underlying this article are available from the public NRAO archive,⁶ project codes are BK128 and BK133. The fully calibrated visibility and the image fits files for the target source are made available online by us in CDS.⁷ We made use of the publicly available Astrogeo VLBI fits image database.⁸

Footnotes

References