Issue 
A&A
Volume 629, September 2019



Article Number  A4  
Number of page(s)  24  
Section  Extragalactic astronomy  
DOI  https://doi.org/10.1051/00046361/201834724  
Published online  22 August 2019 
Using evolutionary algorithms to model relativistic jets
Application to NGC 1052
^{1}
Institut für Theoretische Physik, Goethe Universität, MaxvonLaueStr. 1, 60438 Frankfurt, Germany
email: cfromm@th.physik.unifrankfurt.de
^{2}
MaxPlanckInstitut für Radioastronomie, Auf dem Hügel 69, 53121 Bonn, Germany
^{3}
Mullard Space Science Laboratory, University College London, Holmbury St. Mary, Dorking, Surrey RH5 6NT, UK
^{4}
Anton Pannekoek Institute for Astronomy, University of Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands
^{5}
Departament d’Astronomia i Astrofísica, Universitat de València, Dr. Moliner 50, 46100 Burjassot, València, Spain
^{6}
Observatori Astronòmic, Parc Científic, Universitat de València, C/ Catedràtic José Beltrán 2, 46980 Paterna, València, Spain
^{7}
Frankfurt Institute for Advanced Studies, RuthMoufangStrasse 1, 60438 Frankfurt, Germany
Received:
27
November
2018
Accepted:
18
May
2019
Context. Highresolution very long baseline interferometry (VLBI) observations of NGC 1052 show a two sided jet with several regions of enhanced emission and a clear emission gap between the two jets. This gap shrinks with increasing frequency and vanishes around ν ∼ 43 GHz. The observed structures are due to both the macroscopic fluid dynamics interacting with the surrounding ambient medium including an obscuring torus and the radiation microphysics. In order to model the observations of NGC 1052 via stateofthe art numerical simulations both the fluiddynamical and emission processes have to be taken into account.
Aims. In this paper we investigate the possible physical conditions in relativistic jets of NGC 1052 by directly modelling the observed emission and spectra via stateoftheart specialrelativistic hydrodynamic (SRHD) simulations and radiative transfer calculations.
Methods. We performed SRHD simulations of overpressured and pressurematched jets using the specialrelativistic hydrodynamics code Ratpenat. To investigate the physical conditions in the relativistic jet we coupled our radiative transfer code to evolutionary algorithms and performed simultaneous modelling of the observed jet structure and the broadband radio spectrum. During the calculation of the radiation we consider nonthermal emission from the jet and thermal absorption in the obscuring torus. In order to compare our model to VLBI observations we take into account the sparse sampling of the uv plane, the array properties and the imaging algorithm.
Results. We present for the first time an endtoend pipeline for fitting numerical simulations to VLBI observations of relativistic jets taking into account the macrophysics including fluid dynamics and ambient medium configurations together with thermal and nonthermal emission and the properties of the observing array. The detailed analysis of our simulations shows that the structure and properties of the observed relativistic jets in NGC 1052 can be reconstructed by a slightly overpressured jet (d_{k} ∼ 1.5) embedded in a decreasing pressure ambient medium
Key words: galaxies: active / galaxies: jets / radio continuum: galaxies / radiation mechanisms: nonthermal / radiative transfer / hydrodynamics
© ESO 2019
1. Introduction
The giant elliptical galaxy NGC 1052 (z = 0.005037) harbours a lowluminosity 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 twosided 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 eastwest 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 eastwest direction and are oriented at ∼60° in the sky (measured from north to the east). Multifrequency 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 toruslike structure, so that not only synchrotron selfabsorption is present, but also freefree absorptionThe extent of the torus can be studied via highresolution multifrequency 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 xrays which provide an estimate for the column density of 10^{22} cm^{−2}–10^{24} cm^{−2} (Kadler et al. 2004b, a). Kadler et al. (2002, 2004b) find indications of extended xray 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 parsecscales can be studied. The typical apparent speeds observed in NGC 1052 are subluminal, 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^{−1} 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^{−1}), 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 counterjet, 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 microphysics and the jet launching mechanism.
To investigate the interplay between the nonthermal emission produced in the jets and the thermal absorption provided by the obscuring torus we performed 2Daxisymmetric specialrelativistic hydrodynamic (SRHD) simulations of jets in a decreasingpressure 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 doublehumped spectrum (Fromm et al. 2018). In this paper we combine our jettorus 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 idealfluid equation of state , where p is the pressure, ρ the restmass 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 frequencydependent calibration error σ_{cal} ∼ 0.1−0.14 and the offsource rms: σS_{ν} = σ_{cal}S_{ν} + rms, where S_{ν} is the flux density. We define the spectral index, α, computed between two frequencies, ν_{1} and ν_{2} as: α_{ν1, ν2} = log_{10}(S_{ν1}/S_{ν2})/log_{10} (ν_{1}/ν_{2}). 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^{−2} pc, and a proper motion of 1 mas yr^{−1} 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^{1} (Lister et al. 2009) and the FGamma program^{2} (Fuhrmann et al. 2016). These surveys provide an excellent data base for both highresolution radio images and densely sampled multifrequency 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 steadystate 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^{3}. 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. 
Image parameters for the stacked VLBA images.
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 pixelbypixel 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 zdirection at the left edge of the grid (z = 0) and we assume a decreasingpressure ambient medium. The ambient medium follows a KingProfile which is characterised by the coreradius, z_{c}, and the exponents n and m (see Eq. (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 pressurematched (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 overpressured (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 restmass density for different ambient medium configurations as indicated by the exponent m and d_{k} (d_{k} = 1 corresponds to a PM jet and d_{k} > 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/R_{j} < 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 nonthermal 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.
Parameters for the emission simulations.
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 nonthermal 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 nonthermal 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:
where n_{0} 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:
with m_{p} the mass of the proton. The energy of the nonthermal particles is related to the energy of the thermal particles via the parameter ϵ_{e}:
By performing the integrals in Eqs. (3) and (4) an expression for the lower electron Lorentz factor, γ_{e, min} can be derived:
The last step in the construction of the nonthermal particle distributions assumes that the upper electron Lorentz factor, γ_{e, max}, is a constant fraction ϵ_{γ} of the lower electron Lorentz factor
Given the expressions for the electron Lorentz factors (Eq. (5) and (6)) the normalisation coefficient of the particle distribution, n_{0}, can be written as:
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
Finally, we compute the frequencydependent total intensity, I_{ν}, for each ray by solving the radiative transport equation
where ϵ_{nt} and α_{nt} are the emission and absorption coefficients for nonthermal 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 raytracing. 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 raygrid we would need a high numerical resolution which leads to computationally challenging simulations. A way to overcome these computational limitations is to introduce adaptive nonlinear 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. Restmass density distribution for different grid resolutions (left: uniform cartesian 300^{3}, centre: uniform cartesian 500^{3} and right: adaptive logarithmic 300^{3}). 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 subscript “cart”, with the jet direction. The aligned coordinate system labelled with the subscript “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:
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^{−3}, (x, y, z)_{scale} = 4 and η_{x, y, z} = 1.15. A ray propagating through this nonlinear 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^{−6}. This method guarantees that we are tracing the optical depth cutoff 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 onsource scans on the main target and offsource scans for calibration and focussing on a calibrator source. In addition to the reduced onsource 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 uv 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:

Radiative transfer calculation.

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

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

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 multifrequency 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} ^{4}. 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:
Fig. 4. Location of the VLBA antennas across the USA. 
Used SEFD for the VLBA.
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 onsource and 50 min offsource 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^{5} 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 20170404 0:00 UT to 20170404 24:00 UT using the recipe and array configuration mentioned above.
Fig. 5. Example from the synthetic imaging routine. Top panel: logarithm of the restmass density in g cm^{−3}. Middle panel: reconstructed radio image at 43 GHz as seen by the VLBA and bottom panels: uv plane (left) and the visibility amplitude (right). See text for details on the jet model used. 
4. Constrained nonlinear 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 nonlinear optimisation problem:
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/nonexistence 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, :
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:
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: gradientbased and gradientfree algorithms. Given the high dimensionality of our problem, together with the high computational effort (raytracing and synthetic imaging) a gradientbased 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 gradientfree 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:
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 socalled 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 (ForemanMackey 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 burnin times during the MCMC simulations, significantly reducing the computational effort and speeding up the calculation. We typically apply 400 random walkers and perform 10^{3} iterations, which leads to a total number of 4 × 10^{5} 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:

Start the optimisation using GA or PSO, typically between 10^{4} and 10^{5} 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 endtoend 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:
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^{4} 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 plateaulike 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)
Best fit parameters obtained from nonlinear optimisation.
Fig. 6. Results of the nonlinear 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, S_{max}, 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^{9} Hz < ν < 10^{12} 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 χ^{2} values in Table 5). Both models exhibit the trend towards lower highfrequency emission (ν > 10^{11} Hz) and higher low frequency emission (ν < 5 × 10^{9} Hz). This behaviour can be attributed to the idealised treatment of the nonthermal particles. In our simulations we neglect radiative losses and reacceleration acting on the nonthermal 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 ν_{1} and ν_{2} is related to the energy distribution of the nonthermal 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 reacceleration of nonthermal 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: 2Ddistribution 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 dashdotted 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^{6}. 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 selfabsorption 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^{18} cm and R_{out} = 1.5 × 10^{18} 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:
where the exponents k_{ρ} and l_{ρ} model the decrease in density in radial and Θ directions. Using the values obtained from the nonlinear optimisation (see Table 5) we calculate a number density of:
Both values are in agreement with number density in NGC 1052 derived from Xray observations of 0.6 × 10^{22} cm^{−2} ⩽ N_{H, obs} ⩽ 0.8 × 10^{22} cm^{−2} (Kadler et al. 2004b).
6. Discussion
6.1. Multifrequency VLBA observations and coreshifts
The inclusion of the broadband radio spectrum (see Fig. 7) in our nonlinear optimisation process enables us to compute various synthetic images in addition to 22 GHz and 43 GHz. Therefore, we can model the multifrequency 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 multifrequency 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 nonthermal 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 multifrequency 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 multifrequency 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 coreshift. 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 nonthermal particle distribution, α_{nt}, along the line of sight:
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 microphysics (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^{7}. 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 frequencydependent variation of the coreshift 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 ∝ ν^{−2}.
Fig. 11. Results of the coreshift 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:
Inserting the temperature and density profiles used for the obscuring torus (and omitting the angular dependence for simplicity):
the variation of the core position with frequency can now be written as:
Inserting the values for k and l obtained from the nonlinear optimisation in Eq. (29) leads to:
With increasing frequency, the absorption due to the thermal particle distribution is decreasing and the nonthermal particles start to dominate the opacity. Following Lobanov (1998) the core position can be written as:
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 nonthermal particles, n_{0} ∝ r^{−n}, with distance. In the case of equipartition (kinetic energy density equals magnetic energy density), Eq. (32) leads to: r_{core, nt} ∝ ν^{−1} using m = 1 and n = 2. If the nonthermal particle distribution dominates the absorption along a line of sight, we can derive estimates for the evolution of the magnetic field and the nonthermal particles using Eqs. (8) and (7). Both equations depend on the pressure: B ∝ p^{1/2} and n_{0} ∝ 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^{−2} 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 coreshift behaviour in the case of nonthermal particles dominating the absorption:
Using the developed estimates for the variation of the core position with frequency, Eqs. (30)–(34), we can explain the coreshift 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 nonthermal particles in the jet starts contributing to the opacity. Therefore, we obtain for 8 GHz < ν < 22 GHz a steepening of the coreshift from ∝ν^{−0.3} to ∝ν^{−1}. 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 coreshift in this region. This behaviour is clearly visible in Fig. 11 for ν > 22 GHz for the eastern jets.
Due to the geometry of the jettorus 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 coreshift 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 frequencydependent position of the τ = 1surface 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 restmass density for the OP jet model (see Fig. 12).
Fig. 12. Distribution of the rest mass density for the OP jet in the zx 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 raytracing and regions of the jet within the dashed area are obscured by the torus and not visible in the corresponding radio images. 
6.2. Overpressured vs. pressurematched 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 spacebased 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 spacebased 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 generalrelativistic 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 overpressured 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 overpressured 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} (pressuremismatch) 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 rerun 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 coreradius 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 coreradius, z_{c}, induces an earlier transition from p_{a}(z)∝z^{−1.3} to p_{a}(z)∝z^{−2} (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.
Best fit parameters obtained from OP jet and OP jet refined model.
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 frequencydependent 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 wellreproduced. 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 reacceleration 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 nonthermal 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 shockacceleration and shear acceleration are able to reenergise the relativistic particles during their propagation through the jet (Gonzalez & Parker 2016; Vaidya et al. 2018; Liu et al. 2017). The reacceleration 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 reacceleration processes.
However, such a selfsimilar treatment of the nonthermal particles during the optimisation process is currently computationally too demanding.
7. Conclusions and outlook
In this work we present an endtoend pipeline for the modelling of relativistic jets using stateofthe art SRHD and emission simulations coupled via evolutionary algorithms to high resolution radio images. We use this newlydeveloped 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 Xray observations. The detailed comparison with available data leads to the conclusion that NGC 1052 is best described by an overpressured 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 timedelays (slowlight 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 reacceleration mechanisms of the nonthermal particles along the flow.
To overcome recent limitations with respect to the magnetisation of the jets we will couple in a followup 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.
Including radiative losses similar to Mimica et al. (2009), Fromm et al. (2016) requires the injection and propagation of Lagrangian particles which is currently numerically too expensive within our modelling algorithm.
This value is derived from the observed peak flux density and flux density of the innermost gaussian component for different frequencies using values given in Table 1 of Kadler et al. (2004a).
For more details see https://www3.mpifrbonn.mpg.de/div/vlbi/globalmm/
Due to the constrain handling in ALPSO the iterations are split into inner (unconstrained) and outer (constrained) iterations see Sect. 4 in Jansen & Perez (2011) for details.
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 AYA201348226C0302P 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
 Aloy, M. A., & Rezzolla, L. 2006, ApJ, 640, L115 [NASA ADS] [CrossRef] [Google Scholar]
 Baczko, A.K., Schulz, R., Kadler, M., et al. 2016, A&A, 593, A47 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Baczko, A.K., Schulz, R., Kadler, M., et al. 2019, A&A, 623, A27 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Böck, M. 2013, Ph.D. Thesis, Dr. Karl RemeisSternwarte Bamberg Erlangen Centre for Astroparticle Physics, Germany [Google Scholar]
 Braatz, J. A., Wilson, A. S., Henkel, C., Gough, R., & Sinclair, M. 2003, ApJS, 146, 249 [NASA ADS] [CrossRef] [Google Scholar]
 Chael, A. A., Johnson, M. D., Narayan, R., et al. 2016, ApJ, 829, 11 [NASA ADS] [CrossRef] [Google Scholar]
 Claussen, M. J., Diamond, P. J., Braatz, J. A., Wilson, A. S., & Henkel, C. 1998, ApJ, 500, L129 [NASA ADS] [CrossRef] [Google Scholar]
 Deb, K., Pratap, A., Agarwal, S., & Meyarivan, T. 2002, Trans. Evol. Comp., 6, 182 [Google Scholar]
 Doeleman, S. S. 2017, Nat. Astron., 1, 646 [NASA ADS] [CrossRef] [Google Scholar]
 Engelbrecht, A. P. 2007, Computational Intelligence: An Introduction (Hoboken, NJ: John Wiley & Sons Inc.) [Google Scholar]
 Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2019a, ApJ, 875, L5 [NASA ADS] [CrossRef] [Google Scholar]
 Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2019b, ApJ, 875, L6 [Google Scholar]
 ForemanMackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306 [CrossRef] [Google Scholar]
 Fromm, C. M., Ros, E., Perucho, M., et al. 2013, A&A, 557, A105 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fromm, C. M., Perucho, M., Mimica, P., & Ros, E. 2016, A&A, 588, A101 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fromm, C. M., Perucho, M., Porth, O., et al. 2018, A&A, 609, A80 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fuhrmann, L., Angelakis, E., Zensus, J. A., et al. 2016, A&A, 596, A45 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gabel, J. R., Bruhweiler, F. C., Crenshaw, D. M., Kraemer, S. B., & Miskey, C. L. 2000, ApJ, 532, 883 [NASA ADS] [CrossRef] [Google Scholar]
 Gómez, J. L., Marti, J. M., Marscher, A. P., Ibanez, J. M., & Alberdi, A. 1997, ApJ, 482, L33 [NASA ADS] [CrossRef] [Google Scholar]
 Gonzalez, W., & Parker, E., eds. 2016, in Magnetic Reconnection, Astrophys. Space Sci. Lib., 427 [CrossRef] [Google Scholar]
 Högbom, J. A. 1974, A&AS, 15, 417 [Google Scholar]
 Jansen, P. W., & Perez, R. E. 2011, Comput. Struct., 89, 1352 [CrossRef] [Google Scholar]
 Kadler, M., Ros, E., Kerp, J., et al. 2002, in Proceedings of the 6th EVN Symposium, eds. E. Ros, R. W. Porcas, A. P. Lobanov, & J. A. Zensus, 167 [Google Scholar]
 Kadler, M., Ros, E., Lobanov, A. P., Falcke, H., & Zensus, J. A. 2004a, A&A, 426, 481 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kadler, M., Kerp, J., Ros, E., et al. 2004b, A&A, 420, 467 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kameno, S., Inoue, M., Wajima, K., SawadaSatoh, S., & Shen, Z.Q. 2003, PASA, 20, 134 [NASA ADS] [CrossRef] [Google Scholar]
 Kardashev, N. S., Khartov, V. V., Abramov, V. V., et al. 2013, Astron. Rep., 57, 153 [NASA ADS] [CrossRef] [Google Scholar]
 Lee, S.S., Lobanov, A. P., Krichbaum, T. P., et al. 2008, AJ, 136, 159 [NASA ADS] [CrossRef] [Google Scholar]
 Lister, M. L., Aller, H. D., Aller, M. F., et al. 2009, AJ, 137, 3718 [NASA ADS] [CrossRef] [Google Scholar]
 Liu, R.Y., Rieger, F. M., & Aharonian, F. A. 2017, ApJ, 842, 39 [NASA ADS] [CrossRef] [Google Scholar]
 Lobanov, A. P. 1998, A&A, 330, 79 [NASA ADS] [Google Scholar]
 Lu, R.S., Roelofs, F., Fish, V. L., et al. 2016, ApJ, 817, 173 [NASA ADS] [CrossRef] [Google Scholar]
 Martí, J. M., Perucho, M., & Gómez, J. L. 2016, ApJ, 831, 163 [NASA ADS] [CrossRef] [Google Scholar]
 Mimica, P., Aloy, M. A., Agudo, I., et al. 2009, ApJ, 696, 1142 [NASA ADS] [CrossRef] [Google Scholar]
 Mizuno, Y., Gómez, J. L., Nishikawa, K.I., et al. 2015, ApJ, 809, 38 [NASA ADS] [CrossRef] [Google Scholar]
 Mizuno, Y., Younsi, Z., Fromm, C. M., et al. 2018, Nat. Astron., 2, 585 [NASA ADS] [CrossRef] [Google Scholar]
 Pogge, R. W., Maoz, D., Ho, L. C., & Eracleous, M. 2000, ApJ, 532, 323 [NASA ADS] [CrossRef] [Google Scholar]
 Pushkarev, A. B., Kovalev, Y. Y., Lister, M. L., & Savolainen, T. 2017, MNRAS, 468, 4992 [NASA ADS] [CrossRef] [Google Scholar]
 Rezzolla, L., & Zanotti, O. 2013, Relativistic Hydrodynamics (Oxford: Oxford University Press) [CrossRef] [Google Scholar]
 Shepherd, M. C. 1997, Astron. Data Anal. Softw. Syst. VI, 125, 77 [Google Scholar]
 Sironi, L., Keshet, U., & Lemoine, M. 2015, Space Sci. Rev., 191, 519 [NASA ADS] [CrossRef] [Google Scholar]
 Vaidya, B., Mignone, A., Bodo, G., Rossi, P., & Massaglia, S. 2018, ApJ, 865, 144 [NASA ADS] [CrossRef] [Google Scholar]
 Vermeulen, R. C., Ros, E., Kellermann, K. I., et al. 2003, A&A, 401, 113 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wang, Z., Bovik, A. C., Sheikh, H. R., & Simoncelli, E. P. 2004, IEEE Trans. Image Process., 13, 600 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Wrobel, J. M. 1984, ApJ, 284, 531 [NASA ADS] [CrossRef] [Google Scholar]
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^{3} 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^{3} cells). The DSSIM is computed between the reference image and an image with lower numerical resolution for example n = 700^{3} 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^{9} Hz and 5 × 10^{10} Hz is due to the torus which leads to higher absorption. The red curve indicates the results for the adaptive linearlogarithmic 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^{3} grid.
A.2. Spectral convergence
In Fig. A.2 we show the singledish spectrum between (10^{7}−10^{13}) 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^{8} Hz and 5 × 10^{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^{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 raygrid exhibits the native resolution of the SRHD simulations, which is 800^{3}.
Fig. A.1. Convergence test for the synthetic images. The synthetic images produced with 800^{3} 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 χ^{2} 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 raytraced 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. 
Used telescopes and SEFD for the GMVA observations at 86 GHz.
A.4. Global Millimetre VLBI Array (GMVA)
The GMVA^{8} 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 onsource 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 uv plane together with the amplitude and phase for the OP jet model are displayed in Fig. A.5.
Used telescopes and SEFD for the Radioastron observations at 22 GHz.
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 uv plane, middle panel the visibility amplitude and the bottom panel the phase with uvdistance. 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. SpektrR and ground array at 22 GHz
The space based satellite SpektrR 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 sourcescans and 95 min offsource 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 onsource 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 uv 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^{9}λ 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 endtoend pipeline and investigate the ability to recover the injected parameters. For the PSO we use 400 particles, 200 outer iterations and 5 inner iterations^{9} 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^{5}. 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 uv plane, middle panel: visibility amplitude and bottom panel: phase with uvdistance. For reasons of clarity only every 50th data point is plotted. 
Injected parameters for the reference model and recovered parameter using PSO and MCMC.
MCMC estimates for the OP jet model.
MCMC estimates for the PM jet model.
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 dashdotted 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 dashdotted 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 dashdotted 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
Injected parameters for the reference model and recovered parameter using PSO and MCMC.
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 restmass density for different ambient medium configurations as indicated by the exponent m and d_{k} (d_{k} = 1 corresponds to a PM jet and d_{k} > 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/R_{j} < 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. Restmass density distribution for different grid resolutions (left: uniform cartesian 300^{3}, centre: uniform cartesian 500^{3} and right: adaptive logarithmic 300^{3}). 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 restmass density in g cm^{−3}. Middle panel: reconstructed radio image at 43 GHz as seen by the VLBA and bottom panels: uv plane (left) and the visibility amplitude (right). See text for details on the jet model used. 

In the text 
Fig. 6. Results of the nonlinear 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, S_{max}, 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: 2Ddistribution 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 dashdotted 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 multifrequency 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 multifrequency radio images. 

In the text 
Fig. 11. Results of the coreshift 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 zx 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 raytracing 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 800^{3} 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 raytraced 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 uv plane, middle panel the visibility amplitude and the bottom panel the phase with uvdistance. 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 uv plane, middle panel: visibility amplitude and bottom panel: phase with uvdistance. 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 dashdotted 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 dashdotted 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 dashdotted 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 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.