A small scale structure model of jet based on the observation of microvariability

It is assumed that the emission from blazar jet consists of the stationary emission and flaring events (Kirk et al., 1998; Potter, 2018; Bhatta, 2021). The phenomenon of optical variation can be attributed to various factors, such as alterations in the Doppler factor $\delta$ resulting from changes in visual angle, and variations in the radiation flux stemming from the evolution of particle distribution, among other contributing factors. The model considers that flares emanate from the radiation cooling of particles within the jet, subsequent to their acceleration.

The sketch of jet geometry and emission zones is presented in Figure 2. Based on the light travel time argument, the observed light curves at different timescales may originate from excitation sources at different scales. We assume that the STV/LTV originate within the kinematic blobs or large scale extensions in the jet, where small scale turbulence components may be present. Within jet-extended structures or large spheres, localized regions of perturbation known as turbulent cells may be present, characterized by parameters like particle number density or magnetic field that diverge from those of the surrounding environment. The emission emitted by these cells excited by the shock would correspond to the numerous flare events in the IDV. Overall, local inhomogeneous structures in these blobs will form “flare-in-flare” patterns, and the variability light curves can be simulated as the result of multiple flares superimposed on the LTV. The jets-in-a-jet model proposed by Giannios et al. (2009) based on high-energy fast variability is an alternative interpretation of the “flare-in-flare” patterns. However, the applicability and comparison of these models will require further work for elaboration.

Marscher et al. (1992) studied a turbulence model of the plasma in radio jets and interpreted radio flickering as a result of a shock encountering this turbulence. Bhatta et al. (2013) based on the previous works and explained the microvariability phenomenon with a model of shock passing through turbulent cells, which was first proposed by Webb et al. (2010). The formation of the turbulent cells is a stochastic process, thus the microvariability light curve does not show repeatable periodicity. Each independent turbulent cell may have a different density, size, and magnetic field direction. The model assumes that the microvariations are the result of a shock propagating down a turbulent jet and energizing individual turbulent cells that cool by emitting synchrotron radiation in the form of pulses.

We improved the model by adding the assumption that there may be no causal connection between flares, that is, adjacent flares in time may also come from distant emission regions. It means that there may be more than one excitation source in the jet, and the difference between their excitation moments is similar to the light-travel-time for the spatial distance between them. So, each flaring event at different times is independent of each other and originates from the region of the blob corresponding to the dominant radiation. By breaking down the IDV profile into separate flares, we can gain insight into the turbulence distribution and physical parameters within the jet. For this purpose, it is imperative to obtain a pulse model derived from the turbulent cells.

For the pulse profile in the microvariability light curves, we used the equation of Kirk et al. (1998), hereafter KRM. The KRM equation shows that when a plane shock encounters a cylindrical density enhancement region, particles with various magnetic field directions and particle density are accelerated in front of the shock. An increase in magnetic field strength or particle number density will give rise to a corresponding increment in particle injection rate. Once the shock has fully traversed the turbulent cell, the particles contained within it cease to undergo energy gain. Whenever the cooling efficacy of particles surpasses their acceleration capability, the radiative flux will gradually diminish to the degree of the underlying emission region. The observed microvariability curves are the result of the turbulent cells interacting with propagating shock waves in turn.

Here, we develop the time-domain multi-region radiation model in order to explain more observational characteristics. Based on the geometry proposed above, the acceleration process and radiation mechanism of particles in a single emission region, where the magnetic field is uniform, are analytically studied. The model is computed by assuming the observer lies in the direction of the normal to the shock front. Predecessors have made substantial research on the analysis of KRM equation. In this section, we perform a semi-analytical study on the effect of various parameters on the pulse profile, based on the analytical work of the KRM equation and following the results of previous research (Xu et al., 2019). We use the derived baseline model as a comparison to establish a theoretical foundation for examining the effects of different parameters.

Based on the model assumptions mentioned above, kinetic equations are solved for the time, space and energy dependence of the particle distribution functions in the acceleration region and the emission region respectively. In the acceleration region, particles are accelerated at a shock front and cool by synchrotron radiation in the emission region behind it. We assumed that the electrons accelerated by the shock obey the particle distribution function given by Ball & Kirk (1992):

