© ESO, 2016

1. Introduction

Blazars comprise the radio-loud minority of active galactic nuclei (AGN) which have energetic jets seen at small viewing angles. Observationally, they exhibit strong variability across the whole electromagnetic spectrum. However, its origin and extreme characteristics (broadband flux density outbursts with time scales as short as minutes) are still a matter of intense debate. Moreover, uncertainty still enshrouds the exact location of their high-energy component, observed at MeV–TeV energies.

Numerous mechanisms have been proposed to explain this intense activity. Among these are: shocks generated and propagating along the relativistic jet (e.g., Marscher & Gear 1985; Valtaoja et al. 1992), colliding relativistic plasma shells (Spada et al. 2001), precession of the jet (e.g., Lister et al. 2003), helical motion in the jet (e.g., Camenzind & Krockenberger 1992), and Kelvin-Helmholtz instabilities developing in the outflow (e.g., Lobanov & Zensus 2001). Testing the validity of such models means quantitatively contrasting the observational signatures of the mechanisms they invoke with observations. Multi-frequency flux density monitoring can shed light on blazar physics by accessing time and thus spatial scales which are unreachable even by interferometric imaging. Dense, long-term, multi-frequency monitoring offers a unique opportunity for putting blazar variability models to the test and for comparing with the predictions of theoretical frameworks that attempt to explain the emission and overall observed characteristics of these sources.

Time delays between flare maxima are often seen and usually attributed to opacity effects. Flare onsets and maxima appear first at higher frequencies and successively progress towards the lower end of the observed frequency range (e.g., Valtaoja et al. 1992). This is usually connected to the motion of disturbances whose optical depth decreases, while moving downstream the jet (e.g., Marscher & Gear 1985) and expanding adiabatically (van der Laan 1966). Consequently, lower frequency (hence energy) photons are able to escape the synchrotron-emitting region.

The powerful flat-spectrum radio quasar (FSRQ) PKS 1502+106 (OR 103, S3 1502+10) at a redshift of z = 1.8385, D_L = 14 176.8 Mpc (Adelman-McCarthy et al. 2008) harbours a supermassive black hole of ~10⁹M_⊙ (Abdo et al. 2010, and references therein). The focus of this paper, is the intense 2008/2010 outburst seen in PKS 1502+106 (Ciprini 2008). This pronounced broadband flare, observed from radio up to γ-ray energies, triggered the first Fermi-GST multi-frequency campaign covering the electromagnetic spectrum in its entirety, both with filled-aperture instruments and interferometric arrays (Karamanavis et al. 2016; see also Karamanavis 2015, for first results).

The scope of this paper is to parameterize the observed outburst, and extract its relevant light curve parameters, such as the flare amplitude at each frequency, and the cross-band delays. First, using the observed cross-band and frequency-dependent time lags, we estimate the core shift – that is the frequency-dependent position of the core – with an approach other than through traditional multi-frequency very-long-baseline interferometry (VLBI) measurements. With the time-lag core shift at hand, an opacity profile of the source and, under the assumption of equipartition, the magnetic field strength along the synchrotron-self-absorbed jet are obtained. Furthermore, the distance of each core to the vertex of the jet, in combination with the observed delay between radio and γ rays, allows for decisively constraining the location of the high-energy emission.

The paper is structured as follows. In Sect. 2 we present the data set and the time series analysis techniques used. In Sect. 3 our results are presented. Section 4 discusses and Sect. 5 concludes our findings. In the paper, we adopt S ∝ ν^{+ α}, where S is the radio flux density, ν the frequency at which observations are made, and α the spectral index, along with the following cosmological parameters H₀ = 71 km s^-1 Mpc^-1, Ω_m = 0.27, and Ω_Λ = 0.73.

2. Observations and data analysis

Data used in this work have been collected within the framework of the Fermi-GST AGN Multi-frequency Monitoring Alliance (F-GAMMA¹) program (Fuhrmann et al. 2007; Angelakis et al. 2010), which includes observations with the Effelsberg 100-m telescope (EB) in eight frequency bands from 2.64 to 43.05 GHz, and with the IRAM 30-m telescope (PV) at 86.24 and 142.33 GHz. Monthly observations from the EB and PV telescopes were performed in a quasi-simultaneous manner with a typical separation of days ensuring maximum spectral coherency. A detailed description of the observational setup and data reduction is provided elsewhere (Fuhrmann et al. 2014; Angelakis et al. 2015; Nestoras et al. 2016). We also employ data from the Owens Valley Radio Observatory (OVRO) 40-m blazar monitoring program² at 15 GHz (Richards et al. 2011) and at 226.5 GHz from the Submillimeter Array (SMA) observer center database (Gurwell et al. 2007).

2.1. Gaussian process regression

