© The Authors 2025

1. Introduction

Active galactic nuclei (AGNs) are known to produce powerful relativistic jets that can propagate on scales of up to thousands of kiloparsec (see, for example Blandford et al. 2019, and references therein). One key parameter that regulates the jet launching and propagation is the magnetic field that permeates the entire system, from the central region down to the extended jets. The relativistic jets are expected to be launched in the form of Poynting-flux dominated outflows, and their initial magnetic field strengths lie in the range of hundreds to thousands of Gauss. They decrease to milliGauss on parsec scales (see, for example, O’Sullivan & Gabuzda 2009). On (sub-)parsec scales, magnetic fields play a crucial role in collimating and accelerating the outflows. The magnetic energy is converted into kinetic energy of the bulk flow until equipartition is reached (Vlahakis & Königl 2003a,b, 2004; Komissarov et al. 2007; Komissarov 2012; Ricci et al. 2024). On compact scales, multiple studies highlighted helical magnetic field configurations (see, for example, Hovatta et al. 2014; Gabuzda 2018), which indicates that the toroidal and poloidal field components both play an important role on these scales. Whether the magnetic field structure is universal in AGN jets or varies across different types of sources is an ongoing investigation, however (see, for example, Boccardi et al. 2017; Hodge et al. 2018, and references therein).

When the event horizon scales are considered, their magnetization and the hot and radiatively inefficient accretion disks can be divided into the two main models of standard and normal evolution disks (Narayan et al. 2012) and magnetically arrested disks (MAD, Bisnovatyi-Kogan & Ruzmaikin 1974, 1976; Narayan et al. 2003; Tchekhovskoy et al. 2011). The first model represents the low-magnetization scenario, in which particles lose angular momentum through magnetorotational instability. The MADs, in contrast, are highly magnetized disks that are permeated by strong poloidal magnetic fields that can saturate and disrupt the accretion flow. They form a magnetosphere that surrounds the central engine. The MAD model has gained increasing visibility in recent years because it can successfully model the launching of powerful jets from hot accretion disks (see, for example, Narayan et al. 2022). It is currently proposed for an increasing number of low-luminosity AGNs (Zamaninasab et al. 2014; Event Horizon Telescope Collaboration 2021; Yuan et al. 2022; Ricci et al. 2022).

In this context, the mechanisms described above are expected to produce a pair of largely similar, parallel, and opposite outflows, so-called jets. It remains an open question, however, whether the jets that are launched from AGNs are actually symmetric. Alongside the differences in the surrounding medium between the two regions that are crossed by the outflows, asymmetric jet production was theoretically shown to be possible as a consequence of persistent asymmetries in the accretion disk (see, for example, Wang et al. 1992; Fendt & Sheikhnezami 2013; Nathanail et al. 2020). Recently, new works have shown that even in symmetrically launched jets, asymmetries in the observed physical parameters can arise from the passage of jet instabilities (e.g., shocks) downstream of both jets as a consequence of the different viewing angles between the two sides (Fromm et al. 2018; Saiz-Pérez et al. 2025). From an observational point of view, it is unknown whether the triggered asymmetries vary significantly in time.

The magnetic field content of the relativistic jets and their central regions, along with possible asymmetries in the jet production mechanisms, can be investigated through very long-baseline interferometry (VLBI) observations of jets on subparsec to parsec scales. For example, VLBI images may highlight differences in the velocity and collimation profiles, and in the overall morphology of the two jets (Baczko et al. 2019; Saiz-Pérez et al. 2025). The best sources to explore the rise of jet asymmetries are double-jetted radio galaxies, that is, AGNs seen at large viewing angles. The reduced projection effects enable us to detect both jets.

We focus on NGC 1052, which is a frequently studied double-sided radio galaxy for which an asymmetric jet production was previously suggested (Baczko et al. 2019). NGC 1052 is a nearby source, with a redshift of z = 0.005 (Koss et al. 2022) and the mass of the central black hole (BH) is estimated to be M_BH = 10^8.2 M_⊙ (Woo & Urry 2002). The source was frequently studied by means of VLBI observations, which allowed us to gain a deep insight into the properties of its jets, such as their collimation¹, propagation speed, and magnetic field content (Kameno et al. 2001; Vermeulen et al. 2003; Kadler et al. 2004b; Baczko et al. 2016, 2019, 2022; Baczko et al. 2024). The accretion disk in NGC 1052 was proposed to be an advection-dominated accretion flow (ADAF, see, for example, Falocco et al. 2020). At larger scales, the nuclear region of NGC 1052 is embedded in a geometrically thick, clumpy molecular torus with a radius of 2.4 ± 1.3 pc (Kameno et al. 2020). The column density of the neutral hydrogen (HI) was constrained to be between 10²² cm⁻²–10²⁴ cm⁻² (Weaver et al. 1999; Guainazzi & Antonelli 1999; Kadler et al. 2004a), while more recent measurements constrained higher values for the molecular hydrogen (H2) of (3.3 ± 0.7)×10²⁵ cm⁻² (Kameno et al. 2020). The torus presents an ionized component that extends 0.1 pc (∼1 mas) toward the eastern jet and 0.7 pc (∼7.5 mas) toward the western jet (Kameno et al. 2001). The latter was argued to cause free–free absorption of the jet-emitted radiation, which is largely seen from VLBI observations at relatively low frequencies (see, for example, Kadler et al. 2004b). Previous studies inferred a free–free absorption coefficient of τ_1GHz ∼ 300–1000 (Kameno et al. 2001; Sawada-Satoh et al. 2008) and τ_8GHz ∼ 3 (Kadler et al. 2004b).

We expand the previously published VLBI results on NGC 1052 by adding new observations at 43 GHz and 86 GHz with the goal of exploring whether the asymmetries between the two jets persist. We also provide further insights into their magnetic fields. The paper is divided as follows. In Sect. 2 we present our data. In Sect. 3 we present our results and focus on the jet collimation and velocity profiles together with the (43–86) GHz spectral index map. In Sect. 4 we discuss the results, with a focus on the time variability of the jet asymmetry and magnetic field strength of the jets. In Sect. 5 we present our conclusions. We assumed a ΛCDM cosmology with H₀ = 71 h km s⁻¹ Mpc⁻¹, Ω_M = 0.27, and Ω_Λ = 0.73 (Komatsu et al. 2009). The luminosity distance of NGC 1052 is D_L = 19.23 Mpc (Tully et al. 2013), and 1 mas corresponds to 0.093 pc and 6134 R_S². As the viewing angle for NGC 1052, we assumed θ > 80°, following Vermeulen et al. (2003), Baczko et al. (2019, 2022).

2. Data analysis

2.1. Original data

We present new observations of NGC 1052 that were performed with the global millimeter-VLBI array (GMVA) (Ros et al. 2024) with a six-month cadence in the period April 2021–April 2022 at a frequency of 86 GHz. The observations were accompanied by observations at 43 GHz in intervals performed by the very long baseline array (VLBA). Unfortunately, we are able to only present one image at 86 GHz, namely the one performed in April 2022, given that a large number of problems afflicted the April 2021 and October 2021 epochs. For a detail description of the problems afflicting the 86 GHz observations alongside the amplitude scaling problems at 43 GHz, we refer to Appendix A.

The data were calibrated using the common astronomy software application (CASA)-based pipeline rPicard (Janssen et al. 2019). In rPicard, several parameters must be specified a priori, such as the signal-to-noise ratio thresholds for the various fringe-fitting stages and the choice of reference antennas. Consequently, to achieve an optimal calibration following this procedure, multiple runs on each dataset were performed with different combinations of the prior parameters. In the final step, the imaging and self-calibration steps on the final calibrated dataset were performed using the standard procedure for VLBI data in difmap (Shepherd et al. 1994).

For the 43 GHz data, due to their relatively high quality, only two runs were necessary in rPicard to achieve science-ready calibrated data; the resultant images are shown in Fig 1.

Fig. 1. 43 GHz VLBA observations of NGC 1052 performed in i) April 2021 (top panel), ii) October 2021 (middle panel), and iii) April 2022 (lower panel). The contours start at 3σ, and the value is reported in Table 1 for each map. The gray ellipse represents the naturally weighted beam.

The situation at 86 GHz was more complicated, and epoch April 2022 required five runs to find a set of parameters that led to a properly calibrated dataset. The final image is presented in Fig. 2: the left panel shows the uniformly weighted image, while the right panel displays the same data convolved with the 43 GHz – April 2022 beam. As mentioned earlier, the low data quality for epochs April and October 2021 prevented the generation of a properly calibrated data set. In the final runs to try to recover as many baselines as possible, the S/N threshold for the science fringe fitting was lowered to 2.9 and, despite that, many antennas had a flag rate of 100%.

