A connection between $\gamma$-ray and parsec-scale radio flares in the blazar 3C 273

TL;DR

This study uses multi-frequency VLBA observations to locate the gamma-ray emission zone in blazar 3C 273 near the jet apex, revealing correlations with radio flares, core shifts, and jet dynamics.

Contribution

It provides the first detailed localization of gamma-ray emission in 3C 273 close to the jet base using combined kinematic and spectral analysis.

Findings

Gamma-ray emission zone is 2-7 pc upstream of the 7 mm core.

Brightness temperature in the core reached at least 10^13 K during the flare.

Jet bulk flow speed estimated at 12c, higher than apparent kinematic estimates.

Abstract

We present a comprehensive 5-43 GHz VLBA study of the blazar 3C 273 initiated after an onset of a strong $\gamma$-ray flare in this source. We have analyzed the kinematics of new-born components, light curves, and position of the apparent core to pinpoint the location of the $\gamma$-ray emission. Estimated location of the $\gamma$-ray emission zone is close to the jet apex, 2 pc to 7 pc upstream from the observed 7 mm core. This is supported by ejection of a new component. The apparent core position was found to be inversely proportional to frequency. The brightness temperature in the 7 mm core reached values up to at least $10^{13}$ K during the flare. This supports the dominance of particle energy density over that of magnetic field in the 7 mm core. Particle density increased during the radio flare at the apparent jet base, affecting synchrotron opacity. This manifested itself as an apparent core shuttle along the jet during the 7 mm flare. It is also shown that a region where optical depth decreases from $\tau\sim1$ to $\tau<<1$ spans over several parsecs along the jet. The jet bulk flow speed estimated at the level of 12c on the basis of time lags between 7 mm light curves of stationary jet features is 1.5 times higher than that derived from VLBI apparent kinematics analysis.

Tables

Table 1. Table 1: Epochs of VLBA multi-frequency observations carried out by the authors. HN antenna dropped out during the last epoch. Otherwise no major losses. Visibilities are counted with frequency and time (5 s) averaging.

Date	Number of antennas	Visibilities $\times {10}^3$ per band C, X, U, K, Q
2009 Aug 28	10	5.5, 7.5, 7, 8.5, 12
2009 Oct 25	10	7.5, 8, 9, 9, 15
2009 Dec 05	10	5, 6.5, 6, 5.5, 12
2010 Jan 26	9	6, 5.5, 6.5, 7, 10

Table 2. Table 2: Frequency setup used in our multi-frequency VLBA observations.

Band	Wavelength (cm)	Frequency (GHz)	Bandwidth (MHz)
Q1	0.7	43.2	32
K1	1.3	23.8	32
U1	2	15	32
X2	3.6	8.4	16
X1	3.7	8.1	16
C2	6.0	5.0	16
C1	6.5	4.6	16

Table 3. Table 3: Epochs and data quality of 7 mm VLBA observations, both from our multifrequency campaign and from the Boston University blazar monitoring. N eq subscript 𝑁 eq N_{\mathrm{eq}} – roughly estimated number of antennas after accounting for all flagged data. Number of visibilities is calculated after IF averaging, integration time 20 sec. “-ANTENNA” in comments stands for total loss. Epochs marked with an asterisk were excluded from subsequent analysis. A full version of the table is presented in online-only materials (Table 11).

Date	MJD	$N_{\mathrm{eq}}$	Number of visibilities	Comments
2008-01-17	54482	8	1829	-HN, NL 60% lost
2008-02-29	54525	8	2433	-HN, -FD
*2008-06-12	54629	6	1168	-BR, - LA, -PT
*2008-07-06	54653	7.5	1956	-BR, -OV
2008-08-15	54693	9.5	2408
$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$

Table 4. Table 4: Modelfit parameters for all epochs. 4.6 GHz data. Component 0 represents the apparent core, other components are numbered according to their distance from the core. Rad and theta are polar coordinates of the component with respect to the map center. Major axis is the FWHM of the major axis of a Gaussian component. Axial ratio is a ratio of minor axis to major axis of a two-dimentional Gaussian component and is unity for circular ones. Phi is an inclination angle of the major axis with respect to the South–North direction. The data for all epochs and frequencies are presented in online-only materials (Tables 12–18).

Epoch	Component	Flux density (Jy)	Rad (mas)	Theta (deg)	Major axis (mas)	Axial ratio	Phi (deg)
2009_08_28	0	6.05426	0.56779	50.40940	0	1.0	0.0
	1	5.71666	0.64843	-146.68800	0.39853	1.0	0.0
	2	5.30306	2.23402	-141.43400	1.1442	1.0	0.0
	3	1.15141	3.80385	-134.86500	0	1.0	0.0
	4	5.43772	5.46867	-126.17600	2.5023	1.0	0.0
	5	0.22616	7.37588	-111.46000	0	1.0	0.0
$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$

Table 5. Table 5: Parameters of the newborn components. For components c1 and c2 , which exhibited accelerated motion, radial speed refers to the middle epoch of the observations. Radial speed at the first occurrence of components c1 and c2 in the jet structure was 0.85 ± 0.02 plus-or-minus 0.85 0.02 0.85\pm 0.02   mas ​ yr − 1 mas superscript yr 1 \mathrm{mas\ yr^{-1}} and 0.73 ± 0.04 plus-or-minus 0.73 0.04 0.73\pm 0.04   mas ​ yr − 1 mas superscript yr 1 \mathrm{mas\ yr^{-1}} respectively. Brightness temperature is calculated in the source frame.

Component	Ejection epoch (MJD)	Radial speed ($\mathrm{mas}{\mathrm{yr}}^{-1}$)	Radial acceleration ($\mathrm{mas}{\mathrm{yr}}^{-2}$)	max($T_b$) (K)
c1	$54838\pm 19$	$1.17\pm 0.02$	$0.26\pm 0.04$	2.2e+11
c2	$54885\pm 43$	$0.96\pm 0.02$	$0.17\pm 0.05$	2.1e+11
c3	$55112\pm 17$	$0.83\pm 0.02$	–	4.5e+11
c4	$55189\pm 14$	$0.81\pm 0.02$	–	1.1e+12
c5	$55136\pm 36$	$0.53\pm 0.03$	–	4.6e+10
c6	$54823\pm 381$	$0.20\pm 0.07$	–	4.1e+10

Table 6. Table 6: Core parameters. T b core superscript subscript 𝑇 b core T_{\mathrm{b}}^{\mathrm{core}} is the brightness temperature of the 7 mm apparent core. Other parameters are estimated from the core-shift analysis. k r subscript 𝑘 r k_{\mathrm{r}} is the power law index; r 7 ​ m ​ m subscript 𝑟 7 m m r_{\mathrm{7mm}} – deprojected distance from the true base of the jet to the apparent core at 7 mm; B 1 subscript 𝐵 1 B_{1} – magnetic field strength at a distance of 1 pc along the jet from the true jet base ; B core 7 ​ mm subscript superscript 𝐵 7 mm core B^{\mathrm{7~{}mm}}_{\mathrm{core}} – magnetic field at the position of the apparent core at 7 mm; N 1 subscript 𝑁 1 N_{1} - particle density at 1 pc from the jet origin. B 1 subscript 𝐵 1 B_{1} , B core 7 ​ mm subscript superscript 𝐵 7 mm core B^{\mathrm{7~{}mm}}_{\mathrm{core}} , and N 1 subscript 𝑁 1 N_{1} are estimated assuming equipartition between the magnetic field and particles energy density. Italic font denotes formally estimated values for epochs with high 7 mm core brightness temperature (Fig.  8 ) which violates the equipartition assumption (see Sect.  4.3 ).

Epoch	$T_b^{\mathrm{core}}$ (${10}^{12}$ K)	$k_r$	${\Omega }_{r\nu }$ (pc GHz)	$r_{7\mathrm{mm}}$ (pc)	$B_1$ (G)	$B_{\mathrm{core}}^{7\mathrm{mm}}$ (G)	$N_1$ (cm-3)
2009 Aug 28	1.1	$0.85\pm 0.19$	$21.87\pm 7.66$	$2.43\pm 2.58$	$0.47\pm 0.12$	$0.19\pm 0.21$	1500
2009 Oct 25	4.7	$0.82\pm 0.25$	$15.00\pm 7.68$	$1.47\pm 2.17$	$0.34\pm 0.10$	$0.23\pm 0.35$	$800$
2009 Dec 05	6.1	$1.04\pm 0.29$	$9.63\pm 3.15$	$2.50\pm 2.64$	$0.55\pm 0.18$	$0.22\pm 0.24$	$2100$
2010 Jan 26	5.3	$1.17\pm 0.27$	$11.42\pm 2.62$	$4.37\pm 3.40$	$0.97\pm 0.30$	$0.22\pm 0.19$	$6700$

Table 7. Table 7: Magnetic field at the position of the apparent 7 mm core B core 7 ​ mm subscript superscript 𝐵 7 mm core B^{\mathrm{7~{}mm}}_{\mathrm{core}} (G) for different values of the particle density relative to N = 1500 𝑁 1500 N=1500   cm − 3 superscript cm 3 \mathrm{cm}^{-3} estimated assuming equipartition for the epoch 2009 Aug 28 (Table  6 ).

