© ESO 2019

1. Introduction

The giant elliptical galaxy NGC 1052 (z  =  0.005037) harbours a low-luminosity active galactic nucleus (AGN) in its centre. The strong optical emission lines in its spectrum classifies NGC 1052 as a LINER (see e.g., Gabel et al. 2000). The radio structure at arcsecond scales reveal a two-sided structure, but the central core dominates the emission with about 85% of the flux (Wrobel 1984). The two radio lobes have indications of hot spots, with an east-west orientation covering about 3 kpc (projected). The central radio source has a relatively flat spectrum at cm wavelengths with typical flux densities of 1−2 Jy. Very long baseline interferometry (VLBI) images show a milliarcsecond jet and counterjet extending over 15 mas at cm wavelengths (see e.g., Kameno et al. 2003; Kadler et al. 2004a). The jets are propagating in an east-west direction and are oriented at ∼60° in the sky (measured from north to the east). Multi-frequency VLBI studies of NGC 1052 display an emission gap between the eastern and the western jet, whereas the extent of this emission gap is decreasing with frequency and disappears for ν ≥ 43 GHz (Kadler et al. 2004a). These observations are interpreted to be caused by a torus-like structure, so that not only synchrotron self-absorption is present, but also free-free absorptionThe extent of the torus can be studied via high-resolution multi-frequency VLBI observations, providing estimates between 0.1 pc and 0.7 pc (Kameno et al. 2003). Remarkably, the source also displays as well water maser emission (see e.g., Braatz et al. 2003) observed at positions aligned with the radio jet (Claussen et al. 1998).

The source has shown prominent emission in x-rays which provide an estimate for the column density of 10²² cm⁻²–10²⁴ cm⁻² (Kadler et al. 2004b, a). Kadler et al. (2002, 2004b) find indications of extended x-ray emission in addition to the nucleus in a Chandra image, in agreement with radio emission and optical emission from the Hubble Space Telescope (Pogge et al. 2000). Using VLBI monitoring of the source at 15 GHz, 22 GHz and 43 GHz the dynamics and kinematics at parsec-scales can be studied. The typical apparent speeds observed in NGC 1052 are sub-luminal, with ranges between 0.25 c and 0.5 c. Böck (2013) reports a detailed analysis of the MOJAVE observed images between 1995 an 2012 with average sky motions of 0.74  ±  0.06 mas yr⁻¹ which corresponds to projected speeds of 0.230  ±  0.011 c, consistent with the values reported by Vermeulen et al. (2003; 0.73  ±  0.06 mas yr⁻¹), with no significant differences between the values measured in jet and counter jet.

Based on VLBA observations at 43 GHz between 2004 and 2009 Baczko et al. (2019) derived mean speeds of 0.343  ±  0.037 c in the western jet and 0.561  ±  0.034 c in the eastern jet. Using the obtained kinematics the viewing angle, ϑ, can be computed. Vermeulen et al. (2003) report values for the viewing angle of ϑ ≥ 57°, whereas the values of Böck (2013) constrain it to ϑ ≥ 78.8° for a v = 0.238 c (see their Sect. 3.1.3). From the speeds measured in the jet and counter-jet, the upper limit of ϑ ≤ 85° is consistent with the lower limit of 78.8°. Due to significant differences between the eastern and western jet in flux density ratio as well as speed, the viewing angle based on the 43 GHz VLBA observations could only be derived to lie between 60 ≤ ϑ ≤ 90 which is consistent with the aforementioned upper limit.

The study of the 43 GHz structure of the jets shows a significant difference between the eastern and western jet, and so the question arises as to whether the jets in NGC 1052 appear asymmetric due to the presence of an obscuring torus or if the jets are intrinsically asymmetric that is, asymmetries in the jet launching and formation process close to the central black hole. Given a viewing angle of nearly 90°, NGC 1052 is a perfect laboratory to investigate the influence of the surrounding medium including the obscuring torus, the radiation micro-physics and the jet launching mechanism.

To investigate the interplay between the non-thermal emission produced in the jets and the thermal absorption provided by the obscuring torus we performed 2D-axisymmetric special-relativistic hydrodynamic (SRHD) simulations of jets in a decreasing-pressure ambient medium and compute their radiative signatures. Depending on the parameters of the torus that is torus density, temperature and dimensions, flux density asymmetries in the jets can be produced. In addition, spectral indices, α >  2.5 are obtained within the region which is covered by the obscuring torus. The imprint of the torus can also be found in the broadband radio spectrum, either as a flat or as a double-humped spectrum (Fromm et al. 2018). In this paper we combine our jet-torus model with evolutionary algorithms (EA) and address the question of which kind of jet and torus configuration is required to best model the observed radio images and broadband spectrum of NGC 1052.

The organisation of the paper is as follows. In Sect. 2 we introduce the radio galaxy NGC 1052 and present stacked radio images and the average broadband spectrum of the source. The numerical setup for the jet and the torus are introduced in Sect. 3.1 and a summary on the emission simulation is provided in Sect. 3.2. The mathematical and numerical methods used during the optimisation process are explained in Sect. 4. The obtained results and their corresponding discussion can be found in Sect. 5 and in Sect. 6. Throughout the paper we assume an ideal-fluid equation of state , where p is the pressure, ρ the rest-mass density, ϵ the specific internal energy, and the adiabatic index (see e.g., Rezzolla & Zanotti 2013).

The flux density errors, σS_ν, on the synthetic images are computed using a frequency-dependent calibration error σ_cal  ∼  0.1−0.14 and the off-source rms: σS_ν  =  σ_calS_ν + rms, where S_ν is the flux density. We define the spectral index, α, computed between two frequencies, ν₁ and ν₂ as: α_ν1, ν2  =  log₁₀(S_ν1/S_ν2)/log₁₀ (ν₁/ν₂). At its distance of 17.72  ±  1.26 Mpc, an angular separation of 1 mas corresponds to a projected length of about 9.5  ×  10⁻² pc, and a proper motion of 1 mas yr⁻¹ corresponds to an apparent speed of 0.31 c.

2. Observations of NGC 1052

NGC 1052 is a frequently observed source within the MOJAVE survey¹ (Lister et al. 2009) and the FGamma program² (Fuhrmann et al. 2016). These surveys provide an excellent data base for both high-resolution radio images and densely sampled multi-frequency single dish observations. NGC 1052 frequently ejects new components into the jets which are propagating outwards, and during the passage along the jet their flux density fades away. Modelling the variation on top of the underlying steady-state flow requires the injection of perturbation and increases the computational costs. This complication can be circumvented by producing stacked radio images which smear out the variation on top of the underlying flow and provide an average image of the source. Given the wealth of available information on NGC 1052, we produced stacked radio images at 22 GHz and 43 GHz together with the mean broadband radio spectrum³. In Fig. 1 we show the stacked radio images at 22 GHz and 43 GHz (see Table 1 for image parameters) as well as the average broadband radio spectrum. In the following we provide a short description of the data used in this work.

Fig. 1. Stacked VLBA images of NGC 1052 for 22 GHz (top) and 43 GHz (middle). The average radio spectrum for NGC 1052 between 2.6 GHz and 142 GHz provided by FGamma is plotted in the bottom panel. For details on the image stacking see text.

NGC 1052 was observed 29 times between 2004 and 2009 with the VLBA at 22 GHz and 43 GHz. Details on data calibration and the imaging process can be found in Baczko et al. (2019). The typical VLBA beam sizes are about (0.5  ×  0.2) mas at 43 GHz and (0.8  ×  0.3) mas at 22 GHz, the total flux densities during this interval are between 0.5 and 2.0 Jy at both frequencies. As described in Baczko et al. (2019) several properties as for example dynamics and distribution of Gaussian model fit parameters, including brightness temperatures and component sizes, at 43 GHz have been analysed. Based on these properties the jets appear to be asymmetric. To produce stacked images the individual maps at 22 GHz and 43 GHz had been aligned on an optically thin feature at around −2 mas distance to the west of the 22 GHz map peak. The maps were restored with a common beam for all observations and stacked images for each frequency were produced by performing a pixel-by-pixel average over all maps, similar to the procedures described in Pushkarev et al. (2017).

3. Numerical simulations

3.1. SRHD and torus simulations