To extract the light curve parameters, such as flare amplitude and time scale at each frequency, we employ a Gaussian process (GP) regression scheme (Rasmussen & Williams 2006). This is a non-parametric approach to the generalized regression and prediction problem, widely used in machine learning applications and lately exploited in astronomy (e.g., Gibson et al. 2012; Aigrain et al. 2012). While traditional fitting assumes a specific functional form f for the model, a second way to learn f from the data is by assigning a prior probability to all functions which one considers more likely, for example, on the grounds of their smoothness, in other words, their infinite differentiability. A Gaussian process is a probability distribution over functions. It constitutes the generalization of the Gaussian distribution of random variables or vectors, into the space of functions (e.g., Rasmussen & Williams 2006; Ivezić et al. 2014).

In the case of AGN light curve fitting, the problem takes the form of a one-dimensional regression where the observables are flux density levels at certain times. For instance, given a set of three observations, or training data points, , our problem reduces to finding the functions, which are drawn from the prior distribution, that pass through all these training points. This distribution of functions is the posterior distribution and its mean, being the expected value, can be considered as the “best-fit function”. A Gaussian process is defined by its mean, μ, and covariance, k. Without loss of generality, one can always assume that μ = 0, since the data can always be shifted to accommodate this assumption.

In the context of Gaussian processes, the covariance function, covariance kernel, or simply kernel, Cov [ f(x_i),f(x_j) ] = k(x_i,x_j), is used to define the covariance between any two function values at points x_i and x_j; in other words, the similarity between data points. It is chosen on the basis of the prior beliefs for the function to be learned. Essentially, the kernel defines the class of functions likely to appear in the prior distribution, which in turn, determine the kind of structure that the specific GP is able to model correctly. Kernels have their own set of parameters, called hyperparameters, since they define the properties of the prior distribution over functions, instead of the functions themselves. A very useful kernel property is that addition or multiplication between two or more of them also produces valid covariance kernels. As such, there is always the option of constructing a covariance kernel that fits the characteristics of the modeling problem³. Here we employ the squared exponential kernel $\mathit{k}\mathrm{\left(}{\mathit{x}}_{\mathit{i}}\mathit{,}{\mathit{x}}_{\mathit{j}}\mathrm{\right)}\mathrm{=}{\mathit{\sigma }}_{\mathit{f}}^{\mathrm{2}}\exp \left[\frac{\mathrm{-}\mathrm{(}{\mathit{x}}_{\mathit{i}}\mathrm{-}{\mathit{x}}_{\mathit{j}}{\mathrm{)}}^{\mathrm{2}}}{\mathrm{2}{\mathit{l}}^{\mathrm{2}}}\right]\mathit{,}$(1)where the hyperparameters are

• l

  the characteristic length scale, denoting the distance in the x dimension after which the function changes significantly, and

• ${\mathit{\sigma }}_{\mathit{f}}^{\mathrm{2}}$

  the variance, mapping the mean distance of the function from its mean. Here, it serves only as a scaling factor.

This is a widely used kernel, owing to its generality and smoothness. The latter is the only assumption when using it, justified by the fact that most time series arising from physical processes have no reason not to be smooth. In the case of data with uncertainty, one only has to add it to the noiseless covariance kernel. In the general case, that is k_n(x_i,x_j) = k(x_i,x_j) + σ²I, where I is the identity matrix and σ² is the hyperparameter describing the noise variance (e.g., Roberts et al. 2012).

Selecting the best set of hyperparameters using the data themselves is referred to as training the GP or, more generally, Bayesian model selection. In essence, one aims to update the prior knowledge in the light of a training data set. One way of doing so is by maximizing the marginal likelihood (Williams & Rasmussen 1996). Letting θ denote the vector of hyperparameters, in the case of the squared exponential kernel (Eq. (1)) we have ${\theta }\mathrm{=}\mathrm{\left\{}\mathit{l,}{\mathit{\sigma }}^{\mathrm{2}}\mathrm{\right\}}\mathit{.}$(2)Then, the probability (or evidence) of the training data, y, given the hyperparameters vector θ, is $\mathit{p}\mathrm{(}{y}\hspace{0.17em}\mathrm{|}\hspace{0.17em}{x}\mathit{,}{\theta }\mathrm{)}\mathit{,}$(3)and the log marginal likelihood is given by $\mathrm{ℒ}\mathrm{=}\log \mathit{p}\mathrm{(}{y}\hspace{0.17em}\mathrm{|}\hspace{0.17em}{x}\mathit{,}{\theta }\mathrm{)}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{{y}}^{\mathit{T}}{{k}}_{\mathit{n}}^{-1}{y}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\log \mathrm{|}{{k}}_{\mathit{n}}\mathrm{|}\mathrm{-}\frac{\mathit{n}}{\mathrm{2}}\log \mathrm{2}\mathit{\pi }\mathit{.}$(4)In the general case of a hyperparameter vector θ = { θ_j | j = 1,...,n }, the derivatives of the log marginal likelihood with respect to each θ_j are $\frac{\mathit{\partial }\mathrm{ℒ}}{\mathit{\partial }{\mathit{\theta }}_{\mathit{j}}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}{{y}}^{\mathit{T}}\frac{\mathit{\partial }{{k}}_{\mathit{n}}}{\mathit{\partial }{\mathit{\theta }}_{\mathit{j}}}{{k}}_{\mathit{n}}^{-1}{y}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{Trace}\left({{k}}_{\mathit{n}}^{-1}\frac{\mathit{\partial }{{k}}_{\mathit{n}}}{\mathit{\partial }{\mathit{\theta }}_{\mathit{j}}}\right)\mathrm{·}$(5)Equation (5) can be used with any numerical gradient optimization algorithm in order to maximize the log marginal likelihood and return the set of best hyperparameters for the problem at hand. A schematic algorithm is given by Rasmussen & Williams (2006).