Epoch	$10N$	$100N$	$1000N$
2009 Oct 25	0.05	0.016	0.005
2009 Dec 05	0.08	0.025	0.008
2010 Jan 26	0.14	0.045	0.014

Table 8. Table 8: Average value of the spectral-index distribution taken along the ridge line for different frequency pairs.

${\alpha }_{43-24}$	${\alpha }_{24-15}$	${\alpha }_{15-8}$	${\alpha }_{8-5}$
$-0.73\pm 0.06$	$-0.81\pm 0.04$	$-0.81\pm 0.07$	$-0.68\pm 0.05$

Table 9. Table 9: Comparison of models for α 43 − 24 subscript 𝛼 43 24 \alpha_{43-24} spectral index transition zone. w AICc subscript 𝑤 AICc w_{\mathrm{AICc}} shows the relative probability of the best fit ramp function in comparison with the best fit step function, width (mas) is intrinsic extent of a region with linear descend of spectral index in the fitted ramp function. The last epoch is not included in the table because both models fail to adequately describe the spectral index slice.

Epoch	$w_{\mathrm{AICc}}$	Width (mas)
2009 Aug 28	${10}^5$	0.49
2009 Oct 25	${10}^8$	0.45
2009 Dec 05	$18$	0.33

Table 10. Table 10: Pearson correlation coefficient r 𝑟 r and the permutation test p − limit-from 𝑝 p- value for clusters K1 and K2 for the period of the 7 mm core flare B 7 subscript 𝐵 7 B_{7} ; for the period outside the flare; and for all the data together.

Cluster	Flare		Non-flare		All
	$r$	$p$	$r$	$p$	$r$	$p$
K1	$-0.9$	$2\times {10}^{-4}$	$-0.12$	$0.36$	$-0.57$	${10}^{-3}$
K2	$-0.93$	${10}^{-5}$	$-0.13$	$0.31$	$-0.51$	$4\times {10}^{-3}$

Equations

$$\mathbf{r}=\frac{\Sigma F_{\mathrm{i}}{\mathbf{r}}_{\mathrm{i}}}{\Sigma F_{\mathrm{i}}}$$

