© ESO, 2013

1. Introduction

Blazar S5 0716+714 (DA 237, HB89) is a BL Lac with a redshift of z = 0.3 ± 0.08 (Nilsson et al. 2008). It is located at high declination (+71) and is usually bright which makes it an ideal candidate for micro-variability observations. Correlated intraday variations in optical and radio frequencies have been reported for this source (Quirrenbach et al. 1991). It is also a gamma-ray emitting object associated with an FR II radio source (Wagner et al. 1996; von Montigny et al. 1995; Chen et al. 2008; Villata et al. 2008). Observations over a period of ten years indicate it experiences nearly continuous micro-variability activity (Humrickhouse & Webb 2008; Webb et al. 2010), with duty cycle of about 95.3%.

Attempts to find periodicities or significant correlations have been done with extensive variability data. Nesci et al. (2005) used long term light curves from 1953 to 2005 and reported finding evidence of a long-term variability consistent with a precessing jet. Raiteri et al. (2003) reported a  3.3 year period in their data. Examination for periodicity in the short-term variability time frames yielded high significance in peaks of 25 and 73 min in twenty-two well sampled micro-variability curves (Gupta et al. 2009).

Among the numerous microvariability studies of this object in the literature, Villata et al. (2000) performed a very densely sampled WEBT observation and obtained 635 data points over 72 h. They found that the steepest rise and the steepest decline were both roughly 0.002 magnitudes per minute over several hours with no apparent difference between the rise and decline rates. Wu et al. (2005) performed nearly simultaneous multi-frequency monitoring over a span of seven days and time series analysis of their microvariability observations yielded no repeatable periodicities common in the light curves. Wu et al. (2007) used an objective prism set up to obtain truly simultaneous multicolor observations over several nights. Montagni et al. (2006) monitored the source for microvariability between 1996 and 2003 finding the most rapid variations were on the order of 0.1 mag per hour over a period of two hours. They also found no difference between the rise and decline rates. Fourier transform and wavelet analysis were performed on twenty-three independent microvariability curves obtained at the SARA Observatory between 1998 and 2005, but analysis failed to yield any of the previously reported periods in the data (Dhalla & Webb 2010). In an attempt to determine if the variations in the SARA data more accurately represented by noise, Dhalla & Webb (2010) applied the statistical methods of Vaughan et al. (2003) which involved comparing the rms variations of the flux over sub-samples of the data. The results proved inconclusive and simulations showed that even if the variations were due to noise, more data was necessary to effectively determine the noise content of the microvariations.

Since both period and noise analysis techniques suggested that the individual microvariability curves were not of sufficient length to determine the nature of the microvariations, we organized an international campaign through the Whole Earth Blazar Telescope (WEBT) to observe S5 0716+714 over a three day period. The WEBT is a group of observatories that collaborate on blazar projects ranging from satellite back-up to campaigns such as this¹. We requested WEBT observers around the world observe S50716+714 during this period using standard photometric techniques, a common set of comparison stars, and in a common filter. We report here on the data acquired during this observation (Sect. 2), the resulting time series analysis (Sect. 3), and we interprete the results in Sect. 4.

2. Observations

We selected a February observation date when S5 0716+714 would be readily accessible the entire night from most northern hemisphere sites. The response to the call for observers through WEBT was excellent and plans were made to carry out the observations. The primary observing band was chosen to be R since most observatories own and regularly use that filter for microvariability observations. We divided up the observatories into longitude regions around the globe to help coordinate the observations. It was also decided that if there were multiple telescopes in each longitude range, then at least one telescope would be assigned a different filter so some simultaneous color information could be acquired as the observation progressed. However, the continuity of the R light curve was the highest priority in this campaign. In addition to the common filter system, we also selected four comparison stars from the sequence of Villata et al. (1998). We chose to use stars 3, 4, and 5 as comparison stars and measuring the magnitude of star 6 as a check star in addition to the object. Stars 1 and 2 frequently are overexposed due to their brightness and cause problems when trying to get accurate photometry, so they were left out of the comparison sequence. Table 1 lists the longitude zone in Col. 1, the code assigned to each observatory in Col. 2, the location of the observatory and observatory name and longitude in Cols. 3–5 respectively. The telescope aperture and the filters used are in Cols. 6 and 7.

	Fig. 1 Raw light curve of S5 0716+714 obtained by the compilation of some of the high quality data by major contributors. The light curves contributed from each observatory is plotted together with different symbols indentifying the observatory according to the codes given in Table 1.