with

where $N$ is the number density of the electrons in the energy space represented by Lorentz factor $\gamma$. In the equation, the parameters $t_{acc}^{-1}$, $t_{esc}^{-1}$, $B$, $\sigma_{T}$, $m_{e}$, $\mu_{0}$ and $c$ are acceleration rate, escape rate, magnetic field strength, Thomson’s scattering cross-section, mass of an electron, permeability of the free space and the speed of light respectively. Particles are assumed to be injected into the acceleration process with Lorentz factor $\gamma_{0}$ at a rate of $Q$ particles per second. After acceleration, the particles escape to the downstream plasma, where they are cooled by synchrotron radiation. The time-domain evolutionary pattern of the particle population distribution can be obtained by solving the semi-analytic numerical equations.

We consider the superposition of the total emission obtained by the shock front passing through these large-scale plasma blobs as the global emission. For the flaring event under the STV, LTV described above can be regarded as static emission, that is, the background laminar flow component. This paper focuses on the flare characteristics under the STV, while the detailed discussion of the unified model under the multi-timescale, such as “flare-in-flare”, will be elaborated in future work.

Possible causes of microvariations or flares are diverse, including the existence of a perturbative component in a localized region of the magnetic field or particle number density. The observed emission is considered to be the superposition of static emission and local flaring events, and an increase of the injection rate by a factor $1+\eta_{f}$ for a time $t_{f}$ is assumed. Here we consider the Eq.(24) of KRM as a superposition of these components. After the semi-analytical numerical calculation of the above model, the flare profiles at different frequencies can be obtained.

Figure 3(a) shows the baseline model used in this article, the parameters used in the model are listed in Table 1. The correlation coefficients between $I$-band and other bands are about $99.71^{+0.25}_{-0.37}\%$, which shows that the multi-band flare profiles obtained by this model under the same parameters have high similarity. These patterns show that there are obvious differences in the amplitude, duration and peak time of flare in different bands. The frequencies used in this study to calculate the flare profiles are the typical frequencies that are commonly used by telescopes for observations in optical bands. To obtain the theoretical flare profile for each observed frequency, we perform a Doppler transformation and convert the observed frequencies into the frequencies within the plasma rest frame.

From the theoretical analysis, it can be known that under the condition of the same acceleration region and emission region, the flare amplitude, duration and peak time at different frequencies are not exactly the same. Differences in flare amplitudes at different frequencies will lead to time-domain evolution patterns of spectral indices, which will be analyzed in Section 3.3. As shown in the right panel of Figure 3, different peak time of flare at different frequencies will lead to time delay of light curves at different frequencies, which is caused by different acceleration efficiency of particles with different energies. Although there is little difference in flare duration scales at different frequencies of optical bands, it can still be concluded theoretically that flare duration at lower frequencies are longer. This characteristic indicates that for low-frequency flares, their duration scales usually correspond to the range of large-scale emission, so they are limited by the physical properties of the large-scale emission, and are less affected by the local small-scale structure. In order to further analyze the influence of these acceleration and radiation parameters on flare, the influences of different parameters on multifrequency flare characteristics will be analyzed in details below.

Based on theoretical analysis and numerical simulation, the effects of kinematic parameters (e.g., $\gamma_{0}$, $\gamma_{max}$, and $t_{esc}/t_{acc}$) on flux level, profile and spectral hardness evolution of flare at multifrequencies were investigated. The influences of different parameters on flare observation characteristics at different frequencies, such as flux level, center time and duration, are analyzed below. On this basis, flare symmetry, peak amplitude ratio and time delay between different frequencies, time domain evolution of spectral index ($\alpha$; we assume that the spectral energy distribution (SED) follows $\nu F_{\nu}\propto\nu^{-\alpha}$) and other characteristics are further compared.