For our modelling of the jet in NGC 1052 we use the SRHD simulations presented in Fromm et al. (2018) and included additional values for the pressure mismatch, namely d_k = 1.5,  3.5,  4.5 and 5.0. For clarity we provide a short summary here. The simulations are performed using 2D axisymmetric SRHD and the jets are injected in a numerical grid which consists of 320 cells in the radial direction and 400 cells in the axial direction. Using four cells per jet radius (R_j) the numerical grid covers a range of 80 R_j  ×  100 R_j. We inject jets with a velocity of v_j = 0.5 c in the z-direction at the left edge of the grid (z = 0) and we assume a decreasing-pressure ambient medium. The ambient medium follows a King-Profile which is characterised by the core-radius, z_c, and the exponents n and m (see Eq. (1)):

(1)

For the simulations presented in Fromm et al. (2018) we applied z_c = 10 R_j, n = 1.5 and m = 2, 3, 4. Depending on the ratio between the pressure in the jet and the ambient medium, d_k, and the ambient pressure profile different jet morphologies are established (see Fig. 2). If the jet is in pressure balance at the nozzle, nearly featureless conical jets are established, so called pressure-matched (PM) jets. On the other hand, a mismatch between the jet and ambient medium pressure at the nozzle leads to the formation of a series of recollimation shocks and a pinching of the jet boundary is obtained. These jets are referred to as over-pressured (OP) jets (see e.g., Gómez et al. 1997; Mimica et al. 2009; Aloy & Rezzolla 2006; Fromm et al. 2016).

Fig. 2. Stationary results for the jet simulations. The panels show the 2D distribution of the rest-mass density for different ambient medium configurations as indicated by the exponent m and dk (dk = 1 corresponds to a PM jet and dk >  1 to an OP jet). The upper row spans the entire simulation grid whereas the bottom row shows a magnified view of the nozzle region (−7 <  z/Rj <  7). The white lines in the bottom row show stream lines visualising the direction of the flow and the bold dashed lines correspond to the inward travelling and reflected shock wave.

After roughly 5 longitudinal grid crossing times a steady state is recovered and we insert a steady torus into our simulations. This torus acts only as an absorber of the non-thermal emission and is not dynamical evolved with the RHD simulations. This torus is modelled via several parameters which describe its geometry and the distribution of the torus density and temperature (see Fromm et al. 2018, for details). In Table 2 we present an overview of the jet and torus parameters.

3.2. Emission calculations

In order to compare our SRHD models to the observations that is, VLBI images and single dish spectra, we need to compute the non-thermal and thermal emission. For the computation of the emission we follow the recipe presented in Fromm et al. (2018). For completeness we introduce the basic steps for the emission simulation below. Since we are not evolving the non-thermal particles during the SRHD simulations we have to reconstruct their distribution from the SRHD variables that is, pressure, p and density ρ. We assume a power law distribution of relativistic electrons:

(2)

where n₀ is a normalisation coefficient, γ_e is the electron Lorentz factor, γ_e,  min and γ_e,  max are the lower and upper electron Lorentz factors and s is the spectral slope. In the next step we relate the number density of relativistic particles to the number density of thermal particles in the jet via the scaling parameter ζ_e as:

(3)

with m_p the mass of the proton. The energy of the non-thermal particles is related to the energy of the thermal particles via the parameter ϵ_e:

(4)

By performing the integrals in Eqs. (3) and (4) an expression for the lower electron Lorentz factor, γ_e, min can be derived:

(5)

The last step in the construction of the non-thermal particle distributions assumes that the upper electron Lorentz factor, γ_e, max, is a constant fraction ϵ_γ of the lower electron Lorentz factor

(6)

Given the expressions for the electron Lorentz factors (Eq. (5) and (6)) the normalisation coefficient of the particle distribution, n₀, can be written as:

(7)

Given that information on the magnetic field cannot be obtained from our hydrodynamical numerical simulations, we assume its magnitude is a fraction ϵ_B of the equipartition magnetic field

(8)

Finally, we compute the frequency-dependent total intensity, I_ν, for each ray by solving the radiative transport equation

(9)

where ϵ_nt and α_nt are the emission and absorption coefficients for non-thermal emission and α_th is the absorption coefficient for thermal emission (for details on the computation of the emission and absorption coefficients see Fromm et al. 2018).

3.3. Modification of the emission calculations

In our modelling we will compare VLBI images and single dish spectra at various frequencies. Assuming that the emission at higher frequencies is generated mainly near the apex of the jet (due to the higher density and magnetic field within this region as compared to the more extended outer regions), we have to ensure that we sufficiently cover this region within our numerical grid during the ray-tracing. On the other hand, the low frequency emission zones are not restricted to a region close to the apex of the jet. Therefore a minimum numerical resolution has to be provided to resolve the larger extent of the jet. In order to fulfil our above stated criteria on the ray-grid we would need a high numerical resolution which leads to computationally challenging simulations. A way to overcome these computational limitations is to introduce adaptive non-linear grids. The advantage of these grids, as compared to uniform Cartesian grids, is that we can increase the numerical resolution at the apex of the jet and at the same time provide enough resolution at the larger scales. In this way we can keep the computational effort to a minimum (see Fig. 3).

Fig. 3. Rest-mass density distribution for different grid resolutions (left: uniform cartesian 3003, centre: uniform cartesian 5003 and right: adaptive logarithmic 3003). Top row: entire simulation box and bottom row: zoom into the central region indicated by the dashed boxes in the top row. The adaptive grid is created using Δ(x, y, z)align,  min = 10−3, (x, y, z)scale = 4 and ηx, y, z = 1.15.

In order to focus numerical resolution on the jet we first have to align the initial Cartesian coordinate system, indicated by the sub-script “cart”, with the jet direction. The aligned coordinate system labelled with the sub-script “align” is obtained from the Cartesian system via two rotations using the angles ϑ (viewing angle) and φ (orientation of the jet in plane of the sky). After the rotation into the aligned system the coordinates are modified according to:

(10)

(11)

(12)

where Δ(x, y, z)_{align,  min} is the smallest cell spacing in the different directions, (x, y, z)_scale sets the extent of the linear scale that is, number of cells covered with the highest numerical resolution and η_x, y, z scales the exponential growth of the grid. Once this ray grid is established we interpolate the SRHD parameters to each cell using a Delaunay triangulation. The parameters for the numerical grid used in this work are Δ(x, y, z)_{align,  min} = 10⁻³, (x, y, z)_scale = 4 and η_x, y, z = 1.15. A ray propagating through this non-linear grid will encounter cells with different cell sizes and we therefore have to take special care in the computation of the optical depth along ray paths. Consequently, we interpolate the SRHD values required for the computation of the emission along the ray and include a bisection method for the calculation of the optical depth with an accuracy of Δτ ≃ 10⁻⁶. This method guarantees that we are tracing the optical depth cut-off with high precision and leads to converged spectra and images (see Appendix A for a detailed study on the spectra and image convergence).

3.4. Synthetic imaging

A typical VLBI experiment consists of series of on-source scans on the main target and off-source scans for calibration and focussing on a calibrator source. In addition to the reduced on-source time due to calibration of the experiment, the limited number of telescopes participating in the observations lead to a sparse sampling of the source brightness distribution in Fourier space (hereafter u-v plane). Both effects reduce the number of data points in Fourier space (hereafter termed visibilities). During the standard emission calculations these effects are not considered and the obtained images are typically blurred with the observing beam to mimic real observations. However, if we want to compare our emission simulations directly to VLBI observations we have to take the above mentioned effects into account. The calculation of the synthetic observations can be divided into four main steps:

1. Radiative transfer calculation.

2. Setup of the observing array and the observation schedule (duration of scans, integration time).

3. Fourier transformation of the computed intensity and sampling with the projected baselines of the observing array.

4. Reconstruction of the image.

The solution to the radiative transfer problem is presented in Sect. 3.2 and below we provide some details regarding the remaining steps listed above.

Array setup and observation schedule. Throughout this work we use the Very Long Baseline Array (VLBA) as the observing array. The VLBA consists of ten identical 25 m radio telescopes scattered across the USA (see Fig. 4). The antennas are equipped with several receivers allowing multi-frequency observations between 1.6 GHz and 43 GHz (eight of the ten telescopes can also observe at 86 GHz). The sensitivity of a telescope is usually given in the system equivalent flux density (SEFD) which is computed from the system temperature, T_sys, and the effective area of the telescope, A_eff ⁴. In Table 3 we list the SEFD used for our synthetic imaging. The SEFDs of a baseline (telescope 1 and telescope 2) together with the bandwidth, Δν, and the integration time, t_int, determine the thermal noise in the synthetic images according to:

