Inferring AGN jet parameters using Bayesian analysis of VLBI data with non-uniform jet model

TL;DR

This paper introduces a Bayesian method to directly infer physical parameters of AGN jets from VLBI data using an inhomogeneous jet model, avoiding traditional Gaussian template fitting.

Contribution

It presents a novel Bayesian analysis approach for inhomogeneous jet modeling directly from VLBI data, improving parameter inference accuracy.

Findings

Jet is well described by an inhomogeneous conical model.

Results favor an electron-positron jet composition.

Detected counter jet component suggests an external absorber.

Abstract

Physical parameters of AGN jets observed with Very Long Baseline Interferometry (VLBI) are usually inferred from the core shift measurements or flux and size measured at a peak frequency of the synchrotron spectrum. Both are preceded by modelling of the observed VLBI jet structure with a simple Gaussian templates. We propose to infer the jets parameters using the inhomogeneous jet model directly - bypassing the modelling of the source structure with a Gaussian templates or image deconvolution. We applied Bayesian analysis to multi-frequency VLBA observations of radio galaxy NGC 315 and found that its parsec-scale jet is well described by the inhomogeneous conical model. Our results favour electron-positron jet. We also detected a component in a counter jet. Its position implies the presence of an external absorber with a steep density gradient at close ($r=0.1$ pc) distance from the central engine.

Figures

Tables

Table 1. Table 1: Fit statistics for different models

model	${\chi }^2$	$N_{\mathrm{param}}$	DoF	${\chi }_{\mathrm{red}}^2$
CLEAN	25969.8	3537	22311	1.16$\pm$0.005
Gaussians	27433.6	136	25712	1.07$\pm$0.004
Non-uniform jet	29747.9	25	25823	1.15$\pm$0.004

Table 2. Table 2: Parameters of the counter-jet component.

Frequency	Flux density	$r$
GHz	mJy	mas
8.1	47.9$\pm$2.9	-0.173$\pm$0.015
8.4	49.7$\pm$3.1	-0.168$\pm$0.013
12.1	38.3$\pm$1.2	-0.220$\pm$0.008
15.4	30.9$\pm$1.1	-0.185$\pm$0.007

Table 3. Table 3: Parameters of the jet component.

Frequency	Flux density	$r$
GHz	mJy	mas
8.1	11.6$\pm$0.6	4.35$\pm$0.05
8.4	9.5$\pm$0.6	4.39$\pm$0.05
12.1	9.3$\pm$0.6	4.49$\pm$0.04
15.4	3.6$\pm$0.6	4.83$\pm$0.06

Equations

$$T(r_{\mathrm{obs}},d)=\frac{c^2{C}_1(\alpha )D^{0.5}{\nu }_{\mathrm{obs}}^{0.5}{r}_{\mathrm{obs}}^{0.5m}}{8\pi k_B{C}_2(\alpha )(1+z{)}^{0.5}{B}_1^{0.5}{r}_1^{0.5m}(\sin \theta {)}^{0.5m}}(1-e^{-\tau (r_{\mathrm{obs}},d)})$$

$$\begin{gathered} \tau (r_{\mathrm{obs}},d)=C_2(\alpha )(1+z{)}^{-(\alpha +2.5)}{D}^{1.5+\alpha }{K}_1{B}_1^{1.5+\alpha }{\nu }_{\mathrm{obs}}^{-(2.5+\alpha )} \\ \times 2r_1^{n+m(1.5+\alpha )}{\left(\sin \theta \right)}^{n+m(1.5+\alpha )-1}\sqrt{{\varphi }_{\mathrm{app}}^2-{\left(\frac{d}{r_{\mathrm{obs}}}\right)}^2}{r}_{\mathrm{obs}}^{-(n+m(1.5+\alpha )-1)} \end{gathered}$$

$$P(\theta |{\mathbb{V}}_{\mathrm{obs}})\sim P({\mathbb{V}}_{\mathrm{obs}}|\theta )\cdot P_{\mathrm{pr}}(\theta )$$

$$P_{\mathrm{pr}}(\theta )=\prod _{i=1}^{i=N_{\mathrm{params}}}{P}_{\mathrm{pr}}({\theta }_i).$$

$$\begin{array}{r}P({\mathbb{V}}_{\mathrm{obs}}|\theta )=\prod _{b=X,Y,J,U}\prod _{i=1}^{N_{\mathrm{vis},b}}\frac{1}{2\pi {\sigma }_{b,i}^2}\\ \exp \left(-\frac{|g_b{V}_{\mathrm{model},b,i}(\theta )-V_{\mathrm{obs},b,i}{|}^2}{2\cdot {\sigma }_{b,bl(i)}^2}\right)\end{array}$$

$$\begin{gathered} {\beta }_{\mathrm{p}}=\frac{u_e}{u_B} \\ {u}_e=n_{\mathrm{pl}}{m}_{\mathrm{e}}{c}^2⟨\gamma -1⟩(1+k_{\mathrm{u}}) \\ {u}_B=\frac{B^2}{8\pi } \end{gathered}$$

$$\sigma =\frac{B^2}{8\pi (n_{\mathrm{e}}{m}_{\mathrm{e}}+n_{\mathrm{p}}{m}_{\mathrm{p}})c^2}$$

$$\begin{gathered} P_{\mathrm{jet}}=P_{\mathrm{e}}+P_{\mathrm{B}}+P_{\mathrm{p}}= \\ =\pi c\beta (\Gamma r{)}^2(\eta u_{\mathrm{e}}+{\eta }_B{u}_B)+\pi c\beta \Gamma (\Gamma -1)f_{\mathrm{pl}}{n}_{\mathrm{pl}}(m_{\mathrm{p}}+m_{\mathrm{e}})c^3{r}^2 \end{gathered}$$

$$r_{\mathrm{obs}}({\nu }_{\mathrm{obs}})={\left(\frac{\tau }{C_{\tau }(\alpha )A_2{\varphi }_{\mathrm{app}}}\right)}^{-\frac{1}{n+m(1.5+\alpha )-1}}{\nu }_{\mathrm{obs}}^{-\frac{1}{k_r}}$$

$$\eta ({\nu }_{\mathrm{obs}})=(4\pi {)}^{-1}⟨C_1(\alpha )⟩(1+z{)}^{-\alpha }{D}^{2+\alpha }K{B}^{1+\alpha }{\nu }_{\mathrm{obs}}^{-\alpha }$$

$$\kappa ({\nu }_{\mathrm{obs}})=⟨C_2(\alpha )⟩(1+z{)}^{-(\alpha +2.5)}{D}^{1.5+\alpha }K{B}^{1.5+\alpha }{\nu }_{\mathrm{obs}}^{-(2.5+\alpha )}$$

$$⟨C_1(\alpha )⟩=\sqrt{\frac{\pi }{4}}{C}_1(\alpha )\frac{\Gamma \left(\frac{3+\alpha }{2}\right)}{\Gamma \left(\frac{4+\alpha }{2}\right)}$$

$$⟨C_2(\alpha )⟩=\sqrt{\frac{\pi }{4}}{C}_2(\alpha )\frac{\Gamma \left(\frac{7+2\alpha }{4}\right)}{\Gamma \left(\frac{9+2\alpha }{4}\right)}$$

$$C_1(\alpha )=\frac{\sqrt{3}{e}^3}{m{c}^2(2\alpha +2)}{\left(\frac{3e}{2\pi mc}\right)}^{\alpha }\Gamma \left(\frac{\alpha }{2}+\frac{11}{6}\right)\Gamma \left(\frac{\alpha }{2}+\frac{1}{6}\right)$$

$$C_2(\alpha )=\frac{\sqrt{3}{e}^3}{8\pi m^2{c}^2}{\left(\frac{3e}{2\pi mc}\right)}^{\alpha +0.5}\Gamma \left(\frac{6\alpha +5}{12}\right)\Gamma \left(\frac{6\alpha +25}{12}\right)$$

$$L=2r\sqrt{{\varphi }^2-{\left(\frac{d}{r}\right)}^2}/\sin \theta$$

$$\begin{gathered} \tau (r_{\mathrm{obs}},d)=C_2(\alpha )(1+z{)}^{-(\alpha +2.5)}{D}^{1.5+\alpha }K{B}^{1.5+\alpha }{\nu }_{\mathrm{obs}}^{-(2.5+\alpha )} \\ \times 2r\sqrt{{\varphi }^2-{\left(\frac{d}{r}\right)}^2}/\sin \theta \end{gathered}$$