For the application of the method to the radio light curves of PKS 1502+106, a variant of the algorithm presented by Pedregosa et al. (2011) has been used, developed specifically for this purpose. The complete suite of machine learning programs can be accessed online⁴. Here, we consider the full length of each available light curve. Since consequent flares start without necessarily reaching a quiescent level prior to the onset of a new one, before proceeding with the Gaussian process regression, we subtract the minimum flux density level observed at each frequency, S₀, from the corresponding light curve. To ensure the best unbiased result we perform the process with 100 random initializations of the length scale parameter, l. Upon completion, the posterior mean is returned with a robust uncertainty estimate in the form of a 95% confidence interval for the flux density, along with the set of hyperparameters maximizing the log marginal likelihood. Consequently, we are able to extract the flare amplitude, S_m, the time of maximum flux density, t_m, and the cross-band delay, τ_GP. Additionally, the flare rising and decay times (t_r and t_d, respectively) are extracted from the times of the two flux-density minima bracketing S_m. In Fig. 1 the results of regression at each frequency are visualized. The values of S_m, t_m, t_r, and t_d characterizing the flare are visible in Table 1.

2.2. The discrete cross-correlation function

Putative correlated variability across observing bands, which are based on sampling-rate-limited time series, can be investigated and quantified using the discrete cross-correlation function (DCCF, Edelson & Krolik 1988).

The cross-correlation function (CCF), as a function of the time lag τ, for two discrete and evenly sampled light curves, x(t_i) and y(t_i), is given by $\mathrm{CCF}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathit{N}}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{N}}\frac{\left[\mathit{x}\mathrm{\left(}{\mathit{t}}_{\mathit{i}}\mathrm{\right)}\mathrm{-}\mathit{x̅}\right]\left[\mathit{y}\mathrm{(}{\mathit{t}}_{\mathit{i}}\mathrm{-}\mathit{\tau }\mathrm{)}\mathrm{-}\mathit{y̅}\right]}{{\mathit{\sigma }}_{\mathit{x}}{\mathit{\sigma }}_{\mathit{y}}}\mathit{,}$(6)where

• 

  is the mean value of the light curve x(t_i);

• σ _x

  its standard deviation;

• 

  is the mean of the light curve y(t_i); and

• σ _y

  is the standard deviation of light curve y(t_i).

Uneven sampling of light curves introduces the need for the discrete cross-correlation function (Edelson & Krolik 1988; see also Larsson 2012; Robertson et al. 2015). Contrary to linear interpolation methods, the DCCF takes only the data points themselves into consideration. For two irregularly sampled light curves, the unbinned CCF (UDCF) is ${\mathrm{UDCF}}_{\mathit{ij}}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\mathrm{=}\frac{\mathrm{(}{\mathit{x}}_{\mathit{i}}\mathrm{-}\mathit{x̅}\mathrm{)}\left({\mathit{y}}_{\mathit{i}}\mathrm{-}\mathit{y̅}\right)}{\sqrt{\mathrm{(}{\mathit{\sigma }}_{\mathit{x}}^{\mathrm{2}}\mathrm{-}\mathit{\epsilon ̅}\begin{array}{c}\mathrm{2}\\ \mathit{x}\end{array}\mathrm{)}\mathrm{(}{\mathit{\sigma }}_{\mathit{y}}^{\mathrm{2}}\mathrm{-}\mathit{\epsilon ̅}\begin{array}{c}\mathrm{2}\\ \mathit{y}\end{array}\mathrm{)}}}\mathit{,}$(7)where is the mean measurement uncertainty of each data set. For calculating the mean and standard deviation, only data within each time lag bin are considered. Averaging the UDCF over the number of data pairs, M, yields the DCCF $\mathrm{(}\mathrm{DCCF}\mathrm{=}\frac{\mathrm{1}}{\mathit{M}}\sum {\mathrm{UDCF}}_{\mathit{ij}}$) and an uncertainty associated with each time lag can be determined as ${\mathit{\sigma }}_{\mathrm{DCCF}}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathit{M}\mathrm{-}\mathrm{1}}\sqrt{\begin{array}{c}\sum \end{array}{\mathrm{[}{\mathrm{UDCF}}_{\mathit{ij}}\mathrm{-}\mathrm{DCCF}\mathrm{(}\mathit{\tau }{\mathrm{)}}^{\mathrm{]}}}^{\mathrm{2}}}\mathit{.}$(8)Positive DCCF values indicate a positive correlation with an average time shift τ, while negative values imply anti-correlation.

