Using evolutionary algorithms to model relativistic jets: Application to NGC 1052

TL;DR

This paper develops an end-to-end modeling pipeline combining relativistic hydrodynamics, radiative transfer, and evolutionary algorithms to interpret VLBI observations of relativistic jets in NGC 1052, capturing both fluid dynamics and emission physics.

Contribution

It introduces a novel integrated approach coupling simulations with observational data fitting for relativistic jets, including thermal and non-thermal emission modeling.

Findings

Jet structure in NGC 1052 can be modeled by a slightly over-pressured jet ($d_k extasciitilde1.5$).

The ambient medium pressure decreases with distance, matching observed jet features.

The pipeline accurately reproduces VLBI observations considering array sampling and imaging effects.

Abstract

High-resolution Very-Long-Baseline Interferometry 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 $\nu\sim43$ 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 this paper we investigate the possible physical conditions in relativistic jets of NGC 1052 by directly modelling the observed emission and spectra via state-of-the-art special-relativistic hydrodynamic (SRHD) simulations and radiative transfer calculations. 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 both thermal and non-thermal emission. In order to compare our model to VLBI observations we take into account the sparse sampling of the u-v plane, the array properties and the imaging algorithm. We present for the first time an end-to-end 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/non-thermal 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 over-pressured jet ($d_k\sim1.5$) embedded in a decreasing pressure ambient medium

Tables

Table 1. Table 1: Image parameters for the stacked VLBA images

$\nu$	RMSa	Speak	Stotal	${\Theta }_{\mathrm{maj}}$	${\Theta }_{\min }$	P.A.
GHz	[mJy/beam]	[Jy/beam]	[Jy]	[mas]	[mas]	[$\mathrm{deg}$]
22	0.39 – 1	0.31	0.96	1.25	0.42	$-$13.5
43	0.63 – 1	0.32	0.79	0.72	0.27	$-$11.9
a RMS values are determined in a region of the final map
without significant source flux

Table 2. Table 2: Parameters for the emission simulations. The x 𝑥 x indicates values which will be optimised during the modelling (see Sect. 4 ).

Symbol	Value	Description
scaling parameter
$d_k$	1,1.5, 2.5, 3.5,4.5,5	pressure mismatch at nozzle
$R_j$	$3\times {10}^{16}$ cm	jet radius
$z_c$	$10R_j$	core radius
$n$	1.5	pressure gradient in ambient medium
$m$	2, 3, 4	pressure gradient in ambient medium
$z$	0.005	redshift
$v_j$	$0.5c$	jet velocity
$\hat{\gamma }$	13/9	adiabatic index
${\rho }_a$	$x$ $g{\mathrm{cm}}^{-3}$	ambient medium density
emission parameter
${ϵ}_B$	$x$	equipartition ratio
${ϵ}_e$	$x$	thermal to non-thermal energy ratio
${\zeta }_e$	$x$	thermal to non-thermal number density ratio
${ϵ}_{\gamma }$	1000	ratio between e- Lorentz factors
$s$	$x$	spectral index
$\vartheta$	80∘	viewing angle
torus parameter
$L_{\mathrm{AGN}}$	$1\times {10}^{43}\mathrm{erg}{s}^{-1}$	bolometric luminosity
$R_{\mathrm{out}}$	$x$ $\mathrm{cm}$	torus outer radius
$\theta$	$x$ ∘	torus thickness
$\rho \left(R_{\mathrm{in}}\right)$	$x$ $g{\mathrm{cm}}^{-3}$	torus density at $R_{\min }$
$T_{\mathrm{sub}}$	$x$ K	dust sublimation temperature
$k$, $l$	$x$	exponents for $\rho$ and $T$ distributions

Table 3. Table 3: Used SEFD for the VLBA

frequency [GHz]	1.6	4.8	8	15	22	43
SEFD [Jy]	289	278	327	543	640	1181
data taken from NRAO

Table 4. Table 4: Parameters used for the calculation of the synthetic image presented in Fig.  5 .

Symbol	Value	Description
scaling parameter
$d_k$	2.5	pressure mismatch at nozzle
$R_j$	$3\times {10}^{16}$ cm	jet radius
$z_c$	$10R_j$	core radius
$n$	1.5	pressure gradient in ambient medium
$m$	2	pressure gradient in ambient medium
$z$	0.005	redshift
$v_j$	$0.5c$	jet velocity
$\hat{\gamma }$	13/9	adiabatic index
${\rho }_a$	$1.67\times {10}^{-21}g/{\mathrm{cm}}^3$	ambient medium density
emission parameter
${ϵ}_B$	0.2	equipartition ratio
${ϵ}_e$	0.4	thermal to non-thermal energy ratio
${\zeta }_e$	1.0	thermal to non-thermal number density ratio
${ϵ}_{\gamma }$	1000	ratio between e- Lorentz factors
$s$	2.2	spectral index
$\vartheta$	80∘	viewing angle
torus parameter
$L_{\mathrm{AGN}}$	$1\times {10}^{43}\mathrm{erg}{s}^{-1}$	bolometric luminosity
$R_{\mathrm{out}}$	$1.5\times {10}^{18}\mathrm{cm}$	torus outer radius
$\theta$	40∘	torus thickness
$\rho \left(R_{\mathrm{in}}\right)$	$8.3\times {10}^{-20}g/{\mathrm{cm}}^3$	torus density at $R_{\min }$
$T_{\mathrm{sub}}$	1400 K	dust sublimation temperature
$k$, $l$	1,2	exponents for $\rho$ and $T$ distributions