$$\begin{gathered} \tau (r_{\mathrm{obs}},d)=C_2(\alpha )(1+z{)}^{-3}{D}^2{K}_1{B}_1^2{\nu }_{\mathrm{obs}}^{-3} \\ \times {\left(\frac{r}{r_1}\right)}^{-4}2r\sqrt{{\varphi }^2-{\left(\frac{d}{r}\right)}^2}/\sin \theta \end{gathered}$$

$$\begin{gathered} \tau (r_{\mathrm{obs}},d)=C_2(\alpha )(1+z{)}^{-3}{D}^2{K}_1{B}_1^2{\nu }_{\mathrm{obs}}^{-3} \\ \times 2r_1^4(\sin \theta {)}^2\sqrt{{\varphi }^2-{\left(\frac{d}{r}\right)}^2}{r}_{\mathrm{ob}}^{-3} \end{gathered}$$

$$\begin{gathered} \tau (r_{\mathrm{obs}},d)=C_2(\alpha )(1+z{)}^{-3}{D}^2{K}_1{B}_1^2{\nu }_{\mathrm{obs}}^{-3} \\ \times 2r_1^4(\sin \theta {)}^3\sqrt{{\varphi }_{\mathrm{app}}^2-{\left(\frac{d}{r_{\mathrm{obs}}}\right)}^2}{r}_{\mathrm{obs}}^{-3} \end{gathered}$$

$$I_{\nu }({\nu }_{\mathrm{obs}})=\frac{\eta ({\nu }_{\mathrm{obs}})}{\kappa ({\nu }_{\mathrm{obs}})}(1-e^{-\tau })$$

$$I_{\nu }(r_{\mathrm{obs}},d)=\frac{C_1(\alpha )D^{0.5}{\nu }_{\mathrm{obs}}^{2.5}}{4\pi C_2(\alpha )(1+z{)}^{0.5}{B}^{0.5}}(1-e^{-\tau (r_{\mathrm{obs}},d)})$$

$$T(r_{\mathrm{obs}},d)=\frac{c^2{C}_1(\alpha )D^{0.5}{\nu }_{\mathrm{obs}}^{0.5}}{8\pi k_B{C}_2(\alpha )(1+z{)}^{0.5}{B}^{0.5}}(1-e^{-\tau (r_{\mathrm{obs}},d)})$$

$$T(r_{\mathrm{obs}},d)=\frac{c^2{C}_1(\alpha )D^{0.5}{\nu }_{\mathrm{obs}}^{0.5}{r}_{\mathrm{obs}}^{0.5}}{8\pi k_B{C}_2(\alpha )(1+z{)}^{0.5}{B}_1^{0.5}{r}_1^{0.5}(\sin \theta {)}^{0.5}}(1-e^{-\tau (r_{\mathrm{obs}},d)})$$

$$\begin{gathered} \tau (r_{\mathrm{obs}},d)=C_2(\alpha )(1+z{)}^{-(\alpha +2.5)}{D}^{1.5+\alpha }{K}_1{B}_1^{1.5+\alpha }{\nu }_{\mathrm{obs}}^{-(2.5+\alpha )} \\ \times 2r_1^{n+m(1.5+\alpha )}{\left(\sin \theta \right)}^{n+m(1.5+\alpha )-1}\sqrt{{\varphi }_{\mathrm{app}}^2-{\left(\frac{d}{r_{\mathrm{obs}}}\right)}^2}{r}_{\mathrm{ob}}^{-(n+m(1.5+\alpha )-1)} \end{gathered}$$

$$T(r_{\mathrm{obs}},d)=\frac{c^2{C}_1(\alpha )D^{0.5}{\nu }_{\mathrm{obs}}^{0.5}{r}_{\mathrm{obs}}^{0.5m}}{8\pi k_B{C}_2(\alpha )(1+z{)}^{0.5}{B}_1^{0.5}{r}_1^{0.5m}(\sin \theta {)}^{0.5m}}(1-e^{-\tau (r_{\mathrm{obs}},d)})$$

$$\tau =\log (1+\frac{n+m(1.5+\alpha )-1}{0.5m}\tau )$$

$$A_1={\left(\frac{D}{B_1(\sin \theta {)}^m}\right)}^{0.5}$$

$$A_2=K_1(D{B}_1{)}^{1.5+\alpha }(1+z{)}^{-(\alpha +2.5)}(\sin \theta {)}^{n+m(1.5+\alpha )-1}$$

$$T=C_T(\alpha )A_1{\left(\frac{{\nu }_{\mathrm{obs}}{r}_{\text{obs}}^m}{(1+z)}\right)}^{1/2}(1-e^{-\tau })$$

Full text

††pubyear: 2019††pagerange: Inferring AGN jet parameters using Bayesian analysis of VLBI data with non-uniform jet model–D

1 Introduction

Jets observed in AGNs with Very Long Baseline Interferometery (VLBI) typically have one-sided structure with an optically thick base (called the VLBI core) that has a flat or even inverted spectrum. The position of core was found to be changing with frequency (Marcaide & Shapiro, 1984). The most plausible explanation is that the core at a given frequency is the surface of the unit optical depth to synchrotron self-absorption ($\tau_{\rm ssa}\approx 1$) as predicted by the non-uniform jet model (Blandford & Königl, 1979; Konigl, 1981). However, other explanations do exist (e.g. standing shock, Marscher, 2009). The non-uniform (or inhomogeneous) jet model was initially invoked to explain flat spectra and variability of radio sources. It assumes conical relativistically moving jet with bulk motion Lorenz factor $\Gamma$ and half-opening angle $\phi$ observed at viewing angle $\theta$. Magnetic field (assumed to be tangled) and amplitude of the particle density $K$ (both measured in a plasma rest frame) depend on the distance from cone apex as $K=K_{\rm 1}(r/r_{\rm 1})^{-n}$ and $B=B_{\rm 1}(r/r_{\rm 1})^{-m}$, where $r_{\rm 1}=1$ pc, $B_{\rm 1}$ and $K_{\rm 1}$ are the values at the $r_{\rm 1}$. Energy distribution of particles in the plasma rest frame is assumed to be a power law $N(\gamma)=K\gamma^{-s}$ where $\gamma$ is the particle Lorenz factor. Such steady state distribution along the jet neglects effects of the radiating cooling on the shape of distribution. This could correspond to constant re-acceleration process acting across the jet volume (Konigl, 1981). At the same time radially decreasing magnetic field can affect the radiative cooling rate of emitting particles at some part of the jet. This could lead to $n\approx 1$ (Malyshev et al., 2013) that is different from the canonical value $n=2$ expected from the constant speed conical isothermal jet (Blandford & Königl, 1979). Exponent $s$ is related to the optically thin spectral index $\alpha$ as $\alpha=(s-1)/2$ taking the following definition of the spectral flux density $S_{\nu}\approx\nu^{-\alpha}$. One of the most successive applications of non-uniform model is the explanation of the aforementioned effect of the core position depending on the observing frequency — the core shift effect (e.g. Lobanov, 1998; Sokolovsky et al., 2011).

The brightness temperature distribution in a general case of arbitrary $n$, $m$ and $\alpha$ is (see Appendix B for details):

[TABLE]

where the optical depth is:

[TABLE]

Here the position in a sky is described by $r_{\rm obs}$ and $d$ — distances along and perpendicular to the jet direction, $D$ – Doppler factor, $\phi_{\mathrm{app}}$ is the observed jet half-opening angle, $z$ – redshift, $C_{1}(\alpha)$ and $C_{2}(\alpha)$ – functions of optically thin spectral index $\alpha$, $c$ – speed of light and $k_{\rm B}$ – Boltzmann constant.

Couple of additional assumptions were employed to derive these relations (Blandford & Königl, 1979) (see also Appendix A). The first one is purely geometrical approximation of the ray path length through the jet. The second is that the absorption coefficient in the jet is assumed to be constant along the line of sight. Both approximations are reasonable when the ratio of jet viewing angle to jet half-opening angle $\theta/\phi\sim\phi_{\mathrm{app}}$ is small.

