Variability studies and modeling of the blazar PKS 2155-304 in the light of a decade of multi-wavelength observations

TL;DR

This study analyzes a decade of multi-wavelength data on blazar PKS 2155-304, revealing log-normal variability, potential periodicity, and modeling its behavior with a synchrotron self-Compton approach.

Contribution

It provides a comprehensive variability analysis over ten years and introduces a successful time-dependent model to explain observed behaviors and potential periodicity.

Findings

Optical and X-ray variability follow a log-normal distribution.

A tentative 700-day periodicity is identified.

A one-zone synchrotron self-Compton model reproduces variability evolution.

Abstract

The variability of the high-frequency peaked BL Lac object PKS 2155-304 is studied using almost 10 years of optical, X-ray and gamma-rays data. Publicly available data have been gathered and analyzed with the aim to characterize the variability and to search for log-normality or periodic behavior. The optical and X-ray range follow a log-normal process; a hint for a periodicity of about 700 days is found in optical and in the high energy (100 MeV <E<300 GeV) range. A one zone, time-dependent, synchrotron self-Compton model is successfully used to reproduce the evolution with energy of the variability and the tentatively reported periodicity.

Tables

Table 1. Table 1: F var subscript 𝐹 var F_{\mathrm{var}} values for each energy range of the PKS 2155 − - 304  data set presented in this work.

Energy range	$F_{\mathrm{var}}$
0.2–10 TeVa	$0.657\pm 0.008$
0.1–1 GeVa	$0.36\pm 0.04$
1–300 GeVa	$0.43\pm 0.02$
0.3–2 keV	$0.591\pm 0.004$
2–4 keV	$0.716\pm 0.003$
4–7 keV	$0.796\pm 0.004$
J	$0.383\pm 0.005$
R	$0.369\pm 0.003$
V	$0.371\pm 0.004$
B	$0.378\pm 0.004$
aFrom H.E.S.S. Collaboration et al. (2017)
bFrom Chevalier et al. (2015)

Table 2. Table 2: Left: Values of χ 2 superscript 𝜒 2 \chi^{2} for the Gaussian fit of Φ Φ \Phi and log ⁡ ( Φ ) Φ \log(\Phi) , with values for the significance σ 𝜎 \sigma . Right: Values of the reduced χ 2 superscript 𝜒 2 \chi^{2} of the constant and linear fits of the scatter plots shown in Fig.  4 for each light curve. ρ 𝜌 \rho is the correlation factor. The corresponding data set is named in the first column.

	$\Phi$	$\log \left(\Phi \right)$		constant	linear increase
	${\chi }^2/$d.o.f.	${\chi }^2/$d.o.f.	$\sigma$	${\chi }^2/$d.o.f.	${\chi }^2/$d.o.f.	$\rho$
X-ray	$80.0/12$	$12.5/9$	$7.69$	$782/10$	$260/9$	$0.85\pm 0.02$
SMARTS (J)	$56.6/13$	$5.1/12$	$7.18$	$3077/13$	$865/12$	$0.81\pm 0.01$
SMARTS (R)	$29.9/13$	$8.7/11$	$4.22$	$22462/16$	$7858/15$	$0.93\pm 0.02$
SMARTS (V)	$65.1/12$	$9.1/11$	$7.48$	$3800/15$	$1746/14$	$0.76\pm 0.01$
SMARTS (B)	$30.2/13$	$15.4/12$	$3.85$	$3676/15$	$2234/14$	$0.72\pm 0.01$

Table 3. Table 3: SSC model parameters of the steady state of PKS 2155 − - 304.

Normalisation	$N_0$	$2.7\times {10}^{47}$ electrons
Index	$\alpha$	$2.3$
Cut off energy	$\log \left({\gamma }_{\mathrm{cut}}/1\mathrm{eV}\right)$	$5.3$
Magnetic field	$B$	$0.1$ G
Radius	$R$	$2\times {10}^{16}$ cm
Doppler factor	$\delta$	35

Equations

$$\frac{\mathrm{d}{N}_e(E,t)}{\mathrm{d}t}=\frac{\partial }{\partial E}\left[\left(\dot{{\gamma }_{\mathrm{S}}}+\dot{{\gamma }_{\mathrm{IC}}}\right)N_e(E,t)\right]+Q(E,t),$$

