3.1. Gamma-ray ( $\gamma-ray$ ) Light curve

As shown in the light curve (Fig. 1), three major flares are observed in Ton 599 from 2019 to 2024. We have studied the temporal evolution of all the flares separately. To illustrate the temporal evolution, we fitted the peaks with a sum of exponential functions, providing the rise and decay times for each peak visible in the light curve plots. The functional form of the sum of exponentials is given by Abdo et al. (Reference Abdo2010a)

(1) \begin{equation}F(t) = 2F_0 \left[ \exp\left( \frac{t_0 - t}{T_r} \right) + \exp\left( \frac{t - t_0}{T_d} \right) \right]^{-1}\end{equation}

where $F_0$ is the flux at time $t_0$ representing the approximate flare amplitude, and $T_r$ and $T_d$ are the rise and decay times of the flare.

Fig. 3 shows the light curve of flare 1 and flare 2 with 1-day binning corresponding to the flaring activity during MJD 59370-59420 and 59440-59500. The flare 1 shows six peaks $P_1$ , $P_2$ , $P_3$ , $P_4$ , $P_5$ and $P_6$ and maximum flux during the flare 1 was 2.18 $\times$ 10 $^{-6}$ ph cm $^{-2}$ s $^{-1}$ , which is the flux of peak $P_5$ . In flare 2, no Fermi-LAT data is available in the time range MJD 59394-59412, so the sum of exponentials does not fit it. So, we have tried to fit it using the rise and decay function of the sum of exponentials. The maximum flux during flare 2 was 3.06 $\times$ 10 $^{-6}$ ph cm $^{-2}$ s $^{-1}$ but it is not fitted finely so that showing the highest flux 5.13 $\pm$ 0.16 $\times$ 10 $^{-6}$ ph cm $^{-2}$ s $^{-1}$ . The light curve of flare 3 with 6-h binning corresponds to its flaring activity during MJD 59920-60030.

The result after fitting is shown in Table 3, where rise time $T_r$ and decay time $T_d$ for each flare are mentioned. Flare 1 shows six bright peaks, and their rise and decay time varies from 1.06 to 2.35 h. Peak $P_2$ , $P_3$ , $P_4$ , and $P_6$ is symmetric in nature with almost equal rise and decay time (within error bars) suggesting it may occur when a perturbation in the jet flow or a blob of denser plasma passes through a standing shock present in the jet (Blandford & Königl Reference Blandford and Königl1979). At peak, $P_1$ and $P_5$ rise is comparatively slower than the decay, suggesting a slow injection of electrons into the emission region. In flare 2, the peak rises very fast with a rise time of around 2.78 $\pm$ 0.12 h, and the flux decays very slowly, with a decay time of 7.83 $\pm$ 1.13 h. The fast rise and slow decay in Flare 2 suggest a longer cooling time for the electrons through the various processes. A similar trend is also seen in flare 3, where the decay is slower compared to the rise, suggesting slower cooling of the electrons. It has been argued that any physical process faster than the light travel time will not be detectable from the light curve, and hence, the rise and decay times will always be higher than the light crossing time, and all three scenarios of rise and decay time are possible.

3.2. Multi-wavelength light curves

In Fig. 2, we show the broadband light curves collected from various telescopes. In panel 1, we present a 1-day binned gamma-ray light curve where various flares and quiet periods have been marked in different patches of colours. During the flaring events, the object was also monitored in X-rays with Swift-XRT. In the second panel, we show the XRT light curve for energy 0.3–10 keV. We observed high flux variability during flare 3, but during flare 1 and the quiet period, the flux did not vary much. During flare 2, we do not have any X-ray observations. The flaring behaviour in optical-UV is much clearer compared to X-ray, and the UVOT light curve is shown in panel 3. A close temporal correlation between the Fermi light curve and the UVOT light curve is seen, suggesting that gamma-ray and optical emissions are highly correlated and produced at the same time. A detailed study of the correlation is presented later in this paper. The archival optical light curves from ZTF (g, i, r-bands) and ASAS-SN (g-band) are shown in panels 4 and 6. They clearly follow the UVOT light curve and are in temporal correlation with the gamma-ray emission. In panel 5, we show the archival near-infrared emission from the WISE observatory. Unfortunately, we do not have a clear temporal correlation of gamma-ray with WISE band emission, but we see some correlation with optical band emission. The overall suggestion is that the broadband light curves are correlated temporally and might have been produced at the same location and at the same time. This information is essential for broadband SED modelling to derive the main physical mechanism responsible for their emissions. We have also shown the photon index variation with gamma-ray and x-ray fluxes, respectively in Fig. 4.

Figure 2.

Broadband light curves of Ton 599 during the flaring episodes. The top panel shows the Fermi-LAT light curve, followed by the X-ray and UVOT light curves in panels 2 and 3. The archival ZTF, WISE, and ASAS-SN light curves are shown in panels 4,5, and 6.