Based on theoretical analysis, the evolution of particle number density distribution, and ultimately, the profile of flare, are impacted by several factors including the acceleration rate ($t_{acc}^{-1}$), escape rate ($t_{esc}^{-1}$), initial ($\gamma_{0}$) and final energy ($\gamma_{max}$) of particles, flare timescale ($t_{f}$), and the time it takes for a shock to pass through the plasma blob ($t_{b}$). The results of the numerical simulations align with the aforementioned theoretical expectations. Examples of flare curves for different parameter sets are illustrated in Figure 3. The distribution of particles with different energies is not uniform during the acceleration process, so the flare profile responding to different frequencies is not exactly the same, which will lead to systematic differences in flare at different frequencies.

The simulation results reveal high correlation between flare profiles at different frequencies, which indicates that frequency is not a main parameter affecting the shape of the flare profile. However, there are differences in flare amplitude, peak time and duration at different frequencies, which we believe will lead to differences in observed characteristics at different frequencies. Different combinations of parameters will lead to the asymmetry of the ascending and descending trend of flare profile, the lag or advance of flare center time, and the sharpening or flattening of flare shape. The flare characteristics caused by these different sets of parameters can provide a theoretical foundation for the observed variability curves.

When the escape rate is much faster than the acceleration rate, i.e., $t_{esc}/t_{acc}\ll 1$, the radiation intensity will be too low to observe the substantial variability characteristics. In contrast, the particles do not have enough time to cool before leaving the source. When $t_{esc}$ and $t_{acc}$ are comparable, the simulation shows that the flare profile’s rising and declining times have a similar range of values, so the values used here are similar. Numerical simulations demonstrate that as the ratio of $t_{esc}/t_{acc}$ increases, the acceleration becomes more efficient compared to the escape, resulting in a shorter time to reach the peak. This indirectly affects the frequency at which the flare peaks at, and therefore the time delay between each frequency peak.

Simulations conducted under the baseline model for 1-hour flares in the observer frame demonstrate that the peak times for $I$-band flares occur ahead of the $R$, $V$, $B$, and $U$-band by $0.65$ min, $1.05$ min, $1.74$ min, and $1.81$ min, respectively. The time delay between bands is not an intrinsic constant value, as it is influenced by various other parameters. Figure 3 shows that low-frequency flares can lead high-frequency flares for some parameters, while for other parameters the opposite may occur. Therefore, we further analyze the influence of different parameters on the time delay between different bands through simulation. Here, some quantitative simulation results are presented based on the comparison under the baseline model.

Figure 4 presents sample simulation results of the effects of parameters such as $\gamma_{0}$, $\gamma_{max}$, $t_{b}$, $t_{f}$ and $t_{esc}/t_{acc}$ on the time delays or amplitude ratios of flares in different bands. In these examples, all parameters align with the baseline model (Figure 3a) except for those utilized for quantile testing. The numerical simulation results show that the time delay of flares in different optical bands depends on the various sets of parameters. This implies that flares originating from diverse physical environments may manifest as either low-frequency or high-frequency lag modes, contingent upon the time-containing evolution of the particle distribution. If a variability curve contains multiple flares with similar origins, their superposition will produce a curve with high similarity between bands. In contrast, if the flares have distinct origins, the inter-band time delays, amplitude ratios, and other characteristics possess non-systematic differences, leading to weak correlation between the bands. Hence, this dissimilarity serves to differentiate the source of the flare.

The simulation results indicate that the amplitude ratio of the flare between the specific bands is dependent on the parameter sets. Illustrated in the right panel of Figure 4, the diverse parameters have varying impacts on the amplitude ratio, with the most pronounced influence seen from the $t_{esc}/t_{acc}$, $\gamma_{0}$ and $\gamma_{max}$. This suggests that the initial and maximum energies, as well as the acceleration and escape rates, are critical factors that impact the distribution and evolution of particles, leading to changes in the energy spectrum index over time. However, parameters such as $\delta$, $t_{f}$ and $t_{b}$ have little effect on the amplitude ratio between the optical bands. In addition, consider that the peak frequency of synchrotron radiation grows with increasing $\gamma_{max}$. Thus, if the synchrotron peak frequency is near the optical observation frequency, observed high-frequency flares may grow more than low-frequency flares, but it could also be the opposite. Essentially, the nonlinear evolution of the particle number distribution under different parameter sets results in a dissimilar growth rate of the flare fluxes across each band. As a consequence, there is a non-linear development of the spectral indices in the SED. The time-dependent evolution of the spectral indices will help to distinguish the origins of the variability, which will be analyzed and discussed in detail in the next subsection.

