Gaussian Process Modeling Blazar Multiwavelength Variability: Indirectly Resolving Jet Structure

Blazar jet structure can be indirectly resolved by analyzing the multiwavelength variability. In this work, we analyze the long-term variability of blazars in radio, optical and X-ray energies with the Gaussian process (GP) method. The multiwavelength variability can be successfully characterized by the damped-random walk (DRW) model. The nonthermal optical characteristic timescales of 38 blazars are statistically consistent with the $\gamma$-ray characteristic timescales of 22 blazars. For three individuals (3C 273, PKS 1510-089, and BL Lac), the nonthermal optical, X-ray, and $\gamma$-ray characteristic timescales are also consistent within the measured 95$\%$ errors, but the radio timescale of 3C 273 is too large to be constrained by the decade-long light curve. The synchrotron and inverse-Compton emissions have the same power spectral density, suggesting that the long-term jet variability is irrelevant to the emission mechanism. In the plot of the rest-frame timescale versus black hole mass, the optical-$\gamma$-ray timescales of the jet variability occupy almost the same space with the timescales of accretion disk emission from normal quasars, which may imply that the long-term variabilities of the jet and accretion disk are driven by the same physical process. It is suggested that the nonthermal optical-X-ray and $\gamma$-ray emissions are produced in the same region, while the radio core which can be resolved by very-long-baseline interferometry locates at a far more distant region from the black hole. Our study suggests a new methodology for comparing thermal and nonthermal emissions, which is achieved by using the standard GP method.


INTRODUCTION
Flat spectrum radio quasars (FSRQs) and BL Lac objects (BL Lacs) are classed into a special class of active galactic nuclei (AGNs) called blazars, whose jets nearly point to the Earth. Blazars are highly variable over the entire electromagnetic bands. One popular scenario is that the accretion onto a supermassive black hole is the central engine, driving relativistic jet. But the detailed process is still unclear. Thanks to the high variability of blazars, one can investigate the physical process close to the central engine (e.g., Rieger 2019), such as the location of the emitting region and the jet-disk connection (e.g., Ackermann et al. 2016;Meyer et al. 2019;Zhang et al. 2022).
Using advanced interferometric instruments, blazar radio jet can be directly resovled on ∼parsec-scale (see Hovatta & Lindfors 2019, for a recent review). This provides a calibrator for multi-band variability analysis. There have been lots of works attempting to investigate the underlying physical process of blazar jet with multi-band variability (e.g., Chatterjee et al. 2012;Nakagawa & Mori 2013;Xiong et al. 2017;Goyal et al. 2018Goyal et al. , 2022. Max-Moerbeck et al. (2014) investigated the time-domain relationship between radio and γ-ray emission of blazars, and found the correlations only exist in a minority of the sources over a 4 yr period. They found radio variations lagging the γ-ray variations, suggesting that the γ-ray emission originates upstream of the radio emission. This result is further verified by Liodakis et al. (2018) who concluded that the radio variation is usually substantially delayed to the other wavelengths for blazars. Bhatta (2021) analyzed the correlation between optical (V -band) and γ-ray variabilities for blazars and found that the optical variability is highly correlated with the γ-ray variability except for 3C 273, however, no significant lagging is found. The multi-band variability analysis can be considered as an indirect approach for resolving blazar jet.
The GP method becomes popular in modern time-domain astronomy (e.g., Ryan et al. 2019;Burke et al. 2021;Yang et al. 2021;Griffiths et al. 2021;Covino et al. 2022;Rueda et al. 2022;Stone et al. 2022;Zhang et al. 2022). The GP method enables us to effectively extract information from astronomical variability. For example, Zhang et al. (2022) used a GP method to characterize the γ-ray variability of AGNs with stochastic process. It is found that the DRW model can successfully fit the γ-ray variability, which is similar with the optical variability of AGN accretion disk (Kelly et al. 2009;Li & Wang 2018;Burke et al. 2021). Moreover, Zhang et al. (2022) suggested a connection between the jet and the accretion disk by comparing the rest-frame γ-ray timescales of blazars with the optical accretion disk timescales of quasars.
In this work, we analyze the multi-band variability of blazars with the GP method, which is independent of the temporal correlation analysis. We hope to extract additional information from the variability. Using the data from Fermi-Large Area Telescope (Fermi-LAT), we carried out systematic research of γ-ray variability of AGNs recently (Zhang et al. 2022). So far, the Small and Moderate Aperture Research Telescope System (SMARTS) monitoring program 1 ) and the Steward Observatory (SO) spectropolarimetric monitoring project 2 (Smith et al. 2009) can provide almost ten years' (from 2008 to 2018) optical data of Fermi blazars. RXTE AGN Timing & Spectral Database 3 (Rivers et al. 2013) provides long-term X-ray data, and the Owens Valley Radio Observatory (OVRO) 40 m program (Richards et al. 2011) provides radio light curves (LCs) from 2008 to 2020 4 . Using these public data, we analyze the radio, optical and X-ray variability of three individual blazars, as well as optical variability for a sample including 38 Fermi blazars. The format of this paper is as follows. In Section 2, we describe the data as well as the GP method. The modeling results of the three individual sources and 38 blazars are shown in Section 3. We give discussions and physical interpretations of the results in Section 4. In Section 5, we conclude the paper with a brief summary.