Lobanov (1998) have shown that the inhomogeneous jet model can be used to estimate the physical parameters of the jets, e.g. magnetic field and particle density, using the measured core shift. However the inference based on the inhomogeneous model needs some further assumptions because core shift measurements alone do not constrain the model parameters. The most widely used is the assumption of equipartition between particle and magnetic field energy densities (Lobanov, 1998). Most of the magnetic field estimates using core shift measurements are made under this assumption (Lobanov, 1998; O’Sullivan & Gabuzda, 2009; Sokolovsky et al., 2011; Pushkarev et al., 2012). Zdziarski et al. (2015) used optically thick approximation to avoid the equipartition assumption and to tie particles number density with the observed flux density (for a similar approach see also Nokhrina, 2017). They found that magnetic field estimates obtained without equipartition assumption significantly differ from corresponding equipartition estimates. Both papers used fixed $s=2$, that corresponds to optically thin spectral index $\alpha=0.5$. Considering strong dependence of the model predictions on this particular parameter (Appendix B) it is desirable to simultaneously infer it from the observed data. At the same time even these simplified approaches need measurements of the core shifts. However core shift estimates obtained using the fitting of the VLBI core with a Gaussian are biased with bias that depends on resolution (Plavin et al., 2019) and observed jet parameters (Pashchenko et al., in prep.).

Finke (2019) used analytical expressions for flux density and core shift of the Blandford & Königl (1979) model to estimate physical jets properties. They used fixed exponents of the magnetic field and particle density radial dependence and used jet power as an additional observational constrain. The resulting estimates have large uncertainties due to low accuracy of the available core shift estimates and unknown value of the particle energy spectrum exponent $s$, which is poorly constrained by the data they used. Fromm et al. (2019) used special-relativistic hydrodynamic (SRHD) simulations as a jet model to fit multifrequency VLBI observations of the radio galaxy NGC 1052. They used images stacked over several epochs to smooth out numerous inhomogeneities travelling along the flow and minimised the difference between the observed stacked and model images at several frequencies. Although they used artificially created VLBI images obtained with the same VLBI array and visibility noise as in the observed data, comparing in the image plane – where the noise is correlated – complicates the uncertainty analysis of the fit. Furthermore, the optimization approach could suffer from the local minima and miss separated high probability parameter regions.

In this paper we use the inhomogeneous jet model to make inference from the radio interferometric visibility data directly – bypassing the modelling of the source structure with a Gaussian templates or image deconvolution.

Throughout the paper, we adopt the $\Lambda$CDM cosmology with $\Omega_{m}=0.287$, $\Omega_{\Lambda}=0.7185$ and $H_{0}=69.3$ km s*-1* Mpc*-1* (Hinshaw et al., 2013).

2 Data

We searched for the source with a featureless straight jet and not small jet viewing angle. The first requirement is needed because the model of Blandford & Königl (1979) does not account for inhomogeneities in a jet, e.g. shocks. However this is not a principal obstacle for using our approach as shocks can be described as additional components over the top of the smooth underlying jet. The model we propose does not handle jet curvature as-is, so for now we need to limit our selection to straight jets. The last geometric requirement of relatively high jet viewing angle implies that the chosen object should be a radio galaxy (Urry & Padovani, 1995). It is due to one of the approximations in the original paper and can be completely abandoned by numerically solving corresponding radiative transfer instead of using the analytical solution (1).

We choose radio source NGC 315 (Bridle et al., 1976) which is a close ($z=0.0165$) giant FR I (Fanaroff & Riley, 1974) radio galaxy. It resides in elliptical galaxy and optically classified as broad-lined LINER (Ho et al., 1997). Flux density monitoring of NGC 315 with Westerbork synthesis radio telescope at 5 GHz from 1974 to 1980 revealed constant flux (Ekers et al., 1983). Cotton et al. (1999) noted the flare that was observed using VLA in 1990–1995 that increased the arcsecond core flux 1.5x (from 588 mJy in 1990 Sep to 746 mJy in 1994 June, Venturi et al., 1993; Cotton et al., 1999). OVRO (Richards et al., 2011) 15 GHz light curve and MOJAVE (Kovalev et al., 2005) 15.4 GHz VLBA-scale fluxes are presented in Figure 1. Canvin et al. (2005) modelled kpc-structure of this source with a relativistic decelerated jet model and estimated viewing angle $\theta=38\pm 2\deg$ with initial speed $\beta\approx 0.9$ (in units of speed of light $c$). On epoch (2006–02–12) the source was observed at VLBA at several frequencies (15.4, 12.1, 8.4 and 8.1 GHz – u, j, y and x bands) in MOJAVE survey (Lister et al., 2009a). Corresponding self-calibrated interferometric visibilities were taken from the MOJAVE database111http://www.physics.purdue.edu/astro/MOJAVE/sourcepages/0055+300.shtml.

3 Bayesian inference and probabilistic model

3.1 Re-parametrization and additional assumptions

The model of the observed brightness distribution from (1) besides jet parameters $\theta$, $\phi$, $B_{1}$, $K_{1}$, $m$, $n$, $\alpha$ includes $dx$ and $dy$ – the location of the jet apex in the image, and $\theta_{\rm rot}$ – jet position angle on the celestial sphere.

Due to degeneracies between model parameters (e.g. $B_{1}$ and $K_{1}$, $\theta$ and $\phi$) we introduce a different parametrization (29 and 30) and use the observed half-opening angle $\phi_{\rm obs}$. Thus our re-parametrized model has the following parameters: $dx$, $dy$, $\theta_{\rm rot}$, $\phi_{\rm app}$, $A_{1}$, $A_{2}$, $m$, $n$, $\alpha$. There is still an unidentifiable combination of parameters $n+m(1.5+\alpha)-1$ in the exponent of the observed position in the expression for the optical depth. This, together with $A_{2}$, defines the characteristic size of the core and thus the core shift effect. However the optically thin spectral index $\alpha$ can be constrained by fitting data at several frequencies simultaneously. The total flux spectrum also constrains the exponents $n$ and $m$ (Konigl, 1981). However we found that this weakly bounds parameters $m,n$ and $\alpha$. Fitting the model with arbitrary exponents $n$ and $m$ to the multifrequency data shows that the degeneracy between them still remains. This could results in poor exploration of their distribution so this item deserves additional investigation. Thus we further assume that the particles to the magnetic field energy density ratio is constant, i.e. $n=2m$. This assumption is weaker than the equipartition assumption. However, it could be violated during flares (Lobanov & Zensus, 1999; Lisakov et al., 2017; Plavin et al., 2019), when region of e.g. higher plasma density is travelling along the jet. NGC 315 light curve presented in Figure 1 shows only a slowly rising VLBA-scale flux increasing by 20% from epoch 2001 to 2009. Thus, the proposed assumption seems reasonable. Moreover, as shown in (Hovatta et al., 2014) distribution of the core spectral indices in MOJAVE sample peaks at $\alpha=0$ thus supporting the assumption of a constant particle-to-magnetic energy density ratio ($n=2m$).

We stress that the re-parametrization (29 and 30) is not necessary, i.e. its usage does not bring any new information, nor bias the estimates of the physical parameters (e.g. $B_{1}$ and $K_{1}$) in any way. One could make the inference using original model with parameters that include $B_{1}$ and $K_{1}$. But using re-parametrized model without degeneracies significantly speeds up the inference (i.e. sampling high dimensional posterior distribution of the model parameters - see Section 3.2). Also sampling of high dimensional distributions with pronounced degeneracies is still a challenge for most of the algorithms in wide usage.

As we fit the model to self-calibrated data, we assume that gains of the individual antennas are accounted for during self-calibration process (Thompson et al., 2017). However the remaining amplitude scale factors should be taken into account, because they affect the spectrum and thus our inference. The corresponding absolute calibration uncertainty was estimated by Hovatta et al. (2014) as $\sigma_{\rm b}=5\%$ for x, y and u-bands and 7.5$\%$ for j-band. The scale factors are degenerate as their simultaneous multiplication by the same number is equivalent to scaling $A_{1}$ parameter with this number (see 29 and 30). To break this degeneracy we put an additional constrain on amplitude scale factors: $\prod_{b}g_{b}=1$. However, this leaves an uncertainty of the product not being exactly unity. The remaining uncertainty due to scaling each of the $g_{b}$ to the same factor can be accounted for using priors for $g_{b}$ after the sampling (see Appendix D for details of our implementation). Another option is to fix one amplitude scale factor to 1 (Natarajan et al., 2017). However this makes it difficult to account for the uncertainty of other model parameters due to this chosen scale factor not being exactly 1.