Open with DEXTER

Figure 1 shows the raw observations plotted together. The observations covered the time period between JD 2 454 886.1 (2/23/2009) and JD 2 454 889.5 (2/26/2009). There was nearly continuous coverage between JD 2 454 887.3 and 2 454 888, with overlap from several observatories during many time intervals. The code given in Table 1 is the key to the observatories responsible for each data segment on the plot. The data were reduced at each individual observatory and sent as magnitude files for final compilation and analysis. Overlaying each contributed light curve revealed some offsets due to instrumental/filter differences, but most of the data was of sufficient quality and had sufficient overlap so that minor zero-point adjustments could be made to obtain a consistent continuous light curve. In all cases, the data with the lowest noise and longest overlap with other data sets were used to determine the offsets for the other light curves. Exposure times for individual images ranged from 30 s to 120 s depending on the observatory and telescope.

	Fig. 2 Flux curve of 0716+714 smoothed using the algorithm discussed in the text.
Open with DEXTER

In order to prepare the light curve for time series analysis, we implemented a smoothing algorithm. The data were smoothed based on the assumption that given a time scale as short as two minutes, a data point cannot be too different from the average of the previous and the following point using the algorithm that if any data point is x_i > (x_i − 1 + x_i + 1)/2 + 0.005 then, x_i = (x_i − 1 + x_i + 1)/2. The resulting smoothed data set consisted of 2613 high quality data points. The smoothing algorithm only affected timescales on the order of a few minutes, much shorter than any possible periods we could find from a time series analysis. In order to ensure there was no unexpected bias introduced by the smoothing, we also preformed all of the time series analysis reported below on the unsmoothed data the results were identical to those performed on the smoothed data in the frequencies of interest. The magnitudes were converted to flux using standard flux conversions for the R filter as given by Johnson (1966) using a redshift of 0.30 and Galactic absorption of 0.031 mag. Figure 2 shows the complete smoothed flux curve. Smoothing kept the major trends of light curve intact, while cutting out some of the high frequency noise in the data. The total length of the light curve was 78.88 h.

We analyzed individual segments of the flux curve to determine the maximum climb and decline rates. Twelve individual rapid excursions were noted in the flux curve and we fit a line to each of those segments to determine the maximum climb rates and decline rates, concentrating on segments which had a large number of data points. The fastest rate was an increase of 0.089 mag per hour over a range of 0.15 mag. The correlation coefficient for that fit was r² = 0.997 and it contained 99 data points. Table 2 lists the slopes and correlation coefficients for each of the segments we examined. Column 1 is the segment index, Col. 2 the start and finish time of the segment in hours, and Col. 3 the slope in magnitudes per hour. We calculated the correlation coefficient for each fit and listed them in Col. 4 along with the number of points in Col. 5. Column 6 gives the probability that a random sample would show such a large correlation coefficient. The final column denotes whether the slope is a rise or a decline in magnitude. The average decline rate was 0.042 mag/h (standard deviation of 0.022) while the average rise rate was 0.043 with a standard deviation of 0.028. Thus overall, the rise and decline rates are similar in this segment of light curve. Although the rates are different, the fact that the slopes for the rise and decline are the same agree with the results found by Villata et al. (2000) and Montagni et al. (2006).

3. Time series analysis

3.1. Fourier transform analysis

