On the origin of multifrequency temporal and spectral variability in Ton 599

Abstract

In this work, we studied the broadband temporal and spectral properties of the flat-spectrum radio quasar Ton 599. We collected the long-term data from January 2019 to August 2024 when the source was in a long flaring episode. We used the Bayesian block methodology to identify the various flux states, including three flares. The broadband fractional variability is estimated during two flaring states. The F$_{\text{var}}$ variation with respect to frequency shows a nearly double hump structure similar to broadband SED. The power spectral density shows a pink-noise kind of stochastic variability in the light curve, and we do not see any break in the power spectrum, suggesting a much longer characteristic timescale is involved in gamma-ray variability. The flux distribution is well-fitted with a double log-normal flux distribution, suggesting the variability of non-linear in nature. The gamma-ray, optical, and X-ray emissions were found to be highly correlated with a zero time lag, suggesting a co-spatial origin of their emissions. We used the one-zone leptonic model to reproduce the broadband spectrum in the energy range from the IR to very high-energy gamma rays. The increase in the magnetic field and the Doppler factor were found to be the main causes for high flux states. The XMM-Newton spectra taken during one of the flaring durations exhibit a signature of thermal black body emission from the accretion disc, suggesting a possible disc-jet coupling. This has also been indicated by the gamma-ray flux distribution, which shows the distribution as non-linear in nature, which is mostly seen in galactic X-ray binaries or active galactic nuclei, where the accretion disc dominates the emission.

1. Introduction

Active galactic nuclei (AGNs) are the nuclei of galaxies hosting supermassive black holes (SMBHs) at the center, which are active, meaning mass is accreting onto the black hole. AGN consists of three main parts: a central SMBH, a thermal accretion disc, and the relativistic jets in pairs produced perpendicular to the disc plane. There is also observational evidence of gas clouds hanging above the accretion disc, known as the broad-line region (BLR), and dusty molecules having a torus-like structure to hide the central part of the AGN. AGNs are randomly oriented in the Universe, where the jet points randomly. AGNs with one of the jets pointing toward the Earth are known as blazars (Urry & Padovani 1995). The jets are highly relativistic in nature, which boosts all the emissions produced along the jet axis. The observational results suggest that they emit across the entire electromagnetic spectrum, ranging from low-energy radio waves to very high-energy $\gamma$ -rays (Ulrich, Maraschi, & Urry 1997; Hovatta & Lindfors 2019). Some blazars were also found to be in a good temporal and spatial correlation with the neutrino emissions or events detected by the IceCube observatory (IceCube Collaboration et al. 2018; Prince et al. 2023). Blazars exhibit a very strong flux variability on the timescale from minutes to years across the entire electromagnetic spectrum (Ulrich et al. 1997; Aharonian et al. 2007; Raiteri et al. 2013; Hovatta & Lindfors 2019; Goyal et al. 2022). The broadband information is used to construct the broadband spectral energy distribution (SED) which shows two broad emission components (Marscher 1980; Konigl 1981; Sambruna, Maraschi, & Urry 1996; Fossati et al. 1998; Abdo et al. 2010b), with the low-energy component peaking in the optical to X-ray range, which is explained by synchrotron emission from relativistic electrons travelling in the magnetic fields of the jet (Rawes, Worrall, & Birkinshaw 2015; Urry & Mushotzky 1982). The high energy component peaks at MeV to TeV energy range and it is believed to be the result of inverse Compton (IC) scattering of low-energy photons by the relativistic electrons within the jet (synchrotron-self Compton; SSC; Sikora et al. 2009) or outside the jets (external Compton; EC; Dermer, Schlickeiser, & Mastichiadis 1992; Sikora, Begelman, & Rees 1994). The observation suggests that the peak of the low and high energy components in the SED also changes from source to source and sometimes also with the various flux states of the source. Blazars also have a sub-category of flat spectrum radio quasars (FSRQs) and BL Lacerate objects (BL Lacs) depending upon the presence and absence of broad emission lines in their optical spectra. Further, based on the location of the synchrotron peak blazars are classified into different sub-classes (Padovani & Giommi 1995) such as low synchrotron peak blazar (LSP for peak frequency below 10 $^{14}$ Hz), high synchrotron peak blazar (HSP for peak frequency above 10 $^{15}$ Hz), and intermediate synchrotron peak blazar (ISP for peak frequency between 10 $^{14}$ to 10 $^{15}$ Hz) (Abdo et al. 2010c). Separately, the BL Lacs objects were also classified into various categories depending on the location of the synchrotron peak. The classes are low BL Lac (LBL for peak frequency below 10 $^{14}$ Hz), intermediate BL Lacs (IBL for peak frequency between 10 $^{14}$ –10 $^{15}$ Hz), high BL Lacs (HBL for peak frequency above 10 $^{15}$ Hz), and extreme high BL Lacs (EHBL where peak frequency is above 10 $^{16}$ Hz).