Figure 3.

The local peak in the $\gamma$ -ray light curve for each flare is fitted with the Sum of the exponential function. In flare 2, there was not sufficient observation, so we tried to fit it with the rising and decay function. The reduced $\chi^2/\text{ndf}$ values are calculated to estimate the goodness-of-fit and are mentioned in each plot.

3.3. Fractional variability

Blazars show strong flux variability at all frequencies. Fractional variability allows for the comparison of variability amplitudes across the entire electromagnetic spectrum and can be calculated using the relation given in Vaughan et al. (Reference Vaughan, Edelson, Warwick and Uttley2003).

(2) \begin{equation}F_{\text{var}} = \sqrt{\frac{S^2 - \sigma^2}{r^2}},\end{equation}

(3) \begin{equation}\text{err}(F_{\text{var}}) = \sqrt{\left(\sqrt{\frac{1}{2N}}\cdot \frac{\sigma^2}{ r^2 F_{\text{var}}} \right)^2 + \left( \sqrt{\frac{\sigma^2}{N}}\cdot\frac{1}{r} \right)^2},\end{equation}

where $\sigma^2_{\text{XS}} = S^2 - \sigma^2$ , is called excess variance, $S^2$ is the sample variance, $\sigma^2$ is the mean square uncertainties of each observation, and r is the sample mean.

The plot of $F_{\text{var}}$ versus frequency is shown in Fig. 5 to illustrate its overall shape. We have created plots for three segments total dataset, Flare 1 and Flare 3. However, due to the limited number of data points, Flare 2 has not been plotted. All segments show a similar shape; $\gamma$ -ray shows the highest variability, followed by optical and UV bands, and then X-ray. A similar result is also obtained for other blazars by Schleicher et al. (Reference Schleicher2019). In literature, studies have also shown the increasing or decreasing F $_{var}$ with respect to frequency, suggesting either the large number of particles injected in the jet (Prince Reference Prince2020) or the presence of steady thermal emission (Bonning et al. Reference Bonning2009).

Table 3.

Results of temporal fitting with sum of exponentials.

3.4. X-ray spectral fitting

During the selected period for study in this work, we found 3 observations of XMM-Newton and 2 observations of NuSTAR, which were simultaneous to two of the XMM observations. The spectra were extracted following the standard procedure for both telescopes and were fitted with a single power-law model. The best-fit parameters are shown in Table 5, and the corresponding plots are shown in Fig. 6. The XMM-Newton observation performed in 2019 happens to be in a low-flux state, which is also visible in Fig. 1. The reduced chi-square estimated for these observations suggests that the power-law is a best-fit model. However, for observation done in June 2021, which happens to be in the high gamma-ray flux state, a single power-law does not provide a better fit; the reduced chi-square is very high. We tried different combinations to fit this XMM-Newton spectrum, such as a single log-parabola, a combination of power-law and log-parabola, a combination of black-body (bbody) and log-parabola, and a combination of power-law and black-body. Finally, we found that the power-law + bbody produces a better fit with a better-reduced chi-square. The estimated photon index is $\Gamma$ = 1.51 $\pm$ 0.05, black-body temperature = 0.10 $\pm$ 0.01 keV, and the chi-square is 422.57/435. This suggests that during the gamma-ray flaring state and the X-ray flaring state (Flare-1), the X-ray spectrum has some influence on black-body emission from the accretion disc, linking to a possible accretion-disc connection. A possible accretion disc connection is also suggested from the flux distribution and the PSD analysis in the next few sections. We have also produced the NuSTAR spectrum taken during Flare-1 of the gamma-ray and in the low flux state. These two spectra are used in broadband SED modelling to guide the model better in order to produce the best fit for the data and derive the physical parameters.

Figure 4.

$\gamma$ -ray Photon Index ( $\alpha$ ) and curvature parameter ( $\beta$ ) vs gamma-ray flux and X-ray photon index vs X-ray flux.

Figure 5.

Fractional Variability estimated for various wavebands for flaring state, flare 1, 2, and total dataset.

Table 4.

Fractional variability.

3.5. Gamma ray emission region

The gamma-ray light curves are produced with 6 h of binning, which still provides a significant TS value for each data point. We used this light curve to derive the fastest variability time scale using the following expression

(4) \begin{equation}F(t_2) = F(t_1) \cdot 2^{{(t_2 - t_1)}/{t_d}}\end{equation}

Where $F({t_1})$ and $F({t_2})$ are the fluxes measured at time ${t_1}$ and ${t_2}$ respectively, and ${t_d}$ represents doubling timescale or variability time of flux. A range of variability time is found, from a few hours to a few days. The shortest one is recorded as 2.5 h, which is used as the fastest variability time to estimate the size of the emission region, using the following relation

(5) \begin{equation}R \leq c \, t_{\text{var}} \,\delta (1 + z)^{-1},\end{equation}