$$Q(E,t)=N_0(t)E^{-\alpha (t)}\exp \left(\frac{-E}{{\gamma }_{\mathrm{cut}}(t)}\right),$$

Full text

1 Introduction

Blazars are a subclass of active galactic nuclei with a relativistic jet pointed towards the Earth. The observed emission, from radio up to TeV energies, is dominated by their jet. Nevertheless, the precise composition of the jet as well as the acceleration and emission processes involved are not known. If the jet is dominated by leptons ($e^{\pm}$ pairs), leptonic models such as the synchrotron self-Compton (Band & Grindlay, 1985, SSC) or the external Compton (Dermer & Schlickeiser, 1993) are invoked to reproduce the electromagnetic emission. These models differ essentially in the target photon field for the inverse Compton emission. Alternatively, hadronic models have been considered, where $\gamma$-rays are emitted through photo-production of pions (e.g. Mannheim, 1993) or synchrotron emission of protons (e.g. Aharonian, 2000). Unfortunately, to disentangle between these models, fits of the time-averaged spectral energy distribution (SED) are insufficient.

However, one of the striking properties of blazars is their variability111Another diagnostic tool is provided by multimessenger studies, notably in high-energy neutrinos, which till now is however limited to the single spectacular case of the blazar TXS 0506+056 (Aartsen et al., 2018).. Their brightness can vary at time scales ranging from minutes to years, and this behavior has been observed at all wavelengths, from radio waves to X-rays and gamma rays. The two classes of models predict different variability patterns. Hence, long-term observations and statistical studies of the variability are key tools in the comprehension of these objects.

PKS 2155$-$304 (redshift $z=0.116$; Falomo et al., 1993) is a very well known blazar detected at TeV energies for the first time in 1999 (Chadwick et al., 1999) and later confirmed by Aharonian et al. (2005). It has been classified as high-frequency peaked (HBL) thanks to X-ray observations from the HEAO-1 satellite (Schwartz et al., 1979). At TeV energies, this object exhibits large flares on minutes timescale (Aharonian et al., 2007) but also variation on longer timescales (H.E.S.S. Collaboration et al., 2017). At lower energy, Fermi-LAT reported variability on monthly timescale (Acero et al., 2015) as well as much more rapid flares (Cutini, 2014, 2013). In X-ray, the source exhibits variability (see e.g. Gupta, 2015) and even intra-day variability was reported (Pandey et al., 2017).

PKS 2155$-$304 variability in optical and X-ray can be deeply studied thanks to systematic observations by SMARTS, RXTE, Swift/XRT and XMM-Newton. With the impressive dataset recorded by H.E.S.S. in the TeV range over 9 years (H.E.S.S. Collaboration et al., 2017) and the constant monitoring of the Fermi-LAT, such studies can be extended to the $\gamma$-ray band as well.

This article is structured as follows: Section 2 presents the multiwavelength dataset used. Section 3 details our analyses on the variability of PKS 2155$-$304. A time-dependent modeling is presented in Section 4. In Section 5, we summarize our results and conclude.

2 Observations and analysis

2.1 $\gamma$-rays datasets

In the very high energy (VHE, $E\gtrsim 100$ GeV) range, this work makes use of the H.E.S.S. data presented in H.E.S.S. Collaboration et al. (2017), reporting on 9 years of observations of PKS 2155$-$304. In the high energy (HE, 100 MeV $<$ E $\lesssim$ 300 GeV) range, the Fermi data presented in Chevalier et al. (2015), where the light curves have been computed in two energy bands, are also used.

The possibility to extend the Fermi (as well as X-ray, and optical) light curve to further data taking periods was considered since—contrarily to H.E.S.S. data for an article external to the collaboration as this one—the former are available. Nevertheless, we deemed very important the role of the H.E.S.S. instrument for the phenomenological interpretation, since probing a unique spectral range, and thus more helpful than the added value of a data taking extension limited to the lower energies.

2.2 X-ray datasets

