© The Authors 2025

1. Introduction

Radio-loud active galactic nuclei (AGNs) contain relativistic jets whose formation and evolution are not yet fully understood (for a review of AGN radio jets, see Blandford et al. 2019). Both the ergosphere (Blandford & Znajek 1977) and the accretion disk (Blandford & Payne 1982) of the central black hole of an AGN are thought to play an important role in the formation and launching of relativistic jets. Astrophysical processes are expected to affect the collimation of a jet as it travels through the surrounding medium. This collimation is described by a width profile, w(d)∝d^k, where w is the jet width, d is the distance from the central engine, and k is a slope whose value describes a jet that is staying cylindrical (k = 0), growing parabolically (k = 0.5), or growing conically (k = 1) (see Hada et al. 2013; Kovalev et al. 2020; Boccardi et al. 2021; Ricci et al. 2022).

Ever since the computational power became available, numerical special relativistic hydrodynamic (SRHD) simulations have been used as a tool to study the parsec and sub-parsec scales of relativistic jets. Initial studies focused on the collimation and expected radio emission of steady-state jets (Komissarov & Falle 1996), as well as on the effect of the ambient medium on their standing shock structure (Gomez et al. 1995). Later studies found that perturbing the jet inlet gave rise to traveling shocks, resulting in bright components, observable at apparent superluminal speeds at radio frequencies (Gómez et al. 1997; Komissarov & Falle 1997). These traveling shocks affect the dynamics of the jet, particularly in their interactions with the standing shocks. These works found that when a traveling shock passes through a standing shock, this standing shock is dragged downstream before going back to its equilibrium position. Furthermore, the initially bright and compact shocks become extended and much less bright upon crossing a standing shock, and their width affects the collimation profile of the jet downstream. For all of this to be studied, they determined that the aberration seen by an observer due to the finite speed of light was absolutely necessary. This was confirmed by later studies, which also discovered that sufficiently strong perturbations will cause trailing shocks (Agudo et al. 2001). These shocks travel behind a main component with a slower velocity, and may not be observable at low radio frequencies due to the large blurring beams (Mimica et al. 2009). Fromm et al. (2016) studies in detail the interaction between traveling shock waves and recollimation shocks using SRHD simulations and radiative transfer calculations. The study focuses on the impact of traveling perturbations on the spectral evolution of AGN jets at small inclination angles, so-called blazars, at gigahertz frequencies, as they induce observational signatures that could be detected in very long baseline interferometry (VLBI) observations.

Due to its large viewing angle of θ > 80° (Baczko et al. 2022), meaning that the jet is closely aligned with the plane of the sky, as well as the fact that both its jet and counterjet are visible, the low-luminosity AGN in the galaxy NGC 1052 is a source of particular interest for studying jet properties. The source is located at a redshift of z = 0.005037 (Jensen et al. 2003) (D ≈ 19.5 Mpc) and is classified as a low-ionization nuclear emission line region (LINER) galaxy due to its nuclear properties (Mayall 1939; Fosbury et al. 1978; Ho et al. 1997). It hosts a supermassive black hole with a mass of ≈10^8.2 M_⊙ (Woo & Urry 2002), though more recent works estimate a greater value of ≈10^9.3 M_⊙ (Kameno et al. 2020).

Observations of the source done at frequencies around 1.5 GHz revealed a double-sided jet structure spanning nearly 3 kpc, with a central, flat-spectrum compact core containing around 85% of the emission flux (Wrobel 1984; Kadler et al. 2004a). Very long baseline interferometry studies estimate a position angle¹ of the source in the plane of the sky between 10° and 25° (Kadler et al. 2004a; Baczko et al. 2016; Nakahara et al. 2020). In this work, we set this position angle to 15°. Images by the Very Long Baseline Array (VLBA) and the Global mm-VLBI Array (GMVA) also show a prominent emission gap between the jet and the counterjet at centimeter wavelengths, decreasing in size with increasing frequency and vanishing for ν ≥ 43 GHz (Vermeulen et al. 2003; Kameno et al. 2003; Kadler et al. 2004b; Baczko et al. 2016). Modeling studies show this sort of gap could be caused by a torus-like structure presenting free-free absorption (Fromm et al. 2018).

Sub-parsec-scale studies of the source conducted at 15 GHz between 1995 and 2001 by the VLBA estimated the outward motions in the jets to have a roughly equal apparent velocity of 0.26 ± 0.04 c (Vermeulen et al. 2003). These velocities are in agreement with a detailed 15 GHz analysis of MOJAVE data between 1995 and 2012 that revealed an average velocity of the components of 0.230 ± 0.011 c (Böck 2013). Baczko et al. (2019) calculate faster mean apparent velocities through VLBA observations done at 43 GHz between 2004 and 2009: 0.529 ± 0.038 c for the eastern jet and 0.343 ± 0.037 c for the western counterjet. An asymmetry was also found in the brightness of the jet and counterjet, starting in 2007, while further studies (see Baczko et al. 2022) find asymmetries between the collimation profiles.

Such apparent asymmetry could arise from a symmetric source due to the effects of an obscuring torus, as well as due to relativistic beaming of the misaligned jet-counterjet system. Baczko et al. (2019) could not find consistent parameters for the intrinsic velocities and viewing angle of all observations, pointing to an intrinsic asymmetry between the jet and the counterjet. While the majority of models for relativistic jet launching result in symmetric outflows, several works have investigated intrinsically asymmetric launching. One such proposed setup consists of the torus around a rotating black hole being seeded with a multiloop magnetic field configuration (see Parfrey et al. 2015; Mahlmann et al. 2020). Nathanail et al. (2020), Chashkina et al. (2021), and Jiang et al. (2023) ran general relativistic magnetohydrodynamic simulations of such setup, and found that an alternating polarity in the magnetic loops results in intermittent jets. This effect may be severe enough that the launching of the counterjet is inhibited. Fendt & Sheikhnezami (2013) studied several alternative methods of asymmetric jet launching. They found that a sophisticated prescription of magnetic diffusivity in the torus, accounting for position and local sound speed, resulted in the longest-lasting asymmetry between the jet and the counterjet.

The aim of this work is to study the possibility and extent of inducing perceived asymmetries on a symmetric jet-counterjet system at parsec and sub-parsec scales. We carried out a 2D-axisymmetric SRHD simulation of a mirrored jet in an ambient medium with a decreasing pressure gradient. Once the jets reached a steady state, we injected perturbations into the jet nozzle and had them evolve alongside the jet spine. We then performed radiative transfer calculations on the hydrodynamic results, incorporating both thermal absorption by a torus-like structure and nonthermal emission from the jets. To understand if an asymmetry between both jets could be caused by observational effects, we computed the effect that a finite speed of light has on the emission maps. In addition, to account for the effects of limited resolution and observing time on observational data, we generated synthetic images that recreate real observing conditions. We have built on the work carried out by Fromm et al. (2019), incorporating into our method the injection of perturbations into the jet nozzle as well as the multifrequency functionalities of the imaging algorithm eht-im (Chael et al. 2023).

To scale our results to NGC 1052, we used the redshift-independent distance of D = 19.23 ± 0.14 Mpc, which translates to a linear scale of 0.093 pc/mas (Baczko 2003). This paper is organized as follows. In Sect. 2, we introduce our numerical setup, breaking it down into the jet-torus setup and SRHD simulation in Sect. 2.1 and the radiative transfer in Sect. 2.2. In Sect. 3, we explain our method of generating synthetic images that mimic observational conditions. In Sect. 4, we analyze the final images, and we discuss the results in Sect. 5.

2. Numerical simulations

2.1. Numerical simulations