Where $z=0.72$ is redshift (Hewett & Wild Reference Hewett and Wild2010; Schneider et al. Reference Schneider2010) and $\delta $ is the Doppler factor. The size of the emission region is found to be 5.0 $\times 10^{15}$ cm, for $\delta $ = 18.2 (Rajput et al. Reference Rajput2023). The observational constraint on the size of the emission region is important because it helps to derive the best-fit model for the broadband SED. As expected, the shortest variability time is produced by a smaller region close to the base of the jet.

Table 5.

The XMM-Newton spectra are fitted with a simple power-law model. The flux is in units of 10 $^{-12}$ erg cm $^{-2}$ s $^{-1}$ .

Figure 6.

The XMM-Newton spectra fitted with a simple power-law model. To fit the observation during a high gamma-ray flux state, disc blackbody is required (lower panel).

3.6. $\gamma$ -ray spectral analysis: Locating the gamma-ray emission region

The location of the gamma-ray emission region is important to constrain the broadband SED, and this helps us to decide which external photon fields need to be used for the inverse-Compton scattering. To have an idea about the location of the emission region, we produced the gamma-ray spectral data points as shown in Fig. 7. As seen by the naked eye, the spectrum estimated for various flux states is curved in nature, and we fitted all the spectra with the log-parabola function to derive the best fit. We found that the log parabola fits the spectrum very well, and the best-fit parameters are shown in Table 6. The datapoints with arrows are the upper limits and have not been included in the fitting. In almost all the cases, the spectrum seems to bend around 10 GeV, suggesting that photons above 10 GeV are getting absorbed and hence reduced in number. A similar spectrum has been seen across various FSRQs and BL Lacs during the flaring events and also during the 2018 flares of Ton 599 (Prince Reference Prince2019).

Figure 7.

The gamma-ray SED derived from various states. Upper panel: Flare 1 and Quiet State. Lower panel: Flare 2 & 3.

The presence of curvature in the gamma-ray SED plays a crucial role in constraining the location of the emission region. A curvature or break in the gamma-ray spectrum is interpreted as a signature of photon-photon absorption (pair-production) where a gamma-ray photon interacts with a low-energy photon from the BLR, suggesting the emission region is possibly within the BLR. It has been shown that at the base of the jet, the medium is quite opaque for photons having energy above 20 GeV (Liu & Bai Reference Liu and Bai2006), and hence a curvature or break is seen in that energy range.

The curvature or the cut-off in the spectrum can happen because of other reasons as well when there is already a cut-off in the energy distributions of the particles. The shape of the initial particle distribution injected in the jet remains unknown, and it can have a shape like power-law/broken power-law/log-parabola.

Table 6.

Results of gamma-ray SEDs fitted with spectral type Log-Parabola(LP).

While performing the broadband SED modelling, we keep the emission region close to BLR or within the range of the inner and outer radii of the BLR to maximise the BLR’s contribution.

3.7. Flux distribution

The blazar light curve shows a range of variability on shorter to longer time scales. The cause of the variability can be accessed through the flux distribution study. The flux distribution helps us to probe the nature of variability, whether it is caused by additive or multiplicative processes. A Gaussian flux distribution leans toward an additive process, suggesting stochastic and linear variations (McHardy Reference McHardy and Belloni2010). On the other hand, if the stochastic variations are non-linear, then the log-normal flux distribution is expected, which is nothing but the Gaussian distribution of the logarithmic flux values (Uttley, McHardy, & Vaughan Reference Uttley, McHardy and Vaughan2005; McHardy Reference McHardy and Belloni2010). The log-normal distributions are very common in AGNs, Gamma-ray Bursts, and the galactic X-ray binaries (Uttley & McHardy Reference Uttley and McHardy2001; Negoro & Mineshige Reference Negoro and Mineshige2002; Shah et al. Reference Shah2018a; Khatoon et al. Reference Khatoon, Shah, Misra and Gogoi2020; Prince, Khatoon, & Stalin Reference Prince, Khatoon and Stalin2021) where the emission is dominated by the accretion disc. However, in the case of blazars, the situation is completely different, where most of the emission is dominated by the jets. However, in some of the studies, people have shown that the disc and jet can have some possible connection in the case of a blazar when explored systematically and carefully. In this case, the idea is that the fluctuations or variability produced in the accretion disc can somehow travel to jets and modulate the jet variability, leaving its imprint on the emissions produced in the jets. Following this scenario, it is expected that the gamma-ray emission, which is surely produced in the jet, can have a log-normal flux distribution as expected in the case of emission produced in AGN or X-ray binaries. In a gamma-ray light curve of blazars Shah et al. (Reference Shah2018a) & Prince et al. (Reference Prince, Khatoon and Stalin2021) have found that a single log-normal flux distribution is quite prevalent.

