© The Authors 2025

1. Introduction

Blazars are a class of active galactic nuclei (AGNs) whose powerful relativistic jets are oriented very close to our line of sight Urry & Padovani 1995. This alignment causes significant Doppler boosting of the emitted radiation, leading to extreme observational properties such as high luminosity, rapid and large-amplitude variability across the entire electromagnetic spectrum, and strong, variable polarization Urry & Padovani 1995; Kinman et al. 1966. The variability of blazars occurs on a wide range of timescales, from minutes (intraday variability, IDV) to years (long-term variability, LTV), and is generally stochastic in nature, with larger amplitudes on longer timescales (e.g., Chatterjee et al. 2008, 2009; Kushwaha et al. 2017).

Blazars are traditionally classified into two main categories based on their optical spectra: BL Lacertae objects (BL Lacs) that exhibit weak or no emission lines, and flat-spectrum radio quasars (FSRQs) that are characterized by strong, broad emission lines Blandford & Rees 1978; Stocke et al. 1991. Their spectral energy distribution (SED) typically shows two broad humps Fossati et al. 1998. The low-energy component, from radio to optical/UV frequencies, is dominated by synchrotron radiation from relativistic electrons in the jet Kinman et al. 1966. The high-energy component, from X-rays to γ-rays, is most commonly explained by inverse Compton (IC) scattering of lower-energy photons Marscher & Gear 1985; Sikora et al. 1994; Dermer et al. 1992; Ghisellini et al. 1996.

Based on the position of the synchrotron peak frequency of the SED, blazars are divided into three categories: low-synchrotron peaked (LSP) blazars (ν_peak < 10¹⁴ Hz), intermediate-synchrotron peaked (ISP) blazars (10¹⁴ < ν_peak < 10¹⁵ Hz), and high-synchrotron peaked (HSP) blazars (10¹⁵ < ν_peak < 10¹⁷ Hz) Ackermann et al. 2011. The FSRQs are mostly LSPs, while BL Lacs can be found in all three categories and are, accordingly, called low-, intermediate-, and high-frequency-peaked BL Lac objects Abdo et al. 2010a. There are blazars for which the synchrotron peak is at ν_peak > 10¹⁷ Hz; these are called extremely high-synchrotron-peaked BL Lacs (EHBL; Foffano et al. 2019).

The subject of this study, Ton 599 (also known as 1156+295), is a luminous FSRQ Healey et al. 2007 at a redshift z = 0.72469 Hewett & Wild 2010. It is powered by a supermassive black hole with a mass of approximately 9 × 10⁸ M_⊙ Keck 2019. VLBI observations of its parsec-scale jet have revealed highly superluminal motion, allowing for the determination of a Doppler factor of δ = 12 ± 3 and a jet viewing angle of θ ≤ 2.5° Jorstad et al. 2017. At γ-ray energies, the source was first detected by the Compton Gamma-Ray Observatory (CGRO) and listed in the second Energetic Gamma Ray Experiment Telescope (EGRET) catalog Thompson et al. 1995. Ton 599 was confirmed as a gamma-ray source by the Fermi Large Area Telescope (LAT), a detection made possible by the instrument’s improved sensitivity within the first three months of its mission Abdo et al. 2010b. Ton 599 has a rich history of dramatic activity. Early observations documented rapid flux changes, such as a 0.05 mag brightening in the B band in just 29 minutes Grauer 1984. A long-term study by Fan et al. 2006 found large variations of nearly 5 magnitudes and suggested periodicities of 1.58 and 3.55 years. More recently, a major multiwavelength outburst in late 2017 led to its first detection at very high energies by VERITAS and MAGIC Mirzoyan 2017; Mukherjee & VERITAS Collaboration 2017. Detailed studies of this event focused on constraining the jet’s physical parameters and the location of the emission region, revealing a complex connection between the nonthermal continuum and the broad emission lines Hallum et al. 2022; Patel et al. 2018; Prince 2019. The source’s activity continued, with a record-breaking γ-ray flare in January 2023, which has also been a subject of detailed broadband spectral modeling Patel & Chitnis 2020; Manzoor et al. 2024. Other recent works have investigated its long-term radio kinematics and magnetic field properties Ramakrishnan et al. 2014; Kang et al. 2021.

Given this history of complex and powerful outbursts, continuous and detailed long-term monitoring is essential to fully understand the physical processes at play. In this paper, we present one of the most extensive optical photometric studies of Ton 599 to date, utilizing a rich dataset collected from November 2011 to September 2023. The core of our analysis is based on BVRI photometry from the Whole Earth Blazar Telescope (WEBT) collaboration, which we complement with crucial long-term photometric and spectroscopic data from the Steward Observatory. Our primary goal is to perform a comprehensive characterization of the blazar’s multiband optical variability on diverse timescales in order to investigate the physical parameters and processes in the central part of this AGN. To this end, we test the flux distributions against normal and log-normal models and explore the root mean square (RMS)–flux relation. We quantify intranight changes using χ² tests and fractional variability to estimate the emitting region size and magnetic field strength. Finally, we trace the source’s complex color evolution using color–magnitude and color–time diagrams. Furthermore, we leverage the low-flux optical spectra to measure the properties of the Mg II emission line, after correcting for Fe II contamination, to provide an independent estimate of the central black hole mass via single-epoch scaling.

In Sect. 2, we describe the extensive optical photometric data collected from various observatories worldwide under the WEBT consortium. In Sect. 3, we present our results related to flux, color, and spectral studies. In Sect. 4, we discuss our findings and in Sect. 5 we summarize our conclusions.

2. Optical data

For the analysis of the emission behavior of the blazar Ton 599, we used optical photometric data obtained within the framework of the Whole Earth Blazar Telescope (WEBT)¹ project (e.g., Villata et al. 2002, 2008; Raiteri et al. 2017; Larionov et al. 2020; Jorstad et al. 2022). The analysis covers November 2011 to September 2023 (from JD 2455875.6141 to JD 2460192.30862). Photometric observations were made in B, V, R, and I Johnson-Cousins filters. The total number of data points in the above filter bands included in the final processed light curves are 1489, 2089, 6625, and 1851, respectively. Table A.1 lists the observatories, the telescope diameters that were used for the observations, and the total number of data points provided by individual instruments included in theanalysis.

Instrumental magnitudes were reduced using master bias, dark (where appropriate) and flat-field calibration images. Cosmetic corrections included masking hot and dead pixels, removal of cosmic rays and fringing in the I filter where fringing is most pronounced. Photometric calibration was done with standard stars listed on the official WEBT website².

The optical light curves were carefully processed to ensure homogeneity and precision across all four BVRI bands. The following steps were applied: (i) removal of outliers, including those within intranight sequences; (ii) binning of closely spaced and noisy measurements from the same dataset; (iii) inspection of data points with large uncertainties – when supported by independent observations taken on the same or nearby nights, the uncertainties were revised accordingly, otherwise the data points were discarded; (iv) imposition of a minimum photometric uncertainty of 0.01 mag and exclusion of all measurements with uncertainties larger than 0.15 mag; and (v) correction for systematic offsets between different datasets, applied when necessary to minimize inter-dataset discrepancies, while avoiding artificial time-dependent shifts.

Throughout the paper we use magnitudes corrected for Galactic extinction. Galactic extinction values were taken from the NASA/IPAC Extragalactic Database³, with values of 0.072, 0.054, 0.043, and 0.030 for the B, V, R, and I bandpasses, respectively. For the conversion of magnitudes into fluxes, we use zero-magnitude flux densities provided by Bessell et al. 1998. Throughout the paper, times are expressed as JD–2455000. Assuming a flat Universe (H₀ = 68 km s⁻¹ Mpc⁻¹, Ω_m = 0.3 and Ω_Λ = 0.7), the luminosity distance at redshift z = 0.72469 is D_L ≈ 4.6 Gpc.

3. Analysis and results

Figure 1, viewed from top to bottom, shows the light curves of Ton 599 blazar in the B, V, R, and I filters, respectively, for the whole time period under consideration. Data points from different observatories and instruments are marked with different colors and symbols, as described in Table A.1.

Fig. 1. Flux density evolution of the blazar Ton 599 from November 2011 to September 2023 in B, V, R, and I filters. Colors and symbols in the plots are explained in Table A.1.

3.1. Global flux properties

3.1.1. Flux distribution.

The distribution of the blazar radiation flux can provide clues about the underlying nature of the radiation process. A normal flux distribution typically suggests additive processes, where the total radiation arises from the sum of contributions from multiple independent sources (e.g., independent emission regions or individual shocks within the blazar). In contrast, a log-normal flux distribution implies multiplicative processes, where the radiation is shaped by a chain of factors that multiplicatively influence the total flux (e.g., variability driven by cascading turbulent processes, relativistic beaming, or interactions in a jet). In some blazars, neither of these distributions is suitable to describe the flux distribution (see, e.g., Shah et al. 2018, where the flux distribution in the γ domain suggests the existence of two log-normal distributions).

We performed a flux distribution analysis for the entire light curve across all filters. As an example, the left panel of Figure 2 displays the flux distribution in the R band⁴. Prior to creating the histogram, the data were averaged into one-day bins to reduce noise. The histogram was then fit with both normal and log-normal functions, represented by the smooth red and blue lines, respectively. The reduced χ² values for the normal and log-normal fits indicate that the flux distribution is better described by the log-normal function in all filters. Specifically, the reduced χ² values for the normal (log-normal) fits to the distributions for the B, V, R, and I bands are 24.43 (4.54), 20.62 (2.85), 29.76 (3.10), and 11.01 (2.13), respectively.

Fig. 2. Left panel: Flux distribution from ∼12 years of R-band observations. The histogram (blue bars) shows the flux counts. The red curve is the best-fitting normal function, and the blue curve is the best-fitting log-normal function. Right panel: RMS (i.e., the square root of the excess variance) as a function of the mean flux, computed in (1) bins containing an equal number of data points (blue markers), and (2) bins of fixed time length (red markers). The solid blue and red lines indicate the best linear fits for each binning scheme.

The results of the histogram-based analysis depend on the choice of bin size. To address this, we performed fitting on the unbinned flux values using the maximum likelihood estimation method. To evaluate the goodness-of-fit, we utilized the Bayesian information criterion (BIC) and Akaike information criterion (AIC) statistical measures. The lower values of both BIC and AIC for the log-normal function compared to the normal function further confirm that the flux distribution is better described by a log-normal model.

In addition to the previous analysis, we also performed the Kolmogorov-Smirnov (KS) test to further assess the goodness-of-fit for the flux distributions. The KS test compares the empirical cumulative distribution of the data with theoretical distributions (normal and log-normal), providing two key statistics: D, the maximum deviation between the observed and expected cumulative distributions, and p, the p-value, which indicates the significance of the fit. The results of the KS test are given in Table 1. These results reinforce the conclusion that the flux distribution of the blazar is better described by a log-normal distribution in all observed filters.

