2.1. Swift-XRT

We conducted an X-ray spectral analysis of the source S5 0716+714 within the 0.3–10 keV range, utilising all the data recorded from April 2005 to January 2023, totalling 361 observations. The X-ray spectra were obtained using the automated online tool ‘Swift-XRT data products generator’ (Evans et al. Reference Evans2009). This tool generates X-ray light curves, spectra, and images for any point source within the Swift-XRT field of view. It automatically selects source and background regions based on the count rate and correct the products for instrumental artefacts, such as photon pile-up and CCD bad columns. The spectra were binned into groups of 20 photons using grppha to improve $\chi^2$ statistics.

Table 1. Table showing best-fit parameters of power-law spectral fitting of Swift-XRT observations.

(Table is available in its entirety in the supplementary material.)

Table 2. Table showing best-fit parameters of selected Swift-XRT observations using log-parabola model.

(Table is available in its entirety in the supplementary material.)

All spectra were fitted using the power-law (PL) and log-parabola (LP) models within XSPEC. The neutral hydrogen column density was fixed at $N_H = 3.11\times 10^{20} \, {\textrm{cm}}^{-2}$ (Kalberla et al. Reference Kalberla, Burton, Hartmann, Arnal, Bajaja, Morras and Pöppel2005). We discarded all observations with degrees of freedom less than 4. We have also excluded all spectra which rendered a reduced chi-square ( $\chi_{red}^2$ ) value greater than 2. With these criteria, we obtained 302 spectra. F-tests were performed to determine the statistical significance of the log-parabola model compared to the power-law model for each individual observation. The best-fit parameters for the power-law model and F-test details are provided in Table 1. In Table 2, we present the best-fit parameters of the log-parabola model for the 35 spectra for which the F-statistic value (F-value) was $ \gt $ 5 and the null hypothesis probability (P-value) $ \lt $ 0.05 (these observations are not used in the correlation studies).

2.2. Swift-UVOT

We downloaded all Swift-UVOT (Ultra-Violet and Optical Telescope) observations of S5 0716+714 from the HEASARC archive. UVOT has three optical filters; u (3465 Å), b (4392 Å), and v (5468 Å), and three UV filters; uvw1 (2600 Å), uvm2 (2246 Å), and uvw2 (1928 Å). Standard data reduction proceduresFootnote ^a were followed to obtain the spectral files. Multiple images in each filter were summed over the extensions using the uvotimsum tool for each observations. A circular region with a 5 arcsec radius centred at the source was used to extract the source counts. A circular region with a radius of about 20 arcsec, free from source contamination, was used for background estimation. The spectral products for each filter were then generated using the uvot2pha tool.

To perform a power-law spectral fit in XSPEC, we considered only the observations for which images were available in at least four filters (249/361). The optical/UV data were corrected for Galactic reddening using the UVRED model, with the parameter $E_{B-V}$ fixed at 0.027 (Schlafly & Finkbeiner Reference Schlafly and Finkbeiner2011). We found that some observations required the addition of systematic error to achieve better fit statistics. We discarded those data which demands more than 5% of the systematic error to obtain $\chi_{red}^2$ less than 2. After these selection criteria, we got 235 UVOT observations, and the best-fit parameters of the spectral fit are provided in Table 3.

Table 3. Table showing best-fit parameters of power-law spectral fitting of Swift-UVOT observations.

(Table is available in its entirety in the supplementary material.)

2.3. NuSTAR

The source S5 0716+714 was observed by NuSTAR (Harrison et al. Reference Harrison2013) on 2015 January 24 (MJD 57 046.11) and on 2022 April 4 (MJD 59673.39). This epochs are marked with solid orange vertical lines in Fig. 1. The observations were downloaded through NASA’s HEASARC interface, and the data processing utilised the NuSTARDAS software package (Version 2.1.1) within the HEASoft environment (Version 6.29). For ObsID 90002003002 (MJD 57 046.11), the source spectrum was extracted from a circular region with 45 arcsec radius, while for ObsID 60701037002 (MJD 59 673.39), a radius of 20 arcsec was used. Background estimation was performed using circular regions with radii of 70 arcsec and 50 arcsec, respectively, for ObsIDs 90002003002 and 60701037002. The nuproduct (Version 0.3.3) tool was employed to derive the source and background spectra after running nupipeline (Version 0.4.9) on each observation. Subsequently, the FPMA and FPMB source spectra were individually binned using the grppha tool to ensure a minimum of 30 counts per bin. The final spectra were fitted with an absorbed log-parabola model and the unabsorbed fluxes in the 3–79 keV energy range were used for the broadband spectral study.

