Localizing the γ-ray emitting region in the blazar TXS 2013+370

Aims. The γ-ray production mechanism and its localization in blazars are still a matter of debate. The main goal of this paper is to constrain the location of the high-energy emission in the blazar TXS 2013+370 and to study the physical and geometrical properties of the inner jet region on sub-pc scales. Methods. TXS 2013+370 was monitored during 2002–2013 with VLBI at 15, 22, 43, and 86 GHz, which allowed us to image the jet base with an angular resolution of ≥0.4 pc. By employing CLEAN imaging and Gaussian model-fitting, we performed a thorough kinematic analysis at multiple frequencies, which provided estimates of the jet speed, orientation, and component ejection times. Additionally, we studied the jet expansion profile and used the information on the jet geometry to estimate the location of the jet apex. VLBI data were combined with single-dish measurements to search for correlated activity between the radio, mm, and γ-ray emission. For this purpose, we employed a cross-correlation analysis, supported by several significance tests. Results. The high-resolution VLBI imaging revealed the existence of a spatially bent jet, described by co-existing moving emission features and stationary features. New jet features, labeled as A1, N, and N1, are observed to emerge from the core, accompanied by flaring activity in radio/mmbands and γ-rays. The analysis of the transverse jet width profile constrains the location of the mm core to lie ≤2 pc downstream of the jet apex, and also reveals the existence of a transition from parabolic to conical jet expansion at a distance of ∼54 pc from the core, corresponding to ∼1.5×106 Schwarzschild radii. The cross-correlation analysis of the broad-band variability reveals a strong correlation between the radio-mm and γ-ray data, with the 1 mm emission lagging ∼49 days behind the γ-rays. Based on this, we infer that the high energy emission is produced at a distance of the order of ∼1 pc from the jet apex, suggesting that the seed photon fields for the external Compton mechanism originate either in the dusty torus or in the broad-line region.


Introduction
Matter accretion onto a supermassive black hole (SMBH) -one of the most efficient energy production mechanisms in the Universe (Davis & Laor 2011) -is thought to power Active Galactic Nuclei (AGN). This physical process, combined with the presence of strong magnetic fields accumulated in the accretion disk, is linked to the formation of collimated plasma jets emanating Reduced images of Figs. 1-3 are only available at the CDS via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5)  from the nuclear regions and propagating up to large distances (Blandford & Znajek 1977;Blandford & Payne 1982). Powerful AGN jets, especially those oriented close to our line of sight (i.e., blazars), are observed to radiate over the entire electromagnetic spectrum and to display extreme variability (e.g., Aharonian et al. 2007). A persistent and intriguing question, specifically relevant to the high-energy part of the emission, concerns the physical processes driving the jet γ-ray emission and the location of its production site. Observational findings and theoretical models suggest that the high energy emission can be produced in regions close to the central engine as well as further downstream in the jet (e.g., Jorstad et al. 2001;Rani et al. 2014;Madejski & Sikora 2016, and references therein).
In the framework of leptonic models, the broadband AGN emission is thought to be produced by leptons (e − and e + ) through synchrotron and inverse Compton (IC) processes. The photons that are scattered up to γ-ray energies can either be the same synchrotron photons radiated by the jet (Synchrotron-Self-Compton, Maraschi et al. 1992) or can originate in the jet surroundings (External Compton) (e.g., Ghisellini & Tavecchio 2008;Dermer et al. 2009). With increasing distance from the black hole, possible reservoirs of these seed photons are the accretion disk (Dermer et al. 1992), the broad-line region (BLR, Dermer & Schlickeiser 1993;Sikora et al. 1994;Poutanen & Stern 2010;Dotson et al. 2012), or the dusty torus (Błażejowski et al. 2000;Kataoka et al. 1999). At larger distances, Cosmic Microwave Background (CMB) photons may play a role (Celotti & Ghisellini 2008). If ultra-relativistic protons, as well as leptons, compose the jet, γ-rays could also be produced through proton-synchrotron or photo-pion production (Mannheim & Biermann 1992;Aharonian 2000).
Very-long-baseline interferometry (VLBI) observations at millimeter wavelengths are ideally suited to investigate the origin of γ-ray emission in blazars as they provide a very sharp view of the innermost jet regions, generally affected by synchrotron opacity effects at longer wavelengths. When used in combination with monitoring of the broadband emission variability, mm-VLBI can enable us to effectively pinpoint the location of the γ-ray production site (see Boccardi et al. 2017, and references therein).
In this study, we focus on the compact radio source TXS 2013+370, which is associated with a γ-ray loud object (Mukherjee et al. 2000;Halpern et al. 2001;Kara et al. 2012;Lico et al. 2016) of an uncertain type (Massaro et al. 2015) at redshift z = 0.859 (Shaw et al. 2013) and hosting a SMBH with a mass of 4 × 10 8 M (Ghisellini & Tavecchio 2015). Based on ultra-high-resolution VLBI and single-dish observations performed over 10 years, we obtained a detailed description of the morphological evolution and the variability properties of the radio jet. By combining a kinematic and geometrical analysis of the jet base with the investigation of correlated flux density variability in the cm-and mm-radio and in the γ-ray bands, we were able to constrain the location of the γ-ray production region.
The article is organized as follows. In Sect. 2, we present the multi-frequency data set and the data reduction techniques; in Sect. 3, we report on the results from the VLBI study and the multi-band variability analysis; in Sect. 4, we discuss their implications for the γ-ray production; in Sect. 5, we summarize our conclusions. For our calculations, we adopt the following cosmological parameters: Ω M = 0.27, Ω Λ = 0.73, H 0 = 71 km s −1 Mpc −1 (similar to those used by Lister et al. 2016, and references therein), which result in a luminosity distance D L = 5.489 Gpc and a linear-to-angular size conversion of 7.7 pc mas −1 for the redshift of z = 0.859.  (Kardashev et al. 2013), in combination with the VLBI ground array performed simultaneous observations at two frequencies, 22 GHz and 5 GHz, to facilitate fringe search at the space-ground baselines. In this article we consider the 22 GHz data.