Finally, we give the parameters that best fit a log-normal distribution of the form

$$\begin{array}{c}\hfill f(x)=\frac{1}{xs\sqrt{2\pi }}exp[-\frac{{(ln(x)-m)}^2}{2s^2}],\end{array}$$

where m and s are the mean location (in the natural logarithm space) and shape parameter (standard deviation in the natural logarithm space): B: m = −0.05 ± 0.10, s = 1.24 ± 0.10

V: m = 0.18 ± 0.06, s = 1.17 ± 0.06

R: m = 0.20 ± 0.04, s = 1.24 ± 0.04

I: m = 0.81 ± 0.07, s = 1.17 ± 0.08.

3.1.2. RMS-flux relation.

To further characterize the flux distribution and variation, we analyzed the dependence of the degree of flux variation on the mean radiation flux (known as the RMS-flux relation). We binned the data in two ways to calculate both the degree of variation and the mean flux within each bin: (a) using fixed time intervals of 50 days, and (b) ensuring that each bin contained exactly 20 data points. To quantify the amplitude of the flux variation in each bin, we employed the excess variance ${\sigma }_{\mathit{XS}}=\sqrt{S^2-⟨{\sigma }_{\mathrm{err}}^2⟩}$, where S² is the variance of flux in the light curve and ⟨σ_err²⟩ the average of squared measurement error.

The right panel of Figure 2 shows the RMS-flux relation for the two binning modes for the R-band light curve. Red dots correspond to binning with a fixed number of data, and blue dots correspond to binning with a fixed time interval. We fit a straight line using the least-squares method; the fit is shown in the figure as a straight line with the same color as the corresponding data. The analysis shows a significant positive correlation for all filters, and the results are presented in Table 2.

3.1.3. PSD distribution.

The power spectral density (PSD) of blazars in the optical domain often follows a power law, P(ν)∝ν^−β, with typical slopes of β ∼ 1 − 3. This indicates that flux variability of blazars is a correlated colored-noise stochastic process, with the amplitude of flux variations decreasing toward shorter timescales (specifically, β ∼ 0 is a uncorrelated white noise process, β ∼ 1 is a correlated pink (or flickering) noise process, β ∼ 2 is a red (or random-walk) noise process, β ≥ 3 is a black noise process). Determining the slope and normalization of the PSD, as well as the inflection points, is of great importance because they carry important information about the system in which variations occur, such as the size of the emission region or the timescale of particle cooling in the emission region.

We adopted the method described by Vaughan 2005 to compute the periodogram of our light curves. Concretely, if x_j (j = 1, …, N) denotes the time series of flux values, then we define the periodogram, P(ν_k), as

$$\begin{array}{c}\hfill P({\nu }_k)=\frac{2T}{{⟨x⟩}^2{N}^2}|\mathrm{DFT}({\nu }_k){|}^2,\end{array}$$

where DFT(ν_k) is the discrete Fourier transform of the light curve, $T\equiv \frac{N(t_N-t_1)}{N-1}$ is the effective duration of the time series, and ⟨x⟩ is the mean flux. We computed the DFT at frequencies ν_k = k/T (k = 1, …, N/2), so that ν_min = 1/T and ν_Ny = 2/(N T) is the Nyquist frequency. By integrating the positive part of this periodogram, one recovers the fractional variability (rms/⟨x⟩).

Because our data are not uniformly sampled, we interpolated the data onto a uniform grid with a time step of T/N. Next, we subtracted the mean flux (so that the time series has zero mean) and then applied a Hanning window to reduce leakage effects (see Max-Moerbeck et al. 2014 for the importance of centering the light curve and using a Hanning window before calculating the periodogram). The resulting periodogram, log P(ν_k) versus log ν_k, for the R bandpass is shown in Fig. 3. The vertical dashed line corresponds to the Nyquist frequency. The horizontal dashed line corresponds to the power that comes from the uncertainty of measurement, calculated as $P_{\mathrm{noise}}=(2T)/(N{⟨x⟩}^2)({\displaystyle \sum _{j=1}^N}{\sigma }^2/N)$, where σ is the flux error.

Fig. 3. R-band periodogram of Ton 599.

Following the method of Vaughan 2005, we fit the raw (unbinned) periodogram with a linear function in log–log space using the least-squares method, obtaining both the best-fit parameters and their uncertainties, as well as model-based error estimates. In Figure 3, the red line represents the best-fitting power-law model, while the shaded region indicates the range of model errors around that fit. For the R band, we find a power law slope of −1.40 ± 0.05. We also performed a KS test on the ratio of the observed periodogram to the model, confirming that it follows a χ² distribution with two degrees of freedom (χ₂²). The KS statistic (and corresponding p value) is 2.37 × 10⁻² (7.12 × 10⁻¹).

We calculated the 99.73% confidence limit by evaluating the quantity

$$\begin{array}{c}\hfill {\gamma }_{ϵ}=-2ln[1-(1-{ϵ}_{n^{\prime }}{)}^{1/n^{\prime }}],\end{array}$$

where ϵ represents the false-alarm probability, and n′=n − 1 indicates that the Nyquist frequency is excluded from the calculation. We then added log(γ_ϵ/2) to the periodogram to define the threshold level. The 99.73% confidence limit is drawn as dashed line above the periodogram, and it shows that there is no statistically significant peak indicating periodicity in our light curve. We obtained a similar result for the light curves for the B, V, and I bands. Here we list only the slopes of the obtained periodograms: −1.79 ± 0.07, −1.53 ± 0.06 and −1.69 ± 0.07, respectively.

3.2. Intranight and intraday flux variations

In order to study the shortest variations in brightness of the blazar Ton 599, we extracted the datasets that were taken from the same observatory within a single observing night. For more precise statistical analysis, we considered only datasets with more than 2 hours of observation and more than 20 data points. With these criteria, 31 light curves were extracted.

We used a χ² test to determine whether the light curves exhibit statistically significant variability on a given night. The χ² statistic was computed as

$$\begin{array}{c}\hfill {\chi }^2={\displaystyle \sum _{i=1}^N}\frac{{(O_i-E_i)}^2}{{\sigma }_i^2},\end{array}$$(1)

where O_i is the observed flux, E_i is the expected value (here, the mean flux), σ_i is the measurement uncertainty for the ith observation, and N is the total number of data points. We then compared the resulting χ² value to a χ² distribution with N − 1 degrees of freedom to obtain a p-value (p), which reflects the likelihood that the observed variations are solely due to measurement errors. We adopted p < 0.05 as the threshold to reject the null hypothesis of a constant source flux. Consequently, p < 0.05 indicates significant variability, whereas p ≥ 0.05 suggests that any apparent flux changes could be attributed to random noise rather than genuine variability. As a result of the test, we found statistically significant variability in 14 intranight light curves across multiple filters, corresponding to 8 distinct nights of observation..

We adopted the fractional variability measure from Peterson 2001 to characterize the amplitude of flux variations in the intranight light curves. It is defined by

$$\begin{array}{c}\hfill F_{\mathrm{var}}=\frac{\sqrt{{\sigma }^2-{\delta }^2}}{⟨F⟩},\end{array}$$(2)

where σ² is the sample variance of the fluxes, δ² is the mean square measurement error, and ⟨F⟩ is the mean flux. Errors for F_var were calculated as follows (e.g., Bhatta & Dhital 2020):

$$\begin{array}{c}\hfill {\sigma }_{F_{\mathrm{var}}}=\sqrt{F_{\mathrm{var}}^2+\sqrt{\frac{2}{N}\frac{{⟨{\sigma }_{\mathrm{err}}^2⟩}^2}{{⟨F⟩}^4}+\frac{4}{N}\frac{⟨{\sigma }_{\mathrm{err}}^2⟩}{{⟨F⟩}^2}{F}_{\mathrm{var}}^2}}-F_{\mathrm{var}}.\end{array}$$(3)

This definition accounts for measurement errors, thereby enabling the straightforward identification of light curves in which measurement uncertainties dominate the sample variance.

Table 3 shows the result of the χ² test (column 4) and the F_var measurement (column 5) for observation nights when the blazar was variable according to the χ² test. Considering the observation dates (column 2) and filters (column 3), the source was variable on 8 nights, with a total of 14 intranight light curves showing significant variability.

Figure 4 shows the intranight light curves for which the χ² test showed to be variable. In the reading direction, the subplots are arranged in the same order as in the Table 3. The timescale is expressed in hours relative to the starting time of observation for a given night. In the legends of each subplot, we give the filter and the name of the observatory.

Fig. 4. Intranight light curves displaying statistically significant variability based on the χ2 test. The vertical axis shows the measured magnitude, the legend indicates the filter used for each observation, and the horizontal axis represents the relative time in hours.

For each light curve listed in Table 3, we calculated the variability timescale following the prescription of Burbidge et al. 1974:

$$\begin{array}{c}\hfill {\tau }_{\mathit{ij}}=\frac{\mathrm{\Delta }t}{ln(F_i/F_j)},\end{array}$$

where Δt is the time difference between two flux measurements F_i and F_j. We calculated τ_ij for each pair (i, j) of measurements that satisfies |F_i − F_j|> (σ_Fi + σ_Fj). The characteristic timescale of variability is then taken to be the minimum of all computed τ_ij values. We derived the corresponding uncertainty via standard error propagation from the flux measurement errors. The timescales of variability for each epoch are given in column 6 of Table 3.

The upper limit on the size of the emission region can be estimated using the relation

$$\begin{array}{c}\hfill R\le c{\tau }_{\min }\frac{\delta }{1+z},\end{array}$$

where δ is the Doppler boosting factor. For our calculations, we adopted a mean Doppler factor of δ = 23.1, derived from multiple estimates reported in the literature. Specifically, Jorstad et al. 2017 obtained δ = 11.8 based on Very Long Baseline Array (VLBA) observations, utilizing measurements of the apparent speeds of moving knots and variability timescales from light curves. Savolainen et al. 2010 determined δ = 28.2 using a combination of VLBA and Metsähovi data. Chen 2018 estimated δ = 29.3 from broadband spectral energy distribution modeling. The adopted value represents an average of these independent estimates, providing a robust constraint on the size of the emission region.

To estimate the magnetic field strength, we applied the condition that the observed minimum variability timescale, derived from intranight observations, must be at least as long as the synchrotron cooling timescale Hagen-Thorn et al. 2008:

$$\begin{array}{c}\hfill t_{\mathrm{syn}}\approx 4.75\times {10}^2(\frac{1+z}{\delta \nu B^3}{)}^{0.5},\end{array}$$

where ν is the observation frequency in Hz, t_syn is in days, and B is in Gauss. The results of the analysis for the size of the emission region and B are given in columns 7 and 8 of Table 3.