Using the DCCF, Fuhrmann et al. (2014) searched for possible correlations between the 3 mm radio waveband and the Fermi γ-ray observations for the F-GAMMA blazar sample. Their study has shown the presence of significantly correlated 3 mm/γ-ray variability in an average sense (whole sample). Specifically, for PKS 1502+106 a 3 mm/γ-ray correlation is detected at a significance level above 99%, with the mm-radio emission lagging behind γ rays by 14 ± 11 days. The aforementioned difference in the time of arrival translates into a de-projected relative separation of the corresponding emitting regions of 2.1 pc, with the γ-ray active region upstream of the 3 mm (or 86.24 GHz) core.

Here, we enhance these findings and quantify the correlated variability beyond the mm/sub-mm band, using all available radio observations (see Fig. 1 for the light curves). The DCCFs have been calculated between the 15 GHz light curve of the OVRO 40-m blazar monitoring program and all available light curves from 226.5 down to 2.64 GHz. For each DCCF, we employed a time lag bin width of 30 days within a window of ± 750 d. The better sampling of the OVRO 15 GHz light curve led to its selection as reference. Hereafter, we refer to the cross-correlated light curves as the reference (OVRO 15 GHz), and to all the rest as the target. A positive lag τ_DCCF implies that the target waveband lags behind the reference frequency by τ_DCCF, measured in days. We use a weighted Gaussian fit of the form $\mathit{G}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}\mathit{\alpha }\hspace{0.17em}\exp \left[\frac{\mathrm{-}\mathrm{(}\mathit{x}\mathrm{-}\mathit{\mu }{\mathrm{)}}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}}\right]$, to obtain the DCCF peak and its position. Since the value of interest is the peak of the DCCF, the Gaussian fit is performed in the time lag range that ensures better peak determination. The mean of the Gaussian is the value of time lag τ_DCCF reported in Table 2. The left-hand panels of Figs. 2 and 3 show all DCCFs with respect to our selected reference frequency along with the Gaussian fit.

To estimate the time-lag uncertainty we employ the model-independent Monte Carlo (MC) approach introduced by Peterson et al. (1998; see also Raiteri et al. 2003). The two main sources of uncertainty in cross-correlation analyses – in addition to limited data trains and a consequently limited number of events – are (i) flux uncertainties; and (ii) uncertainties associated with the uneven sampling of the light curves. The method comprises two parts; the flux redistribution (FR), and a random sample selection (RSS), accounting for (i) and (ii) respectively. The algorithm is essentially as follows:

• 1.

  Introduce Gaussian deviates to the data, constrained by the maximum flux density uncertainty.

• 2.

  Use a re-sampling with replacement scheme to randomly select data points and create an equally long light curve; in other words, a bootstrap-like procedure.

• 3.

  Calculate the DCCF and determination of the time lag.

• 4.

  After a given number of MC realizations, obtain the cross-correlation peak distribution (CCPD; see right panels of Figs. 2 and 3).

The uncertainty, ± Δτ_68%, or simply ± Δτ, is obtained directly from the CCPD, and corresponds to the time lag interval of (τ_median − Δτ) and (τ_median + Δτ), such that 68% of the realizations yield results within this interval. This estimate is equivalent to the 1σ error, in case of a Gaussian distribution (Peterson et al. 1998). Here, 5000 MC realizations were performed.

3. Results

3.1. Flare parameters

The flare amplitudes, S_m, show an increasing trend with frequency, up to about 43 GHz whereafter they start decaying again, as is seen for 86.24, 142.33, and 226.5 GHz (see Fig. 4). As discussed later, there is a clear tendency for the flare to be visible earlier at higher frequencies, an effect attributed to synchrotron opacity (see also Fig. 5). The flare rising times, though they initially increase with frequency, at frequencies higher than 32.00 GHz show a decreasing trend (see Table 1). Finally, flare decay times seem to be following a slightly increasing trend with observing frequency in the whole range employed here (Table 1).

3.2. Frequency-dependent time lags

The DCCFs shown in Figs. 2 and 3 confirm the presence of a correlation between the outburst seen across observing frequencies. In essentially all cases, a single strong peak is clearly visible. In the following, we use the obtained time lag where the DCCF peaks (Table 2), referenced to the highest-frequency, well-sampled light curve at 142.33 GHz. The final τ_DCCF are reported in Table 3.