PKS 2155$-$304 has been regularly observed by the X-ray observatories RXTE, Swift/XRT  and XMM-Newton. Preliminarily, we have gathered both monitoring and target of opportunity (ToO) observations of PKS 2155$-$304. However, since ToO observations bias the dataset towards high flux values, only monitoring data have been considered. Further, we applied a correction for the Galactic absorption with $n_{\rm H}=1.48\times 10^{20}$ cm*-2*, according to Kalberla et al. (2005).

RXTE data consist in publicly available222http://cass.ucsd.edu/rxteagn/ pre-analyzed light curves in four energy ranges: the full range from 2 to 10 keV and three subranges 2–4, 4–7 and 7–10 keV.

Swift/XRT  data in the energy bands 0.3–2, 2–4, 4–7, 7–10 and 2–10 keV have been analyzed using the package HEASOFT 6.16. The data were recalibrated using the last update of CALDB and reduced using the standard procedures xrtpipeline and xrtproducts. Caution has been taken to properly account for pile-up effects for corresponding affected exposures, and spectral fits were performed using Xspec 12.8.2 assuming a power-law spectrum.

XMM-Newton  public EPIC (European Photon Imaging Camera) data in the energy ranges 0.3–2, 2–4, 4–7, 7–10 and 2–10 keV have been reduced using the SAS software package (version 14.0) and analyzed following Tatischeff et al. (2012).

Figure 1 presents the X-ray light curves in the total 2-10 keV range as well as in the 4 energy sub-ranges.

2.3 Optical dataset

SMARTS (Small and Moderate Aperture Research Telescope System, Bonning et al., 2012) data are publicly available333http://www.astro.yale.edu/smarts/glast/home.php. Magnitudes have been corrected for the absorption of the Galactic foreground following Schlafly & Finkbeiner (2011) and converted to a spectral flux density using the zero flux values of Cohen et al. (1992). The light curves, shown in Fig. 2, are taken in four bands (J, R, V and B) in the same time windows as Fermi (MJD 54603 to 56622).

3 Results

3.1 Fractional variability

To study the variability of PKS 2155$-$304 , the fractional variability $F_{\mathrm{var}}$ as defined in Vaughan et al. (2003) has been computed for each energy bin (Table 1). Figure 3 presents the evolution of $F_{\mathrm{var}}$ as a function of the energy (hereafter variability energy distribution). The value in the TeV range computed in H.E.S.S. Collaboration et al. (2017) is also reported. There is a clear trend, with $F_{\mathrm{var}}$ increasing with energy from the optical range to X-ray. In $\gamma$-rays, $F_{\mathrm{var}}$ is lower than in X-ray but also increases with energy. This bimodal structure has been already reported for this source in Aharonian et al. (2009b) but is shared also by other objects, e.g. Mrk421 (Giebels et al., 2007; Ahnen et al., 2016).

3.2 Log-normality of the flux

We fit the X-ray and optical flux distribution $\Phi$ and its logarithm $\log(\Phi)$ with a Gaussian. The results are summarized in Table 2, also reporting the $\chi^{2}$ values. All light curves present a $>3\sigma$ preference for a log-normal distribution, i.e. the distribution of $\log(\Phi)$ is better described by a Gaussian than the distribution of $\Phi$.

The excess variance $\sigma_{\rm XS}$ as defined in Vaughan et al. (2003) vs. the average flux $\overline{\Phi}$ is shown in Figure 4. Each light curve is split in several bins with at least 20 points to ensure sufficient statistics to compute $\sigma_{\rm XS}$ and $\overline{\Phi}$. These two quantities are clearly correlated, and linear fit is found to better reproduce the data than a constant fit (see Table 2). Although insufficient statistics has prevented to reach firm conclusions in the HE range, H.E.S.S. Collaboration et al. (2017) has reported a similar behavior in the VHE range. Overall, these results are suggestive of multiplicative processes as main responsible for the variability of PKS 2155$-$304 through the whole spectrum. In such a stochastic process, a high flux leads to an increased variability which, in turn, possibly leads to higher flux.