Finally, our preliminary fits revealed a significant component in a counter-jet side that is fundamentally not accounted by the model described above. Emission on the counter-jet side was detected by examining the residuals in terms of visibilities and reconstructed CLEAN maps. To account for this emission we extended our model by including a point-like component located on the jet axis at some distance from the apex. We used a point component because the sizes of the circular Gaussian components if used are much less than a beam size and weakly constrained. We treat this as component being unresolved. Thus our extended model has 25 parameters in total: 6 non-uniform jet model parameters, jet position angle, location of the jet apex (2 for each of 4 bands), 3 parameters of the amplitude scale and 8 parameters of the counter-jet component (2 for each band).

3.2 Bayesian inference and sampling

Posterior distribution of the model parameters $\mathbb{\theta}$ is expressed as the product of the likelihood and prior distribution (Bayes Theorem, Bayes & Price, 1763; Laplace, 1774):

[TABLE]

where $\mathbb{V}_{\rm obs}$ – the observed data and the prior distribution is itself a product of priors for individuals model parameters:

[TABLE]

We fit the model to the self-calibrated visibility and assume the gaussian noise in its real and imaginary parts. Thus, the likelihood:

[TABLE]

where $g_{b}$ – amplitude scale factor for frequency band $b$, constrained as described in Section 3.1, $V_{{\rm model},b,i}(\mathbb{\theta})$ – Fourier transform of a model at frequency band $b$ with parameter vector $\mathbb{\theta}$ at $i$-th ($u$, $v$)-point, $V_{{\rm obs},b,i}$ – observed visibility function at band $b$ at $i$-th (u, v)-point, $\mathbb{V}_{{\rm obs},b}=(V_{{\rm obs},b,i}),i=\overline{1,N_{{\rm vis},b}}$ – vector of the observed visibilities at frequency band $b$, $\sigma^{2}_{b,bl(i)}$ – dispersion of normal noise on $bl(i)$-th baseline for real/imaginary part of the Stokes I visibility222Stokes I visibilities are calculated as half sum of parallel hand correlations, i.e. $0.5\cdot(<RR*>+<LL*>)$ where L and R are voltages from left and right circular polarized feeds, angular brackets means correlation and star denotes complex conjugation (Thompson et al., 2017) at frequency band $b$. We estimated them with the successive differences approach (Briggs, 1995) for each baseline.

We used non-uniform pixel size for calculating the model image with pixel size increasing from 1 $\mu$as close to the jet apex to 0.1 mas at 20 mas along the jet. This is justified by the fact that the model we use has smooth structure at large scales and more compact details in the core region. Also this significantly speeds up calculation of the model image and Fourier Transform to the visibility space, compared to using the same small pixel size for the whole jet.

We used normal priors for amplitude scale factors $g_{b}$, positions $dx_{b},dy_{b}$, angles $\phi_{\mathrm{app}},\theta_{\mathrm{rot}}$, uniform for positions of the counter-jet component $r_{b}$ and lognormal for $m$, $\alpha$, $A_{1},A_{2}$, fluxes of the counter-jet component $F_{b}$ and checked that our sampled posterior is not bounded by the corresponding priors. We used wide uninformative priors for amplitude scale factors $g_{b}\sim N(1,0.1)$ in sampling and more informed priors $N(1,\sigma_{\rm b})$ in post-sampling step of accounting for the error due to unknown mean scale factor (Appendix D).

The sampling from posterior distribution was done using a nested sampling (Skilling, 2004) algorithm as implemented in PoLyChord (Handley et al., 2015a, b). It implements constrained sampling from prior distribution using slice sampling and allows effectively sample multimodal high dimensional posterior distributions with linear degeneracies. The resulting posterior distribution of the parameters is presented in Figure 2. Position of the cone apex and parameters of counter-jet component are only shown for 15.4 and 8.1 GHz for compactness. Here the prior distributions of individual parameters are plotted as the orange lines on diagonal plots together with histograms of the marginalized posterior distributions.

Finally, we note that contrary to fitting in the image plane (e.g. Fromm et al., 2019), ($u$, $v$)-plane fitting allows to make inference even if the image deconvolution is not feasible (e.g. in Space VLBI observations with sparse coverage of the ($u$, $v$)-plane). Moreover, sampling the full multidimensional posterior distribution of the model parameters ensures to detect possible high-probability modes of the posterior that could be missed by optimisation algorithms. This also provides justified uncertainty estimates of the model parameters.

4 Model evaluation

4.1 Comparing with the observed visibilities and residuals images

Visibility amplitude and phase dependence on radial distance in ($u$, $v$)-plane is presented for the observed and model data for 8.1 and 15.4 GHz in Figure 3. To account for the uncertainty in estimated model parameters we plotted model visibilities not only for the single “best” value of the parameter vector (e.g. Maximum Likelihood) but for a sample of the parameters (24 samples) drawn from the posterior distribution (Figure 2). We also added the observed noise to the model visibilities for each sample from the posterior to aid the comparison with the observed data. This step is necessary because noise contributes significantly to the amplitude of the observed visibilities at the largest baselines with low signal-to-noise ratio. We also accounted for the uncertainty of the mean amplitude scale not being exactly 1.

The original CLEAN images superimposed with the CLEAN images of the differences between the observations and our best model333We use Maximum Likelihood parameters here. are plotted in Figure 4. The most obvious is a residual component in the jet at DEC$\approx$4 at all bands. We discuss it in Section 6.3 further. Also it is apparent that the model fits the inner jet region better at lower frequencies.

4.2 Comparing with CLEAN and Gaussians models

As a part of the model evaluation procedure we compared our model with a CLEAN model and a model consisting of several circular Gaussian components. CLEAN and model fitting were done in Difmap package (Shepherd, 1997). Our model has 25 parameters and describes the observed data at 4 frequency bands simultaneously. To describe the same data using model with circular Gaussian components one needs $\sim 100$ parameters. Corresponding CLEAN model has $\sim$1000 parameters at each frequency band (position and flux of $\delta$–functions). Various fit statistics, including sum of the $\chi^{2}$, reduced chi-squared $\chi_{\rm red}^{2}$ over all 4 bands for CLEAN, Gaussians and our final model are presented in Table 1444Note the critique of using $\chi_{\rm red}^{2}$ (Andrae et al., 2010). Despite the model consisting of Gaussians has slightly lower $\chi^{2}_{\mathrm{red}}$, it has unclear physical interpretation considering the smooth jet structure. The inherently unphysical but highly flexible CLEAN model has the same $\chi^{2}_{\mathrm{red}}$ as the non-uniform jet model for this source.

4.3 Comparing with the observed data - Posterior Predictive Check

Adequate model should create data similar to the observed data (Davies, 1995; Lindsay & Liu, 2009). Posterior predictive check (PPC, Gelman et al., 2013) can be used as a measure of discrepancy that includes the uncertainty associated with the estimated model parameters. PPC consists in simulating the data under the fitted model and then comparing these with the observed data to discover possible systematic differences.

We made the posterior predictive check for the total CLEAN model flux by sampling parameters from the posterior (Figure 2) 200 times, each time constructing the model image, transforming it to ($u$,$v$)-plane, applying corresponding amplitude scaling factors, adding the noise estimated from the observed visibilities and CLEANing with the same parameters (image and pixel size) as for the observed data using Difmap package. All observed values lie in middle quartiles of the posterior predictive distributions. The corresponding percentiles of the observed CLEAN model flux among the posterior predictive distributions are 37%, 65%, 45% and 64% for frequencies 15.4, 12.1, 8.4 and 8.1 GHz correspondingly. Thus, the model describes the total flux well.