Figure 1. Multi-wavelength light curve of S5 0716+714. From top to bottom: panel 1 shows 3-day binned $\gamma$ -ray light curve using Fermi-LAT observations in the energy range 0.1–300 GeV. Second panel shows 2-day binned $\gamma$ -ray light curve simultaneous with Swift observations. X-ray light curve in 0.3–10 keV energy range using Swift-XRT observations is shown in panel 3 (solid circles in red and open circles in magenta colour corresponds to the observations with power-law and log-parabola models, respectively) and optical/UV light curve in the 2–7 eV energy band using Swift-UVOT observations is given in the fourth panel. The dashed blue lines represent the average flux in each energy bands. The vertical lines in orange colour represent the epochs where NuSTAR observations were available.

Table 4. Table showing best-fit parameters of spectral fitting using a power-law model (along with TS values) for selected $\gamma$ -ray observational data (in 2-day bins) simultaneous with Swift observations.

(Table is available in its entirety in the supplementary material.)

2.4 Fermi-LAT

LAT (Large Area Telescope) is the primary instrument on board Fermi $\gamma$ -ray space telescope. In this study, we have determined 2-day binned Fermi-LAT data in the energy range of 0.1–300 GeV to match the Swift observations. To perform data reduction, we employed Fermitools and followed the standard proceduresFootnote ^b. Photons within a circular region of interest (ROI) with a 10 degree radius, centred on the position of S5 0716+714 were chosen. Only photon-like events falling under the category evclass=128 and evtype= 3 were included. Additionally, a zenith angle cut off (less than 90 degree) was employed to eliminate background $\gamma$ -ray contamination originating from the Earth’s limb. We conducted a binned likelihood analysis method to fit the data across the entire time interval. We used the standard templates iso_P8R3_SOURCE_V3_v1.txt and gll_iem_v07.fits as the isotropic background model and Galactic diffuse-emission model, respectivelyFootnote ^c. We included all $\gamma$ -ray sources within a circular region of 20 degree radius from the central point in the fitting process, and their spectral characteristics were adopted from fourth Fermi Large Area Telescope (4FGL) catalogue. The parameters of all sources within the ROI were left as free variables, whereas those of sources outside the ROI were set to their catalogue values. To evaluate the significance of detecting each source within the ROI, we employed the test statistic (TS), defined as TS=2 $\log (\mathcal{L})$ , where $\mathcal{L}$ represents the likelihood parameter of the analysis. We froze the spectral parameters of all sources with TS $ \lt $ 25. The resulting output file was used as the input sky model for the unbinned likelihood analysis, which was then employed to determine the $\gamma$ -ray flux and spectral index for the 2-day time bins that were simultaneous with Swift observations. We opted for the PowerLaw2 function to model the $\gamma$ -ray spectrum of S5 0716+714 within each selected time bin, as it provided a good fit for these short time intervals. We started with the ‘DRMNFB’ optimiser, and sources with a TS value less than 9 were removed when they did not converge during the fitting process. Furthermore, for sources with TS values between 9 and 25, all parameters except the normalisation were held constant. The obtained fits were then optimised using the ‘Newminuit’. The $\gamma$ -ray spectrum of the source in the 4FGL catalogue is described by a log-parabola model. In order to determine whether there is a significant curvature in the $\gamma$ -ray spectrum within the selected time bins, we computed the test statistics for both power-law (TS $_{\textrm{pl}}$ ) and log-parabola (TS $_{\textrm{lp}}$ ) functions. The significance of spectral curvature (TS $_{\textrm{c}}$ ) is then calculated as TS $_{\textrm{c}}$ = TS $_{\textrm{lp}}$ - TS $_{\textrm{pl}}$ . The best-fit parameters along with the TS values are provided in Table 4. A significant spectral curvature will result in a large positive value for TS $_{\textrm{c}}$ . Further, we produced a 3-day binned light curve for the entire period (till MJD 59975), and the obtained fluxes and indices were then used for the flux and index distribution study (Section 3.4).