However, to counter that, there have been various suggestions that the log-normal flux distribution in the case of a blazar can be produced by the mini-jets model. The mini-jets model suggests that the jets can have many mini-jets and that the total emission is the combination of these mini-jets, which leads to the multiplicative nature of the total emission and results in a log-normal flux behaviour (Biteau & Giebels Reference Biteau and Giebels2012). On the other hand, Scargle (Reference Scargle2020) argued that the log-normal flux distribution does not necessarily mean the production of a multiplicative process.

We estimated the flux histogram of the total light curve of Ton 599 from 2020 to 2024. In the analysis, we choose only the flux data points that have TS $\gt$ 9 to account for the highly significant data points in order to achieve the correct representation. To investigate the behaviour of flux distribution, we performed the Anderson-Darling(AD) test on the logarithm of flux. We have performed the AD test for both the normal and log-normal distribution functions. The p-value for normal distribution is found to be 0.007 with AD statistics as 1.098, and for the log-normal distribution p-value is found to be 8.928 $\times$ 10 $^{-22}$ with AD statistics as 8.979. So, these two distributions are not suitable since the p-value is below 0.05. Next, we fit the histogram with the double Gaussian probability density functions (PDF), and the functional form is given by Sharma et al. (Reference Sharma, Kamaram, Prince, Khatoon and Bose2023) and Khatoon et al. (Reference Khatoon, Shah, Misra and Gogoi2020).

(6) \begin{multline}f(x) = \frac{a}{\sqrt{2\pi \sigma_1^2}} \exp \left( - \frac{(x - \mu_1)^2}{2 \sigma_1^2} \right)\\+ \frac{(1 - a)}{\sqrt{2\pi \sigma_2^2}} \exp \left( - \frac{(x - \mu_2)^2}{2 \sigma_2^2} \right)\end{multline}

Figure 8.

The figure shows double log-normal fits of the $\gamma$ -ray flux histogram. The blue line represents a double log-normal fit. The parameter for fit is shown in Table 7.

The observed flux distribution is fitted with the double Gaussian function, and it is shown in Fig. 8, and the best-fitted parameters are shown in Table 7. We found that a bi-model flux distribution can be well-fitted with a log-normal flux distribution, suggesting the variability of non-linear in nature. As we expect, the gamma-ray emissions are produced in the jet far from an accretion disc, but if the observed jet variability is non-linear in nature or best represented with a log-normal flux distribution, it suggests a possible connection between the accretion disc and the jet. It has been shown that it is quite possible, and in a few of the blazars (Sharma et al. Reference Sharma, Kamaram, Prince, Khatoon and Bose2023), it has been shown that the variation or fluctuations produced in the accretion disc can travel to the jet and modify the jet variability accordingly. In the literature, it has also been noticed that some blazars show a double log-normal flux distribution instead of a single log-normal distribution (Shah et al. Reference Shah2018a; Kushwaha et al. Reference Kushwaha2016). Our study concludes that there is a possibility of disc-jet coupling in this source. The nature is clearer when we look at the XMM-Newton spectra, where the thermal disc component dominates over the non-thermal emission.

Table 7.

Results of the double Gaussian fit on the log-flux data for flux distribution.

3.8. Colour-Magnitude (CM) variation

The colour-magnitude (CM) diagram can be measured between various optical-UV and IR filters. It can be used as a tool to study the IR-optical-UV emissions of blazars. The most common trend that has been seen among blazars is ‘redder-when-brighter’ or ‘bluer-when-brighter’ depending upon their types, but sometimes, a complex nature has also been seen where the trend is not clear. Bonning et al. (Reference Bonning2012) have shown that mostly FSRQs follows ‘redder-when-brighter’ whereas Ikejiri et al. (Reference Ikejiri2011) have shown that the BL Lac in general show ‘bluer-when-brighter’ trend.

We can also derive the optical spectral index by following the equation (Wierzcholska, Alicja et al. Reference Wierzcholska and Ostrowski2015).

(7) \begin{equation}\alpha_{UB} = 0.4 \cdot \frac{U - B}{\log_{10} \left( \nu_U / \nu_B \right)},\end{equation}

Here, $(U-B)$ represents the colour index derived from the magnitudes in the U and B bands, while $\nu_U$ and $\nu_B$ denote the effective frequencies for these respective bands. The scaling factor in the numerator accounts for the differences between the bands (Bessell, Castelli, & Plez Reference Bessell, Castelli and Plez1998).