In Fig. 5 the delays obtained with both methods are shown. The GP regression and the DCCF analysis yield conclusive results that agree very well. A systematic trend is clearly seen, in the sense that the time lags become larger towards lower frequencies. This trend is described well by a power-law. In order to obtain its index, the frequency-dependent delays are fitted by a power-law of the form aν^{− 1 /k_r} using a least-squares procedure. The results for index k_r are for the GP regression k_{r, GP} = (1.4 ± 0.1) and k_{r, DCCF} = (1.6 ± 0.1) for the DCCF results. Both are shown, along with the best-fit curves, in Fig. 5.

The rather conservative uncertainty estimates of the FR/RSS method (Peterson et al. 1998) and the good agreement between τ_DCCF and τ_GP allow an average time lag, ⟨τ⟩, and its standard error to be obtained (Table 3). Hereafter, this average value is used for all calculations, and all uncertainties are formally propagated. The longest delay of (205 ± 2.0) days is seen between 142.33 and 2.64 GHz. A decreasing trend is clearly visible towards higher frequencies with a time lag of (136.9 ± 8.1) d at 4.85 GHz and ⟨τ⟩ = (96.2 ± 3.8) d at 8.35 GHz. Between the two highest frequencies of 142.33 GHz and 86.24 GHz, it is ⟨τ⟩ = (24.8 ± 0.2) days (see Table 3).

3.3. Time-lag core shifts and synchrotron opacity structure

Like most blazars, PKS 1502+106 features a flat-spectrum core, indicative of the superposition of optically-thick synchrotron components (e.g., Kaiser 2006). In addition to the scenario that the observed VLBI core is the first stationary hotspot (recollimation shock) of the jet (Daly & Marscher 1988; Gómez et al. 1995), a widely accepted interpretation is that this core is the point of the continuous flow where the photons at a given frequency have sufficient energy to escape the opacity barrier. In this scenario, the core is the surface where the optical depth is close to unity at a given frequency (e.g., Königl 1981). The standard relativistic jet paradigm (Blandford & Königl 1979; Königl 1981) predicts that the apparent position of this unit-opacity surface depends on observing frequency. Therefore, according to the interpretation of the VLBI core discussed above, its position must also be frequency-dependent. This “core-shift effect” was first observed by Marcaide & Shapiro (1984) and ever since has been studied using multi-frequency VLBI (e.g., Lobanov 1998; Kovalev et al. 2008; Sokolovsky et al. 2011; Pushkarev et al. 2012). The importance of core-shift measurements lies in the fact that they can provide critical insights into the structure and physics of ultracompact jets, such as the position of the core at each frequency relative to the jet base, and the magnetic field in different locations along the relativistic flow (e.g., Lobanov 1998; Hirotani 2005).

Relative timing of total flux density outbursts, as seen in single-dish radio observations, presents a viable alternative to VLBI core shifts. Returning our attention to the standard jet model, here, a disturbance originating at the jet base will propagate outwards along the jet. As this disturbance (or shock) moves away from the base it crosses the sequence of physically separated regions where the jet flow becomes optically thin – that is, the VLBI cores – as seen at a given sequence of observing frequencies. When the disturbance reaches a given unit-opacity surface for a given frequency ν_i, the synchrotron flare emission most likely becomes optically thin at ν_i. This appears in single-dish light curves as a total intensity outburst, while in VLBI maps these traveling shocks manifest as regions of enhanced emission after the core, the so-called VLBI knots (Marscher & Gear 1985). Since the times of maxima in single-dish light curves correspond to the flare emission becoming optically thin at the unit-opacity surface, for a given frequency, it follows that these times are frequency-dependent as well, and carry information on the opacity of each core. In other words, through relative timing of the same flaring event seen at different bands, we can obtain a tomography of the jet in terms of nuclear opacity, similar to VLBI core-shift measurements (e.g., Kudryavtseva et al. 2011; Kutkin et al. 2014).

An important consideration to be taken into account is that the method is based on the absence of any optically thick component downstream of the core. If there was any, it would sooner or later induce an opacity-driven enhancement of flux density due to its expansion and subsequent drop of density at a region further away from the core, the position of which we try to constrain. Essentially, the emission needs to become optically thin when it reaches each respective unit-opacity surface. In this respect and given the spectra of VLBI components presented by Karamanavis et al. (2016), component C3, used for our calculations, is optically thin between 43 and 86 GHz at least (spectral index α_{43/86 GHz} = −0.8 ± 0.3).

A “time-lag core-shift” between two frequencies ν₁ and ν₂, with ν₂ > ν₁, can be defined as follows, based on Lobanov (1998; see also Kudryavtseva et al. 2011) ${\mathrm{\Omega }}_{\mathrm{r}\mathit{\nu }}\mathrm{=}\mathrm{4.85}\mathrm{\times }{\mathrm{10}}^{-9}\frac{\mathit{\mu }\hspace{0.17em}\mathrm{\Delta }\mathit{t}\hspace{0.17em}{\mathit{D}}_{\mathrm{L}}}{\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}{\mathrm{)}}^{\mathrm{2}}}\left(\frac{{\mathit{\nu }}_{\mathrm{1}}^{\mathrm{1}\mathit{/}{\mathit{k}}_{\mathrm{r}}}\mathrm{·}{\mathit{\nu }}_{\mathrm{2}}^{\mathrm{1}\mathit{/}{\mathit{k}}_{\mathrm{r}}}}{{\mathit{\nu }}_{\mathrm{2}}^{\mathrm{1}\mathit{/}{\mathit{k}}_{\mathrm{r}}}\mathrm{-}{\mathit{\nu }}_{\mathrm{1}}^{\mathrm{1}\mathit{/}{\mathit{k}}_{\mathrm{r}}}}\right)\mathit{,}$(9)in units of pc GHz for k_r = 1, where