Fig. 2. April 2022 86 GHz GMVA image. Left panel: Uniform weighting. Right panel: Image convolved with the 43 GHz April 2022 beam (0.677 mas × 0.235 mas, −11.8°). The contours start at 3σ.

Table 1 summarizes the important parameters of the final maps used for the analysis. In Sect. 3.1 we discuss extensively the quality and validity of our new 86 GHz image.

2.2. Gaussian modelfit components

To compute the jet kinematics (see Sect. 3.2.2), we fit the visibilities with circular Gaussian components following the modelfit task within difmap. The obtained components are reported in Appendix B. As uncertainties, we assume conservative values of 10% error on the flux densities and a 25% error on the size. The uncertainty on the radial distance for each component is calculated by summing in quadrature 25% of their size together with 25% of the circular equivalent beam. To determine whether a certain component is below the resolution limit, and so has to be neglected, we used (Kovalev et al. 2005)

$$\begin{array}{c}\hfill {\theta }_{\lim }=b_{\mathrm{\Theta }}\sqrt{\frac{2\ln 2}{\pi }\ln (\frac{\mathrm{S}/\mathrm{N}}{\mathrm{S}/\mathrm{N}-1})},\end{array}$$(1)

in which b_Θ is the full-width at half maximum (FWHM) of the beam along a certain direction Θ and S/N is the signal-to-noise ratio of each component, calculated as the component flux over the thermal noise in a nearby region which does not include jet emission. The FWHM of the beam is assumed to be the beam size along the position angle of both jets.

As discussed in the previous section, the 86 GHz 2021 observations were not imaged due to their low data quality. We performed a model fit on the limited available visibilities using a single-Gaussian component to model the core, however. The resultant core component parameters are listed in Appendix B, and they are found to be above the resolution limit, according to Eq. (1).

2.3. Map alignment

When comparing jet emission at different frequencies, it is crucial to consider the so-called core-shift effect. Due to the opacity of the nuclear region, the VLBI core, defined as the optically thick region in the jet, moves downstream of the jet with decreasing observing frequency (Blandford & Königl 1979). Therefore, to perform a multifrequency analysis, the radio images need to be aligned to a common origin, which is ideally the position of the central black hole.

The core displacement of the eastern and approaching jet between 1 GHz and 43 GHz has been evaluated in Baczko et al. (2022, see their Table 3), while the shift between the 43 GHz and 86 GHz cores is taken from Baczko et al. (2024). Recent studies (for example, Plavin et al. 2019) suggest the displacement to be time variable, a consequence of the varying underlying jet physical properties. Therefore, in principle, simultaneous data are needed to adequately describe the core shift across the different frequencies. As a consequence, in Sect. 3.3 we infer the 43 GHz core shift from the new observations presented here, and we discuss the impact of the core-shift variability. In Fig. 3, we plot the core sizes of NGC 1052 as a function of the frequency and we perform a fit using the equation z_core = c_s + a ⋅ ν^−1/k_r, in which c_s is the core-shift value at infinite frequency, a the scaling factor, and k_r the inverse of the power-law index. As best-fit parameters, we obtain c_s = 0.05 ± 0.16 mas, a = 14.94 ± 7.67, and k_r = 0.73 ± 0.38. The parameter k_r is of high importance since it reflects the physical properties of the underlying jet. In detail, in the case of a conically expanding jet in equipartition between particle and magnetic field energies, its value is expected to be ∼1 (Blandford & Königl 1979). On the contrary, k_r < 1 is expected for a non-conical and/or accelerating jet (see, for example, Kravchenko et al. 2025).

Fig. 3. Core shift of the eastern jet as a function of the frequency. The orange data points are taken from Baczko et al. (2022), and the 43 GHz blue data point is taken from Baczko et al. (2024). The continuous orange line represents the best fit, and the dashed lines represent the corresponding upper and lower limits. The gray area shows the possible best-fit line within 1σ.

The obtained value of k_r = 0.73 ± 0.38 is consistent with the theoretical expectations (but it is also consistent with one within 1σ). Indeed, only the data point at 1.5 GHz lays in the conical jet (according to Baczko et al. 2024), while the other core shift measurements are within the collimating jet region. We performed a power-law fit that also excluded the point at 1.5 GHz, in order to consider only the purely collimating jet region. We obtained k_r = 0.46 ± 1.00. While the best-fit value is much lower, the high uncertainty makes it comparable with the estimate obtained using all the data points.

At frequencies ν ≲ 43 GHz, namely, when the core region is heavily affected by the torus opacity, Baczko et al. (2022) attributes the first component of the western jet as the core, being the dynamical center. We made the same choice, as shown in Appendix B. On the contrary, at 86 GHz, where the obscuration effects of the torus are expected to be minimal – if not null – the core is associated with the brightest component, which lies on the position of the brightest pixel in our model fitting. Throughout the analysis, the images are centered on the respective core components. In Sect. 3.3 we discuss whether the core shift between the 43 GHz and 86 GHz cores obtained in Baczko et al. (2024) is also valid for the April 2022 data presented here.

3. Results

3.1. 86 GHz observations

Due to the relatively low surface brightness of the extended emission, imaging the jet in NGC 1052 is a non-trivial task. The interpretation of the visibilities, and thus the resulting image, is not straightforward and may be influenced by human biases during the cleaning process. To mitigate this issue, multiple imaging attempts with slightly different approaches in difmap were carried out to obtain the final image shown in Fig. 2. The selected image is the one with the lowest σ_rms across all the imaging attempts and with a jet structure compatible with the 43 GHz simultaneous observations. In Appendix C, we show analternative version of the imaging at 86 GHz to prove that, even by accounting for some differences in the final reconstructed jet emission, the physical discussion and conclusions we extrapolate from the 86 GHz data remain unaltered.

The image chosen for our analysis has a total flux of ∼0.745 Jy with a brightness peak of 0.62 Jy/beam, largely in agreement with the properties of the 86 GHz image shown in Baczko et al. (2024) from 2017 observations. The core flux density is ∼0.650 Jy corresponding to ∼87% of the total flux detected. From the ALMA monitoring, in April 2022 NGC 1052 showed a flux of 0.994 Jy at 91.5 GHz³, ∼33% higher than our obtained VLBI flux. We expect most of this missing flux to come from the extended jet since the lack of short baselines of the GMVA makes it less sensitive to extended jet emission. When performing the spectral analysis (see next section), we consider and discuss the possible implications of the missing flux.

The situation differs for the 2021 epochs. As reported in Appendix B, the April and October observations have central Gaussian components with fluxes of 1.12 Jy and 1.81 Jy, respectively. From the ALMA monitoring, NGC 1052 was observed to have a total flux of 1.43 Jy on 28.03.2021 and 1.214 Jy on 23.10.2021. In the first case, the flux of the core component corresponds to ∼80% of the source flux, leading to a core compactness in agreement with the 86 GHz observation presented here and in Baczko et al. (2024). In the second one, the flux is higher than the actual flux density observed, a sign of improper amplitude calibration. This is not surprising, considering the numerous problems which afflicted the October 2021 session (see Sect. 2.1).

We calculated the brightness temperature of the 86 GHz core components using Eq. (5) in Kadler et al. (2004b). For the October 2021 epoch, we use 80% of the ALMA monitoring flux as the component flux, consistent with the compactness of the April 2021 observation. We obtained brightness temperatures of T_b = (1.25 ± 0.51)×10¹¹ K, T_b = (1.30 ± 0.54)×10¹¹ K, and T_b = (1.30 ± 0.47)×10¹³ K for the April 2021, October 2021, and April 2022 epochs, respectively. Brightness temperatures of ∼10¹¹ K at 86 GHz are consistent with the values reported at the same frequency in Baczko et al. (2024) and are larger than the value at 230 GHz (T_b = 4.7 × 10⁹ K). The exceptionally large T_b obtained in April 2022, which exceeds the inverse Compton limit, is most likely unrealistic. This large value is due to the very small core component, which is most probably caused by improper sampling of the nuclear region.

3.2. Jet structure

The 43 GHz images of NGC 1052 shown in Fig. 1 lack a central, bright core component at the location of the jet nozzle, which, in contrast, is clearly visible in the 86 GHz image (see Fig. 2). This feature, which is unusual for AGN jets observed with VLBI, is consistent with previous 43 GHz observations reported by Baczko et al. (2022). Specifically, the source exhibited a prominent central component at 43 GHz in data obtained between 2005 and 2007, whereas more recent images from 2008–2009 show a progressive disappearance of this bright core. The absence of a bright, central component at 43 GHz in certain epochs is reproduced in the relativistic hydrodynamic simulations of an NGC 1052-like source presented in Saiz-Pérez et al. (2025). In their work, the phenomenon is attributed to free–free absorption caused by the torus and the injection of traveling shocks downstream the jets. Additionally, all three 43 GHz images show a bridge of fainter brightness and apparent thinner width between the eastern jet core and the western jet. This structure is well known at lower frequencies (see, for example, Kadler et al. 2004b) to be a consequence of the impact of the torus, which therefore remains the main candidate at 43 GHz as well. Further evidence in this direction is given by the presence of free–free absorption caused by the ionized torus in the 43 GHz images and its impact on the observed jet morphology. We discuss this in Sects. 3.3 and 4.1.