We performed Fourier transform analysis on the entire smoothed light curve by removing the linear trend of slope  −2.5 × 10^-5 mJy/h and using a discrete Fourier transform (DFT) algorithm (Deeming 1975). The results of this analysis are shown in Fig. 3 and Table 3. The DFT results yielded some of the large features at periods of 40.00, 21.05 and 13.19 hours which correspond to the first, second and fourth peaks respectively in Fig. 3. Some of the periods which are well above the noise level are listed in Col. 4 of Table 3. The corresponding time-scales in the rest frame in Col. 5 were calculated using (1)where D = 1/Γ(1 − βcosθ) is Doppler factor with βc bulk speed and Γ bulk Lorentz factor. θ is the orientation of the jet axis with respect to the line of sight and z is the red-shift of the object. Using the values Γ = 17 , θ = 2.6° and z = 0.3, a Doppler factor of 21.32 was used to transform the time-scales in the rest frame (see Celotti & Ghisellini 2008).

	Fig. 3 Power spectrum of the smoothed Light curve. The horizontal line at low power indicates the maximum power level of simulated random noise light curves. The numbers associated with the peaks in the power spectrum correspond to the frequencies and periods listed in Table 3.
Open with DEXTER

The horizontal line across the bottom indicates the average power of 100 light curves generated in Interactive Data Language (IDL) assuming Gaussian distributed noise with sigma equal to the sigma of the flux curve. The peaks in the S5 0716+714 DFT are thus many sigmas above the noise level. We then pre-whitened the flux curve by subtracting the derived period, adjusting the amplitude and phase for best fit, and then reanalyzed the resulting light curve using the DFT. The amplitude of peaks in the DFT did not decrease drastically as we pre-whitened the data, indicating that the periods are not indicative of a true period running through the extended data set. Neither did we find the periods at approximately 25 and 73 min in this flux curve seen previously by Gupta et al. (2009). Thus we could not confirm any of the previously detected periods seen in the S50716 + 714 microvariability curves or propose significant new periods. We attribute the seemingly significant peaks seen in the DFT to particular features in pieces of the flux curve, not cyclical oscillations throughout the entire 72-h observation.

	Fig. 4 Wavelet transform of the smoothed light curve linearized by removing a slope of  − 2.5 × 10-5 mJy/h and a Y-intercept of flux 0.01 mJy. Morlet kernel of 6th order was used from the wavelet application in IDL.
Open with DEXTER

We performed wavelet analysis (Torrence & Compo 1998) using the wavelet application in IDL to compute the wavelet transform of the data and then compared with the DFT results. The Morlet kernel of order 6 was used to do the analysis. Figure 4 shows the resulting wavelet transform. The peak of the wavelet transform corresponds to a sinusoidal oscillation with a period of 30.72 h, but is clearly only significant in the center of the data run, fitting extremely poorly at the start of the light curve. The inconclusiveness of both the DFT and wavelet analysis, along with not reproducing any of the previously reported periods in this object led us to examine the flux curve in terms of a noise model.

3.2. Noise analysis

Dhalla & Webb (2010) used the analysis methods of Vaughan et al. (2003) to examine twenty-one single-night microvariability light curves of S5 0716+714. This analysis consisted of dividing the flux curves up into individual bins and calculating the rms deviation within the each bin given by: (2)According to Vaughan et al. (2003), if light curves are strictly the result of Gaussian noise processes, each independent flux curve would be one realization of the underlying stochastic process. If the process is a stationary noise process, the realizations should exhibit similar statistical properties e.g. a linear relationship between excess rms and average flux. Dhalla & Webb (2010) found that the available microvariability curves were too short to reliably recover the noise characteristics of the data and concluded that the individual microvariability curves needed to be longer to determine whether the variations are the result of a deterministic stochastic process. The WEBT observation presented here should be long enough to identifiy the nature of the microvariations if in fact the microvariability is due to purely noise processes. We repeated the same noise analysis reported by Dhalla & Webb (2010) with the WEBT 72-h flux curve of S5 0716+714. The data was binned into 45 min bins and only bins with a minimum of 20 data points were used. Figure 5 shows the rms vs. average flux plotted for the bins and it fails to show the expected linear relationship between rms and average flux if the the variations were due to a totally Gaussian noise process. The correlation of the best fit line to the data is 0.0003. In a further effort to deduce the noise content of the data, we re-plotted the DFT in log–log space in Fig. 6. The DFT plotted in this way is best fit to a slope of 1/f² noise, but there are several prominent features at low frequency which deviates from a simple noise distribution.

	Fig. 5 rms variations versus average flux for the entire 3-day observation.