Table 5. Table 5: Best fit parameters obtained from non-linear optimisation

Symbol	OP jets	PM jets
scaling parameter
$d_k$	1.5	1.0
$R_j$	$3\times {10}^{16}$ cm
$z_c$	$10R_j$
$n$	1.5
$m$	2
$z$	0.005
$v_j$	$0.5c$
$\hat{\gamma }$	13/9
${\rho }_a$	$3.0\times {10}^{-21}g{\mathrm{cm}}^{-3}$	$2.1\times {10}^{-21}g{\mathrm{cm}}^{-3}$
emission parameter
${ϵ}_B$	0.39	0.20
${ϵ}_e$	0.27	0.34
${\zeta }_e$	0.35	0.40
${ϵ}_{\gamma }$	1000	1000
$s$	3.8	3.9
$\vartheta$	80∘	80∘
torus parameter
$L_{\mathrm{AGN}}$	$1\times {10}^{43}\mathrm{erg}{s}^{-1}$
$R_{\mathrm{out}}$	$1.4\times {10}^{18}\mathrm{cm}$	$1.5\times {10}^{18}\mathrm{cm}$
$\theta$	73∘	54∘
$\rho \left(R_{\mathrm{in}}\right)$	$1.1\times {10}^{-19}g{\mathrm{cm}}^{-3}$	$2.6\times {10}^{-20}g{\mathrm{cm}}^{-3}$
$T_{\mathrm{sub}}$	1250 K	1160 K
$k$, $l$	2.4 (2.3), 2.5 (1.5)	2.2 (1.6), 3.8 (1.5)
${\chi }^2$ results
${\chi }_{22\mathrm{GHz}}^2$	4.72	3.88
${\chi }_{43\mathrm{GHz}}^2$	7.30	9.90
${\chi }_{\mathrm{SED}}^2$	1.28	1.94
image metrics
${\mathrm{cc}}_{22\mathrm{GHz}}$	0.98	0.98
${\mathrm{cc}}_{43\mathrm{GHz}}$	0.91	0.85
${\mathrm{DSSIM}}_{22\mathrm{GHz}}$	0.22	0.21
${\mathrm{DSSIM}}_{43\mathrm{GHz}}$	0.10	0.11

Table 6. Table 6: Best fit parameters obtained from OP jet and OP jet refined model

Symbol	OP jets	OP jets (refined)
scaling parameter
$d_k$	1.5
$R_j$	$3\times {10}^{16}$ cm
$z_c$	$10R_j$	$8R_j$
$n$	1.5
$m$	2
$z$	0.005
$v_j$	$0.5c$
$\hat{\gamma }$	13/9
${\rho }_a$	$3.0\times {10}^{-21}g{\mathrm{cm}}^{-3}$	$5.6\times {10}^{-21}g{\mathrm{cm}}^{-3}$
emission parameter
${ϵ}_B$	0.39	0.10
${ϵ}_e$	0.27	0.32
${\zeta }_e$	0.35	0.43
${ϵ}_{\gamma }$	1000	1000
$s$	3.8	3.5
$\vartheta$	80∘	80∘
torus parameter
$L_{\mathrm{AGN}}$	$1\times {10}^{43}\mathrm{erg}{s}^{-1}$
$R_{\mathrm{out}}$	$1.4\times {10}^{18}\mathrm{cm}$	$1.5\times {10}^{18}\mathrm{cm}$
$\theta$	73∘	54∘
$\rho \left(R_{\mathrm{in}}\right)$	$1.1\times {10}^{-19}g{\mathrm{cm}}^{-3}$	$1.0\times {10}^{-19}g{\mathrm{cm}}^{-3}$
$T_{\mathrm{sub}}$	1250 K	1170 K
$k$, $l$	2.4 (2.3), 2.5 (1.5)	1.8 (2.8), 2.9 (1.8)
$N_H$	$0.7\times {10}^{22}{\mathrm{cm}}^{-2}$	$1.1\times {10}^{22}{\mathrm{cm}}^{-2}$
${\chi }^2$ results
${\chi }_{22\mathrm{GHz}}^2$	4.72	4.29
${\chi }_{43\mathrm{GHz}}^2$	7.30	6.90
${\chi }_{\mathrm{SED}}^2$	1.28	1.48
image metrics
${\mathrm{cc}}_{22\mathrm{GHz}}$	0.98	0.98
${\mathrm{cc}}_{43\mathrm{GHz}}$	0.91	0.91
${\mathrm{DSSIM}}_{22\mathrm{GHz}}$	0.23	0.23
${\mathrm{DSSIM}}_{43\mathrm{GHz}}$	0.10	0.10

Table 7. Table 7: Used telescopes and SEFD for the GMVA observations at 86 GHz

Name	code	SEFD [Jy]
Fort Davis	FD	3600
Los Alamos	LA	3100
Pie Town	PT	4100
Kitt Peak	KP	4600
Owns Valley	OV	5800
Brewster	BR	3500
Mauna Kea	MK	4100
North Liberty	NL	4900
Effelsberg	EF	1000
Green Bank	GB	137
Plateau de Bure	PdB	818
Pico Veleta	PV	654
Yebes	YS	1667
Onsala	ON	5102
Methsähovi	MET	17647
Large Millimetre Telescope	LMT	1714
ALMA	ALMA	68