We used the high-resolution shock-capturing code Ratpenat (for more details, see Perucho et al. 2010) to compute 2D SRHD jet simulations in cylindrical coordinates (r, z), on which we tested our synthetic imaging and analysis pipeline. The setup of these simulations can be found in Fromm et al. (2018, 2019). During the simulation, we set the speed of light to c = 1 and the units of length, velocity, density, and pressure as R_j, c, ρ_a, 0, and ρ_a, 0c², respectively. R_j is the jet radius and ρ_a, 0 is the initial density of the ambient medium. For the sake of completeness, we repeated the basic setup below. Our jet is embedded in an isothermal ambient medium with decreasing pressure following:

$$\begin{array}{c}\hfill p_a(z)=\frac{p_{b,0}}{d_k}{[1+{\left(\frac{z}{z_c}\right)}^n]}^{-m/n}.\end{array}$$(1)

In the equation above, p_a is the ambient pressure, p_b, 0 is the jet pressure, and d_k is the over-pressure of the jet with respect to the ambient at injection. Furthermore, z_c is the core distance and the exponents n and m control the decay of the ambient pressure with distance, z. We used z_c = 10 R_j, n = 1.5, and m = 2. The density of the ambient medium at z = 0 is ρ_a, 0 = 1.0 and it decreases with distance. The jet is characterized by its pressure, p_b, 0 = 0.003, density, ρ_b, 0 = 0.04, velocity, v_b, 0 = 0.5 c at injection indicated by the subscript ‘0’, and the overpressure, d_k = 2.5. An overpressured jet in a decreasing-pressure ambient medium results in a jet with a series of recollimation shocks, which are identified as the bright knots observed by VLBI techniques (e.g. Gómez et al. 1997; Fromm et al. 2018). Finally, we assumed an ideal equation of state using an adiabatic index of $\widehat{\gamma }=13/9$, the expected one for relativistic electrons and cold protons at the same temperature.

To mimic an obscuring torus, we included a flared disk characterized by its inner and outer radii, R_in, out, angular thickness, Θ, density ρ(ρ_Rin, r, θ), and temperature profile, T(T_Rin, r, θ) (for details see Fromm et al. 2018). In Table 1, we list all parameters used for the simulation.

The numerical grid covers 80 × 100 R_j, using 4 cells per jet radius at the injection nozzle. The initial overpressure leads to a pinching jet and triggers the formation of recollimation shocks (see, e.g., Fromm et al. 2016). The steady state jet structure was obtained after roughly five longitudinal light crossing times and is presented in Fig. 1. The bottom panel of Fig. 1 displays the density evolution along the jet axis (dashed white line in the top panel). As can be seen from the local density maxima, there are several recollimation shocks located at x = ±7, ±21, and ±42 R_j. The distance between the recollimation shocks increases, while their density jumps decrease with distance. This behavior is expected from relativistic jets evolving in an atmosphere of decreasing pressure (see, e.g., Fromm et al. 2018).

Fig. 1. Distribution of the density of the initial steady-state jet model in a slice through the midplane (y = 0) of the jet. The jet is seen under a viewing angle of ϑ = 80°. The top panel shows the 2D density distribution and the bottom panel a density evolution along the jet axis (dashed white line).

In order to study the jet dynamics, including the propagation of shocks, we injected periodic density and pressure variations into the jet nozzle with a maximum amplitude of four times the injection values. These perturbations lead to the formation of traveling shock waves. We generated snapshots from the SRHD simulations every t = 0.25 R_j/c time steps² for later use during the radiative transfer calculations.

2.2. Radiative transfer

We used the produced SRHD snapshots to compute the jets’ radiative output, assuming nonthermal synchrotron emission from the jet (including self-absorption) and free-free absorption in the torus. In addition, we took the finite speed of light into account; in other words, a slow-light approach.

As we were not evolving a nonthermal particle distribution within our numerical models, we instead related them to the density and energy of the thermal particles. We assumed a power-law distribution of nonthermal particles:

$$\begin{array}{c}\hfill n({\gamma }_e)=\{\begin{array}{cc}{n}_0{\left(\frac{{\gamma }_e}{{\gamma }_{e,\min }}\right)}^{-s}\hfill & \text{if}{\gamma }_{\mathrm{e},\min }\le {\gamma }_e\le {\gamma }_{\mathrm{e},\max }\hfill \\ 0\hfill & \text{otherwise}\hfill \end{array},\end{array}$$(2)

where γ_e is the electron Lorentz factor, n₀ is the normalization coefficient, and s is the spectral slope. The normalization coefficient and the boundaries of the particle distribution, γ_e, min/max, can be obtained assuming that the nonthermal particles are a fraction, ζ_e, of the thermal particles, containing a fraction, ϵ_e, of their energy. Since we were not evolving the magnetic field, we assumed a fraction, ϵ_b, of the equipartition magnetic field during the emission calculations (for details, see Fromm et al. 2019). Given that the SRHD simulations are scale-free, we scaled them to centimeter-gram-second units by providing a jet radius, R_j, in centimeters and an ambient density, ρ_a, in grams per cubic centimeter. Moreover, we assumed a viewing angle with respect to the line of sight of ϑ = 80° and a redshift of z = 0.005. In Table 2, we list the parameters used to scale our model, chosen to match the total flux of the reference source NGC 1052 at the central frequency of 22 GHz (Kadler et al. 2004b). Given NGC 1052’s linear scale of 0.093 pc/mas, the jet scale is ≈2 ⋅ 10¹⁶ cm per jet radius. This is compatible with the length scale in Table 2, as well as the estimated upper limit of the jet core size at 86 GHz (Baczko et al. 2016).

The initial 2D axisymmetric simulations were interpolated to a 3D cartesian grid via a Delaunay triangulation, assuming a viewing angle, ϑ. During the Delaunay triangulation, we stored the vertices and their weights. Once they are computed and kept in memory, additional triangulations represent an interpolation of data on the 3D Cartesian grid using the precomputed vertices and weights. This procedure significantly speeds up the radiative transfer calculations (for details on the numerical grid and conversion studies, see Fromm et al. 2018). The emission structure of the initial jet at 43 GHz is shown in Fig. 2. It should be noted that the radiative transfer leads to a projection of the source structure onto the x-y plane, as is seen in Fig. 1. This image was then rotated by ϕ = 15° in the plane of the sky to match the orientation of NGC 1052. At 43 GHz, the obscuring torus is almost entirely optically thin; therefore, the jet nozzle becomes visible in the emission map. However, due to the slight misalignment of the jet with respect to the observer (ϑ = 80°), the absorption along the line of sight is greater for the right jet, which manifests as a dimmer innermost jet structure. A more detailed view of the emission structure is presented in the bottom panel of Fig. 2, which displays the flux density along the jet axis (dashed white line in the top panel). The aforementioned recollimation shocks are clearly visible and can be found at a right ascension of RA = ±0.8, ±2.1, and ±4.5 mas. The last recollimation shock resembles a plateau rather than a local flux density maxima, as would be typical of a strong recollimation shock.

Fig. 2. Emission from the initial steady-state jet at 43 GHz. The top panel displays the jet structure and the bottom panel the evolution of the flux density along the jet axis (dashed white line).

We initialized the slow-light calculations by using the Delaunay triangulation of the last SRHD snapshot and searched for the jet location closest to the observer. From this position on, we integrated backward in time and loaded the corresponding SRHD snapshots including the plasma parameters. Since we kept the vertices and weights of the Delaunay triangulation in memory, the slow-light calculations could be performed in a fast and efficient way. In Fig. 3, we show the density distribution in the midplane (x-z plane for y = 0) seen by the photons for the slow-light (bottom) and for the fast-light (top) approaches. The rotation and broadening of the shocks are caused by the light travel time delay between the near side (front) and the far side (back) of the perturbation. These time delays are also responsible for the asymmetry between the left and right jet, which are, respectively, beamed toward and away from the observer (see also Aloy et al. 2003).

Fig. 3. Distribution of the density in a slice through the midplane (y = 0) of the jet. The jet is seen under a viewing angle of ϑ = 80° and the radiative transfer was performed in a positive z direction. The top panel shows the density without time delays and the middle one including time delays, where the left jet is beamed toward and the right jet away from the observer. The white boxes correspond to the zoom windows shown in the two bottom panels. Notice the rotation of the traveling shock in the slow-light case (bottom right panel).

In Fig. 4, we show the emission structure obtained from the simulations at 43 GHz for several snapshots separated by Δt = 10 R_j/c ≈ 115 days. The torus turns optically thin at ν ∼ 30 GHz, and thus we can clearly see the jet injection point (dotted red line) together with the first recollimation shocks (dotted green lines). Several traveling shock waves are visible and traceable (dotted gray lines). A deeper look into the structures of the traveling shocks reveals trailing shocks propagating in the wake of the main (forward) shock (see, e.g., Agudo et al. 2001; Fichet de Clairfontaine et al. 2022).

Fig. 4. Evolution of the jet structure at 43 GHz for different time frames, Δt ∼ 115 days. The red line indicates the location of the jet origin and two green lines mark the position of the first recollimation shock. The gray lines track the evolution of propagating shock waves.

Besides the temporal and spatial variations of the jet structure for a fixed frequency, our larger frequency coverage allows us to trace the spectral properties; for example, the turnover flux density, S_t, the turnover frequency, ν_t, and the optically thin spectral index, α_thin. Given the variation in the spectral index together with synchrotron self-absorption (SSA) in the jet and free-free absorption in the torus, we computed the optically thin spectral index from our last frequency bin; that is, between 900 GHz and 1 THz. At this frequency range, SSA and free-free absorption do not affect the spectrum. In Fig. 5, we present from top to bottom the turnover flux density, S_t, the turnover frequency, ν_t, and the optically thin spectral index, α_thin. The variation in the turnover position can clearly be seen from the plots. Due to the obscuring torus and high density and pressure at the injection region, the largest values for the turnover flux density and turnover frequency are located within ±1 mas of the origin.

Fig. 5. Spatial evolution of the spectral turnover position for a selected time frame. The panels show, from top to bottom, the turnover flux density, St, the turnover frequency, νt, and the optically thin spectral index, αthin.

In addition, the local maxima in the turnover position nicely trace the moving shocks. Since we do not include radiative cooling in our emission modeling, the increases in the turnover values are due to shock compression of the underlying flow. The median optically thin spectral index (bottom panel) is ⟨α_thin⟩= − 1.25 ± 0.1, which corresponds to the expected values for the used spectral slope of s = 3.5. During the radiative transfer, we numerically solved the integrals over the Bessel functions and did not use approximations. This leads to a deviation from the expected spectral index in the high-frequency emission, which explains the variation in the spectral index.

The unresolved total spectrum between 10⁸ Hz < ν < 10¹² Hz and its temporal variation over several time frames is shown in Fig. 6. The inlet provides a zoom into the turnover location showing the variation in the spectrum during the course of the simulation, which is connected to the propagation of perturbations and their interactions with the recollimation shocks.

Fig. 6. Temporal evolution of the broadband spectrum for our model. The different colors indicate selected time steps.

3. Synthetic data generation

3.1. Setup

To better compare the output of our simulations with observations, we aim to obtain a dataset that incorporates hurdles present in VLBI results. These include a sparse sampling of the u − v plane as a result of a limited number of telescopes, the presence of noise, and the effects of image reconstruction. These effects are not fully considered when blurring the emission maps from a simulation with an observing beam, as is typically done to imitate real observations. To put our dataset through an image reconstruction process, we used the synthetic observation tools of the imaging code eht-im (Chael et al. 2016, 2018).

We first decided on an observing array as well as an observation schedule to generate the projected baselines of a synthetic observation for 190 consecutive outputs. In this work, we have used the characteristics of the VLBA at five different observing frequencies, as well as certain GMVA stations for 86 GHz. The VLBA consists of ten identical 25 m antennas spread over the US, with a longest baseline of around 8600 km. Multifrequency observations can be carried out in different bands ranging from 0.3 GHz to 90 GHz, with eight of the ten telescopes being able to observe in the band centered at 86 GHz. The telescope sensitivities, given in system equivalent flux density (SEFD)³, of the VLBA corresponding to the observing frequencies used in this work are shown in Table 3. We also included GMVA⁴ stations at 86 GHz, as is given in Table 4. For this work, we used the VLBA bandwidth of Δν = 128 MHz.

With the information given in Tables 3 and 4, as well as as the coordinates of the stations, we computed the projected baselines for a synthetic observation with a duration of 15 hours (typical observing time for NGC 1052), a scan integration time of 12 s, and a cadence of 10 minutes between scans. As per Fromm et al. (2019), we furthermore considered that the telescopes could observe at elevations between 10° and 85°. The time elapsed between each consecutive output was Δt = R_j/4c ≈ 3 days. We started on October 9, 2003, for our first set of data and calculated the starting hour of the observation for each time step. For this purpose, we chose the VLBA station St. Croix as a reference antenna. From these parameters, we generated projected baselines, with Fig. 7 showing the resulting u − v coverage of one such synthetic observation at each of our working frequencies. We then generated synthetic visibilities by Fourier transforming and sampling the flux density distribution, I(x, y), obtained in Sect. 2.2 with the aforementioned baselines. This was done through the equation

$$\begin{array}{c}\hfill V_{\mathit{ij}}={\displaystyle \int \int I(x,y)}{e}^{-2\pi \text{i}(ux+vy)}\text{d}x\text{d}y.\end{array}$$(3)

Fig. 7. u − v coverage of a synthetic observation on October 9, 2003, at different frequencies. Note the different baselines for 86 GHz, due to the additional GMVA stations as well as the VLBA stations not operational at 86 GHz. Notice the different x and y scales.

During the calculation of V_ij, we added thermal noise related to the telescope sensitivities, as well as a time-dependent gain calibration error extracted from a Gaussian distribution with a standard deviation of 0.05 and a mean of 1 (for details, see Chael et al. 2018). Figure 8 shows the amplitudes of the final visibilities corresponding to the time step of Fig. 7, at the central frequency of 22 GHz.

Fig. 8. Total visibility amplitude plotted against u − v distance for 22 GHz. For clarity, every second data point is plotted with no error bars.

3.2. Image reconstruction

The last step to obtain our synthetic emission data was the reconstruction of the synthetic visibilities (for details on image reconstruction in radio astronomy, see Thompson et al. 2017), which we performed with the eht-im code adapted to our dataset. In particular, we incorporated the multifrequency functionality described in Chael et al. (2023). This method simultaneously reconstructs images at different frequencies by parametrizing the intensity as a function of the frequency with a Taylor expansion,

$$\begin{array}{c}\hfill log{\text{I}}_i=log\left[{\text{I}}_0\right]+\mathit{\alpha }log\left(\frac{{\nu }_i}{{\nu }_0}\right)+\mathit{\beta }{log}^2\left(\frac{{\nu }_i}{{\nu }_0}\right)+\dots ,\end{array}$$(4)

where α is the spectral index and β the spectral curvature. As these two parameters are a function of the coordinates, the reconstruction process returns a map of α and a map of β. This multifrequency approach to regularized maximum likelihood (RML) imaging allows the sharing of information across frequencies, acting as a regularizer by improving resolution at lower frequencies and aiding in the recovery of low-brightness structures at higher frequencies. As we can see in Fig. 6, the turnover frequency of our spectra is within our range of 8–86 GHz. As such, the spectral curvature cannot be neglected, and we include the second-order term on Eq. (4).

Before the reconstruction, we added station-dependent systematic noise to improve the convergence of the RML method. We used a Gaussian as a prior image, tilted by 15° as the object (see Chael et al. 2016). At each corresponding time step, this Gaussian had the same flux as the total intensity of the image at the central frequency 22 GHz.

The reconstruction process returned, for each time step, the maps of α and β used as parameters in Eq. (4). It also returned the output of the RML reconstruction for each of the five frequencies, and that same output convolved with its corresponding synthetic beam. The parameters of said beam associated with each frequency, averaged over every time step, are shown in Table 5.

Figure 9 shows the final result of the imaging process for all five frequencies on the initial time step of October 9, 2003. The blurred reconstructions (left) are contrasted with the direct output images blurred with the same beam (right). The frequency-dependent emission gap between the jets is clearly seen in both the reconstructed images and the blurred direct output. We note the asymmetry between the left, eastern jet and the right, western counterjet, as well as the locally enhanced emission due to the traveling shocks. The effects of the sparse u − v sampling as well as the limited dynamical range of the observing arrays are best seen at 43 GHz and 86 GHz.

Fig. 9. Reconstruction of a single time step at five frequencies (left column) contrasted with the direct emission maps (right column). Both have been blurred with the beam of their corresponding frequency, plotted in white.

4. Analysis

4.1. Jet kinematics

Identifying bright features, hereafter referred to as components, in jets and tracking their evolution over time allows us to study the jet dynamics including the propagation of shocks and their interactions. Furthermore, by comparing the identified and tracked components with the numerical simulations we can evaluate the impact of the image reconstruction on inferred observed jet dynamics.

As we generated data for 190 time steps (∼1.5 yr) to study over five frequencies, we set up an algorithm to find bright features over time as well as to link them together into trajectories. For this, we first used a “Laplacian of Gaussian”-based method implemented in scikit-image (van der Walt et al. 2014) (for details on feature finding methods, see Lindeberg 1994) on the reconstructed images. These images were previously convolved with a circular beam equal in area to the corresponding elliptical beam for its time step and frequency. On each separate frequency, we then made use of the Trackpy Python package (Allan et al. 2024) over every time step to identify and connect the center of the identified components into trajectories. We removed spurious trajectories from our dataset, identified as appearing for only a few frames. These spurious trajectories were mostly present on the 43 GHz and 86 GHz maps due to their noisier reconstructions and smaller convolving beams. We also ignored the components found in the last ≈1.7 mas on either side, corresponding to the beam size at 8 GHz.

The results for 43 GHz are plotted in Fig. 10, 43 GHz being the frequency that presents the largest number of detected trajectories, with dashed lines representing linear fits to them. The plots corresponding to the lower frequencies can be found in the Figs. A.1, A.2 and A.3. We contrasted the components across all frequencies and identified them according to their kinematics. The jet core, named C0, is visible in most time steps at this frequency but is obscured by the torus at the lower, optically thick frequencies. Moving components that are traveling outward are identified with the letter M, followed by an odd number for the left, eastern jet (top half of the figures) and an even number for the right, western jet (bottom half of the figures). We furthermore find at this frequency two fainter, harder to detect components moving closely behind M5 and M6. We identify these trajectories as trailing components, named T, followed by the number of the moving feature they are trailing with a slower speed.

Fig. 10. Distance of bright components from the central core, identified here as C0, over time for the 43 GHz maps. The components located in the left or right jet are marked with an odd or even number, respectively. Recollimation shocks are indicated by the letter R, moving components by the letter M, and their trailing shocks by the letter T, followed by the same number of the component which they trail. The dashed lines on each component correspond to a linear fit.

The two main recollimation shocks, identified as quasi-stationary components, R1 and R2, are located at roughly 1 mas on either side. While initially staying stationary as expected, they are notably perturbed halfway through our time range. Starting in 2005 in our simulations, two new moving components, M9 and M12, emerge from either recollimation shock. This disturbance is almost entirely gone in the 15 GHz and 22 GHz maps due to the larger beam size blurring adjacent features into each other. In the 8 GHz maps, these traveling perturbations are not detected until they are at least 5 mas away from the core and the recollimation shocks have instead merged with earlier traveling components M5 and M6, shifting both bright features farther out. The secondary recollimation shocks R3 and R4, being fainter than the main R1 and R2, were more difficult to find. However, they were located at roughly 4 mas on either side across most frequencies, particularly visible at 15 GHz and 22 GHz.

In Fig. 11, we show the reconstructions of the May 5, 2005, maps at the five different frequencies, with the central location and cross-identification of the components plotted on top. The primary recollimation shocks, R1 and R2, marked with two green lines, show what could be interpreted as an inward displacement between 8 GHz and 15 GHz. However, contrasting the position of the components in Figs. A.1, A.2, and A.3 shows that the combined effect of the obscuring torus and the larger resolving beams is the likely cause of this, obscuring inner components and blurring recollimation shocks and nearby traveling components into each other. Figure 11 also shows the newly emerging M9 and M12 components, as well as the more distant traveling shocks, M5 and M8, becoming too unresolved to be picked up by our algorithm as the frequency increases. As the figure displays one of the final time steps, the fainter secondary recollimation shocks, R3 and R4, are not visible, having been obscured, respectively, by the passing shocks, M5 and M6. While Fig. 2 shows three maxima in flux on either jet, the missing ones at roughly 2.5 mas may be obscured by the brightness of the moving components and of the nearby brighter shocks at 1 mas.

Fig. 11. Reconstructed and circularly blurred maps at all five frequencies in the May 5, 2005, snapshot. The components detected by our algorithm have been marked. The identified trajectories are connected through the dotted lines: white for the moving shocks, green for the recollimation shocks, and red for the jet core.

Table 6 shows the slope of the trajectories plotted in Fig. 10 with a dashed black line, as well as the trajectories for the other frequencies. As a component reaches the end of the simulated box, it cannot be tracked any further. This gives the impression that the components slow down as they reach the end of the simulated region. This is a limitation of the hydrodynamic simulation, and not a physical phenomenon. We marked with an asterisk the velocities that were computed up to an intermediate time step, determined from the residuals of a linear fitting.

As was expected, the jet core C0 only appears at the optically thin frequencies with a velocity compatible with zero. The recollimation shocks generally have values oscillating around zero, although perturbed by moving shocks and the limited resolution. For example, the faster velocity of R1 in 8 GHz is accompanied by an anomalously high velocity for M5 when compared to all four of the other frequencies. This disagreement is present in most trajectories visible at 8 GHz, most notably M7 and M8. M7, which shows the biggest disagreement, is difficult to pick up at this frequency until it travels enough to be far away from the faint recollimation shock R3. This is due to the combined effects of the large blurring beam and the obscuring torus, and it makes 8 GHz a nonoptimal frequency at which to study jet dynamics on this scale.

The traveling component M6 shows a largely consistent behavior across frequencies, as well as the late-appearing M12. M9 shows a larger disagreement between 43 GHz and 86 GHz data. The velocities computed for M5, M7, and M8 are also in agreement when 8 GHz data is excluded. We thus match them to the perturbations injected at the nozzle in the hydrodynamic simulation with a velocity of 0.5 c.

4.1.1. Positional errors

As was previously mentioned, we used an automated feature-finding algorithm to locate the components in our map. This was done due to the large number of images, 190 time steps over five frequencies. In observational studies, the component modeling and cross-correlation are typically done on the interferometric u − v plane (see Lister et al. 2009). As we performed our analysis on the image plane, we could not calculate positional errors based on the component fitting results as is described in Lister et al. (2009).

We assigned the location of each component an error, with the criteria being that it should depend on both the observing frequency and the flux density of the component. For the former, we related the error to half of the full width at half maximum (FWHM) of the circular blurring beam at each frequency and time step, b/2. For the latter, we assigned a weight to each component, w_i = S_i/S_max. S_i is the flux density at the center of a given component, while S_max is the maximum flux density found in the map at that frequency and time. As can be seen in Fig. 2, this tracks both the presence of shocks and the distance to the core, as the latter relates to the underlying flux profile.

By combining these two criteria, and following the relation between weights and variance, w_i = 1/σ_i², we have assigned to each of our components the positional error

$$\begin{array}{c}\hfill \delta x_i=\sqrt{\frac{S_{\text{max}}}{S_i}}\frac{b}{2}·\end{array}$$(5)

This equation produces a minimum error of b/2, larger than the minimum error of b/10 adopted by Baczko et al. (2019). This more conservative approach is warranted to ensure that adjacent components identified by the automated feature-finding algorithm in the image plane are truly distinct. It can be seen in Fig. 10, and especially in Fig. 12, that for the higher frequencies this error displays more time variability. This is caused by the low flux density at these frequencies affecting the quality of the reconstructions.

Fig. 12. Distance of bright components from the central core over time for the 86 GHz maps. The components are marked by a letter and a number, as is explained in Fig. 10. The dashed lines on each component correspond to a linear fit.

4.1.2. Shock-shock interaction

The maps at 86 GHz contain the lowest flux, which hinders the reconstruction process. This strong impact of the reconstruction process on the data is seen in Fig. A.4. While this means that outer components cannot be properly detected, the low beam size at this frequency does allow us a closer look into the innermost areas of the jet. In Fig. 12, we have plotted the distance of the components from the core for 86 GHz, where it can be observed that features become too unresolved for detection at roughly 4 mas. The trajectories that can be found at this frequency are consistent with those seen in Fig. 10 at 43 GHz.

The resolution at these two frequencies allows us to see the two newly emerging components, M9 and M12, as well as an oscillatory movement of the recollimation shocks. Using the axial emission of the initial steady-state model at 43 GHz (see Fig. 2) as a guide, we would expect three recollimation shocks on either jet. R3 and R4 can be identified with the peaks between 4 mas and 5 mas, but they disappear early on in our simulation due to the passing of traveling shocks. R1 and R2 appear at all frequencies and most time steps, but seem to travel away from and back to their initial position when crossed by new components. Figure 13 shows five different time steps of the reconstructed 86 GHz maps with highlighted components. The last two time steps show M9 and M12 emerging as well as, as was mentioned, R1 and R2 moving inward. We shall discuss later the possible reasons for the behavior of the standing shocks, and the impact of moving shocks on jet structure.

Fig. 13. Five different snapshots of reconstructed and circularly blurred 86 GHz maps. The components detected by our algorithm have been marked. The identified trajectories are connected through the dotted lines: white for the moving shocks, green for the recollimation shocks, and red for the jet core.

4.1.3. Viewing angle and intrinsic speed estimation

The apparent speeds of traveling components allow us to probe the symmetry in an independent way. Assuming that the jets are intrinsically symmetric, the ratio, R, of the total flux density between the approaching (left) and receding (right) jets is described by the equation:

$$\begin{array}{c}\hfill R=\frac{S_{\text{a}}}{S_{\text{r}}}={\left(\frac{1+\beta cos\theta }{1-\beta cos\theta }\right)}^{2-\alpha }.\end{array}$$(6)

β is the intrinsic jet speed in units of c, θ is the viewing angle, and α was set to −1 following Baczko et al. (2019). Furthermore, the apparent speeds of Table 6 relate to the intrinsic jet speed and viewing angle as

$$\begin{array}{c}\hfill {\beta }_{\text{app,a/r}}=\frac{\beta sin\theta }{1\mp \beta cos\theta }·\end{array}$$(7)

The overlap of Eq. (7) with Eq. (6) defines the allowed parameter space for β and θ. If symmetry holds, the two parameter spaces created by both signs of Eq. (7) should furthermore overlap with each other, within the derived uncertainties. Figure 14 shows the result of the analysis for our reconstructed data at the central frequencies of 22 GHz and 43 GHz. We used the time-averaged mean of R and its standard deviation, while the values for β_app were taken from Table 6. We chose the components M5 and M6, as they were the main traveling shocks during our time range and were picked up with minimal scatter for most frequencies. These properties indicate they are likely traveling shocks, and thus faster than the underlying flow, and so the following results should be taken as an approximate limit.

Fig. 14. Allowed parameter space for the intrinsic speed, β, and viewing angle, θ, for both jets, at 22 GHz (left) and 43 GHz (right). The orange and blue lines represent Eq. (7); the green lines Eq. (6). The dotted and dashed lines represent the lower and upper error ranges, respectively. Both cases allow intrinsic symmetry, as is shown by the overlapping shaded areas.

It can be seen that both cases allow for intrinsic jet symmetry, even with the more restrictive errors of 43 GHz. Furthermore, both parameter spaces contain the ground truth of β = 0.5 that was set into our simulations. The viewing angle of 80° is overestimated by roughly 5°. As the allowed viewing angles are smaller for 22 GHz, an optically thicker frequency, it is possible that our flux ratio, R, was overestimated due to the core emission. Still, this result shows that an intrinsically symmetric jet could be retrieved from the data even after undergoing a reconstruction process.

4.2. Jet collimation profiles

In this section, we analyze the profiles of the width – that is, the jet collimation profiles of both the jet and the counterjet – using the following method. As in the previous section, the flux maps resulting from the imaging process were blurred with a circular beam equal in area to the corresponding elliptical beam for its time step and frequency. The image was then rotated back by 15° so that the jet spine was aligned with the x axis, for ease of computing. From there, the peak flux and jet width were evaluated at every slice perpendicular to the jet spine. To compute the jet width, a Gaussian fit was performed on each slice, from which we obtained a FWHM equivalent to the width of the blurred jet. A flux-dependent error was incorporated into the Gaussian fit. This error was computed from taking the mean value of an empty region of the flux map to obtain an average value of the absolute error, and dividing it by the map to obtain a map of relative errors.

The final jet width was calculated by deconvolving the circular blurring beam and the width of the Gaussian fit, obtaining a jet width of $\sqrt{{\text{FWHM}}_{\text{fit}}^2-{\text{FWHM}}_{\text{beam}}^2}$. This was done over each time step of every frequency, and these deconvolved values were the ones we used for the purpose of studying jet collimation.

Following the method used for the study of the collimation of the NGC 1052 jets in Baczko et al. (2022), a broken power-law model was fit to the data at all frequencies simultaneously,

$$\begin{array}{c}\hfill w(d)=W_0{2}^{(k_{}\text{u}-k_{}\text{d})/s}{\left(\frac{d}{d_{}\text{b}}\right)}^{k_{}\text{u}}{[1+{\left(\frac{d}{d_{}\text{b}}\right)}^s]}^{(k_{}\text{d}-k_{}\text{u})/s},\end{array}$$(8)

where w is the jet width as a function of the distance to the core d and d_b is the breakpoint between the upstream and downstream sections of the jet, defined, respectively, as the region between the core and the breakpoint, and the region beyond the breakpoint. W₀ is the width at the breakpoint, and k_u/d are the power-law slopes of the upstream and downstream sections, respectively. s is a free parameter that controls the smoothness of the transition between power laws, set to 10 following Nakahara et al. (2020) and Baczko et al. (2022). The position of the breakpoint is expected to correspond on each side to the primary recollimation shock, an area characterized by its narrow width and high flux. A recollimation shock is present at approximately 1 mas from the core at all time steps, pointed out by green lines in Fig. 4. These can be seen in the reconstructed images in Fig. 9, particularly at 43 GHz where there is no significant obscuring from the torus. It would be expected for d_b in Eq. (8) to match these main recollimation shocks, which mark the transition from the upstream jet to the downstream one. While our algorithm would have an easier time positioning this d_b on our direct emission maps, it is of interest for this paper to examine the ways in which this analysis was affected by the noise and limited resolution of a reconstruction process. Figure 15 shows the result of the broken power-law fitting for one time step. This fitting was done simultaneously for all frequencies, performed separately over the left jet and right counterjet after removing the invalid values of low-flux regions. We focus our analysis on the values of the downstream slope, k_d, for which we obtained single epoch relative errors on the order of 0.1% for most time steps.

Fig. 15. Width and peak flux profiles for the left (eastern) and right (western) jet at a single time step. The solid black line represents the broken power law resulting from fitting all frequencies together.

We performed a kernel density estimation (KDE) analysis on the k_d distribution of both jets for all time steps. This allowed us to estimate a probability density function and to draw the values of the median as well as the quartiles of the distribution. The results are plotted in Fig. 16, where the solid black line drawn over the histogram is the KDE result. We show the median as well as the first and third quartiles, which we used as a range. The resulting $k_{\mathrm{d},\mathrm{l}}=0.{72}_{-0.37}^{+0.20}$ for the left jet and $k_{\mathrm{d},\mathrm{r}}=0.{60}_{-0.13}^{+0.12}$ for the right counterjet point toward the jet and counterjet not being significantly asymmetric on a global scale. As a test, we ran the analysis again, removing k_d, l > 1.5 values from our dataset. These abnormally high values of the downstream slope could correspond to a failure of the fitting algorithm to identify the position of the recollimation shock as the breakpoint of the broken power law. The new analysis constrained the median value further, $k_{\mathrm{d},\mathrm{l}}=0.{67}_{-0.13}^{+0.11}$. This is seen as the orange line in the left panel of Fig. 16.

Fig. 16. Kernel density estimation of the distribution of downstream slopes resulting from the broken power-law fitting.

A more detailed study of the validity of these fittings is found in Appendix B. Carrying out this analysis on the direct emission maps supports our conclusion that the tail present in Fig. 16 is an artifact introduced by failed fittings on reconstructed images. Furthermore, the range provided by our KDE calculation narrows by a factor of ∼2, while the new values of $k_{\mathrm{d},\mathrm{l}}=0.{69}_{-0.11}^{+0.05}$ and $k_{\mathrm{d},\mathrm{r}}=0.{66}_{-0.07}^{+0.14}$ remain statistically equivalent to those of Fig. 16. While it is of interest for this study to focus on the impact of the reconstruction process on the study of jet kinematics, this test is in agreement with the intrinsically symmetric nature of our simulated observation.

Still, downstream profiles with a significant asymmetry were found when looking at specific time steps. Figure 17 shows a snapshot where the downstream jet has a wider collimation profile than the counterjet and vice versa, as well as a case where they are roughly symmetric. This shows that while an asymmetry may arise at specific time steps due to light aberration and the shock structure of the jet caused by a traveling perturbation, a significant asymmetry cannot be claimed on the global profiles over all epochs.

Fig. 17. Emission maps at 43 GHz with their resulting broken power-law fittings in solid white lines, as well as the values for the downstream slopes for eastern and western jets.

5. Discussion

In this paper, we present new results from numerical simulations in which we have incorporated the aberration caused by the finite speed of light. Past studies showed that relativistic beaming may not be enough to properly explain the dynamics of apparent superluminal motion in AGN jets, making light aberration absolutely necessary. Not accounting for these slow-light calculations may lead to an underestimation of the flux of the moving components, and may also affect the overall profile of the jet. Concerning the flux, following the analysis of Komissarov & Falle (1997) with our values for the viewing angle and injection velocity, 80° and 0.5 c, we calculate that without slow-light we would have underestimated the flux of the moving components by around 10%. This is less dramatic than the example given of a typical source with apparent superluminal motion, where the flux could be underestimated by an order of magnitude. This is expected for two-sided sources observed at subluminal speeds such as NGC 1052, whose large viewing angle makes it ideal for observational comparisons. Still, the introduction of light aberrations affects not only the observed brightness of the moving components, but also their size and shape. Figure 3 shows how accounting for light travel time not only introduces a bipolar asymmetry, but also alters the collimation profile of the jets due to the rotation and broadening of the shocks.

We also introduced a synthetic imaging method whereby several VLBI frequencies are reconstructed simultaneously, allowing for information to be shared between the different frequencies. Chael et al. (2023) introduce the method and show the results of computing a map of the spectral index during the reconstruction. They also show the necessity of using nonzero values of the spectral curvature in frequency ranges where the curvature is not negligible, as is the case here. As we covered frequencies from 8 GHz to 86 GHz, this significantly improved our reconstruction.

During our component tracking, we found two trailing shocks, T5 and T6. These were only found in the 43 GHz maps, where the noise in the reconstruction is much larger, to the point where the detection threshold had to be increased, and thus they could be spurious features. To better understand to which extent these features are affected by the reconstruction process, Fig. A.4 shows our trajectories on top of a map where the one-dimensional flux of the jet spine is plotted against time. This is done for the reconstructed images as well as for the blurred direct maps. It does appear visually as if trailing shocks could be accompanying the moving components. This is particularly visible for components M7, M5, and M6. However, it appears as though our tracking algorithm in its current settings has difficulties in separating them from the main components. The trailing component T5 could, however, be retrieved near the beginning and end of the direct maps.

Agudo et al. (2001) and Mimica et al. (2009) carried out SRHD simulations that showed that perturbations traveling along a jet can lead to the formation of trailing shocks emerging behind a sufficiently strong moving component. They furthermore concluded that these trailing components should move more slowly than the main feature that they are trailing. In addition, Fichet de Clairfontaine et al. (2022) identify these trailing components with relaxation shocks that arise from the interaction of a moving shock passing through a standing shock. The slower velocities we obtained for T5 and T6 are consistent with these findings. These works furthermore study the oscillation of standing components after being perturbed by a passing shock. This effect is visible in our results, particularly for the higher-resolution 43 GHz and 86 GHz maps. Rather than staying at the same distance on either side, the recollimation shocks appear to move inward halfway through our simulation. For our selected time range, the recollimation shocks are disturbed by the passing of the shocks M5 and M6, before the starting point of the dataset shown in this work. Starting roughly halfway, these recollimation shocks slide inwards, and a possible outward motion is starting to appear near the end, caused by the newly emerging shocks M9 and M12. This is displayed in Fig. A.4.

Several recent works study the symmetry of NGC 1052, making it a good source to compare our results to. These studies fit a broken power law to the width profile of the jet, yet there is no consensus on whether an asymmetry is present. Using 2000–2001 observations, Nakahara et al. (2020) calculate the upstream values of the slope to be consistent with zero for both jet and counterjet, and a downstream slope value of 0.96 ± 0.04 for the eastern jet and 1.02 ± 0.03 for the western jet, reaching the conclusion of symmetry. On the other hand, using 2017 observations, Baczko et al. (2022) calculate 1.01 ± 0.01 downstream for the eastern jet and jet 1.22 ± 0.02 for the western jet, with an upstream value of 0.17 for both. They carry out the same analysis on stacked VLBA observations from 2005 to 2009, calculating downstream values of 0.80 ± 0.01 for the eastern jet and 1.22 ± 0.05 for the western jet, and upstream values of around 0.20 for both. Kovalev et al. (2020) only calculate the slope values for the eastern jet, obtaining a steeper value of 0.391 ± 0.048 for the upstream jet, and 1.052 ± 0.081 downstream.

By comparison, carrying out the same broken power-law fitting on our reconstructed images yields upstream slopes of k_u ≤ 0.5 for the majority of time steps. 93% and 85% of the time steps show k_u ≤ 0.25 for the left and right jet, respectively, while this percentage rises to 100% when repeating this analysis on the infinite-resolution images. This collimation profile between cylindrical and parabolic (0 ≤ k ≤ 0.5) is consistent with recent studies. Furthermore, the final analysis of our reconstructed images yields downstream slopes of $k_{\mathrm{d},\mathrm{l}}=0.{67}_{-0.13}^{+0.11}$ for the left jet and $k_{\mathrm{d},\mathrm{r}}=0.{60}_{-0.13}^{+0.12}$ for the right counterjet. Testing our analysis on the direct emission maps narrows the range further to $k_{\mathrm{d},\mathrm{l}}=0.{69}_{-0.11}^{+0.05}$ and $k_{\mathrm{d},\mathrm{r}}=0.{66}_{-0.07}^{+0.04}$. While this is significantly less conical than VLBI observations of NGC 1052, the underlying collimation profile is determined by the chosen pressure gradient of the ambient medium. The aim of this work is to determine whether it is possible to induce asymmetries to the extent found in Baczko et al. (2022), when starting from an intrinsically symmetric double-sided jet. Our results point to the conclusion that perturbing the jet and introducing light travel time delays as well as an obscuring torus is, on its own, not enough to induce a global bipolar asymmetry. We carried out a complementary study of the allowed parameter space and found consistent parameters for both jets, supporting our conclusion of symmetry. However, while this is the result of doing the analysis over an extended period of time, asymmetries may still arise at individual time steps. Studies focusing on inducing asymmetries in the jet launching stage could give us further insight into the morphology and evolution of jets, and the astrophysical processes behind them. Important new constraints will come from future observational facilities at high radio frequencies such as the next-generation Event Horizon Telescope (ngEHT; Doeleman et al. 2023) and the next-generation Very Large Array (ngVLA; Murphy et al. 2018; Kadler et al. 2023).

In the presented work, we have used 2D hydrodynamic simulations. Besides the suppression of 3D instabilities due to the assumed axisymmetry of the jets, there are other limitations to our model, such as the lack of magnetic fields in the hydrodynamic simulations. The effects of magnetic fields can be split into two categories: the impact on the dynamics and the impact on the computed emission. If we assume that the jets are not magnetically dominated on parsec scales, the sound speed will be larger than the Alfvén speed. In this case, the oscillation of the jet surface will be determined by the sound speed and we expect no significant alteration of the jet shape compared to pure SRHD simulations. On the other hand, if the jets are magnetically dominated the Alfvén speed could be larger than the sound speed. In this case, the wavelength of the oscillations will decrease and a more cylindrical collimation profile (k = 0) could be obtained (see, e.g., Porth & Komissarov 2015). Given the cylindrical collimation profile found in some observations, this scenario could be analyzed in future studies but is out of the scope of this work. On the emission side, magnetic fields can introduce a top-bottom asymmetry in relativistic jets due to the aberration of light and pitch angle of the magnetic field (see, e.g., Fuentes et al. 2018). While this aberration is typically small for large inclinations and small velocities, and has not been studied in this work, the presence of magnetic fields could introduce a general asymmetry in a double-sided jet system.

Acknowledgments

It is a pleasure to thank L. Rezzolla for useful comments and suggestions. This research is supported by the DFG research grant “Jet physics on horizon scales and beyond” (Grant No. 443220636) within the DFG research unit “Relativistic Jets in Active Galaxies” (FOR 5195). Support also comes from the European Research Council Advanced Grant “JETSET: Launching, propagation and emission of relativistic jets from binary mergers and across mass scales” (grant no. 884631). The numerical simulations and calculations have been performed on MISTRAL at the Chair of Astronomy at the JMU Wuerzburg and on Iboga and Calea at the Institute for Theoretical Physics (ITP) at University of Frankfurt. ASP and MP acknowledge support by the Spanish Ministry of Science through Grant and PID2022-136828NB-C43. MP acknowledges support from the Generalitat Valenciana through grant CIPROM/2022/49, and from the Astrophysics and High Energy Physics programme supported by the Ministry of Science and Innovation and Generalitat Valenciana with funding from European Union NextGenerationEU (PRTR-C17.I1) through grant ASFAE/2022/005. YM is supported by the National Key R&D Program of China (grant no. 2023YFE0101200), the National Natural Science Foundation of China (grant no. 12273022), and the Shanghai municipality orientation program of basic research for international scientists (grant no. 22JC1410600).

References

Appendix A: Tracking results

Figures A.1, A.2 and A.3 show the component trajectories for 8 GHz, 15 GHz and 22 GHz, respectively. These have been computed with the same method described in Section 4.1.

Fig. A.1. Distance of bright components from the central core over time for the 8 GHz maps. The components are marked by a letter and a number as explained in Fig. 10. The dashed lines on each component correspond to a linear fit.

Fig. A.2. Distance of bright components from the central core over time for the 15 GHz maps. The dashed lines on each component correspond to a linear fit.

Fig. A.3. Distance of bright components from the central core over time for the 22 GHz maps. The dashed lines on each component correspond to a linear fit.

Fig. A.4 compares the results of applying our tracking algorithm on the reconstructed images (bottom) and the direct, blurred ones (top) for 43 GHz and 86 GHz. The minimum detection threshold was increased for the reconstructed images due to the noise present at these frequencies.

Fig. A.4. Evolution over time of the flux of the jet spine. Each vertical slice shows the flux of the jet spine for a given time step, following the positional angle of 15°, showing the evolution of the jets’ shock structure over time. The direct emission maps (top) have been convolved with the same beam as the corresponding reconstructed maps (bottom). Overplotted in cyan are the components detected by separately applying our tracking algorithm to each of the maps. Note the effect of adjacent shocks in the tracking algorithm, such as the oscillation of the recollimation shocks when a moving shock passes through, or the jet core temporarily going undetected in the direct maps when the recollimation shock at ∼1 mas is dragged outward and expands. Note also the faint trailing components visible in the direct images. The direct emission maps also show the purely numerical effect of the moving shocks “slowing down" as they reach the end of the simulated box.

Appendix B: Fitting results

Here we provide a more careful study of the results of fitting Eq. 8 to our width profiles.

Fig. B.1 shows the resulting distribution of all four fitting parameters, for the left and right jets. The downstream slopes, k_d, are the same ones found in Fig. 16. Before the breakpoint, the fitting results in a less than parabolic profile (k_u < 0.5) in > 90% of the time steps, for either jet or counterjet, while k_u < 0.25 for > 80% of time steps. After the breakpoint, a downstream value of k_d > 1 was found in 5% of the time steps for the right jet and in 15% for the left jet. This results in a tail following the central distribution.

Fig. B.1. Distribution of parameters resulting from fitting a broken power-law to our reconstructed maps. Note the tail present on the distribution of the breakpoints, db.

Fig. B.2. Emission map at 43 GHz showing the overshot fitting of the downstream eastern jet. The solid white line represents the broken power-law fit of this specific time step.

These greater than conical values could have resulted from the breakpoint being placed further out than the primary recollimation shocks during the fitting process. We checked this by plotting the result of our fitting on top of its corresponding emission map at 43 GHz. Figure B.2 shows an example of such a reconstruction. On the left jet, the breakpoint was placed at a distance of d_b, l = 3.94 mas from the core. It appears likely that missing information on particularly faint outer features leads to the breakpoint being identified with a cylindrical region at this distance. When extended farther out, we clearly observe the failed fitting overestimates the width at the final few mas of the jet, caused by a greater than conical expansion (k_d > 1). For these reasons, we decided to re-run the analysis in Section 4.2 removing all values k_d > 1.5.

To further validate our analysis, we ran our fitting algorithm on the direct emission maps. The fitting parameters are shown in Fig. B.3. Note the difference in x-axis values when compared to Fig. B.1. The resulting fitting values were constrained to a smaller range for every parameter. The breakpoints do not show a tail in their distribution, with their maximum distance found to be d_b, l = 1.2 mas. This is accompanied by every value of the downstream slope being less than conical, k_d < 1. In Fig. B.4, we ran the same analysis shown in Fig. 16 with this data. The resulting median values support our conclusion that there is no global asymmetry present, in agreement with our intrinsically symmetric simulation. This analysis also shows the impact which a finite-resolution array has on the study of jet kinematics.

Fig. B.3. Distribution of parameters resulting from fitting a broken power-law to our direct, infinite resolution maps. Note the narrower range of values when compared to Fig. B.1.

Fig. B.4. Kernel density estimation of the distribution of downstream slopes resulting from the broken power-law fitting.

Finally, Fig. B.5 shows some of the parameters of the infinite-resolution fitting of Fig. B.3 as a function of time. We focus on the infinite-resolution maps as the reconstruction process introduces significant scatter and makes it challenging to find such clear trends, though the overall distributions point to similar results as seen in Figs. 16 and B.4. The results of Fig. B.5 show the dependence between the collimation profile and the breakpoint location, where an inwards displacement of d_b seems to be accompanied by a narrower collimation profile. Furthermore, when contrasting this to Fig. A.4, the dates where both values of d_b drop roughly correspond to the dates where the recollimation shocks at ±1 mas display inward motion. As previously discussed, this oscillation is caused by the recollimation shock being dragged downstream by a moving shock. This emphazizes the importance of traveling components on the overall jet profile, as they may introduce asymmetries on specific epochs when combined with observational effects.

Fig. B.5. Parameters resulting from fitting a broken power-law to our infinite resolution maps, as a function of time. The downstream slopes for both jets (top panel) and the ratio between them (bottom panel) are contrasted to the breakpoint location (central panel). A dashed line in the bottom panel indicates a ratio of 1.

All Tables

All Figures

	Fig. 1. Distribution of the density of the initial steady-state jet model in a slice through the midplane (y = 0) of the jet. The jet is seen under a viewing angle of ϑ = 80°. The top panel shows the 2D density distribution and the bottom panel a density evolution along the jet axis (dashed white line).
In the text

	Fig. 2. Emission from the initial steady-state jet at 43 GHz. The top panel displays the jet structure and the bottom panel the evolution of the flux density along the jet axis (dashed white line).
In the text

Fig. 3. Distribution of the density in a slice through the midplane (y = 0) of the jet. The jet is seen under a viewing angle of ϑ = 80° and the radiative transfer was performed in a positive z direction. The top panel shows the density without time delays and the middle one including time delays, where the left jet is beamed toward and the right jet away from the observer. The white boxes correspond to the zoom windows shown in the two bottom panels. Notice the rotation of the traveling shock in the slow-light case (bottom right panel).

In the text

	Fig. 4. Evolution of the jet structure at 43 GHz for different time frames, Δt ∼ 115 days. The red line indicates the location of the jet origin and two green lines mark the position of the first recollimation shock. The gray lines track the evolution of propagating shock waves.
In the text

	Fig. 5. Spatial evolution of the spectral turnover position for a selected time frame. The panels show, from top to bottom, the turnover flux density, St, the turnover frequency, νt, and the optically thin spectral index, αthin.
In the text

	Fig. 6. Temporal evolution of the broadband spectrum for our model. The different colors indicate selected time steps.
In the text

	Fig. 7. u − v coverage of a synthetic observation on October 9, 2003, at different frequencies. Note the different baselines for 86 GHz, due to the additional GMVA stations as well as the VLBA stations not operational at 86 GHz. Notice the different x and y scales.
In the text

	Fig. 8. Total visibility amplitude plotted against u − v distance for 22 GHz. For clarity, every second data point is plotted with no error bars.
In the text

	Fig. 9. Reconstruction of a single time step at five frequencies (left column) contrasted with the direct emission maps (right column). Both have been blurred with the beam of their corresponding frequency, plotted in white.
In the text

Fig. 10. Distance of bright components from the central core, identified here as C0, over time for the 43 GHz maps. The components located in the left or right jet are marked with an odd or even number, respectively. Recollimation shocks are indicated by the letter R, moving components by the letter M, and their trailing shocks by the letter T, followed by the same number of the component which they trail. The dashed lines on each component correspond to a linear fit.

In the text

	Fig. 11. Reconstructed and circularly blurred maps at all five frequencies in the May 5, 2005, snapshot. The components detected by our algorithm have been marked. The identified trajectories are connected through the dotted lines: white for the moving shocks, green for the recollimation shocks, and red for the jet core.
In the text

	Fig. 12. Distance of bright components from the central core over time for the 86 GHz maps. The components are marked by a letter and a number, as is explained in Fig. 10. The dashed lines on each component correspond to a linear fit.
In the text

	Fig. 13. Five different snapshots of reconstructed and circularly blurred 86 GHz maps. The components detected by our algorithm have been marked. The identified trajectories are connected through the dotted lines: white for the moving shocks, green for the recollimation shocks, and red for the jet core.
In the text

	Fig. 14. Allowed parameter space for the intrinsic speed, β, and viewing angle, θ, for both jets, at 22 GHz (left) and 43 GHz (right). The orange and blue lines represent Eq. (7); the green lines Eq. (6). The dotted and dashed lines represent the lower and upper error ranges, respectively. Both cases allow intrinsic symmetry, as is shown by the overlapping shaded areas.
In the text

	Fig. 15. Width and peak flux profiles for the left (eastern) and right (western) jet at a single time step. The solid black line represents the broken power law resulting from fitting all frequencies together.
In the text

	Fig. 16. Kernel density estimation of the distribution of downstream slopes resulting from the broken power-law fitting.
In the text

	Fig. 17. Emission maps at 43 GHz with their resulting broken power-law fittings in solid white lines, as well as the values for the downstream slopes for eastern and western jets.
In the text

	Fig. A.1. Distance of bright components from the central core over time for the 8 GHz maps. The components are marked by a letter and a number as explained in Fig. 10. The dashed lines on each component correspond to a linear fit.
In the text

	Fig. A.2. Distance of bright components from the central core over time for the 15 GHz maps. The dashed lines on each component correspond to a linear fit.
In the text

	Fig. A.3. Distance of bright components from the central core over time for the 22 GHz maps. The dashed lines on each component correspond to a linear fit.
In the text

Fig. A.4. Evolution over time of the flux of the jet spine. Each vertical slice shows the flux of the jet spine for a given time step, following the positional angle of 15°, showing the evolution of the jets’ shock structure over time. The direct emission maps (top) have been convolved with the same beam as the corresponding reconstructed maps (bottom). Overplotted in cyan are the components detected by separately applying our tracking algorithm to each of the maps. Note the effect of adjacent shocks in the tracking algorithm, such as the oscillation of the recollimation shocks when a moving shock passes through, or the jet core temporarily going undetected in the direct maps when the recollimation shock at ∼1 mas is dragged outward and expands. Note also the faint trailing components visible in the direct images. The direct emission maps also show the purely numerical effect of the moving shocks “slowing down" as they reach the end of the simulated box.

In the text

	Fig. B.1. Distribution of parameters resulting from fitting a broken power-law to our reconstructed maps. Note the tail present on the distribution of the breakpoints, db.
In the text

	Fig. B.2. Emission map at 43 GHz showing the overshot fitting of the downstream eastern jet. The solid white line represents the broken power-law fit of this specific time step.
In the text

	Fig. B.3. Distribution of parameters resulting from fitting a broken power-law to our direct, infinite resolution maps. Note the narrower range of values when compared to Fig. B.1.
In the text

	Fig. B.4. Kernel density estimation of the distribution of downstream slopes resulting from the broken power-law fitting.
In the text

	Fig. B.5. Parameters resulting from fitting a broken power-law to our infinite resolution maps, as a function of time. The downstream slopes for both jets (top panel) and the ratio between them (bottom panel) are contrasted to the breakpoint location (central panel). A dashed line in the bottom panel indicates a ratio of 1.
In the text