A bias in VLBI measurements of the core shift effect in AGN jets
MMNRAS , 1–11 (2020) Preprint 23 October 2020 Compiled using MNRAS L A TEX style file v3.0
A bias in VLBI measurements of the core shift effect in AGN jets
I. N. Pashchenko, ★ A. V. Plavin, , A. M. Kutkin, , Y. Y. Kovalev , , Astro Space Center, Lebedev Physical Institute, Profsouznaya 84/32, Moscow 117997, Russia Moscow Institute of Physics and Technology, Institutsky per. 9, Moscow region, Dolgoprudny, 141700, Russia ASTRON, The Netherlands Institute for Radio Astronomy, Oude Hoogeveensedijk 4, 7991 PD, Dwingeloo, The Netherlands Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, D-53121 Bonn, Germany
Accepted 2020 October 7. Received 2020 October 7; in original form 2020 July 27
ABSTRACT
The Blandford and Königl model of AGN jets predicts that the position of the apparent opaque jet base — the core — changeswith frequency. This effect is observed with radio interferometry and is widely used to infer parameters and structure of theinnermost jet regions. The position of the radio core is typically estimated by fitting a Gaussian template to the interferometricvisibilities. This results in a model approximation error, i.e. a bias that can be detected and evaluated through simulations ofobservations with a realistic jet model. To assess the bias, we construct an artificial sample of sources based on the AGN jetmodel evaluated on a grid of the parameters derived from a real VLBI flux-density-limited sample and create simulated VLBIdata sets at 2.3, 8.1 and 15.4 GHz. We found that the core position shifts from the true jet apex are generally overestimated. Thebias is typically comparable to the core shift random error and can reach a factor of two for jets with large apparent openingangles. This observational bias depends mostly on the ratio between the true core shift and the image resolution. This impliesthat the magnetic field, the core radial distance and the jet speed inferred from the core shift measurements are overestimated.We present a method to account for the bias.
Key words: techniques: interferometric – galaxies: active – galaxies: nuclei – galaxies: jets
The Very Long Baseline Interferometry (VLBI) technique allows ex-ploring Active Galactic Nuclei (AGN) jets with the highest angularresolution. Synthesised images of most of the radio-loud AGN ob-served with VLBI are dominated by an unresolved or barely-resolvedfeature called the core (or VLBI core, e.g. Zensus 1997; Kovalev et al.2005). The core has a flat or inverted radio spectrum, characteris-tic of optically-thick synchrotron emission (Hovatta et al. 2014).Its frequency-dependent position was discovered by Marcaide &Shapiro (1984) and can be explained within the relativistic jet model(Blandford & Königl 1979; Konigl 1981). The core appears in thatmodel as a surface in a continuous flow where the optical depth ofjets synchrotron radiation at a given frequency 𝜏 𝜈 ≈ ★ E-mail:[email protected] (e.g. Hada et al. 2011; Algaba et al. 2017). They were also used toconstrain the spins of the supermassive black holes (Agarwal et al.2017) and the nature of the central engine (Zamaninasab et al. 2014).If a source demonstrates flaring activity, it is possible to derive thejets innermost regions kinematics based on the measured core shiftand multi-frequency time delays (Kudryavtseva et al. 2011; Kutkinet al. 2014, 2018, 2019; Plavin et al. 2019a). Differences in corepositions should be taken into account for image registration whenconstructing the VLBI spectral index (e.g. Marr et al. 2001; Kovalevet al. 2008; Hovatta et al. 2014) and Faraday rotation maps (e.g.Hovatta et al. 2012). The core shift effect is also important for theradio (VLBI) and optical (Gaia) reference frame alignment (Porcas2009; Kovalev et al. 2017; Plavin et al. 2019b), as it can partiallyexplain the observed offsets.However, conventionally the position of the core and, thus, thecore shift is estimated by fitting a set of Gaussian templates to theinterferometric visibility. This could result in a biased estimates ofthe core shift if the emission profile is non Gaussian. Plavin et al.(2019a) assumed a linear dependence between the measured coreshift and the synthesised beam for a multi-epoch multi-frequency(2.3 and 8.4 GHz) sample of 40 sources and found that the measuredcore shift is positively correlated with resolution, i.e. the beam size.Their results could imply that core shifts are typically overestimated.However, Sokolovsky et al. (2011) did not find larger core shiftsbetween 1.4 and 15.4 GHz in the direction of the elongated beamin their sample of 20 sources. In this paper we study the possiblebias directly, using the jet model (Blandford & Königl 1979; Konigl1981) that is used to explain the core shift effect. © a r X i v : . [ a s t r o - ph . GA ] O c t I. N. Pashchenko et al.
Throughout the paper, we adopt the Λ CDM cosmology with Ω 𝑚 = . Ω Λ = . 𝐻 = . − Mpc − (Hinshaw et al.2013). The spectral index is defined via the dependence of the fluxdensity on the spectral index 𝑆 ∝ 𝜈 − 𝛼 . The jet model by Blandford & Königl (1979); Konigl (1981) assumesa conical geometry of the outflow and the particle density amplitudeand the magnetic field dependence on the distance 𝑟 from the jetapex: 𝐾 = 𝐾 ( 𝑟 / 𝑟 ) − 𝑛 and 𝐵 = 𝐵 ( 𝑟 / 𝑟 ) 𝑚 , where the referenceposition is 𝑟 = 𝐵 and 𝐾 are the values at the 𝑟 . In thegeneral case of arbitrary exponents 𝑚 , 𝑛 and optically thin spectralindex 𝛼 analytical expressions for the distribution of optical depthand brightness temperature can be found in, e.g. Pashchenko & Plavin(2019).Lobanov (1998) has shown that the relativistic jet model can beused to estimate the jets physical parameters, e.g. magnetic fieldand particle density, using the measured core shift. However, thisrequires some further assumptions because core shift measurementsalone do not constrain the parameters of the jet model. More on de-generacies between model parameters can be found in Pashchenko& Plavin (2019). The most widely used is the assumption of anequipartition between emitting particles and magnetic field energydensities (Lobanov 1998; O’Sullivan & Gabuzda 2009; Pushkarevet al. 2012). Zdziarski et al. (2015) used the optically thick approxi-mation to avoid the equipartition assumption (for a similar approach,see also Nokhrina 2017). They obtained magnetic field estimateswhich significantly differ from the corresponding equipartition val-ues. Both approaches needed measurements of the core shifts be-tween two frequencies (cid:52) 𝑟 . In the equipartition case, the magneticfield scales with the core shift as 𝐵 eq1 ∝ (cid:52) 𝑟 . 𝑘 r (Lobanov 1998),where 𝑘 𝑟 = 𝑛 + 𝑚 ( . + 𝛼 )− . + 𝛼 . However, for a more general case, the de-pendence on the measured core shift is much stronger: 𝐵 noeq1 ∝ (cid:52) 𝑟 (Zdziarski et al. 2015).There are several methods to measure the core shift effect (for abrief review, see Kutkin et al. 2014). In general, one has to knowthe relative position of the core at several frequencies. The absoluteposition information is lost during self-calibration of the VLBI data(Pearson & Readhead 1984). Thus in the absence of a phase cali-brator, the first step in deriving the core shift is matching imagesat different frequencies using a reference point that is stable withfrequency. The optically thin structure of a jet is the most obviousreference for core shift measurements. Depending on the morphol-ogy of the extended jet structure, two methods are conventionallyused. In the first one, the core position is measured with respect to anindividual jet feature using a brightness distribution model (Lobanov1998; Kovalev et al. 2008; Sokolovsky et al. 2011; Kutkin et al.2014). The second measure core position with respect to a large jetsection rich in structure using image cross-correlation (Walker et al.2000; Croke & Gabuzda 2008; O’Sullivan & Gabuzda 2009; Algabaet al. 2012; Hovatta et al. 2012; Pushkarev et al. 2012; Fromm et al.2013; Plavin et al. 2019a).Regardless of the method used to align images, the final step inderiving the core shift is the measurement of the core position it-self. Conventionally, the source structure is modelled with a set ofGaussian templates (e.g. Lister et al. 2009b). More recently, Plavinet al. (2019a) employed subtraction of the extended structure fromthe original data set and considered the resulting data as a contribu-tion from the VLBI core only. In both approaches, the core brightness distribution is fitted with a simple model, i.e. usually a circular Gaus-sian. This could result in a bias of the measured core position if thetrue brightness distribution is significantly non-Gaussian. Here wedefine the bias as the difference between the observed and true values.However, before exploring the bias of the measurement, it is worthfinding out what is supposed to be measured and what is actuallymeasured. 𝜏 = REGION BEING OFFSET FROMTHE INTENSITY PEAK
To find the analytical expression for the position of the maximumintensity at some frequency, one has to differentiate the equation ofthe brightness temperature radial distribution. The brightness tem-perature distribution in a general case of arbitrary 𝑛 , 𝑚 and 𝛼 is(Pashchenko & Plavin 2019): 𝑇 ( 𝑟 obs , 𝑑 ) = 𝑐 𝐶 ( 𝛼 ) 𝐷 . 𝜈 . 𝑟 . 𝑚 obs 𝜋𝑘 𝐵 𝐶 ( 𝛼 )( + 𝑧 ) . 𝐵 . 𝑟 . 𝑚 ( sin 𝜃 ) . 𝑚 ( − 𝑒 − 𝜏 ( 𝑟 obs ,𝑑 ) ) , (1)where the optical depth is 𝜏 ( 𝑟 obs , 𝑑 ) = 𝐶 ( 𝛼 )( + 𝑧 ) −( 𝛼 + . ) 𝐷 . + 𝛼 𝐾 𝐵 . + 𝛼 𝜈 −( . + 𝛼 ) obs × 𝑟 𝑛 + 𝑚 ( . + 𝛼 ) ( sin 𝜃 ) 𝑛 + 𝑚 ( . + 𝛼 )− √︄ 𝜙 − (cid:18) 𝑑𝑟 obs (cid:19) 𝑟 −( 𝑛 + 𝑚 ( . + 𝛼 )− ) obs . (2)Here the position in the sky is described by 𝑟 obs and 𝑑 — distancesalong and perpendicular to the jet direction, 𝐷 — Doppler factor, 𝜙 app is the observed jet half-opening angle, 𝑧 — redshift, 𝐶 ( 𝛼 ) and 𝐶 ( 𝛼 ) — functions of optically thin spectral index 𝛼 , 𝑐 — speed oflight and 𝑘 B — Boltzmann constant.To this end, the following relation for the optical depth at themaximum intensity is obtained: 𝜏 𝑟 max = log (cid:18) + 𝑛 + 𝑚 ( . + 𝛼 ) − . 𝑚 𝜏 𝑟 max (cid:19) . (3)For 𝑛 = 𝑚 = 𝛼 = . 𝜏 𝑟 max ≈ .
92. This stays closeto 𝜏 ≈ 𝑛, 𝑚 and spectral index 𝛼 . Thus to obtain the position of the maximum 𝑟 max for some frequency, one has to solve the equation of the opticaldepth radial distribution for 𝜏 = . 𝜏 =
1. These twoapproaches are not equivalent, and combining equations from papersadopting different conventions can become confusing and lead tosubtle biases.Both definitions of the core consider the on-axis distributions. Anexample of the on-axis intensity profile of the Blandford-Königl jetis given in Fig. 1. The ratio of the distances from the jet apex tothe 𝜏 = 𝑟 𝜏 = and to the intensity maximum 𝑟 max can becalculated by substituting 𝜏 𝑟 max into the equation for optical depthdistribution along the jet axis (Pashchenko & Plavin 2019): 𝑟 max 𝑟 𝜏 = = 𝜏 − /( 𝑛 + 𝑚 ( . + 𝛼 )− ) 𝑟 max . (4)For a wide range of parameters, this ratio is 𝑟 max / 𝑟 𝜏 = ≈ .
7, and 𝑟 max differs from 𝑟 𝜏 = by ≈ MNRAS , 1–11 (2020) ias of the core shift in AGN jets Figure 1.
Radial slice of the emission pattern from a Blandford-Königl jetmodel for viewing angle 𝜃 = ◦ and half-opening jet angle 𝜙 = . ◦ showingthe difference between the positions of the VLBI core according to its differentdefinition as the region of 𝜏 = core positions measured by Gaussian model fitting were unbiased,then the magnetic field estimated from the core shifts would bebiased downwards by a factor ≈ . ≈ Blandford & Königl (1979) model has a significantly asymmetricradial profile (Fig. 1). To assess the bias due to its approximationwith a Gaussian template, simulations with a realistic brightnessdistribution model are necessary. Here we describe a study of thesystematics in core shift measurements using brightness distributionof the Blandford-Königl jet model, constructing artificial visibilitydata sets and processing them with conventional methods.In this section, we use the results of MOJAVE observations toconstruct and study the properties of our artificial sample. We notethat through the years MOJAVE group has constructed and publishedresults for slightly different samples, including the MOJAVE-I sam-ple (Lister et al. 2009a; Pushkarev et al. 2012; Clausen-Brown et al.2013), its extension in (Pushkarev et al. 2017), as well as the MOJAVE1.5 Jy Quarter-Century sample (1.5JyQC, Lister et al. 2019). Allof them are VLBI flux-density-limited samples of Doppler-boostedAGN jets. Their properties are comparable and do not differ signifi-cantly for our goals.
The brightness distribution of the Blandford & Königl (1979) modeldepends on several parameters. Only some combinations of theseparameters are identifiable on the basis of the VLBI observations(Pashchenko & Plavin 2019). Using a uniform or random grid inthe parameter space will result in sample parameter values, whichare not typical for the sources observed in flux-density-limited sam-ples (e.g. 1.5JyQC). Thus we construct a grid of model parameters conditioning them on the observed data in the following way. First,we draw a source with the parameters ( 𝐿, 𝑧, 𝜃, Γ ) and the flux den-sity > . 𝐿 — luminosity , 𝑧 —redshift, 𝜃 — viewing angle and Γ — bulk motion Lorentz factor.We draw 𝑧 from constant co-moving density between 𝑧 min = . 𝑧 max = .
4. Thus we obtained the power-law low-luminosity cut-offof the unbeamed luminosity function from Lister et al. (2019). Thenwe draw 𝐿 from the obtained power-law distributed luminosity func-tion and Γ from the power-law with limits Γ min = .
25 and Γ max = − .
4, obtained from the Monte-Carlo modelling of the1.5JyQC sample in Lister et al. (2019). The viewing angle 𝜃 wasdrawn from the isotropical (i.e. uniform in cos 𝜃 ) distribution. Fromthe unbeamed luminosity and Doppler factor 𝛿 ( 𝜃, Γ ) we calculatedthe observed luminosity 𝐿𝛿 𝑝 + 𝛼 , where we assume beaming expo-nent 𝑝 = 𝛼 =
0. If theobserved flux corresponding to the parameters ( 𝐿, 𝑧, 𝜃, Γ ) was largerthan the threshold value 1.5 Jy, the source was accepted to the sim-ulated sample. The distributions of redshifts and Lorentz factors ofthe simulated sources are presented in Fig. 2.As expected, the redshift distribution is close to that of the 1.5JyQCsample, while the distribution of Γ corresponds to its best fit Monte-Carlo simulation (see fig. 11 in Lister et al. 2019). To obtain thevalue of the half-opening angle 𝜙 / , we used relation 𝜙 / ≈ 𝜌 / Γ derived from statistical modelling of the opening angles distributionby Clausen-Brown et al. (2013), where 𝜌 is drawn from the normaldistribution 𝑁 ( 𝜇 = . , 𝜎 = . ) . The resulted distribution ofapparent opening angles 𝜙 app in the simulated sample is consistentwith those obtained by Pushkarev et al. (2017): compare our median 𝜙 app = . ◦ with their 𝜙 app = . ◦ .We assumed a conical jet geometry with a radial velocity fieldand exponents of the magnetic field and particle density amplitude 𝑚 = 𝑛 = Γ and isother-mal particle distribution. The second implies a dominant transversemagnetic field component. Together they result in a constant ratioof the emitting particles to the magnetic field energy densities. Thischoice is justified because most of the core shift measurements wereused in deriving the physical properties of the jet under the equipar-tition assumption. Note that the departures from the conical shapeare seen close to the jet base up to the distance of 10 gravitationalradii (Kovalev et al. 2020). This affects mostly nearby sources withthe given limited VLBI resolution. We employed synchrotron radia-tion transfer coefficients for a tangled magnetic field model (Longair1994) and assumed the optically thin spectral index 𝛼 = .
5, whichcorresponds to 𝑠 = . 𝑁 ( 𝛾 ) ∝ 𝛾 − 𝑠 .To limit the visible size of the jet, we included a steepening of thespectral index with the value 0.0001 pc − . With this settings, themodel spectrum is flat (Blandford & Königl 1979). We also assumedthat the low- and high-frequency cut-offs in the spectrum fall outsideof the chosen frequency range due to an external jet scale and breakin the energy distribution of emitting particles in the magnetic field(Konigl 1981).To obtain the sample of the core shifts, we generated Stokes I im-ages for frequencies 2.3, 8.1 and 15.4 GHz (S, X and U radio bands,respectively) numerically solving the radiative transfer equation us-ing an adaptive Dormand-Prince 5 type of RungeKutta algorithm(Dormand & Prince 1980). As we assume radial velocity field, theDoppler factor 𝛿 changes from the near to the far side of the jet.To obtain the synthetic VLBI data sets, we substituted the observedvisibilities of a single template data set with their predicted modelvalues and added a Gaussian noise estimated from the scatter of the MNRAS , 1–11 (2020)
I. N. Pashchenko et al.
Figure 2.
Left : distribution of redshifts for the simulated sample ( blue ) and sources with 𝑧 > . orange ). Right : distributionof the Lorentz factors for the simulated sample. observed visibilities. We used real VLBI observations of the source1458+718 from Lister et al. (2009b) and Sokolovsky et al. (2011)as a template for ( 𝑢, 𝑣 )-coverage and noise estimation. To estimatethe thermal noise from the observed visibilities, we employed thesuccessive differences approach (Briggs 1995). We choose a sourcewith a high declination to reduce possible effects due to an elongatedVLBA beam. We used major axes of naturally weighted CLEANbeams for calculations: 0.81, 1.61 and 6.51 mas for 15.4, 8.1 and 2.3GHz. As will follow in subsection 4.3, for the straight jet model itis the core shift in units of the beam size that determines the valueof the core shift bias. For elongated beams (i.e. sources with a lowdeclination), we expect qualitatively the same behaviour with the ef-fective resolution determined by the beam size along the core shiftdirection (subsection 4.4).Following the common practice, we modelled the resulting datasets in the ( 𝑢, 𝑣 )-plane at each frequency using simple Gaussian tem-plates in
Difmap (Shepherd 1997). We automated the modelling pro-cedure routine using several stopping and filtering criteria to choosethe “best” number of components, thus mimicking the modellingprocess by a human. At each step, the procedure analyses the resid-ual map and suggests new components estimating their parametersby the method of moments (Wasserman 2010). We also employed theapproach of Plavin et al. (2019a) of obtaining VLBI core parameters.It consists of estimating the extended emission by a list of CLEANcomponents outside of the core region and subtracting them from thedata set in the ( 𝑢 , 𝑣 )-plane. The obtained residuals are then modelledwith a single Gaussian component. We estimated the true core shiftvalue using the model brightness distribution and determined theobserved shift of the VLBI core in the corresponding Difmap modelfitted to the synthetic VLBI data. Note that our brightness distribu-tions are perfectly aligned, i.e. the jet cone apex is always positionedin the phase center of the VLBI images in the model.Finally our simulated sample includes 1000 sources. We filtered22 sources with core component flux densities larger than 10 Jy andless than 100 mJy. We also excluded 3 outliers with the ratio of theobserved to true core shifts between 8.1 and 15.4 GHz less than0.2 and one source with the observed core shift between 8.1 and15.4 GHz larger than 1 mas.
The median of the observed core shifts between 15.4 and 8.1 GHz forthe MOJAVE-I sample is 0.114 mas (see tab. 1 in Pushkarev et al.2012) for 𝑧 > 0.1. This agrees with the value of our simulated sampleof 0 . ± .
003 mas for a circular Gaussian core model (hereafter,CG). Hereafter the error corresponds to the 95% confidence interval.The confidence interval for the median was obtained by applying theBinomial distribution (Conover 1999). This agreement is expectedbecause we run our simulations assuming population properties in-ferred from the MOJAVE 1.5JyQC sample. The median core shift forthe elliptical Gaussian core model (hereafter, EG) is 0 . ± . . ± .
013 mas. This is close to the value 0.61mas extrapolated from the typical core shift between 15.4 and 8.1GHz, assuming that the core shift Δ 𝑟 between frequencies 𝜈 and 𝜈 scales with ( 𝜈 − 𝜈 )/( 𝜈 𝜈 ) (Lobanov 1998). Plavin et al. (2019b)estimated the typical core shift between 2.3 and 8.6 GHz as 0.53 masfor a sample of 40 sources observed at 10 epochs or more, mainlyduring geodetic VLBI sessions. However, a different method of coreparameters extraction was employed (subsection 4.1) there. In oursimulated sample, this method provides the median core shift 0.91mas. Note that both studies required specific criteria on the core shiftmeasurement to be reliable, e.g. the same source structure at bothbands (Kovalev et al. 2008), prominent extended structure (Plavinet al. 2019b), etc. These criteria are potentially capable of loweringobserved typical core shifts relative to core shifts estimated for oursimulated sample. Kovalev et al. (2008) estimated the theoreticalmean core shift from 0.13 to 0.30 mas for a flux-density-limitedsample between 2.3 and 8.6 GHz for the Blandford & Königl (1979)model. The specific value depends on the distribution of the Lorentzfactors and assumes a mean synchrotron luminosity 10 erg/s anda mean redshift 𝑧 =
1. We note that this is significantly lower thanin our simulated sample based on different underlying assumptionsabout the source population. This also could be the reason for a lowertypical core shifts.In Fig. 3 the observed and true core shifts obtained from the
MNRAS , 1–11 (2020) ias of the core shift in AGN jets Figure 3.
Distribution of the observed and true core shifts between 15.1 and 8.4 GHz ( left , 2 (true), 2 (CG), 8 (EG) points are not shown) and between 8.1 and2.3 GHz ( right , 1 (true), 2 (CG) and 14 (EG) points are not shown). Here CG means a circular Gaussian, and EG — an elliptical Gaussian core model. Theobserved core shifts are measured with Circular Gaussian (CG, blue colour) and Elliptical Gaussian (EG, orange colour) core model. simulated sample are compared. The positive bias of the core shiftestimates is evident — the measured values are systematically largerthan the true ones. The median of the true core shifts between 15.4and 8.1 GHz is 0.09 mas, while for 2.3 and 8.1 GHz it is 0.43 mas.Pushkarev et al. (2012) indicate a value of 0.05 mas as the randomerror of the core shift measurements between 15.4 and 8.1 GHz ob-tained using three approaches. This estimate includes a contributionfrom thermal noise, calibration error and mask size error. Plavin et al.(2019b) estimated the typical uncertainty of 8.6 - 2.3 GHz core shiftsas 0.2 mas. The median bias of the estimated core shift between 8.1and 15.4 GHz is 0.024 and 0.078 mas for circular and elliptical coremodels. For 8.6 - 2.3 GHz the corresponding median biases are 0.09and 0.37 mas. These values are comparable to the random errors.However, even more important is the possible bias dependence onthe source properties.
The core shift is determined by the difference in core positions attwo frequencies. The dependence of the ratio of the observed 𝑟 obs tothe true 𝑟 true projected core distance on the true core distance 𝑟 true expressed in units of the restoring beam FWHM for our simulatedsample is shown in Fig. 4 for 15.4, 8.1 and 2.3 GHz. Apparently,the core position is overestimated for the elliptical Gaussian coretemplate. However, using circular Gaussians for core fitting maylead to an underestimation of 𝑟 true . In our simulations, this occursat high jet viewing angles when the core is modelled with two closecomponents. The position of the first one (more compact with a flatterspectra, thus, treated as the core) lays upstream of the maximum ofthe model brightness.We found that the relative angular resolution is the main factorwhich determines the fractional core position bias. This dependencein Fig. 4 can be well described by: 𝑟 obs 𝑟 true = 𝐶 (cid:16) 𝑟 true FWHM (cid:17) − 𝑘 (5)with 𝐶 = . ± .
008 and 1 . ± .
013 for CG and EG. Exponent 𝑘 = . ± .
02 and 0 . ± .
01 for CG and EG. Note that themeasurements at all frequencies are shown in Fig. 4. Assuming thatthe true position of the core for a given source is 𝑟 true ( 𝜈 ) = 𝐶 𝜈 − / 𝑘 r and FWHM ( 𝜈 ) = 𝐶 𝜈 − , the observed core position of this sourceat some frequency 𝜈 yields: 𝑟 obs ( 𝜈 ) = 𝐶 𝐶 (cid:18) 𝐶 𝐶 (cid:19) − 𝑘 𝜈 − / 𝑘 r , eff (6)and the ratio of the observed and true core shifts between two fre-quencies 𝜈 and 𝜈 : Δ 𝑟 obs Δ 𝑟 true = 𝐶 (cid:18) 𝐶 𝐶 (cid:19) − 𝑘 𝜈 − / 𝑘 r , eff − 𝜈 − / 𝑘 r , eff 𝜈 − / 𝑘 r − 𝜈 − / 𝑘 r (7)where the effective coefficient of the frequency dependence is: 𝑘 r , eff = 𝑘 r /( + 𝑘 ( 𝑘 r − )) (8)equals to the true 𝑘 r only for 𝑘 r =
1. In this case, the ratio ofthe observed to the true core shift is frequency independent. As Δ 𝑟 true ∝ 𝐶 and FWHM ∝ 𝐶 , this case also can be written via thetrue core shift in units of the mean restoring beam FWHM: Δ 𝑟 obs Δ 𝑟 true = 𝐶 Δ 𝑟 (cid:18) Δ 𝑟 true (cid:104) FWHM (cid:105) (cid:19) − 𝑘 (9)where 𝐶 Δ 𝑟 = (cid:18) 𝜈 + 𝜈 ( 𝜈 − 𝜈 ) (cid:19) − 𝑘 𝐶 (10)Here (cid:104) FWHM (cid:105) — mean size of the beam for two frequencies and 𝐶 — constant from Equation 5. For 15.4 and 8.1 GHz we obtain 𝐶 Δ 𝑟 = .
41 and for 8.1 and 2.3 GHz 𝐶 Δ 𝑟 = .
53. This relation shouldhold at any pair of frequencies provided that the true core shift andthe resolution scales linearly with the wavelength. The correspondingdependence is shown in Fig. 5 for the core shift between 15.4 and8.1 GHz with its fit by Equation 9, which gives 𝐶 Δ 𝑟 = . ± . 𝑘 = . ± .
02 for CG. This agrees with our estimates of 𝐶 Δ 𝑟 from Equation 10 and fits of 𝐶 from Equation 5. For EG, weobtained 𝐶 Δ 𝑟 = . ± .
03 and 𝑘 = . ± .
02. For both CG andEG core templates, the overestimation of the core shift is the mostpronounced for sources with the largest apparent opening angles 𝜙 app (Fig. 6). In our artificial sample, these are the sources with thesmallest jet viewing angles, i.e. blazars. Radio galaxies are observedat larger viewing angles, have smaller 𝜙 app (Pushkarev et al. 2017)and, consequently, practically unbiased core shifts. MNRAS000
02. For both CG andEG core templates, the overestimation of the core shift is the mostpronounced for sources with the largest apparent opening angles 𝜙 app (Fig. 6). In our artificial sample, these are the sources with thesmallest jet viewing angles, i.e. blazars. Radio galaxies are observedat larger viewing angles, have smaller 𝜙 app (Pushkarev et al. 2017)and, consequently, practically unbiased core shifts. MNRAS000 , 1–11 (2020)
I. N. Pashchenko et al.
Figure 4.
Dependence of the ratio between the observed 𝑟 obs and true 𝑟 true core projected distance on the ratio of the 𝑟 true to FWHM of the beam for circular( left ) and elliptical ( right ) core templates. The points represent data at 15.4, 8.1 and 2.3 GHz (marked by blue , orange and green colour). The thin red linesrepresent the samples from the posterior distribution of the parameters of the fitted model (Equation 5). These lines are close and merging on this plot, thus thewidth of the red curve shows the spread of the model prediction. For the core shift between 8.1 and 2.3 GHz, the correspondingvalues are 𝑘 = . ± .
02 and 𝐶 Δ 𝑟 = . ± .
02 for CG. This isslightly different than estimated from Equation 5 and Equation 10.It could be due to a different VLBI array configuration at 2.3 GHzrelative to other bands used as a template in our simulations.
A recent extensive study of the core shift effect variability (Plavinet al. 2019a) analysed possible systematic biases as well. However,there are several differences in the analysis and the assumptions withour work. First, they focused on 2.3 and 8.6 GHz data obtained fromdifferent configurations of VLBI arrays. However, more importantly,is that they have assumed that the bias in core shift measurements isproportional to the beam FWHM averaged between the two frequen-cies, whereas we chose the dependence that fits simulated data well.Moreover, they employed a different approach to the estimation ofthe core parameters.Plavin et al. (2019a) estimated the proportionality coefficient be-tween core shift bias and mean beam FWHM to be 0 .
14. Assumingthat the effect at each frequency is proportional to its beam size, thisimplies that the bias of the core position is 0 . · (cid:104) FWHM (cid:105) . Apply-ing the same approach to the core parameter estimation as in Plavinet al. (2019a) to our simulated sample, we obtain the median bias of0 . · (cid:104) FWHM (cid:105) . Note that there are other possible sources of uncer-tainties (see Plavin et al. 2019a), which are out of the scope of thiswork. However, if the close agreement does not occur by chance, itsuggests that the remaining measurement uncertainties only increasethe statistical variance and do not make the results more biased. Italso implies that significant beam elongation present in Plavin et al.(2019a) does not significantly change the effect.The assumption adopted in Plavin et al. (2019a) leads to differentdependence of the fractional core shift bias on the true core shiftbetween 2.3 and 8.1 GHz: Δ 𝑟 obs Δ 𝑟 true = + . (cid:18) Δ 𝑟 true (cid:104) FWHM (cid:105) (cid:19) − . (11)This is more consistent with our results obtained for the elliptical Gaussian core model (Fig. 5, left), however, with a steeper depen-dence on the true core shift − − . 𝛿 -functions, alwayscontribute to a correlated flux density at the longest baselines whichcould significantly influence the subsequent estimation of the coreparameters. If only the measurements of the core shift between two frequencies 𝜈 and 𝜈 are available, one can assume 𝑘 r = 1 and obtain the relationfor compensating the bias from Equation 9 and Equation 10: Δ 𝑟 true (cid:104) FWHM (cid:105) = (cid:18) 𝜈 + 𝜈 ( 𝜈 − 𝜈 ) (cid:19) 𝑘 /( − 𝑘 ) 𝐶 /( 𝑘 − ) (cid:18) Δ 𝑟 obs (cid:104) FWHM (cid:105) (cid:19) /( − 𝑘 ) . (12) MNRAS , 1–11 (2020) ias of the core shift in AGN jets Figure 5.
Dependence of the ratio between the observed Δ 𝑟 obs and true Δ 𝑟 true core shifts between 15.4 and 8.1 GHz on the ratio of the Δ 𝑟 true to mean FWHMof the beams for circular ( left ) and elliptical ( right ) core templates. The thin red lines represent the samples from the posterior distribution of the parameters ofthe fitted model (Equation 9). Here these lines are close and merging, thus the width of the red curve shows the spread of the model prediction. Figure 6.
Dependence of the ratio between the observed and true core shiftsbetween 15.4 and 8.1 GHz on the apparent jet opening angle for circular( blue ) and elliptical ( orange ) core templates.
Here constants 𝐶 and 𝑘 are estimated in subsection 4.3 for bothcircular and elliptical core models. Note that Equation 12 can beemployed for any pair of frequencies provided both true core shiftand the resolution scale linearly with the wavelength. In Fig. 7, weplot the measured and true core shifts between 15.4 and 8.1 GHztogether with those corrected for the bias using Equation 12. Herethe residuals between the corrected and true values correspond to 𝜎 MAD = .
02 and 0.01 mas for circular and elliptical Gaussians(Fig. 8), where 𝜎 MAD is the robust estimate of standard deviation(Leys et al. 2013). The biases of both compensations are 0.010 and0.005 mas for CG and EG, respectively.If observations at more than two frequencies are analysed, thenEquation 6 can be fitted to the measured core shift frequency de-pendence. Here 𝐶 and 𝑘 are estimated in this work (see subsec-tion 4.3) and 𝐶 can be estimated from fitting the inverse frequencydependence to beam sizes. After 𝑘 r , eff and 𝐶 are found, the canthen be used to estimate true core shifts using the original relation 𝑟 true ( 𝜈 ) = 𝐶 𝜈 − / 𝑘 r . Finally, we warn readers that the estimates of 𝑘 and 𝐶 (thus 𝐶 Δ 𝑟 ) could depend on the VLBI array configuration, as shown insubsection 4.3. This could introduce uncertainty in the core shiftcorrection procedure described here. Thus the most rigorous wayto compensate for the core shift bias for a particular VLBI arrayconfiguration would be to estimate coefficients 𝑘 and 𝐶 Δ 𝑟 performinga series of simulations using a jet model (e.g. Blandford & Königl1979) and a specific ( 𝑢 , 𝑣 )-coverage. The measured shifts are widely used in estimating the referencevalues (at 1 pc from the apex) of the magnetic field 𝐵 and particledensity 𝐾 assuming equipartition (Lobanov 1998; Hirotani 2005;O’Sullivan & Gabuzda 2009; Pushkarev et al. 2012): 𝐵 eq1 ∝ ((cid:52) 𝑟 obs ) . 𝑘 r . (13)Here we should account for both biases: one is due to defining thecore shift as the 𝜏 = (cid:52) 𝑟 in different directions and compensateeach other for (cid:52) 𝑟 true / FWHM = 0.1 for CG and 0.4 for EG accordingto Fig. 5 and fits with Equation 9. Note that the magnetic field value 𝐵 core at the core itself turns out less biased. The distance of the corefrom the jet apex 𝑟 core ∝ (cid:52) 𝑟 obs for 𝑘 r = 1 (Lobanov 1998). Thus themagnetic field in the core weakly depends on the measured core shiftas 𝐵 eqcore = 𝐵 eq1 𝑟 − ∝ ((cid:52) 𝑟 obs ) − / .O’Sullivan & Gabuzda (2009) have estimated the 1 𝜎 uncertaintyof the magnetic field strength changing from 15% at 15.4 GHz to100% and larger at 5 GHz, assuming equipartition and only consider-ing the errors in the core shift measurements. More robust estimatesshould include others sources of errors as those of Doppler factors,jet opening and viewing angles, and spectral indices. Thus, the biasesconsidered here should be treated as a lower limit of the accuracy ofthe magnetic field estimation. However, as averaging across multipleobservations does not decrease a bias, it has certain astrophysicalimplications discussed below. MNRAS000
02 and 0.01 mas for circular and elliptical Gaussians(Fig. 8), where 𝜎 MAD is the robust estimate of standard deviation(Leys et al. 2013). The biases of both compensations are 0.010 and0.005 mas for CG and EG, respectively.If observations at more than two frequencies are analysed, thenEquation 6 can be fitted to the measured core shift frequency de-pendence. Here 𝐶 and 𝑘 are estimated in this work (see subsec-tion 4.3) and 𝐶 can be estimated from fitting the inverse frequencydependence to beam sizes. After 𝑘 r , eff and 𝐶 are found, the canthen be used to estimate true core shifts using the original relation 𝑟 true ( 𝜈 ) = 𝐶 𝜈 − / 𝑘 r . Finally, we warn readers that the estimates of 𝑘 and 𝐶 (thus 𝐶 Δ 𝑟 ) could depend on the VLBI array configuration, as shown insubsection 4.3. This could introduce uncertainty in the core shiftcorrection procedure described here. Thus the most rigorous wayto compensate for the core shift bias for a particular VLBI arrayconfiguration would be to estimate coefficients 𝑘 and 𝐶 Δ 𝑟 performinga series of simulations using a jet model (e.g. Blandford & Königl1979) and a specific ( 𝑢 , 𝑣 )-coverage. The measured shifts are widely used in estimating the referencevalues (at 1 pc from the apex) of the magnetic field 𝐵 and particledensity 𝐾 assuming equipartition (Lobanov 1998; Hirotani 2005;O’Sullivan & Gabuzda 2009; Pushkarev et al. 2012): 𝐵 eq1 ∝ ((cid:52) 𝑟 obs ) . 𝑘 r . (13)Here we should account for both biases: one is due to defining thecore shift as the 𝜏 = (cid:52) 𝑟 in different directions and compensateeach other for (cid:52) 𝑟 true / FWHM = 0.1 for CG and 0.4 for EG accordingto Fig. 5 and fits with Equation 9. Note that the magnetic field value 𝐵 core at the core itself turns out less biased. The distance of the corefrom the jet apex 𝑟 core ∝ (cid:52) 𝑟 obs for 𝑘 r = 1 (Lobanov 1998). Thus themagnetic field in the core weakly depends on the measured core shiftas 𝐵 eqcore = 𝐵 eq1 𝑟 − ∝ ((cid:52) 𝑟 obs ) − / .O’Sullivan & Gabuzda (2009) have estimated the 1 𝜎 uncertaintyof the magnetic field strength changing from 15% at 15.4 GHz to100% and larger at 5 GHz, assuming equipartition and only consider-ing the errors in the core shift measurements. More robust estimatesshould include others sources of errors as those of Doppler factors,jet opening and viewing angles, and spectral indices. Thus, the biasesconsidered here should be treated as a lower limit of the accuracy ofthe magnetic field estimation. However, as averaging across multipleobservations does not decrease a bias, it has certain astrophysicalimplications discussed below. MNRAS000 , 1–11 (2020)
I. N. Pashchenko et al.
Figure 7.
Raw observed (blue) and corrected (orange) core shifts vs true core shifts between 15.4 and 8.1 GHz for circular ( left ) and elliptical Gaussian ( right )according to Equation 12. The dashed line represents the equality between the observed and true values.
Figure 8.
Residuals between the true and measured but bias-corrected coreshifts between 15.4 and 8.1 GHz for the circular ( blue ) and elliptical Gaussiancore models ( orange ). Nalewajko et al. (2014) found that magnetic fields inferred fromEquation 13 are smaller than the 𝐵 estimates from the modellingof the spectral energy distribution (SED) by a factor of 2. This re-sult is based on a sample of blazars from Zamaninasab et al. (2014)with a core shift estimated by Pushkarev et al. (2012). These authorspropose that jets could have a non-uniform distribution of the mag-netic field across the jet. An alternative explanation is that magneticfields inferred from the core shift measurements are overestimated.Nalewajko et al. (2014) associate this with radio cores not being pho-tospheres due to the synchrotron self-absorbtion (SSA) process, butdue to a break in a lower energy cutoff of the energy distribution ofthe emitting particles. In our simulations, the sources with the small-est viewing angles are the most biased in terms of the core shift. Thiscould reconcile the lower measured 𝐵 with the corresponding SEDestimates within the SSA model of radio cores. Zamaninasab et al. (2014) found observational evidence for themagnetically arrested disc (MAD) (Bisnovatyi-Kogan & Ruzmaikin1974; Narayan et al. 2003) scenario in which magnetic fluxes ofthe jet Φ j ∝ 𝐵 p 𝑅 should be comparable to magnetic fluxes Φ MAD threading a black hole. Here 𝐵 p is the poloidal magnetic field in the jetat some distance where the jet radius equals to 𝑅 j . As Φ j is estimatedusing magnetic field amplitudes (Equation 13), it is overestimated forblazars. Another assumption lowering Φ j was that Γ 𝜃 ∝
1. Jorstadet al. (2005); Pushkarev et al. (2009); Clausen-Brown et al. (2013)estimated Γ 𝜃 ∼ . Φ j on the order of the magnitude for blazars (see also Finke 2019). A magnetic field can be also estimated from the core shift withoutequipartition assumptions using equation (8) from Zdziarski et al.(2015): 𝐵 noeq ( 𝑟 ) (cid:39) . × − 𝐷 𝐿 𝛿 ((cid:52) 𝑟 obs ) tan 𝜙𝑟 ( 𝜈 − − 𝜈 − ) [( + 𝑧 ) sin 𝜃 ] 𝐹 𝜈 (14)where 𝑟 is distance in pc, 𝐹 𝜈 — flux density of an inner jet integratedover both its optically thick and thin parts in Jy, 𝐷 𝐿 — luminositydistance in Mpc, 𝜈 < 𝜈 — frequencies of observations in GHz,and the numerical coefficient is calculated for the exponent of theemitting particles energy distribution 𝑠 =
2. Strong dependence on (cid:52) 𝑟 obs implies that the found systematic error of the core shift couldbias the magnetic field estimates up to an order of magnitude.Zdziarski et al. (2015) estimated magnetic fields from core shiftmeasurements using both equipartition (Equation 13) and a moregeneral approach (Equation 14). They found (cid:104) 𝐵 noeq / 𝐵 eq (cid:105) = 1.6for blazars and (cid:104) 𝐵 noeq / 𝐵 eq (cid:105) = 0.09 for radio galaxies. This ratiois modified in the presence of the core shift bias as (cid:104) 𝐵 noeq / 𝐵 eq (cid:105) = ((cid:52) 𝑟 obs /(cid:52) 𝑟 true ) − . ×(cid:104) 𝐵 unbiasednoeq / 𝐵 unbiasedeq (cid:105) . Here 𝐵 unbiased — mag-netic field estimate which would be obtained if core shift measure-ments were unbiased. From our sample, the mean of (cid:52) 𝑟 obs /(cid:52) 𝑟 true MNRAS , 1–11 (2020) ias of the core shift in AGN jets for sources with 𝜃 > ◦ (i.e. radio galaxies), corrected for biasestimated in section 3, is 0.73. Then, for blazars to have the sameratio of the unbiased estimates (cid:104) 𝐵 unbiasednoeq / 𝐵 unbiasedeq (cid:105) , the mean ra-tio (cid:104)(cid:52) 𝑟 obs /(cid:52) 𝑟 true (cid:105) should be 1.44. In our simulated sample, this isachieved for sources with 𝜃 < . ◦ . Zdziarski et al. (2015) attributethe effect of different (cid:104) 𝐵 noeq / 𝐵 eq (cid:105) for blazars and radio galaxies todifferent viewing angles 𝜃 and the corresponding difference in rela-tivistic beaming. However, it can also be explained by the systematicsof the core shift measurements.Pjanka et al. (2017) found that jet power 𝑃 j values computed fromthe core shift estimated according to Zdziarski et al. (2015) are anorder of magnitude lower than those estimated through the extendedradio lobe luminosity. As the jet power 𝑃 j ∝ 𝐵 and 𝐵 ∝ (cid:52) 𝑟 ,this disagreement could be explained if the measured core shifts areoverestimated by a factor of 1.26, consistent with our findings. Plavin et al. (2019a) found that core shifts vary significantly duringflares. Their results imply that flares are due to an increase in theemitting particle density and a decrease of the magnetic field. Asfollows from the flare model presented in Plavin et al. (2019a), seetheir fig. 13, when a flare reaches the core at some frequency, itshifts the core downstream. This increases the core separation fromthe jet base and, according to subsection 4.2, decreases its fractionalbias. Thus the quiescent state has the most fractionally biased coreshift estimate. Contrary to expectations, it could be that the mostaccurate core shift measurements in terms of a fractional bias arepossible during the flaring state. Not accounting for this effect couldpotentially underestimate departures from the equipartition in coresduring flares, as found in Plavin et al. (2019a).Algaba et al. (2017) found evidence for non-conical jet shapes inradio-loud AGN jets in the core region using core sizes 𝑅 measured atat least four frequencies. The distances of radio cores from the centralengines were estimated using measured core shifts from Pushkarevet al. (2012) and assuming 𝑟 core ∝ (cid:52) 𝑟 obs / 𝜈 (i.e. 𝑘 r = 1). For 56sources with fit statistic 𝑅 < .
85, the obtained outflow geometry ischaracterised by the coefficient 𝜖 of the radial core size dependence 𝑅 ∝ 𝑟 𝜖 , which is consistent with a quasi-parabolic outflow withmean (cid:104) 𝜖 (cid:105) = 0.87. From Equation 9, the true core shift is related to theobserved one as Δ 𝑟 true ∝ ( Δ 𝑟 obs ) /( − 𝑘 ) . Then, in the presence of abias, the underlying dependence 𝑅 ∝ 𝑟 𝜖 ∝ Δ 𝑟 𝜖 true transforms to 𝑅 ∝( Δ 𝑟 obs ) 𝜖 /( − 𝑘 ) , changing the effective exponent to 𝜖 eff = 𝜖 /( − 𝑘 ) .In subsection 4.3, we estimated 𝑘 ≈ .
42 in the case of fitting the corewith a circular Gaussian template. Thus the exponent (cid:104) 𝜖 eff (cid:105) = . (cid:104) 𝜖 (cid:105) = .
50, implying parabolic jet shape. However, Algaba et al. (2017)also employed the elliptical core model in their estimation of (cid:104) 𝜖 (cid:105) .According to subsection 4.3, modelling the core with an ellipticalGaussian results in a lower 𝑘 and, consequently, higher (cid:104) 𝜖 (cid:105) . Thus ourestimate of the intrinsic (cid:104) 𝜖 (cid:105) should be considered as a lower bound.The overestimation of 𝜖 due to the core shift bias could be inter-preted as follows. In the parabolic accelerating outflow, the magneticfield and the particle density gradients are lower than those in a con-ical constant speed jet. This leads to a steeper than 𝜈 − frequencydependence of the core shift for distances closer to the jet base (i.e.higher frequencies), as demonstrated by Porth et al. (2011). Since res-olution scales linearly with wavelength, and the bias is determinedby the ratio of the true core shift to the resolution (subsection 4.3),core shift estimates at higher frequencies will have a larger fractionalbias, effectively increasing the measured 𝜖 . However, this assumes that the model brightness distribution is close to that of Blandford& Königl (1979) for a conical jet used in our simulations. In fact,Abellán et al. (2018) found more elongated cores at 43 GHz relativeto 15 GHz for a fraction of sources in their complete S5 polar capsample (Eckart et al. 1986). This also could be a manifestation of aradio core at 43 GHz residing in the parabolic accelerating part ofthe jet.If a source demonstrates flaring activity, it is possible to derivethe kinematics of the jet innermost regions based on the measuredcore shift and multi-frequency time delays (Kudryavtseva et al. 2011;Kutkin et al. 2014, 2018, 2019). Assuming that the peak of a flare at agiven frequency is observed when a moving jet feature passes throughthe region of the core at that frequency, it is possible to estimate theapparent projected speed of the component by comparing the lightcurve time delay, (cid:52) 𝑡 , with the core shift for a pair of frequencies, (cid:52) 𝑟 : 𝜇 app = (cid:52) 𝑟 / (cid:52) 𝑡 . Kutkin et al. (2014) obtained an estimate ofthe apparent angular velocity 𝜇 app = 0.7 mas yr − , which is 2-8times larger than the speed estimated using VLBI kinematic data(Lister et al. 2013; Jorstad et al. 2010). These papers propose severalexplanations for this discrepancy, e.g. truly fast components or VLBInot measuring true plasma speed. Systematically overestimated coreshifts could at least partially explain this controversy. This especiallyholds for blazars, inclined at small angles to the line of sight anddemonstrating the largest fractional bias. Precise positions of AGNs are measured with VLBI using eithersingle-band phase delays or multi-band group delays of the arrivingwaves. The effects of the source structure and its frequency depen-dence on these positions were extensively studied: see, e.g. Porcas(2009), Plank et al. (2016), Xu et al. (2017), Anderson & Xu (2018).Single-band measurements are sensitive to the brightest part posi-tion itself; thus, they are affected by the full magnitude of the coreposition bias addressed in this paper.Commonly used multi-band astrometry is sensitive to the fre-quency dependence of the core offset (Porcas 2009): the more it de-parts from 𝑟 ∼ 𝜈 − , the larger the position measurement errors are.Such departures often arise due to intrinsic physical properties ofAGNs: see Lobanov (1998) for theoretical modelling and Kudryavt-seva et al. (2011); Algaba et al. (2012); Kutkin et al. (2014); Kara-manavis et al. (2016) for specific observational results. Plavin et al.(2019a) have shown that this frequency dependence is also affectedby strong flares which commonly happen in AGN jets. Our analysissuggests that systematic errors of astrometric positions in all thesecases are effectively amplified. Indeed, as illustrated in Fig. 9, anydeparture of the exponent in 𝑟 ∼ 𝜈 − 𝑘 r from 𝑘 r = MNRAS000
50, implying parabolic jet shape. However, Algaba et al. (2017)also employed the elliptical core model in their estimation of (cid:104) 𝜖 (cid:105) .According to subsection 4.3, modelling the core with an ellipticalGaussian results in a lower 𝑘 and, consequently, higher (cid:104) 𝜖 (cid:105) . Thus ourestimate of the intrinsic (cid:104) 𝜖 (cid:105) should be considered as a lower bound.The overestimation of 𝜖 due to the core shift bias could be inter-preted as follows. In the parabolic accelerating outflow, the magneticfield and the particle density gradients are lower than those in a con-ical constant speed jet. This leads to a steeper than 𝜈 − frequencydependence of the core shift for distances closer to the jet base (i.e.higher frequencies), as demonstrated by Porth et al. (2011). Since res-olution scales linearly with wavelength, and the bias is determinedby the ratio of the true core shift to the resolution (subsection 4.3),core shift estimates at higher frequencies will have a larger fractionalbias, effectively increasing the measured 𝜖 . However, this assumes that the model brightness distribution is close to that of Blandford& Königl (1979) for a conical jet used in our simulations. In fact,Abellán et al. (2018) found more elongated cores at 43 GHz relativeto 15 GHz for a fraction of sources in their complete S5 polar capsample (Eckart et al. 1986). This also could be a manifestation of aradio core at 43 GHz residing in the parabolic accelerating part ofthe jet.If a source demonstrates flaring activity, it is possible to derivethe kinematics of the jet innermost regions based on the measuredcore shift and multi-frequency time delays (Kudryavtseva et al. 2011;Kutkin et al. 2014, 2018, 2019). Assuming that the peak of a flare at agiven frequency is observed when a moving jet feature passes throughthe region of the core at that frequency, it is possible to estimate theapparent projected speed of the component by comparing the lightcurve time delay, (cid:52) 𝑡 , with the core shift for a pair of frequencies, (cid:52) 𝑟 : 𝜇 app = (cid:52) 𝑟 / (cid:52) 𝑡 . Kutkin et al. (2014) obtained an estimate ofthe apparent angular velocity 𝜇 app = 0.7 mas yr − , which is 2-8times larger than the speed estimated using VLBI kinematic data(Lister et al. 2013; Jorstad et al. 2010). These papers propose severalexplanations for this discrepancy, e.g. truly fast components or VLBInot measuring true plasma speed. Systematically overestimated coreshifts could at least partially explain this controversy. This especiallyholds for blazars, inclined at small angles to the line of sight anddemonstrating the largest fractional bias. Precise positions of AGNs are measured with VLBI using eithersingle-band phase delays or multi-band group delays of the arrivingwaves. The effects of the source structure and its frequency depen-dence on these positions were extensively studied: see, e.g. Porcas(2009), Plank et al. (2016), Xu et al. (2017), Anderson & Xu (2018).Single-band measurements are sensitive to the brightest part posi-tion itself; thus, they are affected by the full magnitude of the coreposition bias addressed in this paper.Commonly used multi-band astrometry is sensitive to the fre-quency dependence of the core offset (Porcas 2009): the more it de-parts from 𝑟 ∼ 𝜈 − , the larger the position measurement errors are.Such departures often arise due to intrinsic physical properties ofAGNs: see Lobanov (1998) for theoretical modelling and Kudryavt-seva et al. (2011); Algaba et al. (2012); Kutkin et al. (2014); Kara-manavis et al. (2016) for specific observational results. Plavin et al.(2019a) have shown that this frequency dependence is also affectedby strong flares which commonly happen in AGN jets. Our analysissuggests that systematic errors of astrometric positions in all thesecases are effectively amplified. Indeed, as illustrated in Fig. 9, anydeparture of the exponent in 𝑟 ∼ 𝜈 − 𝑘 r from 𝑘 r = MNRAS000 , 1–11 (2020) I. N. Pashchenko et al.
Figure 9.
Blue line shows the dependence of the observed vs the true 𝑘 r . Theline of equality is shown by the orange color. In this paper, we examined and found significant biases in measure-ments of the frequency dependent core shift effect in AGN jets andsuggested a method to compensate for them. To estimate the modelapproximation bias, we performed simulations using the artificiallygenerated multi-frequency VLBI data sets based on the brightnessdistributions derived from the Blandford & Königl (1979) jet modeland the parameters of the complete MOJAVE sample. We employedtwo approaches to modelling the core with a Gaussian template:fitting source structure with a number of Gaussians and fitting thecore with a single Gaussian after subtracting the extended emission.Our results show that the estimated core positions and core shiftsare positively biased. For example, for the 15.4-8.1 GHz frequencyinterval, the typical bias is ≈ 𝑘 r ofthe core shift frequency dependence 𝑟 ( 𝜈 ) ∼ 𝜈 − 𝑘 r is biased for the 𝑘 r ≠ 𝑘 r = 𝜏 = DATA AVAILABILITY
The datasets were derived from sources in the public do-main: , http://astrogeo.org/vlbi_images/ . ACKNOWLEDGEMENTS
Authors thank the anonymous referee for comments and sugges-tions which helped to significantly clarify the paper. Authors thankEduardo Ros for extensive comments on the manuscript and ElenaBazanova for language corrections. The Russian Science Foundationgrant 16-12-10481 has supported the core shift systematic errorsanalysis. The Russian Basic Research Foundation grant 19-32-90141has supported estimates of the AGN inner jet parameters. This re-search has made use of data from the MOJAVE database that ismaintained by the MOJAVE team (Lister et al. 2018). This researchhas made use of NASA’s Astrophysics Data System. This researchhas made use the Astrogeo VLBI FITS image database.This research made use of
Astropy , a community-developed corePython package for Astronomy (Astropy Collaboration et al. 2013),
Numpy (Van Der Walt et al. 2011),
Scipy (Jones et al. 2001), scikit-learn (Pedregosa et al. 2011) and
PyMC3 (Salvatier et al. 2016).
Matplotlib
Python package (Hunter 2007) was used for generatingall plots in this paper.
REFERENCES
Abellán F. J., Martí-Vidal I., Marcaide J. M., Guirado J. C., 2018, A&A, 614,A74Agarwal A., et al., 2017, MNRAS, 469, 813Algaba J. C., Gabuzda D. C., Smith P. S., 2012, MNRAS, 420, 542Algaba J. C., Nakamura M., Asada K., Lee S. S., 2017, ApJ, 834, 65Anderson J. M., Xu M. H., 2018, Journal of Geophysical Research (SolidEarth), 123, 10,162Astropy Collaboration et al., 2013, A&A, 558, A33Bisnovatyi-Kogan G. S., Ruzmaikin A. A., 1974, Ap&SS, 28, 45Blandford R. D., Königl A., 1979, ApJ, 232, 34Briggs D. S., 1995, PhD thesis, The New Mexico Institue of Mining andTechnology, Socorro, New MexicoCharlot P., 1990, AJ, 99, 1309Clausen-Brown E., Savolainen T., Pushkarev A. B., Kovalev Y. Y., ZensusJ. A., 2013, A&A, 558, A144Conover W., 1999, Practical nonparametric statistics, 3. ed edn.Wiley series in probability and statistics, Wiley, New York,NY [u.a.], http://gso.gbv.de/DB=2.1/CMD?ACT=SRCHA&SRT=YOP&IKT=1016&TRM=ppn+24551600X&sourceid=fbw_bibsonomy
Croke S. M., Gabuzda D. C., 2008, MNRAS, 386, 619Dormand J., Prince P., 1980, Journal of Computational and Applied Mathe-matics, 6, 19Eckart A., Witzel A., Biermann P., Johnston K. J., Simon R., Schalinski C.,Kuhr H., 1986, A&A, 168, 17Finke J. D., 2019, The Astrophysical Journal, 870, 28Fromm C. M., Ros E., Perucho M., Savolainen T., Mimica P., Kadler M.,Lobanov A. P., Zensus J. A., 2013, A&A, 557, A105Hada K., Doi A., Kino M., Nagai H., Hagiwara Y., Kawaguchi N., 2011,Nature, 477, 185Hinshaw G., et al., 2013, ApJS, 208, 19Hirotani K., 2005, ApJ, 619, 73Hovatta T., Lister M. L., Aller M. F., Aller H. D., Homan D. C., KovalevY. Y., Pushkarev A. B., Savolainen T., 2012, AJ, 144, 105Hovatta T., et al., 2014, AJ, 147, 143Hunter J. D., 2007, Computing In Science & Engineering, 9, 90Jones E., Oliphant T., Peterson P., et al., 2001, SciPy: Open source scientifictools for Python,
Jorstad S. G., et al., 2005, AJ, 130, 1418Jorstad S. G., et al., 2010, ApJ, 715, 362Karamanavis V., et al., 2016, A&A, 590, A48Konigl A., 1981, ApJ, 243, 700Kovalev Y. Y., et al., 2005, AJ, 130, 2473Kovalev Y. Y., Lobanov A. P., Pushkarev A. B., Zensus J. A., 2008, A&A,483, 759MNRAS , 1–11 (2020) ias of the core shift in AGN jets Kovalev Y. Y., Petrov L., Plavin A. V., 2017, A&A, 598, L1Kovalev Y. Y., Pushkarev A. B., Nokhrina E. E., Plavin A. V., Beskin V. S.,Chernoglazov A. V., Lister M. L., Savolainen T., 2020, MNRAS, 495,3576Kudryavtseva N. A., Gabuzda D. C., Aller M. F., Aller H. D., 2011, MNRAS,415, 1631Kutkin A. M., et al., 2014, MNRAS, 437, 3396Kutkin A. M., et al., 2018, MNRAS, 475, 4994Kutkin A. M., Pashchenko I. N., Sokolovsky K. V., Kovalev Y. Y., Aller M. F.,Aller H. D., 2019, MNRAS, 486, 430Leys C., Ley C., Klein O., Bernard P., Licata L., 2013, Journal of ExperimentalSocial Psychology, 49, 764Lister M. L., et al., 2009a, AJ, 137, 3718Lister M. L., et al., 2009b, AJ, 138, 1874Lister M. L., et al., 2013, AJ, 146, 120Lister M. L., Aller M. F., Aller H. D., Hodge M. A., Homan D. C., KovalevY. Y., Pushkarev A. B., Savolainen T., 2018, ApJS, 234, 12Lister M. L., et al., 2019, ApJ, 874, 43Lobanov A. P., 1998, A&A, 330, 79Longair M. S., 1994, High Energy Astrophysics, 2 edn. Vol. 2, CambridgeUniversity Press, doi:10.1017/CBO9781139170505Marcaide J. M., Shapiro I. I., 1984, ApJ, 276, 56Marr J. M., Taylor G. B., Crawford III F., 2001, ApJ, 550, 160Marscher A. P., 2008, in T. A. Rector & D. S. De Young ed., AstronomicalSociety of the Pacific Conference Series Vol. 386, Extragalactic Jets:Theory and Observation from Radio to Gamma Ray. p. 437Nalewajko K., Sikora M., Begelman M. C., 2014, ApJ, 796, L5Narayan R., Igumenshchev I. V., Abramowicz M. A., 2003, PASJ, 55, L69Nokhrina E. E., 2017, MNRAS, 468, 2372Nokhrina E. E., 2020, in Asada K., de Gouveia Dal Pino E., Giroletti M.,Nagai H., Nemmen R., eds, IAU Symposium Vol. 342, IAU Symposium.pp 197–200, doi:10.1017/S1743921318006087Nokhrina E. E., Beskin V. S., Kovalev Y. Y., Zheltoukhov A. A., 2015,MNRAS, 447, 2726O’Sullivan S. P., Gabuzda D. C., 2009, MNRAS, 400, 26Pashchenko I. N., Plavin A. V., 2019, MNRAS, 488, 939Pearson T. J., Readhead A. C. S., 1984, ARA&A, 22, 97Pedregosa F., et al., 2011, Journal of Machine Learning Research, 12, 2825Pjanka P., Zdziarski A. A., Sikora M., 2017, MNRAS, 465, 3506Plank L., Shabala S. S., McCallum J. N., Krásná H., Petrachenko B.,Rastorgueva-Foi E., Lovell J. E. J., 2016, MNRAS, 455, 343Plavin A. V., Kovalev Y. Y., Pushkarev A. B., Lobanov A. P., 2019a, MNRAS,485, 1822Plavin A. V., Kovalev Y. Y., Petrov L. Y., 2019b, ApJ, 871, 143Porcas R. W., 2009, A&A, 505, L1Porth O., Fendt C., Meliani Z., Vaidya B., 2011, ApJ, 737, 42Pushkarev A. B., Kovalev Y. Y., Lister M. L., Savolainen T., 2009, A&A,507, L33Pushkarev A. B., Hovatta T., Kovalev Y. Y., Lister M. L., Lobanov A. P.,Savolainen T., Zensus J. A., 2012, A&A, 545, A113Pushkarev A. B., Kovalev Y. Y., Lister M. L., Savolainen T., 2017, MNRAS,468, 4992Salvatier J., Wiecki T. V., Fonnesbeck C., 2016, PeerJ Computer Science, 2,e55Shepherd M. C., 1997, in G. Hunt & H. Payne ed., Astronomical Societyof the Pacific Conference Series Vol. 125, Astronomical Data AnalysisSoftware and Systems VI. p. 77Sokolovsky K. V., Kovalev Y. Y., Pushkarev A. B., Lobanov A. P., 2011,A&A, 532, A38Van Der Walt S., Colbert S. C., Varoquaux G., 2011, Computing in Science& Engineering, 13, 22Walker R. C., Dhawan V., Romney J. D., Kellermann K. I., Vermeulen R. C.,2000, ApJ, 530, 233Wasserman L., 2010, All of Statistics: A Concise Course in Statistical Infer-ence. Springer Publishing Company, IncorporatedXu M. H., Heinkelmann R., Anderson J. M., Mora-Diaz J., Karbon M., SchuhH., Wang G. L., 2017, Journal of Geodesy, 91, 767 Zamaninasab M., Clausen-Brown E., Savolainen T., Tchekhovskoy A., 2014,Nature, 510, 126Zdziarski A. A., Sikora M., Pjanka P., Tchekhovskoy A., 2015, MNRAS, 451,927Zensus J. A., 1997, ARA&A, 35, 607This paper has been typeset from a TEX/L A TEX file prepared by the author.MNRAS000