It has also been shown that the jet emission is not entirely produced in the jet but also has some influence from the accretion disc in some of the AGNs and blazars. In Chatterjee et al. (2018a), authors have shown the possible accretion disc origin of X-ray variability based on the X-ray PSD breaks derived from the combination of soft-X-ray telescope (SXT), LAPXC, and Swift-XRT. A moderate and strong correlation of disc and jet power is also obtained for a sample of blazars in Rajguru & Chatterjee (2022, 2024), suggesting disc-jet coupling in blazars is more common than we think. Based on the gamma-ray power spectra Sharma et al. (2023) have shown a possible disc-jet connection in three bright blazars (S4 0954+65, PKS 0903-57, and 4C +01.02). During the current flaring episode in Ton 599, we see a possible hint of disc-jet coupling, which has been explored in detail.

In Section 2, we elaborate on the analysis procedure of various satellites. In Section 3, we show and discuss our results, followed by a summary in Section 4.

Table 1. Log of the observations during the flaring state (MJD 59436-59509).

2. Multiwavelength observations and data analysis

2.1. Fermi-LAT

Fermi-LAT is a $\gamma$ -ray telescope launched by NASA in 2008 into a near-Earth orbit. It works using the pair-conversion method to detect high-energy $\gamma$ -rays between the energy range of 20 MeV to a Few hundred GeV. When a gamma ray enters the telescope, it interacts with the material in the detector, and if the photon has sufficient energy, it may convert into an electron-positron pair through a process called pair production. A colorimeter at the base of the detector is placed to record the energies and the direction of these charged particles. It is sensitive to photon energies between 20 MeV to higher than 500 GeV and has a field of view of roughly 2.4 Sr (Atwood et al. 2009). The LAT’s field of view covers approximately 20% of the sky at any given moment, allowing it to scan the entire sky every three hours. Ton 599 has been continuously monitored by the Fermi-LAT since 2008. We analysed Fermi-LAT data from June 2020 to August 2024. The data is downloaded from the archive for a particular duration of time, and a circular region of 20 degrees is chosen around the source. We followed the standard analysis procedure of Fermi Tools and used Fermipy to derive the light curve and the spectrum. We followed the documentation of the Fermipy. ^a To account for the galactic and extragalactic background, we have used the background files (gll_iem_v07.fits, iso_P8R3_SOURCE_V3_v1.txt) provided by the Fermi support team. ^b

We produced the light curve for June 2020 to August 2024 and applied the Bayesian block method (Wagner et al. 2022) to identify the flares. Three major flares were identified, and the remaining time source was in a quiescent state. We named them Flare 1,2, and 3. Flare 1 was observed during MJD 59370-MJD 59430, Flare 2 was observed during MJD 59440-MJD 59500, and Flare 3 was observed during MJD 59920-MJD 60030. Maximum recorded flux during the Flare 1 was 2.31 $\times$ 10 $^{-6}$ ph cm $^{-2}$ s $^{-1}$ , during the Flare 2 was 3.06 $\times$ 10 $^{-6}$ ph cm $^{-2}$ s $^{-1}$ and during the Flare 3 was 3.63 $\times$ 10 $^{-6}$ ph cm $^{-2}$ s $^{-1}$ .