Otero-Santos et al. 2024 demonstrated that the flux distributions derived from TESS observations on IDV timescales are well fit by a Gaussian model. They argue that the success of the Gaussian representation suggests that the rapid variability is generated in a compact region of the jet through linear (additive) processes rather than by disk-related mechanisms. In our study, we applied two normality tests – the Shapiro–Wilk Shapiro & Wilk 1965 and Anderson–Darling Anderson & Darling 1952 tests – to the light curves presented in Table 3. Our results indicate that the IDV flux distributions in all cases are consistent with a Gaussian profile.

3.3. Long-term variability

As can be seen from Figure 1, over a period of 12 years, the Ton 599 blazar goes through alternating phases of weak and strong activity. The average flux density for the entire time series in the BVRI filters is 2.50, 3.16, 3.31, and 5.73 mJy, respectively. For the same order of filters, the flux density range (minimum–maximum) is 0.11−9.49, 0.13−12.75, 0.11−23.25, and 0.14−20.39 mJy, respectively. We note here that, due to the difference in sampling in the BVRI filters, the maximum flux peaks do not coincide in time. For instance, in the R filter, the maximum flux is recorded around JD–2455000 = 4958.5, which is not at all covered by the observations in the other filters.

The time interval from November 2011 to September 2017, i.e., from JD–2455000 = 900 to 3000, is marked by several average-intensity flares (Period I). The maximum fluxes in this period are 3.82, 5.42, 6.49, and 10.10 mJy for the BVRI filter, which, due to different sampling, do not occur in the same epoch.

This period is followed by two large outbursts that occur in about 400 days from about JD–2455000 = 3000 to 3400 (Period II). Their structure is much more complex than it can be seen in Figure 1 and consists of a series of smaller-larger flares (see Figure 5). The peak of the first outburst, which in all filters occurs nearly simultaneously around JD–2455000 = 3098 (see Table 4), reaches the value of 9.48, 12.75, 14.71, and 19.89 mJy for the BVRI filters, respectively (see the inset plot of the Figure 5). According to the R-band light curve, this peak is the second largest for the entire observation period. Assuming that the spectrum can be described with a power law of the form F_ν ∼ ν^−α, where α is the spectral index, we calculated α for the 3098 peak using the relation

$$\begin{array}{c}\hfill \alpha =\frac{log(F_1/F_2)}{log({\nu }_1/{\nu }_2)},\end{array}$$(4)

Fig. 5. Zoomed-in R-band light curve from JD–2455000 = 3050 to 3300 (Period II) and from JD–2455000 = 4150 to 5120 (Period III) in the R band. Complex structure of the outbursts is visible. The inset highlights the second-largest peak, well sampled in all filters.

where the subscripts 1 and 2 denote any combination of filters. By averaging spectral indices obtained from all possible combinations of filters (BV, BR, BI, VR, VI, RI), the mean value is α = 1.2 ± 0.2, where the error is given as the standard deviation of individual α values. After correction to the cosmological redshift, using equation F_corr = F(1 + z)^3 + α Nilsson et al. 2018, the logarithms of luminosity for the BVRI filters are 48.21, 48.24, 48.23, and 48.27 erg/s, respectively, which is on average for the optical equal to log(νL_ν) = 48.24 ± 0.02 [erg/s]. It is worth noting that the maximum flux density for the entire time series in the R band is 23.25 mJy, which is about 1.6 times greater than this peak value.

The peak of the second outburst in period II, which also occurs almost simultaneously (within one day around JD–2455000 = 3248), is 5.09, 7.07, 9.59, and 14.12 mJy for BVRI, respectively. With the same procedure as for the first peak, for this peak we find that the logarithm of the luminosity is 48.06, 48.11, 48.17, and 48.25 erg/s in the BVRI filters; that is, in the average log(νL_ν) = 48.15 ± 0.07 [erg/s].

After these two outbursts, there follows a period of about 2.2 years of relatively weak activity, marked by only two smaller flares (maximum flux in the R filter about 2.2 mJy). This phase is interrupted by an active phase lasting about 2.8 years (from about JD–2455000 = 4200 to the end of the time series; Period III) which is marked by a series of large flares. This period also saw the biggest peak in the blazar light curve, with flux density 23.25 mJy in the R band, which is a record for the period we are analyzing. As we mentioned earlier, the maximum flux of 23.25 mJy in the R band was not observed in other filters, so it is not possible to determine the luminosity as for the peaks in the period II. However, we can make an estimate using the average spectral indices that we calculate in Sect. 3.4. The mean spectral indices for the B-V, B-R, B-I, V-R, V-I, and R-I color indices (see Figure 8) are 1.38, 1.30, 1.43, 1.21, 1.45, and 1.64, respectively, from which it follows that the mean luminosity for the highest peak in the R band is log(νL_ν) = 48.48 ± 0.03 [erg/s], where the error is calculated as the standard deviation of the luminosity values calculated with different mean α.

Table 4 summarizes some of the variability parameters mentioned above. The table is divided into four large groups, which present the parameters for the entire time series, period I, period II, and period III, respectively, as indicated in their headers. Each of these groups list four parameters: the JD–2455000 of the maximum flux density, the average flux density in mJy, the maximum flux density in mJy, and the amplitude of variation calculated using Equation (2) in percentage.

Figure A.1 provides an overview of the entire dataset, divided into 12 time segments, each highlighted in light blue in the top panel. The lower panels offer a closer view of each segment, arranged sequentially from left to right.

3.4. Color analysis

Blazars show specific spectral characteristics. BL Lac blazars often show a bluer-when-brighter (BWB) long-term trend. Some authors found that long-term flux variations are quasi-achromatic and likely due to changes in the orientation of the jet, while the strongly chromatic short-term changes are attributed to energetic processes taking place in the optical jet (magnetic reconnection, shock) (e.g., Raiteri et al. 2021a,b, 2023).

Spectral behavior of the FSRQ type of blazars is somewhat different. In the case of CTA102, Raiteri et al. 2017 found RWB trend in the faint state, switching to a mild BWB trend above R ∼ 15 mag. A similar behavior was described by Villata et al. 2006 for the blazar 3C 454.3 and was interpreted as the influence of radiation from the disk in the faint state, which becomes less pronounced with the increase in blazar luminosity, where synchrotron radiation is dominant and the blazar exhibits BWB trend.

To study the spectral behavior of the blazar Ton 599, we used color indices and optical spectra. Color indices were calculated by combining data from the same observatory that were acquired within a time interval of 15 minutes. In the analysis we used only magnitudes with measurement uncertainties smaller than 0.03 mag. The number of color indices, which were obtained based on this criteria, ranges from ∼500 to ∼800 and is listed in the second column of Table 5. Magnitudes in these 15-min bins were calculated as weighted average of all magnitudes in the corresponding bins. For weights we used w = 1/magErr², where magErr denotes the magnitude uncertainty of each data point. Magnitude errors in bins were calculated according to the relation $\sqrt{1/{\displaystyle \sum {\mathrm{magErr}}^2}}$. To obtain errors for the color indices, we added the magnitude errors in bins in quadrature.

Figure 6 shows the change in various color indices over time. The names of the color indices are shown in the legends of the subplots. For comparison, the first subplot from the top shows the light curve in the R-band. The red line is the best linear fit to the data obtained by the least-squares method. The results of the best linear fit, along with the fitting errors for each parameter, are given in Table 5. The table also lists the Pearson’s linear correlation coefficient and the p value that tests the null hypothesis that the slope of the linear function is zero. These results suggest that there is no significant long-term color evolution for the blazar Ton 599.

Fig. 6. Evolution of the color indices with time. The first panel from the top shows the light curve in the R filter. The rest of the panels present the CI evolution. The red line is the best fit of a linear function to the data. Parameters of the linear function are provided in the Table 5.

The variation in the CI on smaller timescales is much more pronounced (see Figure 6). Furthermore, it is striking that the amplitude of the color variation decreases with average flux density. To test this claim, we divided the data into 12 segments of about one year using the light curve in the R band. As a measure of the color variation, we use the standard deviation of the colors in the segments and correlate this parameter with the corresponding mean flux density. The results of the correlation are shown in the Figure 7 for B-V, B-R, and V-R color indices. The upper subplot shows the R-band light curve, where the shaded parts of the plot indicate the 12 segments in which we perform the correlation analysis. The red lines are the best linear fit to the data obtained by the least square method. The Pearson’s correlation coefficients are −0.63, −0.59, and −0.58 for the B-V, B-R, and V-R color indices, respectively, as is indicated in the subplot legends. The p values, which tests the null hypothesis that the slope of the regression line is zero, are 0.03, 0.04, and 0.05 for the same order of the color indices. This result indicates that there is a significant negative correlation between the colors variations and the blazar brightness.

Fig. 7. Correlation between variation in the color indices and the mean flux density in 12 segments (shaded in the top R-band light curve). Red lines show the best linear fits.

Color-magnitude plots for various color indices are shown in the Figure 8. The dashed blue lines on the subplots indicate the mean values of the color indices and they are 0.45, 0.84, 1.49, 0.40, 1.04, and 0.65 for the B-V, B-R, B-I, V-R, V-I, and R-I colors, respectively. The right side of the y axis shows the spectral index, which was calculated from the Equation (4)

$$\begin{array}{c}\hfill \alpha =\frac{0.4CI-K}{log({\nu }_1/{\nu }_2)},\end{array}$$(5)

where subscripts 1 and 2 correspond to filter bands, K is a constant calculated from zero flux densities ZP₁ and ZP₁ according to relation K = log(ZF₁/ZF₂), and CI is the appropriate color index. The mean values of spectral slopes of the B-V, B-R, B-I, V-R, V-I, and R-I color indices are 1.38, 1.30, 1.43, 1.21, 1.45, and 1.64, respectively. The mean values of the color indices and α are summarized in the third and fourth columns of Table 5

Figure 8 shows the dependence of the color indices and spectral slopes on magnitude. As can be seen, in 12 years of observation, the magnitude of the blazar Ton 599 has changed by about 5 mag in all filters. In the phase of weak blazar brightness, we can see a RWB trend for all color indices except for R-I. For these color indices, RWB switches to an achromatic trend above a certain brightness. On the other hand, color indices that include the I filter (B-I, V-I, and R-I), show transition from achromatic to BWB trend above a certain magnitude level. In order to quantitatively describe these trends, we fit linear curves to certain segments of the color-magnitude plots showing RWB, achromatic and BWB trends. We determined the boundaries that describe the transition from one trend to another visually. The results of the linear fitting are presented in the Table 6. The table is divided into three groups corresponding to the RWB, achromatic, and BWB trends, which can be clearly discerned from the Figure 8. For each of these groups, we have three columns that describe the Pearson correlation coefficient, the p value that tests the null hypothesis that there is no linear relationship between the two parameters we are correlating (p < 0.05 rejects the null hypothesis), and the magnitude at which one trend changes to another. The results shown in the table suggest that the RWB, achromatic, and BWB trend are statistically significant.