Observations, data calibration and imaging
The space-VLBI data were correlated with a special RadioAstron-enabled version  of DiFX software correlator running on a desktop computer. The fringe search performed with PIMA (Petrov et al. 2011) resulted in space-ground fringe detection at baselines up to 1.7 Earth diameters. The full description of the space-VLBI experiment will be given by Sokolovsky et al. (in prep.). At 22 GHz, we re-imaged the source based on the data originally presented by Sokolovsky (2014) and Kardashev et al. (2015), who also show a 5 GHz image.
At 15 GHz we re-analyzed fifteen epochs of Very Long Baseline Array (VLBA, Napier 1994) observations that cover a period from 2002 to 2012 and are publicly available at the MOJAVE data archive 2 (see Lister et al. 2009aLister et al. , 2011Lister et al. , 2018, and references therein).

VLBI data calibration
The data reduction was performed using the National Radio Astronomy Observatory's (NRAO) Astronomical Image Processing System (AIPS, Greisen 1990). The calibration of the GMVA data at 86 and 43 GHz was performed in the standard manner (e.g., Nair et al. 2019): after an initial parallactic angle correction of the phases, we determined the inter-band phase and delay offsets between the intermediate frequencies using some high signal-to-noise ratio scans (manual phase calibration, Martí-Vidal et al. 2012). After the phase alignment across the observing band, the global fringe fitting was performed (Schwab & Cotton 1983), correcting for the residual delays and phases with respect to a chosen reference antenna. Finally, the visibility amplitudes were calibrated, taking into account corrections for atmospheric opacity computed using the measured system temperatures and gain-elevation curves of each telescope. The RadioAstron data were analyzed in a similar fashion as described in Gómez et al. (2016) and Bruni et al. (2017). The calibration of the 15 GHz data was carried out by the MOJAVE team, following the procedure described in Lister et al. (2009b).

VLBI imaging and model-fitting
Frequency and time-averaged data were imported to DIFMAP (Shepherd et al. 1994) for the imaging and Gaussian modelfitting. Before imaging, the visibilities were carefully inspected and spurious data points were flagged. Using the CLEAN algorithm (Högbom 1974) and SELFCAL procedures, which are implemented in DIFMAP, we created the radio images of the source. While imaging the ground-based data at 15, 43, and 86 GHz was straightforward thanks to a large number of stations, the imaging of the RadioAstron 22 GHz data was more challenging due to the limited uv-coverage. Five stations (including the space telescope) provided useful data (Table B.1). Robledo 70 m recorded only right-hand circular polarization while the space telescope recorded only left-hand circular polarization at 22 GHz (other telescopes recorded both polarizations), which resulted in no space-ground fringes to Robledo 70 m. The amplitude calibration of the Jodrell Bank Mark II telescope was corrupted for an unknown technical reason. A Gaussian source model derived from the near-in-time 43 GHz data was used to correct the dimensional circular Gaussian components model the flux density distribution along the jet. The data were imaged under a uniform weighting scheme and without uv-tapering (all the visibilities weighed equally, independent of their uv-distance). The contour levels are set to 0. 25,0.5, 1, 2, 4, 8, 16, 32, and 64% Jy beam −1 of each image peak flux density (see Table B.1). All the images are convolved with a common beam of 0.8 × 0.5 mas at PA −16 • . The time corresponding to each image is indicated in the x-axis.
amplitude of Mark II telescope data before producing the image with several iterations of CLEAN and SELFCAL.
For all data sets, the jet brightness distribution was then parameterized by fitting two-dimensional Gaussian components to the fully calibrated visibility data by using the MODELFIT algorithm, which is implemented in DIFMAP. The uncertainty in the component positions was set to one-fifth of the beam size if the component's size was smaller than the equivalent circular beam b = (b max b min ) 1/2 , otherwise, we assumed one-fifth of the component full width at half maximum (FWHM) as an estimate of the uncertainty. The uncertainty in the position angle, PA, was calculated based on the positional uncertainty ∆X using the trigonometric formula ∆PA = arctan(∆X/r), with r the radial separation in mas. For the component flux density and FWHM, we adopt an uncertainty of 10%, following Lister et al. (2009bLister et al. ( , 2013 and Karamanavis et al. (2016) respectively. The resulting images are shown in Figs. 1-3, while the parameters of the Gaussian components are reported in Tables B.2-B.5.
Since the source is located very close to the Galactic plane, where the column density of the interstellar medium is relatively high, there is the possibility for the images to be affected by interstellar scattering. In order to verify the impact of scattering in TXS 2013+370, we investigated two possible observational effects that could be introduced by it: angular broadening and fast flux density variations. By analyzing the dependence of the angular size with frequency and based on the variability properties inferred from an Effelsberg monitoring at 5 GHz, we conclude that none of these effects play a dominant role at the frequencies considered in our study. The detailed analysis is presented in Appendix A.