(13)

Fig. 4. Location of the VLBA antennas across the USA.

We use a bandwidth of 256 MHz and an integration time of ten seconds for all frequencies listed in Table 3. In order to include the calibration gaps in the synthetic radio observations we assume ten minutes on-source and 50 min off-source per observing hour. Together with an integration time of ten seconds this transforms into ∼60 data points per scan per baseline (the exact number depends on rise and set time of the source). Additionally we include the effects of 10% gain calibration errors following Chael et al. (2016).

Fourier transformation and sampling. The computed intensity distribution is Fourier transformed and sampled with the projected baselines for the VLBA array. Given that each telescope has certain elevation constraints together with the coordinates of the source in the sky (RA and DEC) and the date of the observation, the source is not always visible for the entire array. We apply a lower elevation cutoff of 10° an upper limit of 85°. The synthetic data is generated using the EHTim library⁵ which we modified to suit our needs.

Image reconstruction. The simulated visibilities are imported into DIFMAPShepherd (1997) and the image is reconstructed using the CLEAN algorithm Högbom (1974) combined with the MODELFIT algorithm. The final image is stored in FITS format for further analysis for example comparison between synthetic and observed images via normalised cross correlation coefficients and other image metrics. In Fig. 5 we show an example from the synthetic imaging routine. The underlying SRHD simulation and emission model parameters are summarised in Table 4. The location of the source in the sky is 2h41m4.799s in RA and −8d15′20.752″ in Dec. We observe the source from 2017-04-04 0:00 UT to 2017-04-04 24:00 UT using the recipe and array configuration mentioned above.

Fig. 5. Example from the synthetic imaging routine. Top panel: logarithm of the rest-mass density in g cm−3. Middle panel: reconstructed radio image at 43 GHz as seen by the VLBA and bottom panels: u-v plane (left) and the visibility amplitude (right). See text for details on the jet model used.

4. Constrained non-linear optimisation

For applying our model to the observational data we need to find the best values for the different parameters listed in Table 2. This task can be considered as a constrained non-linear optimisation problem:

(14)

where x is an m dimensional vector including the model parameters, for example x=[j_mod, ρ_a, ϵ_b, ϵ_e, ζ_e, s, R_out, θ, ρ(R_in), T_sub, k, l]^T and f(x) is the objective function (minimisation function), g_j(x) are the constraints and x_L, i and x_R, i are lower and upper boundaries for the model parameters.

The minimisation function and the constraints can be constructed in such a way that key properties of the data are used to guide the optimisation process which will speed up the convergence of the algorithm. In the case of NGC 1052 we use the number of local flux density maxima along the jet axis and the existence/non-existence of an emission gap between the western jet and the eastern jet.

For the minimisation function we use the least squares computed from the flux density along the jet axes, and the least squares computed from the observed radio spectrum and simulated one, :

(15)

with q denoting the number of images (frequencies) included in the data set with weighting factors w_i. The above mentioned constraints can be expressed as:

(16)

(17)

(18)

The first constraints (Eq. (16)) set the required minimum cross correlation coefficient between the observed and the synthetic radio images, that is, rejecting or accepting solutions based on the structural agreement between the images. Equation (17) increases the agreement between the number of the flux density peaks, , in the radio maps and Eq. (18) is enforcing a consistence between the images regarding a possible emission gap and its extent, . During the optimisation process we allow for all constraints (Eqs. (16)–(18)) a tolerance of 0.01 that is, a cross correlation coefficient of cc_p(x) = 0.79 is accepted as a solution.

The numerical handling of the constraints depends on the implementation of the optimisation algorithm. A common approach includes the addition of penalty functions to the minimisation function, f(x). More details on the constraint implementation can be found in Deb et al. (2002) and Jansen & Perez (2011).

4.1. Optimisation algorithms

The optimisation problem can be solved by two kind of algorithms: gradient-based and gradient-free algorithms. Given the high dimensionality of our problem, together with the high computational effort (ray-tracing and synthetic imaging) a gradient-based algorithm will most likely get stuck in a local minimum and requires a lot of computational resources for mapping out the gradient with sufficient resolution. Therefore, we apply a gradient-free search algorithm. Among these classes of algorithm we select a genetic algorithm (GA) and a particle swarm optimisation (PSO). In the next Section we provide a short introduction to GA and PSO algorithms and refer to Engelbrecht (2007) for further details.

4.2. Genetic algorithm (GA)

Genetic algorithms are motivated by Darwin’s theory of the survival of the fittest. The main steps of a GA are ranking, crossover (mating), and mutation of the individuals for several generations. An individual can be seen as set of parameters, in our case x = [j_mod,ρ_a,ϵ_b,ϵ_e,ζ_e,s,R_out,θ,ρ(R_in),T_sub,k,l]^T, also referred to as a chromosome, and each entry for example θ is labelled as a gene. During the initial step, N random chromosomes are produced and their fitness is computed (in our case the fitness function is given by Eq. (15)). Based on their fitness the chromosomes are ranked and selected for crossover. During the crossover new chromosomes are created from the parent and the fitness of the offsprings is computed. If the fitness of the offspring is improved compared to the parent, the worst parent with respect to the minimisation function can be replaced (details of crossover and replacing methods depend on the implementation of the GA).

In addition to the crossover process, mutation is used to enforce diversity into the population. Mutation is applied to the offspring of the crossover process at a certain rate, typically ≪1, which guarantees that good offsprings are not overwritten. During the mutation, one or more genes within a chromosome are selected and replaced, where again details on the mutation depend on the implementation of the GA. To summarise, the main parameters of a GA are the number of chromosomes, the number of generations (iterations), the fitness function, the rate of crossover and the rate of mutation. In this work use the Non Sorting Genetic Algorithm II (NSGA2; Deb et al. 2002).

4.3. Particle swarm optimisation (PSO)

A Particle swarm optimisation mimics the behaviour of animal swarms, for example birds searching for food. A swarm particle can be described by its position vector x_i and velocity v_i in the parameter space, and the food can be related to the minimisation function, here Eq. (15). The optimisation is driven by the velocity v_i, which reflects both the individual experience of the particle, commonly referred to cognitive component, and the social experience, that is, exchange of information between neighbouring particles. Thus, the update on the position of each particle can be written as:

(19)

Depending on the variant of the PSO, different recipes for the calculation of the velocity exist (see Chap. 16 in Engelbrecht 2007, for details). In the so-called global best PSO, the velocity update of a particle, v_i(t + 1), is computed from the best position found by the entire swarm at the time t and the best position visited so far by the individual particle (i). The weighting between the best position of the swarm and the best position of the individual particles is set by the social and cognitive weights. The PSO terminates after a maximum number of iterations or after a convergence of the solution, that is, Δf(x) <  ϵ, where f is the minimisation function. In this work we make use of the augmented lagrangian particle swarm optimisation (ALPSO; Jansen & Perez 2011).

4.4. Markov chain Monte Carlo (MCMC) simulation

To further investigate the uncertainties of optimal solutions found by the GA or the PSO we perform MCMC simulations using emcee (Foreman-Mackey et al. 2013). As initial positions for the random walkers we employ a Gaussian distribution with a mean equal to the solution found by the GA or PSO and 50% standard deviation. The large standard deviation allows the random walkers to spread out and sufficiently sample the parameter space around the PSO/GA position. This hybrid approach is advantageous in that we don’t have to use large burn-in times during the MCMC simulations, significantly reducing the computational effort and speeding up the calculation. We typically apply 400 random walkers and perform 10³ iterations, which leads to a total number of 4  ×  10⁵ iterations, similar to the number of iterations used by the GA and the PSO runs.

4.5. Optimisation strategy

In order to model the observations of NGC 1052 we employ the following strategy:

• Use different SRHD models presented in Table 2 and Fig. 2.

• Start the optimisation using GA or PSO, typically between 10⁴ and 10⁵ iterations

• Use the broadband radio spectrum to guide the optimisation.

• Explore the parameter space around the best position provided by GA or PSO via MCMC simulations.