Evidence for log-normality on different time scales has been reported for different sources: BL Lacertae in X-ray (Giebels & Degrange, 2009), in the VHE for the BL Lac Mrk5 01 (Tluczykont et al., 2010; Chakraborty et al., 2015) or for PKS 2155$-$304 , during flaring event in VHE (Aharonian et al., 2009a). This behavior was also observed for non-blazar objects, such as for the Seyfert 1 galaxy IRAS 13244$-$3809 in X-rays (Gaskell, 2004). There are growing evidences that this behavior is a common feature of accreting systems. In the context of galactic X-ray binaries, where log-normal flux variability has first been established, such a behavior is thought to be linked to the underlying accretion process (Uttley & McHardy, 2001). The detection of log-normality in a Seyfert 1 galaxy, a class of radio quiet AGN lacking a relativistic jet and whose emission line emission correlates with the amount of gas surrounding the central engine, strengthens the link between accretion disk and log-normal behavior.

3.3 Search for periodicity

Sandrinelli et al. (2014) reported a possible periodic behavior in the optical and HE light curves of PKS 2155$-$304 and an intriguing coincidence of a period in HE roughly double the one in optical was noted. To study further such periodic features, the multi-wavelength light curves are analyzed with a Lomb Scargle periodogram (LSP; Lomb, 1976; Scargle, 1982). The LSP is a method to estimate the Power Spectrum Density (PSD) of a time series based on a least squares fit of sinusoids to the data sample. The advantage of the LSP compared to a classical Fourier analysis is that it takes into account the uneven spacing of the data. The standard LSP was used (section 5 of VanderPlas, 2018), as implemented in the astropy package (Astropy Collaboration et al., 2018).

In this work, the light curves are not evenly sampled and binned differently. Moreover, gaps between observations - appearing from the impossibility to observe the source during some period of the year - have to be taken into account. For all the analyses, the maximum frequency is set following section 4.1.3. of VanderPlas (2018) to $f_{\mathrm{min}}=1/(2\delta t)$ where $\delta t$ is the typical integration time. To have uniform results, this has been set according to the most constraining data set, i.e. the Fermi-LAT integration time (10 days).

The LSPs of the SMARTS (in the R band only for the sake of clarity), X-ray, Fermi-LAT and H.E.S.S. light curves are shown in Fig. 5 along with the 1$\sigma$ and 2 $\sigma$ confidence intervals. The X-ray and H.E.S.S. light curves do not show any periodic feature. In optical, the B,R,V and J bands exhibit a significant periodicity, with the best fit period ranging from $715$ to $733$ days depending on the band, while the HE light curve is found to have a periodicity of 685$\pm 9$ days. Sandrinelli et al. (2014) found a similar period in the HE range as the one reported here, but a T$\approx$315 day period for the optical light curve. While we do confirm the presence of a peak in the LSP of optical data around 300 days, the most intriguing excess of the power is at $\simeq$700 days (see Fig. 5), since it is found both in optical and HE light curves.

In order to assess the significance of our results, light curves without periodicity have been simulated and rebinned according to the observational time binning. This allows one to factor out instrumental effects such as windowing due to sparse observation and/or binning in time due to limited sensitivity. Each simulated light curve has been obtained by inverse Fourier transform of power-law noise, without adding a constant term, adopting a different spectral index for each energy ranges. In VHE, an index of 1.1 has been taken while in HE, the used value is 1.2. (H.E.S.S. Collaboration et al., 2017). In optical and X-ray, an index of 1.2 has been used. We have also checked that the results are robust with respect to slight ($\sim 10-20\%$) changes of the spectral index, representative of typical fit uncertainties.Then, for each period, the range of LSP power spanned by 68% (95%) of the realizations is used to derive the 1 $\sigma$ (2 $\sigma$) contours shown in Fig. 5.

One limitation of the visual inspection of Fig. 5 is that it cannot obviously account for the trial factor effect, coming from a scan over different frequencies tested. Note that the astropy python package used here can provide an estimate of a false alarm probability (FAP) taking into account the trial factor; however, it implicitly assumes non-varying data with Gaussian noise, while the real data follow a red noise behavior. As such, this estimator cannot be taken at face value. We merely use it to perform some sanity checks, e.g. to verify that the FAP computation following Baluev (2008) yields more conservative results than the method (VanderPlas, 2018), as expected.