Table 8. Table 8: Used telescopes and SEFD for the Radioastron observations at 22 GHz

Name	code	SEFD [Jy]
Fort Davis	FD	640
Hancock	HN	640
Los Alamos	LA	640
Pie Town	PT	640
Kitt Peak	KP	640
Owns Valley	OV	640
Brewster	BR	640
Mauna Kea	MK	640
North Liberty	NL	640
St. Croix	SC	640
Effelsberg	EF	90
Green Bank	GB	20
Yebes	YS	200
SpektR	RS	30000

Table 9. Table 9: Injected parameters for the reference model and recovered parameter using PSO and MCMC

type	$d_k$	$m$	$\log {\rho }_a$	$\log \rho \left(R_{\mathrm{in}}\right)$	$R_{\mathrm{out}}$	$\Theta$	$T_{\mathrm{sub}}$	${ϵ}_e$	${ϵ}_B$	$k_T$	$k_{\rho }$	$l_T$	$l_{\rho }$	${\zeta }_e$	s
			[g/cm3]	[${\rho }_a$]	$\left[R_j\right]$	[∘]	[K]
ref. model	2.5	2	-20.78	1.70	50.00	40.00	1400.00	0.40	0.20	1.00	1.00	2.00	2.00	1.00	2.20
PSO/GA	2.5	2	-20.67	1.54	35.00	47.00	1283	0.47	0.10	1.20	1.20	1.70	1.70	0.98	2.30
MCMC (mean)	2.5	2	-20.89	1.41	44.18	43.35	1300	0.44	0.14	1.19	1.22	2.05	1.91	0.98	2.41
MCMC ($1\sigma$)	–	–	${}_{}{}^{+0.43}{}_{-0.35}{}^{}$	${}_{}{}^{+0.35}{}_{-0.59}{}^{}$	${}_{}{}^{+4.53}{}_{-12.99}{}^{}$	${}_{}{}^{+15.07}{}_{-17.03}{}^{}$	${}_{}{}^{+208.48}{}_{-276.06}{}^{}$	${}_{}{}^{+0.05}{}_{-0.13}{}^{}$	${}_{}{}^{+0.05}{}_{-0.03}{}^{}$	${}_{}{}^{+0.38}{}_{-0.15}{}^{}$	${}_{}{}^{+0.37}{}_{-0.16}{}^{}$	${}_{}{}^{+0.56}{}_{-0.42}{}^{}$	${}_{}{}^{+0.54}{}_{-0.52}{}^{}$	${}_{}{}^{+0.16}{}_{-0.24}{}^{}$	${}_{}{}^{+0.81}{}_{-0.25}{}^{}$

Table 10. Table 10: MCMC estimates for the OP jet model

type	$d_k$	$m$	$\log {\rho }_a$	$\log \rho \left(R_{\mathrm{in}}\right)$	$R_{\mathrm{out}}$	$\Theta$	$T_{\mathrm{sub}}$	${ϵ}_e$	${ϵ}_B$	$k_T$	$k_{\rho }$	$l_T$	$l_{\rho }$	${\zeta }_e$	s
			[g/cm3]	[${\rho }_a$]	$\left[R_j\right]$	[∘]	[K]
PSO/GA	1.5	2	-20.52	1.55	50	74	1249	0.27	0.39	2.44	2.50	2.28	1.54	0.36	3.86
MCMC (mean)	1.5	2	-20.48	1.54	48	72	1260	0.27	0.40	2.45	2.49	2.28	1.54	0.36	3.75
MCMC ($1\sigma$)	–	–	${}_{}{}^{+0.81}{}_{-1.08}{}^{}$	${}_{}{}^{+0.30}{}_{-0.32}{}^{}$	${}_{}{}^{+11}{}_{-10}{}^{}$	${}_{}{}^{+14}{}_{-15}{}^{}$	${}_{}{}^{+220}{}_{-260}{}^{}$	${}_{}{}^{+0.06}{}_{-0.06}{}^{}$	${}_{}{}^{+0.08}{}_{-0.08}{}^{}$	${}_{}{}^{+0.47}{}_{-0.51}{}^{}$	${}_{}{}^{+0.48}{}_{-0.49}{}^{}$	${}_{}{}^{+0.46}{}_{-0.46}{}^{}$	${}_{}{}^{+0.32}{}_{-0.31}{}^{}$	${}_{}{}^{+0.07}{}_{-0.08}{}^{}$	${}_{}{}^{+0.85}{}_{-0.72}{}^{}$

Table 11. Table 11: MCMC estimates for the PM jet model

type	$d_k$	$m$	$\log {\rho }_a$	$\log \rho \left(R_{\mathrm{in}}\right)$	$R_{\mathrm{out}}$	$\Theta$	$T_{\mathrm{sub}}$	${ϵ}_e$	${ϵ}_B$	$k_T$	$k_{\rho }$	$l_T$	$l_{\rho }$	${\zeta }_e$	s
			[g/cm3]	[${\rho }_a$]	$\left[R_j\right]$	[∘]	[K]