3.1. Flux-index correlation

The scatter plots of the integrated fluxes in optical/UV, X-ray, and $\gamma$ -ray energies with their corresponding photon indices $\alpha_{\textrm{O/UV}},\alpha_{\textrm{X}}\, \textrm{and} \, \alpha_{\gamma}$ , respectively, are shown in Fig. 3. We observe a significant negative correlation between the optical/UV flux and $\alpha_{\textrm{O/UV}}$ . A Spearman rank correlation study between these quantities yield a correlation coefficient, $r=-0.66$ with a null hypothesis probability, $p \lt 0.0001$ . This negative correlation suggests a harder when brighter behaviour. Whereas, a mild positive correlation is observed between the X-ray flux in the 0.3–10 keV range and $\alpha_{\textrm{X}}$ ( $r=0.38$ , $p \lt 0.0001$ ), which implies a softer when brighter behaviour of the source. These correlation results are consistent with the broadband SED of the blazar shifting towards blue end during flux enhancement as shown in Fig. 6. The spread of $\alpha_{\textrm{X}}$ around 2 suggests the presence of both steep synchrotron and hard inverse Compton components in the X-ray band. Also, the number of observations with $\alpha_{\textrm{X}}$ $ \gt $ 2 is high compared to those with $\alpha_{\textrm{X}}$ $ \lt $ 2, suggesting that the synchrotron spectral component often dominates in this energy range (Fig. 6). The observed softer when brighter trend may then indicate a shift of the spectrum to high energies during flux enhancement (bluer when brighter) (Giommi et al. Reference Giommi1999). Similar behaviour of the source in the X-ray region was also reported in earlier studies by measuring the hardness ratio (Giommi et al. Reference Giommi1999; Zhang Reference Zhang2010; Wierzcholska & Siejkowski Reference Wierzcholska and Siejkowski2015). In the $\gamma$ -ray region, we did not find any significant correlation between flux and $\alpha_{\gamma}$ ( $r=-0.001$ , $p=0.987$ ). However, Geng et al. (Reference Geng2020) reported a significant spectral hardening of the $\gamma$ -ray spectrum with an increase in flux during certain outbursts. The absence of such a correlation over a long period, regardless of the flux states, could also suggest that the source displays either softer or harder when brighter behaviour during different flaring events. Alternatively, the presence of a peak or a break in the $\gamma$ -ray energy region can destroy the possible correlation. The latter inference is further asserted by the distribution of the $\gamma$ -ray spectral index around 2 (bottom panel of Fig. 3).

Figure 3. Scatter plots of the fluxes and indices in optical/UV, X-ray and $\gamma$ -ray energies. The upper left panel is for optical/UV, the upper right is for X-ray (data points in orange colour filled circles correspond to the X-ray observations with log-parabola model), the lower left is for $\gamma$ -ray in the 0.1–300 GeV energy range for 2-day binned data simultaneous with Swift observations, and the lower right is for 3-day binned $\gamma$ -ray data for the entire epoch. The red and blue coloured lines represent the flux-index relation in each energy-band with the variation in Lorentz factor and magnetic field, respectively. The shaded regions represent the maximum attainable flux at a viewing angle of $\theta$ = 2 degrees in each energy band (see Section 4).

Figure 4. Scatter plots between fluxes in optical/UV, X-ray and $\gamma$ -ray energies. The upper panel is for optical/UV (2–7 eV) and X-ray (0.3–10 keV), the lower left is for optical/UV (2–7 eV) and $\gamma$ -ray (0.1–300 GeV), and the lower right is for X-ray (0.3–10 keV) and $\gamma$ -ray (0.1–300 GeV). Data points in violet colour filled circles correspond to the X-ray observations with the log-parabola model. The red and blue colour lines represent the flux-flux relation with the variation in Lorentz factor and magnetic field, respectively.

Figure 5. Scatter plots between the spectral indices. The upper panel is for optical/UV ( $\alpha_{\textrm{O/UV}}$ ) and X-ray ( $\alpha_{\textrm{X}}$ )(filled triangles in dark-grey colour and filled inverted triangles in dark-yellow colour represent the observations with low and high X-ray fluxes, respectively), the lower left is for optical/UV ( $\alpha_{\textrm{O/UV}}$ ) and $\gamma$ -ray ( $\alpha_{\gamma}$ ), and the lower right is for X-ray ( $\alpha_{\textrm{X}}$ ) and $\gamma$ -ray ( $\alpha_{\gamma}$ ). The red and blue coloured lines represent the index-index relation with variation in Lorentz factor and magnetic field, respectively.