We have used the synchrotron radiation model (Kirk et al., 1998) to simulate time-evolving optical spectral variation patterns in a specific physical environment. As described in Section 3.1, we consider that the jet radiation component originates from the superposition of various components such as the global radiation and the local flare radiation. The global radiation is assumed to be constant at STV, and the flare model parameters here are based on the baseline model parameters (Figure 3a). We utilize the aforementioned simulation as the fundamental framework for analyzing how varying parameter sets impact the time-contained evolution of the SED spectral index and the intensity hardness pattern. The parameter set in the baseline model was configured to achieve $\alpha\sim 0.9$, based on calculations derived from observation data. Several hypotheses attempt to explain the steep SED in the optical regime, including its proximity to the exponential cut-off region, specific acceleration processes (Summerlin & Baring, 2012), and the fine-tuned superposition of components, etc. Nevertheless, further studies are needed to establish the exact origin of the steep spectrum.

From the Eq.(24) of KRM, it is shown that the ratio of the flux of the flare to the flux of the static emission ($\eta_{f}$; $\eta_{f}$ is defined in relation to the I-band flux, noting that the ratio of the flux increment to the basal flux differs among various bands.) significantly affects the fluctuation level of the amplitude in the active state, which in turn affects the specific value of the spectral index of the spectral energy distribution. Figures 5(a) and 5(b) illustrate the time-contained evolution of the $\alpha$ and $\alpha$-intensity maps for different values of $\eta_{f}$, respectively. Figure 5(a) shows that $\alpha$ can decrease or increase during the flare event, and the flare’s intensity only affects the range of the $\alpha$ variation, not its overall trend over time. In Figure 5(b), the $\alpha$-intensity loop pattern suggests multiple transitions between RWB and BWB behavior during the flare event. Similarly, the comparable loop patterns of other $\eta_{f}$ parameters indicate that changes in only flare intensity do not affect the hardness-intensity trend. This points to radiation intensity changes not being the primary cause of the shift between the BWB and RWB phenomena. It is worth noting that this does not mean that parameters such as $\eta_{f}$ do not affect RWB and BWB phenomena, and they control the significance level of BWB and RWB variation during the flare period: if $\eta_{f}\ll 1$, i.e., the static emission is significantly larger than the flare flux, then $\alpha$ is almost constant at this time; if $\eta_{f}\sim 1$, then significant BWB and/or RWB phenomena can be observed.

Based on the parametric analysis in Section 3.3, the evolution of the particle number distribution and multi-band flare amplitude ratio are significantly impacted by the ratio of $t_{esc}/t_{acc}$, which makes it a critical factor in the time-containing evolution with $\alpha$. As illustrated in Figure 5(c), the softening trend in the evolution of $\alpha$ becomes more apparent as the ratio of $t_{esc}/t_{acc}$ decreases, and vice versa. Figure 5(d) shows that the $t_{esc}/t_{acc}$ ratio’s variation affects not only the $\alpha$ profile over time but also the $\alpha$-intensity trend. For a specific flare event, the evolution of $\alpha$ may demonstrate a monotonic trend of either BWB or RWB.

In essence, we consider that the evolution of the spectral index over time at STV is due to the non-uniformity and non-simultaneity of the evolution of the particle number between different frequencies. Assuming that the ratio of flare amplitudes between the variability curves at different frequencies remains constant, it can be inferred that the $\alpha$ does not change with time. Thus, the observed variability of the $\alpha$ implies that the ratio of the intensity between the flares at different frequencies is changing with time. Furthermore, assuming similar flare profiles at different frequencies, if the acceleration and cooling processes of particles between different frequencies occur simultaneously, i.e., there is zero time delay between the flares at different frequencies, the pattern of the $\alpha$ evolution with time in a single flaring event will only show a monotonic RWB/BWB trend, while the presence of time delay will make it possible for both BWB and RWB to co-exist in the same flaring event. In summary, the amplitude ratios and time delays between multi-band flares can be used as measurements to distinguish differences in flares and track particle distribution evolution in accelerated and radiated regions.