Despite this limitation, our results remain intriguing: While taken separately each of the peaks found in the right panels of Fig. 5 might not be statistically very significant, by interpreting e.g. the Fermi-LAT results as a test search to suggest the most interesting periods to search a priori in the optical bands, the $\sim 3\,\sigma$ excess found in SMARTS data sample at comparable period of $\simeq$700 days can be taken more or less at face value, since (most of) the trial factor is basically accounted for in the Fermi-LAT sample search. Albeit heuristic, this argument is also what singles out this period compared to others, for which no matching is seen in the multiwavelength comparison.

Besides the above-mentioned technical difficulties, these kinds of long-term periodicity analyses suffer from physical limitations, such as the fact that only a few periods are probed. It is also worth noting that in a recent article, Covino et al. (2018) warned that for none of the blazars whose periodicity in Fermi-LAT  band is reported in (Prokhorov & Moraghan, 2017; Sandrinelli et al., 2017; Zhang et al., 2017a, b) a strong statistical case can be made. However, Covino et al. (2018) only used Fermi data; one comforting cross-check from the multiwavelength data presented here is that cross-correlations between LAT and optical data are significant both at times $\tau=0$ (no delay) and at a timescale approximately equal to the reported periodicity. We take the simultaneous hints for a periodicity around 700 days in both Fermi and optical datasets as the most intriguing indication of our analysis. While awaiting a definitive confirmation in a priori searches in future independent datasets, in the next section we will tentatively considering the implications of including or not this periodicity for the interpretations in the context of a simple SSC model.

4 Time-dependent modeling

4.1 The synchrotron self-Compton model and the steady-state of PKS 2155-304

In order to reproduce the results found in this work and especially the variability evolution with energy, a one-zone synchrotron SSC model has been considered (Band & Grindlay, 1985). In this model, the first bump of the SED is produced by the synchrotron radiation of electrons spinning into the uniform magnetic field $B$ of the jet. The second bump of the SED is explained by inverse Compton scattering of the same electrons population on the previous synchrotron photon field. The emission is assumed to be produced by a homogeneous region of radius $R$ propagating in the jet with a Doppler factor $\delta$.

If electrons are injected with a time dependent function $Q(E,t)$ and radiate their energy via synchrotron or inverse Compton processes, the electron density $N_{e}$ is given by the diffusion equation:

[TABLE]

where $\dot{\gamma_{\rm S}}$ and $\dot{\gamma_{\rm IC}}$ are the synchrotron and inverse Compton cooling rates of the electrons. The escape of the electrons is not taken into account in this model. If the escape time scale is larger than the cooling time scale, this has no effect in the model. Lower time scales will lead to an achromatic decrease of the variability. The injection of the electrons $Q(E,t)$ is chosen to be a power-law with exponential cut-off:

[TABLE]

where $N_{0}$ is the injection normalization, $\alpha$ the power-law index and $\gamma_{\rm cut}$ the energy of the exponential cut-off.

The equation is solved numerically using the algorithm of Chang & Cooper (1970) for each time step, which allows us to follow the evolution of the electron density and hence of the emitted flux. The time-averaged SED has been modeled by reaching the steady state of the diffusion equation for $Q(E,t)\equiv Q(E)$. The parameters used are given in Table 3.

4.2 Simulation of the variability

To introduce the variability in the model used in this work, one of the parameters of the model was chosen to vary with time. Giebels et al. (2007) modeled the emission of Mrk 421 with a similar model and an injection function being a relativistic Maxwellian function. They found that two flux states (high flux and low flux) can be reproduced by merely changing the characteristic energy. Most interesting is that they predicted small flux variations in optical and GeV ranges and high variations in X-ray and TeV range. In the modeling presented here, $\gamma_{\rm cut}$ is similar to their Maxwellian characteristic energy. Since the number of injected particles increases exponentially with $\gamma_{\rm cut}$, the relevant parameter is then $\log(\gamma_{\rm cut})$.

AGN light curves generally show power-law noise of the form $1/f^{\beta}$ (Lawrence & Papadakis, 1993). To reproduce this property, $\log(\gamma_{\rm cut})$ is varied during the simulations following a power-law noise of index $\beta$ and total variance $\sigma$, and a mean of $\log(\gamma_{\rm cut})=5.3$. Simulations of the variation of $\log\gamma_{\rm cut}$ are drawn following Timmer & Koenig (1995). Note that time series of $\log\gamma_{\rm cut}$ were constructed on a timescale ten times longer than the needed amount, to ensure that long term variations are well reproduced using this technique.