3.2. Flux-Flux correlation

The scatter plots between the integrated fluxes of optical/UV, X-ray, and $\gamma$ -ray energies are shown in Fig. 4. A significant positive correlation is observed between the X-ray flux and optical/UV flux, with $r=0.74$ and $p \lt 0.0001$ , whereas a moderate correlation is observed between the X-ray and $\gamma$ -ray flux ( $r=0.42$ , $p \lt 0.0001$ ). A moderate correlation is also observed between the optical/UV and $\gamma$ -ray fluxes ( $r=0.49$ , $p \lt 0.0001$ ). The latter correlation was also reported before, though it was performed for short time periods (Larionov et al. Reference Larionov2013; Reference RaniRani et al. 2013b). This study identifies for the first time a strong correlation between optical/UV and X-ray flux over a long duration for the source S5 0716+714. These flux correlations are also consistent with the bluer when brighter behaviour of the source which can be visualised from Fig. 6. The correlations between the $\gamma$ -ray flux and the low-energy fluxes are often interpreted as supporting a leptonic origin for the $\gamma$ -rays, as the same electron population is responsible for both low-energy photons (through synchrotron processes) and high-energy photons (through inverse Compton processes). However, the observed correlations do not eliminate the possibility of hadronic models. Additionally, a correlation between hard X-ray and neutrino emission in blazars was reported recently (Plavin et al. Reference Plavin, Burenin, Kovalev, Lutovinov, Starobinsky, Troitsky and Zakharov2024). This, in turn, may suggest a possible correlation between low- and high-energy emissions. The observed correlations between the fluxes thus do not allow us to confirm the origin of $\gamma$ -ray emission mechanism.

3.3. Index-Index correlation

The scatter plots illustrating the relationships between the optical/UV, X-ray, and $\gamma$ -ray indices are presented in Fig. 5. In contrast to the positive correlation between X-ray and optical/UV fluxes, the X-ray index shows a significant negative correlation with the optical/UV index. The Spearman rank correlation analysis yielded a coefficient of $r=-0.69$ , with $p \lt 0.0001$ . On the other hand, a mild correlation is found between the optical/UV and the $\gamma$ -ray indices ( $r=0.32$ , $p \lt 0.0001$ ), while a mild negative correlation between the X-ray and $\gamma$ -ray indices ( ${r=-0.32}$ , $p=0.0001$ ). The observed anti-correlation between the optical/UV and X-ray indices may be associated with the harder when brighter behaviour in optical/UV and softer when brighter behaviour in the X-ray region. This is again consistent with the spectrum shifting towards the blue end during high flux states (Fig. 6).

To investigate whether the correlation between X-ray and optical/UV indices depends upon the flux state of the source, we repeated the study for those observations with X-ray flux above and below the average flux. Interestingly, we found that both these states support the negative correlation; however, it is more prominent for low-flux states, the correlation coefficient being $r=-0.72$ ( $p \lt 0.0001$ ) while that for the high flux states is $r=-0.42$ ( $p \lt 0.0001$ ).

The power-law fit to the optical/UV spectrum mostly results with index $\alpha_{\textrm{O/UV}}$ $\gtrsim$ 2. This indicates that the optical/UV emission mostly falls on the decaying part of the synchrotron spectral component. The X-ray spectral indices, on the other hand, have both cases with $\alpha_{\textrm{X}}$ $ \lt $ 2 and $\alpha_{\textrm{X}}$ $ \gt $ 2. This suggests that the X-ray energy region may fall either on the raising part of the Compton spectral component ( $\alpha_{\textrm{X}}$ $ \lt $ 2) or decaying part of the synchrotron spectral component ( $\alpha_{\textrm{X}}$ $ \gt $ 2). If we assume the underlying electron distribution responsible for the broadband emission as a broken power-law, then $\alpha_{\textrm{X}}$ $ \gt $ 2 maps the high energy end of the particle distribution, while $\alpha_{\textrm{X}}$ $ \lt $ 2 maps the low energy end of the particle distribution. Since the condition $\alpha_{\textrm{O/UV}}$ $\gtrsim$ 2 was true for majority of the observations, the optical region maps the particle distribution immediately after the break energy. In Fig. 5 (upper panel), we marked this region by vertical and horizontal dashed lines.