As seen in Fig. 1, the launched jets show some degree of asymmetries, with high-intensity regions whose positions change across the epochs. Focusing on the eastern jet, in April 2021, the jet shows two bright spots connected by a relatively faint bridge of radio emission. In the consecutive epoch, October 2021, this bridge greatly increased its luminosity, matching that of the two neighboring regions. Finally, in April 2022, we recovered a morphology similar to that observed one year before, with the first bright spot closer to the core remaining relatively steady while the second moved downstream the jet at a distance of ∼1.8 mas in RA The western jet also shows a high level of variability, with two bright regions within 2 mas from the central engine in the October 2021 epoch and only one bright spot in April 2021 and 2022. The variability described could, in principle, be associated with the time-variable torus and the traveling of jet components downstream of both jets (Saiz-Pérez et al. 2025). We analyze and discuss in detail the knots kinematics in Sect. 3.2.2.

In Fig. 4 we present the flux density profiles of the stacked images at 43 GHz and 86 GHz, which are shown in Appendix D. The 43 GHz profiles (blue data points) show asymmetries between the two jets. Within the first milliarcsecond, the western jet appears dimmer with respect to the eastern counterpart, while at a larger distance from the core the two jets switch, with the recessing one showing a region of enhanced brightness at ∼ − 4 mas. In Sect. 4.1 we discuss whether this asymmetry is a consequence of intrinsic differences in the jet launching. On the contrary, the 86 GHz observations reveal a different scenario with the flux density profiles (orange data points) showing an high level of symmetry within ∼0.3 mas.

Fig. 4. Flux density profile for the 43 GHz and 86 GHz stacked images for the eastern and western jet. The double-sided jets are highly symmetrical at 86 GHz, and a dimmer western jet on submilliarcsecond scales is visible at 43 GHz.

3.2.1. Collimation profile

In Fig. 5, we present the collimation profile shown in Baczko et al. (2024) updated with the 43 GHz and 86 GHz stacked jet width. The new data points have been extrapolated using the pixel-by-pixel method presented in Ricci et al. (2022). For the details on the code, we refer to Sect. 2 in the paper. Here we highlight how the code allows us to extrapolate the integrated flux, jet width, and opening angle for each consecutive jet slice perpendicular to the jet direction both in the jet (the eastern one) and in the counter-jet (the western one). The de-projected width was computed as $w=\sqrt{w_m^2-{\theta }^2}$ in which w_m is the width of the fitted Gaussian profile and θ is the FWHM of the beam. For the uncertainties on the jet width, we assume a conservative value of 20% of the extrapolated jet width, which also accounts for the fitting errors, summed in quadrature with 10% of the equivalent circular beam.

Fig. 5. Deconvolved jet width as a function of the distance from the core for the eastern and western jet. This plot, originally published in Baczko et al. (2024), is updated with the new profiles at 43 GHz and 86 GHz we extrapolated. We show three different fits: i) black continuous line, best-fit line for the high frequency data (43, 86, and 230 GHz); ii) red dashed line, best-fit line for all the data; iii) brown dot-dashed line, best-fit line reported in Baczko et al. (2022). In the tables, W0 are the initial jet widths, ku and kd are the upstream and downstream power law indexes, respectively, and zb are the transition distances.

As seen in Fig. 5, the new profiles coincide with the previous estimates of the jet width (shown as gray points). On the one hand, the new 43 GHz jet widths are largely in agreement with the previous data points at the same frequency and match within 1σ. On the other hand, the 86 GHz data points show a slightly thinner jet upstream ∼0.1 mas than the one recovered at 43 GHz (although consistent within 1σ with the previously obtained 43 GHz data points) and are in agreement with the width recovered from the EHT observation at 230 GHz (Baczko et al. 2024). Downstream ∼0.1 mas the 86 GHz data points highly match the ones recovered at 43 GHz. We notice how the jet widths from the modelfit of the core components of the 2021 86 GHz observations are consistent with the measurement obtained from the stacked image at the same frequency (see Appendix D). On submilliarcsecond scales, the jet becomes progressively thinner at increasing observing frequencies, which we suggest to be a consequence of a progressively less important scattering caused by the torus. Indeed, the wider jet at relatively low frequencies (ν ≤ 22 GHz) could be a consequence of the emitted jet radiation being highly scattered by the torus, which leads, in turn, to the cylindrical geometry observed and the deviation from the parabolic expansion, as previously hypothesized in Baczko et al. (2022). On the contrary, at higher frequencies (ν ≥ 43 GHz), the scattering caused by the torus becomes progressively negligible and the true jet width profile is recovered. Overall, the newly added data points confirm the scenario of the jet expanding initially with a parabolic regime.

We fit a broken power-law to the two jets that is described by

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

in which W₀ is the initial jet width, k_u and k_d are the upstream and downstream power law indexes, respectively, s is the hardness of the broken profile, and z_b is the distance at which the break occurs. At first, we fit all the data available for both jets (see Fig. 5). For the eastern jet, we derive an initial jet width of W₀ = 0.23 ± 0.04 mas, with upstream and downstream power-law indexes of k_u = 0.52 ± 0.03 and k_d = 0.75 ± 0.01, respectively. For the western jet, the corresponding values are W₀ = 0.27 ± 0.02 mas, k_u = 0.42 ± 0.02, and k_d = 1.20 ± 0.05. The transition from a parabolic to a conical jet shape is constrained to occur at z_b = 1.2 ± 0.3 mas for the eastern jet and z_b = 2.7 ± 0.2 mas for the western one. The inclusion of new data points allowed us to sample the submilliarcseond scales in more detail, which revealed that the jet geometry steepens more strongly closer to the core than constrained previously. This finding agrees with the results suggested by Baczko et al. (2024), aided by the newly added EHT measurement at 230 GHz. On the contrary, in Baczko et al. (2022), upstream indexes of k_u = 0.22 ± 0.06 and k_u = 0.21 ± 0.05 were derived for the approaching and receding jets, respectively (shown as brown dashed dotted lines in Fig. 5). For both the eastern and western jets, both the transition distance and the conical jet geometry obtained here are largely in agreement with those of Baczko et al. (2022). Additionally, we notice that along the western jet, we detected a longer transition distance which leads to a jet width expansion rate steeper than the one inferred for the eastern jet. Interestingly, the theoretical work presented in Saiz-Pérez et al. (2025) shows a very similar behavior: a jet break occurring farther from the core results in a larger downstream expansion rate. We also performed a separate fit using only the high-frequency data presented here (43 GHz and 86 GHz) together with the 230 GHz point (black continuous line in Fig. 5), which yields consistent results with the previous fit. The best-fit parameters for the different scenarios are shown in Appendix B.

Overall, we propose the fit to all data points (dashed red line in Fig. 5) as the best representation for both jet geometries. For the eastern jet, it is the only fit that adequately captures the parabolic expansion on submilliarcsecond scales and the quasi-conical expansion downstream of the transition distance, while for the western, it largely matches with the fit proposed in Baczko et al. (2022) while describing the parabolic region in a more accurate way. Accordingly, throughout the rest of the paper, we will adopt the following parameters. For the eastern jet: k_u = 0.52 ± 0.03, k_d = 0.75 ± 0.01, and z_b = 1.2 ± 0.3 mas; while for the western one: k_u = 0.42 ± 0.02, k_d = 1.20 ± 0.05, and z_b = 2.7 ± 0.2 mas. The expansion rates we extrapolate are in agreement with the theoretical predictions from Saiz-Pérez et al. (2025). Overall, our findings indicate that on submilliarcsecond scales the jets appear symmetric, showing similar expansion rates and initial jet widths, while asymmetries arise further downstream, with the western jet showing a larger transition distance and a much steeper conical profile than the eastern one.

3.2.2. Jet speed