• μ

  is the VLBI jet proper motion inmas yr^-1;

• Δt

  is the observed time lag between the frequency pair ν₁ and ν₂ in yr;

• D _L

  is the luminosity distance in pc;

• z

  is the redshift;

• ν ₁ , ν ₂

  are the first and second observing frequency, both expressed in GHz; and

• k _r

  is the power-law index obtained by fitting the cross-band delays.

Here, k_r = [ (3−2α)m + 2n−2 ] /(5−2α) with m and n denoting the index of the power-law dependence of the magnetic field B(r) ∝ r^{− m} and the electron number density N_e(r) ∝ r^{− n}, at distance r along the jet, respectively. Finally, α denotes the optically thin spectral index. In case of external density gradients and/or free-free absorption, k_r> 1. For synchrotron self-absorption (SSA) dominated opacity and equipartition between particle energy and magnetic field energy, k_r = 1 for m = 1 and n = 2 independently of α. The angular offset in units of mas is given by Δr = μ⟨τ⟩. That is, the distance traveled by the disturbance within the time interval ⟨τ⟩ (Table 3). Using Ω_rν and following Lobanov (1998), Hirotani (2005), and Kudryavtseva et al. (2011), we obtain the distance between the core and the base of the assumed-conical jet. The frequency-dependent distance of the core – at any frequency ν – is given by ${\mathit{r}}_{\mathrm{core}}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}\mathrm{=}\frac{{\mathrm{\Omega }}_{\mathrm{r}\mathit{\nu }}}{\sin \mathit{\theta }}{\mathit{\nu }}^{\mathrm{-}\mathrm{1}\mathit{/}{\mathit{k}}_{\mathrm{r}}}\mathit{.}$(10)We employ Eqs. (9) and (10) to infer the separation of the radio core at each frequency from the jet base. For this, we use the average k_r = (1.5 ± 0.1) calculated based on the DCCF analysis and the Gaussian process regression. The two agree very well with each other indicating that k_r> 1. This implies the possible presence of jet external density gradients and/or foreground free-free absorption. Other parameters used are μ = 0.09 mas yr^-1 which corresponds to the proper motion of component C3 at 86 GHz, responsible for the flare, and θ = 2.6° the critical viewing angle (Karamanavis et al. 2016).

Our results reveal the nuclear opacity structure of PKS 1502+106 (see Table 4). The largest distance of about 10 pc is seen for the 2.64 GHz core, as expected, and distances diminish towards higher frequencies. Between 23.05 and 43.05 GHz the light curves become slightly more sparsely sampled because, at those bands, weather effects and poorer system performance of the Effelsberg 100-m telescope play an increasingly important role. At the highest frequency of 86.24 GHz, the core lies at a distance of (4.0 ± 1.1) pc away from the jet base.

3.4. Equipartition magnetic field

We obtain an expression for the magnetic field at 1 pc from the jet base, B_{1 pc}, by evaluating Eq. (43) from Hirotani (2005; see also O’Sullivan & Gabuzda 2009; Kutkin et al. 2014). Under equipartition between the energies of the magnetic field and that of radiating particles, and with a spectral index α = −0.5, this yields ${\mathit{B}}_{\mathrm{1}\hspace{0.17em}\mathrm{pc}}\mathrm{\approx }\mathrm{0.014}{\left[\frac{{\mathrm{\Omega }}_{\mathrm{r}\mathit{\nu }}^{\mathrm{3}{\mathit{k}}_{\mathrm{r}}}\hspace{0.17em}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}{\mathrm{)}}^{\mathrm{2}}\ln \mathrm{(}{\mathit{\gamma }}_{\max }\mathit{/}{\mathit{\gamma }}_{\min }\mathrm{)}}{{\mathit{\delta }}^{\mathrm{2}}\mathit{\phi }{\sin }^{\mathrm{3}{\mathit{k}}_{\mathrm{r}}\mathrm{-}\mathrm{1}}\mathit{\theta }}\right]}^{\mathrm{1}\mathit{/}\mathrm{4}}\mathrm{\left(}\mathrm{G}\mathrm{\right)}\mathit{,}$(11)where

• γ _max , γ _min

  are the maximum and minimum Lorentz factors of the emittinge⁻,

• δ

  is the Doppler factor,

• φ

  is the half-opening angle of the jet, and