The correlation study between $\alpha_{\textrm{O/UV}}$ and $\alpha_{\textrm{X}}$ for those observation having $\alpha_{\textrm{X}}$ $ \lt $ 2, resulted in a moderate negative correlation with $r=-0.46$ ( $p=0.001$ ). This result is consistent with earlier study for the blazar Mkn 421 (Baheeja et al. Reference Baheeja, Sahayanathan, Rieger, Jagan and Ravikumar2022). Such a correlation strongly disagrees with the radiative loss origin of the broken power-law electron distribution (Kardashev Reference Kardashev1962; Rybicki & Lightman Reference Rybicki and Lightman1986). We also performed the $\alpha_{\textrm{O/UV}}$ – $\alpha_{\textrm{X}}$ correlation study for the observations with $\alpha_{\textrm{X}}$ $ \gt $ 2. Again we found a moderate negative correlation with $r=-0.49$ ( $p \lt 0.0001$ ). This observed correlation may possibly be associated with a shift in the spectrum towards blue or red side. For instance, shift towards the bluer end will harden the optical/UV index (optical/UV flux close to the peak) while softening the X-ray index (the X-ray flux falling at the synchrotron tail) as shown in Fig. 6.

Our correlation analysis suggests that the dominant spectral changes encountered during different flux states may be associated with the shift in the SED towards bluer/redder end. Interestingly, this also indicates the spectral variation in blazar due to the changes in the underlying electron distribution may not be substantial. We explore the possible scenario under which these correlation results can be inferred in Section 4. The correlation between different observed quantities are summarised in Table 5.

3.4. Distribution of fluxes and indices

A long term flux distribution study has the potential to understand the process responsible for the flux variations in blazars. Studies conducted on blazar light curves across different energy bands largely support a log-normal variability (Vaughan et al. Reference Vaughan, Edelson, Warwick and Uttley2003; Romoli et al. Reference Romoli, Chakraborty, Dorner, Taylor and Blank2018; Sinha et al. Reference Sinha, Khatoon, Misra, Sahayanathan, Mandal, Gogoi and Bhatt2018; Rieger Reference Rieger2019; Khatoon et al. Reference Khatoon, Shah, Misra and Gogoi2020). Such a variability can be associated with the moving perturbation in the accretion disc (multiplicative process) or a large ensemble of mini-jets buried with in the blazar jet (Narayan & Piran Reference Narayan and Piran2012). Alternatively, for a power-law spectrum, a log-normal variability in flux can be an outcome of normal variation in the photon index (Sinha et al. Reference Sinha, Khatoon, Misra, Sahayanathan, Mandal, Gogoi and Bhatt2018; Khatoon et al. Reference Khatoon, Shah, Misra and Gogoi2020). The flux distribution of certain blazar light curves also support a double log-normal behaviour (Sinha et al. Reference Sinha, Khatoon, Misra, Sahayanathan, Mandal, Gogoi and Bhatt2018; Khatoon et al. Reference Khatoon, Shah, Misra and Gogoi2020; Thekkoth et al. Reference Thekkoth, Sahayanathan, Shah, Paliya and Ravikumar2023).

To investigate the behaviour of the flux and spectral variability, we performed Anderson-Darling (AD) tests on the logarithm of fluxes and the spectral indices (Press et al. Reference Press, Teukolsky, Vetterling and Flannery1992). Under the AD test, the deviation of a given distribution from normality is prominent when the test statistics exceed a critical value. In Table 6, we present the details of the AD test for spectral indices and the logarithm of fluxes, with critical values calculated at a 5% significance level for the null hypothesis. From the AD test, we find that only the optical/UV and X-ray index distributions satisfy the normal behaviour, while the distributions of the logarithm of fluxes in optical/UV, X-ray, and $\gamma$ -ray energies, as well as the $\gamma$ -ray index are inconclusive.

We extended the study on the distribution of indices and the logarithm of fluxes by analysing the corresponding histograms. These histograms were first fitted using a single Gaussian probability density function (PDF) defined as