In Fig. 9, we show the possible colour-index variations of combinations of WISE W1 and W2 bands, ZTF g and r bands, and Swift optical (U, B, V) bands. In the case of WISE bands, a very small number of observations are available, but they cover both the high and low flux states, as can be seen in Fig. 2, and as a result, we see a mild trend of ‘redder-when-brighter’ with a positive correlation of $r = 0.12$ . The mild trend could be because of the lack of data points in the W1 and W2 bands of WISE. Swift optical bands $(B-V)$ and $(U-V)$ also show similar characteristics of ‘redder-when-brighter’ but with a much stronger positive correlation of $r = 0.65$ and $0.99$ , respectively. In the case of Swift $B-V$ , we see a mild ‘bluer-when-brighter’ with a negative correlation coefficient of $r = -0.17$ . A similar trend of ‘bluer-when-brighter’ with a strong correlation coefficient ( $r =-0.81$ ) is observed in ZTF g-r. This shows that the object behaves differently in different wavebands or shows a complex nature, as also reported in (Safna et al. Reference Safna, Stalin, Rakshit and Mathew2020).

Figure 9.

Colour index variations of possible IR and optical combinations from WISE, ZTF, and Swift.

In the community, the understanding is that the emission from the accretion disc is mostly bluer (bluer than the optical-UV of the synchrotron), and the jet emission produced by the Compton scattering is mostly redder (Ghisellini Reference Ghisellini2013; Sarkar et al. Reference Sarkar2019).

Our observation of mixed trends in optical and IR suggests that it is difficult to separate the optical-IR emission from the synchrotron and the accretion disc, and therefore, we do not see a clear trend of one behaviour.

3.9. Correlation

The broadband correlation is an important way to understand if the multi-wavelength emissions are produced simultaneously at the same location or if they are produced at separate locations with some time delays (Fuhrmann et al. Reference Fuhrmann2014). This information is important when we do the broadband SED modelling, which helps us to decide if one should choose a one-zone or multi-zone emission region (Patel et al. Reference Patel2018). This also helps us to understand what exactly is happening inside the jets (Pushkarev, Kovalev, & Lister Reference Pushkarev, Kovalev and Lister2010; Cohen et al. Reference Cohen2014).

The broadband coverage allows us to explore the correlation among various emissions coming out of the jet. We used the traditional z-transformed discrete correlation function (zDCF) to find out the correlations among various wavebands and the possible time lags, if any are present. Details about the code can be found in (Alexander Reference Alexander2013). A Fortran version of the code is available online.Footnote ^j The main reason for using zDCF over the discrete correlation function of Edelson & Krolik (Reference Edelson and Krolik1988) is that it corrects several biases using equal population binning and Fisher’s z-transform. These lead to a more robust and powerful estimate of cross-correlations among sparse light curves. If the two light curves, say LC1 and LC2, are cross-correlated, then the positive time lags indicate that the LC1 light curves lead the LC2 and vice-versa. A zero time lag suggests the co-spatial origin of the emissions (Pushkarev et al. Reference Pushkarev, Kovalev and Lister2010; Raiteri et al. Reference Raiteri2013; Cohen et al. Reference Cohen2014; Sarkar et al. Reference Sarkar2019).

The gamma-ray, optical from ZTF, and the X-ray light curves are well sampled, and therefore, we choose these three frequencies for the correlation study. We focus on the total light curve rather than the individual flaring episode because we do not have a sufficient number of data to derive any significant correlations in individual flares. The correlation results are presented in Fig. 10. The correlation coefficient estimated between gamma-ray and the optical-g (ZTF) band at zero time lag is $\sim$ 0.7 ( $\sim$ 70%), suggesting the gamma-ray and optical emissions are highly correlated and are produced at the same location mostly simultaneously. The gamma-ray versus X-ray correlation shows a similar behaviour with zero time lag, and the correlation coefficient is estimated as $\sim$ 0.6 ( $\sim$ 60%), suggesting gamma-ray and X-ray are also highly correlated and are produced at the same location. Combining both results, we conclude that the gamma-ray, X-ray, and optical emissions are produced mostly simultaneously at the same location. Therefore, a single-zone emission model is sufficient to model the broadband SED. Following these arguments, we proceed with one-zone emission modelling of the broadband SED as discussed in section 3.11. Similar results were also seen in Prince (Reference Prince2019) during the flare of 2018 and in the recent flare by Rajput et al. (Reference Rajput2023).

Figure 10.

The cross-correlation estimated for various combinations among gamma-ray, optical, and X-ray emissions.

3.10. Power spectral density

The variability in the source light curve can also be quantified by the power spectral density (PSD), which determines the amplitude of variation in the temporal light curve as a function of Fourier frequency or variability time scales (Ryan et al. Reference Ryan, Siemiginowska, Sobolewska and Grindlay2019). PSD is important to understand the average properties of the variability, whereas the source light curve could be thought of as only a single realisation of an underlying stochastic process, as shown by Vaughan et al. (Reference Vaughan, Edelson, Warwick and Uttley2003).