2.2. Swift-XRT/UVOT

Neil Gehrels Swift Observatory is a space-based observatory launched by NASA in 2004 to catch the brightest gamma-ray bursts (GRBs) in the sky. It consists of three main telescopes, namely: X-ray Telescope (XRT), Ultraviolet and Optical Telescope (UVOT), and Burst Alert Telescope (BAT). XRT and BAT are X-ray telescopes working in the energy range 0.3–10 (soft) keV and 3–150 (hard) keV, respectively. Most of the bright Fermi blazars detected in gamma rays have been monitored with the Swift. ^c The advantage of Swift is that it can provide simultaneous observations of any object in soft and hard X-ray, along with optical-UV. Particularly, the simultaneous observations across the wavebands are important to understand the physical mechanism producing occasional broadband spectacular flares in blazars. Blazar Ton 599 has been observed with Swift on multiple occasions during the gamma-ray flaring as well as non-flaring states. We have analysed all the swift observations available between May 2020 to September 2024. We present the outcome in this paper (see Table 1).

XRT: We first downloaded the raw data from the HEASARC Archives and ran the XRTPIPELINE to create the cleaned events files. Most of the observations are done with the photon counting (PC) mode, and hence, we took the cleaned event files of the PC mode to create the source and background regions. We chose the 20 arcsec radius around the source and away from the source for source and background regions and extracted the source and background spectra using the tool XSELECT. We then collected the redistribution matrix file (RMF) from the command ‘quzcif’ using the XRTMKARF; we then created the ancillary response file. All these files are loaded in GRPPHA to merge and then are used in XSPEC (Arnaud 1996) for modelling. We chose the simple power law to model the X-ray spectrum and also consider the galactic absorption ( $N_H$ =1.63 $\times$ 10 $^{20}$ cm $^{-2}$ ; HI4PI Collaboration et al. 2016) and estimated the spectral index and the flux in the 0.3–10.0 keV energy regime. The estimated photon index and the flux are quoted in Table 1, along with their error bars. To create an average spectrum for the SED modelling, we add all the observations together during that particular period using ADDSPEC and then model them in XSPEC to derive the SED data points.

UVOT: During Swift’s observation, the UVOT Telescope (Roming et al. 2005) observed the TON 599 in its three optical (U, B, and V) as well as in three ultraviolet(UVW1, UVM2, and UVW2) filters. We have downloaded the UVOT data from the HEASARC Data archive ^d and performed the data analysis. We begin by accessing the uvot/image within the specified obsid, where images for all wavebands are located. Opening the filename-sk.img.gz file in DS9, we use SIMBAD for object identification to generate the source and background region file. We then run uvotsource for each filter separately to obtain AB magnitude. The source magnitudes were obtained from a circular region with a 5 arcsec radius centered on the source, while the background magnitudes were measured from a nearby, source-free circular region with a radius of 20 arcsec. For SED modelling, we summed up all the observation IDs during the particular period using uvotimsum and extracted the source magnitude using uvotsource. Magnitude is corrected for galactic extinction (Schlafly & Finkbeiner 2011) and converted into flux using zero-points and conversion factors (Larionov et al. 2016).

2.3. NuSTAR

The Nuclear Spectroscopic Telescope Array (NuSTAR) mission, launched on 2012 June 13, is the first focusing high-energy X-ray telescope in orbit (Harrison et al. 2013). It consists of two co-aligned X-ray detectors paired with their corresponding focal plane modules, called FPMA and FPMB, and it operates over a wide energy range from 3 to 79 keV. It recorded two observations for Ton 599; we have analysed one observation (Obsid 60463037004) for 2021-06-25 (MJD 59390.53). We have downloaded the NuSTAR data from the HEASARC Data archive and run the nupipeline, so two cleaned event files were produced. We extracted the source and background spectra using a circular region and used the tool nuproducts. We grouped the spectra using grppha with 30 counts per bin and then used XSPEC for modelling.