Data and Sources
We use photometric data of blazars from the SMARTS and SO monitoring projects. The SMARTS program gives photometric data at five wavelength bands (B, V, R, J, K), which were taken from the 1.3 m telescope at the Cerro Tololo Inter-American Observatory. SO is a long-term optical program to support the Fermi Telescope, utilizing both the 2.3 m Bok Telescope on Kitt Peak and the 1.54 m Kuiper Telescope on Mt.Bigelow in Arizona. The campaign of the SO program spanned almost a decade from 2008 November to 2018 July. The X-ray data can be gained from RXTE observation which provided us with 16 yr (1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012) data in 2-10 keV. OVRO 40 m program gives radio data of blazars from 2008 to 2020, which is a large-scale, fast-cadence 15 GHz radio monitoring program. We select sources having long-term continuous observations and a good sampling. For the source with a large gap in the LC, we only use the data covering a longer period before or after the gap for analysis. Finally, We have 38 blazars in the optical band, including 23 FSRQs and 15 BL Lacs. Three blazars (3C 273, BL Lac, and PKS 1510-089) have long-term RXTE X-ray data. They are also in the sample of selected optical sources. Unfortunately, among the three sources, only 3C 273 has the OVRO LC. Table 1 gives the general information of these targets.

Gaussian Process Method
GPs are a class of statistical models, which are widely applied for modeling stochastic processes. For the one who is interested in the stochastic behavior of astronomical variability, GP provides a flexible method to model the LC with stochastic processes. The application of GPs for astronomical time-series is discussed in a recent review (Aigrain & Foreman-Mackey 2022). Considering a data set of y n at coordinates x n , the GP model consists of a mean function µ θ (x) parameterized by θ and a kernel function (covariance function) k α (x n , x m ) parameterized by parameters α (Foreman-Mackey et al. 2017). For time-series data, the GP is one-dimensional, and the coordinates are time, x n =t n . After obtaining the likelihood function with the above information, one can use Bayesian inference to estimate the posterior distribution over the parameter space.  In practical application, the key point is choosing the kernel function. The DRW process (called Ornstein-Uhlenbeck process in physics) is widely used to describe the variability of AGNs (e.g., Burke et al. 2021), and it is defined by an exponential covariance function (e.g., Kelly et al. 2009;Zu et al. 2013), where t nm = |t n − t m | is the time lag between measurements m and n. The amplitude term (σ DRW ) represents the amplitude of the random disturbance, and the damping (characteristic) timescale (τ DRW ) represents the timescale that the system returns to the stability after experiencing a disturbance. Sometimes, an excess white noise term (σ 2 n δ nm where σ n is the excess white noise amplitude and δ nm is the Kronecker δ function) is needed in the situation that there is a white noise in the LC in addition to the quoted measurement errors (Foreman-Mackey 2018; Burke et al. 2021). A more complex kernel is the stochastically driven damped simple harmonic oscillator (SHO), which is described by a second-order differential equation (Foreman-Mackey et al. 2017). The SHO kernel has been used to model the AGN accretion disk (Yu et al. 2022) and jet (Zhang et al. 2022) variability.
Celerite software package is a GP tool for a stationary process (Foreman-Mackey et al. 2017). It uses the semiseparated structure of a special covariance matrix to directly analyze and compute the GP likelihood for large data sets. Yang et al. (2021) and Zhang et al. (2021Zhang et al. ( , 2022 have tested the efficiency of this method for the study of AGN jet variability, and suggested that celerite has a strong potentiality for studying AGN variability (also see Burke et al. 2021). Here, we use the DRW model implemented in celerite package to model the multi-band variability of blazars.
The Markov Chain Monte Carlo (MCMC) sampler emcee 5 is adopted to perform posterior analysis. We assume log-uniform priors on each of the parameters. The MCMC sampler is run for 50,000 iterations with 32 parallel walkers. The first 20,000 steps are taken as burn-in. After modeling the LCs, we should estimate the fitting quality for assessing whether the fitting results are reliable, e.g., whether the standardized residuals follow a Gaussian white-noise sequence.
The power spectral density (PSD) can be constructed by using the fitting results. The DRW PSD is in the form of It is a broken power-law form with slope 0 below the broken frequency (f b ) and slope -2 above the broken frequency. The conversion between the timescale τ DRW and f b is τ DRW = 1/(2πf b ). The LC with large cadence or insufficient length leads to a large bias on the characteristic timescale derived from modeling. If the timescale is larger than the mean cadence of LC and less than 1/10 of the length of LC, the measurement of the damping timescale from the LC is reliable (Burke et al. 2021).

Results of 3C 273, PKS 1510-089 and BL Lac
We first analyze the multi-band variability of the three individual sources, 3C 273, PKS 1510-089, and BL Lac. We present the celerite fitting results of the LC for each source in the following. The measured timescales given in the main text are with errors in 95% confidence intervals.
For 3C 273, the optical data in both B and V bands are available. We show the modeling results in Figure 1, in which the left panel is for B-band LC and the right is for V -band LC. The DRW model can agree well with both LCs. Looking at the distribution of standardized residuals and the auto-correlation function (ACF) of standardized residuals (see details in Zhang et al. 2022), we believe the characteristic of each LC has been captured successfully. Through MCMC sampling, we get the posterior probability density distributions of two parameters (σ DRW and τ DRW ) and show them in Figure 2. The values are listed in Table 2. The results are different between the two bands. The parameters can be constrained by the B-band data but with large uncertainties, e.g., τ DRW = 59 +41 −28 days. Comparing the timescale with the cadence and the length of the LC, we believe that the B-band timescale is reliable. A broken frequency corresponding to the characteristic timescale is shown in the B-band PSD ( Figure 3). While the V -band timescale is ≈ 3200 days, very close to the length of the LC. This means that the V -band timescale is unreliable, which is also confirmed by the single power-law PSD ( Figure 3). We show the modeling results of X-ray LC, the posterior probability density distribution of parameters, and the PSD in the right panel in Figure   respectively. The values of the parameters can be found in Table 3. It is shown that the DRW model can describe the X-ray variability of 3C 273. The parameters are well constrained. The X-ray PSD presents a broken frequency that corresponds to a timescale of τ DRW = 28 +7 −6 days. We give the radio results together with the X-ray results. The radio PSD (the left panel of Figure 6) of 3C 273 is a single power law. The radio timescale is too large to be reliable.  For PKS 1510-089, the V and B-band LCs can be described by the DRW model ( Figure 7 and Figure 8). The V -band τ DRW of 39 +18 −14 days is larger than the B-band τ DRW of 11 ± 3 days ( Table 4). As expected, we get a smaller value of f b in V -band PSD (Figure 9). The X-ray LC of PKS 1510-089 also can be fitted well by the DRW model ( Figure 10). The parameters are well constrained (Table 3), and the PSD is in the form of typical DRW PSD. A trusted timescale τ DRW = 26 ± 3 days is obtained.
For BL Lac, only the V -band and X-ray data are available. For the X-ray data, there are two large gaps in the first 2800 days of LC, we then take the following 2500 days of LC for analysis. When modeling the two LCs of BL Lac, we get poor fitting (ACF of residuals deviating from the white noise) with the two-parameter DRW model. An excess white noise term is then added to the DRW model, and we model the LCs with the three-parameter DRW model again. The modeling of LC, the posterior distribution of parameters, and the broken power-law PSDs are shown in Figure   are shown in the right panels. One can see that the three-parameter DRW can fit the LCs well. Note that the highest flux point in the LC is poorly fitted. After removing the highest flux point, the modeling results are unchanged. We obtain the X-ray timescale of τ DRW = 63 +49 −30 days and V -band τ DRW of 47 +26 −19 days (Table 4). We applied the SHO model to the optical and X-ray data of BL Lac. The fitting is not improved significantly, and the parameters cannot be constrained. This suggests that the SHO model is not a good choice. The DRW including an additional white noise can describe the variability behavior. The value of σ 2 n (0.01 for the V -band LC; 0.04 for the X-ray LC) is larger than the squared of the mean size of the light curve uncertainties (σ y 2 ) where σ y 2 =0.0001 for V -band and 0.0036 for X-ray data. We have σ 2 DRW >σ 2 n +σ y 2 for both the V -band and X-ray data, which ensures that the fitted DRW amplitude is reasonable (Burke et al. 2021). It is possible that the quoted measurement errors are underestimated, and the excess white noise can account for excess measurement noise.
The γ-ray variability of the three sources has been analyzed in our previous work (Zhang et al. 2022) with the same method. We give the multi-band timescales with the errors in 95% confidence intervals of the three sources in Table 4. For 3C 273, the B-band, X-ray, and γ-ray timescales are consistent within the errors. The V -band and radio PSDs   are single power-law, having no corresponding characteristic timescales. For PKS 1510-089, the V -band, X-ray and γ-ray timescales are consistent within the errors but the B-band one has a smaller value. For BL Lac, the V -band, X-ray and γ-ray timescales are also consistent within the errors.

Optical Results of 38 Blazars
The DRW model can successfully fit the long-term optical LCs of the 38 blazars. Based on the criteria of selecting reliable measurements of the damping timescale, we get reliable optical timescales for the 38 blazars. The basic information of the 38 blazars and the modeling results are given in Table 1 and Table 2, respectively. Except for Note-(1) source name, (2) data source, S is for SMARTS and SO is for Steward Observatory blazar data archive,    3C 273 and PKS 1510-089 which are analyzed in Section 3.1, the timescales for different optical bands are consistent for the other 36 sources. This indicates that the optical emission of the 36 blazars has the same origin, i.e., the jet emission. In Table 2, we only list one optical band result for these sources. The timescale is between 10 days and 200 days. Some notes should be given on PKS 2052-47 and Ton 599. The fitting to the LC of PKS 2052-47 needs an additional white noise, and the relation σ 2 DRW (0.16) > σ 2 n (0.026) + σ y 2 (0.0016) still holds. Ton 599 has big gaps and few data in the first half of its V -band LC, and we select the second half of the LC to analyze.  Note-(1) source name, (2)(3)(4)(5) multi-band damping timescales in the observed frame. The uncertainties of the damping timescales represent 95% confidence intervals of the distribution of the parameter.

Origin of the optical emission from 3C 273 and PKS 1510-089
The optical emission of 3C 273 and PKS 1510-089 is complicated. Blue bump can be seen in their multi-band spectral energy distributions (SEDs; e.g., Abdo et al. 2010;Nalewajko et al. 2012;Castignani et al. 2017). SED modeling results showed that the accretion disk has a significant contribution to the optical emissions of 3C 273 and PKS 1510-089 (e.g., Nalewajko et al. 2012;Yan et al. 2012;Castignani et al. 2017). In addition, Zhang et al. (2019) found that a long-term variation trend in the optical continuum LC of 3C 273 does not appear in the emission-line variation. This suggests that the long-term variation trend is not contributed by the accretion disk, and it could originate from the jet. Li et al. (2020) quantitatively decoupled the optical emissions from the jet and accretion disk in 3C 273 and found that the jet emission accounts for 10%-40% of the total optical emission. Pandey et al. (2022) studied the correlation between V -band flux and polarization degree (PD) variations using SO observation during 2008-2018. They found a significant positive correlation only in two of the ten observing cycles. Note that the PD is quite small, and it changes from 0.04% to 1.58% during 2008-2018. The V -band single power-law PSD we obtained here is different from the typical PSD of the accretion disk (Suberlak et al. 2021;Burke et al. 2021) and jet variability (Zhang et al. 2022). The complicated mixture of the jet and accretion disk emissions at the V -band may result in the single power-law PSD. The mixed emission also results in the weak correlation between V -band and Fermi γ-ray variabilities reported by Bhatta (2021). We find no significant correlation between B-band variability and γ-ray variability for 3C 273 and PKS 1510-089. Looking at the location of the blue bump in SED (Roy et al. 2021), we suggest that the B-band emission of 3C 273 is dominated by the accretion disk photons.
For PKS 1510-089, the V and B-band timescales are clearly different, indicating different origins for the two bands' emissions. The V -band polarization of PSK 1510-089 is averagely greater than that of 3C 273, varying from 0.2% to 25.82% (Pandey et al. 2022). Among the ten observing cycles during 2008-2018, a significant positive correlation Note-(1) waveband, (2) the range of black hole mass in solar mass, (3) the mean damping timescale (redshift-corrected) with unit day. The uncertainties of timescales represent 1σ confidence intervals.
between V -band flux and PD variations is found in 5 cycles. Moreover, Castignani et al. (2017) found a good correlation between the long-term SO V -band and γ-ray LCs. These results suggest that the V -band emission is dominated by jet contribution. Also looking at the location of the blue bump in SED (Nalewajko et al. 2012), the B-band emission with a smaller timescale of 11 days is suggested as the accretion disk contribution.

Comparing Optical and γ-ray results
Long-term Fermi γ-ray LCs of 22 blazars have been analyzed by Zhang et al. (2022) with the same GP method. The optical timescale in this work is generally consistent with the γ-ray timescale (Figure 14). We examine the consistency of the timescales in the two energy-bands by using a statistical significance test (T-test). We get t-statistic=1.1 and p-value=0.28 (>0.05), which means that in statistic there is little difference between the two groups of timescales. The optical amplitude term σ DRW is less than one, and the γ-ray σ DRW can be greater than 10. This means that γ-ray variability can be more energetic than optical variability.
We separated the sources into two groups with M BH < 10 9 M and M BH > 10 9 M . The mean timescales (redshiftcorrected) in different ranges of black hole mass are listed in Table 5. It is found that the mean timescale of the sources in the mass range of 10 9 -10 10 M is smaller in both γ-ray and optical energies. However, we have a few sources with the mass of 10 9 -10 10 M , therefore this result may be tentative.
In Figure 15, we plot the relationship between the damping timescale in the rest frame (τ rest damping ) and the black hole mass of blazars along with the results of normal quasars from Burke et al. (2021). The timescales should be modified into the rest frame with the following formula: An average Doppler factor of δ D =10 is used here and the redshift z for each source is given in table 1. We show the optical, X-ray, and γ-ray results in the plot. It is found that the nonthermal optical τ rest damping of blazars and the thermal optical timescale of normal quasars occupy the same space in the plot of τ rest damping − M BH . The X-ray results for the three individual blazars are also in the same area as the optical results. The B-band timescale of 3C 273 is a typical value of accretion disk timescale. The B-band timescale of PKS 1510-089 is an outlier value among the accretion disk timescales. This value significantly deviates from the relation between damping timescale and black hole mass reported by Burke et al. (2021).   Figure 15. Plot of the rest-frame timescale versus black hole mass. The gray data, lines, and area represent the optical accretion disk results for normal quasars taken from Burke et al. (2021). Red data are γ-ray results for blazars taken from Zhang et al. (2022), and the purple and blue data respectively represent the optical and X-ray results for blazars obtained in this work.
It is difficult to directly resolve the inner jet structure of the blazar 6 . Especially, the location of the high-energy emission region is still a hot open question (e.g., Madejski & Sikora 2016;Böttcher 2019). Multi-band variability analysis provides an indirect approach to resolve the emission regions. The cross-correlation method is frequently used in multi-band variability analysis (e.g., Liodakis et al. 2018;Bhatta 2021).
GP method has been wildly used to characterize the AGN accretion disk variability (Kelly et al. 2009;Zhang et al. 2018;Lu et al. 2019;Burke et al. 2021). In blazar science, it becomes popular in recent several years (e.g., Goyal et al. 2018;Ryan et al. 2019;Covino et al. 2020;Tarnopolski et al. 2020;Yang et al. 2021;Zhang et al. 2022). In this work, we use the GP method to study the multi-band variability of the blazar. This provides results independent of the cross-correlation method.
The γ-ray variability of the blazar has been studied by Zhang et al. (2022) with the GP method. Here we focus on the X-ray and optical variability of the blazar. Multi-band emission from the blazar is dominated by the nonthermal jet contribution. Two special blazars are 3C 273 and PKS 1510-089. An optical-ultraviolet bump appears in their SED, which is associated with their thermal accretion disk emission (e.g., Nalewajko et al. 2012;Yan et al. 2012;Castignani et al. 2017).
We fit the long-term optical LCs from the database of SO and SMARTS with the DRW model. Finally, 38 blazars with a reliable characteristic timescale are selected. Except for 3C 273 and PKS 1510-089, the timescales in different optical colors agree with each other for the remaining 36 blazars. This indicates that the emissions in different optical colors of the 36 blazars have the same origin, i.e., the jet emission. Ruan et al. (2012) modeled the optical LCs covering from 2002 December through 2008 March of 51 blazars using the DRW model. They found that the observed damping timescale peaks at ∼80 days, and the intrinsic timescale τ rest damping peaks at ∼800 days 7 . The distribution of the optical timescale obtained in this work is flat (Figure 14), and the average optical τ rest damping is ∼400 days, which is smaller than the result of Ruan et al. (2012). All blazars in our sample are Fermi-detected γ-ray sources. While the sample studied by Ruan et al. (2012) would be dominated by the blazars of non-Fermi detection. Therefore, the results indicate that the optical timescale of the blazar of non-Fermi detection may be longer than that of the blazar of Fermi detection. Xiong et al. (2015) found that the two population blazars indeed have different physical properties, for example, the blazar of non-Fermi detection has a smaller Doppler factor (Paliya et al. 2017).
In the reverberation mapping studies of 3C 273 and PSK 1510-089, a nonechoed long-term trend is found in the optical continuum LC Li et al. 2020;Rakshit 2020). This reveals the mixed origin of their optical emission. New clues on the origin of the optical emission can be found in our results. The V and B-band timescales of PSK 1510-089 are different. Its long-term V -band variability is correlated with the γ-ray variability (Castignani et al. 2017), suggesting that the V -band emission is dominated by jet contribution. The long-term polarization variation (Pandey et al. 2022) also supports that the nonthermal component is dominated at V -band. The V -band emission of 3C 273 seems to be more complicated. The jet contribution to V -band emission may be strongly time-dependent and may vary in a large range. This complicated mixture of jet and accretion disk emission results in a single power-law PSD. For the two sources, no significant correlation is found between B-band and γ-ray variabilities in our analysis. The B-band emission is naturally considered as the accretion disk contribution. For 3C 273, the B-band timescale of ≈ 60 days is a typical value for the accretion disk emission of normal quasars. While the B-band timescale of ≈ 11 days of PKS 1510-089 is significantly smaller, and it deviates from the τ rest damping − M BH relation of Burke et al. (2021) (Figure 15). This short timescale may imply special properties of its accretion disk.
The nonthermal optical, X-ray and γ-ray variabilities all have the typical DRW PSD. Namely, the PSD of synchrotron emission is the same as that of inverse-Compton (IC) emission, consistent with the simulations with a time-dependent one-zone leptonic blazar emission model (Thiersen et al. 2022). In other words, the long-term jet variability is irrelevant to the underlying emission mechanism. Burke et al. (2021) suggested that the DRW damping timescale measured from the accretion disk variability of normal quasars could be associated with the thermal instability timescale expected in the AGN standard accretion disk theory. Zhang et al. (2022) measured the γ-ray DRW damping timescale of AGNs from the Fermi-LAT data, and found that the γ-ray timescales of 23 AGNs occupy almost the same space with the optical variability timescales of normal quasars in the plot of τ rest damping − M BH . In this work, we add the nonthermal optical timescale of blazars in this plot. The nonthermal optical timescale of blazars also locates at the same region with the thermal optical timescale of normal quasars in the plot (Figure 15). This implies that the jet variability is relevant to the accretion disk. The thermal instability in accretion disk may not only cause the accretion disk variability but also the jet multi-band variability. Statistically, the nonthermal optical τ rest damping of 38 blazars are consistent with the γ-ray τ rest damping of 22 blazars. Individually (3C 273, PKS 1510-089, and BL Lac), the damping timescales of the jet variability in optical, X-ray, and γ-ray energies are consistent within the measured errors. Our results indicate that multi-band jet emissions are produced in the same region. However, we still cannot know the distance from the emission region to the central black hole. The radio observation is helpful to constrain this distance (Max-Moerbeck et al. 2014). We modeled the OVRO radio LCs covering over ∼ten years, and we obtain a single power-law PSD. In this work, we only show the radio result for 3C 273 as an example. We also modeled the 30-yr radio LCs of 3C 279 and 3C 454.3 obtained from Aalto University Metsähovi Radio Observatory, and we still get an unconstrained timescale. The results indicate the radio timescale is very large and may be larger than 10 years. Through the very long baseline interferometry (VLBI) observation, one can determine the distance from the radio core to the central black hole. Comparing the optical/X-ray/γ-ray timescale and the radio timescale, we can infer that the optical/X-ray/γ-ray emission region is far upstream from the radio core.

SUMMARY
We analyze the blazar's radio, optical, and X-ray variabilities using the GP tool celerite. The DRW model can successfully fit the jet multi-band variabilities. The multi-band characteristic timescale is used to probe the structure of the emission region in the blazar jet. Our main results are as follows.
(i) The synchrotron and IC emissions have the same PSD, i.e., the typical DRW PSD. This indicates that the jet's long-term variability is irrelevant to the underlying emission processes. In the plot of τ rest damping −M BH , the jet timescales locate at almost the same space as the accretion disk timescales of normal quasars, implying that the jet and accretion disk variability is driven by the same physical process (Zhang et al. 2022).
(ii) The nonthermal optical, X-ray, and γ-ray variability has a consistent characteristic timescale. The radio characteristic timescale is very long which cannot be constrained by decades-long LC. The results indicate that the nonthermal optical-X-ray-γ-ray emission is produced in the same region, which is upstream and far from the radio core. This supports the basic hypothesis of the standard Synchrotron-Self-Compton jet model.
The GP method provides a flexible approach to understand the variability pattern of AGN in the framework of stochastic process. Adopting the standard GP tool (Foreman-Mackey et al. 2017), we build the link between accretion disk (thermal emission) and the jet (nonthermal emission), i.e., Figure 15. This is a new methodology for comparing thermal and nonthermal emissions, additional to the comparison between the thermal and nonthermal luminosities (e.g, Ghisellini et al. 2011;Sbarrato et al. 2012;Ghisellini et al. 2014).