Fig. 8. Dependence of the color indices and spectral slopes on brightness. The dashed blue line corresponds to the mean value of the color indices. The right y axis shows the spectral slope, α, which was calculated under the assumption that the optical flux can be described by the relation Fν ∼ ν−α (see the text).

In order to better understand the color-magnitude trends, we show several spectra of the blazar Ton 599 at different flux levels in the Figure 9. The spectra were taken from the publicly available archive of the Steward Observatory⁵, which were acquired as part of a monitoring program designed to support gamma-ray observations with the Fermi Gamma-ray Space Telescope Smith et al. 2009. Unfortunately, the Steward spectra were recorded from 2011 to 2018, so they do not cover Period III in which the most powerful outbursts occurred, but they cover the Period II well. For illustration, we selected representative spectra at different flux levels. In the observation frame, the spectra range from ∼4000 Å to ∼7600 Å, which corresponds to ∼2300 Å to ∼4400 Å in the rest frame of the blazar. Therefore, the broad emission line, which is pronounced in the low activity phase, is Mg II λ2798 originating from the broad line region (BLR) of the AGN. In addition to the strong Mg II emission line, broad Fe-complex features are clearly visible, which also originate from the AGN (e.g., Hallum et al. 2022). Absorption lines in the infrared part of the spectrum are atmospheric lines⁶. In addition to the spectra, the figure shows the transmission curves of the Johnson-Cousin filters that were used for photometric observations and are marked with transparent blue, green, orange, and red colors for the B, V, R, and I filters, respectively.

Fig. 9. Steward optical spectra of the blazar Ton 599 at different activity states.

From the spectra shown in Figure 9, several things can be noted. The first is that the contribution of thermal radiation is significant in the phases of weak blazar brightness. As the mean flux increases, the thermal component weakens and synchrotron radiation becomes increasingly dominant. This spectral behavior can explain the RWB trend we observe in the color-magnitude plots. Then, we see a diverse behavior of the spectra for (approximately) the same mean flux level, which explains the large dispersion of the color indices for a certain magnitudes that can be seen in the color-magnitude plots. Colors involving the I band show a large dispersion and this makes the BWB trend suspect.

In order to test if the thermal radiation may be responsible for the RWB trend of the color-magnitude diagram, we use Steward spectra. In the first step, we fit the spectra with a power law function F = F_o(λ/λ_o)^α through the points representing the continuum, where F_o is the normalization parameter at λ_o = 3000 Å in the rest frame and α is the spectral index. Both F_o and α are free parameters. The obtained best fit represents the spectrum without Mg II and Fe emission lines; that is, the continuum. Then we do convolution of both spectra, the continuum and the original, with the BVRI filter transparency curves to determine the corresponding magnitudes. The top two panels of the Fig. 10 show two spectra from the Steward archive, one corresponding to the low flux state of the blazar (left), and the other to a relatively high flux state (right). Data points we used for fitting the power law function, i.e., the continuum windows, are marked in red colors. The two middle continuum windows (5226−5329 Å, 6105−6209 Å in the observed frame) are taken from Hallum et al. 2022. The remaining two continuum windows were not taken from the literature, but were defined with the aim of eliminating the emission and absorption lines in the B and I filters as best as possible. We have also marked the BVRI transparency curves with the colors indicated in the plot legends. The black line represents the best fit power law functions for these two spectra. The lower left panel shows the B-R color versus R magnitude plot. Observations are marked with transparent gray circles. Transparent red and green circles represent calculated colors/magnitudes; red ones correspond to filter convolution with the original Steward spectra, and green ones correspond to convolution with the continuum. For better visualization, with red and green lines we present the binned values of colors and magnitudes with bin widths of 1 magnitude. The lower right panel compares the calculated color indices. The diagonal red line is the equality line y = x. The points are color-coded so that one can see how color indices depend on the brightness of the blazar. Both lower panels indicate that the color indices tend to be redder when the contribution from the thermal component–likely associated with the accretion disk emission in the optical–UV range–is not included, and that this trend becomes more pronounced at fainter flux levels (i.e., higher magnitudes).

Fig. 10. Top: Steward spectra in low (left) and high (right) flux states, with continuum fits (black) and BVRI filter curves. Bottom left: (B-R) vs. R color–magnitude diagram from observed (gray) and model-based (red/green) data. Bottom right: comparison of observed and continuum-based color indices. The color bar encodes B magnitude.

Vanden Berk et al. 2001 built a composite quasar spectrum from over 2200 SDSS spectra (rest frame 800−8555 Å) and showed that, between ∼1300 Å and 5000 Å, the continuum follows F_ν ∝ ν^−0.44. We adopt this as the “disk” index α_disk = 0.44. Jorstad et al. 2013 analyzed the blazar 3C 454.3 and measured the continuum spectral index α_cont as a function of the V-band flux S​_V using simultaneous photometry (within ∼1 h). They showed that the blue (disk) component remains constant while the red (jet) component varies, fit a linear relation over the range 0 ≤ S​_V ≤ 4 mJy, and inverted it at α_disk = 0.44 to estimate the disk flux as S​_disk ≃ 0.85 ± 0.15 mJy in V (and 0.91 ± 0.16 mJy in R). We followed a similar procedure on our two-filter photometry using the BV, BR, and BI pairs. The magnitudes were first converted to flux densities using standard zero-flux values. Errors were propagated using the relation σ_F = 0.4ln(10)Fσ_m. For each filter pair, we then computed the spectral index as α = −ln(F₁/F₂)/ln(ν₁/ν₂). The resulting α values were binned in log₁₀F₂, and we performed a weighted linear fit restricted to log₁₀F₂ ≤ 0. Finally, we solved for F₂ at a fixed value of α_disk = 0.44. As an example, Figure 11 shows the binned spectral index α versus log₁₀F_V for the BV pair. Blue symbols are the bin means with 1σ errors; gray points are individual measurements. The dashed red line is the weighted linear fit for log₁₀F_V ≤ 0, and the vertical red line marks the derived log₁₀F_V at α_disk = 0.44 (i.e., F_disk = 0.152 mJy). Table 7 summarizes the results for the three B-anchored pairs.

Fig. 11. Binned spectral index α vs. log10FV for the BV data. Blue symbols are the bin averages with 1σ errors; gray points are individual measurements. The dashed red line shows the weighted linear fit for log10FV ≤ 0, and the vertical red line indicates log10FV at αdisk = 0.44 (i.e., Fdisk = 0.152 mJy).

All three estimates of F_disk (0.13−0.20 mJy) lie well below the minimum observed flux (≳1 mJy), confirming that the thermal disk (big blue bump) contributes only a few tenths of a millijansky, while the variable nonthermal jet dominates at higher brightness. The slight increase in F_disk from BV to BR is consistent with a redder jet spectrum varying on top of a nearly constant blue disk component, in agreement with the two-component model of Jorstad et al. 2013 and Raiteri et al. 2007.

Fig. 12. Example of the Mg II + UV Fe II fit in the 2650−3050 Å range for the Steward spectrum observed on February 7, 2011. The core and wing Gaussian components of the Mg II line are plotted with solid lines, while the UV Fe II template from Popović et al. 2019 is plotted with a dashed line. The observed data are represented by dots, and the red line indicates the best fit.

3.5. Estimation of the black hole mass

We determined the black hole mass of Ton 599 using the Steward spectra and the empirical relation that connects the black hole mass with the measured parameters of the Mg II emission line Kong et al. 2006:

$$\begin{array}{c}\hfill \frac{M_{\mathrm{BH}}}{M_{\odot }}=2.9\times {10}^6{\left(\frac{L_{\mathrm{MgII}}}{{10}^{42}\mathrm{erg}/\mathrm{s}}\right)}^{0.57\pm 0.12}{\left(\frac{{\mathrm{FWHM}}_{\mathrm{MgII}}}{{10}^3\mathrm{km}/\mathrm{s}}\right)}^2,\end{array}$$

where L_MgII is the line luminosity and FWHM_MgII is the full width at half maximum (FWHM) of the line.

To accurately measure the Mg II spectral parameters needed for M_BH calculation, it is of great importance to extract the profile of the Mg II line, which overlaps with numerous UV Fe II lines. The additional complexity of the Mg II line is that it is actually a doublet (λ2796 Å and λ2803 Å) with an intensity ratio dependent on the optical depth of the emitting region. Given that the doublet separation (∼8 Å) is much smaller than the total line width, Mg II can be treated as a single line, as suggested by Marziani et al. 2013. We performed a multicomponent fit in the 2650−3050 Å spectral range, as described in detail in Popović et al. 2019. After subtracting the continuum, the Mg II emission line was fitted with a two-Gaussian model, where one Gaussian fits the core and one the wings of the Mg II line. In this way, we were able to fit the asymmetry well in the Mg II line. The UV Fe II lines in 2650−3050 Å range were fit with the complex UV Fe II template described in Popović et al. 2019. In that model, the UV Fe II lines are divided into six groups of Fe II lines, which have different intensity parameters, while the widths and the shifts of all UV Fe II lines are taken to be the same. The total flux of the Mg II line, as well as the FWHM(Mg II) used for the M_BH estimation, were obtained using the total Mg II profile. The flux was calculated as the sum of the two Gaussian components from the best fit. The FWHM of the total Mg II profile was measured by normalizing the peak of the summed two-Gaussian fit to unity and determining its half-maximum width. Figure 12 shows an example of the Mg II + UV Fe II fit in the 2650−3050 Å range for the Steward spectrum observed on February 7, 2011.

To determine the black hole mass, we selected ten spectra with the lowest continuum flux, minimizing the jet contribution while maximizing the contribution from the broad-line region (BLR). The measured parameters are listed in Table 8, with associated standard errors from the fit. The uncertainties in the individual black hole mass estimates were obtained by propagating the measurement uncertainties through the calibration relation. For the final black hole mass and its uncertainty, we adopted the mean value from Table 8 and the mean propagated uncertainty, yielding log(M/M_⊙) = 8.77 ± 0.26. The obtained value is consistent with typical SMBH masses of 10⁶ − 10⁹ M_⊙, exceeding the 10⁸ M_⊙ typical of radio-loud AGNs Chiaberge & Marconi 2011. It is also comparable to the result by Keck 2019, which for the black hole mass of Ton 599 obtained ∼9 × 10⁸ M_⊙. The Eddington luminosity is given by (e.g., Wiita 1985)

$$\begin{array}{c}\hfill L_{\mathrm{E}}=1.3\times {10}^{38}\left(\frac{M_{\mathrm{BH}}}{M_{\odot }}\right)[\mathrm{erg}/\mathrm{s}],\end{array}$$

which for our black hole mass estimate gives ∼7.7 × 10⁴⁶ erg/s.

4. Discussion

4.1. Flux distribution and RMS-flux relation.