The computational costs for a single run depends on the number of frequency points included in the broadband radio spectrum and on the required step size to achieve the optical depth accuracy, Δτ. Typically 100 s are required to compute a radio spectrum consisting of 10 frequencies. This leads to a total required computational time of 300–3000 cpus h. Using MPI parallelisation the duration of the computation can be reduced to a few tens of hours.

We perform a parameter recovery test to explore the capabilities of our end-to-end pipeline by inserting the reference model into our code. The obtained uncertainties on average smaller 30% and the large uncertainties occur mainly due the skewed distribution. Notice that the parameters of the reference model are well recovered within 1σ. Given the complexity of the numerical model and the large degeneracy of the parameters, for example a configuration consisting of a small torus with constant density can lead to a similar radio map and broadband spectrum as larger torus with step density gradients, we consider this deviation as acceptable. The radio images and broadband spectra are most sensitive to the density scaling, ρ_a. This dependency can be understood in the following way: The density scaling is used to convert pressure and density from code units to cgs units via the following relations:

(20)

(21)

Using the recipe presented in Sect. 3.2 pressure and density determine shape of the relativistic particle distribution (see Eq. (7)) and the emission and absorption coefficients (see Eqs. (12)–(20) in Fromm et al. 2018). Therefore ρ_a mainly determines the radio maps and broadband spectra. The additional parameters introduce/modify for example the extent of the gap between the radio jets. In addition to the performed parameter recovery test, we compared our method to THEMIS a MCMC based feature extraction framework (Broderick et al. in prep.) and the results of this comparison can be found in Event Horizon Telescope Collaboration (2019a,b).

5. Results

In Table 5 we present the results of the optimisation runs, for both OP (d_k >  1) and PM (d_k = 1) jets. Both optimisation runs used 4  ×  10⁴ iterations and results of the MCMC can be found in the Appendix A. In Fig. 6 we compare the jet structure and the flux density along the jet axes between the models and observations of NGC 1052. Both jet models OP and PM jet are in good agreement with the observed jet structure (panels a–d in Fig. 6) and a more detailed view of the distribution of the flux density along the jet axis is provided in panels e–h. As mentioned in Sect. 3, OP jets generate local maxima in pressure and density, i.e., recollimation shocks, while PM jets show a monotonic decrease in pressure and density. This behaviour is clearly visible in the flux density profiles along the jet axis (panels e–h). In the 22 GHz flux density profiles (panels e and g) the OP jet shows a plateau-like feature around ±2 mas, whereas the flux density profile of the PM jet is continuously decreasing. With increased resolution, in other words, a higher observing frequency, recollimation shocks closer to the jet nozzle can be resolved and additional local flux density maxima appear (see panel f at ±1 mas) in contrast to the PM jet (see panel h)

Fig. 6. Results of the non-linear optimisation for NGC 1052 for 22 GHz (Cols. 1 and 2) and for 43 GHz (Cols. 3 and 4): Panels a–d: flux density contours for the VLBA observation (black) and the jet models (OP jets in red and PM jet green) and panels e–h: flux density profiles along the blue dashed line in panels a–d. The black points correspond to the VLBA observations of NGC 1052, red to the OP jet model and green to the PM jet. In panels a–d, the peak flux density, Smax, is indicated in the middle and the convolving beam is plotted in lower left corner of each panel. The lowest flux density contour is drawn at 1 mJy and increases by a factor 2.

The broadband radio spectrum between 10⁹ Hz <  ν <  10¹² Hz for the different jet models and NGC 1052 can be seen in Fig. 7. Typically the radio emission seen by single dish telescopes includes radiation emerging from large scale structures not observed by VLBI observations. We take this observational effect into account by including 10% flux density variations indicated by the blue and green shaded bands around the simulated spectrum. Both models, OP and PM jets are able to reproduce the observed spectrum, whereas the OP jet model provides a slightly better fit to the observed spectrum of NGC 1052 (see χ² values in Table 5). Both models exhibit the trend towards lower high-frequency emission (ν >  10¹¹ Hz) and higher low frequency emission (ν <  5  ×  10⁹ Hz). This behaviour can be attributed to the idealised treatment of the non-thermal particles. In our simulations we neglect radiative losses and re-acceleration acting on the non-thermal particles. In Sect. 6 we provide a detailed discussion of the impact of the above mentioned physical mechanism on the broadband spectrum.

Fig. 7. Broadband radio spectrum of NGC 1052 including the averaged observations and the simulated OP and PM jet models. The red and blue shaded regions correspond to 10% flux density variations (see text for details).

A glimpse into the acting radiation microphysics in relativistic jets can be obtained via spectral index studies and the computed spectral index between two frequencies ν₁ and ν₂ is related to the energy distribution of the non-thermal particles, s, via s = −(2α − 1). Changes in the spectral index can be attributed to losses including adiabatic (expansion) and radiative (synchrotron and inverse Compton) (see e.g., Mimica et al. 2009; Fromm et al. 2016) and to re-acceleration of non-thermal particles, for example internal shocks, shear and magnetic reconnection (see e.g., Sironi et al. 2015; Liu et al. 2017). For the calculation of a spectral index in each pixel of the radio map we convolve both radio images with a common beam, typically the one of the smaller frequency (here 22 GHz) and align the structure by common optically thin features in the jet (see Fromm et al. 2013, for details). The computed spectral indices between 22 GHz and 43 GHz are presented in Fig. 8.

Fig. 8. Spectral index maps for NGC 1052. Panel a: 2D-distribution of the spectral index for the OP jet model, panel b: VLBA observations of NGC 1052 and panel c: PM jet model. The variation of the spectral index along the jet axis (black dashed line in panels a–c) can be seen in panel d. In panels a–c, the convolving beam is plotted in the lower left corner of each panel. The contours correspond to the 22 GHz flux density distribution and the lowest flux density contour is drawn at 1 mJy and increases by factors of two. The dash-dotted line in panel d corresponds to the optically thin spectral index of the models (α = 1.4).

The spectral index computed from the OP and PM model approximates at large distances the theoretical value of α = 1.4 (OP jet) and α = 1.45 (PM jet). The small variation in the spectral index for the models is an artefact of the image reconstruction algorithm and the convolution with observing beam (see Appendix A for details). The largest difference between our numerical models and the VLBA observations occurs at a distance of x =   ±  2 mas. These locations coincides in the case of the OP jet with the position of a recollimation shock. As shown in Fromm et al. (2016) the spectral index at the location of recollimation shocks can be increased or inverted if radiative losses are taken into account. However due to computational limitations we excluded radiative cooling in our emission calculations⁶. The central region of NGC 1052 is dominated by absorption due the obscuring torus, which is reflected by large spectral index values exceeding the typical value for synchrotron self-absorption of α >  2.5. With increasing distance from the centre, the spectral index is decreasing (less absorption due the obscuring torus) and approaches the already mentioned value of α = 1.4 (OP jet) and α = 1.45 (PM jet).

In Fig. 9 we present the temperature and density distribution within the torus. The top panels (a and b) show the torus for the OP jet model and the bottom panels (c and d) for the PM jet model. The inner and outer radii for both models are very similar R_in = 0.2  ×  10¹⁸ cm and R_out = 1.5  ×  10¹⁸ cm (0.5 pc) which is within the estimates 0.1 pc ⩽ R_torus ⩽ 0.7 pc reported by Kameno et al. (2003). The torus in the OP jet model has a larger opening angle and less steep temperature gradient than the torus of the PM jet (see Table 5).

Fig. 9. Distribution of the temperature (panels a and c) and the density (panels b and d) in the torus for the OP jet (top) and the PM jet (bottom).

A measurable quantity of the obscuring torus is the number density, which can be computed by integrating the density profile along a line sight. In our case the density of the torus is modelled via:

(22)

where the exponents k_ρ and l_ρ model the decrease in density in radial and Θ directions. Using the values obtained from the non-linear optimisation (see Table 5) we calculate a number density of:

(23)

(24)

Both values are in agreement with number density in NGC 1052 derived from X-ray observations of 0.6  ×  10²² cm⁻² ⩽ N_H, obs ⩽ 0.8  ×  10²² cm⁻² (Kadler et al. 2004b).

6. Discussion

6.1. Multi-frequency VLBA observations and core-shifts

