Issue 
A&A
Volume 590, June 2016



Article Number  A48  
Number of page(s)  10  
Section  Extragalactic astronomy  
DOI  https://doi.org/10.1051/00046361/201527796  
Published online  10 May 2016 
What can the 2008/10 broadband flare of PKS 1502+106 tell us?
Nuclear opacity, magnetic fields, and the location of γ rays
^{1}
MaxPlanckInstitut für Radioastronomie,
Auf dem Hügel 69,
53121
Bonn,
Germany
email:
vkaramanavis@mpifr.de
^{2}
Instituto de Radio Astronomía Milimétrica,
Avenida Divina Pastora 7, Local
20, 18012
Granada,
Spain
^{3}
Harvard–Smithsonian Center for Astrophysics,
Cambridge, MA
02138,
USA
Received: 20 November 2015
Accepted: 11 March 2016
Context. The origin of blazar variability, as seen from radio up to γ rays, is still a heavily debated matter, and broadband flares offer a unique testbed for developing a better understanding of these extreme objects. One of these energetic outbursts was detected by Fermi/LAT in 2008 from the blazar PKS 1502+106. The outburst was observed from γ rays down to radio frequencies.
Aims. Through the delay between flare maxima at different radio frequencies, we study the frequencydependent position of the unitopacity surface and infer its absolute position with respect to the jet base. This nuclear opacity profile enables the jet’s magnetic field tomography. We also localize the γray emission region and explore the flare production mechanism.
Methods. The PKS 1502+106 radio flare is studied through singledish flux density measurements at 12 frequencies in the range 2.64 to 226.5 GHz. To quantify the flare, we employ both a Gaussian process regression and a discrete crosscorrelation function analysis.
Results. We find that the light curve parameters (flare amplitude and crossband delays) show a powerlaw dependence on frequency. Delays decrease with frequency, and the flare amplitudes increase up to about 43 GHz, and then decay. This behavior is consistent with the propagation of a shock downstream in the jet. The selfabsorbed radio cores are located between approximately ten and four pc from the jet base, and their magnetic field strengths range between 14 and 176 mG, at the frequencies 2.64 to 86.24 GHz. Finally, the γray active region is located at (1.9 ± 1.1) pc away from the jet base.
Key words: galaxies: active / galaxies: jets / radiation mechanisms: nonthermal / radio continuum: galaxies / gamma rays: galaxies
© ESO, 2016
1. Introduction
Blazars comprise the radioloud 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 highenergy 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 KelvinHelmholtz 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. Multifrequency flux density monitoring can shed light on blazar physics by accessing time and thus spatial scales which are unreachable even by interferometric imaging. Dense, longterm, multifrequency 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 synchrotronemitting region.
The powerful flatspectrum radio quasar (FSRQ) PKS 1502+106 (OR 103, S3 1502+10) at a redshift of z = 1.8385, D_{L} = 14 176.8 Mpc (AdelmanMcCarthy et al. 2008) harbours a supermassive black hole of ~10^{9}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 FermiGST multifrequency campaign covering the electromagnetic spectrum in its entirety, both with filledaperture 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 crossband delays. First, using the observed crossband and frequencydependent time lags, we estimate the core shift – that is the frequencydependent position of the core – with an approach other than through traditional multifrequency verylongbaseline interferometry (VLBI) measurements. With the timelag core shift at hand, an opacity profile of the source and, under the assumption of equipartition, the magnetic field strength along the synchrotronselfabsorbed 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 highenergy 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_{0} = 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 FermiGST AGN Multifrequency Monitoring Alliance (FGAMMA^{1}) program (Fuhrmann et al. 2007; Angelakis et al. 2010), which includes observations with the Effelsberg 100m telescope (EB) in eight frequency bands from 2.64 to 43.05 GHz, and with the IRAM 30m telescope (PV) at 86.24 and 142.33 GHz. Monthly observations from the EB and PV telescopes were performed in a quasisimultaneous 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) 40m blazar monitoring program^{2} 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 nonparametric 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 onedimensional 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 “bestfit 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^{3}. Here we employ the squared exponential kernel (1)where the hyperparameters are

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

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}) + σ^{2}I, where I is the identity matrix and σ^{2} 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 (2)Then, the probability (or evidence) of the training data, y, given the hyperparameters vector θ, is (3)and the log marginal likelihood is given by (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 (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^{4}. 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_{0}, 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 crossband delay, τ_{GP}. Additionally, the flare rising and decay times (t_{r} and t_{d}, respectively) are extracted from the times of the two fluxdensity 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.
Fig. 1 Gaussian process regression curves for the radio light curves in the range 2.64–226.50 GHz, shown here with the minimum flux density, S_{0}, subtracted. Observations are shown in blue, the posterior mean (prediction curve) in red, and the 95% confidence interval is the lighter redshaded area. 

Open with DEXTER 
Results of the Gaussian process regression.
2.2. The discrete crosscorrelation function
Putative correlated variability across observing bands, which are based on samplingratelimited time series, can be investigated and quantified using the discrete crosscorrelation function (DCCF, Edelson & Krolik 1988).
The crosscorrelation 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 (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 crosscorrelation 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 (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 ) and an uncertainty associated with each time lag can be determined as (8)Positive DCCF values indicate a positive correlation with an average time shift τ, while negative values imply anticorrelation.
Using the DCCF, Fuhrmann et al. (2014) searched for possible correlations between the 3 mm radio waveband and the Fermi γray observations for the FGAMMA 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 mmradio emission lagging behind γ rays by 14 ± 11 days. The aforementioned difference in the time of arrival translates into a deprojected 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.
Fig. 2 DCCFs (left) and CCPDs (right) between the 15 GHz light curve and those in the range 2.64 to 23.05 GHz. Solid lines represent the bestfit Gaussian functions. 

Open with DEXTER 
Fig. 3 DCCFs (left) and CCPDs (right) between the 15 GHz light curve and those in the range 32.00 to 226.5 GHz. Solid lines represent the bestfit Gaussian functions. 

Open with DEXTER 
Results of the correlation analysis at all frequencies.
Here, we enhance these findings and quantify the correlated variability beyond the mm/submm 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 40m 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 crosscorrelated 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 , 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 lefthand panels of Figs. 2 and 3 show all DCCFs with respect to our selected reference frequency along with the Gaussian fit.
To estimate the timelag uncertainty we employ the modelindependent Monte Carlo (MC) approach introduced by Peterson et al. (1998; see also Raiteri et al. 2003). The two main sources of uncertainty in crosscorrelation 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 resampling with replacement scheme to randomly select data points and create an equally long light curve; in other words, a bootstraplike procedure.

3.
Calculate the DCCF and determination of the time lag.

4.
After a given number of MC realizations, obtain the crosscorrelation peak distribution (CCPD; see right panels of Figs. 2 and 3).
Fig. 4 Frequency dependence of the flare amplitudes. Solid lines show the wings of the fitted broken power law, while the values denote their slopes. 

Open with DEXTER 
Fig. 5 Frequency dependence of the time lags with respect to the data at 142.33 GHz. Red circles denote values resulting from the GP regression and blue squares those from the DCCF analysis. Solid lines represent the bestfit power laws. 

Open with DEXTER 
Frequencydependent time lags and “timelag core shifts” with respect to the observations at 142.33 GHz.
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. Frequencydependent 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 highestfrequency, wellsampled 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 powerlaw. In order to obtain its index, the frequencydependent delays are fitted by a powerlaw of the form aν^{− 1 /kr} using a leastsquares 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 bestfit 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. Timelag core shifts and synchrotron opacity structure
Like most blazars, PKS 1502+106 features a flatspectrum core, indicative of the superposition of opticallythick 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 unitopacity surface depends on observing frequency. Therefore, according to the interpretation of the VLBI core discussed above, its position must also be frequencydependent. This “coreshift effect” was first observed by Marcaide & Shapiro (1984) and ever since has been studied using multifrequency VLBI (e.g., Lobanov 1998; Kovalev et al. 2008; Sokolovsky et al. 2011; Pushkarev et al. 2012). The importance of coreshift 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 singledish 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 unitopacity surface for a given frequency ν_{i}, the synchrotron flare emission most likely becomes optically thin at ν_{i}. This appears in singledish light curves as a total intensity outburst, while in VLBI maps these traveling shocks manifest as regions of enhanced emission after the core, the socalled VLBI knots (Marscher & Gear 1985). Since the times of maxima in singledish light curves correspond to the flare emission becoming optically thin at the unitopacity surface, for a given frequency, it follows that these times are frequencydependent 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 coreshift measurements (e.g., Kudryavtseva et al. 2011; Kutkin et al. 2014).
Synchrotron opacity structure and magnetic field tomography of the jet of PKS 1502+106.
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 opacitydriven 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 unitopacity 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 “timelag coreshift” between two frequencies ν_{1} and ν_{2}, with ν_{2} > ν_{1}, can be defined as follows, based on Lobanov (1998; see also Kudryavtseva et al. 2011) (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 ν_{1} and ν_{2} in yr;

D
_{L}
is the luminosity distance in pc;

z
is the redshift;

ν
_{1}
,
ν
_{2}
are the first and second observing frequency, both expressed in GHz; and

k
_{r}
is the powerlaw index obtained by fitting the crossband delays.
Here, k_{r} = [ (3−2α)m + 2n−2 ] /(5−2α) with m and n denoting the index of the powerlaw 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 freefree absorption, k_{r}> 1. For synchrotron selfabsorption (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 assumedconical jet. The frequencydependent distance of the core – at any frequency ν – is given by (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 freefree 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 100m 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 (11)where

γ
_{max}
,
γ
_{min}
are the maximum and minimum Lorentz factors of the emittinge^{−},

δ
is the Doppler factor,

φ
is the halfopening angle of the jet, and

θ
is the viewing angle.
Then, the magnetic field strength at the core is given by (12)We adopt ln(γ_{max}/γ_{min}) = 10 and employ the values for the Doppler factor and halfopening 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 nonaccelerating, 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 multiwavelength 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 unitopacity surface of about 2.1 pc (cf. Ramakrishnan et al. 2015). Here, based on the opacitydriven 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, MaxMoerbeck 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 unitopacity 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 mmVLBI imaging.
At (1.9 ± 1.1) pc from the jet base, the highenergy 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 inverseCompton upscattering to γ rays, and as a result other or additional seed photon fields (e.g., torus, synchrotron selfCompton, SSC) need to be involved in the process.
4.2. A shock origin for the 2008/10 flare
Overall, the slopes extracted from the multiwavelength light curves agree well with the expectations of the shockinjet 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 crossband delays) in conjunction with the results of analytical shockinjet model simulations (Fromm et al. 2015). The dependence of flare amplitude on frequency is shown in Fig. 4 and the crossfrequency 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 powerlaw 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 frequencydependent crossband 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, longterm FGAMMA radio light curves in the frequency range 2.64–142.33 GHz obtained with the Effelsberg 100m and IRAM 30m telescopes. The FGAMMA 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 crosscorrelation 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 frequencydependence for the crossband 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 pcscale jet acceleration.
Through the observed opacitydriven time lags, the structure of PKS 1502+106 in terms of synchrotron opacity is deduced – that is, using a “timelag core shift” method. The positions of the ten radio unitopacity 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 frequencydependent 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 unitopacity 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 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 shockinjet 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 mmVLBI images as component C3, travelling downstream of the core at 43/86 GHz maps. This traveling disturbance is associated with the multifrequency flare seen from radio up to γray energies.
Acknowledgments
V.K. is thankful to B. Boccardi, A. P. Lobanov, W. MaxMoerbeck, 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 100m telescope of the MPIfR at Effelsberg and with the IRAM 30m telescope. The singledish 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 40m Telescope Fermi Blazar Monitoring Program is supported by NASA under awards NNX08AW31G and NNX11A043G, and by the NSF under awards AST0808050 and AST1109911. 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
 Abdo, A. A., Ackermann, M., Ajello, M., et al. 2010, ApJ, 710, 810 [NASA ADS] [CrossRef] [Google Scholar]
 AdelmanMcCarthy, J. K., Agüeros, M. A., Allam, S. S., et al. 2008, ApJS, 175, 297 [NASA ADS] [CrossRef] [Google Scholar]
 Aigrain, S., Pont, F., & Zucker, S. 2012, MNRAS, 419, 3147 [NASA ADS] [CrossRef] [Google Scholar]
 Angelakis, E., Fuhrmann, L., Nestoras, I., et al. 2010, ArXiv eprints [arXiv:1006.5610] [Google Scholar]
 Angelakis, E., Fuhrmann, L., Marchili, N., et al. 2015, A&A, 575, A55 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Blandford, R. D., & Königl, A. 1979, ApJ, 232, 34 [NASA ADS] [CrossRef] [Google Scholar]
 Camenzind, M., & Krockenberger, M. 1992, A&A, 255, 59 [NASA ADS] [Google Scholar]
 Ciprini, S. 2008, ATel, 1650, 1 [NASA ADS] [Google Scholar]
 Daly, R. A., & Marscher, A. P. 1988, ApJ, 334, 539 [NASA ADS] [CrossRef] [Google Scholar]
 Edelson, R. A., & Krolik, J. H. 1988, ApJ, 333, 646 [NASA ADS] [CrossRef] [Google Scholar]
 Fromm, C. M., Fuhrmann, L., & Perucho, M. 2015, A&A, 580, A94 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fuhrmann, L., Zensus, J. A., Krichbaum, T. P., Angelakis, E., & Readhead, A. C. S. 2007, in The First GLAST Symposium, eds. S. Ritz, P. Michelson, & C. A. Meegan, AIP Conf. Ser., 921, 249 [Google Scholar]
 Fuhrmann, L., Larsson, S., Chiang, J., et al. 2014, MNRAS, 441, 1899 [NASA ADS] [CrossRef] [Google Scholar]
 Gibson, N. P., Aigrain, S., Roberts, S., et al. 2012, MNRAS, 419, 2683 [NASA ADS] [CrossRef] [Google Scholar]
 Gómez, J. L., Martí, J. M. A., Marscher, A. P., Ibanez, J. M. A., & Marcaide, J. M. 1995, ApJ, 449, L19 [NASA ADS] [CrossRef] [Google Scholar]
 Gurwell, M. A., Peck, A. B., Hostler, S. R., Darrah, M. R., & Katz, C. A. 2007, in From ZMachines to ALMA: (Sub)Millimeter Spectroscopy of Galaxies, eds. A. J. Baker, J. Glenn, A. I. Harris, J. G. Mangum, & M. S. Yun, ASP Conf. Ser., 375, 234 [Google Scholar]
 Hirotani, K. 2005, ApJ, 619, 73 [NASA ADS] [CrossRef] [Google Scholar]
 Ivezić, V., Connolly, A. J., VanderPlas, J. T., & Gray, A. 2014, Statistics, Data Mining, and Machine Learning in Astronomy, stu – student edition edn., Princeton Series in Modern Observational Astronomy (Princeton University Press), 544 [Google Scholar]
 Kaiser, C. R. 2006, MNRAS, 367, 1083 [NASA ADS] [CrossRef] [Google Scholar]
 Karamanavis, V. 2015, Ph.D. Thesis, University of Cologne [Google Scholar]
 Karamanavis, V., Fuhrmann, L., Krichbaum, T. P., et al. 2016, A&A, 586, A60 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Königl, A. 1981, ApJ, 243, 700 [NASA ADS] [CrossRef] [Google Scholar]
 Kovalev, Y. Y., Lobanov, A. P., Pushkarev, A. B., & Zensus, J. A. 2008, A&A, 483, 759 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kudryavtseva, N. A., Gabuzda, D. C., Aller, M. F., & Aller, H. D. 2011, MNRAS, 415, 1631 [NASA ADS] [CrossRef] [Google Scholar]
 Kutkin, A. M., Sokolovsky, K. V., Lisakov, M. M., et al. 2014, MNRAS, 437, 3396 [NASA ADS] [CrossRef] [Google Scholar]
 Larsson, S. 2012, Fermi & Jansky Proceedings – eConf C1111101 [arXiv:1207.1459] [Google Scholar]
 Lister, M. L., Kellermann, K. I., Vermeulen, R. C., et al. 2003, ApJ, 584, 135 [NASA ADS] [CrossRef] [Google Scholar]
 Lister, M. L., Cohen, M. H., Homan, D. C., et al. 2009, AJ, 138, 1874 [NASA ADS] [CrossRef] [Google Scholar]
 Lobanov, A. P. 1998, A&A, 330, 79 [NASA ADS] [Google Scholar]
 Lobanov, A. P., & Zensus, J. A. 2001, Science, 294, 128 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Marcaide, J. M., & Shapiro, I. I. 1984, ApJ, 276, 56 [NASA ADS] [CrossRef] [Google Scholar]
 Marscher, A. P., & Gear, W. K. 1985, ApJ, 298, 114 [NASA ADS] [CrossRef] [Google Scholar]
 MaxMoerbeck, W., Hovatta, T., Richards, J. L., et al. 2014, MNRAS, 445, 428 [NASA ADS] [CrossRef] [Google Scholar]
 Nestoras, I., Fuhrmann, L., Angelakis, E., et al. 2016, A&A, submitted [Google Scholar]
 O’Sullivan, S. P., & Gabuzda, D. C. 2009, MNRAS, 400, 26 [NASA ADS] [CrossRef] [Google Scholar]
 Pedregosa, F., Varoquaux, G., Gramfort, A., et al. 2011, J. Machine Learning Res., 12, 2825 [Google Scholar]
 Peterson, B. M., Wanders, I., Horne, K., et al. 1998, PASP, 110, 660 [NASA ADS] [CrossRef] [Google Scholar]
 Pushkarev, A. B., Hovatta, T., Kovalev, Y. Y., et al. 2012, A&A, 545, A113 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Raiteri, C. M., Villata, M., Tosti, G., et al. 2003, A&A, 402, 151 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ramakrishnan, V., Hovatta, T., Nieppola, E., et al. 2015, MNRAS, 452, 1280 [NASA ADS] [CrossRef] [Google Scholar]
 Rasmussen, C. E., & Williams, C. K. I. 2006, Gaussian Processes for Machine Learning (The MIT Press) [Google Scholar]
 Richards, J. L., MaxMoerbeck, W., Pavlidou, V., et al. 2011, ApJS, 194, 29 [NASA ADS] [CrossRef] [Google Scholar]
 Roberts, S., Osborne, M., Ebden, M., et al. 2012, Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, 371 [Google Scholar]
 Robertson, D. R. S., Gallo, L. C., Zoghbi, A., & Fabian, A. C. 2015, MNRAS, 453, 3455 [NASA ADS] [Google Scholar]
 Sokolovsky, K. V., Kovalev, Y. Y., Pushkarev, A. B., Mimica, P., & Perucho, M. 2011, A&A, 535, A24 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Spada, M., Ghisellini, G., Lazzati, D., & Celotti, A. 2001, MNRAS, 325, 1559 [NASA ADS] [CrossRef] [Google Scholar]
 Valtaoja, E., Terasranta, H., Urpo, S., et al. 1992, A&A, 254, 71 [NASA ADS] [Google Scholar]
 van der Laan, H. 1966, Nature, 211, 1131 [NASA ADS] [CrossRef] [Google Scholar]
 Williams, C. K. I., & Rasmussen, C. E. 1996, in Advances in Neural Information Processing Systems 8, eds. D. S. Touretzky, M. C. Mozer, & M. E. Hasselmo (MIT Press) [Google Scholar]
All Tables
Frequencydependent time lags and “timelag core shifts” with respect to the observations at 142.33 GHz.
Synchrotron opacity structure and magnetic field tomography of the jet of PKS 1502+106.
All Figures
Fig. 1 Gaussian process regression curves for the radio light curves in the range 2.64–226.50 GHz, shown here with the minimum flux density, S_{0}, subtracted. Observations are shown in blue, the posterior mean (prediction curve) in red, and the 95% confidence interval is the lighter redshaded area. 

Open with DEXTER  
In the text 
Fig. 2 DCCFs (left) and CCPDs (right) between the 15 GHz light curve and those in the range 2.64 to 23.05 GHz. Solid lines represent the bestfit Gaussian functions. 

Open with DEXTER  
In the text 
Fig. 3 DCCFs (left) and CCPDs (right) between the 15 GHz light curve and those in the range 32.00 to 226.5 GHz. Solid lines represent the bestfit Gaussian functions. 

Open with DEXTER  
In the text 
Fig. 4 Frequency dependence of the flare amplitudes. Solid lines show the wings of the fitted broken power law, while the values denote their slopes. 

Open with DEXTER  
In the text 
Fig. 5 Frequency dependence of the time lags with respect to the data at 142.33 GHz. Red circles denote values resulting from the GP regression and blue squares those from the DCCF analysis. Solid lines represent the bestfit power laws. 

Open with DEXTER  
In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.