(1) \begin{align} {f(x)} = \frac{a}{\sqrt{2\pi \sigma^2}} \,\,\textrm{exp}\left[- \frac{(x-\mu)^2}{2\sigma^2}\right]\end{align}

where, a is the normalisation factor, $\mu$ and $\sigma$ are the mean and standard deviation of the distribution, respectively. In case of poor fit statistics, the fitting is repeated with a double Gaussian PDF given by

(2) \begin{align} g(x) =& \frac{a}{\sqrt{2\pi \sigma_1^2}}\,\, \textrm{exp}\left[-\frac{(x-\mu_1)^2}{2\sigma_1^2}\right] \nonumber \\ +& \frac{1-a}{\sqrt{2\pi \sigma_2^2}}\,\, \textrm{exp}\left[-\frac{(x-\mu_2)^2}{2\sigma_2^2} \right]\end{align}

where, $\mu_1$ and $\mu_2$ are the means of the distribution with standard deviations $\sigma_1$ and $\sigma_2$ , respectively. The histograms of the logarithm of fluxes and indices along with the best-fit PDF are shown in Fig. 7. The parameters obtained from the fitting are provided in Table 7. Contrary to the AD test, we find the histogram of the logarithm of optical/UV fluxes can be well-fitted with a single Gaussian PDF with reduced chi-square, $\chi_{red}^2$ = 1.14. A plausible reason for this can be the non-availability of high flux states that govern the tail of the distribution and which is evident from the histogram (top left panel of Fig. 7). Since the AD test is sensitive to the tails of the distribution, the test results are inconclusive. The distributions of optical/UV and X-ray spectral indices, on the other hand, support a single Gaussian PDF with $\chi_{red}^2$ values 0.82 and 0.90, respectively. This is consistent with the conclusion drawn from the AD test. Nevertheless, if the results obtained from the histogram analysis were correct then we could extend our inference that the log-normal variabitly of the optical/UV fluxes may be associated with the normal variations in the spectral index.

Figure 6. The shift of the SED towards the bluer end as the source becomes brighter is demonstrated. The hardening of the optical spectrum and the steepening of the X-ray spectrum are associated with the blue shift.

A single Gaussian PDF did not provide a satisfactory fit to the histogram of logarithm of X-ray fluxes ( $\chi_{red}^2$ = 1.99). Instead, the double Gaussian PDF was able to provide a better fit ( $\chi_{red}^2$ = 1.47). The double Gaussian behaviour of the logarithm of X-ray fluxes may be associated with the presence of two emission processes (synchrotron and SSC) active at this energy band. However, such a feature does not appear in the index distribution. Extending this study to hard X-rays (with predominant emission from SSC process) can probably provide a better answer to this unusual behaviour. At the $\gamma$ -ray energies, the histograms of both indices and the logarithm of fluxes deviate from a Gaussian nature but can closely represent a double Gaussian PDF. The best fit reduced chi-squares obtained for a double Gaussian fit to indices and fluxes being 1.83 and 1.42, respectively. These fit results are better than the ones obtained for a single Gaussian fit (5.33 and 2.82, respectively). This is consistent with the conclusion drawn from the AD test as well.

3.5. Broadband spectral energy distribution

To examine the broadband emission mechanisms responsible for the optical/UV–X-ray– $\gamma$ -ray emissions from S5 0716+714, we fitted the observed SED using synchrotron, SSC and EC processes. We considered only those epochs for which NuSTAR observations were available. Among the selected epochs, X-ray and gamma-ray flux during MJD 57045–57048 is large compared to the epoch MJD 59640–59699. To model the broadband SED, we considered the emission region to be a spherical blob of radius R and embedded with a tangled magnetic field B. The emission region is assumed to be populated homogeneously by a broken power-law electron distribution of the form