The inclusion of the broadband radio spectrum (see Fig. 7) in our non-linear optimisation process enables us to compute various synthetic images in addition to 22 GHz and 43 GHz. Therefore, we can model the multi-frequency behaviour of the source and compare the obtained results to the observations of NGC 1052. The most remarking feature in the VLBA observations of NGC 1052 is the emission gap between the eastern (left) and western (right) jet, which is shrinking with increasing frequency (see Fig. 1 in Kadler et al. 2004a). Thus, a valid model of NGC 1052 must reproduce this behaviour. In Fig. 10 we present our synthetic multi-frequency images for both jet models, OP jet and PM jet. Since the absolute position of the jets is lost during the image reconstruction, we align the jets by the centre of the emission gap. The emission gap is clearly visible at lower frequencies and the distance between the jets is decreasing with increasing frequency. As mentioned in Sect. 5, the gap between the jets is produced by the combined absorption of thermal particles in the torus and the non-thermal particles in the jet, where the thermal particles provide the major contribution to the appearance of the emission gap at lower frequencies.

Fig. 10. Synthetic multi-frequency VLBA images from 5 GHz towards 43 GHz computed for the OP jet model (left) and the PM jet model (right). The images are aligned by the centre of the emission gap and in all images the lowest flux density plotted is 1 mJy. The convolving beam is plotted next to the radio images. The green stars mark the position of the flux density maximum and the white stars correspond to the position where the flux density reaches 20% of the peak flux density for the first time. The dashed lines trace the position of the aforementioned locations through the multi-frequency radio images.

A more qualitative discussion on which particle distribution is dominating the absorption process at different frequencies can be obtained by means of the core-shift. The radio core of a relativistic jet is usually defined as the τ = 1 surface, that is, the onset of the jet. The optical depth τ is computed from the sum of absorption coefficients for the thermal distribution, α_th, and non-thermal particle distribution, α_nt, along the line of sight:

(25)

Since the absorption coefficient depends on the physical conditions in the jet, (e.g. density, temperature and magnetic field) and on frequency, the shift in the core position can be used to probe the radiation micro-physics (Lobanov 1998).

In the analysis of VLBI data, the jet is typically modelled via several gaussian components representing the observed brightness distribution, and the innermost component is selected as the core. However, such a detailed modelling of our synthetic radio images is beyond the scope of this work. Based on detailed modelling of the NGC 1052 observations by Kadler et al. (2004a) we define the observed onset of the jet as the innermost location where the flux density reaches ∼20% of the peak flux density in each jet⁷. In Fig. 10 the green stars correspond to the flux density maximum in the eastern and western jet, and the white stars mark the onset of the jets using the method described above. The frequency-dependent variation of the core-shift with respect to the centre of the torus is presented in Fig. 11. At lower frequencies the distance between the jets decreases as Δr ∝ ν^−0.3 and continuously steepens towards Δr ∝ ν⁻².

Fig. 11. Results of the core-shift analysis for the OP jet (top panel) and PM jet (bottom panel). Solid lines correspond to the eastern (left) jet, dashed lines to the western (right) jet.

This change in the slope is an indication of a change in the radiation process responsible for the absorption. On larger scales, which are probed by lower frequencies, the opacity is dominated by the thermal particle distribution in the torus. Therefore, the absorption coefficient (in the Rayleigh–Jeans regime hν ≪ kT) can be written as:

(26)

Inserting the temperature and density profiles used for the obscuring torus (and omitting the angular dependence for simplicity):

(27)

(28)

the variation of the core position with frequency can now be written as:

(29)

Inserting the values for k and l obtained from the non-linear optimisation in Eq. (29) leads to:

(30)

(31)

With increasing frequency, the absorption due to the thermal particle distribution is decreasing and the non-thermal particles start to dominate the opacity. Following Lobanov (1998) the core position can be written as:

(32)

where α is the spectral index and b and n are the exponents of the evolution of the magnetic field, B ∝ r^−b, and the density of the non-thermal particles, n₀ ∝ r⁻ⁿ, with distance. In the case of equipartition (kinetic energy density equals magnetic energy density), Eq. (32) leads to: r_core, nt ∝ ν⁻¹ using m = 1 and n = 2. If the non-thermal particle distribution dominates the absorption along a line of sight, we can derive estimates for the evolution of the magnetic field and the non-thermal particles using Eqs. (8) and (7). Both equations depend on the pressure: B ∝ p^1/2 and n₀ ∝ p. The evolution of the pressure is given by Eq. (1), which is correct for PM jets and differs only in the position of the recollimation shocks in the case of OP jets. The evolution of the pressure can be divided into two different regimes: p ∝ r⁻² for r >  r_c and p ∝ r^−4/3 for r ⩽ r_c. Inserting the values reported in Table 5 into Eq. (32) leads to the following core-shift behaviour in the case of non-thermal particles dominating the absorption:

(33)

(34)

Using the developed estimates for the variation of the core position with frequency, Eqs. (30)–(34), we can explain the core-shift behaviour presented in Fig. 11. At lower frequencies, here between 5 GHz and 8 GHz, the thermal particle distribution in the torus provides the major contribution to the absorption along the line of sight which leads to a slope of −0.3. This value is in agreement with the derived values of −0.27 and −0.22, respectively. As the frequency increases, the torus becomes more transparent and the absorption due to the non-thermal particles in the jet starts contributing to the opacity. Therefore, we obtain for 8 GHz <  ν <  22 GHz a steepening of the core-shift from ∝ν^−0.3 to ∝ν⁻¹. At frequencies ν >  22 GHz we probe the regions close to the jet nozzle and as mentioned above, there is a change in the pressure gradient. Therefore we expect a strong steepening of the core-shift in this region. This behaviour is clearly visible in Fig. 11 for ν >  22 GHz for the eastern jets.

Due to the geometry of the jet-torus system relative to the observer, the path of a light ray through the obscuring torus is longer for a ray emerging from the western jet than for a ray leaving the eastern jet. Thus, the thermal absorption for a light ray from the western jet is larger than for its eastern counterpart. Therefore, we expect the steeping of the core-shift to be shifted to higher frequencies in the case of the western jet. This effect can be seen in the western jets for both jet models (see Fig. 11). To illustrate the frequency-dependent position of the τ = 1-surface and its impact on the visible regions in the radio images we plot the opacity for four different frequencies on top of the z − x slice of the rest-mass density for the OP jet model (see Fig. 12).

Fig. 12. Distribution of the rest mass density for the OP jet in the z-x plane at y = 0. The dashed lines trace the τ = 1 surfaces for 2 GHz, 5 GHz, 8 GHz, 15 GHz and 22 GHz. The yellow arrow indicates the direction of the ray-tracing and regions of the jet within the dashed area are obscured by the torus and not visible in the corresponding radio images.

6.2. Over-pressured vs. pressure-matched jets

So far both models, OP jet and PM jet, provide very similar results and can successfully reproduce the observations of NGC 1052. It is therefore difficult to distinguish both models based on current observations. A major difference between OP jets and PM jets is the generation of recollimation shocks. However, at low frequencies the torus can mimic a recollimation shock, that is, a local flux density maxima, by a drop in opacity (see e.g. the flux density peak −1 mas in panel h of Fig. 6). A way out of this dilemma can be provided by μas resolution VLBI observations. This high resolution can be achieved in two ways: increasing the observing frequency and/or increasing the baselines. The space-based radio antenna of the RadioAstron satellite operates at 1.6 GHz, 5 GHz and 22 GHz, and extends the projected baseline up to 10 Earth radii (Kardashev et al. 2013). In addition to the space-based antenna, there are two more VLBI experiments providing μas – resolution: The Global Millimetre VLBI Array (GMVA) at 86 GHz (Lee et al. 2008) and the Event Horizon Telescope (EHT) at 230 GHz (Doeleman 2017). At 230 GHz we will observe regions close to the central engine, requiring a general-relativistic treatment of the MHD and the radiative transport, which will be addressed in a future work. Therefore, we focus in this paper on synthetic radio images as observed by RadioAstron at 22 GHz and the GMVA at 86 GHz. In the Appendix A we provide details on the observing schedule, the array configuration and the SEFDs for the two arrays. For the reconstruction of the synthetic radio images we apply a maximum entropy method MEM provided by the EHTim package, since the CLEAN algorithm tends to produce patchy structure for smooth flux density distributions.