From the theoretical analysis in Section 3.3, it is shown that the main parameters affecting the time delay and amplitude ratio are $t_{esc}/t_{acc}$, $\gamma_{0}$, and $\gamma_{max}$. Taking Figure 5(e) and 5(f) as examples, we show the variability patterns of the $\alpha$ under different $\gamma_{0}$. The $\alpha$ changes significantly in these cases, potentially due to the dominant presence of the hard component (flare) over the steady component. Besides, according to the Eq.(6) of KRM and $t_{acc}=(\beta_{s}\gamma_{max})^{-1}$, it is known that the variability of parameters such as $B$ as well as $\gamma_{max}$ will lead to the change of $t_{acc}$, which in turn will affect the change of $\alpha$. In addition to the above parameters, we also simulated the effects of other parameters in the Eq.(1) of KRM, and the results are consistent with the above theoretical analysis. To refine the theoretical model, we analyze it in the next section by fitting the observed data with statistical methods.

Before fitting, to avoid the influence of introducing STV/LTV into the IDV fitting, we obtained the curve of basal components by fitting the LTV trend with polynomials. By fitting multiple pulses in numerous microvariability curves, we can obtain the size and the enhanced particle number density of the small-scale inhomogeneous emission region, and the distribution of acceleration parameters corresponding to these flares (assuming that other parameters such as shock velocity and magnetic field strength remain constant locally). These non-uniform regions exist in the form of superimposed background laminar flow, and the decomposition of background radiation components and radiation enhancing components will be further discussed in the following results analysis. The fitting process is described below.

Based on previous successful studies (e.g. Bhatta et al., 2013; Xu et al., 2019; Webb et al., 2021) that fitted independent pulse profile to microvariability curves, we have further refined the fitting method. By fitting the time-domain evolution of the optical multi-wavelength variability curves and spectral indices simultaneously, we obtain their flare parameters and thus calculate the physical parameters of the emission regions, unlike the previous paper where only the baseline model profile was used as a template for the fit. Here, simultaneous multi-band optical data after re-sampling of the original data are used as the sample, excluding the effect of small fluctuations on very short time scales. The lowest value of the intraday flux is taken as the baseline flux, and the average of the enhanced flux above it is taken as the average daily flare growth. Note that the parameters of emission zone obtained only from the observed optical variability features as a distinction are concurrent and thus their solutions are not unique, but for flares with different parameter sets still indicate that their origins may not be the same. Unlike the method described in the previous article, the number of flares by estimating the number of inflection points and extrema of the light curves at several different frequencies as well as the evolution of the $\alpha$, and the initial value of the flare simulation is thus obtained.

Based on the above improvements, the automatic fitting method automatically fits the number of intraday flares, the amplitude and duration of individual flares by iterative random number and least squares methods to approximate the results that minimize the $\chi^{2}$ sum of the original data and the fitted curve. In this research, the light curves with at least 50 data points in each band on a single day were selected as the original data input, and the correlation coefficients between the light curves generated by the fitted flare superposition and the original data were all above 99%. A full version of the fitted light curve plots is available online. The duration, amplitude, and central time of each pulse were directly obtained during the fitting of the flares. The size of the turbulence zone, the particle number density increment, and the relative position of the flares are also calculated.

Some examples of the fitting for $I$, $R$, and $V$-band variability for BL Lacertae are presented in Figure 6. The top and bottom panels give the observed data and fit patterns, the evolution of the spectral indices, and the hardness-intensity diagrams, respectively. The spectral index $\alpha_{opt}$ is derived by fitting optical multi-band observations using a power-law function ($\nu F_{\nu}=A\nu^{-\alpha_{opt}}$). In each panel, the red dots are raw observations, and the blue square dots are the data points of multi-band simultaneity after smooth interpolation. The colored dashed lines indicate the different flares obtained by the fit, the solid black line is the total fit, and the horizontally dashed line indicates the constant background flux for a single day. The lower plot of each fitted pattern indicates the distribution of residual between the fitted results and the original data. The left panel of Fig. 6 shows the fitting results for the single/dual emission region obtained by the general fitting method, while the right panel shows the fitting results for the multi-region obtained with our improved method.