The statistical properties of long-term optical variability in blazars, such as the flux distribution and the RMS–flux relationship, offer profound insights into the underlying emission mechanisms and the physical processes at play within their relativistic jets. Our comprehensive analysis of Ton 599 over a 12-year period reveals characteristics that consistently point toward multiplicative processes governing its variability on these extended timescales.

The flux distribution of Ton 599, when examined over the entire observational period, is best described by a log-normal distribution across all BVRI bands, rather than a normal (Gaussian) distribution. This finding is a strong indicator of multiplicative processes, where variations are compounded. Such log-normal behavior in long-term optical blazar light curves is corroborated by studies such as Bhatta 2021, who analyzed a sample of 12 γ-ray-bright blazars with decade-long optical data and similarly found their flux distributions to be well-characterized by log-normal probability density functions (PDFs). This was also interpreted as evidence of multiplicative nonlinear processes driving the variability. The log-normal nature of blazar flux has also been noted in other bands, such as γ-rays (e.g., Shah et al. 2018; Bhatta & Dhital 2020).

Complementing the flux distribution analysis, the RMS-flux relationship in Ton 599 shows a significant positive linear trend in all optical bands (see Table 2 and Figure 2). This linear scaling, where the amplitude of variability increases with the mean flux, provides further robust evidence for dominant multiplicative processes. Again, this aligns with the findings of Bhatta 2021 for their sample of blazars, where linear RMS-flux relations were observed in the optical domain, supporting the multiplicative scenario. The combination of a log-normal flux distribution and a linear RMS-flux relation is widely considered a strong signature of multiplicative processes in accreting astrophysical systems (e.g., Uttley et al. 2005). Such behavior is thought to arise from various physical scenarios. Bhatta 2021 discusses several such mechanisms, including the general concept of multiplicative coupling of perturbations occurring within the accretion disk and/or the jet. More specifically, models invoking propagating fluctuations within the accretion disk, where uncorrelated variations in viscosity (the α-parameter) at different radii modulate the accretion rate as they propagate outward Lyubarskii 1997, can account for these observations. These disk-borne modulations can subsequently propagate into the blazar jet via a strong disk-jet connection and are then significantly amplified by relativistic beaming and projection effects. Alternatively, processes within the relativistic jet itself, such as the “jets-in-jets” scenario whereby Poynting-flux-dominated jets might generate numerous isotropically distributed mini-jets Giannios et al. 2009, have also been proposed. Such mini-jet emission has been shown to produce highly skewed flux distributions and adhere to the RMS-flux relation. Our findings for Ton 599 are consistent with these general frameworks proposed for blazar variability.

It is noteworthy that a single, geometric interpretation can provide a compelling explanation for both of these findings. As has been proposed by Raiteri et al. 2024, changes in the viewing angle to the jet naturally lead to the observed effects; a decrease in the viewing angle enhances the Doppler beaming, which not only increases the observed flux, contributing to the log-normal distribution, but also amplifies the intrinsic variability, thus producing the linear relationship between RMS and flux that we observe.

In our detailed examination of the RMS-flux relation for Ton 599, the intercept of the linear fit was treated as a free parameter, and its value exhibited some dependence on the specific filter and binning method employed (Table 2 in this work). Using the “size correlation” method, statistically significant positive intercepts on the RMS axis were observed for the V (0.21 ± 0.09), R (0.17 ± 0.07), and I (0.54 ± 0.21) bands. While a literal interpretation of nonzero RMS at zero flux is generally considered nonphysical, such positive intercepts can, as discussed in studies such as Bhatta 2021 (who references Gleissner et al. 2004), indicate the presence of an additional variability component that may become more prominent at very low flux states. Alternatively, it may suggest that the true relationship deviates from perfect linearity at the lowest flux levels probed. The B-band intercept (0.09 ± 0.11) with this binning method was consistent with zero. In contrast, using fixed 50-day time intervals for binning (“time correlation”), the intercepts for the B (−0.02 ± 0.06), V (0.06 ± 0.08), and R (0.00 ± 0.06) bands were consistent with zero, supporting a more direct proportionality between RMS and flux in these cases. The I band showed a marginally positive intercept (0.20 ± 0.14) with this latter method. These nuances in the intercept values might reflect the sensitivity of this parameter to the binning approach and the consequent sampling of low-flux states. However, the robustly determined positive slope of the RMS-flux relation across all filters and both binning methods consistently points to the multiplicative nature of the dominant LTV in Ton 599.

The characterization of LTV in Ton 599 contrasts with its behavior on shorter timescales. Our analysis of intranight light curves indicates that these rapid (hourly) fluctuations are well represented by a Gaussian distribution, suggesting that additive processes are likely dominant on these very short timescales. This finding is in agreement with other studies, such as Otero-Santos et al. 2024, who reported Gaussian flux distributions for IDV in TESS observations of blazars, linking them to linear (additive) processes in compact jet regions.

Bridging these extremes, the study by Pininti et al. 2023 using ∼25-day continuous TESS light curves provides valuable context for intermediate timescales. While their general finding for a sample of 29 blazars was that optical flux histograms on these day-to-week timescales are often consistent with normal (Gaussian) PDFs (frequently bimodal), their specific result for Ton 599 (Sector 22) indicated that a log-normal distribution provided a better fit than a normal one, based on AIC and BIC criteria (see Table 3 in Pininti et al. 2023). This persistence of log-normality for Ton 599, from ∼25-day segments Pininti et al. 2023 up to the multiyear periods analyzed in this work, is particularly compelling.

Thus, a hierarchical picture of variability emerges for Ton 599: additive processes appear to govern the most rapid, hourly fluctuations (IDV), while multiplicative processes become dominant as the timescale extends from weeks (Pininti et al. 2023, for Ton 599) to years (this work). This makes Ton 599 an interesting object where the signature of multiplicative processes is robustly observed even on relatively short, continuous observational segments, setting it apart from the average behavior seen in the TESS sample by Pininti et al. 2023 and highlighting the persistent nature of these mechanisms in this particular blazar. This complex, timescale-dependent nature of flux statistics underscores the interplay of different emission mechanisms or conditions within blazar jets.

It is important to note, however, that Scargle 2020 has cautioned that certain classes of stationary, linear, additive time-series models can also, under specific conditions (e.g., with a particular choice for the distribution of innovations), reproduce statistical properties such as log-normal flux distributions and linear or near-linear RMS-flux relationships. This suggests that while the combination of these properties observed in Ton 599 strongly points toward a multiplicative origin in line with prevailing interpretations (e.g., Uttley et al. 2005; Bhatta & Dhital 2020), definitive conclusions about the underlying physical processes (i.e., ruling out all possible linear, additive scenarios) solely based on these flux statistics should be made with awareness of these potential mathematical degeneracies. Therefore, the interpretation of multiplicative processes is further strengthened when considered alongside other observational evidence from blazars, such as spectral variability, polarization characteristics, and higher-order timing statistics, which may provide more direct constraints on the nonlinear and coupled nature of the emission mechanisms.

4.2. Power spectral density.

Our PSD analysis of the long-term optical light curves of Ton 599 reveals that its variability is characteristic of a red-noise process across all BVRI bands. The PSDs are described well by a single power-law model, P(ν)∝ν^−β, with derived slopes β ranging from approximately 1.40−1.79 (e.g., β_R = 1.40 ± 0.05). Such slopes, typically falling in the range 1 ≤ β ≤ 3 for blazar optical variability, indicate that the flux variations are strongly correlated over time, with more power concentrated on longer timescales. This is a hallmark of stochastic processes where past variations influence future ones.

The specific slopes we find for Ton 599 are broadly consistent with those reported for other blazars in the optical domain based on long-term monitoring. For instance, our R-band slope of β_R ≈ 1.40 is remarkably similar to the average optical PSD slope of β_R ≈ 1.42 reported by Nilsson et al. 2018 for a large sample of blazars. Other studies of long-term optical variability have also found comparable slopes; for example, Goyal et al. 2017 reported a long-term optical (R-band) PSD slope of β = 2.0 ± 0.1 for PKS 0735+178 spanning 23 years. Earlier work on PKS 0735+178 by Ciprini et al. 2007, using structure function analysis, also found PSD slopes between 1.5 and 2.0 for long-term optical data. Analyses of Kepler optical data for blazars yielded β ∼ 1.5 − 2.0 Edelson et al. 2013; Revalski et al. 2014. The PSDs in this work were derived by fitting the periodogram of linearly interpolated light curves, applying mean subtraction and Hanning windowing to minimize leakage, a methodological approach also adopted in other optical blazar studies confronted with unevenly sampled data (e.g., Goyal et al. 2017).

A key result from our PSD analysis is the absence of any statistically significant peaks indicative of quasi-periodic oscillations (QPOs), or distinct breaks signifying characteristic timescales, within the frequency range probed by our decade-long dataset. This lack of pronounced features reinforces the interpretation that the variability of Ton 599 is predominantly stochastic. This finding is consistent with many blazar studies; for instance, Pininti et al. 2023 also reported no significant PSD breaks or QPOs in their analysis of short (∼25-day) TESS light curve segments for their blazar sample. Similarly, Goyal et al. 2017 did not detect any low-frequency flattening or high-frequency cutoffs in their composite optical PSD for PKS 0735+178, which extended down to minute timescales. While some blazars do exhibit QPOs in their optical light curves on year-like timescales (e.g., Bhatta 2021) for sources such as S5 0716+714 and Mrk 421), Ton 599 does not show such persistent periodic features in our current, extensive dataset. We also note that earlier work by Fan et al. 2006 suggested potential periodicities of 1.58 and 3.55 years for Ton 599 from R-band data spanning 1974−2002. Our more recent and extensive dataset does not provide evidence supporting these specificperiodicities.

The derived PSD slopes for Ton 599 can also be considered in the context of its blazar classification and spectral characteristics. Ton 599 is a flat spectrum radio quasar (FSRQ), a class of blazar whose synchrotron peak frequency (${\nu }_{\mathrm{peak}}^{\mathrm{syn}}$) typically lies in the low-synchrotron peaked (LSP) range (i.e., < 10¹⁴ Hz). While this provides a general classification, recent detailed multiwavelength spectral energy distribution (SED) modeling of Ton 599 by Manzoor et al. 2024 during its exceptionally bright flare in January 2023 has revealed a more dynamic picture of its synchrotron peak. Analysis of their fit synchrotron components across different flux states during this flare (visually inferred from Figure 5 in Manzoor et al. 2024) suggests that ${\nu }_{\mathrm{peak}}^{\mathrm{syn}}$ can shift to higher frequencies, reaching or even exceeding 10¹⁴ Hz, particularly during periods of heightened activity. This indicates that Ton 599 can transition toward or into the intermediate-synchrotron peaked (ISP) regime during such active phases, highlighting the dynamic nature of its SED. Our long-term optical PSD slopes (e.g., β_R ≈ 1.40) are derived from data spanning over a decade, thus representing an average characteristic of Ton 599’s variability across its diverse activity states, which would encompass both quiescent periods where it may behave as a typical LSP, and more active phases where it might exhibit ISP characteristics.