• θ

  is the viewing angle.

Then, the magnetic field strength at the core is given by ${\mathit{B}}_{\mathrm{core}}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}\mathrm{=}{\mathit{B}}_{\mathrm{1}\hspace{0.17em}\mathrm{pc}}\hspace{0.17em}{\mathit{r}}_{\mathrm{core}}\mathrm{(}\mathit{\nu }{\mathrm{)}}^{-1}\mathit{.}$(12)We adopt ln(γ_max/γ_min) = 10 and employ the values for the Doppler factor and half-opening angle, deduced by Karamanavis et al. (2016), of δ = 10 and φ = (1.90 ± 0.25)°, respectively corresponding to components traveling within the inner jet of PKS 1502+106. Resulting values for B_{1 pc} and B_core are summarized in Table 4. We also report a conservative range for these values, calculated by combining the uncertainties of the input parameters in a way that maximizes this range. Inferred magnetic field strengths are higher at higher frequencies, when approaching the jet base. The values range between 14 and 176 mG for the 2.64 GHz (outermost) and 86.24 GHz (innermost) cores, respectively.

The results presented here are in reasonable agreement with the core shift measurements of PKS 1502+106 derived by Pushkarev et al. (2012) who employ standard VLBI techniques to obtain a separation of ~8 pc between the 15 GHz VLBI core and the vertex of the jet. Some discrepancies can be seen in the physical parameters of the jet, which are inferred from VLBI core shifts, such as the magnetic field at a distance of 1 pc and in the 15 GHz core (B_{1 pc} and B_core), but these can be explained through the different kinematical parameters used. Specifically, Pushkarev et al. (2012) make use of the fastest non-accelerating, component’s apparent speed from data covering the period between August 1994 and September 2007, as reported by Lister et al. (2009). However, here we use the proper motion, μ, of the knot most likely to be responsible for the multi-wavelength flare of 2008/2010 (knot C3 in Karamanavis et al. 2016), whose apparent angular velocity – while traveling in the inner jet of PKS 1502+106 – is about half that used by Pushkarev et al. (2012).

4. Discussion

4.1. Localizing the γ-ray emission

Our ultimate goal of pinpointing the γ-ray emission region is reached by combining the findings of the present study with those of Fuhrmann et al. (2014), who deduce a relative distance between the γ-ray active region and the 86.24 GHz unit-opacity surface of about 2.1 pc (cf. Ramakrishnan et al. 2015). Here, based on the opacity-driven time lags across the eleven observing frequencies, the distance of the 86.24 GHz core from the jet base is (4.0 ± 1.1) pc. From that combination it follows that the γ-ray emission region is best constrained to a distance of (1.9 ± 1.1) pc away from the jet base. In a similar study, Max-Moerbeck et al. (2014) obtain a time lag of (− 40 ± 13) d between the leading γ rays and the lagging 15 GHz radio emission. Following their reasoning, this time lag translates into a relative distance of (2 ± 1) pc between the 15 GHz unit-opacity surface and the γ-ray site. Combining our r_core at 15 GHz (Table 4) with the aforementioned value, we obtain a separation of (1.8 ± 1.3) pc between the jet base and the γ-ray active region, comparable with our previous result based on the 86.24 GHz time lag. We reiterate that in the present work we use the VLBI kinematical parameters of the knot connected to the 2008/10 flare of PKS 1502+106 (C3, see Sect. 3.4). Finally, our results are in agreement with the upper limit of ≤ 5.9 pc reported by Karamanavis et al. (2016), based on mm-VLBI imaging.

At (1.9 ± 1.1) pc from the jet base, the high-energy emission originates well beyond the bulk of the BLR material of PKS 1502+106, placed at R_BLR ≈ 0.1 pc (Abdo et al. 2010). A γ-ray emission region placed so far away from the BLR almost negates the notion that BLR photons alone can be used as the target photon field for inverse-Compton up-scattering to γ rays, and as a result other or additional seed photon fields (e.g., torus, synchrotron self-Compton, SSC) need to be involved in the process.

4.2. A shock origin for the 2008/10 flare

Overall, the slopes extracted from the multi-wavelength light curves agree well with the expectations of the shock-in-jet model (Marscher & Gear 1985), indicating a shock origin for the 2008/2010 flare of PKS 1502+106. In this section, we discuss the observed light curve parameters (flare amplitudes and cross-band delays) in conjunction with the results of analytical shock-in-jet model simulations (Fromm et al. 2015). The dependence of flare amplitude on frequency is shown in Fig. 4 and the cross-frequency delays in Fig. 5. Flare amplitudes rise with frequency until they culminate at a frequency of about 43 GHz, after which the amplitude of the flare drops. The rise and decay wings follow a trend that is described well by a broken power law, in very good agreement with the shock model expectations. The slopes are ϵ_{flare amp rise} = (0.7 ± 0.1) for the rising part of the flare amplitude and ϵ_{flare amp decay} = (−0.3 ± 0.05) for the decaying part. The power-law slope of the observed time lags, ϵ_time−lag = −1 /k_r (see Sect. 3) is ϵ_time−lag ≈ −0.7 (see Fig. 5). This value is obtained by averaging over the two values arising from the GP regression and the DCCF analysis and is in accord with the pattern expected from the shock model. By comparing the slope of the frequency-dependent cross-band delays with the findings of Fromm et al. (2015), we can safely exclude the scenario of a decelerating jet. Our findings, in fact, suggest the presence of a certain degree of acceleration taking place within the relativistic jet of PKS 1502+106.