Hovatta et al. (2014) found spectral index in the jet region to be $\approx 1.1$ with a typical error $\approx 0.2$. Notably they used the component we interpret being on the counter-jet side as the core to align the images at different frequencies. However this should not influence the jet spectral index. Our modelling results in the median value for $\alpha=0.95$ with 68%-credible interval – (0.91, 0.99), consistent with theirs. We conducted posterior predictive check for spectral index in the following way. We chose 200 random samples from the posterior distribution of model parameters (Figure 2). Then for each sample we simulated the observed data in the same way as for the total flux PPC. To align CLEAN images at different frequencies we shifted the phases of model visibilities and put the position of the jet apex at each frequency to the phase centre. The same shifts were also applied to the real data. After CLEANing of all data sets with the same pixel size and restoring beam we obtained many realizations of the spectral index maps for the model and observed data (conditioning on the current model parameters). The difference between the observed and model spectral index normalized by the corresponding standard deviation is presented in Figure 5. Apparently, the model describes the observed spectral index distribution well besides the barely significant difference in the region of inner jet at $\approx 2$ mas from the cone apex. Here the observed spectral index is flatter than the model one. To show the magnitude of the effect we plot the posterior distribution of the difference between the observed and model spectral index map along the model jet axis in Figure 6. Corresponding difference in the on-axis spectral index is only $\approx$0.05.

(Pushkarev et al., 2017) obtained the median value of half-opening angle $\phi_{\rm app}=3.45\pm 0.05$ degrees for 15.4 GHz VLBA images stacked over 14 epochs, while our estimate for epoch 2006-02-12 is nearly two times larger – median 6.61 with 95% credible interval (6.27, 6.94). There are several possible causes of such inconsistency. (Pushkarev et al., 2017) used deconvolved FWHM of the Gaussians fitted to the transverse slices of the observed brightness distribution to calculate jet width and half-opening angle, while our estimate concerns intrinsic jet geometry. Also they possibly used the counter-jet component as a core, making $\phi_{\rm app}$ smaller and analyzed stacked image.

(Pushkarev et al., 2017) found exponent of the jet width $d$ radial dependence $k=0.86\pm 0.01$, where $d\propto r^{k}$. Although we can not compare to their results directly as they used an image stacked over many epochs, we conducted PPC using the jet shape at single epoch used in our analysis. We generated artificial data from the posterior distribution of model parameters (Figure 2) and added noise as in the original data. Then we CLEANed each artificial data set and convolved corresponding CLEAN models with the same circular beam. For each of the obtained image we applied the method of (Pushkarev et al., 2017) to calculate shape parameter $k$. The posterior predictive distribution together with the observed value are presented in Figure 7. The width of the distribution is large but the observed value lies at its low tail (with $p$-value 0.06 for the corresponding two-sided statistical test). Interesting that the PPC distribution is well below the model value $k=1$. This controversy deserves further investigation. It could imply that the conclusions from such jet shape measurements to the outflow geometry should be treated with caution.

5 Inferred physical properties of the jet

5.1 Spectral index and particle energy distribution

We obtained the median value of the optically thin spectral index $\alpha=0.95$ and 68% credible interval – (0.91, 0.99). This results in a power law exponent of the emitting electrons $s=2.90$ with 95% credible interval (2.75, 3.07). It is significantly higher than $s=2.0$ assumed in the original paper (Blandford & Königl, 1979) and also significantly differs from both $s=2.5$ attributed to weakly magnetized shock acceleration and ultra-relativistic limit $s=2.2$ (Sironi et al., 2015).

However in strong magnetic field the emitting particles are cooled efficiently (Kardashev, 1962) and the break in their energy spectrum occurs at some $\gamma_{\rm br}$. Above the break $\gamma>\gamma_{\rm br}$ the power law steepens from the original (acceleration) value $s_{\rm acc}$ to $s_{\rm acc}+1$. If we observe the steepen part of corresponding synchrotron spectrum then $s_{\rm acc}=s-1=1.90$ with 95% credible interval (1.75, 2.07). This also significantly differs from the ultra-relativistic shocks universal value. One can estimate the break frequency $\nu_{\rm br}$ at given location $r$ by equating the jet travel time to a distance $r$ with synchrotron cooling time (Blandford & Königl, 1979; Konigl, 1981). We obtain the median $\nu_{\rm br}=387$ GHz with 68%-credible interval (217, 701) [GHz] at the position of the VLBI core at 15.4 GHz. This implies that the observed emission resides well below the break and, thus, our model is self-consistent.

Venturi et al. (1993) observed optically thin spectral index $\alpha_{5}^{8.4}\approx 0.5$ between 5 and 8.4 GHz and sharp steepening at $r\approx 4$ mas from the core using VLBI observations at close epochs (1989 Apr – 1990 Nov). They conclude that the break frequency $\nu_{\rm br}\approx 5$ GHz at distance $r$ = 4 mas from the core. Thus at our frequency bands (i.e. above $\nu_{\rm br}$) we should see a steepen spectrum at this location. Interesting that the steepening is observed at 6 mas in our residuals between data and model (Figure 6). However, Hovatta et al. (2014) showed that for steepening with distance the jet should be collimating with radius $R\propto r^{k}$ where $k<2/3$ assuming magnetic field with dominating transverse component. Another explanation is a cutoff in synchrotron spectrum due to high energy cutoff $\gamma_{\rm max}$ of the electron energy spectrum (Hovatta et al., 2014). This could imply $\gamma_{\rm max}=370\pm 20$ at this distance assuming that emission at given frequency 8.1 GHz is dominated by electrons with $\gamma_{\rm rad}\approx\sqrt{(\nu_{\rm 8.1GHz}/D)/\nu_{B_{\rm 6mas}}}=370$, where $\nu_{B}$ - Larmor frequency for magnetic field at a given radial distance.

5.2 Magnetic field and particle density

We obtained values of the exponents $m$ and $n$, that determine the radial dependence of the magnetic field and particle number density, close to the canonical (Blandford & Königl, 1979): $m=0.94\pm 0.01$ and $n=1.88\pm 0.02$ (68% credible intervals). The value $m\approx 1$ implies the dominant transverse component of the magnetic field.

To infer the value of the fields we need an independent estimate of the Doppler factor as it is not constrained in our model. Velocity $\beta$ could be estimated to lie in a range 0.7–0.96 of component speeds observed in pc-scale radio structure (Cotton et al., 1999). Together with estimate of the jet viewing angle from Canvin et al. (2005) we used the obtained posterior distribution of model parameters (Figure 2) to calculate magnetic field and particle density at $r=1$ pc (Figure 8). Here the lines correspond to the equipartition and are shown for $\gamma_{\rm min}$ = 10 and 100. When using relation between the jet opening angle and Lorentz factor derived in Clausen-Brown et al. (2013), the median and 95% credible intervals are -1.41 and (-1.53, -1.29) for $\lg{B_{1}{\rm[G]}}$ and 3.50 and (3.0, 4.18) for $\lg{K_{1}[{\rm cm^{-3}]}}$.

It is interesting to compare this estimate of the magnetic field with the traditional estimate that assumes equipartition (e.g. Lobanov, 1998; Hirotani, 2005; O’Sullivan & Gabuzda, 2009). Pushkarev et al. (2012) obtained the following values of the core shift projected on the jet median position angle: 0.179 mas for 15.4–8.1 GHz555The posterior distribution of the model core shift (see Section 6.1 and Appendix B) has median 0.165 mas and 95% credible interval (0.159, 0.171), 0.064 mas for 15.4–8.4 GHz and 0.052 mas for 15.4–12.1 GHz. To calculate the magnetic field value one would assume equipartition (thus $k_{r}=1$, Lobanov, 1998) and $\alpha=0.5$. With core shift value for 15.4–8.4 GHz that is consistent with one at higher frequency, we use equation (5) from (O’Sullivan & Gabuzda, 2009) to obtain $\lg(B_{1}^{\rm eq}[{\rm G}])=$ -1.63 and -1.43 for $\beta=$ 0.7 and 0.96. For value of the core shift between 15.4 and 8.1 GHz the magnetic field is (-1.25, -1.05) for corresponding speeds. As discussed in O’Sullivan & Gabuzda (2009) using any $\alpha$ other than 0.5 in derivation of (Hirotani, 2005) results in $B_{1}^{\rm eq}$ dependence on $\gamma_{\rm min}$. Accounting for the estimated spectral index lowers $B_{1}^{\rm eq}$ with a factor $\approx$2 for values of $\gamma_{\rm min}$ from 10 to 100.