This observed dynamic behavior of ${\nu }_{\mathrm{peak}}^{\mathrm{syn}}$ in Ton 599 is relevant when considering studies such as that by Nilsson et al. 2018, who investigated the relationship between long-term optical PSD slopes (β) and ${\nu }_{\mathrm{peak}}^{\mathrm{syn}}$ for their blazar sample and found no significant correlation. The PSD characteristics of Ton 599, an FSRQ exhibiting such spectral peak variability, are therefore consistent with this general finding from Nilsson et al. 2018 of a lack of a simple, direct dependence of the long-term optical PSD slope on a single, fixed synchrotron peak frequency. Instead, the PSD slope likely reflects more fundamental properties of the energy dissipation and particle acceleration processes that persist across these varying spectral states.

A comparison of our optical PSD slopes with those obtained at radio and γ-ray energies provides additional context. In the radio domain, Kankkunen et al. 2025 analyzed a multi-decade 15 GHz OVRO light curve of Ton 599 and reported a clear break in the PSD at ∼1500 days, interpreted as a transition from a steep red-noise process (β ≈ 2) on shorter timescales to a flatter slope (β ≈ 1) on longer timescales. In this work, based on ∼4300 days of optical monitoring, we find that the BVRI light curves are well described by single power-law PSDs with slopes β_opt ≈ 1.4 − 1.8, without statistically significant evidence for a break. Although hints of a steepening above ∼1000 days may be present, the limited baseline prevents us from confirming a robust feature comparable to the radio break. At γ-ray energies, Wang et al. 2024 performed a systematic Fermi-LAT analysis of TeV blazars including Ton 599, and derived PSD slopes for monthly binned light curves (their Table 6). For FSRQs the average slope is close to unity (β_γ ≈ 1.0), consistent with flicker noise. This is flatter than both our optical slopes (β_opt ≈ 1.4 − 1.8) and the radio slope below the 1500-day break (β_radio ≈ 2), suggesting a wavelength-dependent steepening from γ-ray to optical to radio frequencies. The overall picture indicates that variability processes in Ton 599 become progressively more correlated (i.e., longer memory, steeper PSD) when moving from high to low energies. This trend may reflect different dominant timescales for particle acceleration and cooling in the jet: γ-ray variability is driven by fast processes in compact regions, optical variability (this work) arises from intermediate-scale synchrotron zones, while the radio emission integrates over larger jet volumes, where long-term jet-dynamics and propagation effects dominate.

Finally, while a power-law PSD provides a useful mathematical description of Ton 599’s observed variability spectrum, it is important to acknowledge the cautionary notes raised in the literature. As has been pointed out by Nilsson et al. 2018 and more broadly by Scargle 2020, this statistical descriptor does not unequivocally define the underlying physical process as simple power-law noise. The observed PSD could potentially arise from a superposition of discrete events or flares, or even from certain classes of linear, additive models under specific conditions. Thus, the interpretation of the PSD slopes in terms of specific physical mechanisms such as turbulence or propagating accretion instabilities is best supported when considered alongside other observational characteristics of Ton 599 discussed in this work.

4.3. IDV analysis.

Our analysis of IDV in Ton 599, based on observations where significant flux changes were detected on hourly timescales (Sect. 3.2 of this work), provides estimates for key physical parameters of the emitting regions. The characteristic variability timescales (τ_var) observed range from approximately 0.3 hours to 12.5 hours. Using these timescales and an adopted mean Doppler factor of δ = 23.1, we derived upper limits on the size of the emission region (R) and the strength of the magnetic field (B).

The IDV-derived emission region sizes (R) for Ton 599 range from 9.8 × 10¹⁵ cm to 4.3 × 10¹⁷ cm, with a mean value around 2 × 10¹⁷ cm. These can be compared with estimates of the emission region size from broadband SED modeling. For instance, Chen 2018 derived an emission region radius of R ≈ 7.9 × 10¹⁶ cm for Ton 599. While Ghisellini et al. 2010 primarily list the dissipation distance from the central black hole (R_diss) for Ton 599 as 1.14 × 10¹⁷ cm, their underlying model assumes that beyond an initial acceleration zone, the jet becomes conical with a semi-aperture angle ψ ≈ 0.1 radians. This allows for an estimation of the transverse size (radius) of the emission region (r_diss) in their model as r_diss ≈ ψ × R_diss ≈ 1.14 × 10¹⁶ cm. Our IDV-derived region sizes, particularly at the lower end of the range, are thus remarkably comparable to this estimated size from the Ghisellini et al. 2010 model parameters, and the size derived by Chen 2018 also falls within our IDV estimates. This broad agreement suggests that the rapid flares observed on intraday timescales might originate from compact regions whose dimensions are consistent with those implied by one-zone leptonic models. It is plausible that IDV probes these entire compact dissipation zones during coherent variations, or perhaps more localized, very dense subregions within a larger emitting volume, aligning with commonly inferred blazar emission region sizes of 10¹⁶ − 10¹⁷ cm (e.g., Ghisellini & Madau 1996; Raiteri et al. 1999). The range of our derived emission region sizes, R ≈ 9.8 × 10¹⁵ − 4.3 × 10¹⁷ cm can be further contextualized by a focused comparison with studies of the major November 2017 flare (p4) by Patel et al. 2018; Patel & Chitnis 2020. Our estimate for the most compact region, R ≈ 9.8 × 10¹⁵ cm, derived from the fastest optical IDV timescale of 0.3 h shows a remarkable agreement with the size of R′≈5.9 × 10¹⁵ cm estimated by Patel et al. 2018 using a timing argument on the gamma-ray flare data. This consistency, despite the different wavelengths and variability characteristics (stochastic IDV versus a discrete flare), suggests that simple timing arguments are effective at probing the scale of the most compact, active sites within the jet. In a subsequent work, Patel & Chitnis 2020 performed detailed SED modeling of the same flare and concluded that a more complex, two-component model was necessary. Their refined model includes a compact inner region with a size of R′ = 1 × 10¹⁶ cm and a much larger outer region of R′ = 1.26 × 10¹⁷ cm. It is noteworthy that the size of theirmodeled inner component remains highly consistent with our estimate for the most compact IDV region. This suggests that the rapid optical IDV likely originates from physical regions with dimensions comparable to the compact, primary energy dissipation sites of major gamma-ray flares. The necessity of a second, larger component in the Patel & Chitnis 2020 model to explain the full SED further strengthens the case for a structured jet, where our IDV analysis effectively isolates the scale of the core emission zone, while a larger volume contributes to the overall flare energetics.

The magnetic field strengths (B) derived from the IDV synchrotron cooling timescale arguments in our study range from 0.14 G to 0.50 G, with a mean of approximately 0.2 G. This range is notably consistent with the magnetic field of B ≈ 0.10 G derived for Ton 599 by Chen 2018 from their SED fit. However, these values are significantly lower than other estimates derived from broadband SED modeling, such as the B = 4.0 G reported by Ghisellini et al. 2010 for an average activity state. A particularly relevant comparison can also be made with the studies of the major November 2017 flare. An initial estimate by Patel et al. 2018 based on an equipartition argument yielded a much stronger field of B_e ≈ 6.9 G. Subsequent, detailed two-component SED modeling of the same event by Patel & Chitnis 2020 refined these values to B = 1.3 G for an inner component and B = 2.3 G for an outer component. While the specific values from different SED models vary, a consistent picture emerges where our IDV-derived magnetic field is weaker than estimates based on the overall energetics of the source. It is important to stress that one-zone models commonly adopted to fit blazar SEDs can be degenerate, which means that multiple choices of the parameters can lead to similarly good results. The differences between the various SED-based estimates (e.g., Ghisellini et al. 2010; Patel & Chitnis 2020) can therefore stem not only from different source activity states or distinct modeling assumptions, but also from this inherent model degeneracy. Despite this, the overarching trend suggests that our IDV analysis consistently probes a weaker magnetic field. This likely reflects the different physical conditions being probed: our IDV-derived B is specific to the rapidly varying, compact subregions, while the global SED-fitting approaches are constrained by the energetics of a larger emission zone. This may imply that our IDV analysis probes localized zones of weaker magnetic fields or enhanced turbulence within an inhomogeneous jet, while the flare and average-state models characterize the stronger, more organized magnetic field of the main jet body.

In addition to constraining the size of the emission region, the rapid variability timescales derived in our work can be used to estimate its distance (d) from the central engine. Following the methodology of Patel et al. 2018, which assumes d ∼ 2cΓ²t_var/(1 + z), our shortest observed timescale of τ_var ≈ 0.3 h (Table 3), combined with a Lorentz factor of Γ = 10 Jorstad et al. 2017, places the origin of this fast optical variability at a distance of d ≈ 0.001 pc. This finding can be contrasted with the results for major gamma-ray flares in Ton 599 presented by Patel et al. 2018; Patel & Chitnis 2020. Their initial analysis, based on a timing argument for a major flare, placed the emission region at d ≈ 0.051 pc, a location inside the blazar’s BLR, which is estimated to be at R_BLR ≈ 0.08 pc Patel et al. 2018. However, in a subsequent, more detailed study involving spectral energy distribution (SED) modeling, the same authors concluded that a simple one-zone model inside the BLR could not self-consistently explain the observed flare properties. Their refined two-component model required emission regions located entirely outside the BLR, at distances of d ≈ 0.097 pc and d ≈ 1.22 pc, with the external Compton process relying on photons from the dusty torus Patel & Chitnis 2020. Our result is therefore highly complementary; it suggests that the most rapid, low-amplitude optical variability we observe originates in a very compact region significantly closer to the central engine than the sites of the more energetic, large-scale gamma-ray flares. This supports a structured jet scenario, where different variability phenomena probe distinct physical locations and scales.

4.4. Long-term analysis.

Our extensive 12-year monitoring of Ton 599 (2011−2023) reveals a complex pattern of long-term optical variability, characterized by alternating phases of quiescence and pronounced activity, including several major outbursts (see Figures 1 and 5). The fractional variability (F_var) calculated for the entire period is substantial across all bands, reaching 99% in the R band for example (Table 4), which is characteristic of highly variable FSRQs (cf. Bhatta 2021, who reports F_var values for a sample of γ-ray bright blazars, with FSRQs such as 3C 279 and 3C 454.3 showing F_var of ∼80%).