Simulations with $\beta\in$ [1.0, 1.5, 2.0] and $\sigma\in$ [10%, 15%, 20%, 25%] were performed to find the couple that best reproduces the variability energy distribution. In total, 200 simulations of 10 years each with a binning of 9.5 minutes in the observer rest frame (to ensure that small variations are simulated) were performed. Values of $\beta=1.0$ (i.e. flicker noise) and $\sigma$= 20% were found to best match the variability energy distribution. The results of the simulations are found to be mostly sensitive to $\sigma$. Indeed, increasing this parameters increases the measured variability mainly in optical and Fermi energy ranges.

Figure 6 shows again the variability energy distribution (black points as from Table 1) but also reports the simulated $F_{\mathrm{var}}$ as an orange band. The variability increases from the lowest energy up to the X-ray domain, then drops in the HE range, and eventually increases towards TeV energies, following the same trend as the data, albeit quantitative discrepancies are noticeable.

To take into account the differences between each instrument, each simulated light curve is rebinned in time and energy following the observations of each instrument. The cyan points/boxes represent the variability of the simulated light curves after rebinning. Such binning does not change much the pattern in optical, X-ray and TeV energies. However, the variability in the Fermi-LAT ranges is sensibly reduced, now matching the data within errors. This is likely due to the important time bins (10 days) used for the Fermi-LAT analysis. The only band which appears in clear disagreement with the model is the optical range, where the variability of the SMARTS data is not reproduced; another source of variability has to be invoked.

4.3 Power spectral density

In H.E.S.S. Collaboration et al. (2017), both the HE and VHE power spectral density were computed. It was found that they can be quite well characterized by a flicker noise, the best-fit power-law indexes being $\beta_{\rm HE}=1.20^{+0.21}_{-0.23}$ and $\beta_{\rm VHE}=1.10^{+0.10}_{-0.13}$, respectively. For completeness, we report in Figure 7 a comparison of PSD obtained in the simulations (orange-shaded areas) vs. the multiwavelength data (VHE, HE, X-ray and optical/R band). The cyan band accounts for the correct binning in time, according to the actual observations: notice how this effect is crucial for the simulation to match the data. It is noticeable that this is an important consistency check, since nothing in the simulation was tuned to reproduce the PSD.

4.4 Adding a periodic component

The hinted periodicity of PKS 2155$-$304 in the optical band (and, to some extent, at HE) could be explained by different physical effects.

The most frequently discussed culprit for quasi-periodic behavior in blazars are quasi-periodic modifications of the Doppler effect. A fascinating possibility is that a binary supermassive black hole (SMBH) system could be at the center of these AGN, instead of just one SMBH as assumed in the general AGN picture (Begelman et al., 1980). This binary system could cause a periodic change in the accretion rate of the matter coming from the disk and even misalign the accretion disk (Doğan et al., 2015, and references therein). A similar outcome may be due to the Lense-Thirring effect, breaking the central regions of tilted accretion disks around spinning black holes, see e.g. (Doğan et al., 2015, , and references therein). These scenarios however face the difficulty that jet precession is expected to happen on too long time scales, $\sim 10^{6}$ years according to the analyses of King et al. (2008) and Nixon & King (2013). A recent study (Sandrinelli et al., 2018) also points out the tension that a binary SMBH origin associated to (the relatively common) blazar periodicity may have with pulsar timing array limits on the gravitational wave emission of such close SMBH binaries. In Raiteri & Villata (2017), optical-to-radio monitoring of the blazar CTA 102 has been argued to support a scenario where magnetohydrodynamic instabilities or the rotation of a twisted jet cause different jet regions to change their orientation444however see e.g. Zacharias et al. (2017) for an alternative interpretation of the recent variability exhibited in CTA 102., hence their relative Doppler factors. Other observational evidence in the AGN BL Lacertae and M 81 suggesting a precession motion of their jets looking at radio knots with VLBI observations can be found in Stirling et al. (2003); Caproni et al. (2013); Marti-Vidal et al. (2013).