Open with DEXTER

The result of not verifying any of the previously detected periods and the inconclusiveness of the noise analysis has led the authors propose a new model for the interpretation of microvariability. We propose the variations are due to stochastic pulses generated by a shock encountering a turbulent jet flow. This model is described in the following section then fit to the 72-h microvariability curve.

	Fig. 6 log–log plot of the power spectrum of 0716+714.
Open with DEXTER

4. Modeling the microvariations

Since time series analysis of the microvariability flux curve did not show repeatable periodicities nor a strictly noise characteristics, and since close inspection of many microvariability curves show that the microvariations can be resolved into pulses or shots, we investigate a model where individual synchrotron cells are energized by a plane shock propagating down the jet and result in an increase in flux resembling a pulse.

4.1. Turbulence in the plasma jet

The most natural model for the generation of stochastic synchrotron cells is the presence of turbulence in the jet. Jones (1988) numerically simulated a turbulent relativistic jet using nested cells to investigate the circular de-polarization properties of the VLBI radio jets. Marscher et al. (1992) investigated a turbulent model of the plasma in the radio jets discussed the interpretation of radio flickering seen in jets as a result of a shock encountering this turbulence. We are following these seminal works by assuming a turbulent jet is responsible for the microvariations and that as the shock accelerates particles in each cell, the particles cool by synchrotron emission. The individual vortices each could have different densities, sizes and magnetic field orientations. We further assume that as the strong shock hits each cell, rapid particle acceleration and subsequent cooling by synchrotron emission produces a pulse in the flux curve. The convolution of these individual pulses leads to the observed microvariability. The formation of the turbulent cells is a stochastic process, thus producing no strict periodicities in the data. We suppose that by de-convolving the microvariability curve into individual pulses, we can gain understanding into the underlying turbulence in the jet. To do this, we need to have a reasonable idea of the pulse profiles expected from the turbulent cells.

4.2. Calculation of the synchrotron pulse profiles

Lehto (1989) first investigated shot models to describe the 1/f nature of X-ray light curves. He studied shots that had a delta function rise and exponential decline and was able to show that a random distribution of such shots yields a 1/f power spectrum. The shot profiles studied by Lehto (1989) did not resemble the microvariability seen in 0716+714, but visual inspection of over 100 microvariability light curves between 1989 and 2009 (Montagni et al. 2006; Humrickhouse & Webb 2008) led Webb et al. (2010) to propose that the microvariability curves were indeed shots, but shots with a light curve profile given by Kirk et al. (1998; hereafter KRM). KRM calculated the particle acceleration in the shock front for various magnetic field orientations and particle densities assuming a plane shock encounters a cylindrical density enhancement. In this model, the ratio of the acceleration time t_acc and escape time t_esc, in addition to constraints on the cooling length L control the pulse shape. The amplitude is given by the parameter Q and is related to the magnetic field strength B and orientation θ in addition to the enhanced electron density. The profiles described by KRM look similar to the microvariability pulse profiles seen in the 0716+714 microvariability curves and microvariability curves of other objects (Bhatta et al. 2011). The general picture of this model is schematically shown in Fig. 7. We re-evaluated the KRM model for our case assuming that the electrons accelerated by the shock obey the particle distribution function given by Ball &Kirk (1992). (3)with (4)where N is the number density of the electrons in the energy space represented by Lorentz factor γ. t_acc and t_esc are the time of acceleration and time of escape respectively. The term β_sγ² in Eq. (2) represents energy loss of an electron due to synchrotron emission where σ_T, μ₀, m_e and c are Thomson’s scattering cross-section, permeability of the free space, mass of an electron and the speed of light respectively.

	Fig. 7 Example of geometry used in solving the Kirk diffusion equations.
Open with DEXTER