$$S(t)=\{\begin{array}{ll}{S}_{\max }{e}^{(t-t_{\mathrm{peak}})/({\tau }_{\uparrow })},& t<t_{\mathrm{peak}},\\ S_{\max }{e}^{(t_{\mathrm{peak}}-t)/({\tau }_{\downarrow })},& t>t_{\mathrm{peak}},\end{array}$$

$$T_{\mathrm{b}}=\frac{2\ln 2}{\pi k_{\mathrm{B}}}\frac{(1+z)S{\lambda }^2}{{\theta }_{\min }{\theta }_{\mathrm{maj}}},$$

$${\Omega }_{r\nu }=4.85\cdot 10^{-9}\frac{\Delta r_{1,2\mathrm{mas}}{D}_{\mathrm{L}}}{(1+z{)}^2}\frac{{\nu }_1^{1/k_{\mathrm{r}}}{\nu }_2^{1/k_{\mathrm{r}}}}{{\nu }_2^{1/k_{\mathrm{r}}}-{\nu }_1^{1/k_{\mathrm{r}}}}$$

$$r_{7\mathrm{mm}}={\Omega }_{r\nu }{\left({\nu }^{1/k_r}\sin \theta \right)}^{-1}.$$

$$B_1=0.014{\left(\frac{{\Omega }_{r\nu }^{3k_r}(1+z{)}^3\ln ({\gamma }_{max}/{\gamma }_{min})}{\varphi {\delta }^2{\sin }^{3k_r-1}\theta }\right)}^{0.25}$$

$$N_1=\frac{B_1^2}{8\pi \ln ({\gamma }_{\max }/{\gamma }_{\min })m_{\mathrm{e}}{c}^2}$$

$${B^{\prime}}^{\mathrm{7~{}mm}}_{\mathrm{core}}=B^{\mathrm{7~{}mm}}_{\mathrm{core}}\frac{\nu^{1/{k^{\prime}}_{r}}}{\nu^{1/{k}_{r}}}\frac{\big{(}{\Omega^{\prime}}_{r\nu}/\sin(\theta)\big{)}^{\frac{{k^{\prime}}_{r}}{k_{b}}-m}}{\big{(}{\Omega}_{r\nu}/\sin(\theta)\big{)}^{\frac{{k}_{r}}{k_{b}}-m}}\Bigg{(}\frac{{N^{\prime}}_{1}}{{N}_{1}}\Bigg{)}^{\frac{-1}{k_{b}(3/2-\alpha)}},$$

$$\Delta r=\frac{\delta \Gamma \beta c\Delta t_{7\mathrm{mm}-\gamma }^{\mathrm{obs}}}{1+z}=\frac{{\beta }_{\mathrm{app}}c\Delta t_{7\mathrm{mm}-\gamma }^{\mathrm{obs}}}{(1+z)\sin \theta }$$

Full text

††pubyear: 2017††pagerange: A connection between $\gamma$-ray and parsec-scale radio flares in the blazar 3C 273–References

1 Introduction

A connection between $\gamma$-ray activity in blazars with evolution of their parsec-scale structure and flux-density variations is now claimed by many authors (e.g. Pushkarev et al., 2010; Agudo et al., 2011a; Jorstad et al., 2012) . However it is still debated from where the $\gamma$-ray emission originates: on sub-parsec scales or on a few parsec scales. The former scenario implies that the $\gamma$-ray emission zone is located close to the true beginning of the jet and the black hole (e.g. Tavecchio et al., 2010; Pushkarev et al., 2010; Rani et al., 2014). Arguments for this point are the following: short timescale of variability in $\gamma$-rays implies that the emitting region should be small and cross-correlation between light curves in $\gamma$-rays and radio statistically puts $\gamma$-ray flares before radio flares. $\gamma$-ray spectral breaks, if caused due to pair production, also argue in favor of the former scenario.

The latter scenario claims that $\gamma$-ray photons should inevitably be absorbed if emitted close to the central black hole (e.g. Ghisellini & Madau, 1996; Ghisellini & Tavecchio, 2008). Moreover some $\gamma$-ray flares are reported to coincide with notable changes in flux and polarization of features located many parsecs away from the jet apex (i.e. Agudo et al., 2011b; Jorstad et al., 2011; Schinzel et al., 2012). The inconsistency of rapid variability with this scenario is usually dealt with by introducing a fine structure inside a jet with some parts moving faster or closer to the line of sight and emitting $\gamma$-rays (e.g. Ghisellini & Tavecchio, 2008; Marscher, 2014).

If a $\gamma$-ray event occurs in a part of a relativistic jet that is optically thin for radio emission, then radio band variability is expected to be simultaneous at different radio frequencies and be simultaneous or even precede the $\gamma$-ray variability (e.g. Ghisellini & Tavecchio, 2008; Agudo et al., 2011a; Jorstad et al., 2013). But if a $\gamma$-ray flare arises in a region that is opaque to radio waves, then one should expect a time lag between $\gamma$-ray event and corresponding event at radio frequencies as well as a delay between activity at different radio frequencies (e.g. Pushkarev et al., 2010; Fuhrmann et al., 2014; Ramakrishnan et al., 2015). At the same time, one should take into account that the emitting region at radio frequencies is substantially larger than that at $\gamma$-rays. This may introduce delays in the observed light curves even if flares are ignited co-spatially and simultaneously. To distinguish between the two scenarios we have conducted a series of multi-frequency VLBI observations that covered a 5 month period during and after a major $\gamma$-ray flare in 3C 273 (Abdo et al., 2010b).

3C 273 is one of the most frequently observed and studied sources. It is the first discovered quasar (Schmidt, 1963) and also the first quasar detected in $\gamma$-rays (Swanenburg et al., 1978). 3C 273 was proven to have extremely high values of brightness temperatures in excess of $T_{\mathrm{b}}\approx 10^{13}$ K and features as small as 26 $\mu$as (Kovalev et al., 2016; Johnson et al., 2016). There are already several published papers on the $\gamma$-ray activity of 3C 273 (e.g. Abdo et al., 2010c, b; Rani et al., 2013), some with single-dish radio data and $\gamma$-ray data (e.g. Fuhrmann et al., 2014; Ramakrishnan et al., 2015), and others using radio VLBI observations (e.g. Pushkarev et al., 2010; Jorstad et al., 2012).

In this paper, we aim to use the capabilities of multi-frequency VLBA observations to pinpoint the location of the $\gamma$-ray–emitting zone in the jet of 3C 273. We used measurements of the frequency-dependent apparent core position, light-curve cross-correlation, and component kinematics. Moreover, several by-products are discussed: a gradient of spectral index in the vicinity of the apparent jet beginning, apparent core movement, and we introduce light-curve-based method for jet speed estimation.

We adopt the $\Lambda$CDM cosmology with the following parameters: $H_{\mathrm{0}}=71$ $\mathrm{km\ s^{-1}}$, $\Omega_{\mathrm{m}}=0.27$, $\Omega_{\mathrm{\Lambda}}=0.73$ (Komatsu et al., 2009). 3C 273 is at a redshift of $z=0.1583$ (Strauss et al., 1992), which corresponds to $2.7\,\mathrm{pc}\,\mathrm{mas^{-1}}$ linear scale.

2 Observations and data processing

2.1 Gamma-rays

Fermi/LAT (Large Area Telescope, Atwood et al., 2009) is a great facility to study temporal variability of blazars because of its sensitivity and all-sky coverage. Fermi revealed that many of the blazars are powerful and variable in $\gamma$-rays.

For the analysis presented here, we used Fermi/LAT Pass 8 data to produce the 100 MeV - 300 GeV $\gamma$-ray light curve using the ScienceTools version v10r0p5. In the event selection we followed the LAT team recommendations for the Pass 8 data111http://fermi.gsfc.nasa.gov/ssc/data/analysis/

documentation/Pass8_usage.html. We modeled a $20^{\circ}$ region around 3C 273 using the instrument response function P8R2_SOURCE_V6, Galactic diffuse model “gll_iem_v06.fits”, and isotropic background model “iso_P8R2_SOURCE_V6_v06.txt”. The light curves were binned using the adaptive binning method by Lott et al. (2012), with estimated 20 $\%$ statistical flux uncertainty in each bin. The flux in each bin was then estimated using the unbinned likelihood analysis and the tool gtlike. All sources within $20^{\circ}$ of 3C 273 that are listed in the 3FGL catalog (Acero et al., 2015) were included in the likelihood model. In order to ensure the convergence of the fits in the weekly bins, we froze all the spectral parameters of all sources (including 3C 273) to the values reported in 3FGL. While in most cases the variability in the spectral index is modest compared to flux variations (Abdo et al., 2010a), we note that in bright sources, such as 3C 273, the spectral index can also vary during flares (e.g. Rani et al., 2013). Because we are not modeling the flaring in detail and only use the light curve as a guidance to activity periods in $\gamma$-rays, this will not affect any of our conclusions.

2.2 Radio

After the $\gamma$-ray flux rose by a factor of 3 above its quiescent level we started our multiwavelength VLBA observations (proposal code S2087A). We conducted four observations which covered a five month long period (see Table 1). The $\gamma$-ray flux rose even more throughout our observations and reached $\approx 20$ times its average level in late September 2009 (Abdo et al., 2010b). Thus our multiwavelength observations covered different stages of the $\gamma$-ray flaring activity.

We observed in 5 bands (C, X, U, K, Q), with four IF channels per band. To equalize sensitivity and extend frequency coverage, the C and X bands were divided into 2 sub-bands thus resulting in 7 frequency bands (denoted further as stated in the Table 2 ). The data from each of these bands were then reduced individually. Initial calibration was performed in AIPS (Greisen, 1990), imaging, self-calibration and model-fitting in Difmap (Shepherd, 1997). Following Sokolovsky et al. (2011), we paid special attention to the accuracy of amplitude calibration. After preliminary calibration, we performed overall gain correction per telescope/frequency/polarization/IF. We checked all available auxiliary information to find potential problems. If the gain correction exceeded 10 $\%$, we applied it to the data. Then we ran self-calibration and imaging again on the gain-corrected data. This approach allowed us to estimate amplitude calibration accuracy of 10 $\%$ in K1 and Q1 bands and 5 $\%$ in C1,C2, X1,X2, U1 bands.

For flux calibration purposes, we utilized several AGN monitoring programs: MOJAVE at 2 cm222http://www.physics.purdue.edu/MOJAVE/allsources.html, Boston University blazar monitoring program at 7 mm333http://www.bu.edu/blazars/VLBAproject.html, UMRAO single dish monitoring at 6.3, 3.8 and 2 cm444http://dept.astro.lsa.umich.edu/datasets/umrao.php, VLA calibrator survey at 6, 3.6 1.3 cm and 7 mm555http://www.aoc.nrao.edu/$\sim$smyers/calibration.

To improve our temporal coverage at 7 mm, we used data from the Boston University blazar monitoring group. We took calibrated uv-data and performed further analysis ourselves using the same techniques as for our data. The 7 mm data, used for the kinematics analysis, spanned a period of 5 years from 2008 to 2012. For the brightness temperature analysis we added two more years of data to cover the RadioAstron detection dates. We excluded epochs with insufficient data quality (i.e. with number of antennas less than 9) from our analysis (Table 3).

3 Modelfitting and components parameters

3.1 Model fitting

The first step of the analysis we performed on our VLBA data was model fitting of source structure with a set of two dimensional Gaussian components (Fig. 1). We used circular Gaussian components for all bright features of the jet structure except sometimes the core. In cases when a fitted size of a component approached zero, it was replaced by a $\delta$-function.

3C 273 has a very broad and complex jet and the number of components in models usually exceeds 15 for 7 mm and 20 for longer wavelengths (Table 4). However, cross-identification of the brightest components between epochs is unambiguous and could be used for analysis of kinematics and evolution of the components parameters. Fig. 1 illustrates models for all frequencies at a single multi-frequency epoch. Each circle represents the full width at half maximum (FWHM) of a circular Gaussian component. A full set of models is available in online materials, Fig. 16-18.

For the analysis of the brightness temperature we calculated a resolution limit for each component and used that, if the fitted size of the model component was smaller. The resolution limit was calculated following Kovalev et al. (2005). Noise root mean square (rms), used in the signal-to-noise ratio estimation, was calculated in a $3\times 3$ beam area around the component centroid position.

3.2 Kinematics analysis at 7 mm

We have performed an extensive kinematics analysis at 7 mm searching for possible correlation of flares in $\gamma$-rays with components ejection and/or flux density variations in components. Owing to its better sampling, the shortest wavelength (7 mm) data are chosen for the jet kinematics analysis.

We cross-identified components between epochs to trace their motion. Sometimes for a given component we had to discard several epochs from our analysis when its position became uncertain due to its low brightness or possible interaction or even merging with adjacent components. We used the following criteria to check the quality of component fitted parameters. First, position of a component should not change abruptly. Second, we checked that flux and brightness temperature follow a general trend to decrease with increasing separation from the core, while component size tends to increase.

As a reference point for all kinematics study we choose the VLBI core – the first component at the north-eastern part of the jet structure. There were already published results on the presence of another stationary feature at $0.15$ mas to $0.16$ mas (Savolainen et al., 2006; Jorstad et al., 2005) in the 7 mm jet of 3C 273. Our model fitting analysis revealed two stationary features (ST1 and ST2) besides the core at a distance of $r_{\mathrm{ST1}}=0.14\pm 0.01$ mas and $r_{\mathrm{ST2}}=0.30\pm 0.01$ mas from the core. We found no components between either core and ST1 or ST1 and ST2, hence we tracked the motion of the components after passing the most downstream stationary feature ST2.

Our model fitting results disagree with ones presented by Jorstad et al. (2012, J12), since we used more components in our models. In particular, their core component A0 is likely to be a sum of two innermost components in our models: the core and ST1; a short-living moving I4 from J12, which reported to be the brightest, according to light curve could be our stationary feature ST2. This affects kinematics results and cross-identification of components with the $\gamma$-ray events.

In the kinematic analysis, we preferred to use a linear function to fit the separation of a component from the core versus time. However there were 2 components that have shown significant deviation from a linear fit and thus were fitted with a parabola. We used a parameterization introduced by Homan et al. (2001) to easily derive angular speed and acceleration.

There is a known issue with parabolic fit that the accuracy of an extrapolation towards the core in such case would not be of a high precision. This complicates the calculation of the ejection (from the core) epoch.

We have tested two approaches. First method is just to take the first few measurements of a component separation from the core and fit them with a straight line. Second is to use a tangent line to a parabolic fit at the position of the first data point. Both methods show consistent results within errors for both components that show accelerated motion.

Another difficulty arises when the trajectory of a component in R.A.–Dec plane does not cross the core, which is the case for the component c1. The minimal separation does not exceed $0.1$ mas, thus we have taken the epoch of minimal separation from the core as an estimate of ejection epoch.

We have found 6 new components during the observed period, labeled c1 to c6 according to the order of their appearance on the images. Three of them passed through the core during the period of $\gamma$-ray activity, MJD 55100 to 55300. For the last detected new component c6 the fitted parameters are not well constrained, therefore we have not used it in the analysis. Components c1 and c2 exhibited accelerated motion, whilst for components c3, c4, and c5 acceleration was not significant, see Fig. 2. Trajectories of all new components except c1 pass the core within the errors. So for the first component we used the moment of minimal separation from the core position as an estimate of the ejection epoch; for the second component we fit the innermost data points, spanning 400 days (11 measurements), with a straight line and use it for estimating the ejection epoch; for others we used simply linear fits. The derived kinematics parameters for all newborn components are summarized in Table. 5.

3.3 Cluster kinematics

Due to the great complexity of the jet in 3C 273, the same region of the extended structure could be model fitted differently at different epochs subject to the image sensitivity, uv-coverage, and slight physical changes in the jet structure. This increases the scatter of the kinematic data for each component, which is not well isolated from adjacent components. Therefore, to improve the accuracy of kinematic measurements, one could study kinematics of a cluster of components that describe a complex region of the jet structure, if such region is isolated from the adjacent parts of the jet.

We found two well isolated clusters, K1 and K2, in the structure of the 3C 273 jet at 7 mm during the flaring period MJD $\approx 55000-55400$, marked in Fig. 1. First, we derived the position of the centroid of all components in the cluster. Coordinates were calculated as

[TABLE]

where $F_{\mathrm{i}}$ is flux density, $\bm{r_{\rm i}}$ - R.A. and Dec coordinates of every single component in the cluster relative to the core at this epoch. Summation was made over all components in the cluster. Depending on the observing epoch, each cluster contained 2 to 6 components. We traced their position over the period and found that both clusters move with constant speed $\mu_{\mathrm{K1}}=0.92\pm 0.01$ $\rm mas~{}yr^{-1}$ and $\mu_{\mathrm{K2}}=0.89\pm 0.02$ $\rm mas~{}yr^{-1}$. Kinematics plot for both clusters is presented in Fig. 3. To search for fine effects we subtracted a linear fit from the kinematic data and scrutinized the residuals.

Fig. 3 shows the difference between the measured separation of the clusters from the core and the best linear fit for both clusters. We found that this difference for both clusters is correlated (chance probability $p=0.008$) and have contemporaneous minima. It should be noted that these clusters are $2$ mas (50 pc) away from each other, so due to causality arguments we reject the idea of their connection. This leads to the conclusion that these correlated residuals reflect not the movement of the clusters but movement of the reference point i.e. the core.

Estimation of uncertainties of the cluster centroid position is a bit tricky. One may not use formulas suitable for isolated components because this is definitely not the case. Instead we estimated uncertainties as root mean square (rms) of the residuals of the linear fit to cluster kinematics. It should be noted that these estimates are upper limits, because measured residuals, as we know, are not randomly distributed.

We have detected a steady movement of the 7 mm core downstream the jet by $\approx 0.17$ mas in 5 month period which corresponds to $\Delta r=4.4$ pc along the jet for the adopted viewing angle $\theta=6^{\circ}$.

4 CROSS-CORRELATION of light curves

Light curves of the 7 mm core and two stationary features along with the $\gamma$-ray light curve are plotted in Fig. 4. Five most prominent $\gamma$-ray flares are labeled $A_{\gamma}$, $B_{\gamma}$, $C_{\gamma}$, $D_{\gamma}$, $E_{\gamma}$ following naming in Rani et al. (2013). Two flaring periods that occurred before $A_{\gamma}$ are labeled $Y_{\gamma}$ (several peaks around MJD 54750) and $Z_{\gamma}$ (around MJD 54950). We aimed to cross-identify flares in $\gamma$-rays and in the 7 mm radio band. Almost the whole period of $\gamma$-ray activity is covered by flares in the 7 mm radio band, which makes one-to-one cross-identification challenging. For further discussion we labeled the 7 mm core flares as $Y_{7}$ (peaked at MJD 54855) and $B_{7}$ (peaked at MJD 55206).

First, we compared the $\gamma$-ray light curve with the 7 mm light curves of the core and both stationary components. These components are usually the brightest among all components of 3C 273 at this wavelength. Each of these components underwent 2 major flares with peak flux densities of $4$ to $7$ Jy. It is worth to note that none of these 7 mm radio flares coincided with the major flare $B_{\gamma}$. The only two peak-to-peak coincidences were: $B_{7}$ in the core coincided with $D_{\gamma}$, and second flare in ST2 ($MJD=55261$) – with $E_{\gamma}$.

Second, we ensured that no other moving component except the newborn components shows any notable flares during the period analyzed. Generally the flux density of all moving components tended to decrease monotonically with distance. Only the brightest ones exhibited flux density of an order of 1 Jy.

Third, we analyzed light curves of new components c1 to c6. Component c1 had its flux density maximum coinciding with the $\gamma$-ray flare $B_{\gamma}$. However, this maximum was the first data point for this component after separation from ST2 (it was already $0.6$ mas from the core and was $0.15$ mas across) and its flux density was $\approx 2.5$ Jy, lower than peaks of the core or stationary features.

4.1 Gamma-ray to 7 mm core time lag

We have analyzed the cross-correlation of the $\gamma$-ray and 7 mm radio core light curves using the discrete correlation function (DCF, Edelson & Krolik, 1988). DCF analysis curves are presented in Fig. 6.

The statistical significance of the correlation was estimated by simulation of artificial light curves. In case of $\gamma$-ray data we created a set of simulated light curves with the same power spectral density (PSD) and probability distribution functions (PDF) as the original data (Emmanoulopoulos et al., 2013; Connolly, 2015).

The 7 mm data contain too few measurements to obtain good PSD and PDF distributions so we applied another approach in this case. For the 7 mm data, each simulated light curve was made from the original data as follows: in an array ($t_{i}$, $S_{i}$) of time and flux density we took a random number $k$ ($1<k\leq N$, where $N$ is the total number of 7 mm data points) of flux density measurements from the end of the array and placed them at the beginning, thus transforming initial array $(t_{1},t_{2}\ldots t_{N};S_{1},S_{2}\ldots S_{N})$ into $(t_{1},t_{2}\ldots t_{N};S_{k},S_{k+1}\ldots S_{N},S_{1},S_{2}\ldots S_{k-1})$.

We performed 1000 simulations calculating DCF between an artificial $\gamma$-ray and 7 mm core light curves. This resulted in a distribution of 1000 correlation functions of two unrelated signals. We then calculated a $95\,\%$ confidence level as a value such that a probability of DCF to be higher by chance than this value is $5\,\%$. Such a one-tailed test can also be applied to the negative correlation coefficients (note its absolute value is much lower due to the properties of the PDF distribution of the $\gamma$-ray data).

Based on these tests we find that the correlation between $\gamma$-ray and 7 mm light curves is significant for a time lag in a range 90 to 130 days. A more accurate estimate of uncertainties could be obtained using this information along with Monte Carlo simulations of the data. Now we perform 1000 evaluations of DCF using flux randomization (FR) and random subsample selection (RSS) methods (Peterson et al., 1998). First, at each stage the flux density of every point is taken as a random value normally distributed around the measured value with dispersion corresponding to its error (FR). Second, a random subsample of the same length as the data is selected from the original light curve (RSS – modified bootstrap). The former step (FR) accounts for uncertainties due to flux measurements, while the latter (RSS) allows one to eliminate the contribution of outliers to the correlation results. Finally we obtain a cross-correlation peak distribution (Fig. 7) for time lags in the range $80-140$ days, from which the time lag and its error were estimated as the mean and standard deviation. We obtain $\Delta t_{\mathrm{7~{}mm-\gamma}}^{\mathrm{obs}}=112\pm 9$ days lag with $\gamma$-rays leading the 7 mm core. This value is close to the peak-to-peak lag of $115$ days between $B_{\gamma}$ and $B_{7}$.

For our calculations we used a DCF time bin $\delta t=21$ days, which is the mean time spacing of $\gamma$-ray and radio light curves. We note, however, that the correlation results do not exhibit any dependence on the choice of the $\delta t$ value in the range $5\leq\delta t\leq 50$.

We compared the 7 mm to $\gamma$-ray lag with that derived with single-dish observations at 8 mm by Ramakrishnan et al. (2015) who found a lag of 160 days and at 7 mm by Chidiac et al. (2016) who found a lag of 160 days. These values agree well with our 112 day $\gamma$-ray to 7 mm core lag, considering that the single-dish observations could not disentangle flares of comparable amplitudes in the core and stationary components. We checked this by estimating the time lag from the peak of the $\gamma$-ray flare $B_{\mathrm{\gamma}}$ to the peak of the aggregate core+ST1+ST2 light curve and found a lag of 170 days.

4.2 Subsequent flares in core–ST1–ST2

Fig. 4 shows the light curves of $\gamma$-rays and the core and two stationary components in the 7 mm band. It is clear that a radio flare first occurs in the core and then subsequently in the two stationary components. Taking components’ proper motion to be $\mu=0.8$ $\mathrm{mas\,yr^{-1}}$ – typical for newborn components at that period – and taking the measured distances between stationary features, one would expect the peak of the ST2 light curve to occur at MJD $~{}55400$, significantly later than observed.

This inspired us to calculate speed of a moving disturbance that could produce such flares. Cross-correlation of light curves is inappropriate in cases of such scarcely sampled and small data samples. So to get the moment $t_{\mathrm{peak}}$ of the peak we fitted flares with an exponential rise and decay

[TABLE]

where $\tau_{\uparrow}$ and $\tau_{\downarrow}$ are variability time scales for the rise and the decay of a flare, $S_{\mathrm{max}}$ is the peak flux density, similarly to Valtaoja et al. (1999), but with no assumed relation between $\tau_{\uparrow}$ and $\tau_{\downarrow}$.

This worked well for the core and ST1, but ST2 has even worse sampling, so we chose the position of a maximum in the middle between 2 points with the highest flux density with uncertainty as large as to let the true maximum fall between adjacent points with $95\,\%$ probability.

This resulted in a proper motion $\mu=(1.21\pm 0.06)$ $\mathrm{mas\,yr^{-1}}$ or $\beta=(12.3\pm 0.6)\,c$ which is not extraordinary for 3C 273 (e.g. Savolainen et al., 2006; Lister et al., 2009; Jorstad et al., 2012) but is 1.5 times higher than the speed measured downstream for the components that may be associated with these flares.

This approach is important in two senses. First, this is an alternative way to estimate the jet plasma flow speed. Second, speed is measured very close to the apparent jet beginning where direct component kinematics usually gives poor quality due to the strong dependence between adjacent components parameters.

The difference in speed estimated near the apparent core position using the flux density light curves of the 7 mm core and stationary features and that calculated using direct position measurements farther downstream might be explained by either different jet geometry in the core region, or some kind of a deceleration, or the difference between the flow speed and the pattern speed.

4.3 Brightness temperature of the core and ST1 components

With measured flux density and size of an elliptical Gaussian component one could estimate its brightness temperature in the rest frame of the source as follows (e.g., Kovalev et al., 2005):

[TABLE]

where $k_{\mathrm{B}}$ is the Boltzmann constant, $\lambda$ is the wavelength, $S$ is flux density of the component, and $\theta_{\mathrm{min}}$ and $\theta_{\mathrm{maj}}$ represent either the measured FWHMs of the component or the resolution limits along the component’s axes if the component is not resolved.

We found that the brightness temperature for the most compact and prominent features, the core and ST1, in the parsec-scale structure of 3C 273 could vary by several orders of magnitude, from $10^{10}$ K to $10^{13}$ K in the source frame. Throughout our multi-frequency observations the brightness temperature was of the order of several times $10^{12}$ K at all observed frequencies.

We used the long-term monitoring data at 7 mm to trace the evolution of the brightness temperature on the timescale of years. Figure 8 shows the evolution of the brightness temperature for the most prominent features at the apparent beginning of the jet at 7 mm, the core and ST1. Since both these components were unresolved at some epochs, we used the upper limits on their size for those epochs, following Kovalev et al. (2005), to estimate the lower limit on the brightness temperature.

The observed brightness temperature of the core in the source frame varies from $10^{10}$ K to $>10^{13}$ K (Fig. 8), reaching the highest values simultaneously with the peak of the flare $B_{7}$ on 2010 Jan 10. The highest measured values and lower limits agree with that derived by RadioAstron (Kovalev et al., 2016; Johnson et al., 2016). Meanwhile, the Doppler-boosted equipartition brightness temperature (Readhead, 1994, eq. 4b) in the rest frame of 3C 273 is $T^{\prime}_{\mathrm{eq}}\approx 4\times 10^{11}$ K, for Doppler factor $\delta=6$ (Savolainen et al., 2006). According to, e.g., Readhead (1994) the excess of a brightness temperature over the equipartition value could be ascribed to the dominance of the particle energy density over that of the magnetic field. Measured values of the brightness temperature of the order of $10^{13}$ K require that the conditions at the apparent 7 mm core should significantly deviate from equipartition (see also Kovalev et al. (2016)). During the flare the particle energy density should be greater than the magnetic field energy density. We note, that the huge range of the $T_{\mathrm{b}}$ variations (Fig. 8) could not be explained by the variations of the Doppler factor and is most probably due to variations in the particle density and the magnetic field strength.

5 Frequency-dependent apparent core position

We have taken advantage of our multi-frequency dataset and performed analysis of the frequency dependent apparent core position, the so-called “core shift” effect (e.g. Marcaide & Shapiro, 1984; Lobanov, 1998; Kovalev et al., 2008; O’Sullivan & Gabuzda, 2009; Sokolovsky et al., 2011; Kravchenko et al., 2016). At each observed frequency the core is believed to be a surface of optical depth $\tau=1$ and its location depends on frequency (e.g., inversely, in model of Blandford & Königl, 1979).

With the assumption that optically-thin synchrotron radio emission from the jet originates in the same volume at different frequencies, the position of such optically-thin features should not be affected by opacity effects and thus could be used as a reference point to derive a frequency dependent shift of optically thick parts of the jet. We used two methods to derive the shifts between maps: a cross-correlation of large optically-thin parts of the jet structure and alignment of distinct model components (e.g. Pushkarev et al., 2012; Kutkin et al., 2014).

In the former method a shift between two maps is derived at first using two dimensional cross-correlation (2DCC) on the optically-thin portions of the jet structure. Then a core shift is calculated using the map shift and core coordinates in each map. The latter method uses distinct optically-thin model components at different frequencies as reference points because they should not experience a frequency-dependent shift. A component used in this method ought to be optically thin and preferably isolated to avoid contamination with adjacent components. With our model fitting strategy, all strong components in the jet structure are usually not isolated and hence their position is affected by adjacent components. Therefore, we just checked that optically-thin component alignment gave results consistent with the 2DCC method, but for our subsequent analysis we used only results obtained with the 2DCC method.

For 2DCC analysis we used maps with identical map size, pixel size, beam parameters, and range of sampled baselines (uv-range). Therefore we had to cut off data from long baselines at higher frequency and data from short baselines at lower frequency to match the uv-range. We cross-correlated all pairs of maps that exhibited sufficient source structure, not only at adjacent frequencies. Then for each pair of frequencies we calculated the avareage map shift $\langle\Delta r_{\mathrm{map}}\rangle$ and its rms $\sigma_{\Delta r_{\mathrm{map}}}$.

We inspected spectral-index maps with the map shift applied in order to check the fidelity of the derived shift. In all cases the resulting spectral-index map looked smoother and exhibited an increased spectral index $\alpha$ ($S\propto\nu^{\alpha}$) near the core position indicating that the jet is optically thick in this region.

The core shift comprises the map shift $\langle\Delta r_{\mathrm{map}}\rangle$ and the position of the cores. The total error of the core shift derived using 2DCC consists of pixel size, core position error and $\sigma_{\Delta r_{\mathrm{map}}}$. To decrease the impact of pixel size on the resulting errors in core shift we reduced the pixel size to values small as $0.01$ mas which is about $1/20$th of the beam FWHM at 7 mm and a much lower value for longer wavelengths. Areas on the maps for performing the cross-correlation were chosen manually, so we ran 2DCC several times for each frequency pair per epoch to reduce the effect of a human-driven choice. Errors in the position of the core, $\sigma_{\mathrm{core}}$, were estimated in the image plane following Fomalont (1999).

5.1 Derived physical parameters of the jet

We fitted our core-shift measurements with a power law $\Delta r=a\nu^{-\frac{1}{k_{\mathrm{r}}}}+b$ as shown in Fig. 9. A model of the jet emission (Blandford & Königl, 1979) assuming equipartition between magnetic field and kinetic energy and a conical shape of the jet implies $k_{\mathrm{r}}\approx 1$ (Lobanov, 1998). The value of $k_{r}=((3-2\alpha)m+2n-2)/(5-2\alpha)$ depends on the optically thin spectral index $\alpha$ and power indices of the B-field and particle density decay along the jet, $B(r)=B_{1}r^{-m}$ and $N(r)=N_{1}r^{-n}$ respectively, where $B_{1}$ and $N_{1}$ are measured at 1 pc from the jet beginning.

In the model of Blandford & Königl (1979) the apparent jet base is located farther downstream at lower frequencies due to frequency dependent synchrotron opacity. Thus, if a flare is caused by a disturbance that moves downstream and carries increased values of the particle density or magnetic field, then our single-epoch observation could catch enhanced values of $N$ and $B$ in the apparent cores at higher frequencies and unchanged ones at lower frequencies. This would introduce a steeper gradient of $N(r)$ or $B(r)$ in the region, where the core shift measurement is performed, which would in turn increase the $k_{r}$, as observed (Table 6).

Taking the measured values of the core shift $\Delta r_{\mathrm{1,2\,mas}}$ and $k_{r}$ we calculated $\Omega_{r\nu}$ (Lobanov, 1998; Pushkarev et al., 2012) as follows:

[TABLE]

where $D_{\mathrm{L}}$ is the luminosity distance. A deprojected distance between the VLBI core at 7 mm and the jet nozzle was estimated as:

[TABLE]

The core ($\tau=1$ surface, Blandford & Königl, 1979) position depends on both particle density and magnetic field strength (e.g Lobanov, 1998). Most of researchers assume equipartition between these two quantities (e.g., Lobanov, 1998; Hirotani, 2005; Kutkin et al., 2014). Even though Zdziarski et al. (2015) suggested a method to calculate the jet parameters not assuming equipartition, the accuracy of typical core-shift measurements, including ours, is not sufficient to apply it. Meanwhile, the measured values of the brightness temperature $T_{\mathrm{b}}>10^{13}$ K (Fig. 8) imply that the particles dominate in the total core energy budget. To estimate parameters of the jet base, we first derive the magnetic field and particle density assuming equipartition and then recalculate them for a more plausible set of parameters.

For all four epochs we derive the equipartition magnetic field strength at the distance of 1 pc from the central engine following Kutkin et al. (2014) as

[TABLE]

where $\gamma_{max}/\gamma_{min}$ ratio is adopted to be $10^{3}$, the Doppler factor $\delta=6$, the viewing angle $\theta=6^{\circ}$ (see discussion in sect. 7.1) and the corresponding intrinsic opening angle $\phi=1.1^{\circ}$ (Pushkarev et al., 2009). To estimate the magnetic field strength at the location of the VLBI core we have to assume some reasonable values of power law indices $n$ and $m$. We adopted the equipartition values $n=2$, $m=1$ (Blandford & Königl, 1979).

Equipartition implies the following relation between magnetic field and particle density:

[TABLE]

where $m_{\mathrm{e}}$ is electron mass, and $\gamma_{\mathrm{min}}$ and $\gamma_{\mathrm{max}}$ are the electron energy cutoffs. Equipartition values of $B_{1}$, $B^{\mathrm{7~{}mm}}_{\mathrm{core}}$, and $N_{1}$ are presented in Table 6. The values of electron density at 1 pc are comparable with those presented by Lobanov (1998). However, since we have detected the brightness temperature up to $>10^{13}$ K, the real particle density should be greater than that estimated assuming equipartition. We estimate below the magnetic field strength assuming several increased values of the particle density.

According to Blandford & Königl (1979) and Lobanov (1998) the magnetic field in the 7 mm apparent jet base depends on $\Omega_{r\nu}$ and $N_{1}$. We consider the epoch 2009 Aug 28 as a reference one, because physical conditions in the 7 mm core were likely close to equipartition, as deduced from the brightness temperature (Table 6), and the estimated B-field should be close to its real value. We took the measured values of $k_{r}$ and $\Omega_{r\nu}$ and adopted $\alpha=-0.5$ (hence $k_{b}=2/3$) and $m=1$ to calculate the ${B^{\prime}}^{\mathrm{7~{}mm}}_{\mathrm{core}}$ for other three multi-frequency epochs:

[TABLE]

where primed values refer to the epoch of interest and unprimed refer to 2009 Aug 28. This formula takes into account both the change of the position of the 7 mm core and the change of the particle density and does not assume equipartition. Table 7 presents the magnetic field value at the position of the 7 mm core for different epochs and for 10, 100, and 1000 times increase of the particle density relative to 2009 Aug 28.

We compared our estimates of the apparent 7 mm core B-field with that reported by Savolainen et al. (2008) who have estimated \log_{10}\Bigl{(}B^{\mathrm{7~{}mm}}_{\mathrm{core}}\Bigr{)}=0.3\pm 0.6 G on the basis of synchrotron self-absorbed spectrum turn-over measurements. We estimate the observed brightness temperature of the 7 mm core using parameters provided by Savolainen et al. (2008) to be $T_{\mathrm{b}}\approx 2\times 10^{11}$ K, close to its equipartition value, relativistically boosted. Despite the time difference of 6 years and different model-dependent approaches, $B^{\mathrm{7~{}mm}}_{\mathrm{core}}$ of Savolainen et al. (2008) is comparable to that estimated for our epoch 2009 Aug 28, when conditions in the apparent 7 mm core were also close to equipartition. For other three epochs, exhibiting particle density dominance, according to high $T_{\mathrm{b}}$, the magnetic field in the 7 mm core should have been lower (Table 7).

6 SPECTRAL INDEX analysis

After aligning maps we analyzed the two-point spectral index distribution (see Fig. 10) and its evolution during the flaring period. Spectral-index maps were prepared in the same manner as for the core-shift analysis (see Sect. 5) with a $5\sigma$ cutoff applied to the total-intensity maps. Examples of spectral-index maps made between 43 GHz and 24 GHz are shown in Fig. 10.

Following Hovatta et al. (2014), we performed tests to check the impact of different uv-coverage (visibility space coverage) on the spectral index maps. For each frequency pair we produced artificial data in the following manner. We took a model of brightness distribution (result of the CLEAN algorithm) and simulated data on the uv-coverages of both frequencies using the AIPS task UVMOD. The resulting data were then imaged with parameters identical to those used for real data. Since the brightness distribution the same at both frequencies, spectral index map based on these simulated data should be zero and any discrepancies are due to uneven uv-coverage. This was done separately for high- and low-frequency CLEAN models. We considered a pixel in the spectral-index map as affected by uv-coverage effects if for any of the input CLEAN models the resulting spectral-index value is not zero within the errors. Deviation of this artificial spectral index value from zero tends to increase with distance from the image center. According to our tests uv-coverage-related effects could be neglected up to certain distances from the map center along the jet. For 43 GHz–24 GHz this distance is $4$ mas–$5$ mas, depending on epoch, with one exception, epoch 2010 Jan 26, which was affected strongly because there were only 9 antennas and the uv-coverage hence was more sparse. For 24 GHz–15 GHz it is $5$ mas–$6$ mas, for 15 GHz–8 GHz – up to $25$ mas, for 8 GHz–5 GHz – up to $30$ mas.

6.1 Core spectral index

Measuring core spectral index using the flux density of model components at different frequencies might be misleading because we believe the cores at different frequencies are physically different regions in the jet due to opacity. We used a conservative way to estimate the spectral index at the core position. On a two-frequency spectral-index map we derived spectral index value as an averaged value of the 3x3 pixels area around the center of the VLBI core at a lower frequency.

We studied the evolution of spectral index between adjacent frequencies during the flare at the apparent jet base. Fig. 11 shows that during all stages of the flare $B_{7}$, the spectrum in the core at all frequencies was flat or inverted.As expected, higher frequencies respond to the changing physical conditions more swiftly than lower ones. The spectral index between 43 GHz and 24 GHz exhibits the most remarkable rise from $\alpha_{\mathrm{43-24}}\approx 0$ to $\alpha_{\mathrm{43-24}}\approx 1.5$ following the rapidly rising flux density of the 7 mm core. The flux density of the 1.3 cm core also rose during the flare and was reflected in the 24 GHz to 15 GHz spectral index, but not so quickly and with some delay. The spectral index at lower frequencies was almost unchanged during the period considered. This may also indicate that a flare in the cm-band apparent jet base is due to a disturbance carrying larger $N_{\mathrm{e}}$ that propagates along the jet, passing regions with $\tau\approx 1$ at subsequently lower frequencies. Our measurements did not cover the rise of spectral indexes $\alpha_{\mathrm{15-8}}$ and $\alpha_{\mathrm{8-5}}$ associated with the flare.

6.2 Ridge-line spectral index

We analyzed the distribution of spectral index along the ridge line of the jet. For each two-point spectral index map, the ridge line was calculated at a lower frequency, starting at the position of the VLBI core in a manner similar to that presented by Hovatta et al. (2014). A plot of the spectral index $\alpha_{\mathrm{43-24}}$ vs distance along the ridge line for 2009 Oct 25 is shown in Fig. 12. Locations where the ridge line crosses a component are overplotted as a shaded area to compare the component’s convolved size with the typical scale of spectral-index variations.

Average values of spectral index along the ridge line are presented in Table 8. We restricted the range of distances to calculate the spectral index values to the following: $0.6$ mas to $5$ mas for $\alpha_{\mathrm{43-24}}$, $0.8$ mas to $6$ mas for $\alpha_{\mathrm{24-15}}$, $1$ mas to $25$ mas for $\alpha_{\mathrm{15-8}}$, and $1.7$ mas to $30$ mas for $\alpha_{\mathrm{8-5}}$. The lower bound is set by the edge of the convolved core (and stationary features for 43 GHz), while the higher bound is set by uv-coverage effects.

Because of our model-fitting strategy (many model components) and restricted distances along the ridge line, according to our tests, ridge lines for all frequencies are almost entirely covered by components. This prevents us from detecting significant deviations of spectral-index values at the positions of components and between them.

We detect a decrease of spectral index, $\Delta\alpha$, with respect to the distance along the jet for $\alpha_{\mathrm{15-8}}$ and $\alpha_{\mathrm{8-5}}$ spectral indices on the common range of ridge-line distances $1.7$ mas–$25$ mas. Since uv-coverage effects are taken into account we simply fit the spectral index along the jet with a linear function over a range of distances mentioned above. The average slope of the spectral index along the ridge line is $0.036$ $\mathrm{mas}^{-1}$, which corresponds to decrease of $10^{-3}$ per deprojected parsec. The value of slope does not differ significantly between $\alpha_{\mathrm{15-8}}$ and $\alpha_{\mathrm{8-5}}$.

6.3 Spectral index gradient in the transition zone 43 GHz–24 GHz

We checked whether a transition zone between optically thick ($\tau\sim 1$) and thin ($\tau\ll 1$) emission in the core region could be resolved with our VLBA observations. We analyzed slices of spectral index (43 GHz–24 GHz and 15 GHz–8.4 GHz) in the direction of the jet. Errors of spectral index values along the slice are taken from the spectral-index map.

Hovatta et al. (2014) report that a transition from inverted spectrum in the optically-thick core to a steep one in the optically-thin jet is smooth not only due to convolution with the beam but is intrinsic to the jet. We took maps of spectral index between 43 GHz and 24 GHz for all epochs and sliced a region of 1 mas around the core position. A slice starts behind the core and extends in the direction of the jet. The distribution of spectral index we observe is a smoothed intrinsic distribution.

There are two ways to estimate the shape of the intrinsic spectral index distribution. First, one may deconvolve the measured distribution to get the intrinsic distribution. Second, one may assume a shape for the intrinsic distribution and convolve it with the beam to compare with the observed distribution. We believe the latter is more appropriate for this kind of study. We considered 2 types of intrinsic distributions: a simple step function and a step with linear transition, hereafter ramp. We set values of the functions at the ends of the slice to maximum and minimum values of measured spectral index. Thus the step function has only a single free parameter, ramp has two. We convolved every function with a Gaussian function describing a projection of the beam onto the slice, which resulted in a beam FWHM $\approx 0.35$ mas. By varying the free parameters we fitted the convolved function to the data. An example of the best fitted functions for epoch 2009 Aug 28 is presented in Fig. 13.

The ramp function with intrinsically extended transition from $\alpha>0$ to $\alpha<0$ always fits better than a simple step function with a sharp intrinsic drop of the spectral index. We used corrected Akaike information criterion (AICc) (Akaike, 1974; Burnham & Anderson, 2002) to test whether ramp fits better because it has more free parameters than step or does it really better describe the data. The difference in AICcs is very convincing, see Table 9, favoring the ramp function. The resulting relative probability is calculated as $exp((w_{\mathrm{step}}-w_{\mathrm{ramp}})/2)$, where $w$ stands for the corresponding AICc value.

We should note that with the given pixel size, adjacent points are not completely independent, which in turn could considerably increase the AICc difference. We have tested a dependence of AICc difference on the number of data points along the slice. We used $\alpha_{43-24}$ and $\alpha_{15-8}$ spectral index data for the epoch 2009 Aug 28. We constructed regular subsamples of the dataset taking every first, every second, every third point etc. A relative probability of 0.01 in $\alpha_{43-24}$ is reached when the number of data points is reduced to values small as 10, which corresponds to distance 0.1 mas (1/3 of a beam) between points. For the $\alpha_{15-8}$ spectral index the ramp model is overwhelmingly more probable regardless of the number of data points. Therefore in our analysis we did not use the whole dataset, but a subsample selected with the intent to leave 5-7 data point per beam FWHM. This resulted in 22 points in 43 GHz–24 GHz and 25 points in 15 GHz–8 GHz .

Another test was performed to check the dependency of the relative model probability on the number of data points removed from beginning of the slice, i.e. from the region upstream the core. This dependence could affect the result if a significant amount of data is removed. As shown in Fig. 14, the dependence is not substantial for $\alpha_{15-8}$ until 0.7 mas (0.6 of a beam FWHM) is deleted, which would place a starting point almost onto the upstream core edge. For $\alpha_{43-24}$ the difference between models vanishes when 0.15 mas (0.45 of a beam FWHM) is deleted. Dependency of the determined width of the transition zone on the amount of data taken from behind the core shows a roughly linear trend in data range when step and ramp models became significantly different.

The minimal estimate of the width of the transition zone for $\alpha_{\mathrm{43-24}}$ is $0.3$ mas which corresponds to 1 pc projected distance and 7 pc deprojected distance. For $\alpha_{\mathrm{43-24}}$ most of the transition zone covers a region between the core and the stationary feature ST2.

For $\alpha_{\mathrm{15-8}}$ the estimated width of the transition zone is $2$ mas, and a significant part of the transition zone is located upstream from the 15 GHz core at all epochs.

7 DISCUSSION

7.1 Location of the $\gamma$-ray emission zone

We considered two effects that might result in a time lag between $\gamma$-ray and radio light curves. First, if both $\gamma$-ray and radio emission is caused by the same perturbation moving downstream from the jet beginning and if a $\gamma$-ray flare is located in a region that is opaque to radio waves then the radio emission is absorbed until the perturbation reaches the partially optically-thin apparent core. Second, even if $\gamma$-rays and radio are ignited simultaneously in the same location, transparent for both radiowaves and $\gamma$-rays, then a delay measured between peaks of the flares might be a consequence of a larger size of the radio-emitting region compared to the that of $\gamma$-rays.

Adopting the former scenario and given the delay between light curves from cross-correlation one may calculate a distance between $\gamma$-ray and radio emitting regions as

[TABLE]

following Pushkarev et al. (2010) and the assumptions discussed by them.

The main uncertainty of the $\Delta r$ arises from the poor knowledge of the viewing angle of the jet. Published estimates range from $\theta_{\mathrm{var}}=10\mathring{.}1$ (Lähteenmäki & Valtaoja, 1999) to $\theta_{\mathrm{var}}=3\mathring{.}3$ (Hovatta et al., 2009) and $\theta_{\mathrm{VLBI}}=6\mathring{.}1\pm 0\mathring{.}8$ Jorstad et al. (2005) and $\theta_{\mathrm{VLBI}}=9\mathring{.}8\pm 2\mathring{.}1$ (Savolainen et al., 2006). According to our $\beta_{\mathrm{app}}$ estimates on the basis of core and stationary features variability (Sect. 4.2) we could derive the maximum value of the viewing angle for $\beta_{\mathrm{app}}=12.3\,c$: $\theta_{\mathrm{max}}=\arcsin(2\beta_{\mathrm{app}}/(1+\beta_{\mathrm{app}}^{2}))=9\mathring{.}3$ which in turn gives a lower limit on $\Delta r_{\mathrm{min}}=6.5$ pc. Estimates of $\theta_{\mathrm{max}}$ based on component speed measurements, $\beta_{\mathrm{app}}=8.4\,c$, give $\theta_{\mathrm{max}}=13\mathring{.}6$ and $\Delta r_{\mathrm{min}}=3.0$ pc. We believe the measurements of Jorstad et al. (2005) are most robust given their excellent temporal coverage, hence, we adopt their $\theta=6\mathring{}$. With the measured delay of $112$ days and mentioned above range of the jet parameters we estimate the deprojected distance between $\gamma$-ray emitting zone and 7 mm VLBI core to be $\Delta r=4-12$ pc. For adopted $\theta=6\mathring{}$ the corresponding $\Delta r=6.5$ pc is taken as an estimate of the distance between the $\gamma$-ray emitting region and the 7 mm core. Since $\gamma$-rays lead radio, the $\gamma$-ray emission zone is located closer to the central engine.

The latter scenario implies that both $\gamma$-ray and radio flares start simultaneously. There might be the following associations: $Y_{\gamma}$ starts simultaneously with $Y_{7}$, and both $A_{\mathrm{\gamma}}$ and $B_{\mathrm{\gamma}}$ are consistent with starting contemporaneously with $B_{\mathrm{7}}$ or with the second flare in ST1, and $D_{\mathrm{\gamma}}$ coincides with the start of the second flare in ST2, see Fig. 4.

To check if the measured time lag between flares $B_{\gamma}$ and $B_{7}$ could be fully explained in the latter scenario, we have compared characteristic rise time of the flare $B_{\mathrm{7}}$ with the light travel time across the core extent. The 7 mm core was unresolved in most observations during the flare $B_{\mathrm{7}}$, so its size was adopted from the resolution limit, which resulted in an average Gaussian FWHM of $a=0.035\pm 0.005$ mas during the flare. For further calculation we used the core size of $1.6a$ in case the core is a uniform optically-thick disk. The expected peak-to-peak delay due to the larger size of the emitting region at 7 mm in the observer’s frame is $\Delta t_{\mathrm{obs}}=208\ \delta^{-1}$ days, where $\delta$ is a Doppler factor. For $\delta=4-9$ (Lähteenmäki & Valtaoja, 1999; Jorstad et al., 2005; Savolainen et al., 2006) this implies $\Delta t_{\mathrm{obs}}=52$ days–$23$ days and for $\delta=D_{\mathrm{var}}=17$ (Hovatta et al., 2009) $\Delta t_{\mathrm{obs}}=12$ days. For adopted $\delta=6$, $\Delta t_{\mathrm{obs}}=34$ days.

We fitted 7 mm flares $Y_{\mathrm{7}}$ and $B_{\mathrm{7}}$ with exponentials, Eq. 2. For the flare $B_{\mathrm{7}}$ this results in $\tau_{\uparrow}=82$ days, $\tau_{\downarrow}=94$ days. Comparison of $\tau_{\uparrow}$, $\tau_{\downarrow}$ with the expected $\Delta t_{\mathrm{obs}}$ shows that none of the reported values of the Doppler factor could fully reconcile the light-travel time across the 7 mm core with the variability time scale of the flare $B_{\mathrm{7}}$. Hence we can not treat this effect as the sole factor determining the peak-to-peak lag between 7 mm and $\gamma$-rays, and thus the $\gamma$-ray emitting region should be located upstream from the 7 mm VLBI core.

Fine structure of the $\gamma$-ray flare in 3C 273 (Abdo et al., 2010b; Rani et al., 2013) reveals $\gamma$-ray flux variations on the time scale of less than 1 day, which impose certain constraints on the size of the emitting region. Assuming that the $\gamma$-ray emitting region spans the whole width of a conical jet, the variations imply that the distance between the jet nozzle and $\gamma$-ray-emitting region is less than $10^{-3}$ pc and thus is also placed close to the true base of the jet.

It is worth noting that during the flare $B_{7}$ the 7 mm core moved downstream the jet. Since the peak of the flare $B_{\mathrm{7}}$ occurred when the 7 mm core was at the most downstream position, the cross-correlation analysis gave the maximal value of $\Delta r$. If the core-shuttle amplitude of $4.4$ pc is taken into account, the distance from the $\gamma$-ray emitting region to the 7 mm core calculated using the light-curve lag and derived from the core-shuttle analysis is $2.1$ pc to $6.5$ pc, depending on the observational epoch. This agrees well within errors with the distance from the 7 mm VLBI core to the jet nozzle, $1.5$ pc to $4.4$ pc, derived from the core-shift analysis. This leads us to the conclusion that the $\gamma$-ray emitting region in 3C 273 is located $2-7$ pc upstream from the 7 mm apparent core, close to the true base of the jet.

7.2 Passage of moving components through the core, ST1, ST2 and relation to $\gamma$-ray flares

Since the jet structure of 3C 273 contains 3 stationary features Core, ST1, and ST2, which could be physical structures like standing shocks, we have performed a cross-identification of all events related to these components with the $\gamma$-ray flares. We calculated epochs when components c3–c5 passed through stationary features. The choice of these components is dictated by their appearance in the region of stationary features contemporaneously with the period of high $\gamma$-ray activity. In total there were 9 events that might be examined, see Fig. 15. Except for the coincidence of c5 ejection with $C_{\gamma}$ there is no robust evidence of any correspondence between high-energy flares and transit of components through standing features.

If, in contrast, a $\gamma$-ray flare was caused by a passing component in the vicinity of the jet apex, then solving Eq. 9 for $\Delta t_{\mathrm{7~{}mm}-\gamma}^{\mathrm{obs}}$ could link ejected components with $\gamma$-ray flares. This requires an assumption of a constant component speed from the beginning of the jet to the 7 mm core. With the measured $\Delta r=2-7$ pc and adopted $\theta=6^{\circ}$, c1 and c2 could have produced $Y_{\gamma}$, c4 – either $A_{\gamma}$ or $B_{\gamma}$, c5 – $Z_{\gamma}$, and c3 is left with no $\gamma$-ray counterpart.

The derived ejection epochs for components c1 and c2 allow one to connect their passage through the core with the 7 mm core flare $Y_{7}$. The ejection epoch of component c4 connects it to the 7 mm core flare $B_{7}$ and therefore component c4 could be linked to the preceding $\gamma$-ray flare $B_{\gamma}$. The component appeared after our multi-frequency observations so we have checked parameters of the component derived from only 7 mm data. The component speed is similar to values derived for other components in this research or to values published before (Savolainen et al., 2006; Lister et al., 2009; Jorstad et al., 2005). We found that this component has maximal brightness temperature $T_{\mathrm{b}}=10^{12}$ K in the source frame and a maximal peak flux density $S=6.4$ Jy, the highest among all of the new components ejected between 2008 and 2011.

7.3 Core shuttle

Possible movement of the apparent core, inferred from VLBI kinematics studies, was mentioned previously (e.g. Kovalev et al., 2008; Lister et al., 2009; Hodgson et al., 2017). A recent dedicated astrometric experiment (Niinuma et al., 2015) has proven that the core movement along the jet is significant in Mrk 421 and occurs after a strong X-ray flare. There are two effects that could lead to the shuttle of the core. First, if the core position is constant and a new component emerges from the core but is still unresolved from it, one could observe that the core “moves” downstream until the new feature becomes resolved from it and then the core jumps back. Second, if physical conditions in the jet change, the core should shuttle down- and upstream. This is valid for both interpretations of the mm-core: a $\tau\approx 1$ (e.g. Lobanov, 1998) surface or a standing shock (Daly & Marscher, 1988).

Using the kinematics of the two clusters of components, we found clear evidence that the apparent core of 3C 273 at 7 mm moved downstream during the onset of the powerful flare $B_{\mathrm{7}}$ by $4.4$ pc. The core position was anti-correlated with the core flux density. The Pearson correlation coefficient is $r=-0.57$ with a chance probability $p=10^{-3}$ for the cluster K1, and $r=-0.51$, $p=4\times 10^{-3}$ for K2. To check if the core shuttle was induced by the flare we performed an additional test by considering independently data during the flare $B_{\mathrm{7}}$ (MJD 55059-55361) and outside the flare (MJD 55003-55038, 55409-55763). For both clusters there is no significant correlation left for the data outside the flare. Meanwhile during the flare $B_{7}$ the negative correlation is even more significant. The correlation parameters are summarized in Table 10.

A possible explanation of such a prominent negative correlation during the flare is that the core as a $\tau\approx 1$ surface moves along the jet while physical conditions change. A flare in the core could be caused for example by an increased particle density $N_{\mathrm{e}}$ or increased magnetic field $B$ which also affect the opacity. Higher values of $N_{\mathrm{e}}$ or $B$ act so that $\tau=1$ occurs farther downstream. The core shuttle together with the frequency dependent core shift (Kovalev et al., 2017) could affect the accuracy of the derived offset between radio VLBI and optical Gaia reference frames.

8 Conclusions

Correlated variability of radio and $\gamma$-ray emission, kinematics of a new emerged component and short time scale of $\gamma$-ray variability together with core-shift measurement lead us to the conclusion that the $\gamma$-ray emission site of the most powerful flare $B_{\gamma}$ is located in the nearest vicinity of the true base of the jet in 3C 273, $2$ pc to $7$ pc upstream from the 7 mm VLBI core. This strongest $\gamma$-ray flare might be associated with the brightest moving component c4. 2. 2.

Kinematic analysis shows that the 7 mm core shuttles along the jet during the 7 mm core flare $B_{7}$. Core flux density anti-correlates with the core position $r_{7}$. The total range of the core displacement is $0.17$ mas, or $4.4$ pc along the jet. 3. 3.

Brightness temperature in the most compact and prominent components at 7 mm, the core and ST1, rises up to $>10^{13}$ K in the source frame. This agrees with previous discovery of the extremely high brightness temperatures in 3C 273 with RadioAstron (Kovalev et al., 2016; Johnson et al., 2016). Given the moderate measured speed of the model components, this favors the non-equipartition scenario, with the particle energy density dominance. 4. 4.

Magnetic field strength in the 7 mm core is estimated to be $0.2$ G from the core-shift analysis for the epoch 2009 Aug 28. Equipartition between the magnetic field and particles energy density was assumed for this epoch based on the brightness temperature measurements. This magnetic field value is comparable to estimates of Savolainen et al. (2008) made on the basis of synchrotron self-absorbed spectrum turn-over measurements. Other three epochs (2009 Oct 25, 2009 Dec 05 and 2010 Jan 26) were observed by us during a 7 mm flare for which we believe the magnetic field in the core to drop significantly. 5. 5.

We estimated the jet bulk speed on the basis of 7 mm light curves of three adjacent stationary jet features. Estimated speed during the flare $B_{\mathrm{7}}$ was $12c$. This is $1.5$ times higher than the speed of the moving component c4, which was associated with the flare, as measured downstream the jet. 6. 6.

We detected both temporal variations of spectral index in the core and spatial variations along the jet ridge line. The spectrum of the apparent core at any of the observed frequencies is flat or inverted and becomes harder at high frequencies while a radio flare develops reaching values of $\alpha_{\mathrm{43-24}}=1.5$ between 43 GHz and 24 GHz. The spectral index steepens along the jet ridge line by $0.036$ $\mathrm{mas}^{-1}$, or $10^{-3}$ per deprojected parsec.

Study of the transition zone from the opaque core to the optically-thin jet at high frequencies shows that this transition extends by several parsecs and could be resolved with our observations at both 24 GHz and 8 GHz.

Acknowledgments

We thank the anonymous referee and Alan Roy for thorough reading of the manuscript and constructive suggestions which helped to significantly improve the paper. Authors thank Sebastian Kiehlmann for participation in discussion of light curve cross-correlation, and Andrei Lobanov for review and useful discussion. Core-shift analysis is supported by Russian Science Foundation grant 16-12-10481. Light curve and radio-to-$\gamma$-ray relation analysis is supported by the Russian Foundation for Basic Research (RFBR) grant 13-02-12103. TH was supported by the Academy of Finland project number 267324. TS was funded by the Academy of Finland projects 274477 and 284495. This research has made use of data from the MOJAVE database that is maintained by the MOJAVE team (Lister et al., 2009). This study makes use of 43 GHz VLBA data from the VLBA-BU Blazar Monitoring Program (VLBA-BU-BLAZAR; http://www.bu.edu/blazars/VLBAproject.html), funded by NASA through the Fermi Guest Investigator Program. The VLBA is an instrument of the National Radio Astronomy Observatory. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated by Associated Universities, Inc. This research has made use of data from the University of Michigan Radio Astronomy Observatory which has been supported by the University of Michigan and by a series of grants from the National Science Foundation, most recently AST-0607523.