5.3 Plasma and jet magnetization

The dependence of the energy densities ratio of magnetic field and emitting particles on the jet Lorenz factor $\Gamma$ is presented in Figure 9. The vertical red stripe represents the spread of the component speeds observed within VLBI structure by Cotton et al. (1999). With the sampled posterior distribution of our model parameters (Figure 2) we also show the values of $\Gamma$ obtained using relation between jet opening angle and bulk motion Lorentz factor in Clausen-Brown et al. (2013) (shown as vertical green stripe around $\Gamma\approx 1.9$).

Our fit reveals that emitting plasma has $\beta_{\rm p}\approx 1$ for low energy cutoff $\gamma_{\rm min}=100$, where the plasma parameter $\beta_{\rm p}$ is defined through the energy densities (Zdziarski, 2014):

[TABLE]

where $u_{\rm e}$ and $u_{B}$ — emitting particle and magnetic field energy densities, $n_{\rm pl}$ — number density of electrons in a power-law at a distance $r_{1}=1$ pc from the apex, $B$ — magnetic field at this point (both in a plasma frame), $c$ — speed of light, $k_{\rm u}$=0 — contribution from particles outside of the power-law and ions. For $\gamma_{\rm min}=10$ we obtain $\beta_{\rm p}\approx 10$ with 2$\sigma$ uncertainty 0.5 dex. If the break in particles energy spectrum exists $\gamma_{\rm br}\leq 100$, $\beta_{\rm p}$ becomes lower and consistent with the equipartition.

Although we obtained $\beta_{\rm p}\approx 1$ for $\gamma_{\rm min}=100$ it does not imply that jet magnetization parameter $\sigma\approx 1$, where $\sigma$ is the ratio of the Poynting flux to the kinetic energy flux in the black hole frame (Lyutikov et al., 2005; Nokhrina, 2017):

[TABLE]

where $n_{\rm p}$ is the number density (in jet frame) of protons at $r_{1}$, ($n_{\rm p}=0$ for pairs jet), $m_{\rm p}$ - proton mass. We also assume that number density of electrons in a power-law $n_{\rm pl}$ is some fraction ($1/f_{\rm pl}$) of all electrons $n_{\rm e}$ and in case of normal plasma jet $n_{\rm e}=n_{\rm p}$. Figure 10 shows the dependence of the jet magnetization $\sigma$ on the lower cutoff in the power law distribution of the emitting particles $\gamma_{\rm min}$ for $e^{-}/e^{+}$– and $e^{-}/p$–jet. It is expected that acceleration and collimation of initially highly magnetized ($\sigma\gg 1$) jet takes place till the equipartition $\sigma\propto 1$ (Potter & Cotter, 2013, and references therein). Our results are consistent with this only for the pairs plasma jet and $\gamma_{\rm min}\leq 100$. Normal plasma jet is particle-dominated up to $\gamma_{\rm min}\propto 1000$.

Soft spectral index obtained in Section 5.1 could imply the magnetic reconnection as a particle acceleration process (Sironi & Spitkovsky, 2014). Magnetic reconnection in a both normal and pairs plasma results in various electrons spectra – from $s\approx 1.0$ to to $s\approx 4$ depending on magnetization $\sigma$ and value of the guide magnetic field (Sironi & Spitkovsky, 2014; Guo et al., 2015; Guo et al., 2016; Werner & Uzdensky, 2017; Werner et al., 2018). It however requires $\sigma\gg 1$ to be efficient (Werner et al., 2018). Obtained values (Figure 10) are thus consistent with magnetic reconnection only for $e^{-}/e^{+}$ jet and $\gamma_{\rm min}\geq 100$.

5.4 Jet composition

To constrain the jet composition we plot the dependence of the jet power $P_{\rm jet}$ estimated from our fit on $\gamma_{\rm min}$ for both normal and pairs plasma jet in Figure 11. Here, following Zdziarski (2014) the jet power is presented as sum of components – power in electrons/positrons, magnetic field and bulk-motion kinetic power:

[TABLE]

where $\eta_{B}$ - magnetic adiabatic index (4/3 for tangled field), $\eta$ - average adiabatic index (4/3 < $\eta$ < 5/3). Lines represent median of posterior distribution and shaded regions – corresponding $1\sigma$-interval. Continuous lines show $f_{\rm pl}=10$ and dashed – $f_{\rm pl}=100$.

Bicknell (1994) estimated the jet power of NGC 315 using extended radio structure from $4\cdot 10^{42}$ up to $10^{43}$ erg/s assuming an equal contribution of the lobe internal energy and the work done on the lobe expansion. De Young (2006) shown that the ratio of the work done on the lobe expansion to the lobe internal energy could be higher (with typical ratio values equal 3–10). Cavagnolo et al. (2010) estimated $\lg{P_{\rm jet}}=6.58(+2.48-4.96)\cdot 10^{42}$ erg/s using cavities in a X-ray emitting gas. As the radio lobes of NGC 315 are poorly confined, they assumed that volume of a cavity equals the volume of the corresponding radio lobe and extrapolated pressure profiles beyond the region observed in X-rays. Morganti et al. (2009) used correlation between core radio luminosity and jet power from (Heinz et al., 2007) and obtained $P_{\rm jet}=10^{44.15}$ erg/s with uncertainty 0.5 dex. However, this estimate was obtained for a sample of sources that are biased against beaming (Heinz et al., 2007) and probably is the upper limit for NGC 315.

Using expression for Lorentz factor $\gamma_{\rm rad}$ of particles emitting at synchrotron-self absorption peak $\nu$ ($\gamma_{\rm rad}\approx\sqrt{\nu/\nu_{B}}$, where the observed frequency $\nu$ is corrected for the Doppler factor and $\nu_{B}$ is Larmor frequency for given magnetic field) we can constrain electrons low-energy cutoff as $\gamma_{\rm min}<\gamma_{\rm rad}$. For flat spectrum sources $\nu/\nu_{B}$ is nearly constant and emission is dominated by electrons with the same $\gamma_{\rm rad}$ at all radii (Beckert & Falcke, 2002; Bjornsson, 2019). In our model $\gamma_{\rm rad}$ = 99 for 15.4 GHz core. Thus $\gamma_{\rm min}$ can not be larger than $\approx$99. This corresponds to the vertical red stripe in Figure 11 that shows 1$\sigma$ interval obtained from the posterior (Figure 2).

(Celotti & Fabian, 1993; Reynolds et al., 1996) found that if for $e^{-}/p$–jet lower energy cutoff $\gamma_{\rm min}<100$ than such jet carries more energy than it seems to be dissipated. For pairs plasma jet $\gamma_{\rm min}$ should be $\approx 1$ for jets to carry the amount of the energy that is dissipated. From Figure 11 it is apparent that electron-proton jet fails to explain both the necessary energetic and the lower energy cutoff in a power-law. Note, that possible break in a power-law at $\gamma_{\rm br}$ does not influence the jet power estimate for those $\gamma_{\rm min}$ that are close to the break. In case we are observing the steepen part of the synchrotron spectrum the break $\gamma_{\rm br}$ can not be larger than $\gamma_{\rm rad}\approx 100$ for core region at 15.4 GHz that is well described by the non-uniform model with optically thin steep spectral index $\alpha\approx 1.0$. However $\gamma_{\rm min}$ should not be much less or jet will carry the excessive power. In other words the possible break does not influence the conclusion concerning excess energetic of a normal plasma jet.

If estimate of $P_{\rm jet}$ from Cavagnolo et al. (2010) is underestimated, i.e. value from Morganti et al. (2009) is applicable, possibly with the corresponding correction for beaming (see also Pjanka et al., 2016; Inoue et al., 2017, for discussion of different methods of $P_{\rm jet}$ estimation) than electron-proton jet can be reconciled with this jet power only if $f_{\rm pl}$ is not larger than $\propto$10. This implies quite efficient particles acceleration that is higher than expected for weakly magnetized relativistic shocks ($\propto 1\%$, Sironi et al., 2013).