We assumed the bulk velocity of the jet material as 0.9c with the shock moving at a velocity of 0.1c relative to the jet material. We first determined a linear baseline in flux space for the entire microvariability curve. This is interpreted as the background emission of the laminar flow in the jet. The pulse code written in IDL following KRM required the solution of a transcendental equation for the cooling length and a numerical integration for the total intensity emitted by each turbulent cell. The numerical solution of these equations were for the case where there is a constant injection rate Q₀, before switch-on, which occurs at t = 0. The solution involved the variables t_acc, t_esc, B, and γ which are physical parameters that determine the pulse rise and decay time; the amplitude is determined by changing the Q value for the pulse duration. The magnetic field B and the angle of the field relative to the line of sight θ can be different for each cell in a turbulent medium thus affecting the amplitude of the synchrotron emission and is folded into the Q parameter. Pulses are produced by allowing Q(t) = Q₀ for t < 0 and t > t_f, while the strength of the pulse is Q(t) = (1 + η_f)Q₀ for 0 < t_f where 1 + η_f represents the factor by which the rate of injection increases during t_f. The intensities are then calculated by I(ν,t) = I₁(ν,∞) + I(ν,t_flare) which is similar to Eq. (24) in KRM paper. The value of Q during the flare determines the amplitude of the pulse, the ratio t_acc/t_esc determines the shape, and the duration t_f determines the width of the pulse. After a number of test solutions, we determined that for t_acc/t_esc = 0.5 we get symmetric pulse shapes similar to what we see in the microvariability curves. The pulse shape is shown in Fig. 8. We used this as the standard pulse shape for all subsequent fits.

	Fig. 8 Single pulse emission in the local frame of reference due to the particle density and magnetic field enhancement.
Open with DEXTER

4.3. Modeling the 72-h flux curve

We fit thirty-five pulses to the microvariability curve by varying the width and amplitude of the standard pulse. Figure 9 shows the light curve fitted with the convolved pulses. The blue points are the smoothed data and the red points are the model. The resulting parameters for the pulses used in modeling the light curve are listed in Table 4. The first column of Table 4 shows the shot number in time order, while the second column records the center time of the pulse. The center time is related to the relative location of the cell along the jet as the shock progresses at a velocity of 0.1c. The amplitude and width (τ_flare) of each of pulse are given in the 3rd and 4th columns respectively, and the calculated electron density enhancement in the 5th column. The amplitude of the pulse over the background gives the value of Q, which is related to either the magnetic field enhancement or the density enhancement. The number quoted in Col. 5 is calculated assuming the pulse is due to the density enhancement only. Column 6 gives an estimate of the size for the cell in AU based on the assumed shock speed (u_s = 0.1c) and the duration of the pulse. The final column is an index number determined by sorting the cell sizes in order from the smallest to the largest size. The resulting fit of the 35 convolved pulses over the baseline flux of 7.85 mJy to the 72-h light curve had 2507 degrees of freedom (n_pts − (35 × 3) − 1) and the fit yielded a correlation coefficient of 0.98. Although the fit is not unique, it is representative of how well the model compares to the data.

	Fig. 9 Light curve fitted with the convolutions of the synchrotron pulses of various amplitudes and widths presented in Table 4 – the curve in blue is the data and the one in red is the fit. The correlation coefficient between them is 0.98.
Open with DEXTER

We can estimate the turbulent parameters from the observations by first noting that the very large duty cycle for microvariations in this object suggests the Reynolds number is well above the critical value and the plasma flow is normally turbulent. Since turbulence is a stochastic process, each micro-variability curve is a realization of that stochastic process. There is a large range of length/time scales for the turbulent vortices. The smallest vortex timescale is normally associated with the Kolmogorov scale (where most of the dissipation takes place in non-relativistic plasma) and the largest length scales are associated with either the size of the plasma jet or the correlation length within the plasma.

	Fig. 10 Distribution of the cell sizes in AU constructed using the data from the Table 4. The figure shows that the cell sizes form a continuum up to the size about 60 AU.
Open with DEXTER