PSO/GA	1	2	-20.67	1.08	45	55	1161	0.34	0.20	2.25	3.79	1.58	1.52	0.42	3.96
MCMC (mean)	1	2	-20.72	1.10	46	53	1179	0.35	0.20	2.25	3.80	1.55	1.51	0.41	3.95
MCMC ($1\sigma$)	–	–	${}_{}{}^{+0.96}{}_{-0.95}{}^{}$	${}_{}{}^{+0.20}{}_{-0.22}{}^{}$	${}_{}{}^{+9}{}_{-10}{}^{}$	${}_{}{}^{+11}{}_{-10}{}^{}$	${}_{}{}^{+195}{}_{-240}{}^{}$	${}_{}{}^{+0.07}{}_{-0.07}{}^{}$	${}_{}{}^{+0.04}{}_{-0.04}{}^{}$	${}_{}{}^{+0.42}{}_{-0.46}{}^{}$	${}_{}{}^{+0.75}{}_{-0.79}{}^{}$	${}_{}{}^{+0.30}{}_{-0.28}{}^{}$	${}_{}{}^{+0.31}{}_{-0.26}{}^{}$	${}_{}{}^{+0.09}{}_{-0.07}{}^{}$	${}_{}{}^{+0.70}{}_{-0.76}{}^{}$

Equations

$$p_{\mathrm{a}}(z)=\frac{p_{\mathrm{b}}}{d_{\mathrm{k}}}{\left[1+{\left(\frac{z}{z_{\mathrm{c}}}\right)}^n\right]}^{-\frac{m}{n}}.$$

$$n\left({\gamma }_{\mathrm{e}}\right)=n_0{\left(\frac{{\gamma }_{\mathrm{e}}}{{\gamma }_{\mathrm{e},\min }}\right)}^{-s}\mathrm{for}{\gamma }_{\mathrm{e},\min }\le {\gamma }_{\mathrm{e}}\le {\gamma }_{\mathrm{e},\max },$$

$${\int }_{{\gamma }_{\mathrm{e},\min }}^{{\gamma }_{\mathrm{e},\max }}n\left({\gamma }_{\mathrm{e}}\right)d\gamma ={\zeta }_{\mathrm{e}}\frac{\rho }{m_{\mathrm{p}}},$$

$${\int }_{{\gamma }_{\mathrm{e},\min }}^{{\gamma }_{\mathrm{e},\max }}n\left({\gamma }_{\mathrm{e}}\right){\gamma }_{\mathrm{e}}{m}_{\mathrm{e}}{c}^{2}d\gamma ={ϵ}_{\mathrm{e}}\frac{p}{\hat{\gamma }-1},$$