Figure 1. The $\gamma$ -ray light curve using 1-day binning. Multiple flares have been identified by applying the Bayesian block method.

2.4. XMM-Newton

In gamma-ray, the source was reported to be flaring during June 2021 through a telegram (ATel#14722; Principe 2021). Following this event, we proposed a target-of-opportunity observation in XMM-Newton to study the short-term variability in X-rays. Finally, the observation was done on 2021-06-25 at 18:52:30. We also looked at the archive and collected the observations done during a low flux state. The timeline of the observations is marked in Fig. 1.

We followed the standard analysis procedure to analyse the XMM-Newton data using the Science Analysis Software (SAS) version 18.0.0. We extracted both the light curves and the spectra by selecting a source region of 20 arcsecs around the source and a background region of a circular radius of 50 arcsecs away from the source. The details of the observations are shown in Table 2.

Table 2. The observational log of XMM-Newton and NuSTAR, which are used in this work.

2.5. ZTF

Zwicky Transient Facility (ZTF) is a new optical time-domain survey that uses the Palomar 48-inch Schmidt telescope. ZTF uses a 47 deg $^2$ field with a 600-megapixel camera to scan the entire northern visible sky at rates of $\sim$ $3\,760$ deg $^2$ per hour (Bellm et al. 2018). ZTF has three filters: ZTF-g, ZTF-i, and ZTF-r. We have downloaded ZTF observations from the NASA/IPAC Infrared Science Archive. ^e We get the magnitude(AB) and time (MJD), convert the magnitude into flux using zero-points, and take the center wavelength $\lambda_{\text{cen}}$ from the SVO filter service. ^f

2.6. NEOWISE

The Wide Field Infrared Explorer(WISE) has surveyed the entire sky in four infrared wavelengths with high sensitivity and spatial resolution devices. NEOWISE(Near-Earth Object Wide-Field Infrared Survey Explorer) is the enhancement of WISE. The NEOWISE mission continues to detect, track, and characterise minor planets (Mainzer et al. 2011). The four infrared filters are W1, W2, W3, and W4; we have taken only W1 and W2 in our study. We downloaded the data from the Infrared Science Archive and changed AB magnitude to flux using center wavelength $\lambda_{\text{cen}}$ from the SVO ^g filter service.

2.7. ASAS-SN

All-Sky Automated Survey for Supernovae (ASAS-SN), currently consisting of 24 telescopes around the globe, automatically surveys the entire visible sky every night down to about 18th magnitude in both V and g filters. We have accessed the data using ASAS-SN Sky Patrol V2.0 ^h (Hart et al. 2023; Shappee et al. 2014). We took only the g-filter data with good(G) quality, used it in the light curve, and took the average flux for the particular period to use it in SED modelling.

2.8. Archival

We have downloaded the archival data from SED Builder ⁱ and use it as background data in our broadband SED modelling. This helps us to guide the SED modelling and put better constraints on the total SED.

3. Results and discussions

We have analysed Fermi-LAT, Swift-XRT/UVOT data from 2019 January to August 2024 Jan (MJD 58500–MJD 60542). We have also collected archival data to study the broadband behaviour of the source.

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. (2010a)

(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 1979). 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. (2003).

(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. (2019). 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 2020) or the presence of steady thermal emission (Bonning et al. 2009).

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 2010; Schneider et al. 2010) 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. 2023). 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 2019).

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 2006), 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 2010). 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 2005; McHardy 2010). The log-normal distributions are very common in AGNs, Gamma-ray Bursts, and the galactic X-ray binaries (Uttley & McHardy 2001; Negoro & Mineshige 2002; Shah et al. 2018a; Khatoon et al. 2020; Prince, Khatoon, & Stalin 2021) 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. (2018a) & Prince et al. (2021) 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 2012). On the other hand, Scargle (2020) 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. (2023) and Khatoon et al. (2020).