Using X-ray observation of NGC 315 we can check self-consistency of our model. As the number of particles depends on $s$ and $\gamma_{\rm min}$, with high $s$ and low $\gamma_{\rm min}$ we can get too many particles, which will produce excessive X-rays via Synchrotron-Self Compton (SSC) mechanism (Reynolds et al., 1996). For particles density and size of the emitting region corresponding to the 15.4 GHz core we obtain 1 keV X-ray flux $1.2\cdot 10^{-7}$ Jy assuming high-energy cutoff in the synchrotron spectrum $\nu_{\rm b}=\exp(10)\cdot\nu_{\rm 15.4GHz}$ (Ghisellini et al., 1992; Reynolds et al., 1996). Worrall et al. (2007) used Chandra X-ray observations and obtained flux of the power law component of the central core region $(1.2\pm 0.2)\cdot 10^{-7}$ Jy which they attributed to a jet.

6 Components besides inhomogeneous model

6.1 Counter-jet component

As we noted in Section 3 the analysis of the residuals between the data and best fitted inhomogeneous model revealed the presence of the component in a counter-jet. Parameters of this component are presented in Table 2. Mean values of the posterior distributions are presented as well as corresponding 1$\sigma$ credible intervals.

The spectrum of the counter-jet component and its fit by an optically thin synchrotron model are presented in Figure 12. Optically thin spectral index is $0.73\pm 0.09$ and agrees within 2$\sigma$ with the spectral index of our jet model. Note that for a counter-jet the same observed frequencies correspond to intrinsic frequencies $D^{2}$ times ($\approx$3–4 times) higher than intrinsic frequencies for the approaching jet. This component can hardly be explained by a local brightening in a counter jet as the approaching jet is completely featureless. Thus it should be the position of a $\tau\approx 1$ surface, where optical depth $\tau$ is determined by intrinsic and(or) external absorption.

The position of this component relative to the cone apex of the jet changes little with frequency (see Table 2) and is consistent with being constant. This is shown in Figure 13 together with the same dependence for VLBI-core of the approaching jet according to our model. The intensity peak of the model emission was used as the VLBI core. As shown in Appendix B to find the core position one has to solve an exponential equation for parameters $n$, $m$ and $\alpha$, which gives $\tau\approx 3.1$ for maximum brightness. Then using our re-parametrization the corresponding dependence of the position of true core at some frequency $\nu_{\rm obs}$ is calculated as follows:

[TABLE]

where $k_{r}=(n+m(1.5+\alpha)-1)/(2.5+\alpha)$. From the posterior (Figure 2) we obtain the median $k_{r}=0.89$ and 95$\%$ credible interval for – $(0.86,0.91)$.

We also plot the expected position of the counter-core for a range of jet to counter-jet Doppler factor ratio. Here we assumed that synchrotron self-absorption is the only source of opacity. In Figure 14 green line shows the dependence of jet to counter-jet Doppler factor ratio on the jet Lorenz factor $\Gamma$ for viewing angle $\theta=38\pm 2\deg$ obtained in (Canvin et al., 2005). The red wide vertical stripe indicates range of the jet speeds 0.7–0.96 observed in pc-scale radio structure in Cotton et al. (1999). Purple vertical lines show values of $\Gamma$ obtained from the joint posterior of $\phi_{\rm app}$, $\theta$ from (Canvin et al., 2005) and relation between the jet opening angle $\phi$ and bulk motion Lorenz factor $\Gamma$ obtained in (Clausen-Brown et al., 2013).

As shown in (Lobanov, 1998)666In Section 2.4 of that paper $r_{\rm core}\approx\nu^{-1/2.5}$ should be instead of $\nu^{-2.5}$. See also (Kadler et al., 2004). in the environment with a steep density gradients $k_{r}$ as high as 2.5 are possible. In the presence of strong external density gradients values of $k_{r}$ are generally even higher. In our case position of the counter-jet component is nearly constant within the errors and requires $k_{r}\gg 1$. We thus attribute the counter-jet component to “core” for which the dominating opacity is due to the external absorber with extremely high ionized gas density gradient and size $\approx 0.1$ pc. This is close to the size of the Broad Line Region (BLR). Barth et al. (1999) detected broad polarized H$\alpha$ and H$\beta$ lines in optical spectrum of NGC 315 suggesting an obscured BLR to exists in its nucleus.

The core of counter-jet and its frequency dependent shift were detected before in e.g. Haga et al. (2013) in radio galaxy NGC 4261 using multi-frequency (7 bands from 1.4 to 43 GHz) VLBA data. Authors concluded that neither pure synchrotron self-absorption (SSA) nor pure external (free-free) absorption (FFA) can explain the shift. Kadler et al. (2004) obtained for NGC 1051 the value of $k_{r}$ as high as 6.8 for counter jet and from 2.1 to 4.1 for the approaching jet. They attributed this to dominating contribution of FFA to the opacity on the counter-jet side. Haga et al. (2013) found that lower and higher frequency position of the counter-jet core are explained by SSA and its position at medium frequency band needs contribution from the external absorption. They invoke model of outer disk and inner radiation inefficient accretion flow (RIAF) to explain such behaviour. In our data with only 3 well separated relatively high-frequency bands the counter-jet core is even not apparent in the reconstructed CLEAN-images and in simplistic Gaussians models due to its weakness and closeness to the bright core of approaching jet. Thus it can be recovered only through the detailed modelling of the source structure.

6.2 Constrain on plasma velocity

We can constrain the ratio of the Doppler factors $R$ for approaching jet and counter-jet by the position of the counter-jet component at the lowest frequencies. Indeed, visible component can not be closer to jet apex than the region with optical depth $\tau\sim 1$. The depth is determined by the opacity – both internal and external. In Figure 13 the orange stripe corresponds to $R\in(2.5,7.0)$ or equivalently $\beta\in(0.7,0.96)$ as estimated by Canvin et al. (2005). It turns out that independently from these estimates we can constrain $R>2.5$, i.e. $\beta>0.7$ at $r\approx$0.1 pc for viewing angle $\theta=38\deg$. We stress that this constrains the plasma velocity, while VLBI kinematics could measure the pattern speed velocity (e.g. Cohen et al., 2014). Note that possible parabolic shape of the jet at its origin only strengthens the constrain. In this case the apex of the jet is closer to the counter-jet side than it follows from the conical geometry of our model (Kovalev+ in prep). Thus the observed counter-jet core is also closer to the jet apex. This implies that the ratio of the Doppler factors and corresponding lower limit on the plasma velocity should be even higher.

6.3 Component in the approaching jet

The residuals between data and model imply the existence of a component at distance $\approx 4$ mas from core at all observed frequency bands. The parameters of this component obtained with Difmap modelling of the residual visibilities are listed in Table 3. Component shows optically thin spectrum with $\alpha_{8.1}^{12.1}$ = $0.55\pm 0.20$ and $\alpha_{12.1}^{15.4}$ = $4.0\pm 0.8$. The first agrees with our model spectral index at $2\sigma$ level, but the second is significantly steeper. Velocity $\beta=0.7$ corresponds to angular apparent speed $\mu_{\rm app}$ = 0.85 mas/year for the adopted viewing angle 38∘. Thus, the distance of 4 mas corresponds to travelling time from the core 4.7 years for $\beta=0.7$ and 3.6 years for $\beta=0.9$ assuming constant propagation speed. There are no signs of flaring activity near epoch 2001 in Figure 1, that is for 5 years preceding the observations we use. However the light curve cadence is too low to exclude flares with duration less than a year.

7 Conclusions

For the first time we fitted a physical jet model directly to the interferometric observables — visibilities. We used a non-uniform conical jet model which was generalized from Blandford & Königl (1979) by allowing arbitrary spectral index and exponents in power-law dependencies of the magnetic field and particle density on the distance. To break the degeneracies in model parameters we used multi-frequency VLBA data of a radio galaxy NGC 315 and assumed a constant magnetic field to particle energy density ratio. The observed data at all frequencies is well described by the model. We found that electron-positron jet is consistent with an independently derived jet power. The emitting plasma is consistent with the equipartition between magnetic field and emitting particles. The jet magnetization $\sigma$ is 0.1–1 depending on the fraction of pairs in a power-law. Electron-proton jet could be reconciled with the data only if the jet power inferred from cavities in a X-ray emitting gas is underestimated for NGC 315 and particles acceleration is highly efficient. Such jet is particle dominated with $\sigma\approx 10^{-3}$–$10^{-2}$.