Fermi-LAT data analysis
The Large Area Telescope (LAT) is a pair-conversion detector sensitive to γ-rays from below 20 MeV to more than 300 GeV. It was launched on June 11, 2008 as the main scientific instrument on board the Fermi Gamma-ray Space Telescope (Atwood et al. 2009). We used the Python package Fermipy (Wood et al. 2017) throughout the analysis and assumed as a starting sky model the Fermi-LAT third source catalog (3FGL) (Acero et al. 2015). We consider a Region of Interest (ROI) of 10 • around the target position and include in the model all point sources from the 3FGL within 15 • of the ROI center, together with the corresponding model for the Galactic and isotropic diffuse emission (gll_iem_v06.fits and iso_P8R2_SOURCE_V6_v06.txt, respectively). We performed a binned analysis with 10 bins per decade in energy and 0.1 • binning in space, in the energy range 0.1-300 GeV. We first performed a likelihood analysis over the full-time range considered here, i.e., 2008.58 15:43:36.000 UTC to 2017.17 15:43:36 UTC. We fit the ROI with the initial 3FGL model, freeing all the parameters of the target source and the normalization of all sources within 5 • of the ROI center. Since our data set more than doubles the integration time with respect to the 3FGL catalog, we look for new sources with an iterative procedure. We produce a map of test statistic (TS). The TS is defined as  1.2, 2.4, 4.8, 9.6, 19.2, 38.4, and 76.8% of the peak flux density (see Table B.1); at 86 GHz, images are restored with a beam of 0.16 × 0.08 mas, 0 • , and the contour levels are set to 0.5, 1, 2, 4, 8, 16, 32, and 64% of the peak flux density (see Table B.1). The two beams are displayed on the left and right corners respectively. The time stamp for each image is indicated on the x-axis. source. A value of TS = 25 corresponds to a significance level of 4.2σ (Mattox et al. 1996). A TS map is produced by inserting a test source at each map pixel and evaluating its significance over the current model. We look for TS > 25 peaks in the TS map, with a minimum separation of 0.3 • and add a new point source to the model for each peak, assuming a power-law spectrum. We then fit the ROI again and produce a new TS map. This process is iterated until all significant excesses are modeled out 3 . We also perform a localization analysis on the target source and all new sources with TS > 25 found in the ROI.
For the variability analysis, we perform a likelihood fit in each time bin, using the average model as a starting point. We first attempt a fit leaving the full spectrum of the target source, which is described by a LogParabola, free to vary. If the statistics do not allow the fit to converge or results in a non-detection (TS < 25), we fix all parameters except the target source's normalization. We consider the target source to be detected if TS > 10 in the corresponding bin and the signal-to-noise ratio (i.e., flux over its error) in that bin is larger than two. If this is not the case, we report a 95% confidence upper limit.

Single-dish radio light curves
Single-dish radio data, contemporaneous to the VLBI observations, were provided by the 40-m telescope of the Owens Valley Radio Observatory (OVRO) at 15 GHz (Richards et al. 2011). The data at 235 GHz were obtained with the 8-element Submillimeter Array (SMA) (Gurwell et al. 2007) for the period 2008-2017.
As shown in Fig. 4, where the OVRO, SMA, and Fermi light curves are presented, the source was quite active during our monitoring period, and several flaring episodes occurred both at low and high energies. The investigation of possible correlated variability between the considered energy bands is presented in Sect. 3.3. 3 In our case, the iterative source finding procedure resulted in the addition of 19 new point sources in total. This relatively high number is justified by the increased integration time of our analysis with respect to the 3FGL, and by the fact that the source lies close to the Galactic plane. The latter implies that some of the point sources added might correspond to un-modeled diffuse emission. However, none of these sources has a flux comparable to the target source.