$$\gamma_{\rm e,\,min}=\left\{\begin{array}[]{ll}\displaystyle\frac{p}{\rho}\frac{m_{\rm p}}{m_{\rm e}c^{2}}\frac{(s-2)}{(s-1)(\hat{\gamma}-1)}\frac{\epsilon_{\rm e}}{\zeta_{\rm e}}&\textrm{if }s>2\,,\vspace*{3mm}\\ \displaystyle\left[\frac{p}{\rho}\frac{m_{\rm p}}{m_{\rm e}c^{2}}\frac{(2-s)}{(s-1)(\hat{\gamma}-1)}\frac{\epsilon_{\rm e}}{\zeta_{\rm e}}\gamma_{\rm e,\,max}^{s-2}\right]^{1/(s-1)}&\textrm{if }1<s<2\,,\vspace*{3mm}\\ \displaystyle\frac{p}{\rho}\frac{\epsilon_{\rm e}}{\zeta_{\rm e}}\frac{m_{\rm p}}{m_{\rm e}c^{2}(\hat{\gamma}-1)}\left[\ln\left(\frac{\gamma_{\rm e,\,max}}{\gamma_{\rm e,\,min}}\right)\right]^{-1}&\textrm{if }s=2\,.\end{array}\right.$$

$${\gamma }_{\mathrm{e},\max }={ϵ}_{\gamma }{\gamma }_{\mathrm{e},\min }.$$

$$n_0=\frac{{ϵ}_{\mathrm{e}}p(s-2)}{\left(\hat{\gamma }-1\right){\gamma }_{\mathrm{e},\min }^2{m}_e{c}^2}{\left[1-{\left(\frac{{\gamma }_{\mathrm{e},\max }}{{\gamma }_{\mathrm{e},\min }}\right)}^{2-s}\right]}^{-1}.$$

$$B=\sqrt{8\pi {ϵ}_{\mathrm{B}}\frac{p}{\hat{\gamma }-1}}.$$

$$\frac{d{I}_{\nu }}{ds}={ϵ}_{\nu ,\mathrm{nt}}-\left({\alpha }_{\nu ,\mathrm{nt}}+{\alpha }_{\nu ,\mathrm{th}}\right)I_{\nu },$$

$$x_{\mathrm{align}}$$

$$y_{\mathrm{align}}$$

$$z_{\mathrm{align}}$$

$$\sigma =\frac{1}{0.88}\sqrt{\frac{\mathrm{SEF}{\mathrm{D}}_1\times \mathrm{SEF}{\mathrm{D}}_2}{2\Delta \nu t_{\mathrm{int}}}}.$$

$$\displaystyle\begin{array}[]{ll@{}ll}\mathrm{minimize}&\displaystyle f(\@vec{x})&\\ \mathrm{subject\,to}&\displaystyle g_{j}(\@vec{x})\leq 0,&&j=1,...,n\\ &\displaystyle x_{L,i}\leq x_{i}\leq x_{R,i},&&i=1,...,m,\end{array}$$

$$f\left(\text{@vec}x\right)=-{\left\{\left[\sum _{p=1}^q{w}_i{\chi }_{\mathrm{ridge},\mathrm{i}}^2(\text{@vec}x)\right]+w_{2q+1}{\chi }_{\mathrm{spec}}^2(\text{@vec}x)\right\}}^{-2},$$

$$g_{i=1,\cdots ,q}(\text{@vec}x)$$

$$g_{i=q+1,\cdots ,2q+1}(\text{@vec}x)$$

$$g_{i=2q+1,\cdots ,3q+1}(\text{@vec}x)$$

$$\text{@vec}{x}_i(t+1)=\text{@vec}{x}_i(t)+\text{@vec}{v}_i(t+1).$$

$$\rho =\rho \left(R_{\mathrm{in}}\right){\left(\frac{r}{R_{\mathrm{in}}}\right)}^{-k_{\rho }}{e}^{-l_{\rho }|\cos \Theta |},$$

$$N_{\mathrm{H}}$$

$$\tau =\int \left({\alpha }_{th}+{\alpha }_{nt}\right)ds.$$

$${\alpha }_{\mathrm{th},\mathrm{RT}}\propto T^{-3/2}{\rho }^2{\nu }^{-2}.$$

$$T\propto r^{-l},$$

$$\rho \propto r^{-k},$$

$$r_{\mathrm{core},\mathrm{th}}\propto {\nu }^{-1/(0.75l+k-0.5)}.$$

$$r_{\mathrm{core},\mathrm{th}}\propto {\nu }^{-0.27}\mathrm{OP}\mathrm{jet},$$

$$r_{\mathrm{core},\mathrm{th}}\propto {\nu }^{-0.22}\mathrm{PM}\mathrm{jet}.$$

$$r_{\mathrm{core},\mathrm{nt}}\propto {\nu }^{-(5-2\alpha )/\left[\left(2\alpha -3\right)b-2n-2\right]},$$

$$r_{\mathrm{core},\mathrm{nt}}$$

Full text

*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 state-of-the-art special-relativistic hydrodynamic (SRHD) simulations and radiative transfer calculations.

*Methods. *We performed SRHD simulations of over-pressured and pressure-matched jets using the special-relativistic 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 both thermal and non-thermal emission. In order to compare our model to VLBI observations we take into account the sparse sampling of the u-v plane, the array properties and the imaging algorithm.

*Results. *We present for the first time an end-to-end 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/non-thermal 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 over-pressured jet ($d_{k}\sim 1.5$) embedded in a decreasing pressure ambient medium

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. 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. 2004b). The jets are propagating in an east-west direction and are oriented at $\sim 60^{\circ}$ in the sky (measured from North through 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 $\nu\geqslant 43$ GHz (Kadler et al. 2004b). 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\,\mathrm{pc}$ and $0.7\,\mathrm{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^{22}\,\mathrm{cm}^{-2}$ to $10^{24}\,\mathrm{cm}^{-2}$ (Kadler et al. 2004a, b). Kadler et al. (2002, 2004a) 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). The X-ray spectra reveal a strongly absorbed continuum (see, e.g.,  Guainazzi et al. 2000) and hints of a relativistically broadened Fe K$\alpha$ emission originating from a few gravitational radii of the central object (Brenneman et al. 2009). However, Böck (2013) shows that the data of NGC 1052 can also be modelled consistently without a broad Fe K$\alpha$ line.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$\pm$0.06 mas yr*-1* which corresponds to projected speeds of 0.230$\pm$0.011 c, consistent with the values reported by Vermeulen et al. (2003) (0.73$\pm$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\pm 0.037$ c in the western jet and 0.561$\pm$0.034 c in the eastern jet. Using the obtained kinematics the viewing angle, $\vartheta$, can be computed. Vermeulen et al. (2003) report values for the viewing angle of $\vartheta\geq 57^{\circ}$, whereas the values of Böck (2013) constrain it to $\vartheta\geq 78.8^{\circ}$ for a v=0.238 c (see its Sect. 3.1.3). From the speeds measured in the jet and counter-jet, the upper limit of $\vartheta\leq 85^{\circ}$ 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\leq\vartheta\leq 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 i.e., asymmetries in the jet launching/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 perform 2D-axisymmetric SRHD simulations of jets in a decreasing-pressure ambient medium and compute their radiative signatures. Depending on the parameters of the torus i.e., torus density, temperature and dimensions, flux density asymmetries in the jets can be produced. In addition, spectral indices, $\alpha>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 $p=\rho\epsilon(\hat{\gamma}-1)$, where $p$ is the pressure, $\rho$ the rest-mass density, $\epsilon$ the specific internal energy, and $\hat{\gamma}$ the adiabatic index (see, e.g., Rezzolla & Zanotti 2013).

The flux density errors, $\sigma S_{\nu}$, on the synthetic images are computed using a frequency- dependent calibration error $\sigma_{\rm cal}\sim 0.1-0.14$ and the off-source rms: $\sigma S_{\nu}=\sigma_{\rm cal.}S_{\nu}+\rm{rms}$, where $S_{\nu}$ is the flux density. We define the spectral index, $\alpha$, computed between two frequencies, $\nu_{1}$ and $\nu_{2}$ as:

$\alpha_{\nu_{1},\nu_{2}}=\log_{10}\left(S_{\nu_{1}}/S_{\nu_{2}}\right)/\log_{10}\left({\nu_{1}/\nu_{2}}\right)$.

At its distance of $17.72\pm 1.26$ Mpc, an angular separation of 1 mas corresponds to a projected length of about 9.5$\times$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 111https://www.physics.purdue.edu/MOJAVE/ (Lister et al. 2009) and the FGamma program222https://www3.mpifr-bonn.mpg.de/div/vlbi/fgamma/fgamma.html (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 spectrum333We only considered frequencies where at least 20 measurements are available. In the following we provide a short description of the data used in this work.

NGC1052 was observed 29 times between 2004 and 2009 with the VLBA at 22GHz and 43GHz. 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\times 0.2)\,$mas at 43 GHz and $(0.8\times 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 43GHz 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\,\mathrm{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 4 cells per jet radius $(R_{j})$ the numerical grid covers a range of $80\,R_{j}\times 100\,R_{j}$. We inject jets with a velocity of $v_{j}=0.5\,\mathrm{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):

[TABLE]

For the simulations presented in Fromm et al. (2018) we applied $z_{c}=10\,R_{\rm 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)

After roughly 5 longitudinal grid crossing times a steady state is recovered and we insert a torus into our 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 i.e., 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 i.e., pressure, $p$ and density $\rho$. We assume a power law distribution of relativistic electrons:

[TABLE]

where $n_{0}$ is a normalisation coefficient, $\gamma_{\rm e}$ is the electron Lorentz factor, $\gamma_{\rm e,\,min}$ and $\gamma_{\rm 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 $\zeta_{\rm e}$ as:

[TABLE]

with $m_{\rm 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 $\epsilon_{\rm e}$:

[TABLE]

By performing the integrals in Eq. 3 and Eq. 4 an expression for the lower electron Lorentz factor, $\gamma_{\rm e,min}$ can be derived:

[TABLE]

The last step in the construction of the non-thermal particle distributions assumes that the upper electron Lorentz factor, $\gamma_{\rm e,max}$, is a constant fraction $\epsilon_{\gamma}$ of the lower electron Lorentz factor

[TABLE]

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:

[TABLE]

Given that information on the magnetic field cannot be obtained from our hydrodynamical numerical simulations, we assume its magnitude is a fraction $\epsilon_{\rm B}$ of the equipartition magnetic field

[TABLE]

Finally, we compute the frequency-dependent total intensity, $I_{\nu}$, for each ray by solving the radiative transport equation

[TABLE]

where $\epsilon_{\rm nt}$ and $\alpha_{\rm nt}$ are the emission and absorption coefficients for non-thermal emission and $\alpha_{\rm 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

3.3.1 Ray-tracing through an adaptive grid

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).

In order to focus numerical resolution on the jet we first have to align the initial Cartesian coordinate system, indicated by the subscript “$\mathrm{cart}$”, with the jet direction. The aligned coordinate system labeled with the subscript “$\mathrm{align}$” is obtained from the Cartesian system via two rotations using the angles $\vartheta$ (viewing angle) and $\varphi$ (orientation of the jet in plane of the sky). After the rotation into the aligned system the coordinates are modified according to:

[TABLE]

where $\Delta(x,y,z)_{\rm align,\,min}$ is the smallest cell spacing in the different directions, $(x,y,z)_{\rm scale}$ sets the extent of the linear scale i.e., number of cells covered with the highest numerical resolution and $\eta_{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 $\Delta(x,y,z)_{\rm align,\,min}=10^{-3}$, $(x,y,z)_{\rm scale}=4$ and $\eta_{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 $\Delta\tau\simeq 10^{-6}$. This method guarantees that we are tracing the optical depth cut-off with high precision and leads to converged spectra and images (see Appendix 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:

Radiative transfer calculation. 2. 2.

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

Fourier transformation of the computed intensity and sampling with the projected baselines of the observing array. 4. 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 10 identical 25 metre radio telescopes scattered across the USA (see Fig. 4). The antennas are equipped with several receivers allowing multi-frequency observations between 1.6 GHz to 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_{\rm sys}$, and the effective area of the telescope, $A_{\rm eff}$444the effective area is the product of the geometrical telescope area and the aperture efficiency. 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, $\Delta\nu$, and the integration time, $t_{\rm int}$, determine the thermal noise in the synthetic images according to:

[TABLE]

We use a bandwidth of 256 MHz and an integration time of 10 s for all frequencies listed in Table 3. In order to include the calibration gaps in the synthetic radio observations we assume 10 minutes on-source and 50 minutes off-source per observing hour. Together with an integration time of 10 s this transforms into $\sim 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^{\circ}$ an upper limit of $85^{\circ}$. The synthetic data is generated using the EHTim library555https://github.com/achael/eht-imaging which we modified to suit our needs.

Image reconstruction:

The simulated visibilities are imported into DIFMAP Shepherd (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 e.g., 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 R.A. 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.

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:

[TABLE]

where $\@vec{x}$ is an $m$ dimensional vector including the model parameters, e.g., $\@vec{x}=\left[j_{\rm mod},\rho_{a},\epsilon_{b},\epsilon_{e},\zeta_{e},s,R_{\rm out},\theta,\rho(R_{\rm in}),T_{\rm sub},k,l\right]^{T}$ and $f(\@vec{x})$ is the objective function (minimisation function), $g_{j}(\@vec{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, $\chi^{2}_{\rm{ridge,i}}$ and the least squares computed from the observed radio spectrum and simulated one, $\chi^{2}_{\rm spec}$:

[TABLE]

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:

[TABLE]

The first constraints (Eq. 19) set the required minimum cross correlation coefficient between the observed and the synthetic radio images, i.e., rejecting or accepting solutions based on the structural agreement between the images. Equation 20 increases the agreement between the number of the flux density peaks, $n(\@vec{x})_{\rm peak,p}^{\rm obs\,sim}$, in the radio maps and Eq. 21 is enforcing a consistence between the images regarding a possible emission gap and its extent, $\Delta r(\@vec{x})_{\rm gap,p}^{\rm obs,\,sim}$. During the optimisation process we allow for all constraints (eqns. 19–21) a tolerance of 0.01 i.e., a cross correlation coefficient of ${\rm cc_{p}(\@vec{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(\@vec{x}$). More details on the constraint implementation can be found in (Deb et al. 2002; 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)

GAs 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 $\@vec{x}=\left[j_{\rm mod},\rho_{a},\epsilon_{b},\epsilon_{e},\zeta_{e},s,R_{\rm out},\theta,\rho(R_{\rm in}),T_{\rm sub},k,l\right]^{T}$, also referred to as a chromosome, and each entry e.g., $\theta$ 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. 18). 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 $\ll 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 PSO mimics the behaviour of animal swarms, e.g., birds searching for food. A swarm particle can be described by its position vector $\@vec{x_{i}}$ and velocity $\@vec{v_{i}}$ in the parameter space, and the food can be related to the minimisation function, here Eq. 18. The optimisation is driven by the velocity $\@vec{v_{i}}$, which reflects both the individual experience of the particle, commonly referred to cognitive component, and the social experience, i.e., exchange of information between neighbouring particles. Thus, the update on the position of each particle can be written as:

[TABLE]

Depending on the variant of the PSO, different recipes for the calculation of the velocity exist (see Chapter 16 in Engelbrecht 2007, for details). In the so-called global best PSO, the velocity update of a particle, $\@vec{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, i.e., $\Delta f\left(\@vec{x}\right)<\epsilon$, 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^{3}$ iterations, which leads to a total number of $4\times 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:

• •

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

• •

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 accuray, $\Delta\tau$. Typically 100 s are required to compute a radio spectrum consisting of 10 frequencies. This leads to a total required computational time of 300 to 3000 cpus hours. 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 parameters of the reference model are well recovered within $1\sigma$ and the scatter of the solutions provides insights into the uncertainties of the method (see Appendix for details).

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 $\leavevmode\nobreak\ 4\times 10^{4}$ iterations and results of the MCMC can be found in the Appendix.

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 $\pm 2$ mas, whereas the flux density profile of the PM jet is continuously decreasing. With increased resolution, i.e., higher observing frequency, recollimation shocks closer to the jet nozzle can be resolved and additional local flux density maxima appear (see panel f at $\pm 1$ mas) in contrast to the PM jet (see panel h)

The broadband radio spectrum between $10^{9}\,\mathrm{Hz}<\nu<10^{12}\,\mathrm{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 $\chi^{2}$ values in Table 5). Both models exhibit the trend towards lower high-frequency emission ($\nu>10^{11}\,\mathrm{Hz}$) and higher low frequency emission ($\nu<5\times 10^{9}\,\mathrm{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.

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 $\nu_{1}$ and $\nu_{2}$ is related to the energy distribution of the non-thermal particles, $s$, via $s=-(2\alpha-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, e.g., 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.

The spectral index computed from the OP and PM model approximates at large distances the theoretical value of $\alpha=1.4$ (OP jet) and $\alpha=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 for details). The largest difference between our numerical models and the VLBA observations occurs at a distance of $x=\pm 2\,\mathrm{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 calculations666Including 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.. 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 $\alpha>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 $\alpha=1.4$ (OP jet) and $\alpha=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_{\mathrm{in}}=0.2\times 10^{18}\,\mathrm{cm}$ and $R_{\mathrm{out}}=1.5\times 10^{18}\,\mathrm{cm}$ ($0.5\,\mathrm{pc}$) which is within the estimates $0.1\,\mathrm{pc}\leqslant R_{\mathrm{torus}}\leqslant 0.7\,\mathrm{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).

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:

[TABLE]

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

[TABLE]

Both values are in agreement with number density in NGC 1052 derived from X-ray observations of $0.6\times 10^{22}\,\mathrm{cm^{-2}}\leqslant N_{\rm H,obs}\leqslant 0.8\times 10^{22}\,\mathrm{cm^{-2}}$ (Kadler et al. 2004a).

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. 2004b). 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.

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 $\tau=1$ surface, i.e., the onset of the jet. The optical depth $\tau$ is computed from the sum of absorption coefficients for the thermal distribution, $\alpha_{th}$, and non-thermal particle distribution, $\alpha_{nt}$, along the line of sight:

[TABLE]

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. (2004b) we define the observed onset of the jet as the innermost location where the flux density reaches $\sim 20\%$ of the peak flux density in each jet777This 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. (2004b). In Fig. 5 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 $\Delta r\propto\nu^{-0.3}$ and continuously steepens towards $\Delta r\propto\nu^{-2}$.

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\nu\ll kT$) can be written as:

[TABLE]

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

[TABLE]

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

[TABLE]

Inserting the values for $k$ and $l$ obtained from the non-linear optimisation in Eq. 30 leads to:

[TABLE]

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:

[TABLE]

where $\alpha$ is the spectral index and $b$ and $n$ are the exponents of the evolution of the magnetic field, $B\propto r^{-b}$, and the density of the non-thermal particles, $n_{0}\propto r^{-n}$, with distance. In the case of equipartition (kinetic energy density equals magnetic energy density), Eq. 33 leads to: $r_{\rm core,nt}\propto\nu^{-1}$ 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 Eqns. 8 and 7. Both equations depend on the pressure: $B\propto p^{1/2}$ and $n_{0}\propto 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\propto r^{-2}$ for $r>r_{c}$ and $p\propto r^{-4/3}$ for $r\leqslant r_{c}$. Inserting the values reported in Table 5 into Eq. 33 leads to the following core-shift behaviour in the case of non-thermal particles dominating the absorption:

[TABLE]

Using the developed estimates for the variation of the core position with frequency, Eqns. 31–35, 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\,\mathrm{GHz}<\nu<22\,\mathrm{GHz}$ a steepening of the core-shift from $\propto\nu^{-0.3}$ to $\propto\nu^{-1}$. At frequencies $\nu>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 $\nu>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 $\tau=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).

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, i.e., a local flux density maxima, by a drop in opacity (see for example the flux density peak -1 mas in panel h of Fig. 6). A way out of this dilemma can be provided by $\mu$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 $\mu$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 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 $\mu$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.

At 86 GHz the torus is optically thin and a clear view to the central region is obtained, i.e., 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 $\pm 0.5$mas which is not seen in the PM jet. Given the most recent GMVA observations of NGC 1052 there are two bright features next to core at roughly $\pm 0.4$mas (Baczko et al. 2016) which could be interpreted as the recollimation shocks.

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. We can further improve our model for NGC 1052 if we modify the underlying SRHD models with respect to the pressure-mismatch, $d_{k}$, and the gradients in the ambient medium. Therefore, we produce several SRHD simulations changing $d_{k}$ in 0.1 steps and $z_{c}$ in steps of 2. We inject these models in our pipeline and use 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_{\rm j}$ to $z_{c}=8\,R_{\rm j}$, while $d_{k}=1.5$ is not changed. Reducing the core-radius, $z_{\rm c}$, induces an earlier transition from $p_{a}(z)\propto z^{-1.3}$ to $p_{a}(z)\propto 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 $\sim 1$ mas are also better fitted by the 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 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 $\epsilon_{\rm b}<0.5$ during the modelling and optimisation process. Given a viewing angle $\vartheta>80^{\circ}$, 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 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 $\epsilon_{\rm 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 ($\nu>86\,\mathrm{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\approx 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\sim 2.2$, if radiative losses are taken into account. Especially between the jet nozzle and the first recollimation shock, relativistic particles with large $\gamma_{e}$ will suffer radiative losses which will steepen the particle distribution and lead to large spectral slopes. 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, i.e., 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\propto 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 $\epsilon_{b}$ and furthermore enables us to include polarised observations i.e., 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.

Acknowledgements.

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)

Appendix

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.

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 $\mathrm{DSSIM}$, two identical images would have $\mathrm{SSIM}=1$ and $\mathrm{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 e.g., $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. 16. 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\times 10^{9}$ Hz and $5\times 10^{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^{3}$ grid.

Spectral convergence

In Fig. 17 we show the single-dish spectrum between $\left(10^{7}-10^{13}\right)$ 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 $S_{\mathrm{total,n^{3}}}/S_{\mathrm{total,n=800^{3}}}\leq 1\%$. The increased absorption of the obscuring torus leads to a decrease in the flux density between $3\times 10^{8}$ Hz and $5\times 10^{10}$ Hz.

The small variation see in the inset of Fig. 17 is caused by absorption in obscuring torus at high frequencies $\nu>1\times 10^{10}\,\mathrm{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^{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. 18 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\times 0.72\,\mathrm{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. 18 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 $\chi^{2}$ values presented in Tables 5 and 6.

Global Millimetre VLBI Array (GMVA)

The GMVA888for more details see https://www3.mpifr-bonn.mpg.de/div/vlbi/globalmm/ 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 7 and the location of the telescopes can be seen in Fig. 19.

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. 20.

Spektr-R and ground array at 22 GHz

The space based satellite Spektr-R extents the baselines up to 10 Earth radii and provides $\mu as$ resolution. The SEFD of the ground stations and the space antenna are listed in Table 8. 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. 21. 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 8 for station codes). The u-v plane and the visibility amplitude and phase for the OP jet model are plotted in Fig. 22. Notice the large extent up to $12\times 10^{9}\lambda$ in North-South direction which leads to a beam of $160\,{\rm\mu as}\times 20\,{\rm\mu as}$.

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 iterations999Due to the constrain handling in ALPSO the iterations are split into inner and outer iterations see Sect. 4 in Jansen & Perez (2011) for details. refer and for the MCMC we employ 400 random walkers and 1000 iterations. These number lead to similar number of total iterations, namely $4\times 10^{5}$. The results of the parameter recovering test can be seen in Fig. 23 and in Table 9. All histograms show clear maxima and injected model values are within $1\sigma$ of the recovered values.

The spread around the optimal solution found for the OP and PM jet models are presented in Tables 10 and 11 and the histograms for the distribution can be found Fig. 24 – 25.