We found a weak component on a counter-jet side in the residuals between the observed data and best model at all frequency bands. We attribute it to the external absorber with steep density gradient, possibly a border region of the BLR with a significantly ionized plasma. Position of the component at the 15.4 GHz constrains the plasma flow velocity $\beta>0.7$ at distance $r\approx 0.1$ pc from the central engine.

Acknowledgements

We thank the anonymous referee for helpful comments and suggestions. We also thank Marina Butuzova and Mikhail Lisakov for reading the manuscript and helpful comments. This work is supported by Russian Science Foundation grant 16-12-10481. This research has made use of data from the MOJAVE database that is maintained by the MOJAVE team Lister et al. (2009b).

This research has made use of NASA’s Astrophysics Data System.

This research made use of Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration et al., 2013), Numpy (Van Der Walt et al., 2011), Scipy (Jones et al., 2001), numba (Lam et al., 2015), FINUFFT (Barnett et al., 2018), emcee (Foreman-Mackey et al., 2013), PolyChord (Handley et al., 2015a, b). Matplotlib Python package (Hunter, 2007) was used for generating all plots in this paper.

Appendix A Deriving original formulas

Here we derive the essential formula from the original paper (Blandford & Königl, 1979) for completeness. Coefficients of emission and absorption for for synchrotron radiation in the observer frame expressed in terms of the plasma rest frame values of $B$ and $K$ in case of tangled magnetic field, power law energy distribution of the emitting particles and isotropic pitch-angles distribution (Rybicki & Lightman, 1979):

[TABLE]

[TABLE]

where $z$ is the redshift, $\alpha$ is optically thin spectral index related to the emitting particles powe-law exponent $\alpha=(s-1)/2$, $D$ – Doppler factor, $\langle C_{1}(\alpha)\rangle$ and $\langle C_{2}(\alpha)\rangle$ are averaged over isotropic pitch-angle distributions777For tangled magnetic field one has to integrate $C_{1}(\alpha)(\sin{\theta})^{\alpha+1}$ and $C_{2}(\alpha)(\sin{\theta})^{\alpha+1.5}$ over all directions values of the coefficients, expressed through Gamma functions as follows:

[TABLE]

[TABLE]

where

[TABLE]

[TABLE]

Using (11) and expressing geometrical length $L$ along line of sight at a map point parametrized by the observed distances along the jet $r_{\rm obs}=r\sin{\theta}$ and in perpendicular direction $d$ as:

[TABLE]

where $\phi$ is half-opening angle of the cone and $\theta$ – jet viewing angle, one can derive the expression for optical depth $\tau=\int\kappa dl$. The last can be approximated as $\tau=\kappa L$, where $\kappa$ is the value of the absorption coefficient at the apex of the jet for a given line of sight:

[TABLE]

Here $B$, $K$ and $D$ are values at point $r=r_{\rm obs}/\sin{\theta}$, $C_{1}(\alpha)$ and $C_{2}(\alpha)$ are already averaged over isotropic direction of the tangled magnetic field coefficients, i.e. it is (12) and (13). Substituting values of fields $B$ and $K$ at point $r$ (i.e. $B=B_{1}(r/r_{1})^{-1}$, where $r_{1}$ is 1 pc, and $K=K_{1}(r/r_{1})^{-2}$) and assuming $\alpha=0.5$, $m=1$, $n=2$:

[TABLE]

where $r_{1}$ is 1 pc. Substituting $r_{\rm ob}$ instead of $r$:

[TABLE]

This coincides with (26) from Blandford & Königl (1979). Using cgs values of the constants888These are already averaged values, i.e. it is (12) and (13) $C_{1}(0.5)=3.6\cdot 10^{-19}$ and $C_{2}(0.5)=9.24\cdot 10^{10}$ we obtain coefficient in (19) $\approx 570.4$. That is close to the original 500 in equation (26) in (Blandford & Königl, 1979). To completely move from true half-opening angle and distance from cone apex ($\phi$, $r$) to the observed values ($\phi_{\rm app}$, $r_{\rm obs}$) getting $\sin{\theta}$ out of the brackets gives:

[TABLE]

For brightness distribution using the solution of the radiative transfer equation:

[TABLE]

expressions for coefficients (10)–(11) and Lorentz invariant $\frac{I_{\nu}}{\nu^{3}}$ we obtain:

[TABLE]

To compare with Blandford & Königl (1979) results we express the brightness temperature:

[TABLE]

Substituting the magnetic field coordinate dependence and using $r_{\rm obs}$ instead of $r$ we obtain:

[TABLE]

where $r_{1}$ is 1 pc.

Using cgs values of the constants $C_{1}(0.5)=3.6\cdot 10^{-19}$ and $C_{2}(0.5)=9.24\cdot 10^{10}$ we obtain coefficient $\approx 3.186\cdot 10^{10}$K that coincides with that from equation (25) in Blandford & Königl (1979).

Appendix B Extending for arbitrary $\alpha$, $m$, $n$

Generalization to arbitrary $m$, $n$, $\alpha$ is obvious. One should use $(r/r_{1})^{-m}$ instead of $(r/r_{1})^{-1}$ (similar for $K$) when going from (17) to (18) and from (23) to (24) and keep $\alpha$ after (17). For the optical depth we obtain:

[TABLE]

For canonical $n=2$, $m=1$, $\alpha=0.5$ we obtain $\tau(r_{\rm obs},d)\propto r^{-3}_{\rm obs}$ in accordance with (20). The numerical coefficient strongly depends on $\alpha$ through $C_{2}(\alpha)$. E.g. for $\alpha=0.55$ it is $\approx 1237$ and for $\alpha=0.75$ it is $\approx 27980$.

For brightness temperature:

[TABLE]

where $\tau(r_{\rm obs},d)$ is from (25) and numerical coefficient depends only on $\alpha$. E.g. for $\alpha=0.75$ it is $\approx 2.37\cdot 10^{10}$.

To find position of the core (i.e. maximum intensity) at some frequency one has to differentiate Eq. 26 with respect to $r_{\rm obs}$. This way the following relation for the optical depth at the maximum intensity is obtained:

[TABLE]

For $n=2$, $m=1$ we obtain $\tau\approx 2.92$. Thus to obtain position of maximum $r_{\rm max}$ for some frequency one has to solve equation Eq. 25 for $r_{\rm obs}$ given $\tau=2.92$,

Appendix C “Natural” parametrization

With (25) and (26) the essential parametrization:

[TABLE]

[TABLE]

With this parametrization:

[TABLE]

[TABLE]

where the observed frequency $\nu_{\rm obs}$ is in GHz and projected distance from cone apex $r_{\rm obs}$ is in pc and $C_{T}(\alpha$) and $C_{\tau}(\alpha)$ are:

[TABLE]

[TABLE]

expressed in cgs units with values that depend on $\alpha$ only. For example, $C_{\tau}(0.5)\approx 570.4$ in accordance with that from (19) and $C_{T}(0.5)\approx 3.186\cdot 10^{10}$ in accordance with formula (26) and (25) from (Blandford & Königl, 1979).

Recall that for intrinsically symmetrical counter jet the parameters of are related to that of the approaching jet:

[TABLE]

[TABLE]

where $R=(1+\beta\cos{\theta})/(1-\beta\cos{\theta})$ is the ratio of the Doppler factors of the approaching and receding jets.

Physical parameters can be obtained in “natural” parametrization as:

[TABLE]

[TABLE]

[TABLE]

Appendix D Accounting for the uncertainty in the amplitude scale factors

After sampling of the posterior distribution with the constrained ($\prod_{b}g_{b}=1$) amplitude scale factors $g_{b}$ we can infer the uncertainty resulting from their product not being exactly 1. As we noted in Section 3.1, multiplication of all scale factors by the same number is equivalent to multipliction of $A_{1}$ on this number. Thus, corresponding uncertainty of the $A_{1}$ is the uncertainty of the mean of the scale factors. We used priors $N(1,\sigma_{b})$ for $g_{b}$ at each frequency band, where $\sigma_{b}$ – absolute calibration uncertainty estimated by Hovatta et al. (2014) and samples from the obtained posteriors of $A_{1}$, $A_{2}$:

[TABLE]

where $\sigma_{\rm gains}/n_{b}$ – uncertainty of the mean of scale factors deduced from their priors and $n_{b}$ – number of frequency bands.