The correlation coefficients between observations and calculated curves for the $I$, $R$, and $V$-band obtained from the single/dual emission region are $99.24\%$, $99.21\%$, and $98.44\%$, respectively, and the correlation coefficient between the fitted result and the interpolated data of the hardness-intensity diagram is $74.22\%$. For comparison, the correlation coefficients of the light curves obtained by fitting the improved multi-region radiation model were $99.86\%$, $99.78\%$, and $99.72\%$, respectively, and the correlation coefficient of the hardness-intensity diagram was $95.58\%$. It is obvious that the flares in this IDV curve is more likely to originate from multiple non-uniform components than from similar emission regions, and the improved fitting method can explain the observed features more extensively.

The various flare durations, amplitudes, and center times obtained from the simulated IDV curves in Figure 6 are shown in Table 2. The full version of table 2 listing the flare parameters obtained from the BL Lacertae fitting is available online. Columns 1 and 2 of the table list the index of the flares and the observed bands, as well as columns 3 and 4 give the flare centers and amplitudes. The amplitude of the flare is then converted to the enhanced electron number density as listed in column 5. Column 6 gives the width of the flare in the observer’s coordinate system ($\tau_{flare}^{obs}$). Column 7 indicates the size of the local turbulent cells, which correspond to the radiation timescale of the flare according to the relation

where $\beta_{s}$c, $\Gamma$, $\beta$, and $\theta$ represent shock velocity, bulk Lorentz factor, normalized speed of the emission region, and angle with the line of sight of the jet, respectively.

After fitting and decomposition of the multi-band light curves for 89 nights selected for the BL Lacertae from $2009$ to $2021$, we obtained $378$ multi-wavelength flares, including $368$ in the $I$-band, $371$ in the $R$-band, and $369$ in the $V$-band. Using our fitting method and procedure, we directly obtained the pulse amplitude, duration, center time, and we calculated the enhanced particle number density, scale and relative position of the flare corresponding to the emission region. Since the multi-wavelength microvariability curves are obtained by simultaneous fitting, we can also infer the acceleration parameters based on the differences in the time delay and amplitude ratio of the flares at different wavelengths, and thus distinguish the composition of different emission regions.

Statistical distributions of the parameters that can be obtained directly from the fitting procedure are shown in Figure 7, where the left panel shows the distribution of the parameters for $378$ flares, with the colored triangular dots indicating the fitting results for different bands, and the right panel indicates the probability distributions obtained statistically for these fitted parameters. As can be seen from the statistical results, the distributions of the parameters of flare amplitude, duration, particle injection rate, and emission region scale conform a log-normal distribution, whose best fitting distribution is shown as a dashed line, with the expectation and variance indicated in the panel. Although the amplitudes and durations of the flares obtained by fitting curves in different bands are different, the particle injection rates and the scales of the radiation cells calculated from the parameters obtained in these different bands are very similar, which is consistent with the theoretical assumption that the radiation sources at different frequencies are in the same region.

The size distributions of the radiation cells calculated from the fitting results for the $I$, $R$, and $V$-band are $18.30^{+7.58}_{-5.25}$ AU, $18.03^{+7.63}_{-5.25}$ AU, and $17.77^{+7.44}_{-5.14}$ AU, i.e., statistically $99.73\%$ of the turbulent cell sizes are from $2.40$ AU to $33.66$ AU. Note that the cell size obtained here is sensitive to the assumed parameters, including shock velocity, flow speed, and inclination. The energy spectrum fitting of turbulent cells obeys a power-law distribution $E(k)\propto k^{-p}$ obtained from the statistics of their turbulent cell distribution, where $p=1.54^{+0.22}_{-0.21}$ and $k$ represents the wave number. This statistic is similar to the distribution obtained in relativistic turbulence simulations (Zrake & MacFadyen, 2012), as well as with the $E(k)\propto k^{-5/3}$ distribution derived from Kolmogorov theory. Based on the theoretical parameters from previous studies (e.g. Baring et al., 1997; Summerlin & Baring, 2012; Baring et al., 2017) and the statistical analysis of turbulence scales, the Reynolds number is roughly estimated to be $R_{e}\sim 10^{3}-10^{5}$. In addition, we find from the statistics that the distribution of particle injection rates may be a superposition of two Gaussian distributions, which predicts that the 378 flares we fitted may be located mainly in at least two large-scale plasma blobs at different locations during the observation time (about 10 years), a conjecture that will be analyzed and verified in further LTV studies.