5. Conclusions

In this work we exploited the dense, long-term F-GAMMA radio light curves in the frequency range 2.64–142.33 GHz obtained with the Effelsberg 100-m and IRAM 30-m telescopes. The F-GAMMA data were complemented by light curves at 15.00 GHz (OVRO) and 226.5 GHz (SMA).

A detailed section was devoted to the two independent approaches used in order to quantify the observed flare of PKS 1502+106 in the time domain. These are: (i) a Gaussian process regression and (ii) a discrete cross-correlation function analysis. This is the first time that GP regression is used in the field of blazar variability and it has shown to perform well. It is a viable approach to the problem of modeling discrete, unevenly sampled, blazar light curves and extracting the relevant parameters.

The flare of PKS 1502+106 is a “clean”, isolated outburst seen across observing frequencies. Our findings suggest a frequency-dependence for the cross-band delays and the flare amplitudes. The flare amplitudes first increase – up to the frequency of ~43 GHz, where they peak – after which they decrease (Fig. 4). The times at which the flare peaks are characterized by increasing time delay towards lower frequencies (Fig. 5). The characteristic dependence of these parameters on the observing frequency can be approximated by power laws, indicative of shock evolution and pc-scale jet acceleration.

Through the observed opacity-driven time lags, the structure of PKS 1502+106 in terms of synchrotron opacity is deduced – that is, using a “time-lag core shift” method. The positions of the ten radio unit-opacity surfaces with respect to the jet base were derived, with distances in the range of (10.2 ± 1.2) pc to (4.1 ± 1.1) pc for the core between 2.64 and 86.24 GHz, respectively.

The frequency-dependent core positions also enable a tomography of the magnetic field along the jet axis. We obtained the equipartition magnetic field at the position of each core, B_core, and also at a distance of 1 pc from the jet base, B_{1 pc}. The former, between the 2.64 and 86.24 GHz unit-opacity surfaces, are in the range B_{core, 2.64} ~ 14 mG and B_{core, 86.24} ~ 176 mG. For the latter, our figures indicate an average value of ${}^{\mathrm{⟨}}\mathit{B}_{\mathrm{1}\hspace{0.17em}\mathrm{pc}}{}^{\mathrm{⟩}}\mathrm{=}\mathrm{(}\mathrm{174}\mathrm{\pm }\mathrm{55}\mathrm{)}$ mG.

The γ-ray emission region is constrained to (1.9 ± 1.1) pc away from the jet base, well beyond the bulk of the BLR material of PKS 1502+106. This yields a contribution of IR torus photons and/or SSC to the production of γ rays in PKS 1502+106, while almost negating the contribution of the accretion disk or BLR photons alone.

The flare complies with the typical shock evolutionary path and its overall behavior in the time domain is in accord with the shock-in-jet scenario. When seen in the light of the VLBI findings discussed by Karamanavis et al. (2016), our results corroborate the scenario that the flare of PKS 1502+106 was induced by a shock, seen in mm-VLBI images as component C3, travelling downstream of the core at 43/86 GHz maps. This traveling disturbance is associated with the multi-frequency flare seen from radio up to γ-ray energies.

---

Acknowledgments

V.K. is thankful to B. Boccardi, A. P. Lobanov, W. Max-Moerbeck, B. Rani, S. Larsson, and A. Karastergiou for informative and fruitful discussions. We thank the anonymous referee for useful suggestions. V.K., I.N., and I.M. were supported for this research through a stipend from the International Max Planck Research School (IMPRS) for Astronomy and Astrophysics at the Universities of Bonn and Cologne. Based on observations with the 100-m telescope of the MPIfR at Effelsberg and with the IRAM 30-m telescope. The single-dish mm observations were coordinated with the flux density monitoring conducted by IRAM, and data from both programs are included here. IRAM is supported by INSU/CNRS (France), MPG (Germany), and IGN (Spain). The OVRO 40-m Telescope Fermi Blazar Monitoring Program is supported by NASA under awards NNX08AW31G and NNX11A043G, and by the NSF under awards AST-0808050 and AST-1109911. The Submillimeter Array is a joint project between the Smithsonian Astrophysical Observatory and the Academia Sinica Institute of Astronomy and Astrophysics who also fund it. This research has made use of NASA’s Astrophysics Data System and the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with NASA.

References

All Tables

All Figures