The blazar light curve shows aperiodic flux variations across the wavebands over the short as well as the longer time scales. The PSD can be used as a tool to derive the characteristic time scale in the aperiodic light curves, which will correspond to the variability time scale in the system, and further, it can help us to constrain the size of the emitting zone. The characteristic time scale in the system represents the breaking time in the PSD when the PSD is best fitted by the bending or broken power law rather than a simple power law. The breaking time scale in the light curves is identified as the time scale of variability in the source or the particle cooling or escape time scales (Kastendieck, Ashley, & Horns Reference Kastendieck, Ashley and Horns2011; Sobolewska et al. Reference Sobolewska, Siemiginowska, Kelly and Nalewajko2014; Finke & Becker Reference Finke and Becker2014; Chen et al. Reference Chen, Pohl, Böttcher and Gao2016; Kushwaha et al. Reference Kushwaha, Sinha, Misra, Singh and de Gouveia Dal Pino2017; Chatterjee et al. Reference Chatterjee, Roychowdhury, Chandra and Sinha2018a; Ryan et al. Reference Ryan, Siemiginowska, Sobolewska and Grindlay2019; Bhattacharyya et al. Reference Bhattacharyya, Ghosh, Chatterjee and Das2020).

In general, the observed PSD can be fitted with a single power law model, which can be defined as $P(\nu) \propto \nu^{-\beta}$ , where $\beta$ is the slope. Depending upon the value of the slope parameter, the earlier results suggest that the variability in high energy bands (X-ray and gamma-ray) is characterised by pink or flicker noise ( $\beta=1$ ) (Abdo et al. Reference Abdo2010a; Isobe et al. Reference Isobe2014; Abdollahi et al. Reference Abdollahi2017) and in lower energy (radio and optical) by damped/red-noise type processes ( $\beta=2$ ) (Max-Moerbeck et al. Reference Max-Moerbeck2014; Nilsson et al. Nilsson, Lindfors, Takalo, Reinthal, Berdyugin, Sillanpä 2018). The $\beta$ $\sim$ 0 interpreted as an uncorrelated white-noise-type stochastic process (Press Reference Press1978). However, Press (Reference Press1978) also interpreted $\beta$ $\sim$ 1 & 2 as a flicker (or pink-noise) and damped random-walk (or red-noise) type correlated stochastic processes, respectively.

To probe the characteristic variability time scale and the type of variability in Ton 599, we produced the gamma-ray light curve from August 2008 to August 2024 ( $\sim$ 16 yr). We have produced the power spectrum using the discrete Fourier transform following the Vaughan et al. (Reference Vaughan, Edelson, Warwick and Uttley2003) and then fit it with a simple power law. The observed and fitted PSD is shown in Fig. 11. The best-fitted slope ( $\beta$ ) is estimated as $\beta=1.16$ , suggesting a pink-noise kind of stochastic variability in the light curve of Ton 599. Pink noise or flickering represents that the light curve has more power in the short-term variability (Vaughan et al. Reference Vaughan, Edelson, Warwick and Uttley2003). We do not see any signature of curvature or break in the power spectrum, suggesting a much longer characteristic time scale is involved in the gamma-ray variability of Ton 599 and, therefore, a much longer baseline light curve is needed to probe that. With the Fermi, which is still in operation, we hope Ton 599 will be monitored continuously for a much longer duration, and possibly the break in the PSD can be observed.

Figure 11.

The Power Spectral density (PSD) derived for the total gamma-ray light curve from August 2008 to August 2024. The continuous red line shows the best fit to the PSD with slope, $\unicode{x03B2}=1.16$ .

3.11. disc-Jet coupling

In the literature, Blandford & Znajek (Reference Blandford and Znajek1977) & Blandford & Payne (Reference Blandford and Payne1982) theorise that the jet can be collimated or launched via the extraction of rotational energy from the supermassive black hole and via the extraction of power from the accretion disc in the presence of magnetic fields, which shows that the jet and the accretion disc in some way are connected. The observational evidence suggests that a blazar hosts mostly a thin accretion disc and strong relativistic jets. Indeed, most of the emissions that we observe in the blazar are highly dominated by the jets, and the disc emissions mostly hide behind them. However, in some cases, especially in FSRQs, people have observed strong disc emission (blackbody emission) along with high-energy gamma rays. The concept of disc-jet coupling describes the relationship between the inflow of material through an accretion disc around a black hole and the simultaneous outflow in the form of relativistic jets (Fender & Belloni Reference Fender and Belloni2004). This coupling plays a critical role in understanding black hole X-ray binaries, AGNs, and gamma-ray bursts.

Observationally, the disc-jet coupling can be understood via various studies such as detecting the break in PSD (McHardy Reference McHardy2008; Chatterjee et al. Reference Chatterjee, Roychowdhury, Chandra and Sinha2018a) and comparing the time scale of PSD with the accretion disc time scale as done in Sharma et al. (Reference Sharma, Kamaram, Prince, Khatoon and Bose2023). The accretion disc variability is expected to produce the log-normal flux distribution (Uttley & McHardy Reference Uttley and McHardy2001), and hence, the gamma-ray flux distribution can be investigated to obtain an indication of whether the jet emission is somehow linked with the disc variability.