The time delay of each flare at different frequencies are calculated based on the difference in the center time of the flares obtained by fitting the curves of different bands. The time delay distribution for 378 flares is presented at the bottom panel of Figure 7. It shows that the time delay distribution obtained by fitting approximately follows the normal distribution. The time delay distributions between $I$ and $R$-band, $I$ and $V$-band, and $R$ and $V$-band are $-0.84\pm 8.40$ min, $-0.53\pm 10.97$ min, and $1.22\pm 6.95$ min, respectively. The obtained minute-scale time delays from the fitting align with the theoretical predictions in Section 3.3. This suggests that the minute-scale time delays between IDV flares may originate from accelerating or radiating processes. Furthermore, the parameter-dependent time delays can cause fluctuations around the zero-delay, where co-spatiality plays a role. For comparison, we also calculated the discrete correlation function (DCF; Edelson & Krolik, 1988). The Figure 8 gives an example of the DCF obtained from the single-night IDV curve, which shows a clear positive peak. In order to estimate the most likely value of time delay, we fitted a Gaussian function over a limited range of time delays (red dashed line). Statistical time delay distributions were obtained for all 89 nights of independent DCF results between the $I$ and $R$-band, and between the $I$ and $V$-band, producing $\bar{\tau}=-0.33\pm 5.81$ min and $\bar{\tau}=-1.19\pm 5.88$ min, respectively. The results from the DCF calculations indicate that the time delays converge to the order of minutes, consistent with the time delays obtained by fitting flares.

From the parametric analysis, it can be seen that the ratio of the amplitudes of the flares at different frequencies depends on the physical parameters of the acceleration zone, and considering that the variation of the spectral index and its evolutionary trend originates from the amplitude ratios of the flares at different frequencies, we speculate that the difference between the spectral index and the amplitude ratios can reflect the similarities and differences of the acceleration zone in which the flares are located. According to the model described in Section 3, we decompose the variability into long-term components and short-term flares. Based on the fitting results, we have performed statistics on the observed flux ratios and the flux ratios of the different components obtained from the decomposition, as shown in Figure 9. The statistical results show that the distribution of the amplitude ratios of each component can be approximated as Gaussian distribution. The distribution of the flare flux ratio, the flux ratio of the intraday enhancement, the flux ratio of the basal component, and the observed flux ratio between the $R$ and $I$-band are $0.77\pm 0.21$, $0.78\pm 0.11$, $0.71\pm 0.02$, and $0.71\pm 0.02$, respectively. From the statistical results, the distribution of the intraday enhanced component is similar to that of the decomposed flare component, while the distribution of the flux ratio of the basal component is similar to that of the total flux ratio, which is consistent with the hypothesis of the physical model about the multi-component decomposition.

Based on the flux ratios of various components, we further calculate the distribution of the spectral indices, as shown in the bottom panel of Figure 9. The distribution of the spectral indices demonstrates that the basal component remains consistent with the overall distribution of about $1.84\pm 0.16$, while the distribution of the spectral indices obtained from the intraday enhanced fluxes is in a broad range of $0.83\pm 0.96$. These spectral index distributions illustrate that the dominant trend of the spectral index variation can be determined by the LTV, and that the observed evolutionary trends in the various diversities of intraday spectral indices may be due to the superposition of several different regions with their different acceleration environments. Furthermore, we find that the spectral index distribution obtained from the intraday increments can be approximately decomposed into two Gaussian distributions, which may provide support for a conjecture that the dominant optical emission region may be located in at least two large-scale blobs during the period 2009 to 2021.