Relativistic simulations show that turbulent relativistic extragalactic jets show a similar relationship between vortex length scale and energy (Zrake & MacFadyen 2012) as that seen in non-relativistic plasmas. Figure 10 shows the cell size in AU plotted against index number. The largest cell sizes are  ~160 AU which could either correspond to the correlation length scale or the physical width of the jet, while the smallest cell sizes seen are around 6 AU, could correspond to the Kolmogorov scale length of the turbulent plasma. The distribution of cell sizes is continuous which is what we would expect from a turbulent jet. Further analyses of the turbulent properties of the jet are underway by looking at other microvariability curves of this object in terms of this model. Fitting pulses to other microvariability curves might give us a much better indication of the turbulent nature of the plasma in these sources since each microvariability curve is a realization of the turbulent plasma in the jet.

5. Conclusions

The WEBT provided an excellent means to use longitude distributed ground-based telescopes to continuously monitor a very active blazar with minute time resolution. As a result of that program, we compiled a 72-h nearly continuous light curve of 0716 + 714. This light curve enabled us to do detailed time series analysis on the data and to confirm the maximum rise/decline timescales seen by other observers in the microvariability curves of this object. Based on these observations, we feel we may have resolved the fastest fluctuations. Time series analysis failed to yield reasonable periodicities and noise analysis did not confirm the presence of Gaussian distributed noise, so we have developed a model based on individual rapid synchrotron pulses caused by a shock wave encountering individual turbulent cells in the jet flow. The size of the vortex region and the magnetic filled orientation of the turbulent cells is stochastic, and results in pulses. We can decompose the microvariability curve into individual pulses and can get a picture of the underlying turbulent structure. We were able to fit the 72-h nearly continuous flux curve in terms of this model and obtain an estimate of the range of sizes of the turbulent cells and the density enhancements.

Acknowledgments

I would like to thank Florida International University for the financial support through Dissertation Year Fellowship 2012. The Abastumani Observatory team acknowledges financial support by the Georgian National Science Foundation through grant GNSF/ST09/521 4-320. Shao Ming Hu would like to thank the support by the National Natural Science Foundation of China under grants 11143012, 10778619 and 10778701. Rumen Bachev acknowledges the support by Science Research Fund of the Bulgarian Ministry of Education and Sciences through grants DO 20-85 and Bln 13-09.

References

All Tables

All Figures

	Fig. 1 Raw light curve of S5 0716+714 obtained by the compilation of some of the high quality data by major contributors. The light curves contributed from each observatory is plotted together with different symbols indentifying the observatory according to the codes given in Table 1.
Open with DEXTER
In the text

	Fig. 2 Flux curve of 0716+714 smoothed using the algorithm discussed in the text.
Open with DEXTER
In the text

	Fig. 3 Power spectrum of the smoothed Light curve. The horizontal line at low power indicates the maximum power level of simulated random noise light curves. The numbers associated with the peaks in the power spectrum correspond to the frequencies and periods listed in Table 3.
Open with DEXTER
In the text

	Fig. 4 Wavelet transform of the smoothed light curve linearized by removing a slope of  − 2.5 × 10-5 mJy/h and a Y-intercept of flux 0.01 mJy. Morlet kernel of 6th order was used from the wavelet application in IDL.
Open with DEXTER
In the text

	Fig. 5 rms variations versus average flux for the entire 3-day observation.
Open with DEXTER
In the text

	Fig. 6 log–log plot of the power spectrum of 0716+714.
Open with DEXTER
In the text

	Fig. 7 Example of geometry used in solving the Kirk diffusion equations.
Open with DEXTER
In the text

	Fig. 8 Single pulse emission in the local frame of reference due to the particle density and magnetic field enhancement.
Open with DEXTER
In the text

	Fig. 9 Light curve fitted with the convolutions of the synchrotron pulses of various amplitudes and widths presented in Table 4 – the curve in blue is the data and the one in red is the fit. The correlation coefficient between them is 0.98.
Open with DEXTER
In the text

	Fig. 10 Distribution of the cell sizes in AU constructed using the data from the Table 4. The figure shows that the cell sizes form a continuum up to the size about 60 AU.
Open with DEXTER
In the text