To calculate the jet speed, we cross-correlated the modelfit components across the three different epochs at 43 GHz (see Appendix B for the list of the jet components). In Fig. 6, we present the displacement of each component for the eastern (left panel) and the western jet (right panel). For each individual component, we performed a linear fit of the form z_c = a + μ t, where a is the initial position, z_c is the distance from the core, μ is the jet speed in units of milliarcseconds/years and t is the time in years. The best-fit value of μ for each component, expressed as a fraction of the speed of light, is reported in Table 2. In both jets, we detect the speed of the components to increase with the distance from the core, as observed for the majority of AGN jets (see, for example Lister et al. 2021). Specifically, the components reach the speed β = μ/c ∼ 0.25 on both jets at ∼4 mas from the core. The only exception is represented by component E5, which shows an apparent speed of β = 0.41 ± 0.09. It is, however, possible that what we call component E5 detected in epoch 2022 is actually a combination of components E5 and E6, and it is located halfway between the two. Indeed, as visible from Fig. 1, bottom panel, the easternmost jet region does not show a clear double-component structure, visible on the contrary in the two previous epochs, likely due to a consequence of a lower S/N in the region or evolving jet morphology. In this scenario, the higher velocity of component E5 would be a consequence of an improper cross-identification and not representative of the actual speed at such distance. We note that the first two components in the eastern jet are not moving and can be associated with recollimation shocks, in agreement with the findings of Saiz-Pérez et al. (2025).

Fig. 6. Distance from the core across the three epochs between April 2021 and April 2022 for the Gaussian components reported in Appendix B. Left panel: Eastern jet. Right panel: Western jet. The differently colored continuous lines show the best-fit curves for the different components.

The jet speeds we inferred are lower than those inferred by Baczko et al. (2019) between 2005 and 2009, in which average speeds of 0.343 ± 0.037 c and 0.529 ± 0.038 c have been inferred for the western and eastern jet, respectively. In contrast, we inferred average speeds of 0.19 ± 0.01 c and 0.26 ± 0.01 c, respectively. Additionally, the previously detected asymmetry in the jet speed between the approaching and receding jet is not recovered between 2021 and 2022. Overall, this analysis might indicate how the components in both jets are moving at lower and more symmetric speeds with respect to the past, even if the result might be biased by the lower number of epochs employed for the analysis.

3.3. Spectral index map

In this section, we present the first (43–86) GHz spectral index map of NGC 1052 (we use the formalism S_ν ∝ ν^α). To obtain it, we constrained the visibilities in both datasets to the same (u, v)-range, we convolved the final 86 GHz image with the beam of the 43 GHz one (see Table 1), and applied a flux cutoff at 5σ. The pixel size was set to 0.041 mas. To properly compute the spectral index map, the correct displacement between the 43 GHz and 86 GHz images in April 2022 had to be determined. Initially, we aligned the two frequencies using the core shift obtained in Baczko et al. (2024). As shown in Appendix E, the two maps appear not to be correctly aligned. As a consequence, we attempted to align the images at the two frequencies by employing an automatic 2D cross-correlation analysis, similar to what has been done in Fromm et al. (2013), Baczko et al. (2022). It was impossible to obtain meaningful results, however, because the jet emission recovered at 86 GHz is poorly extended and weak. To perform the alignment, we had to rely on the high-brightness region at a distance of ∼ − 1 mas along the western jet found at both frequencies. The shift that allows to align the contours at the two frequencies is (x, y) = (1 ± 1, 0 ± 1) pixels = (0.041 ± 0.041, 0.000 ± 0.041) mas; we show the correspondent spectral index map in Fig. 7. The shift along the RA direction is positive, indicating that the 43 GHz map is shifted west. This is unexpected, as the radio core at progressively lower frequencies is typically expected to move downstream along the approaching jet, namely the eastern one. A similar scenario was observed by Baczko et al. (2022) when they calculated the alignment between 15 GHz and 22 GHz, however. A possible explanation might be that the shift, being minimal at these frequencies, could oscillate due to different conditions in the underlying jet. As uncertainty in our alignment method, we adopt a displacement of one pixel, corresponding to the minimum calculable shift. In Appendix E we present alternative spectral index maps obtained by applying all possible shift combinations within a ±1 pixels range. We highlight that the main physical properties discussed in the consecutive paragraphs remain unchanged across the different spectral index maps, being the values of α across the different regions consistent within the errors.

Fig. 7. 43–86 GHz spectral index map of NGC 1052 obtained using the core shift proposed here. The black contours show the 86 GHz emission, and the green contours show the 43 GHz image.

We find a highly inverted core region, with α_thick ∼ (3.0–3.3) with σ = 0.16, according to the different displacements. The measured spectral index exceeds the theoretical limit expected from synchrotron self-absorption (α_thick = 2.5, under the assumption of the radiation emitted from a single emitting region), further suggesting that free–free absorption from the ionized torus plays a role at 43 GHz. This value of the spectral index is in agreement Baczko et al. (2024), who report a spectral index of α_thick = 3.27 between 43 GHz and 86 GHz for the core component. Detecting free–free absorption at such high frequencies is relatively surprising, and we further discuss it in Sect. 4.1.

Along the eastern jet, the spectral index goes down to ∼ − 2 outside of the core with the first component located at ∼0.5 mas, which shows an average spectral index of −1.9 and a gradient perpendicular to the jet direction. The latter could be a consequence of an improper alignment for the eastern jet or the bulk of the jet emission at 86 GHz having a slightly different position angle than the one at 43 GHz. We notice that with different alignments (see, for example, Fig. E.2 first row, left column), such a gradient can strongly attenuate. The two more downstream regions inferred in the spectral index map have average index values of −1.97 ± 0.30 and −3.89 ± 0.29, respectively. Along the western jet, the spectral index goes down to ∼ − 1.4 outside of the core with the inverted spectrum region being more extended than the eastern counterpart, which we associate with the free-free absorption caused by the torus since it coincides with the bridge of absorbed emission mentioned in the previous section. The component recovered at ∼ − 1 mas, shows an average spectral index of α_thin = −0.98 ± 0.24 and a gradient parallel to jet direction. This is likely a consequence of a region in the 43 GHz jet with reduced brightness in between two bright spots (see Fig. 1, lower panel). Concerning the very steep spectral index values obtained (down to ∼ − 3.9), while steep spectra on sub-parsec scales are recovered on similar distances in other sources (down to ∼ − 2.5, see, for example, Ro et al. 2023; Ricci et al. 2025), in NGC 1052 this may instead be a consequence of missing flux at 86 GHz, which, as mentioned in Sect. 3.1, is expected to originate from the extended emission. In the scenario where a fraction of the missing flux comes from the core, the spectrum of the optically thick region would be even larger than α_thick ∼ (3.0–3.3), strengthening the idea of the radiation at 43 GHz being absorbed by the ionized torus. Additionally, in Appendix E we report and discuss the (43–86) GHz spectral index map obtained employing the alternative 86 GHz image shown in Appendix C.

Finally, using the (43–86) GHz displacement of (x, y) = (0.041 ± 0.041, 0.0 ± 0.041) mas, the resulting power-law index for the core shift becomes k_r = 0.87 ± 0.10, consistent with the value of k_r = 0.73 ± 0.38 inferred in Sect. 2.3. We will use k_r = 0.87 ± 0.10 to extrapolate the magnetic field strength in Sect. 4.2.

4. Discussion

4.1. Constraints on the production of the symmetric jet

The flux density profiles of the eastern and western jets shown in Fig. 4, exhibit distinct behaviors at the two frequencies. At 86 GHz, the jets appear symmetric within ∼0.3 mas from the core, with fluxes consistent within uncertainties at each given distance. At 43 GHz, the flux density profiles differ, with the eastern jet being brighter than the western one within the inner ∼1 mas from the core. To investigate whether the observed asymmetry could be attributed to absorption by the ionized torus, we calculated the required optical depth for free–free absorption and the corresponding column density of the ionized material. These values are then compared with those reported in the literature. We work under the assumption that the decrease in brightness in the western jet is only due to the free–free absorption from the torus, and that the jets are intrinsically symmetric (as suggested by the 86 GHz profiles). Between (0.3–0.8) mas, namely the region in which we see the reduced brightness, the average flux densities are 27 mJy and 12 mJy for the eastern and western jet, respectively. Employing S_ν^obs = S_ν^inte^−τ_f, we constrain an optical depth for the free–free absorption of τ_f = 0.8. The free–free optical depth is expected to evolve with the frequency as ∝ν^−2.1 (see, for example, Kadler et al. 2004b). In Kadler et al. (2004b) the authors inferred τ_f ∼ 3.0 at 8.4 GHz, which would lead to τ_f ∼ 0.09 at 43 GHz, one order of magnitude lower than the one we infer here. Such a discrepancy could be attributed to the absorption being generated in two different regions at the two frequencies. In detail, the lower frequency emission is produced in a more extended region, which could suffer absorption by a diffuse screen, whereas the absorption of the more compact area at 43 GHz could be generated by some denser clump toward our line of sight, as justified later in the text. Additionally, the value of τ_f ∼ 3.0 at 8.4 GHz is in principle an average of a much larger region with respect to the one inferred at 43 GHz, due to the reduced resolution at the lower frequency. An alternative explanation could, in principle, involve temporal variations in the absorption caused by the torus, a plausible scenario in the case of a clumpy structure. A clumpy torus has indeed been suggested to surround NGC 1052 from ALMA observations (Kameno et al. 2020).