According to theoretical analysis from Section 3.3 and the parametric simulations from Section 3.4, it is shown that different spectral index evolution patterns are obtained for different acceleration parameters, which means different variable components may have different color behaviors. In order to study the evolution of the spectral indices of various components (e.g., constant and variable components) and their correlation with flux variations, we used the methods of Hagen-Thorn (1997) and Larionov et al. (2008) to plot the flux-flux and hardness-intensity diagrams (HID) of various components, as shown in Figure 10. The flux-flux diagram demonstrates that the basal component (i.e., the constant component) shows almost identical linear dependence to the slope of the observed total flux, while the slope of the flux-flux dependence of the variable component is slightly different compared to them. Meanwhile, according to the HID, it demonstrates an overall BWB trend, while the local spectral index evolution is not a single linear dependence, and the presence of a loop pattern similar to that in the simulation of Figure 5 can be clearly found in numerous HIDs obtained from IDV curves.

Therefore, based on these observed features combined with the theoretical simulation results, we interpret this correlated evolutionary feature of the flux and spectral index as the result of the superposition of multiple components. As the LTV component that accounts for most of the flux, it dominates the overall dependence of the HID and observed flux-flux diagram, while these STV produced by small-scale emission regions in the jet are superimposed on these long-term components and determine the short-term spectral index evolution of the observed objects according to their relative contribution to the total flux. The evolution of the spectral indices of the large scale blobs in the jet will produce similar evolution due to the change of the Doppler factor or other parameters. In contrast, the HID evolution of the small-scale radiation cells will produce various specific patterns depending on the acceleration parameters, which can explain the evolution characteristics of most observed flux-flux diagrams and HID under this model conjecture.

The overall time delay between different bands is obtained by calculating the DCF for the IDV curves, and the definition of the time delay obtained by this method is different from that of the time delay we obtained by fitting the flares. The set of flares obtained during the decomposition of the variability curve is not unique and only one possible solution is provided here, thus there is a selection effect in the method of obtaining time delays by fitting. In order to compare the similarities and differences of the time delays obtained by these two methods, a statistical analysis of the DCF with the time delays obtained by fitting the flares is performed, as shown in Figure 11. As can be seen from the left panel, the time delay of each individual flare is closer to its overall time delay and fluctuates around its value. The right panel shows that, the larger the difference between the overall time delay and the time delay of each flare obtained from model decomposition, the worse the correlation of its IDV tends to be between different bands.

The similarities and differences between the time delay of the flare and the IDV curves can be illustrated by this model, where the overall IDV fluctuations are generated by the superposition of radiation from many different small-scale emission regions. If these areas have similar physical environments, the time delays between the different bands will be similar, and thus the IDVs of different bands produced by the superposition will have a high correlation coefficient. Conversely, because the physical environment differs greatly, the different flares at different frequencies will have different time delays and amplitude ratios, and thus the correlation coefficients will tend to become worse. However, the correlation between IDVs of different bands is also affected by other reasons and the above inference is only taken as a possibility.

In order to investigate whether there is a correlation between the parameters obtained from the fitting results, we have statistically analyzed the parameters such as the turbulence size of the flare, the particle injection rate, the amplitude ratio of different bands and the time delay, whose distributions are shown in Figure 12. From the statistical results, it can be seen that the distribution of the time delay is independent of the scale of the emission regions and the number of enhanced particles. However, the distribution of flare amplitude ratios in different bands differs in regions with different particle number density enhancements, and the amplitude ratio tends to a fixed value for larger particle number enhancement. We find that this value is close to the inferred value obtained from the flux-flux diagram, while the amplitude ratios are spread over a wide range when the injection rate is small. This suggests that the regions with stronger particle number enhancement are closer to or in large scale blobs in the jet, while the small turbulent regions may be widely distributed at different locations in the jet.

In addition, we compared the enhancement of the particle number density and size of the cells obtained by the fitting. Although there is no strict dependence between them, statistically the particle injection rate tends to be larger as the cell scale increases. Our statistical results are in agreement with the study of Zrake & MacFadyen (2012), where it was shown that the magnetic cells become larger as the turbulence grows. The size, location distribution, and degree of turbulence of these non-uniform regions have been a popular issue in microvariability studies in recent years, and the precise parameters of these regions need to be further investigated and verified.