Source structure and jet kinematics
Our imaging at 15, 22, 43, and 86 GHz shows that TXS 2013+ 370 features a bright core and a bent jet extending to 5 mas at the lowest frequency. We determined the kinematics of the individual jet features using the parameters of the Gaussian components derived from the model fits of each data set. The component cross-identification between the observing epochs (and frequencies) was done by comparing their positions, flux densities, and sizes. The angular proper motion µ of the components was computed through linear fits of their radial core separation as a function of time. The apparent speed β app is related to µ and the intrinsic speed β as (Urry & Padovani 1995;Rees 1966): where θ is the viewing angle of the jet, µ is the proper motion in rad s −1 , D L is the luminosity distance in m, z the source redshift and c is the speed of light in m s −1 (e.g., Hogg 1999).
15 GHz. Figure 1 presents the 15 GHz contour images with super-imposed circular Gaussian model-fit components, while the parameters of each component are listed in Table B.2.
At this frequency, the source is well modeled by four circular Gaussian components, a core component (labeled as Core) and three jet features (C3, C2, C1), numbered in order of decreasing distance from the core. During the observing interval 2002-2012, component C3 appears quasi-stationary (Lind & Blandford 1985), oscillating around an average distance of r ∼ 0.2 mas from the center (Fig. 5, top-left panel). Numerous studies have shown that when a moving shock passes through a stationary one, the latter could be displaced in position for a short time and then return to its initial position (see Rani et al. 2015, and references therein). Component C3 seems to exhibit such behavior. The features C2 and C1 are, instead, moving components, separating from the core with apparent superluminal speed (see Table 1). From our kinematic analysis, we infer that C2 is the fastest component of the jet, moving with an apparent speed of β app = 13.8 ± 0.9. We note that a quadratic fit was performed by Lister et al. (2016) Fig. 2) respectively. The vertical lines designate the estimated ejection time, whereas the width of the shadowed areas indicates the uncertainty of this estimation, based on the uncertainty of component A1. We adopt for the N as ejection time the peak of the 43 GHz core flare, whereas for N1, we take the appearance time from the 22 GHz VLBI image. We set for N and N1 the same ejection time uncertainty as the most well-defined feature C2, owing to lack of data. measurement uncertainty. Component C1 is the second fastest feature of the jet, moving at a speed of β app = 7.0 ± 0.8. We also note that a new component, labeled A1, becomes visible at 15 GHz in the last two VLBI epochs (after 2011.53), downstream of C3. See the last two images in Fig. 1, and the component separation as a function of time in Fig. 5 (top left panel). The observation of A1 in only two epochs limits the accuracy of its speed determination. Based on the 15 GHz data points, we obtain an apparent speed of β app = 4.2 ± 11.7, which is in agreement with the value derived for the same component at 43 and 86 GHz (see below).
22 GHz. Figure 2 shows the 22 GHz image of the inner jet region (pc-scale) of TXS 2013+370 as observed by space-VLBI with RadioAstron (details in Sokolovsky 2014). The image shows a core, which is elongated in the East-West direction and along the major axis of the elliptical observing beam. The jet appears to be propagating towards a South-West direction, with a slight curvature towards the southern direction beyond r = 0.4 mas. The overall morphology of the VLBI structure is in good agreement with the 43 and 86 GHz images (Fig. 3). Two components can model the central region, the core and a new feature, labeled as N1 (the second grey area in Fig. 2). The appearance and the propagation of N1 probably is connected with the close-in-time increased activity seen in the radio and γ-ray bands (see Fig. 4). Further downstream, three additional Gaussian components represent the jet. Based on their core separation and flux density, we identify these features with the components C3 and A1 at 15 GHz, and with component N at 86 GHz (see below).
43 and 86 GHz. Figure 3 shows the model-fit images of the source at high radio frequencies. Only the most compact regions in the inner jet are visible. The source structure can be modeled by a compact core and three jet components, the innermost of which, A2, is not visible at lower frequencies. The core separation of each component as a function of time is shown in Core Separation (mas)  by Fermi-LAT (see Fig. 4). Further downstream, component A1, which is resolved at 15 GHz only after 2011, moves with an apparent speed β app = 4.0 ± 0.7. The motion of the outermost component, C3, is not well constrained at 86 GHz, since in the region occupied by C3 the jet becomes faint and partially resolved, and after 2009 there is a possible blending between C3 and A1. However, the 43 GHz data show a well defined C3 component in the first three epochs, from which we infer a proper motion µ = (0.07 ± 0.02) mas yr −1 , which corresponds to an apparent speed of β app = 3.3 ± 1.1. This result is comparable to the speed of the nearest upstream component A1. As discussed above, C3 appears quasi-stationary in the long 15 GHz monitoring, therefore the small displacement observed at 43 GHz in the short time is likely associated with flowing plasma crossing C3 at that time.
In summary, the source TXS 2013+370 shows a bent jet, curving from an east-western to a southern orientation. The VLBI core is much more variable in flux density than the jet components. These move at apparent superluminal speeds of ∼(3-14) c, indicating highly relativistic motion of the jet plasma. Near the core, stationary components, as well as newly ejected features, are observed. The ejection times for several jet components could be estimated and are reported in Table 1. A discussion of possible reasons for the wide range of apparent speeds and on the implications for the jet intrinsic parameters (e.g., Lorentz factor and viewing angle) is reported in Sect. 4.1.

Location of the jet apex
In this section, we aim to constrain the location of the jet apex with respect to the VLBI core. Such an estimate can be obtained through an analysis of the frequency-dependent shift of the VLBI core position, caused, in the simplest scenario, by synchrotron opacity and self-absorption. In the case of TXS 2013+370, however, this approach cannot be used directly, since a formal 2D cross-correlation using close-in-time pairs of images at 15 GHz, 43 GHz, and 86 GHz, yielded no measurable shifts. As we will discuss in Sect. 4 based on the variability time lags, the core shift is indeed much smaller than our resolution limits in the considered frequency regime. In order to obtain an estimate of the jet apex location with respect to the VLBI core, we then follow a geometrical approach based on the investigation of the transverse jet expansion profile. At each frequency, we convolved all images with the average equivalent circular beam (0.08 mas at 86 GHz, 0.16 mas at 43 GHz, and 0.75 mas at 15 GHz) and created stacked images. We have not considered the 22 GHz RadioAstron data since there is only one epoch available. We then measured the jet width as a function of core separation in stacked images. This was A112, page 6 of 16 done by slicing the jet pixel-by-pixel in the direction perpendicular to the jet axis, and by fitting single Gaussian profiles to the transverse intensity distribution to infer the jet width. For the error of the jet width, we used one-tenth of the convolved FWHM.
The expansion profile is presented in Fig. 6, top panel. Distances along the x-axis are relative to the position of the VLBI core, which we assume as being fixed, due to the negligible coreshift. The de-convolved FWHM values inferred for the inner jet region based on the 86 and 43 GHz data smoothly connect to those inferred at 15 GHz for the outer regions. By fitting a power-law of the form d = ar b to the 15 GHz data at distances between ∼0.4 and ∼3 mas from the core, we infer that the jet has a conical shape (d ∝ r (1.02±0.01) ). However, this power law does not describe well the higher frequency data, as in the inner jet we observe a flattening of the expansion profile. By fitting a power law of the same form to the 43 GHz and 86 GHz data, we obtain d ∝ r (0.49±0.04) , i.e., the jet has a parabolic shape in the proximity of the black hole. This is expected based on theoretical models for jet formation (e.g., Meier et al. 2001, and references therein), which predict the jet to be actively collimating and accelerating at its onset. While we have not considered in the fit the 15 GHz data relative to the inner jet, the lower resolution data points also lie on the same profile. The fact that the jet is collimating on the scales probed by mm-VLBI is also evident by examining the evolution of the apparent opening angle with distance (Fig. 6, bottom panel). The angle decreases in the inner ∼0.5 mas, and then reaches a roughly constant value of ∼23 • on the scales probed at 15 GHz, until the recollimation region at a distance of ∼4 mas. As observed in other more nearby jets like M 87 (Asada & Nakamura 2012) or Cygnus A (Boccardi et al. 2016), for which higher spatial resolution is achieved, a single parabolic profile describes the expansion of the jet from its onset up to the parsec scales. Therefore, we take the inferred expansion law to back-extrapolate the location of the jet apex, which should have an approximately zero width, with respect to the 86 GHz core. Such a method has been applied in literature before (Agudo et al. 2011;Karamanavis et al. 2016), considering the core size as a reference and then assuming a conical expansion between the black hole and the core. Similarly, we can compute the distance R it takes for the jet to reach the width of the core, assuming the inferred parabolic expansion rate. Based on the analysis of the stacked image at 86 GHz, the jet width at the location of the emission peak is d c ≤ (0.04±0.01) mas, from which we infer that the jet apex is located at a distance of R ≤ ((0.019 ± 0.009)/ sin θ) mas, or ((0.146 ± 0.07)/ sin θ) pc upstream, where θ is the viewing angle of the jet.