Thanks to the knowledge of τ_f, we can calculate the column density employing the equations shown in Lobanov (1998),

$$\begin{array}{c}\hfill {\tau }_{\mathrm{f}}=5\times {10}^{16}\overline{g}{n}_{\mathrm{e}}{n}_{\mathrm{i}}[\frac{T}{\mathrm{K}}{]}^{-1.5}[\frac{\nu }{\mathrm{Hz}}{]}^{-2}[\frac{L}{\mathrm{pc}}]{\mathrm{cm}}^{-3},\end{array}$$(3)

in which n_e = n_i are the numerical density of the electron and ions, which we assume to coincide, T is the temperature of the torus, L is the path length of the light, and $\overline{g}$ is the free–free Gauge factor. We approximate the latter as $\overline{g}=1$, according to Table 5 shown in van Hoof et al. (2014). As parameters, we use T = 10⁴ K, L = 0.3 pc (following Kadler et al. 2004b), τ_f = 0.8, and ν = 43 GHz, obtaining n_e = 3.1 × 10⁵ cm⁻³ and, consequently, a column density of 2.9 × 10²³ cm⁻². Temperatures up to T = 10⁴ K are expected to be high enough to sublimate the dust in the torus, however, which argues against the idea of a clumpy structure (see, for example, Kishimoto et al. 2007). By assuming a temperature one order of magnitude lower (expected to be smaller than the temperature needed for the sublimation), namely T = 10³ K, we obtain a column density of 5.2 × 10²² cm⁻². Both the column densities constrained are in agreement (although not considering different possible values for the ionization and metallicity of the torus) with X-ray observations for the hydrogen column density n_H ∼ 10²² cm⁻² − ∼10²⁴ cm⁻² (Weaver et al. 1999; Guainazzi & Antonelli 1999; Kadler et al. 2004a), supporting the idea of absorbed emission at 43 GHz.

Observing free–free absorption at the high frequencies of 43 GHz is somewhat unexpected. Nonetheless, in support of this scenario we found: i) the inferred column density is in agreement with the one independently determined with X-ray observations; ii) the (43–86) GHz spectral index map presents values of α_thick > 2.5 indicative of free–free absorption. To explain the reduced brightness along the western jet, alternatives include Doppler de-boosting and/or an intrinsic asymmetric jet production. On the one hand, the large viewing angle of NGC 1052 (θ ≥ 80°) and the low velocities inferred at such distance (see Sect. 3.2.2) cannot lead to such Doppler de-boosting. On the other hand, the hypothesis of asymmetric jet production does not match with the symmetry between the two jets seen at 86 GHz.

In conclusion, we propose that the jets in NGC 1052 are produced intrinsically with a high level of symmetry. Asymmetries in the jet morphologies, as the one observed in Fig. 1 and in Baczko et al. (2019), are suggested to be a consequence of the interplay between a clumpy, time-variable torus (as shown in Kameno et al. 2020), predominantly affecting the western jet, and observational effects. The latter are theorized in Saiz-Pérez et al. (2025), in which the authors demonstrate that intrinsically symmetric jets can appear asymmetric at certain epochs due to the downstream traveling of perturbations injected at the jet nozzle. These apparent asymmetries arise as a consequence of relativistic effects, such as differences in light travel time and viewing angle, which affect the apparent timing and brightness of the features in the two jets.

4.2. Magnetic fields on compact scales

In this section, we use the observational parameters reported in Sect. 3 and interpret them in the context of theoretical models to investigate the magnetic field properties of the jets and assess whether the accretion disk in NGC 1052 has reached a MAD state. To do this, we follow an approach similar to the one shown in Ricci et al. (2022).

4.2.1. Saturation of the magnetic field

First, we calculated the magnetic field strength along the jet, extrapolated it to its base, and then compared it with the expected magnetic field needed to saturate the disk and form a MAD. To calculate the magnetic field, we employed the equations reported in Lobanov (1998), Ricci et al. (2022), Nokhrina & Pushkarev (2024), in which the field strength on (sub-)parsec scales is calculated starting from the knowledge of the core-shift through the parameter

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_{r\nu }=4.85\times {10}^{-9}\frac{\mathrm{\Delta }{r}_{\mathrm{mas}}{D}_{\mathrm{L}}}{{(1+z)}^2}(\frac{{\nu }_1^{1/k_r}{\nu }_2^{1/k_r}}{{\nu }_2^{1/k_r}-{\nu }_1^{1/k_r}})\mathrm{pc}{\mathrm{GHz}}^{1/k_r},\end{array}$$(4)

in which Δr_mas is the core shift between the two frequencies ν₁ and ν₂, and D_L is the luminosity distance of the source.

By assuming a jet expanding with a conical shape and k_r = 1, the magnetic field at one parsec from the core can be computed as (Lobanov 1998)

$$\begin{array}{c}\hfill B_1=0.025[\frac{{\mathrm{\Omega }}_{r\nu }^3{(1+z)}^3}{{\delta }^2\varphi {\sin }^2\theta }{]}^{1/4}\mathrm{G},\end{array}$$(5)

where δ is the Doppler factor, and ϕ is the intrinsic jet half-opening angle.

To compute Ω_rν, we employed the core shift between 1 and 5 GHz, as these frequencies correspond to a displacement close to the distance of one parsec. Using the values reported in Baczko et al. (2022), we obtain Ω_rν = 1.01 ± 0.10 pc GHz. For the jet parameters, we consider average values between the eastern and western jets. We assume a jet speed of β = 0.40  ±  0.20, consistent with the average values extrapolated by Baczko et al. (2016) and the speed of components closer to one parsec we observed, as shown in Sect. 3.2.2. The jet half-opening angle at one parsec is ϕ = (4.0 ± 2.0)° (Baczko et al. 2022). Overall, with such parameters, we obtained an average magnetic field strength of B₁ = (0.018 ± 0.008) G at one parsec on both jets.

We showed in Sects. 2.3 and 3.2.1 that k_r is lower than one, however, and the jets of NGC 1052 show a quasi-parabolic geometry up to a few milliarcseconds. Therefore, for a more precise evaluation of the magnetic fields we need to drop the assumption of a conical jet used to compute Eq. (5). To do so, we used the model shown in Ricci et al. (2022), Nokhrina & Pushkarev (2024), which allows us to compute the magnetic field strength starting from the core-shift parameters in a generic jet geometry. There are three different cases, depending on the interplay between the viewing angle and the inverse of the Lorentz factor, namely θ ∼ Γ⁻¹, θ ≤ Γ⁻¹, and θ ≥ Γ⁻¹. For NGC 1052 the first and third scenarios are the most plausible ones, since θ ≥ 80° ≥1.1 rad and Γ⁻¹ ∼ 0.92 (using an average β ∼ 0.40, see Table 2). In the first scenario, θ ∼ Γ⁻¹,

$$\begin{array}{c}\hfill B_z=0.025[\frac{{\mathrm{\Omega }}_{r\nu }^{6\psi }{(1+z)}^3}{{\mathsf{z}}^{\mathsf{6}\psi }{\mathsf{r}}_{\mathsf{z}}{\delta }^{\mathsf{2}}{\sin }^{\mathsf{6}\psi -\mathsf{1}}\theta }{]}^{1/4}\mathrm{G},\end{array}$$(6)

while in the third, θ ≥ Γ⁻¹,

$$\begin{array}{c}\hfill B_z=0.025[\frac{{\mathrm{\Omega }}_{r\nu }^{8\psi }{(1+z)}^3{(1-cos\theta )}^2}{{\mathsf{z}}^{\mathsf{8}\psi }{\mathsf{r}}_{\mathsf{z}}{\sin }^{\mathsf{8}\psi -\mathsf{1}}\theta }{\left(\frac{r_{\mathsf{z}}}{R_L}\right)}^2{]}^{1/4}\mathrm{G}.\end{array}$$(7)

In these equations, r_z is the jet radius (in parsec) at a distance z (in parsec) from the radio core, ψ is the power-law index of the jet expansion profile (d = z^ψ), and R_L is the light cylinder radius (in parsec). The latter is the intrinsic length scale in MHD models and is expressed as R_L = r_br/d_*(σ_M) where r_br is the jet radius at the transition distance and d_* is the non-dimensional jet radius whose value depends on the Michel magnetization σ_M (see Beskin et al. 2017, and references therein).