(3) \begin{align} n(\gamma)\,d\gamma = K\times\left\{ \begin{array}{ll} \left(\frac{\gamma}{\gamma_b}\right)^{-p}\,d\gamma&\textrm{for}\quad \mbox {$\gamma_\textrm{min} \lt \gamma \lt \gamma_b$} \\ \left(\frac{\gamma}{\gamma_b}\right)^{-q}\,d\gamma&\textrm{for}\quad \mbox {$\gamma_b \lt \gamma \lt \gamma_\textrm{ max}$} \end{array}\quad {{\textrm{cm}}}^{-3} \right.\end{align}

where, p and q are the low and high energy indices, respectively, of the broken power-law electron distribution, $\gamma_b$ the break energy with K the corresponding number density, and $\gamma_\textrm{min}$ and $\gamma_\textrm{max}$ are the available minimum and maximum electron energies, respectively. The emission region moves down the blazar jet at relativistic speed with a bulk Lorentz factor $\Gamma$ at a viewing angle $\theta$ . The electron distribution undergoes energy losses through synchrotron, SSC, and EC processes. The source of external radiation field considered is the thermal IR photons from the dusty torus (EC/IR) at temperature $T\sim 1000$ K (Ghisellini & Tavecchio Reference Ghisellini and Tavecchio2009; Bła˙zejowski et al. 2000).

Table 5. Spearman correlations coefficient ( $r_s$ ) and null hypothesis probability (p) of various relations among fluxes and indices.

Table 6. Anderson Darling test results of index/flux distributions in optical/UV, X-ray and $\gamma$ -ray bands.

The observed synchrotron flux, after accounting for the relativistic and cosmological effects, can be estimated from the single particle emissivity ( $f_s$ ), following Rybicki & Lightman (Reference Rybicki and Lightman1986) as

(4) \begin{align} F_{syn}(\nu)=\frac{4\pi\delta^3(1+z)}{3d_L^2} R^3 \frac{\sqrt{3} e^3 B}{16m_ec^2} \int_{\gamma_\textrm{min}}^{\gamma_\textrm{max}} f_s(\nu/ \kappa\gamma^2)n(\gamma)d\gamma,\end{align}

where, $\delta\, ([\Gamma(1-\beta_\Gamma\, \textrm{cos} \theta)]^{-1})$ is the jet Doppler factor, z the redshift of the source and $d_L$ the luminosity distance. The quantity $\kappa$ is given by

\begin{align}\kappa = \frac{ 3 e \delta B }{16 m_e c (1+z)} \nonumber\end{align}

The observed SSC flux can be estimated (Blumenthal & Gould Reference Blumenthal and Gould1970; Sahayanathan, Sinha, & Misra Reference Sahayanathan, Sinha and Misra2018) as

(5) \begin{align} \begin{split} F_{ssc}(\nu) =\frac{\delta^3(1+z)}{d_L^2} \pi R^3 \sigma_T \nu_s & \int_{\gamma_\textrm{min}}^{\gamma_\textrm{max}} \frac{1}{\gamma^2} \int_{x_1}^{x_2} \frac{I_{syn}(\nu_i)}{\nu_i^2} \\& f_c(\nu_i, \nu, \gamma) d\nu_i n(\gamma)d\gamma \end{split} \end{align}

where, $\nu_s$ is the frequency of the scattered photon, $I_{syn}(\nu_i)$ is the synchrotron intensity at photon frequency $\nu_i$ , and the limits $x_1$ and $x_2$ decide the maximum and minimum frequency, respectively, of the synchrotron photons involved in the scattering process (Blumenthal & Gould Reference Blumenthal and Gould1970; Jones Reference Jones1968)

\begin{align} x_1&= MAX \left[\nu_{syn}^\textrm{min}, \frac{ \nu_s}{4\gamma^2(1-h\nu_s/\gamma m_ec^2)}\right] \quad \textrm{and} \nonumber \\ x_2&= MIN \left[\nu_{syn}^\textrm{max},\frac{\nu_s}{(1-h \nu_s/\gamma m_e c^2)}\right] \nonumber\end{align}

The emissivity function $f_c$ given by

\begin{equation}f_c(\nu_i, \nu, \gamma)= 2r\log r+ (1+2r)(1-r)+\frac{\zeta^2r^2(1-r)}{2(1+\zeta r)} \nonumber\end{equation}

with

\begin{align} r=\frac{\nu_s}{4\nu_i\gamma^2(1- h \nu_s/\gamma m_ec^2)} \nonumber\end{align}

and

(6) \begin{align}\zeta=\frac{4 h \nu_i\gamma}{ m_e c^2}\end{align}

For $\Gamma\gg1$ , the observed EC flux can be obtained as (Dermer & Schlickeiser Reference Dermer and Schlickeiser1993; Dermer & Menon Reference Dermer and Menon2009; Finke Reference Finke2016)

(7) \begin{equation} \begin{split} F_{ec}(\nu) \approx \frac{\delta^3(1+z)}{8d_L^2} R^3c \beta_{\Gamma} \sigma_T \nu & \int_{0}^{\infty} d\nu_i^* \int_{\gamma_\textrm{min}}^{\gamma_\textrm{max}} \frac{U_{ph}^*}{\gamma^2 \,\nu_i^{*2}} \\ & \phi(\gamma,\nu_s, \nu_i') n(\gamma)d\gamma \end{split}\end{equation}

where, $U_{ ph}^*$ is target photon energy density at frequency $\nu_i^*$ measured in the rest frame of the AGN. The emissivity function $\phi$ is given by

(8) \begin{align} \phi(\gamma,\nu_s, \nu_i') =\left[y+\frac{1}{y}+\frac{\nu_s^2}{\gamma^2\nu_i'^2y^2}- \frac{2 \nu_s}{\gamma\nu_i'y}\right]\end{align}

with

\begin{align} y=1-\frac{ h\nu_s}{\gamma m_ec^2} \nonumber\end{align}

Figure 7. Histograms of logarithmic flux (left) and index (right) in optical/UV (top panel), X-ray (middle panel) and $\gamma$ -ray (bottom panel) spectra of S5 0716+714. The solid red and orange curves represent the best-fit single Gaussian function and double Gaussian function, respectively, and dotted curves represent the two components of double Gaussian function.

Table 7. Best fit parameter values of the probability density functions fitted to the logarithm of flux and spectral indices in optical/UV, X-ray, and $\gamma$ -ray.

The broadband spectral fit is performed by coupling the numerical routines used for calculating the observed fluxes due to synchrotron, SSC, and EC processes with the statistical spectral fitting code XSPEC. The main parameters governing the observed broadband SED of the source are K, $\gamma_\textrm{b}$ , p, q, B, R, $\Gamma$ , $\theta$ , $\gamma_\textrm{min}$ , $\gamma_\textrm{max}$ and the energy density of the external target photon field $U_{ph}^*$ . However, due to the limited information available in the narrow optical/UV, X-ray, and $\gamma$ -ray energy bands, all these parameters cannot be simultaneously constrained. Nevertheless, additional constraints are imposed in the form of near equipartition between the non-thermal particle and magnetic field energy densities. We defined the equipartition parameter $\eta$ as the ratio of magnetic field energy density ( $U_m$ ) to electron energy density ( $U_e$ ). The initial fit was performed for the low-flux state by letting the parameters p, q, $\gamma_b$ , B, $\Gamma$ , $\eta$ , and $U_{ph}^*$ as free, while the parameters $\gamma_\textrm{min}$ , $\gamma_\textrm{max}$ , $\theta$ , and R were fixed at values 10, $10^6$ , 1.3 deg (see Section 4), and 4 $\times$ 10¹⁶ cm, respectively. However, the confidence intervals on the fit parameters were obtained only for p, q, $\gamma_b$ , $\Gamma$ , and B while the rest of the parameters were frozen to their best fit values. We found that the broadband SED of S5 0716+714 can be fitted reasonably well by considering synchrotron, SSC, and EC/IR processes and the corresponding source parameters are provided in Table 8. The model fit to the observed fluxes along with the residuals are shown in upper panel of Fig. 8.

Table 8. Best fit values of the model parameters from broadband SED fitting.

We repeated the fitting for the high-flux state with the parameters corresponding to the low-flux state set as an initial guess. The fitting was performed only on the parameters p, q, $\gamma_b$ , B, $\Gamma$ , and $U_{ph}^*$ and the best-fit model parameters are provided in Table 8. The model fit to the observed fluxes, along with the residuals are shown in lower panel of Fig. 8. Interestingly, we find the flux enhancement is predominatly associated with an increase in $\Gamma$ which will blue-shift the spectrum. Hence this result is consistent with our correlation study; however, the absence of low energy information (infrared) and the plausible degeneracy among the source parameters does not let us assert this conclusion. Further, the increase in $\Gamma$ did not cause a significant shift in the peak frequencies of the synchrotron and EC spectral component due to minor reduction in B and $\gamma_b$ during the high flux states.