We collected the XMM-Newton observation during the low and high flux states based on gamma-ray flux as marked in Fig. 1. We analysed and fitted them using a simple power-law model. Out of the three, two spectra are well fitted with a power law; however, the third one, taken exactly during the gamma-ray flaring state, required an additional model of the accretion disc to fit the data. The X-ray spectral fitting of XMM-Newton spectra is best fitted by a combination of power-law + bbody, the X-ray spectrum has the influence of black-body emission from the accretion disc, which suggests that there is a possible connection between the disc-jet which leads to the appearance of disc emission in the X-ray spectra in blazar which is very rare since the jet emission is highly dominant over the disc.

To confirm the disc-jet coupling in Ton 599, we produced the PSD using the longest gamma-ray light curve available, and we found that a single power law can fit the PSD, suggesting that an even higher baseline is required to see the break in the PSD. Next, we produced the flux distribution to see if the disc variability has some influence on the gamma-ray emission. We found a bi-modal flux distribution that can be fitted with a log-normal distribution, suggesting non-linear variability in the gamma-ray emission from the jet. Given that these emissions are typically expected to originate in jets far from the accretion disc, the observed variability is non-linear and suggests that the variability produced in the accretion disc travels to the jet through an unknown process and leads to a possible connection between the accretion disc and the jet. It has been shown earlier that it is quite possible, and in a few of the blazars (Sharma et al. Reference Sharma, Kamaram, Prince, Khatoon and Bose2023), it has been seen that the variation or fluctuations produced in the accretion disc can travel to the jet and modify the jet variability accordingly. Furthermore, it has also been noticed that some blazars show a double log-normal flux distribution instead of a single log-normal distribution (Shah et al. Reference Shah2018a; Kushwaha et al. Reference Kushwaha2016). Based on the modelling of the XMM-spectrum and flux distribution study, we conclude that there is a possibility of disc-jet coupling in the blazar Ton 599.

3.12. Spectral energy distribution (SED)

The emission mechanisms of blazars can be better understood through the modelling of their broadband spectral energy distributions (SED). We have collected broadband data from various space-based and ground-based observatories to construct the broadband SED. SED modelling provides insight into the real physical scenario happening close to the SMBHs. We choose to model all the identified flaring states as well as quiet or low flux states to understand which jet parameters cause the flaring events. To perform the SED modelling, we used the publicly available code JetSet Footnote ^k (Tramacere et al. Reference Tramacere, Giommi, Perri, Verrecchia and Tosti2009; Tramacere, Massaro, & Taylor Reference Tramacere, Massaro and Taylor2011; Tramacere Reference Tramacere2020). JETSET (version 1.3.0) fits numerical models to the data to identify the optimal parameter values that most accurately represent the observed data. We obtained the broadband SED using the data from Fermi-LAT, Swift-XRT, Swift-UVOT, ZTF, ASAS-SN, and NuSTAR for each flare and quiet period as shown in Fig. 2, except for Flare 2 due to lack of observational data in that period.

The broken power law (bkn) model was considered for electron distribution with a lower energy spectral slope be p, a high energy spectral slope be $p_1$ , and a turn-over energy be $\gamma_{\text{break}}$ . The functional form is shown in Rajput et al. (Reference Rajput2023). The JetSet uses this particle distribution to solve the transport equation and derive the photon flux due to synchrotron and various cases of inverse-Compton scattering. The JetSet has a large set of parameters in the model that will be used for modelling the SED. To reduce the number of free parameters in the model, we fixed some of the parameters such as the inner and outer radii of the BLR to standard values ( $R_{BLR, in} = 2\times10^{17}$ cm; $R_{BLR, out} = 1\times10^{18}$ cm). The accretion disc temperature is fixed at the standard value for the source $T_{disc}$ = 1 $\times10^5$ K (Patel & Chitnis Reference Patel and Chitnis2020) and the luminosity is also fixed to the standard value, $L_{disc}=1\times10^{45}$ erg/s (Ghisellini et al. Reference Ghisellini2010). The viewing angle of the jet is fixed as 2 deg (Pushkarev et al. Reference Pushkarev, Kovalev, Lister and Savolainen2009) and the redshift of the source to $z = 0.72$ . Through the modelling, we want to explore if the jet is particle-dominated or magnetic field-dominated. We also choose the number of cold protons in the jet that might be contributing to the production of total jet power, and for that, we consider the ratio of cold protons to relativistic electrons to be 0.1, which is fixed during fitting (Ghisellini Reference Ghisellini2012). We optimize the other parameters such as magnetic field, particle distribution slopes, and energy ( $\gamma_{\min}$ , $\gamma_{\max}$ ), the BLR optical depth, size of the emission region, location of the emission region, and Doppler boosting to achieve a best-fit model for the broadband SED. The best-fit parameters are shown in Table 8, and the corresponding fits are shown in Fig. 12.