We computed Equations (6) and (7) for the transition distance, for which we assume the average between the two jets according to Sect. 3.2.1, namely z_br = (0.18 ± 0.02) pc. The corresponding average jet radius is r_br = (0.011 ± 0.003) pc. As a power-law index for the jet expanding profile, we assume the average between the upstream values in the two jets, namely ψ = 0.47 ± 0.04. The Michel magnetization parametrizes the terminal Lorentz factor reached by the jet flow when all the electromagnetic flux is transferred to the particles (see, for example Nokhrina et al. 2015, and references therein). The highest speed observed in NGC 1052 is β ∼ 0.65 (Baczko et al. 2016, 2019) leading to σ_M = Γ ∼ 1.32. In Nokhrina et al. (2019), the different values of d_*(σ_M) are given starting from a magnetization of five, which we therefore assume as our lower limit. Because of the large jet viewing angle, however, the radio emission detected from NGC 1052 is most likely associated with the slower sheath, prohibiting us to detect the faster central spine. NGC 1052 is a low-power radio galaxy, and its aligned counterpart would correspond to a BL Lac object, with expected bulk Lorentz factors of Γ ∼ 10 (see, for example, Hovatta et al. 2009), leading to a magnetization σ_M = 10. As a consequence, we calculated Eq. (7) using values of d_*(σ_M) assuming magnetizations of five and ten, and then averaging the result. We constrain a magnetic field strength at the transition distance of B = 0.27 ± 0.05 G and B = 2.35 ± 0.59 G, using Equations (6) and (7), respectively. While the former is in agreement with the results from Eq. (5), the latter is roughly one order of magnitude higher. The obtained magnetic field strengths are shown in Fig. 8.

Fig. 8. Average spatial magnetic field evolution for the two jets in NGC 1052. The horizontal black line represents the magnetic field strength estimates Baczko et al. (2016) at 1 RS, and the vertical green, blue, and orange lines show the values (with their uncertainties) from Equations (5), (6), and (7), respectively. The dashed lines show the extrapolation assuming b = 0.47 and b = 0.94 (B ∝ z−b) for the parabolic region and b = 1.0 and b = 2.0 for the conical region. The red triangle is a lower limit estimate for the magnetic field from Baczko et al. (2016), and the red dot is taken from Baczko et al. (2024). The red box represents the magnetic field that is expected to be needed to saturate the accretion disk.

The next step is to extrapolate the magnetic field strength computed downstream of the jets up to their onset. To do that, we need to make assumptions on the radial magnetic field dependence, namely b in B(z)∝z^−b. When the magnetic fields in jets are dominated by their toroidal component, b is equal to the inverse of the power-law index for the jet width, namely −ψ. On the contrary, when the magnetic fields are dominated by their poloidal component, b = −2ψ⁴. Having no constraints on the magnetic field geometry, we extrapolate the magnetic field strength considering both cases using ψ = 0.47. To include all possible magnetic field values and account for uncertainties, we extrapolate the lower limits of the magnetic field estimates (value minus the error) with b = 0.47 and the upper limits with b = 0.94 (b = 1 and b = 2 for the conical region, respectively).

In the following, we compare our estimated B with previous works. Namely, in Baczko et al. (2016) the authors inferred a lower limit on the magnetic field strength of 6.7 G at 1.5 mpc = 99 R_S by means of synchrotron cooling arguments which leads to extrapolate values of 360 G ≤ B ≤ 6.9 × 10⁴ G at 1 R_S. In their work, they assumed the decrease in radio emission with the distance from the core at 86 GHz only to the cooling, without effects from the free–free absorption. This hypothesis is supported by our proposition here. Additionally, in Baczko et al. (2024) the authors constrained a magnetic field strength of B ∼ 1.25 G at 43 μas = 264 R_S by fitting the synchrotron-self absorbed spectrum of the core component between 22 GHz to 230 GHz. The literature data points are consistent with the magnetic field strengths calculated using Eqs. (5) and (6) when evolving with b ≥ 0.47, with the values at 1 R_S in agreement with our estimates with a magnetic field geometry of b ∼ 1.0.

Merging all the information together, as shown in Fig. 8 a magnetic field evolving with an index 0.47 ≤ b ≤ 0.94 is expected, since none of the magnetic field estimates along the jet match the values at 1 R_S when evolving with b ∼ 0.47. The situation is more complicated, however. As shown in Fig. 5 on scales smaller than 10⁴ R_S, the parabolic profiles inferred are the average jet geometries, with the jet alternating regions of steeper and flatter expansion rates. Since the magnetic field evolution follows the same trend (b ∝ −ψ), changes in the jet width profile reflect into changes in the magnetic field one. In particular, we observe regions where the jet expansion profile nearly saturates (ψ ∼ 0), such as near the jet injection point, where, consequently, the magnetic field strength remains approximately constant as well. As a consequence, to reach magnetic field strengths in agreement with the findings of Baczko et al. (2016) at 1 R_S, when the jet is expanding the magnetic field need to evolve with an index close to one. A magnetic field trend with the distance b ∼ 1 has been recovered for M 87 and NGC 315 as well on similar scales (see Ro et al. 2023; Ricci et al. 2025, respectively) and is in agreement with previous findings for the source (Baczko et al. 2016). Values of b ∼ 1 can be a consequence of a purely toroidal magnetic field being highly dissipated to accelerate the jet or a magnetic field with an important contribution from its poloidal component leading to a helical field geometry up to scales of ∼10⁴ R_S at least.

Finally, we can compare the magnetic field strength in the nuclear region with the magnetic field needed to saturate the disk (precisely, magnetic force ≥ gravitational force) and form a MAD. To compute the saturation fields, we used Eq. (17) in Ricci et al. (2022). For the mass accretion rate ($\dot{M}=L_{\mathrm{acc}}/\eta c^2$), we assume the value of L_acc = 10^41.24 erg s⁻¹ from X-ray measurements (Woo & Urry 2002; González-Martín et al. 2009). As radiative efficiencies, we assume η = 0.001, 0.01, and 0.4. On the one hand, NGC 1052 has been proposed to have an ADAF, which is a type of disk with efficiencies lower than the standard thin accretion disks (Yuan & Narayan 2014) and that can go down to η = 0.001, 0.01 (Narayan & Yi 1995; Mahadevan 1997). On the other hand, η = 0.4 has been assumed for NGC 1052 in sample studies including it, and we consider it for consistency (Zamaninasab et al. 2014). Finally, the radial velocity of the gas is assumed to be closer to the speed of light and the scale height of the disk h/r = 0.25 ± 0.10. This value and uncertainty is chosen to include the different possible black hole spin scenarios, as shown in Fig. 7 in Narayan et al. (2022). Considering all the possible sources of uncertainties and the different values of η, we find the saturation magnetic field to vary between 80 G ≤ B ≤ 2438 G at 1 R_S and 3 G ≤ B ≤ 77 G at 10 R_S, chosen as an arbitrary upper limit on the putative magnetosphere. Such ranges are shown by the red box in Fig. 8. Following this calculation, the magnetic field strength required to saturate the disk and reach the MAD state is achieved when b ≳ 0.47, as already suggested for NGC 1052 earlier in this section. This analysis provides a very first indication that the accretion disk in NGC 1052 might have reached a MAD state. Overall, our extrapolated magnetic field strength with b ∼ 1 would lead to values at the jet base in the range of ∼10² − 10⁴ G, consistent with the findings for other radio galaxies, such as NGC 315 (Ricci et al. 2022, 2025) and M 87 (Ro et al. 2023). In the scenario, the real magnetic field strengths at the jet base are in the range of tens of Gauss, however, such as those that were proposed for M 87 from EHT observations (Event Horizon Telescope Collaboration 2021), values of b ∼ 0.50 would best describe the field strength extrapolation. In this case, either the magnetic field evolves up to 1 R_S with a purely toroidal configuration or the fields saturate more downstream the jet, forming a magnetosphere on scales of tens of R_S, as shown in Fig. 8.

4.2.2. Magnetic flux

An additional approach to assess whether the accretion disk in NGC 1052 has reached a MAD state involves comparing the magnetic flux along the jet with the one expected from a MAD. This comparison can be done following the work by Zamaninasab et al. (2014). In particular, under the assumption that the gravitational potential of the black hole is the main regularizer of the width of the optical broad lines coming from the clouds surrounding the central region, the expected MAD flux is

$$\begin{array}{c}\hfill F_{\mathrm{MAD}}=2.4\times {10}^{34}[\frac{\eta }{0.4}{]}^{-1/2}[\frac{M_{\mathrm{BH}}}{{10}^9{M}_{\odot }}][\frac{L_{\mathrm{acc}}}{1.26\times {10}^{47}\mathrm{erg}/\mathrm{s}}{]}^{1/2}\mathrm{G}{\mathrm{cm}}^2.\end{array}$$(8)

We use the same parameters as for calculating the saturation fields, namely L_acc = 10^41.24 erg s⁻¹ and η = 0.001, 0.01, and 0.4 obtaining, respectively F_MAD = 7.45 × 10³¹ G cm², F_MAD = 2.35 × 10³¹ G cm², and F_MAD = 3.72 × 10³⁰ G cm².