The results of the μas resolution imaging for our jet models are presented in Fig. 13. The 22 GHz RadioAstron images show a clear emission gap between the eastern and western jet, similar to the 22 GHz VLBA images (see Fig. 10). The synthetic radio images from both models show a similar flux density distribution and more details can be seen in the flux density along the jet axis (panel c in Fig. 13). The variation of the flux density along the jet axis is comparable in both models and therefore it is not possible to distinguish both models based on RadioAstron images.

Fig. 13. Synthetic radio images for NGC 1052 as seen by the RadioAstron satellite (panels a and b) at 22 GHz. The OP jet is presented in panel a and PM jet in panel b. The convolving beam, 160 μas  ×  20 μas, is plotted in the lower right corner. The flux density along the jet axis, white dashed line in panels a and b, is shown in panel c. See text and Appendix A for details on the observation schedule and the array settings.

At 86 GHz the torus is optically thin and a clear view to the central region is obtained, that is, no absorption due the obscuring torus is seen. The GMVA observations provide, for the first time, a clear detectable difference between the OP jet and PM jet: there is a local flux density maximum at ± 0.5 mas which is not seen in the PM jet (see Fig. 14). Given the most recent GMVA observations of NGC 1052 there are two bright features next to core at roughly ±0.4 mas (Baczko et al. 2016) which could be interpreted as the recollimation shocks.

Fig. 14. Synthetic radio images for NGC 1052 as seen by the GMVA (panels a and b) at 86 GHz. The OP jet is presented in panel a and the PM jet in panel b. The convolving beam, 345 μas  ×  82 μas, is plotted in the lower right corner. The flux density along the jet axis, white dashed line in panels a and b is shown in panel c. See text and Appendix A for details of the observation schedule and the array settings.

Based on this result, together with the variation in the spectral index (see Fig. 8), we favour the over-pressured jet scenario as the most likely configuration of NGC 1052.

6.3. Refining the SRHD model

Based on the discussion in the previous section we conclude that NGC 1052 is best modelled by an over-pressured jet.

To further investigate and refine our obtained solution for NGC 1052 we increase the number of OP jets in our simulation library by generating additional SRHD simulations. Do reduce the computational costs we perform a grid search in the SRHD parameters d_k (pressure-mismatch) and the core radius, r_c (see Eq. (1)). We vary d_k between 1.2 and 2.0 in steps of 0.1 and change r_c within 2 R_j and 10 R_j in steps of 2 R_j. We add these models to our simulation library and re-run the optimisation procedure using the values reported in Table 5 as initial positions for the parameter search.

To avoid biasing effects during the optimisation we add some random scatter on the initial position. This approach has the advantage that computational efforts for the refined study are lower than for an optimisation starting from random positions in the parameter space (as performed in Sect. 5).

The result of the refinement of the underlying SRHD simulation is presented in Fig. 15 and the obtained values are presented in Table 6. The main difference between the initial model and the refined one is the core-radius of the ambient medium, which reduced from z_c = 10 R_j to z_c = 8 R_j, while d_k = 1.5 is not changed. Reducing the core-radius, z_c, induces an earlier transition from p_a(z)∝z^−1.3 to p_a(z)∝z⁻² (see Eq. (1)). Due to the steeper decrease in the ambient pressure, the position where a transversal equilibrium between the pressure in jet and the ambient medium is established is shifted downstream. Therefore, the recollimation shocks will be formed at a larger distance from the jet nozzle as compared to the jet embedded in an ambient medium with a larger core size. This effect can be seen in flux density cuts along the jet axis in panels e–h in Fig. 15. The local flux density maximum at 2 mas is better approximated by the refined model than by the initial one. In addition, the innermost flux density peaks at ∼1 mas are also better fitted by the refined model.

Fig. 15. Same as Fig. 6 including the refined OP jet model. See text for details.

6.4. Limitations of the model

Our current model is able to successfully reproduce several features of the VLBA observations, including the extent of the torus, the number density within the torus, the frequency-dependent emission gap between eastern and western jets, and the distribution of the flux density and spectral index. However, some details of the NGC 1052 observations cannot be reproduced with high accuracy. The observed flux density evolution along the jet axis in the western jet is decreasing faster than in the OP and PM models, while flux density evolution in the eastern jet is very well-reproduced. Since our underlying SRHD jet models are symmetric, this could be an indication of: (i) asymmetries in the ambient medium and/or (ii) asymmetries in the jet launching. These limitations could be addressed by future 3D simulations embedded in a slightly asymmetric ambient medium.

Since we use SRHD simulations to model the jets in NGC 1052, the magnetic field is not evolved as an independent parameter and we compute the magnetic field from the equipartition pressure (see Eq. (8)). Therefore, we restrict ourself to ϵ_b <  0.5 during the modelling and optimisation process. Given a viewing angle ϑ >  80°, no asymmetries are expected even if the dominating component of the field is toroidal. From a dynamical point of view, a strong toroidal magnetic field could reduce the distance between the recollimation shocks for a given jet overpressure (see e.g., Mizuno et al. 2015; Martí et al. 2016), thus allowing for a larger overpressure factor d_k between the jet and the ambient medium in the jets of NGC 1052.

In our current model we ignore the impact of radiative cooling and re-acceleration on the relativistic particles. Depending on the strength of the magnetic field, radiative losses can lead to a steeping of spectral indices and a shortening of the jets at high frequencies (Mimica et al. 2009; Fromm et al. 2016). In our case, with ϵ_b <  0.5 we expect only a small impact of the radiative losses on the large scale structure of the jets. In addition, the magnetic field is decreasing with distance from the jet nozzle which will further reduce its influence on the jet structure. However, at the jet nozzle, at high frequencies (ν >  86 GHz in the case of NGC 1052 Baczko et al. 2016) the magnetic field will become important for both the dynamics and the radiative properties of the jet. Our current model requires a large value for the spectral slope s ≈ 4 to model the structure of NGC 1052 (especially the distribution of the spectral index). Such a large spectral slope can be obtained from a non-thermal particle distribution with s ∼ 2.2, if radiative losses are taken into account. Especially between the jet nozzle and the first recollimation shock, relativistic particles with large γ_e will suffer radiative losses which will steepen the particle distribution and lead to large spectral slopes.

Several processes such as magnetic reconnection, diffuse shock-acceleration and shear acceleration are able to re-energise the relativistic particles during their propagation through the jet (Gonzalez & Parker 2016; Vaidya et al. 2018; Liu et al. 2017). The re-acceleration of relativistic particles leads to an increase (flattening) of the spectral indices. Thus, the steepening of the spectral indices due to radiative cooling can be partially compensated by the various re-acceleration processes.

However, such a self-similar treatment of the non-thermal particles during the optimisation process is currently computationally too demanding.

7. Conclusions and outlook

In this work we present an end-to-end pipeline for the modelling of relativistic jets using state-of-the art SRHD and emission simulations coupled via evolutionary algorithms to high resolution radio images. We use this newly-developed pipeline to model stacked radio images and the broadband radio spectra of NGC 1052. The obtained results, that is, synthetic radio images and broadband spectrum, mimic very well their observed counterparts and the recovered parameters for the obscuring torus are in agreement with derived estimates from radio and X-ray observations. The detailed comparison with available data leads to the conclusion that NGC 1052 is best described by an over-pressured jet (d_k = 1.5) in a decreasing pressure ambient medium (p ∝ r^−1.3).

In a follow up work we will explore the time evolution of the jet in NGC 1052 and the impact of time-delays (slow-light radiative transfer) on the radio structure using the obtained values for the jet and torus as initial parameters. The obtained jet configuration can also be used as a framework for further studies including the improvement of our radiation model including radiative loss and re-acceleration mechanisms of the non-thermal particles along the flow.

To overcome recent limitations with respect to the magnetisation of the jets we will couple in a follow-up work our G(S)RMHD and polarised radiative transfer codes to the presented pipeline. This improvement will allow us to drop the limitations on ϵ_b and furthermore enables us to include polarised observations that is, fraction of polarisation and rotation measures, in the modelling. The inclusion of polarisation will provide an additional independent constraint on the magnetic field strength and its geometry, which will allow us to investigate different magnetic field configurations.

Acknowledgments