Table 8.

The parameters obtained from modelling the multifrequency SEDs of flares 1, 3, and the quiet period using JetSeT. The viewing angle is fixed at $\unicode{x03B8}$ =2 deg from Pushkarev et al. (Reference Pushkarev, Kovalev, Lister and Savolainen2009). The parameters with (^*) are fixed during modelling.

Figure 12.

The broadband SEDs of Flare 1, 3 & the quiet period fitted with one-zone leptonic model. The data and the various colourful lines are self-explanatory.

Our modelling reveals that to fit the broadband SED, an external field is required, such as the BLR and the accretion disc. These two act as external photon fields to provide a sufficient amount of photons to produce the gamma-ray emission. We noticed that higher photon densities are required during the flaring episode compared to the quiet state. We successfully fit the IR, optical, and UV spectra with the synchrotron model, and the required magnetic fields are the following: during flares 1 & 3 the required magnetic field is almost the same (0.83 & 0.87 Gauss respectively) and during the low flux state the magnetic field is a bit lower with the value of 0.59 Gauss. However, the magnetic field estimated in Rajput et al. (Reference Rajput2023), Manzoor et al. (Reference Manzoor, Shah, Sahayanathan, Iqbal and Dar2024) is a bit higher, with values ranging between 1.5–2.0 Gauss. Another interesting thing we noticed is that we need higher Doppler boosting to produce the bright gamma-ray flaring states (Flare 1 & 2) compared to the low flux states. The best-fitted Doppler factor for Flare 1 & 2 is 25 and 28, respectively, while for the quiet state, it is around 19, suggesting a sudden rise in Doppler factor can produce bright gamma-ray flaring episodes. The size of the emission region is smaller during the flaring episodes compared to the quiet state, which is expected because a smaller emission region generally produces short-term flares.

The energy densities in the accretion disc ( $U'_{disc}$ ), BLR ( $U'_{BLR}$ ), electrons ( $U'_e$ ) and magnetic field ( $U'_B$ ) are calculated directly by the model using JetSeT. The BLR and disc photon densities are much higher during the flaring case compared to the low flux state case, therefore supplying more photons for inverse-Compton scattering and eventually producing more gamma-ray emission. The magnetic energy density is higher in the flaring case, producing more synchrotron and SSC. The SSC plays an important role here in fitting the X-ray spectral points in the SED.

The estimated jet power in the electron and magnetic fields appears to be higher during the flaring case compared to the low flux (quiet) state, which is expected. The total jet power is estimated by summing over the electron and the magnetic power carried by the jet, and it is found to be of the order of 10 $^{45}$ erg/s. Comparing this with the Eddington luminosity of the source, which is estimated as L $_{Edd}$ = 10 $^{38}$ (M/M $_\odot$ ) erg/s. In the various estimates, the mass of the BH is estimated between $(0.79-3.47)\times10^8\,{\rm M}_{\odot}$ (Xie, Zhou, & Liang Reference Xie, Zhou and Liang2004; Liu, Zhao, & Wu Reference Liu, Zhao and Wu2006). In our case, we take the black hole mass of Ton 599 to be around 1 $\times$ 10 $^{8}$ M $_\odot$ , which gives the L $_{Edd}$ = 10 $^{46}$ erg/s, which is much higher than the total jet power. We conclude that the flaring events are produced by the sudden increase in the Doppler factor and the magnetic field, and the total jet power produced through this process is below the Eddington luminosity.

Patel & Chitnis (Reference Patel and Chitnis2020) models the broadband SED of the flare of TON 599 during November 2017 using the external component (EC). Their modelling suggests the requirement of a two-component leptonic emission model to reproduce the observed broadband SED. The authors have chosen two specific periods corresponding to low and high gamma-ray flux states to compare the jet parameters. The low flux state, in their case, is well fitted with a single leptonic emission model where the GeV emission is explained by the EC-BLR. However, during the flaring state, the SED is modelled with a two-zone emission model; one is located within the BLR, and the second is outside the BLR but within the DT. But in our study, we found that to explain the broadband SED during the flaring state, a one-zone emission model is sufficient, where the emission region is located within the BLR and the seed photons from the accretion disc and the BLR can act as seed photons for the external Compton process. Their modelling suggests that during the flaring event, the jet has become magnetically dominant and, as a result, has more magnetic power than electrons. Our study suggests the jet is more particle-dominated and hence has more jet power in electrons compared to the magnetic field. The magnetic field value estimated in our case is more or less consistent with their estimate. The main difference between our results is that during the flaring state, the Doppler factor has increased significantly from $\sim$ 19 in the low state to $\sim$ 28 during the flare (i.e., especially flare 3). Other differences are in the modelling as we consider the size of the BLR as a shell rather than a single boundary.