(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. 2023), 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. 2018a; Kushwaha et al. 2016). 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. (2012) have shown that mostly FSRQs follows ‘redder-when-brighter’ whereas Ikejiri et al. (2011) 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. 2015).

(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 1998).

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. 2020).

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 2013; Sarkar et al. 2019).

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. 2014). 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. 2018). This also helps us to understand what exactly is happening inside the jets (Pushkarev, Kovalev, & Lister 2010; Cohen et al. 2014).

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 2013). A Fortran version of the code is available online. ^j The main reason for using zDCF over the discrete correlation function of Edelson & Krolik (1988) 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. 2010; Raiteri et al. 2013; Cohen et al. 2014; Sarkar et al. 2019).

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 (2019) during the flare of 2018 and in the recent flare by Rajput et al. (2023).

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. 2019). 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. (2003).

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 2011; Sobolewska et al. 2014; Finke & Becker 2014; Chen et al. 2016; Kushwaha et al. 2017; Chatterjee et al. 2018a; Ryan et al. 2019; Bhattacharyya et al. 2020).

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. 2010a; Isobe et al. 2014; Abdollahi et al. 2017) and in lower energy (radio and optical) by damped/red-noise type processes ( $\beta=2$ ) (Max-Moerbeck et al. 2014; Nilsson et al. Nilsson, Lindfors, Takalo, Reinthal, Berdyugin, Sillanpä 2018). The $\beta$ $\sim$ 0 interpreted as an uncorrelated white-noise-type stochastic process (Press 1978). However, Press (1978) 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. (2003) 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. 2003). 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 (1977) & Blandford & Payne (1982) 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 2004). 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 2008; Chatterjee et al. 2018a) and comparing the time scale of PSD with the accretion disc time scale as done in Sharma et al. (2023). The accretion disc variability is expected to produce the log-normal flux distribution (Uttley & McHardy 2001), 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. 2023), 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. 2018a; Kushwaha et al. 2016). 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 ^k (Tramacere et al. 2009; Tramacere, Massaro, & Taylor 2011; Tramacere 2020). 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. (2023). 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 2020) and the luminosity is also fixed to the standard value, $L_{disc}=1\times10^{45}$ erg/s (Ghisellini et al. 2010). The viewing angle of the jet is fixed as 2 deg (Pushkarev et al. 2009) 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 2012). 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. (2009). 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. (2023), Manzoor et al. (2024) 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 2004; Liu, Zhao, & Wu 2006). 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 (2020) 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.

4. Summary

The Bayesian block methodology turned out to be the best way to characterise the light curve into various high and low flux states. We identified three bright gamma-ray flares during the years 2019–2024. The rise and decay times of flares are found to be in the range of 1–7 h. The F $_{\text{var}}$ shows a double hump structure resembling the broadband SED. The XMM-Newton spectra taken during the flaring episode reveal the presence of a thermal black-body component, suggesting the possible contribution of the accretion disc in the total jet emission. The gamma-ray spectral points show a clear break around 10 GeV, suggesting the emission region to be located within the BLR. A double-log normal flux distribution is found, which is rarely seen in blazars. The gamma-ray, optical, and X-ray emissions were found to be highly correlated with zero time lag, suggesting their co-spatial origin. A total of 16 yr of Fermi data are used to derive the PSD, which is best fitted with a single power law, suggesting that a long baseline is required to probe the characteristic time scale. The one-zone leptonic model is successfully used to fit the broadband data, and the finding is that the high flux state is most probably caused by the increase in the magnetic field and the rise in the Doppler factor. The results of XMM-Newton spectral modelling, along with the flux distribution, suggest there is a possible disc-jet coupling in this source.