Support comes from the ERC Synergy Grant “BlackHoleCam – Imaging the Event Horizon of Black Holes” (Grant 610058). ZY acknowledges support from a Leverhulme Trust Early Career Fellowship. M.P. acknowledges partial support from the Spanish MICINN grant AYA-2013-48226-C03-02-P and the Generalitat Valenciana grant PROMETEOII/2014/069. CMF wants to thank Walter Alef and Helge Rottmann for useful comments and fruitful discussions on synthetic imaging and image reconstruction. This research has made use of data obtained with the Very Long Baseline Array (VLBA). The VLBA is an instrument of the National Radio Astronomy Observatory, a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. This work has made use of NASA’s Astrophysics Data System (ADS).

References

Appendix A: Convergence studies and array configurations

Here we provide details on the convergence studies performed for both the broadband radio spectrum and synthetic images, and on the parameter recovery test. For the convergence study we use the parameters presented in Table 4 and increased the numerical resolution in 100³ steps.

A.1. Image convergence

In order to quantify the convergence of our computed radio image we follow the approach of Lu et al. (2016), Mizuno et al. (2018) and use the structured dissimilarity (DSSIM) index (for details see Wang et al. 2004).

Using the definition of the DSSIM, two identical images would have SSIM = 1 and DSSIM = 0. For this study we select as reference the images obtained with the highest numerical resolution (here 800³ cells). The DSSIM is computed between the reference image and an image with lower numerical resolution for example n = 700³ cells. For each frequency and numerical resolution we obtain a DSSIM value. The result of this study can be seen in Fig. A.1. The calculated DSSIM varies between 0.05 and 0.0001 which indicates an overall good agreement between the images. The DSSIM decreases with increased numerical resolution as expected. The local maximum between 1  ×  10⁹ Hz and 5  ×  10¹⁰ Hz is due to the torus which leads to higher absorption. The red curve indicates the results for the adaptive linear-logarithmic grid. The obtained DSSIM for this grid is in the low frequency regime similar to the uniform cartesian grid with the same number of cells. However, the strength of this grid becomes obvious at higher frequencies where the DSSIM drops below its uniform cartesian counterpart and approaches the DSSIM of the 500³ grid.

A.2. Spectral convergence

In Fig. A.2 we show the single-dish spectrum between (10⁷−10¹³) Hz for different resolutions and the inset panel shows the variation of the total flux density with respect to the number of grid cells. The calculated flux density converges to a common value with a variation of . The increased absorption of the obscuring torus leads to a decrease in the flux density between 3  ×  10⁸ Hz and 5  ×  10¹⁰ Hz. The small variation see in the inset of Fig. A.2 is caused by absorption in obscuring torus at high frequencies ν >  1  ×  10¹⁰ Hz. The absorption in the torus is calculated from the temperature and the density. If we increase the numerical resolution variations both quantities can be better resolved by the grid. This higher resolution can lead to lower or higher absorption along a ray depending on its path through the torus. As a results, there are some small variations in the emission at high frequencies. These variations vanish once the ray-grid exhibits the native resolution of the SRHD simulations, which is 800³.

Fig. A.1. Convergence test for the synthetic images. The synthetic images produced with 8003 resolution are used as reference images.

Fig. A.2. Convergence test for the broadband radio spectrum for eight different resolutions.

A.3. Impact of the image reconstruction and beam convolution

In order to quantify the impact of the image reconstruction algorithm and the beam convolution on the distribution of the flux density and therefore on the spectral index we performed a test on two different reconstruction algorithms. In Fig. A.3 we show reconstructed images for the OP jet model using the CLEAN algorithm as implemented in DIFMAP (panel a) and MEM algorithm from the EHTim package (panel c). The reconstructed images are convolved with a common beam of 0.27  ×  0.72 mas and flux along the jet axis is compared to the blurred infinite resolution image (panel b). Both reconstruction algorithms provide similar results and are capable of reproducing the distribution of the flux density along the jet axis (panel d). In panel e of Fig. A.3 we compute flux density difference between the theoretical value and the reconstructed values. The average flux density difference is around 10–15%. This error could be added to the error budget on the flux density which will improve the χ² values presented in Tables 5 and 6.

Fig. A.3. Result of the image reconstruction algorithm using CLEAN (panel a) and the MEM (panel c) as compared to the convolved ray-traced image (panel b). The flux density along the black dashed line in panels a–c is shown in panel d. The relative deviation from the theoretical flux density in precent is plotted in panel e.

Fig. A.4. Location of the different GMVA radio antennas (for station code see Table A.1).

A.4. Global Millimetre VLBI Array (GMVA)

The GMVA⁸ operates at 86 GHz and currently consists of the telescopes of the VLBA, the European VLBI Network (EVN) and ALMA. Future observations may also include the Korean VLBI stations. The system equivalent flux densities (SEFD) for the involved telescopes are listed in Table A.1 and the location of the telescopes can be seen in Fig. A.4.

For the synthetic GMVA observations of NGC 1052 we use an integration time of 12 s, a on-source scan length of 10 min and a calibration and pointing gap of 50 min. The observing date is set to October 4th 2004 and we observe NGC 1052 from 22:00 UT until 14:00 UT on the next day. The u-v plane together with the amplitude and phase for the OP jet model are displayed in Fig. A.5.

Fig. A.5. Visibilities for the synthetic RadioAstron observations using the parameters listed in Table A.1 and the observing schedule described in the text. The top panel shows the sampling of the u-v plane, middle panel the visibility amplitude and the bottom panel the phase with uv-distance. For reasons of clarity only every 50th data point is plotted.

Fig. A.6. Array configuration for the synthetic RadioAstron observations. The green line corresponds to the ground track of the RadioAstron satellite and black points indicate the individual scans. The ground stations are indicated by the blue points and the station codes are listed in Table A.2.

A.5. Spektr-R and ground array at 22 GHz

The space based satellite Spektr-R extents the baselines up to 10 Earth radii and provides μas resolution. The SEFD of the ground stations and the space antenna are listed in Table A.2. The observing schedule for our synthetic RadioAstron observations of NGC 1052 is the following: 30 min on source-scans and 95 min off-source with an integration time of 10 s. The observations are performed on the 17th of October from 17:00 UT until 18th of October 10:00 UT. This time span corresponds to the periastron transit of the satellite and the ground track of the satellite together with the ground array can be seen in Fig. A.6. The solid green line is the ground track of RadioAstron and the black points on top of the green line corresponds to the on-source scans. Notice, that due the change in the orbital velocity of the satellite during the observations the spacial extent of the black points on top of the green lines shrinks. The location of the ground stations are marked by blue circles (see Table A.2 for station codes). The u-v plane and the visibility amplitude and phase for the OP jet model are plotted in Fig. A.7. Notice the large extent up to 12  ×  10⁹λ in North–South direction which leads to a beam of 160 μas  ×  20 μas.

A.6. Parameter recovering tests and MCMC results

To test and explore the capabilities of your numerical optimisation we perform a parameter recovering test using the reference model given in Table 4 and methods described in Sect. 4. We inject the reference model into our end-to-end pipeline and investigate the ability to recover the injected parameters. For the PSO we use 400 particles, 200 outer iterations and 5 inner iterations⁹ refer and for the MCMC we employ 400 random walkers and 1000 iterations. These number lead to similar number of total iterations, namely 4  ×  10⁵. The results of the parameter recovering test can be seen in Fig. A.8 and in Table A.3. All histograms show clear maxima and injected model values are within 1σ of the recovered values.

The spread around the optimal solution found for the OP and PM jet models are presented in Tables A.4 and A.5 and the histograms for the distribution can be found Figs. A.9 and A.10.

Fig. A.7. Visibilities for the synthetic RadioAstron observations using the parameters listed in Table A.2 and the observing schedule described in the text. Top panel: sampling of the u-v plane, middle panel: visibility amplitude and bottom panel: phase with uv-distance. For reasons of clarity only every 50th data point is plotted.

Fig. A.8. Results for the parameter recovering test for the MCMC simulation and the PSO. The histogram show the distribution of the parameters obtained from the MCMC run, solid black line indicates mean value and the dash-dotted line refer to 1σ standard deviation. The red dashed lines correspond to the parameters of the injected reference model and the green dashed lines to the best position recovered by the PSO.