Multi-band variability
In order to investigate possible correlated variability across the observing bands, we applied a discrete cross-correlation function (DCCF; Edelson & Krolik 1988). To test the significance of the DCF peak, we simulated the light curves using the Emmanoulopoulos et al. (2013) method implemented by Connolly (2015). We first determined the power spectral density (PSD) slope of the observed γ-ray light curve. Since we ignored the upper limits in the observed γ-ray light curve, the data are not evenly sampled. We therefore interpolated them using the cubic spline method in R 4 . As the data are manipulated  via interpolation, we tested our PSD results using simulations. Specifically, we first calculated the slope and normalization of the original PSD to determine its slope and normalization. The PSD slope and normalization were then used as variables in the simulation algorithm. The next step sampled the simulated light curves at the same bin width as the observed ones. Finally, for each pair of slope and normalization, we simulated a hundred light curves. The PSD slope of each simulated light curve was then averaged to determine its mean and standard deviation. The detailed step-by-step procedure is described by Chidiac et al. (2016). Our analysis suggests that the observed source variability can be well described by a PSD slope of −(0.9 ± 0.2). The method also takes into account the underlying probability density function (PDF) of the given light curve. We simulated a total of 5000 light curves using the best-fit PDF and PSD parameters. The simulated data have been correlated with the observed radio data at 235 GHz, and 15 GHz bands to calculate the 99% confidence levels. A correlation is found between the γ-ray and the 235 GHz light curves, with the radio following the high-energy activity by (45 ± 30) days, whereas between 235 and 15 GHz the time lag is (15 ± 15) days (see Fig. 7). The confidence level of this correlation exceeds 99%. The correlation between 15 GHz and γ-rays exceeds the 95% significance level with a delay of (45 ± 30) days, which is above the 2σ statistical threshold. In order to cross-check the degree of correlation, we applied a Spearman's correlation test (Spearman 1904). The Spearman's rho correlation coefficient (ρ) is a statistical A112, page 7 of 16 A&A 634, A112 (2020) Fig. 7. DCCF results between the γ-ray and 15 GHz light curves (left), 15-235 GHz (middle) and γ-rays-235 GHz (right). Positive time lags indicate that γ-ray activity leads the activity in radio. The significance of the correlations is displayed by a dotted line for the 2σ and by the dashed line for the 3σ level. measure of the strength of a monotonic relationship between paired data, which in this case corresponds to the differences between the fluxes of the light curves. Specifically, we used two light curves at a time, assuming as a reference the delayed one, as established by the DCCF analysis. Thus, we shifted the other data set from 0 to 200 days with a step of 1 day and taking into account a coherence time of 2 days. For each iteration, we calculated the p-value to quantify the significance of this statistical hypothesis. We performed this procedure for the pairs γ-ray-235 GHz, γ-ray-15 GHz, and 15 GHz-235 GHz. Errors in the time-lag estimation have been set equal to the data sampling interval. The time lag between γ-rays-15 GHz was calculated based on the Spearman's Rho results of the other two data sets. The obtained time lag between γ-ray-15 GHz is (56 ± 30) days with p-value = 1.5 × 10 −3 , for γ-ray-235 GHz is (52 ± 30) days with p-value = 2.1 × 10 −4 , whereas the time lag between 235 and 15 GHz was found to be (11 ± 6) days with p-value = 9.5 × 10 −9 . As a sanity check, we searched for correlated activity between the 15 GHz data set and itself. The result was, as expected, ρ = 0.997 with p-value = 0 at 0 days shift. Both methods, the DCCF and the Spearman's rho test, showed that the high-energy emission leads the activity in the radio band (see a summary of the results in Table 2). This indicates that the high-energy event may have occurred in a region that is opaque to radio waves (Fuhrmann et al. 2014;Pushkarev et al. 2010) or, alternatively, that a denser source of seed photons for inverse Compton scattering is present upstream of the main millimeter-wave emission region.
In summary, we conclude that the γ-rays lead the variability. In the following analysis, we will assume the mean time lag for each pair of light curves. Hence, we consider that the γ-ray variability is ahead of 235 GHz by (49 ± 30) days. The 15 GHz data lag behind by another (13 ± 11) days (relative to the 235 GHz data), while the lag between 15 GHz and γ-rays is (51 ± 30) days (see Col. 4 in Table 2).