To calculate the poloidal magnetic flux in the jet, we used Eq. (1) of Zamaninasab et al. (2014)

$$\begin{array}{c}\hfill F_{\mathrm{jet}}=1.2\times {10}^{34}f(a_{\ast })\mathrm{\Gamma }\varphi [\frac{M_{\mathrm{BH}}}{{10}^9{M}_{\odot }}][\frac{B_1}{1\mathrm{G}}]\mathrm{G}{\mathrm{cm}}^2,\end{array}$$(9)

in which f(a_*) = 1/a_*[1 + (1 − a²)], and a_* is the dimensionless spin of the central black hole. This equation is based on the assumption of a jet expanding conically and allows for calculating the magnetic field flux at one parsec downstream of the jet. While we know that the assumption of a conical jet does not hold for NGC 1052 (it does at one parsec, but not upstream, see Sect. 3.2.1), we use this equation to extrapolate an order of magnitude for the magnetic flux. Being the spin of the black hole unknown, we use the two opposite scenarios of a slowly and highly rotating object, namely a_* = 0.10 and a_* = 0.99, obtaining F_jet = (3.01 ± 1.32)×10³² G cm² and F_jet = (1.74 ± 0.8)×10³¹ G cm², respectively.

Assuming magnetic flux conservation from one parsec up to the jet injection point, the calculated flux via Eq. (9) is consistent with the value expected in the MAD scenario under the conditions of low radiative efficiency (in agreement with the disk being an ADAF) and a rapidly spinning central supermassive black hole. Interestingly, the nuclear region of 3C 84 is also compatible with this model, perhaps indicative of such conditions being prevalent in nearby radio galaxies (Paraschos et al. 2023, 2024).

5. Summary

We have analyzed new 43 GHz and 86 GHz observations of NGC 1052 to investigate the symmetry of its jets and their magnetic field properties. Our results are summarized below.

• We obtained three new images at 43 GHz. Their morphology is consistent with the most recent epochs presented by Baczko et al. (2019) and is characterized by the absence of a bright core component that corresponds to the jet-launching region. Additionally, the images reveal asymmetries between the eastern and western jets and a time-variable morphology that is likely associated with the passage of components downstream in the two jets.

• The 2021 observations encountered various issues. We were therefore only able to acquire a single new image at 86 GHz. Based on the stacked 86 GHz image and by combining our new data with the map published by Baczko et al. (2024), we derived the flux density profile of the two jets. The profiles are highly symmetrical.

• We updated the jet width profile published by Baczko et al. (2022, 2024) with the results from the new observations we presented here. We infer symmetric, parabolic jets on submilliarcsecond scales. The eastern jet transitions to a quasi-conical shape at z_b  =  1.2 ± 0.3 mas, and the western jet transitions at z_b  =  2.7  ±  0.2 mas. The jets are asymmetrical downstream. The geometry of the receding jet is steeper.

• The (43–86) GHz spectral index map indicates that the core region is highly inverted, with a spectral index as high as α_thick ∼ 3.3. This suggests that free–free absorption still affects the emission at 43 GHz. This scenario is further supported by the flux density profile at 43 GHz. The brightness of the western jet is approximately three times lower than that of its eastern jet counterpart at the same distance.

• We propose that the bulk flows are overall intrinsically symmetric when launched. The asymmetries observed at 43 GHz are likely a combination of the free–free absorption exercised by the clumpy torus on the western jet and the propagation of jet components downstream of the two jets. This scenario is consistent with the numerical simulations presented by Saiz-Pérez et al. (2025) for NGC 1052: time-dependent asymmetries might arise from symmetrically launched jets caused by the passage of the jet perturbations.

• Finally, we explored the magnetic field strength along the jets in NGC 1052 by employing our knowledge of the core shift and the extent of the collimation region in the two jets. Our analysis is consistent with previous estimates of the magnetic field. This field strength is expected to be strong enough to saturate the central region, implying the possible presence of a MAD. This scenario was further examined by analyzing the magnetic flux at the jet base. The magnetic flux expected in a MAD scenario agrees with the jet flux in the scenario of low radiative efficiency and a highly spinning central supermassive black hole.

In conclusion, the new datasets we presented provide important insights into the jet properties of NGC 1052, from the role of the absorbing torus to the intrinsic jet symmetry and their magnetic field content. Future multifrequency and simultaneous observations will be essential to further constrain these properties and explore the proposed scenarios in greater detail.

Acknowledgments

We would like to thank the referee for the meaningful comments, which improved the manuscript. This research is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – project number 443220636 and by the European Research Council advanced grant “M2FINDERS – Mapping Magnetic Fields with INterferometry Down to Event hoRizon Scales” (Grant No. 101018682). SdP gratefully acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 789410, PI: S. Aalto). MP acknowledges support from MICIU/AEI/10.13039/501100011033 and FEDER, UE, via the grant PID2022-136828NB-C43, from the Generalitat Valenciana through grant CIPROM/2022/49, and from the Astrophysics and High Energy Physics project program supported by the Spanish Ministry of Science and Generalitat Valenciana with funding from European Union NextGenerationEU (PRTR-C17.I1) through grant ASFAE/2022/005. This research has made use of data obtained with the Global Millimeter VLBI Array (GMVA), which consists of telescopes operated by the MPIfR, IRAM, Onsala, Metsahovi, Yebes, the Korean VLBI Network, the Greenland Telescope, the Green Bank Observatory and the Very Long Baseline Array (VLBA). The VLBA and the GBT are facilities of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. The data were correlated at the correlator of the MPIfR in Bonn, Germany. Partly based on observations with the 100-m telescope of the MPIfR (Max-Planck-Institut für Radioastronomie) at Effelsberg. This work made use of Astropy (http://www.astropy.org): a community-developed core Python package and an ecosystem of tools and resources for astronomy (Astropy Collaboration 2013, 2018, 2022).

References

Appendix A: Problems afflicting the observations

Concerning the 3 mm observations, in April 2021, Metsahovi and Pico Veleta were partially afflicted by poor weather, Yebes participated in only a few scans and fringes were detected only for one of them, Green Bank had pointing problems, and the VLBA antennas suffered from problems at the 3 mm receiver ⁵. The October 2021 session suffered from poor weather conditions at several stations (mostly Europeans, such as Effelsberg, Onsala, and Metsahovi), no fringes at Yebes, the data from Green Bank were corrupted and so discarded, and an unusually large number of scans with high system temperature (> 1000 K) at the VLBA stations⁶. The aforementioned problems, combined with the relatively low flux density of NGC 1052 and consequent high uncertainties in the final calibrated data, did not allow us to obtain reliable 86 GHz images from those two observing sessions.

At 43 GHz, the observations performed in April and October 2021 suffered from flux scaling problems (the same problems were reported in Ricci et al. 2025). In April 2021, the output from the correlator at 7 mm contains all autocorrelation records twice, leading to a flux density lower than expected. To account for this issue, the baseline-based fluxes were multiplied by a factor of $\sqrt{2}$, following the suggestions provided by the correlator team. For the October 2021 observations, the gain curves originally provided from the VLBA team were corrupted due to issues with the focus and rotation at the VLBA stations ⁷. This problem was fixed by applying the updated gain curves provided by the VLBA staff to the data.

Appendix B: Modelfit components and best-fit tables

In this section, we report the main parameters of the modelfit components in Table B.1 together with the images with the modelfit component overlaid in Fig. B.1. Additionally, we report the best-fit values (Table B.2) for the different fits shown in Fig. 5.

Fig. B.1. Modelfit components (in red) overlaid to the images presented in this paper. Upper panels: 43 GHz April 2021 (left), 43 GHz October 2021 (right). Lower panels: 43 GHz April 2022 (left), 86 GHz April 2022 (right).

Appendix C: 86 GHz alternative image

In Fig. C.1, we report a 3 mm April 2022 image obtained through an independent imaging run. The main properties of the resultant image remain largely similar of the one shown in Sect. 3.1, namely a total flux of 779 Jy, a brightness peak of 0.66 Jy/beam, and a rms of 1.023 mJy/beam. In this case, the image shows an extended structure within ∼0.6 mas on both sides which is highly symmetric, with one bright component aligned along the same position angle. Additionally, we note how the component around −1 mas in the western jet is not recovered. We produced an 86 GHz stacked image using the alternative result shown in this section and derived the jet width profile following the procedure described in Sect. 3.2.1. In Fig. C.2 we show the difference between the obtained profile and the one employed in Fig. 5. The difference between the vast majority of the data points is consistent with zero within 1σ, with the exception of the last four points in the western jet, which appear to be slightly thinner and consistent with zero within 2σ. Overall, this additional image proves that, despite the slightly different morphology obtained from this cleaning run, the scientific discussion and conclusions we extrapolate from the 86 GHz data remain unaltered even for a non-univocal solution for imaging.

Fig. C.1. 86 GHz image obtained by means of a second and independent imaging run. Left panel: uniform weighting. Right panel: convolved with the 43 GHz April 2022 beam. The contours start at 3σ.

Fig. C.2. Differences in jet width from the two stacked images at 86 GHz created with the images shown in Fig. 2 and Fig. C.1, respectively. The differences are consistent with zero within 1σ except for the four outermost points in the western jet, which are consistent within 2σ.

Appendix D: Stacked images

In this section, we present the two stacked images in Fig. D.1. The stacked image at 43 GHz (left panel) is created by centering the single-epoch observations on their respective core components and by convolving them with the minimum common circular beam of 0.4 mas. The one at 86 GHz (right panel) is obtained by combining the April 2022 observation presented here with the 2017 one published in Baczko et al. (2024) with a convolved beam of 0.144 mas.

Fig. D.1. Stacked images of NGC 1052. Left panel: stacked image at 43 GHz obtained combining the three observations reported in Fig. 1. Each map is aligned with each respective core as described in Sect. 2.3 and convolved with a circular beam of 0.4 mas. The contours start at 3σ (σ = 0.15 mJy/beam) and increase by a factor of two up to 384σ. Right panel: stacked image at 86 GHz obtained by combining the observation reported in Table 1 and the 86 GHz map shown in Baczko et al. (2022). The circular beam is 0.144 mas and the contours start at 3σ (σ = 0.8 mJy/beam).

Appendix E: Alternative spectral index maps

In this section, we present the collection of the alternative spectral index maps.

At first, we present the (43-86) GHz spectral index map obtained employing the core shift reported in Baczko et al. (2024) of 0.098 mas in declination and 0.176 in right ascension, namely 2 and 4 pixels, respectively. The resultant map is shown in Fig. E.1, in which the color scale highlights the spectral index values, the black contours trace the 86 GHz image, and green contours the 43 GHz one. The two images at the two frequencies appear not to be correctly aligned. On the one hand, the spectral index reaches at the onset of the western jet values up to α ∼ 6.5, higher than the summed contribution from the self-absorbed synchrotron emission (up to α ∼ 2.5) and free–free absorption (up to α ∼ 2.0). On the other hand, at both frequencies we recover a high-intensity region laying at ∼ − 1  mas which appears to be misaligned with the 86 GHz counterpart, which is too far downstream in the 43 GHz jet. Since the core shift is a consequence of the physical properties of the jet, such as the strength of the underlying magnetic fields, variability across different epochs is expected (see, for example, Plavin et al. 2019).

Fig. E.1. Spectral index maps of NGC 1052 between 43 and 86 GHz obtained using the shift from Baczko et al. (2024).

Second, in this section, we show in Fig. E.2 the alternative (43-86) GHz spectral index maps obtained employing the map shown in the main text (Fig. 2) but with the different possible alignments discussed in Sect. 3.3.

Finally, we show the alternative (43-86) GHz spectral index image (Fig. E.3) obtained using Fig. C.1. The position angle in both jets at 86 GHz appears to be slightly off, as indicated by the transverse gradient in the spectral index around 0.5 mas downstream of the core on both jets (we highlight how such a gradient persists with different core shifts applied). Nonetheless, the core in the spectral index map remains highly inverted (α_thick > 2.5), consistent with Fig. 7.

Fig. E.2. (43-86) GHz spectral index images with different shifts applied. Starting from the top-left panel the shifts applied are: (x, y) = (2, 1), (x, y) = (2, 0), (x, y) = (2, −1), (x, y) = (1, 1), (x, y) = (1, −1), (x, y) = (0, 1), (x, y) = (0, 0), (x, y) = (0, −1).

Fig. E.3. Alternative (43-86) GHz spectral index map using the map shown in Fig. C.1. The black contours are for the 86 GHz emission, while the green ones for the 43 GHz image.

All Tables

All Figures

	Fig. 1. 43 GHz VLBA observations of NGC 1052 performed in i) April 2021 (top panel), ii) October 2021 (middle panel), and iii) April 2022 (lower panel). The contours start at 3σ, and the value is reported in Table 1 for each map. The gray ellipse represents the naturally weighted beam.