Fig. A.9. Results for the MCMC simulation for the OP jet model with the PSO best position as starting point and a standard deviation of 50%. The histograms show the distribution of the parameters obtained from the MCMC run. The solid black line indicates mean value and the dash-dotted line refer to 1σ standard deviation computed from the MCMC results. The green dashed lines correspond to the best position recovered by the PSO.

Fig. A.10. Results for the MCMC simulation for the PM jet model with the PSO best position as starting point and a standard deviation of 50%. The histograms show the distribution of the parameters obtained from the MCMC run. The solid black line indicates mean value and the dash-dotted line refer to 1σ standard deviation computed from the MCMC results. The green dashed lines correspond to the best position recovered by the PSO.

All Tables

All Figures

	Fig. 1. Stacked VLBA images of NGC 1052 for 22 GHz (top) and 43 GHz (middle). The average radio spectrum for NGC 1052 between 2.6 GHz and 142 GHz provided by FGamma is plotted in the bottom panel. For details on the image stacking see text.
In the text

Fig. 2. Stationary results for the jet simulations. The panels show the 2D distribution of the rest-mass density for different ambient medium configurations as indicated by the exponent m and dk (dk = 1 corresponds to a PM jet and dk >  1 to an OP jet). The upper row spans the entire simulation grid whereas the bottom row shows a magnified view of the nozzle region (−7 <  z/Rj <  7). The white lines in the bottom row show stream lines visualising the direction of the flow and the bold dashed lines correspond to the inward travelling and reflected shock wave.

In the text

	Fig. 3. Rest-mass density distribution for different grid resolutions (left: uniform cartesian 3003, centre: uniform cartesian 5003 and right: adaptive logarithmic 3003). Top row: entire simulation box and bottom row: zoom into the central region indicated by the dashed boxes in the top row. The adaptive grid is created using Δ(x, y, z)align,  min = 10−3, (x, y, z)scale = 4 and ηx, y, z = 1.15.
In the text

	Fig. 4. Location of the VLBA antennas across the USA.
In the text

	Fig. 5. Example from the synthetic imaging routine. Top panel: logarithm of the rest-mass density in g cm−3. Middle panel: reconstructed radio image at 43 GHz as seen by the VLBA and bottom panels: u-v plane (left) and the visibility amplitude (right). See text for details on the jet model used.
In the text

Fig. 6. Results of the non-linear optimisation for NGC 1052 for 22 GHz (Cols. 1 and 2) and for 43 GHz (Cols. 3 and 4): Panels a–d: flux density contours for the VLBA observation (black) and the jet models (OP jets in red and PM jet green) and panels e–h: flux density profiles along the blue dashed line in panels a–d. The black points correspond to the VLBA observations of NGC 1052, red to the OP jet model and green to the PM jet. In panels a–d, the peak flux density, Smax, is indicated in the middle and the convolving beam is plotted in lower left corner of each panel. The lowest flux density contour is drawn at 1 mJy and increases by a factor 2.

In the text

	Fig. 7. Broadband radio spectrum of NGC 1052 including the averaged observations and the simulated OP and PM jet models. The red and blue shaded regions correspond to 10% flux density variations (see text for details).
In the text

Fig. 8. Spectral index maps for NGC 1052. Panel a: 2D-distribution of the spectral index for the OP jet model, panel b: VLBA observations of NGC 1052 and panel c: PM jet model. The variation of the spectral index along the jet axis (black dashed line in panels a–c) can be seen in panel d. In panels a–c, the convolving beam is plotted in the lower left corner of each panel. The contours correspond to the 22 GHz flux density distribution and the lowest flux density contour is drawn at 1 mJy and increases by factors of two. The dash-dotted line in panel d corresponds to the optically thin spectral index of the models (α = 1.4).

In the text

	Fig. 9. Distribution of the temperature (panels a and c) and the density (panels b and d) in the torus for the OP jet (top) and the PM jet (bottom).
In the text

Fig. 10. Synthetic multi-frequency VLBA images from 5 GHz towards 43 GHz computed for the OP jet model (left) and the PM jet model (right). The images are aligned by the centre of the emission gap and in all images the lowest flux density plotted is 1 mJy. The convolving beam is plotted next to the radio images. The green stars mark the position of the flux density maximum and the white stars correspond to the position where the flux density reaches 20% of the peak flux density for the first time. The dashed lines trace the position of the aforementioned locations through the multi-frequency radio images.

In the text

	Fig. 11. Results of the core-shift analysis for the OP jet (top panel) and PM jet (bottom panel). Solid lines correspond to the eastern (left) jet, dashed lines to the western (right) jet.
In the text

	Fig. 12. Distribution of the rest mass density for the OP jet in the z-x plane at y = 0. The dashed lines trace the τ = 1 surfaces for 2 GHz, 5 GHz, 8 GHz, 15 GHz and 22 GHz. The yellow arrow indicates the direction of the ray-tracing and regions of the jet within the dashed area are obscured by the torus and not visible in the corresponding radio images.
In the text

Fig. 13. Synthetic radio images for NGC 1052 as seen by the RadioAstron satellite (panels a and b) at 22 GHz. The OP jet is presented in panel a and PM jet in panel b. The convolving beam, 160 μas  ×  20 μas, is plotted in the lower right corner. The flux density along the jet axis, white dashed line in panels a and b, is shown in panel c. See text and Appendix A for details on the observation schedule and the array settings.

In the text

Fig. 14. Synthetic radio images for NGC 1052 as seen by the GMVA (panels a and b) at 86 GHz. The OP jet is presented in panel a and the PM jet in panel b. The convolving beam, 345 μas  ×  82 μas, is plotted in the lower right corner. The flux density along the jet axis, white dashed line in panels a and b is shown in panel c. See text and Appendix A for details of the observation schedule and the array settings.

In the text

	Fig. 15. Same as Fig. 6 including the refined OP jet model. See text for details.
In the text

	Fig. A.1. Convergence test for the synthetic images. The synthetic images produced with 8003 resolution are used as reference images.
In the text

	Fig. A.2. Convergence test for the broadband radio spectrum for eight different resolutions.
In the text

	Fig. A.3. Result of the image reconstruction algorithm using CLEAN (panel a) and the MEM (panel c) as compared to the convolved ray-traced image (panel b). The flux density along the black dashed line in panels a–c is shown in panel d. The relative deviation from the theoretical flux density in precent is plotted in panel e.
In the text

	Fig. A.4. Location of the different GMVA radio antennas (for station code see Table A.1).
In the text

	Fig. A.5. Visibilities for the synthetic RadioAstron observations using the parameters listed in Table A.1 and the observing schedule described in the text. The top panel shows the sampling of the u-v plane, middle panel the visibility amplitude and the bottom panel the phase with uv-distance. For reasons of clarity only every 50th data point is plotted.
In the text

	Fig. A.6. Array configuration for the synthetic RadioAstron observations. The green line corresponds to the ground track of the RadioAstron satellite and black points indicate the individual scans. The ground stations are indicated by the blue points and the station codes are listed in Table A.2.
In the text

	Fig. A.7. Visibilities for the synthetic RadioAstron observations using the parameters listed in Table A.2 and the observing schedule described in the text. Top panel: sampling of the u-v plane, middle panel: visibility amplitude and bottom panel: phase with uv-distance. For reasons of clarity only every 50th data point is plotted.
In the text

Fig. A.8. Results for the parameter recovering test for the MCMC simulation and the PSO. The histogram show the distribution of the parameters obtained from the MCMC run, solid black line indicates mean value and the dash-dotted line refer to 1σ standard deviation. The red dashed lines correspond to the parameters of the injected reference model and the green dashed lines to the best position recovered by the PSO.

In the text

Fig. A.9. Results for the MCMC simulation for the OP jet model with the PSO best position as starting point and a standard deviation of 50%. The histograms show the distribution of the parameters obtained from the MCMC run. The solid black line indicates mean value and the dash-dotted line refer to 1σ standard deviation computed from the MCMC results. The green dashed lines correspond to the best position recovered by the PSO.

In the text

Fig. A.10. Results for the MCMC simulation for the PM jet model with the PSO best position as starting point and a standard deviation of 50%. The histograms show the distribution of the parameters obtained from the MCMC run. The solid black line indicates mean value and the dash-dotted line refer to 1σ standard deviation computed from the MCMC results. The green dashed lines correspond to the best position recovered by the PSO.

In the text