Intrinsic jet parameters
The diversity of apparent speeds found in TXS 2013+370, with values ranging from moderate (β app ∼ 3) to high (β app ∼ 14, see Table 1), is quite common among extragalactic jets (e.g., Fromm et al. 2013;Rani et al. 2015;Karamanavis et al. 2016;Jorstad et al. 2017) and could result from several effects. The VLBI images show a jet bending from west to south at a distance of 0.1-0.2 mas from the core, and then again towards west on the scales probed at 15 GHz. Therefore a first hypothesis is that the apparent speed variations are geometry-dependent, with the apparent speed increasing when the jet points closer to the line of sight. By considering the position angles of the moving components in Tables B.2, B.4, and B.5, there is indeed an indication that the slowest features are those moving towards south (A1, C3), while the highest speed is observed for C2, which follows a trajectory towards west and south-west. Another possibility is that the Lorentz factor is not constant along the jet, but increases as a function of distance. This hypothesis is supported by the observation of active jet collimation in the inner 0.5 mas of the jet (see Sect. 3.2), indicating that the terminal Lorentz factor is not yet reached on the scales probed by mm-VLBI, where the slowest apparent speeds are measured. Ultimately, it is likely that both geometrical effects and intrinsic variations of the bulk speed cause the observation of this wide range of apparent speeds in TXS 2013+370.
Based on this premise and the results of the kinematic analysis, in the following we estimate some of the intrinsic jet parameters. The observation of a maximum speed of ∼13.9 c for C2 implies that at a distance of ∼1-2 mas from the core the flow has a minimum Lorentz factor Γ min = 14.0 ± 0.8, since Γ min is expressed as A112, page 8 of 16 The viewing angle that maximizes the apparent speed for a given Lorentz-factor is called the critical viewing angle θ c : (3) For C2 we obtain θ c = (4.1±0.2) • . The critical viewing angle of the fastest moving component is often assumed in the literature to be equal to the characteristic jet viewing angle (e.g., Vermeulen & Cohen 1994;Lister & Marscher 1997;Cohen et al. 2007). This angle could vary along the jet and be larger for those regions of the jet which are pointing towards south, and for which the slower speeds are measured (e.g., θ c 16 • for components A1 and C3). However, for this paper we are mostly interested in the viewing angle of the inner jet and of the core region, where the plasma moves in a direction similar to C2. Therefore we will adopt this viewing angle value in the following for the deprojection.
Concerning the bulk Lorentz factors, the high apparent speed measured for C2 indicates that the plasma is fast and highly boosted in the regions probed at 15 GHz (Γ > 14). On the other hand, the low apparent speeds and the increase of the jet opening angle towards the jet apex suggests that the Lorentz factor is lower in the vicinity of the core than further downstream. For a jet orientation at θ c = (4.1 ± 0.2) • and the observed low apparent speed of A1 (β app ∼ 4.2) we can estimate the jet Lorentz-factor to be of the order of 6 in the regions probed by millimeter VLBI. This also implies that the Doppler factor δ = 1/(γ · (1 − β cos θ)) in the core region is not very large (δ ≤ 10), which is in good agreement with the moderate variability of the source.
Having measured the opening angle and the Lorentz factor at different locations along the jet, we can also test how are these two quantities related to each other. Based on hydrodynamical (Blandford & Königl 1979) and magneto-hydrodynamical models (Komissarov et al. 2007) we expect the Lorentz factor Γ and the intrinsic jet opening angle φ to be inversely proportional, with Γφ < 1 for a causally connected jet. This is consistent with our results. In the outer jet, the constant apparent opening angle of ∼23 • implies an intrinsic full opening angle φ ∼ 1.6 • , thus the product Γφ yields ∼0.4 rad for Γ = 14. In the inner jet, speeds should be lower, as discussed, and the opening angle is larger. Assuming as a reference an apparent opening angle of ∼40 • , measured at distances of 0.1−0.3 mas from the core (Fig. 6, bottom panel), and Γ = 6 as estimated above for the jet base, we obtain Γφ ∼ 0.30 rad. Both products in the two regions are in good agreement with the median value obtained for the MOJAVE sample, Γφ = 0.35 rad (Pushkarev et al. 2017), as well as with the results of a statistical modeling considering the same population (Clausen-Brown et al. 2013) 5 .