With the aim to test if a periodic variation of the Doppler factor can account for the observations of PKS 2155$-$304 , we performed simulations analogous to what previously described, but adding on top of the stochastic variation of $\log(\gamma_{\rm cut})$ a sinusoidal time series $\delta(t)=\delta_{\textnormal{steady state}}+5\times\sin(t+T)$. We fix $\delta_{\textnormal{steady state}}=35$ from what is found in Section 4.1. The amplitude is an had-oc value.

This addition yields an increase of the global variability in an achromatic way. To compensate for this effect, the simulations were redone with $\sigma_{\rm cut}=15\%$, keeping $\beta=1$.

The resulting variability energy distribution is shown in Figure 8. The shape of the variability energy profile stays roughly the same, with however a flatter part in the optical range around $F_{\rm var}=0.20-0.25$, rising the variability levels of the simulation in this range compared to the non periodic one (Figure 6).

It is clear that the periodicity can (at least partially) explain why the SMARTS data are more variable than the previously considered model. Within this new scenario, this energy range would be dominated by the variability of the periodicity and not by the one of the stochastic process. This is not surprising, since a modification of the value of the cut-off energy has a small impact on the electrons producing the optical photons. It is also encouraging that a minor deviation between low-energy Fermi data and our simulations present in Figure 6 shrinks to an irrelevant difference in presence of periodicity.

Figure 9 displays the periodograms of the simulated light curves in the SMARTS, X-ray, Fermi-LAT and H.E.S.S. ranges, without applying any temporal binning on the short term simulated light curves. The optical and GeV simulated light curves have a clear and strong bump around 600-700 days, being appreciably (if not dominantly, for the optical) sensitive to the periodic variation of $\delta$, while the X-ray and TeV light curves show a less significant bump, since they are especially sensitive to the stochastic change of $\gamma_{\rm cut}(t)$. This is qualitatively consistent with the lack of detection of periodicity in X-ray or VHE ranges till now.

5 Discussion and conclusions

Multiwavelength data spanning up to 10 years of observations of the blazar PKS 2155$-$304  have been gathered and studied in this work. SMARTS optical, RXTE, Swift/XRT and XMM-Newton have been analyzed, as well as HE and VHE gamma-ray data taken from Fermi and H.E.S.S. Collaboration et al. (2017), respectively. Optical and HE ranges shows only little variability with respect to the variability found in X-ray or in VHE. This seems to be a characteristic shared with other BL Lac (e.g. Mrk 421) and indicate a close link in the population of particles that emits the low and high energy parts of the SED.

As in VHE (H.E.S.S. Collaboration et al., 2017), X-ray and optical bands exhibit a log-normal behavior, indicative of a multiplicative process. This is also an argument for a link in the origin of the variability of both bands. Nevertheless, the tentative indication for periodicity around 700 days found in optical and HE is observed neither in X-ray nor in VHE.

The time-dependent SSC model used in this work explains well the evolution of the variability with energy, except for the optical band. Adding a periodic component in the model helps in better describing this energy range, and also improves the agreement in the HE range. Although model-dependent, this is another encouraging indication that stimulates further studies to confirm the hint of a periodicity reported here with future, independent datasets. Our model also reproduces the non-detection of the periodic behavior in X-ray and VHE bands. Independently of how realistic the models discussed are, our results are a healthy reminder that depending on the energy range of interest, the mechanism(s) dominating the observed variability can be different.

Still, some questions remain on the origin of the variability of the whole spectrum of PKS 2155$-$304. For instance, it is worth keeping in mind that the observed log-normal behavior is not explained in the models discussed in this article, and a consensual quantitative theory of its microscopic origin is still lacking. Definitely, more long-term observations of different blazars would be needed to extract common features and differences between objects, in turn helping refining the theoretical models.

Acknowledgements

This work has been done thanks to the facilities offered by the Université Savoie Mont Blanc MUST computing center, as well as by CC-IN2P3 (cc.in2p3.fr).

This research has made use of NASA’s Astrophysics Data System, as well as data and/or software provided by the High Energy Astrophysics Science Archive Research Center (HEASARC), which is a service of the Astrophysics Science Division at NASA/GSFC and the High Energy Astrophysics Division of the Smithsonian Astrophysical Observatory.