During periods of enhanced activity, such as the major outbursts in Period II (around JD 2458098 and JD 2458248) and Period III, simultaneous multi-filter observations allowed us to derive an average optical spectral index α ≈ 1.2 (where F_ν ∝ ν^−α) during the peak of the first major outburst in Period II. This spectral index, corresponding to a relatively flat optical spectrum, is typical for FSRQs during active states when the synchrotron emission from the jet is dominant and can be indicative of fresh particle injection or acceleration. The redshift-corrected monochromatic luminosities derived for these flares, reaching log(νL_ν)≈48.24 [erg s⁻¹] for the JD 2458098 event, and notably log(νL_ν)≈48.48 [erg s⁻¹] for the R-band peak of 23.25 mJy during the record-breaking activity in Period III (January 2023), underscore the extreme power of Ton 599. These luminosity values are consistent with those observed in other powerful, high-redshift FSRQs during major flaring episodes (e.g., as is discussed for samples in Ghisellini et al. 2010).

The most intense optical activity recorded in our dataset, the R-band peak in January 2023, is particularly noteworthy as it appears to be the optical counterpart to the “brightest ever γ-ray flare” from Ton 599 observed in January 2023 (MJD 59954.5) reported by Manzoor et al. 2024. Manzoor et al. 2024 reported a peak γ-ray flux of (3.63 ± 0.20)×10⁻⁶ photon cm⁻² s⁻¹, and their multiwavelength light curve (their Figure 1) shows simultaneous enhancement in optical/UV and X-ray bands, strongly correlated with the γ-ray emission. Although a detailed multiwavelength correlation is beyond the scope of this particular paper, the temporal coincidence of our observed record optical flux with this extreme γ-ray event strongly suggests a common underlying physical origin, likely related to a major energy dissipation event within the jet, such as the propagation of a strong shock or a significant magnetic reconnection event. The complex, often multi-peaked structure of the major outbursts observed in our optical light curves (e.g., Period II and III, see Figure 5) is also a common feature in blazar flares and can be attributed to multiple emission subregions or evolving shock fronts.

Overall, the pronounced long-term optical variations in Ton 599, characterized by significant changes in both flux (spanning over two orders of magnitude) and fractional variability, and culminating in extreme outbursts with luminosities exceeding 10⁴⁸ erg s⁻¹, provide critical constraints on jet emission processes. These observations hint at episodic, and at times dramatic, enhancements in the accretion rate or, more likely, in the efficiency of energy transfer to and dissipation within the relativistic jet, driving the observed powerful outbursts across the electromagnetic spectrum.

4.5. Color analysis.

The long-term optical color analysis of Ton 599, presented in Sect. 3.4, reveals a remarkably complex spectral evolution that is strongly dependent on the source’s brightness level. Our data, spanning 12 years, show three distinct regimes in the CI versus magnitude diagrams, which offer insights into the interplay of different emission components and physical processes.

At faint states, Ton 599 consistently exhibits a RWB trend. This behavior is a well-documented phenomenon in some FSRQs and is typically interpreted as the increasing relative contribution of a relatively constant, redder thermal emission component as the bluer, nonthermal jet emission fades. This thermal component is likely associated with the accretion disk and the BLR. Steward Observatory spectra (Figure 9) provide direct support for this interpretation, revealing clear thermal emission features, such as the Mg II line and Fe II complexes, particularly during low-flux states. As the source brightens, the nonthermal jet continuum increasingly dominates, diluting these thermal features. A similar RWB trend at faint states, attributed to the thermal emission from the accretion disk (often referred to as the “big blue bump”), was also observed in the powerful FSRQ CTA 102 during its less active phases preceding its major 2016−2017 outburst (see Figure 3 in Raiteri et al. 2017). The study of 3C 454.3 by Villata et al. 2006 also implies that at fainter jet states, an underlying, relatively stable “blue bump” component (accretion disk) would be more prominent.

As Ton 599 brightens into intermediate flux levels, its spectral behavior transitions to being largely achromatic, where significant flux variations are not accompanied by systematic color changes. Such achromatic phases in blazars can be more challenging to interpret uniquely but have been attributed to scenarios where the spectral index of the dominant jet emission remains relatively constant during flux changes. Alternatively, geometric effects, such as variations in the viewing angle to different parts of an inhomogeneous or structured jet, or a twisting jet, could lead to flux changes without strong chromatic effects (e.g., Raiteri et al. 2017). Raiteri et al. 2017 invoked a “twisting inhomogeneous jet” model for CTA 102, where such geometric changes could explain complex variability patterns, including quasi-achromatic periods.

Finally, at its brightest states, Ton 599 displays a clear and statistically significant BWB trend across all color indices. This is the most commonly observed spectral behavior in FSRQs and many BL Lacertae objects during flaring activity (e.g., Wierzcholska et al. 2015). The BWB trend is generally explained by models dominated by synchrotron emission from relativistic electrons within the jet. It can arise from the injection or re-acceleration of electrons to higher energies during outbursts, leading to an increase in the relative contribution of higher-frequency (bluer) synchrotron emission, or a shift of the synchrotron peak frequency toward higher energies. Shock-in-jet models, where shocks accelerate particles and compress the magnetic field, also naturally predict such BWB behavior (e.g., Kirk et al. 1998; Sikora et al. 2001; Mastichiadis & Kirk 2002). This BWB characteristic was prominently observed in the FSRQ 3C 454.3 during its major optical outbursts, where it was interpreted as the synchrotron jet component, whose peak shifts with flux, overwhelming any constant thermal emission from the accretion disk Villata et al. 2006. The FSRQ CTA 102 also exhibited a strong BWB trend during its exceptional 2016−2017 outburst, consistent with these synchrotron-based models Raiteri et al. 2017. The overall sequence observed in Ton 599 – from RWB at faint states, through an achromatic phase, to BWB at bright states – strongly supports a two-component model where a relatively stable, redder thermal component is progressively dominated by an increasingly powerful and spectrally hardening (bluer) jet emission as the source brightens. Our estimate of a small underlying thermal disk contribution of 0.13−0.20 mJy is consistent with this picture of jet dominance at higher brightness levels.

Our analysis also revealed significant short-term color variability in Ton 599 on timescales of days to weeks. An intriguing finding is the observed anti-correlation between the amplitude of these short-term CI changes and the average flux density (Figure 7), suggesting that the jet emission is more spectrally stable (i.e., exhibits smaller rapid color fluctuations) during its brightest phases. While detailed modeling of such short-term color evolution is complex, studies of other blazars using high-cadence TESS and ground-based data, such as S5 0716+714 and S4 0954+65, have also revealed intricate and strongly chromatic short-term color variations, typically exhibiting BWB trends (e.g., Raiteri et al. 2021a,b). More broadly, the blazar literature also documents even more complex phenomena such as “clockwise hysteresis loops” in color-magnitude diagrams during such short-term events. These features are often attributed to the interplay of particle acceleration and cooling timescales within shocks, or the presence of multiple emission components with slight delays.

In summary, the complex and brightness-dependent color evolution of Ton 599 is consistent with a scenario involving contributions from both a thermal accretion disk and BLR component (dominant at faint states) and a powerful, spectrally variable relativistic jet (dominant at bright states). This behavior, including the RWB-achromatic-BWB transitions, mirrors the complexities observed in other well-studied, powerful FSRQs and underscores the dynamic interplay of various emission mechanisms and jet conditions.

4.6. Black hole mass.

Our analysis of Steward Observatory optical spectra, obtained during low-flux states of Ton 599 when thermal contributions from the accretion disk and BLR are most discernible, allowed for a robust determination of the central supermassive black hole (SMBH) mass. We carefully measured the Mg II λ2798 emission line parameters, accounting for complex Fe II blends by fitting the Mg II line with a two-Gaussian model to derive its total flux and FWHM. Using these parameters in the established empirical relation from Kong et al. 2006, our estimate for the black hole mass of Ton 599 is log(M_BH/M_⊙) = 8.77 ± 0.26.

This value is consistent with previous estimates for Ton 599 found in the literature, which employ a variety of methods and datasets, though a range of values exists. A summary of these mass estimates is presented in Table 9. For instance, Keck 2019, using the Mg II line from DCT spectra, also derived a comparable mass of log(M_BH/M_⊙) = 8.95. Broadband SED modeling by Ghisellini et al. 2010 for their sample of bright Fermi blazars listed log(M_BH/M_⊙) = 9.0. Similarly, Chen et al. 2009, using SDSS spectra and averaging estimates from Hβ and Mg II lines via empirical relations (Vestergaard & Peterson 2006 for Hβ; Kong et al. 2006 for Mg II), obtained a value of log(M_BH/M_⊙) = 9.11. More recently, Chen 2018, in a compilation of jet properties for γ-ray loud AGN, listed a black hole mass of log(M_BH/M_⊙) = 8.38 for Ton 599, a value primarily compiled from the literature employing virial methods (e.g., Shaw et al. 2012; Shen et al. 2011).

Other studies based on single-epoch virial methods have reported somewhat lower, yet broadly consistent, masses. Liu et al. 2006, employing the Hβ line and an empirical R_BLR − L₅₁₀₀ relationship, estimated log(M_BH/M_⊙) = 8.54. Using HST UV spectra, Pian et al. 2005 derived log(M_BH/M_⊙) = 8.63 from the C IV line and an R_BLR − L₁₃₅₀ relation. Xie et al. 2005 also reported log(M_BH/M_⊙) = 8.54 ± 0.21 based on a mass-luminosity relation for L5100 derived from their spectrophotometry of Ton 599. It is noteworthy that Xie et al. 2005 also presented a significantly lower mass estimate of log(M_BH/M_⊙) = 7.98 using a method based on the optical variability timescale (Δt_min > 50 min for Ton 599), which relies on different physical assumptions concerning the innermost stable orbit around a Kerr black hole.

The general consistency of our Mg II-based mass estimate (8.77 ± 0.26) with the majority of other single-epoch virial estimates (typically ranging from log(M_BH/M_⊙) ∼ 8.5 to ∼9.1) lends confidence to our result. The observed scatter among different literature values is typical for these methods and can be attributed to factors such as the intrinsic variability of the AGN continuum and broad emission lines, different observational epochs and data quality, the specific emission line and calibration used for the R_BLR − L relationship, and varying assumptions in SED modeling or alternative estimation techniques. Our derived mass firmly places Ton 599 in the category of hosting a very massive black hole, which is characteristic of luminous FSRQs and radio-loud AGNs, where black hole masses are generally found to be above 10⁸ M_⊙. Furthermore, the corresponding Eddington luminosity derived from our mass estimate, L_Edd ≈ 7.7 × 10⁴⁶ erg/s, is in line with expectations for such powerful sources, reinforcing the reliability of our spectral measurements and the applied methodology for mass estimation.

5. Summary and conclusions

In this work, we have presented a comprehensive analysis of the multiband optical variability of the FSRQ blazar Ton 599, using an extensive dataset collected between 2011 and 2023. We investigated its flux and color evolution across a wide range of timescales, from minutes to years, and used the observations to constrain some of its key physical parameters. Our main findings can be summarized as follows:

1. The long-term flux distribution in all observed bands (B, V, R, I) is better described by a log-normal distribution than a normal (Gaussian) one. This finding, supported by a strong linear RMS-flux relation, suggests that multiplicative processes drive the LTV. In contrast, the IDV flux distributions are consistent with a Gaussian profile, pointing to additive processes on the shortest timescales.

2. The PSD of the light curves follows a red-noise power-law trend, with slopes ranging from β ≈ 1.40 to 1.79. This confirms that the variability is a correlated stochastic process. We found no statistically significant periodicities in our12-year dataset.

3. We detected significant IDV on multiple nights. Using the fastest variability timescales, we constrained the size of the emitting region to the range R ≈ 10¹⁶ − 10¹⁷ cm and the magnetic field strength to B ≈ 0.14 − 0.50 G.

4. Based on the shortest observed optical variability timescale of ∼0.3 hours, we estimate that the most rapid, low-amplitude variations originate from a very compact region located at a distance of approximately d ≈ 0.001 pc from the central engine, placing it significantly closer than the sites of major γ-ray flares reported in other studies.

5. Over the monitoring period, the blazar exhibited extreme flaring activity, reaching a peak flux of 23.5 mJy in the R band. This corresponds to an exceptionally high monochromatic luminosity of log(νL_ν) = 48.48 [erg s⁻¹].

6. The blazar’s color evolution is complex and strongly depends on its brightness, revealing three distinct regimes: a RWB trend at low fluxes, which we associate with the thermal emission component; an achromatic phase at intermediate brightness; and a dominant BWB trend at the highest flux states, characteristic of the synchrotron jet emission.

7. Using low-flux optical spectra from the Steward Observatory, we analyzed the Mg II λ2798 emission line. After carefully accounting for Fe II contamination, we estimated the mass of the central supermassive black hole to be log(M/M_⊙) = 8.77 ± 0.26, a value consistent with previous estimates for this blazar.

Overall, our results underscore the complexity of blazar variability, pointing to a combination of different physical processes acting on different timescales. This comprehensive, long-term photometric and spectroscopic study of Ton 599 provides crucial constraints on the models of particle acceleration and emission in relativistic jets.

Data availability

The optical photometric data underlying this article are available at the CDS via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/703/A259

Acknowledgments

We thank the anonymous referee for constructive comments and suggestions that helped improve the manuscript. C.M.R., M.V. and M.I.C. acknowledge financial support from the INAF Fundamental Research Funding Call 2023. P.K. acknowledges support from the Department of Science and Technology (DST), government of India, through the DST-INSPIRE Faculty grant (DST/INSPIRE/04/2020/002586). We acknowledge the Maidanak Observatory for a sufficient amount of the observations. The research at Boston University was supported in part by the National Science Foundation grant AST-2108622, and by several NASA Fermi Guest Investigator grants, the latest is 80NSSC23K1507. This study was based in part on observations conducted using the 1.8m Perkins Telescope Observatory (PTO) in Arizona, which is owned and operated by Boston University. This article is partly based on observations made with the IAC80 operated on the island of Tenerife by the Instituto de Astrofisica de Canarias in the Spanish Observatorio del Teide. Many thanks are due to the IAC support astronomers and telescope operators for supporting the observations at the IAC80 telescope. This article is also based partly on data obtained with the STELLA robotic telescopes in Tenerife, an AIP facility jointly operated by AIP and IAC and on observations made with the LCOGT 0.4 m and 1.0 m telescope networks, one of whose nodes is located in the Spanish Observatorio del Teide. O.V. is partially supported by Chinese Academy of Sciences (CAS) President’s International Fellowship Initiative (PIFI) (grant No. 2024VMB0006). G.D., O.V., M.S. M.D.J., J.D.K. and M.L. acknowledge support by the Astronomical Station Vidojevica and the Ministry of Science, Technological Development and Innovation of the Republic of Serbia (MSTDIRS) through contract no. 451-03-136/2025-03/200002 made with Astronomical Observatory (Belgrade), by the EC through project BELISSIMA (call FP7-REGPOT-2010-5, No. 256772), the observing and financial grant support from the Institute of Astronomy and Rozhen NAO BAS through the bilateral SANU-BAN joint research project “GAIA astrometry and fast variable astronomical objects”, and support by the SANU project F-187. Also, this research was supported by the Science Fund of the Republic of Serbia, grant no. 6775, Urban Observatory of Belgrade – UrbObsBel. L.Č.P. acknowledge funding provided by the University of Belgrade – Faculty of Mathematics (the contract 451-03-136/2025-03/200104) and Astronomical Observatory Belgrade (the contract 451-03-136/2025-03/200002) through the grants by the Ministry of Science, Technological Development and Innovation of the Republic of Serbia. Observations with the SAO RAS telescopes are supported by the Ministry of Science and Higher Education of the Russian Federation. We acknowledge support by Bulgarian National Science Fund under grant DN18-10/2017 and Bulgarian National Roadmap for Research Infrastructure Project D01-326/04.12.2023 of the Ministry of Education and Science of the Republic of Bulgaria. This research was partially supported by the Bulgarian National Science Fund of the Ministry of Education and Science under grants KP-06-H68/4 (2022), KP-06-H88/4 (2024) and KP-06-KITAJ/12 (2024). NRIAG team acknowledges financial support from the Egyptian Science, Technology & Innovation Funding Authority (STDF) under grant number 45779. The Abastumani team acknowledges financial support by the Shota Rustaveli National Science Foundation of Georgia under contract FR-24-515. We acknowledge observations at the Tien Shan Observatory (Kazakhstan). We thank iTelescope for the use of their telescopes that made it possible for us to carry out optical observations. Based on observations made with the Nordic Optical Telescope, owned in collaboration by the University of Turku and Aarhus University, and operated jointly by Aarhus University, the University of Turku and the University of Oslo, representing Denmark, Finland and Norway, the University of Iceland and Stockholm University at the Observatorio del Roque de los Muchachos, La Palma, Spain, of the Instituto de Astrofisica de Canarias. The data presented here were obtained [in part] with ALFOSC, which is provided by the Instituto de Astrofisica de Andalucia (IAA) under a joint agreement with the University of Copenhagen and NOT. James R. Webb is honored to be permitted to conduct astronomical research on Iolkam Du’ag (Kitt Peak), a mountain with particular significance to the Tohono O’odham Nation. Some of the data are based on observations collected at the Observatorio de Sierra Nevada; which is owned and operated by the Instituto de Astrofísica de Andalucía (IAA-CSIC); and at the Centro Astronómico Hispano en Andalucía (CAHA); which is operated jointly by Junta de Andalucía and Consejo Superior de Investigaciones Científicas (IAA-CSIC). The IAA-CSIC co-authors acknowledge financial support from the Spanish “Ministerio de Ciencia e Innovación” (MCIN/AEI/10.13039/501100011033) through the Center of Excellence Severo Ochoa award for the Instituto de Astrofíisica de Andalucía-CSIC (CEX2021-001131-S), and through grants PID2019-107847RB-C44 and PID2022-139117NB-C44. Data from the Steward Observatory spectropolarimetric monitoring project were used. This program is supported by Fermi Guest Investigator grants NNX08AW56G, NNX09AU10G, NNX12AO93G, and NNX15AU81G.

References

Appendix A: WEBT data

The following appendix provides additional information on the data used in this study. Table A.1 presents details of the participating observatories and instruments in the WEBT campaign. Figure A.1 displays the segmentation of the optical light curve into 12 time intervals, which form the basis of the variability analysis discussed in the main text.

Fig. A.1. Light curve segmentation of the blazar Ton 599. The top panel shows the full light curve with 12 segments shaded in light blue, numbered from left to right. The lower panels provide a detailed view of each segment, arranged in a 4 × 3 grid.

All Tables

All Figures

	Fig. 1. Flux density evolution of the blazar Ton 599 from November 2011 to September 2023 in B, V, R, and I filters. Colors and symbols in the plots are explained in Table A.1.
In the text

Fig. 2. Left panel: Flux distribution from ∼12 years of R-band observations. The histogram (blue bars) shows the flux counts. The red curve is the best-fitting normal function, and the blue curve is the best-fitting log-normal function. Right panel: RMS (i.e., the square root of the excess variance) as a function of the mean flux, computed in (1) bins containing an equal number of data points (blue markers), and (2) bins of fixed time length (red markers). The solid blue and red lines indicate the best linear fits for each binning scheme.

In the text

	Fig. 3. R-band periodogram of Ton 599.
In the text

	Fig. 4. Intranight light curves displaying statistically significant variability based on the χ2 test. The vertical axis shows the measured magnitude, the legend indicates the filter used for each observation, and the horizontal axis represents the relative time in hours.
In the text

	Fig. 5. Zoomed-in R-band light curve from JD–2455000 = 3050 to 3300 (Period II) and from JD–2455000 = 4150 to 5120 (Period III) in the R band. Complex structure of the outbursts is visible. The inset highlights the second-largest peak, well sampled in all filters.
In the text

	Fig. 6. Evolution of the color indices with time. The first panel from the top shows the light curve in the R filter. The rest of the panels present the CI evolution. The red line is the best fit of a linear function to the data. Parameters of the linear function are provided in the Table 5.
In the text

	Fig. 7. Correlation between variation in the color indices and the mean flux density in 12 segments (shaded in the top R-band light curve). Red lines show the best linear fits.
In the text

	Fig. 8. Dependence of the color indices and spectral slopes on brightness. The dashed blue line corresponds to the mean value of the color indices. The right y axis shows the spectral slope, α, which was calculated under the assumption that the optical flux can be described by the relation Fν ∼ ν−α (see the text).
In the text

	Fig. 9. Steward optical spectra of the blazar Ton 599 at different activity states.
In the text

	Fig. 10. Top: Steward spectra in low (left) and high (right) flux states, with continuum fits (black) and BVRI filter curves. Bottom left: (B-R) vs. R color–magnitude diagram from observed (gray) and model-based (red/green) data. Bottom right: comparison of observed and continuum-based color indices. The color bar encodes B magnitude.
In the text

	Fig. 11. Binned spectral index α vs. log10FV for the BV data. Blue symbols are the bin averages with 1σ errors; gray points are individual measurements. The dashed red line shows the weighted linear fit for log10FV ≤ 0, and the vertical red line indicates log10FV at αdisk = 0.44 (i.e., Fdisk = 0.152 mJy).
In the text

	Fig. 12. Example of the Mg II + UV Fe II fit in the 2650−3050 Å range for the Steward spectrum observed on February 7, 2011. The core and wing Gaussian components of the Mg II line are plotted with solid lines, while the UV Fe II template from Popović et al. 2019 is plotted with a dashed line. The observed data are represented by dots, and the red line indicates the best fit.
In the text

	Fig. A.1. Light curve segmentation of the blazar Ton 599. The top panel shows the full light curve with 12 segments shaded in light blue, numbered from left to right. The lower panels provide a detailed view of each segment, arranged in a 4 × 3 grid.
In the text