Location of the γ-ray emission
In Sect. 3, we investigated the existence of correlations between the variability observed in the radio band (15 GHz and 235 GHz) and the γ-ray band, and we derived time lags indicating that the high-energy activity leads the one in radio (see Table 2). Having identified the most likely intrinsic parameters of the innermost jet regions, we can now translate these time lags into deprojected physical scales. Following León-Tavares et al. (2011), Pushkarev et al. (2010, the distance between the dominant emission regions in two different bands ∆r is related to the time lag 5 Note that Clausen-Brown et al. (2013) consider the half opening angle, obtaining Γθ = 0.2. Since we consider the full opening angle, there is a factor-of-2 difference in the product. ∆t as: In our calculations, we adopt the mean time lags obtained in the DCCF analysis and through the Spearman's rho test analysis, the apparent speed of the innermost moving component, A1, as the best representative of the plasma speed in the region close to the core, and a viewing angle of 4.1 • . The de-projected distances are in parsecs, for the three considered pairs of variability curves: Both the 2D cross-correlation between pairs of VLBI images and the results of the variability analysis indicate that the core shift in the frequency range between 15 GHz and 86 GHz is well below the resolution limit of our observations. In fact, a de-projected shift of 0.34 pc inferred between 235 and 15 GHz would translate into a projected angular separation of only ∼0.003 mas, which will be even smaller between, e.g., 43 GHz and 86 GHz. This result suggests that the observed VLBI cores may not be associated with the unit radio opacity surfaces (τ = 1), but that the core region may be a stationary shock. Negligible core shifts are derived for several blazars in the MOJAVE sample (Pushkarev et al. 2012) and, especially at higher frequencies, the characteristics of the VLBI core often resemble those of a stationary shock (see Jorstad et al. 2017, and references therein).
Based on the analysis of the jet expansion profile presented in Sect. 3.2, this possibly stationary core feature is located at a distance R ≤ (0.019±0.009) mas from the jet apex, which translates to a de-projected separation of ≤(2.05 ± 0.97) pc for a viewing angle of 4.1 • . Following the results of the variability study for the most strongly correlated pair γ-1 mm (∆r = 1.30 ± 0.81 pc), we estimate the γ-ray emission location to be at a distance of ∼(0.75 ± 1.26) pc downstream from the jet apex. By taking into account the estimated error bars, this result tells us that that the high-energy event occurs on scales ranging from sub-parsec to about ∼2 pc distance from the jet apex.
On such scales, the most likely mechanism leading to highenergy production is highly dependent on the source type. In powerful blazars known as Flat-Spectrum Radio Quasars (FRSQs), intense external photon fields originating in the accretion disk, the Broad Line Region (BLR) or the dusty torus are likely to act as seeds for the inverse Compton radiation, while in BL Lac objects these fields are expected to be less prominent or even absent, and the high-energy emission is often well reproduced by synchrotron self-Compton models.
As mentioned in Sect. 1, the classification of TXS 2013+370 as a BL Lac object or an FSRQ is uncertain in blazar catalogs (e.g., Massaro et al. 2015); this is due to the significant Galactic extinction in the source direction, which prevents a solid determination of its optical properties to be obtained. However, the recently determined, relatively high redshift (z = 0.859, Shaw et al. 2013) is uncommon among BL Lacs (see, e.g., their redshift distribution in the latest Fermi-LAT catalog, The Fermi-LAT Collaboration 2019). Moreover, Kara et al. (2012) have examined the spectral energy distribution (SED) of the source, showing it is characterized by a hard X-ray photon index and by the dominance of the inverse Compton component over the synchrotron one in power output, as usually observed in FSRQs. Indeed, reasonable fits of the SED were only obtained when including an external Compton (EC) contribution.
An estimate of the BLR size in this source can be obtained by considering the bolometric luminosity of the accretion disk L d , as determined by Shaw et al. (2013), and then following the approach of Ghisellini & Tavecchio (2015), which assumes a spherical BLR with radius R BLR = 10 17 L 1/2 d cm. Through this method, we obtain a radius of the order of ∼0.07 pc. The analysis presented in this paper points towards distances larger than this for the location of the γ-ray emission. On scales of 1-2 parsecs from the central engine, the dusty torus is the best candidate for providing a rich seed photon field. Indeed, Kara et al. (2012) showed that the best fit in the SED modeling of TXS 2013+370 required equipartition conditions for the dominant emitting region and an external radiation field with a rather low temperature (T ext ∼ 10 2 K), thus possibly originating from cold dust. Temperatures in this range are found in the outer regions of the torus. For instance, previous studies performed in 3C 454.3 andNGC 1068 (Jaffe et al. 2004;Shah et al. 2017) have derived a temperature of the order 300-600 K at a distance of ∼3.5 pc from the SMBH. Therefore, the distance that we obtain from our analysis can also be considered as a lower limit for the outer radius of the dusty torus.
A physical scenario where the EC process is supplied by infrared photons (IR) from the dusty torus is supported by several studies of powerful blazars (Sahayanathan & Godambe 2012;Sikora et al. 2009). Costamante et al. (2018) showed that in the vast majority of the Fermi FSRQs, the γ-ray emission appears to originate outside the BLR. Moreover, recent findings support a torus geometry that deviates significantly from the standard picture of a "donut"-like structure (Carilli et al. 2019;Asmus 2019;Lyu & Rieke 2018;Hönig & Beckert 2007). The torus is likely to be rather clumpy, with polar molecular clouds providing an even richer photon field available for EC scattering.

Transition from parabolic to conical expansion
In summary, we comment on the transition observed in the jet expansion profile, with the jet switching from a parabolic to a conical shape (Sect. 3.2). As the resolution of radio observations increases, this phenomenon is observed in more and more jets and supports the currently most favored physical models for magnetic jet launching. These predict the jet base to be actively collimated and accelerated along an extended region, up to parsec distances from the central engine (e.g., Vlahakis & Konigl 2004;Komissarov et al. 2007).
In TXS 2013+370 the transition is observed at a separation of ∼0.5 mas from the jet apex, corresponding to a de-projected distance of ∼54 parsecs for θ = 4.1 • . For a black hole mass of 4 × 10 8 M (Ghisellini & Tavecchio 2015), we estimate that the jet collimation stops at 1.5 × 10 6 Schwarzschild radii from the black hole, which is of the same order as the transition distance found for M 87 (Asada & Nakamura 2012) and other sources in the MOJAVE sample (Kovalev et al. 2019).
This result indicates that at millimeter wavelengths, we are probing jet regions where the magnetic field is still important and that the γ-ray emission in this source is produced in the magnetically-dominated part of the jet base.

Conclusions
In this study, we gave a complete picture of the radio morphology and jet evolution of the blazar TXS 2013+370 during the last ten years, and we constrained the location of the γ-ray emission site. For this purpose, we employed state-of-art VLBI observations from 2 cm down to 3 mm, as well as space-VLBI data. This unique data set allowed us to investigate the kinematic properties of the source and the jet expansion profile, from which we could derive the linear separation between the VLBI core and the jet apex. In addition, high cadence single-dish observations allowed us to monitor the flux density variability of the source, and to search for correlated activity between the radio and the γ-ray bands. The results can be summarized as follows: Ultimately, space-VLBI RadioAstron data in 2012 showed the emergence of another knot, labeled as N1. The appearance of knot N1 seems to precede the enhanced emission in radio and γ-ray bands. 2. The study of the jet transverse expansion profile allowed us to quantify the distance between the mm VLBI core and jet apex to be R ≤ (2.05 ± 0.97) pc. This analysis also revealed the existence of a geometrical transition from a parabolic to a conical jet shape, taking place at a projected distance of ∼0.5 mas. This corresponds to a de-projected distance of ∼54 parsecs, or of ∼1.5 × 10 6 Schwarzschild radii from the jet base. 3. The correlated activity between the γ-ray and radio bands enabled us to translate the observed time lags to linear distances. The strongest correlation was found between γ-rays and 1 mm activity, with the γ-rays leading by (49 ± 30) days. The estimated delay, taking into account the plasma speed on scales close to the jet apex and the viewing angle of the jet, corresponds to a de-projected distance of ∆r γ−235 = (1.30 ± 0.81) pc. By combining the knowledge of the distance between the jet apex and the mm VLBI core and of the γ-ray production region with respect to the mm VLBI core, we constrained the location of the high-energy production site at (0.75 ± 1.26) pc downstream of the jet apex. For the determined range of distances, the external seed photon field for inverse Compton emission is most likely to be produced by the dusty torus, although an origin in the BLR is also possible.
TXS 2013+370 is located close to the Galactic plane (b ∼ 2 • ) and in the Cygnus super-bubble region, raising the possibility that the source is affected by interstellar scattering. Since interstellar scattering produces a rich range of observational effects that may influence the results of the current study, we discuss it briefly. One manifestation induced by scattering is angular broadening, where the apparent angular size scales approximately as ν −2 (see for example Lazio et al. 2008, and reference therein). In order to test this size-frequency relation, we combined our VLBI measurements with size measurements from Pushkarev & Kovalev (2015). For the error of the angular sizes, we adopt a conservative 10% of the FWHM. Following Lazio et al. (2008), we fit the measured angular sizes to the functional form where θ s and θ i are the scattering and intrinsic source sizes, respectively. We found the best-fitting values for θ s and θ i using a least-squares fit approach, minimizing χ 2 . For the power-law index of the intrinsic size, we considered both k = 0 (i.e., a frequency-independent intrinsic size, for a flat-spectrum radiojet Blandford & Königl 1979) and k = −1 (i.e., a frequency scaling typical for a single inhomogenous synchrotron source Kellermann & Pauliny-Toth 1981), and selected the value of k that produced the lowest χ 2 . The inferred scattering and intrinsic sizes from the fitting are summarized in Table A.1. In Fig. A.1, we plot the observed size versus frequency. The solid lines indicate the best fit of Eq. (A.1) to the observations, and the dashed lines illustrate the "decomposition" of the apparent source angular size into the two physical effects, scattering and intrinsic synchrotron emission. The data show that the angular broadening is prominent below 10 GHz, which is consistent with earlier findings (Spangler et al. 1986;Fey et al. 1989). Above 10 GHz, however, the intrinsic source size begins to dominate, where the size-frequency relation shows the expected slope for a synchrotron self-absorbed (SSA) jet. Another phenomenon that might be induced by interstellar scattering is the so-called "intra-day" variability (IDV, Heeschen 1984;Witzel et al. 1986), which is present in 30−50 % of all flat-spectrum quasars and BL Lac objects in cmwavelengths (Quirrenbach et al. 1992;Lovell 2008). We note that TXS 2013+370 was observed in several IDV monitoring campaigns with the 100 m telescope at Effelsberg at 5 GHz, as part of the coordinated ground support for the RadioAstron space-VLBI experiments (Liu et al. 2018). According to χ 2 -tests, the source did not show significant fast variability over ∼3 days. This can be explained by quenched refractive interstellar scintillation (Narayan 1992), with a source size larger than  k θ s (mas) θ i (mas) χ 2 0 48.2 ± 10.3 (3.6 ± 1.0) ×10 −2 56.9 −1 37.1 ± 3.0 2.4 ± 0.1 3.5 Notes. Columns from left to right: (1) power-law index, (2) observed angular size, (3) intrinsic angular size, (4) chi-square value.
the scattering size at 5 GHz. The Scattering Measure (SM) in the Cygnus region varies by a factor of 2−5 on angular scales of only a few degrees. Although TXS 2013+370 is scatter broadened, it shows a much lower SM than some other prominent AGN in the same region, like e.g., 2005+403, which is separated by only ∼4 • (cf. Fig. 6 in Fey et al. 1989). We, therefore, conclude that for TXS 2013+370, interstellar scattering is not very strong and is dominant only at the longer cm-wavelengths. Thus it should not affect the mm-flux density and imaging significantly.