In the text

	Fig. 2. April 2022 86 GHz GMVA image. Left panel: Uniform weighting. Right panel: Image convolved with the 43 GHz April 2022 beam (0.677 mas × 0.235 mas, −11.8°). The contours start at 3σ.
In the text

	Fig. 3. Core shift of the eastern jet as a function of the frequency. The orange data points are taken from Baczko et al. (2022), and the 43 GHz blue data point is taken from Baczko et al. (2024). The continuous orange line represents the best fit, and the dashed lines represent the corresponding upper and lower limits. The gray area shows the possible best-fit line within 1σ.
In the text

	Fig. 4. Flux density profile for the 43 GHz and 86 GHz stacked images for the eastern and western jet. The double-sided jets are highly symmetrical at 86 GHz, and a dimmer western jet on submilliarcsecond scales is visible at 43 GHz.
In the text

Fig. 5. Deconvolved jet width as a function of the distance from the core for the eastern and western jet. This plot, originally published in Baczko et al. (2024), is updated with the new profiles at 43 GHz and 86 GHz we extrapolated. We show three different fits: i) black continuous line, best-fit line for the high frequency data (43, 86, and 230 GHz); ii) red dashed line, best-fit line for all the data; iii) brown dot-dashed line, best-fit line reported in Baczko et al. (2022). In the tables, W0 are the initial jet widths, ku and kd are the upstream and downstream power law indexes, respectively, and zb are the transition distances.

In the text

	Fig. 6. Distance from the core across the three epochs between April 2021 and April 2022 for the Gaussian components reported in Appendix B. Left panel: Eastern jet. Right panel: Western jet. The differently colored continuous lines show the best-fit curves for the different components.
In the text

	Fig. 7. 43–86 GHz spectral index map of NGC 1052 obtained using the core shift proposed here. The black contours show the 86 GHz emission, and the green contours show the 43 GHz image.
In the text

Fig. 8. Average spatial magnetic field evolution for the two jets in NGC 1052. The horizontal black line represents the magnetic field strength estimates Baczko et al. (2016) at 1 RS, and the vertical green, blue, and orange lines show the values (with their uncertainties) from Equations (5), (6), and (7), respectively. The dashed lines show the extrapolation assuming b = 0.47 and b = 0.94 (B ∝ z−b) for the parabolic region and b = 1.0 and b = 2.0 for the conical region. The red triangle is a lower limit estimate for the magnetic field from Baczko et al. (2016), and the red dot is taken from Baczko et al. (2024). The red box represents the magnetic field that is expected to be needed to saturate the accretion disk.

In the text

	Fig. B.1. Modelfit components (in red) overlaid to the images presented in this paper. Upper panels: 43 GHz April 2021 (left), 43 GHz October 2021 (right). Lower panels: 43 GHz April 2022 (left), 86 GHz April 2022 (right).
In the text

	Fig. C.1. 86 GHz image obtained by means of a second and independent imaging run. Left panel: uniform weighting. Right panel: convolved with the 43 GHz April 2022 beam. The contours start at 3σ.
In the text

	Fig. C.2. Differences in jet width from the two stacked images at 86 GHz created with the images shown in Fig. 2 and Fig. C.1, respectively. The differences are consistent with zero within 1σ except for the four outermost points in the western jet, which are consistent within 2σ.
In the text

Fig. D.1. Stacked images of NGC 1052. Left panel: stacked image at 43 GHz obtained combining the three observations reported in Fig. 1. Each map is aligned with each respective core as described in Sect. 2.3 and convolved with a circular beam of 0.4 mas. The contours start at 3σ (σ = 0.15 mJy/beam) and increase by a factor of two up to 384σ. Right panel: stacked image at 86 GHz obtained by combining the observation reported in Table 1 and the 86 GHz map shown in Baczko et al. (2022). The circular beam is 0.144 mas and the contours start at 3σ (σ = 0.8 mJy/beam).

In the text

	Fig. E.1. Spectral index maps of NGC 1052 between 43 and 86 GHz obtained using the shift from Baczko et al. (2024).
In the text

	Fig. E.2. (43-86) GHz spectral index images with different shifts applied. Starting from the top-left panel the shifts applied are: (x, y) = (2, 1), (x, y) = (2, 0), (x, y) = (2, −1), (x, y) = (1, 1), (x, y) = (1, −1), (x, y) = (0, 1), (x, y) = (0, 0), (x, y) = (0, −1).
In the text

	Fig. E.3. Alternative (43-86) GHz spectral index map using the map shown in Fig. C.1. The black contours are for the 86 GHz emission, while the green ones for the 43 GHz image.
In the text