We analyze the optical light curve data obtained with the Zwicky Transient Facility ZTF survey for 47 ray blazars monitored by the Large Area Telescope onboard the Fermi Gammaray Space Telescope Fermi. These 47 sources are selected because they are among the Fermi blazars with the largest optical variations in the ZTF data. Two colormagnitude variation patterns are seen in them with one being reddertostablewhenbrighter RSWB in 31 sources and the other being stable when brighter in 16 sources. The patterns fit with the results recently reported in several similar studies with different data. Moreover we find that the colors in the stable state of the sources share similar values for which after being corrected for the Galactic extinction most sources are in a range of 0.40.55. This feature could be intrinsic and may be applied in for example study of the intragalactic medium. We also determine the turning points for the sources showing the RSWB pattern after which the color changes saturate and become stable. We find a correlation between optical fluxes and ray fluxes at the turning points. The physical implications of the correlation remain to be investigated probably better with a sample of highquality ray flux measurements.
['(galaxies:) BL Lacertae objects: general', '(galaxies:) quasars: general', 'gamma-rays: galaxies', 'Astrophysics - High Energy Astrophysical Phenomena']
Fermi Blazars in the Zwicky Transient Facility Survey Properties of Large Optical Variations
{'No Header': 'R esearchin  \nA  \nstronomyand  \nA  \nstrophysics', 'Fermi Blazars in the Zwicky Transient Facility Survey: Properties of Large Optical Variations': 'Si-Si Sun 1 , Zhong-Xiang Wang 1 , 2 and Shun-Hao Ji 1  \n- 1 Department of Astronomy, School of Physics and Astronomy, Yunnan University, Kunming, Yunnan, 650091, China wangzx20@ynu.edu.cn\n- 2 Shanghai Astronomical Observatory, Chinese Academy of Sciences, 80 Nandan Road, Shanghai 200030, China  \nReceived 20xx month day; accepted 20xx month day  \nAbstract We analyze the optical light-curve data, obtained with the Zwicky Transient Facility (ZTF) survey, for 47 γ -ray blazars monitored by the Large Area Telescope onboard the Fermi Gamma-ray Space Telescope (Fermi) . These 47 sources are selected because they are among the Fermi blazars with the largest optical variations in the ZTF data. Two color-magnitude variation patterns are seen in them, one being redder to stable when brighter (RSWB; in 31 sources) and the other being stable when brighter (in 16 sources). The patterns fit with the results recently reported in several similar studies with different data. Moreover, we find that the colors in the stable state of the sources share similar values, which (after corrected for the Galactic extinction) of most sources are in a range of 0.4-0.55. This feature could be intrinsic and may be applied in, for example, the study of intragalactic medium. We also determine the turning points for the sources showing the RSWB pattern, after which the color changes saturate and become stable. We find a correlation between optical fluxes and γ -ray fluxes at the turning points. The physical implications of the correlation remain to be investigated, probably better with a sample of high-quality γ -ray flux measurements.  \nKey words: BL Lacertae objects: general - quasars: general - gamma-rays: galaxies', '1 INTRODUCTION': "As a subclass of the active galactic nuclei (AGNs) that are powered by accretion of matter onto central supermassive black holes of galaxies (Urry & Padovani 1995; Ghisellini et al. 1998), blazars are those with a jet pointing close to our line of sight. They thus exhibit highly variable radiation across the whole electromagnetic spectrum due to the Doppler beaming effect (e.g., Hovatta et al. 2009; Hovatta & Lindfors 2019), since their emissions we observe are mostly from the jets. Blazars have a double-hump structure of broadband spectral energy distributions (SEDs), with the low-energy hump peaking at frequencies from radio to ultraviolet/X-rays and the high-energy hump at energies from hard X-ray to γ -ray. Generally, a so-called leptonic scenario is often considered for the emissions of blazars, in which the low-energy hump arises from synchrotron radiation of relativistic electrons accelerated in a jet, and the high-energy hump from inverse Compton Scattering (ICS) of photons by the same population of the electrons, where the seed photons are provided by either the jet (e.g.,Maraschi et al. 1992; Bloom & Marscher 1996) or/and the other AGN components (e.g., broad-line region, dusty torus, and accretion disk; Sikora et al. 1994; Dermer & Schlickeiser 1994; Bła˙zejowski et al. 2000; Ghisellini & Tavecchio 2009).  \nFig. 1: Left: ∆ zr and F var in the zr band calculated for 810 blazars, which can be described with a fit of F var = -0 . 25∆ zr -0 . 11 (red line). Right: γ -ray variability versus ∆ zr for the 810 blazars. No obvious variability correlation between the two bands is seen.  \n<!-- image -->  \nSince the launch of the Fermi Gamma-ray Space Telescope (Fermi) in 2008 (Atwood et al. 2009), blazars have been fully observed at γ -rays by the Large Area Telescope (LAT) onboard Fermi in energy range of approximately from 100 MeV to 500 GeV. Not only have nearly 4000 blazars been detected (Ballet et al. 2023), but their emissions also have been monitored since 2008. Thus, Fermi -LAT blazars constitute a valuable sample for different studies. In particular, γ -ray light curves binned on three different timescales for more than 1400 blazars are being provided by the Fermi -LAT light curve repository (LCR; Abdollahi et al. 2023). Flux and spectral variation information for these blazars can be extracted for studies.  \nVariability studies of blazars help probe their properties and related radiation mechanisms, jet properties, and emission regions (see, e.g., Abdo et al. 2010a; Bottcher 2019). For example, the short, intraday variations are used to constrain the sizes of the emission regions, while the ultra-short ones (on timescales of a few minutes) may further indicate the detailed emission scenarios (e.g., Aharonian et al. 2017 and references therein). A more common type of variations seen in blazars may be described with a stochastic process with time scales of ∼ 10-100days (e.g., Zhang et al. 2023 and references therein), accompanied with (large) flares (e.g., Hayashida et al. 2015). Study of the variations probe the jets' activity and their potential connection with accretion disks, and multi-band monitoring of flares provide clues for emission scenarios and locations of emission regions through analyses of broadband SEDs and variation correlations between different bands (e.g., Hayashida et al. 2015; Krauß et al. 2016; Cohen et al. 2014; Liodakis et al. 2018; Bhatta 2021; Sun et al. 2023).  \nAmong variability studies of blazars, one is to investigate the color-magnitude changes at optical bands (e.g., Gu et al. 2006; Otero-Santos et al. 2022; Zhang et al. 2022; McCall et al. 2024). The investigations reveal the optical spectral changes during activities of blazars and possibly indicate the emission contribution from accretion disks (e.g., Raiteri et al. 2008; Isler et al. 2017; Fan et al. 2018). Such studies can be widely carried out, with the advent of large surveys of optical transients. A rich amount of data at decent cadences covering large sky areas are provided by these surveys. In this work, we used the data from the Zwicky Transient Facility (ZTF; Bellm et al. 2019) survey to study the optical variations of Fermi blazars. Interesting color-magnitude variation patterns from a sample of 47 blazars were found. In addition, simultaneous γ -ray and X-ray data were analyzed for helping understand the optical variations. We report the results in this paper.  \nBelow in Section 2, we first describe the data we used and the target selection procedure, where 47 blazars with the largest zr -band variations in the ZTF data were selected. In Section 3, we describe the analyses we conducted for the 47 blazars using the ZTF optical and Fermi -LAT γ -ray light-curve data; available Neil Gehrels Swift Observatory (Swift) X-ray data were also analyzed. We present the results and discuss them in Section 4.  \nFig. 2: Redshifts versus γ -ray luminosities of Fermi -LAT detected blazars. Among our 47 blazar targets, 40 have redshifts, and they are marked with pentagons.  \n<!-- image -->", '2.1 Fermi -LAT γ -ray and ZTF optical light-curve data': 'Fermi -LAT γ -ray light curves from its LCR (Abdollahi et al. 2023) were used in our analysis. This public database provides light curves for Fermi -detected variable sources binned on timescales of 3day, 7-day (weekly), and 30-day (monthly). There were 1414 blazars among the Fermi -LAT source catalog that have optical counterparts with given coordinates and are included in Fermi -LAT LCR 1 . We downloaded the monthly light curves for these 1414 blazars, and used the flux data points with test statistic (TS) values greater than 4 (i.e., detection significances of the data points were greater than 2 σ ).  \nWe obtained the coordinates of 1414 blazars from the Fermi LAT source catalog (Ballet et al. 2023), and searched for their magnitude data in the ZTF Data Release 20 with the provided API 2 . The ZTF survey can provide magnitude data points at a cadence of as high as ∼ 2day for sources in its field. To ensure high-quality of the light-curve data, we selected data points with catflags = 0 and chi < 4 . We found 1403 blazars covered by the ZTF survey. There are three bands, zg , zr , and zi , in the ZTF survey, and we only used data at zg and zr bands since there often are limited data at zi band. Given our purpose of studying optical variations, which requires as many data points as possible for individual sources, we excluded those sources with the number of data points less than 100 at zr band. After the exclusion, there were 810 blazars. The time period of the ZTF data was approximately MJD 58190-60200.', '2.2 Target selection': 'We checked variations of the 810 blazars using the zr -band data. First, the variation amplitudes ∆ zr were obtained by subtracting the maximum magnitude from the minimum magnitude of the zr light curve of each blazar. Secondly, we calculated the fractional variability amplitude ( F var ; e.g., Schleicher et al. 2019) of each light curve as the further check. The formula for calculating F var for a light curve of N data points is as follows:  \nF var = √ S 2 -σ 2 err ¯ x 2 , (1)  \nerr ( F var ) = √ √ √ √ √ { √ 1 2 N · σ 2 err ¯ x 2 F var } 2 +    √ σ 2 err N · 1 ¯ x    2 , (2)  \nSun et al. Table 1: Information for 47 blazar targets  \nNotes: Column 7 gives the magnitude of the zr band at the turning point, Column 8 gives the corresponding γ -ray flux (in units of photon cm -2 s -1 ), and NULL indicates no data.  \nσ 2 err = 1 N N ∑ i =1 σ 2 err ,i , (3)  \nwhere ¯ x and S 2 are the mean flux and total variance of a light curve, respectively, σ err is the error of a flux measurement, σ 2 err is the mean square error, and err ( F var ) is the error of F var . In calculating F var , zr magnitudes were converted to flux densities (in units of Jy); the ZTF magnitudes are in the AB photometric system (Bellm et al. 2019).  \nIn the left panel of Figure 1, we show F var versus ∆ zr of the 810 blazars. A relationship of F var = -0 . 25∆ zr -0 . 11 can be obtained by fitting the data points. The p-value of the fitting is extremely small ( /similarequal 10 -16 ), indicating the high significance of the relationship. This comparison shows that the variability of the blazars is approximately proportional to ∆ zr . In addition, we also checked γ -ray variability versus ∆ zr for these 810 blazars, where we used Variability index given in the Fermi -LAT  \nsource catalog (Ballet et al. 2023) to represent γ -ray variability. As can be seen in the right panel of Figure 1, there are no obvious correlated variabilities between the γ -ray and optical band.  \nBased on the optical-variation plot shown in Figure 1, we may use ∆ zr to select targets for simplicity. We found that there were 52 sources with the ∆ zr variation greater than 4 mag (i.e., sources with largest variations, right of -4 mag in the left panel of Figure 1). We decided to focus on these largest-variation sources. However, as we further analyzed the data for them, we realized that five of them had too few zg data points simultaneous to the zr data points in our process of obtaining their color variations. We excluded these five sources and had 47 blazars as the final targets. Information for these targets is provided in Table 1.  \nBlazars are divided into two sub-types based on emission-line features detected in their optical emissions (Stickel et al. 1991; Urry & Padovani 1995): Flat Spectrum Radio Quasars (FSRQs) that show emission lines with an equivalent width > 5 ˚ A and BL Lacertae objects (BL Lacs) that have weak or no emission lines. We checked our targets and provide their sub-type information in Table 1. Among them, 33 are FSRQs and 7 are BL Lacs, while the remaining 7 are blazar candidates of uncertain type (BCUs; e.g., Ballet et al. 2023). In addition, blazars are also classified by the peak frequencies ν sy pk of their synchrotron humps (Abdo et al. 2010b; Fan et al. 2016; Ghisellini et al. 2017): Low Synchrotron Peaked (LSP) blazars ( 10 14 Hz > ν sy pk ), Intermediate Synchrotron Peaked (ISP) blazars ( 10 15 Hz > ν sy pk > 10 14 Hz), and High Synchrotron Peaked (HSP) blazars ( ν sy pk > 10 15 Hz). We note that all the targets belong to LSP.  \nFinally we checked the redshifts ( z ) and γ -ray luminosities ( L γ ) of the targets. Among them, 40 have given redshifts. We show them in Figure 2 along with other Fermi -detected blazars. As can be seen, most of them are among the luminous blazars with L γ in ∼ 10 46 -10 48 erg s -1 , and there are also 8 of them with L γ in ∼ 10 45 -10 46 erg s -1 . The bias towards luminous blazars (or FSRQs) could be set by our use of the ZTF data, which have 5 σ detection limits of 20.8 mag in zg and 20.6 mag in zr (Masci et al. 2019). Thus, relatively more luminous blazars would be well covered by the ZTF survey.', '3.1 Color-magnitude analysis': 'A zg -zr color curve for each target was calculated, where zg and zr magnitudes obtained within one same day were considered as simultaneous measurements and used in the calculation. Color-magnitude diagrams of zg -zr versus zr for the targets were thus obtained. There are three types of colormagnitude variations shown by our targets. The first type is curved changes following the redder when brighter (RWB) pattern, and after the color reaches a point, it stays around a constant line without obvious trending. Following Zhang et al. (2022), we call the latter constant-color part as the stable state. Figure 3 shows a good example of this type. For this type, we used the following function to fit the color-magnitude variations:  \ny = { k ( x -x b ) 2 + c , if x ≥ x b c , if x < x b , (4)  \nwhere x and x b are the magnitude and the break point between the curved changes and the constant c line, respectively, and y is the color (note that we plot magnitudes and color values from large values to small ones). With this fitting (by using curve fit in Python SciPy), we obtained both x b and c .  \nThe second type is similar to the first type, but the RWB changes appear linear. Figure 4 shows an example of this type of changes. We fit the color-magnitude changes similarly to that of the first type, but used a straight line instead for the RWB changing part:  \ny = { k ( x -x b ) + c , if x ≥ x b c , if x < x b . (5)  \nFig. 3: Light curve ( upper left ) and color-magnitude ( upper right ) diagrams of 4FGL J0112.8+3208. The curved and constant parts of the color changes are indicated by blue and red lines respectively in the color-magnitude diagram. The corresponding data points used to construct the color-magnitude changes are respectively marked as blue and red in the light curves. The black data points ( zg or zr , not used) in the light curves are those do not have the same-day measurements in the respective other band. The γ -ray light curve of the source during the ZTF time period is shown in the lower left panel.  \n<!-- image -->  \nFig. 4: Same as Figure 3 for 4FGL J1350.8+3033. The RWB part is fitted with a straight line ( upper right panel).  \n<!-- image -->  \nThe x b and c values were obtained from the fitting. In total, 31 of our targets show these two types of changes, and the first and second type is 16 and 15 respectively.  \nFor the remaining 16 blazars, their color changes appear to be around a constant line without obvious RWB or other types of patterns seen. Figure 5 is an example of this type. For this third type, we fitted the color values with a constant c .  \nBesides the three example sources shown in Figures 3-5, the fits to the color-magnitude variations of the other 44 sources are shown in Figures A.1 & A.2.  \nFig. 5: Same as Figure 3 for 4FGL J1303.0+2434. The color changes of this source can be described with a straight line ( upper right panel).  \n<!-- image -->', '3.2 Determination of the γ -ray flux at x b': "We compared the γ -ray light curves with the optical ones for each target. The comparison indicates two different types of cases. One is that γ -ray flaring activity can be seen to correspond to optical flares (e.g., Figures 3 & 4) and the other is that there is not clear flaring activity in γ -ray with respect to optical flare-like variations (e.g., Figure 5). Moreover, examining the optical light curves and color changes, a common property can be drawn: the constant-color part in the first two types of color changes (Section 3.1) consists of brighter magnitude data points contained in one or several flares. This is clearly shown in Figures 3 & 4 (red data points in the light curves). Therefore it appears that there exists a turning point, x b in our fitting, after which the color changes saturate and become to be around a constant value (e.g., Fan et al. 2018).  \nIn order to connect the optical variation behavior with the γ -ray one, we approximated the corresponding γ -ray flux at x b for the first two types of color changes. The method we used is as follows. For each γ -ray flux data point of a target, we found its closest zr measurement in time, but also required the former's (30-day) time bin covering the latter's date. Thus, a plot of γ -ray flux versus zr flux density was constructed (see Figure 6 as an example). Then based on x b (Figure 6) determined above for this target, we chose two γ -ray fluxes closest to x b in time and averaged the values as the corresponding γ -ray flux at the turning point. The flux values obtained for 31 blazars of the first two types of color changes are provided in Table 1.", '3.3 Swift X-ray data analysis': 'It is interesting to learn the variation behavior at other wavelengths during the constant-color part. We thus searched archival X-ray data from Swift X-ray Telescope (XRT) observations during the time periods when a blazar target was in the stable state. Note that the 16 blazars fitted with a constant are considered to be always in the stable state (Zhang et al. 2022). Observations of six sources were found, and the times of the observations for each source are provided in Table 2.  \nWe employed the online Swift-XRT data products generator tool 3 to obtained the 0.3-10 keV Xray spectra. This tool allows the submission of a large sample of objects, and automatically extract and fit the spectrum in each observation (for details about the online tool, see Evans et al. 2009). In each extraction, we set an absorbed power-law model to fit the spectrum, where the Galactic hydrogen column density N H was automatically fixed at the value from Willingale et al. (2013). In total, 71 observations  \nFig. 6: Example of determining the γ -ray flux at the turning point. γ -ray fluxes and corresponding zr flux densities of 4FGL J0112.8++3208 are shown. The red vertical line is the zr -band flux density of x b . The average of two γ -ray fluxes closest in time to that of x b is taken as the value at the break point.  \n<!-- image -->  \nof the six blazars were analyzed, and unabsorbed fluxes and photon indices ( Γ x ) were obtained from the observations.  \nTable 2: Times of Swift -XRT observations', '4 RESULTS AND DISCUSSION': 'By investigating color-magnitude changes of 47 blazars that had the largest optical variations in the time period of ∼ MJD 58190-60200, we have found 31 of them showing a RWB pattern and the other 16 being around a constant of zg -zr . More specifically, the color changing behavior of the 31 blazars may be described as redder-stable-when-brighter (RSWB; e.g., Zhang et al. 2022; McCall et al. 2024), since they all appeared to have a turning point over the course of brightening, after which the zg -zr color stayed around a constant. The other 16 blazars may be described as stable-when-brighter (SWB; Zhang et al. 2022). These behaviors have been observed from earlier studies of individual sources (e.g., Gu et al. 2006; Villata et al. 2006; Raiteri et al. 2008; Fan et al. 2018) and more recent studies of large samples (e.g., Otero-Santos et al. 2022; Zhang et al. 2022; McCall et al. 2024).  \nThe RSWB pattern is considered to be mostly associated with FSRQs, which has been explained with a two-component scenario (e.g., Gu et al. 2006; Zhang et al. 2022). There is a thermal component arising from the accretion disk and broad-line region in an FSRQ, whose emission has a steep spectrum, and the other one from the jet, whose emission is non-thermal with a less steep spectrum. As the source is brightening, the latter component contributes more to the observed emission, thus inducing the RWB behavior. After the source reaches the point when the latter component becomes dominant, the source enters the stable state. Thus, the SWB behavior seen in the 16 blazars is likely because the sources were in the high-flux, stable state (Zhang et al. 2022), although we note that their luminosities are not particularly different, spreading from 10 45 to 10 47 erg s -1 along with the sources with variable colors (Color const and Color var data points in Figure 2, respectively). These SWB sources would have their RWB part, but the fluxes would be too faint to be detected by ZTF, whose detection limit is ∼ 20.8mag.  \nFig. 7: Left: distribution of observed zg -zr at the stable state. Right: distribution of ( zg -zr ) 0 , which has been corrected with the Galactic extinction.  \n<!-- image -->  \nFig. 8: Dereddened color ( zg -zr ) 0 as a function of redshift z . Downward triangles mark those blazars not showing a turning point in their color-magnitude diagrams. A weak trend of ( zg -zr ) 0 ∝ 0 . 05 z is possibly seen.  \n<!-- image -->  \nThe near-future surveys, such as the Legacy Survey of Space and Time (LSST; with the Vera C. Rubin Observatory), will be able to detect their > 20.8 mag turning points, by providing much deeper lightcurve data.  \nAmongthe 31 RSWB blazars, only two are BL Lacs, which are PKS 0235+164 and TXS 0404+075. The first one has been in the samples of Otero-Santos et al. (2022), Zhang et al. (2022), and McCall et al. (2024), and they all reported the RWB pattern for it, although Otero-Santos et al. (2022) noted a mild bluer-when-brighter (BWB) trend in the brightest spectra. We have not found color-pattern information for the second source in the literature. As pointed out by Zhang et al. (2022), BL Lacs can also show the RSWB behavior, as long as the accretion disks are sufficiently brighter (e.g., disk luminosities are probably as large as 5 × 10 45 -10 46 erg s -1 ).  \nIn our analysis, we noted that the color values in the stable state are similar. We plot the color distribution in the left panel of Figure 7. Most sources are within ∼ 0.4-0.6, but several have larger values. We adopted the estimates of the Galactic extinction E ( B -V ) values to our sources from Schlafly & Finkbeiner (2011), and the values are provided in Table 1. We then dereddened zg -zr with the respective extinction values for the sources (where we approximated with reddening at zg , A zg /similarequal 3 . 303 E ( B -V ) , and that at zr , A zr /similarequal 2 . 285 E ( B -V ) ; Schlafly & Finkbeiner 2011). The dereddend color ( zg -zr ) 0 distribution is shown in the right panel of Figure 7. As can be seen, the sources are slightly more concentrated to a range of 0.4-0.55. A Gaussian fit can approximately describe the distribution, with the peak at /similarequal 0.49 and a standard deviation of /similarequal 0.052. There are three notable sources  \nFig. 9: X-ray photon index ( Γ x ) versus 0.3-10 keV X-ray flux of Ton 599 during its optical stable state. A brighter and harder (lower Γ x values) pattern is seen.  \n<!-- image -->  \nwith ( zg -zr ) 0 > 0 . 7 , two of which are the BL Lacs PKS 0235+164 (at ∼ 0.9) and TXS 0404+075 (at ∼ 0.7). These two are the only known BL Lacs mentioned above that showed the RSWB pattern in our sample. Given this, we might consider that BL Lacs tend to have redder colors in the stable state. However, we note that towards PKS 0235+164 ( z = 0 . 94 ), there is a known absorber at redshift z = 0 . 524 (Junkkarinen et al. 2004; Raiteri et al. 2005). Due to this case, the large color values could likely be caused by extra-absorption along the line of sight towards the sources.  \nWe checked ( zg -zr ) 0 as a function of z with the plot shown in Figure 8. The redshift ranges from 0.3 to 2.6, which shift emission at ∼ 5000 / (1 + z ) ˚ Ato the observed bands. The largest value of z /similarequal 2 . 6 corresponds to an ultraviolet (UV) wavelength of ∼ 1300 ˚ A. Except the three ( zg -zr ) 0 > 0 . 7 sources mentioned above, the sources do not show large color variations or a strong trend. If we consider a power-law spectrum F ν ∝ ν -α for the optical-UV emission at the rest wavelengths, the narrow range of ( zg -zr ) 0 ∼ 0.4-0.55 would suggest that there are no large changes in α from UV to optical in these sources. We tried to fit the colors as a function z (not including the three sources with ( zg -zr ) 0 > 0 . 7 ). Apossible trend of ( zg -zr ) 0 ∝ 0 . 05 z was obtained, although the colors are relatively largely scattered around the fit (Figure 8). In any case, this weak trend could be the indication of the extra-reddening caused by the intragalactic medium (IGM), if we assume the UV-optical spectra of blazars in the stable state have highly similar α values. In X-rays, the IGM has been detected as the extra-absorption in spectra of blazars (and gamma-ray bursts; see Gatuzz et al. 2024 and references therein). A large sample of blazars with the stable-state colors would provide a test for this possibility (we will carry out the study in another work).  \nWe also checked if there were significant X-ray spectral changes during the stable state of the blazars. Swift XRT data for six sources were available (Table 2). Except 4FGL J1159.5+2914 (or Ton 599), no significant changes were detected. However, because the X-ray fluxes of the sources were in a range of ∼ (1-10) × 10 -12 erg cm -2 s -1 and the exposures of the Swift observations were short (generally < 1000 sec), the measurements suffer large uncertainties (see, e.g., Figure 9). In any case, no evidence of spectral changes at X-rays were found for the sources at the optical stable state. In Ton 599, we detected a harder when brighter pattern (Figure 9). However, its SED likely consists of contributions from multiple components of two emission zones in its brightening state: its X-ray emission probably arises dominantly from a self ICS process in an out zone with the seed photons provided by the jet, while the optical emission mainly comes from the synchrotron radiation in an inner zone (Patel & Chitnis 2020). The connection between the optical and X-ray emissions in this source is likely not direct, which may explain the spectral variation at X-rays in the optical stable state.  \nAs pointed out by Zhang et al. (2022), the turning point from the RWB to the stable state would commonly exist, brighter than which the non-thermal, synchrotron emission become dominant in the observed optical emission. For the 31 blazars showing such a turning point, we notice that above it, the data points were in one or several flares, which often have the corresponding γ -ray flares. The γ -ray  \nFig. 10: zr flux density F zr versus γ -ray flux F γ at the turning points of the 31 blazars. A F zr ∝ F 1 . 02 γ relationship is seen. The 16 blazars not showing a turning point, marked with the Downward triangles (i.e., the flux-density upper limits in zr -band), are also included for comparison.  \n<!-- image -->  \nactivity possibly reflects the general state of such a blazar is at. In Figure 10, we show the zr -band flux F zr and the corresponding γ -ray flux F γ of the 31 blazars (i.e., values in the last two columns of Table 1). A correlation between the optical and γ -ray fluxes is seen, and F zr ∝ F 1 . 02 ± 0 . 16 γ . We also investigated if there is some physical condition for the turning point that could be indicated in γ -rays. For example, a comparison of F zr versus F γ , normalized respectively by the maximum zr flux density and the maximum γ -ray flux of a source, was conducted. However, because of large uncertainties in γ -ray fluxes (cf., Figure 10), no clear results were seen. Such investigations could be further conducted by selecting only bright γ -ray blazars in a sample.  \nAs a summary, we have analyzed the ZTF optical data for 47 Fermi -LAT blazars that showed the largest optical variations in the past several years. Among them, 33, 7, and 7 are FSRQs, BL Lacs, and BCUs, respectively, and all are LSPs. Their color changes show either a RSWB pattern or a constant-line pattern, in the former the colors reach a constant-line state (or the stable state) after a turning point. The colors at the stable state share similar values, and the dereddened values were obtained. Except three sources have large color values of > 0.7, which could be caused by the existence of extra-absorbers along the line of sight towards them based on one known case PKS 0235+164, most of the sources have the color values in a range of 0.4-0.55. If this feature of a narrow color range is true for blazars in the stable state, it could be used to explore the IGM. In addition, a correlation between optical flux and γ -ray flux at the turning point is found, but what this correlation implies remains to be investigated.  \nAcknowledgements This work was based on observations obtained with the Samuel Oschin Telescope 48-inch and the 60-inch Telescope at the Palomar Observatory as part of the Zwicky Transient Facility project. ZTF is supported by the National Science Foundation under Grant No. AST-2034437 and a collaboration including Caltech, IPAC, the Weizmann Institute for Science, the Oskar Klein Center at Stockholm University, the University of Maryland, Deutsches Elektronen-Synchrotron and Humboldt University, the TANGO Consortium of Taiwan, the University of Wisconsin at Milwaukee, Trinity College Dublin, Lawrence Livermore National Laboratories, and IN2P3, France. Operations are conducted by COO, IPAC, and UW.  \nWe thank the referee for detailed and helpful comments. This research is supported by the Basic Research Program of Yunnan Province (No. 202201AS070005), the National Natural Science Foundation of China (12273033), and the Original Innovation Program of the Chinese Academy of Sciences (E085021002). S.J. acknowledges the support of the science research program for graduate students of Yunnan University (KC-23234629).', 'References': "```\nAbdo, A. A., Ackermann, M., Ajello, M., et al. 2010a, ApJ, 722, 520 2 Abdo, A. A., Ackermann, M., Agudo, I., et al. 2010b, ApJ, 716, 30 5 Abdollahi, S., Ajello, M., Baldini, L., et al. 2023, ApJS, 265, 31 2, 3 Aharonian, F. A., Barkov, M. V., & Khangulyan, D. 2017, ApJ, 841, 61 2 Atwood, W. B., Abdo, A. A., Ackermann, M., et al. 2009, ApJ, 697, 1071 2 Ballet, J., Bruel, P., Burnett, T. H., Lott, B., & The Fermi-LAT collaboration. 2023, arXiv e-prints, arXiv:2307.12546 2, 3, 5 Bellm, E. C., Kulkarni, S. R., Graham, M. J., et al. 2019, PASP, 131, 018002 2, 4 Bhatta, G. 2021, ApJ, 923, 7 2 Bła˙zejowski, M., Sikora, M., Moderski, R., & Madejski, G. M. 2000, ApJ, 545, 107 1 Bloom, S. D., & Marscher, A. P. 1996, ApJ, 461, 657 1 Bottcher, M. 2019, Galaxies, 7, 20 2 Cohen, D. P., Romani, R. W., Filippenko, A. V., et al. 2014, ApJ, 797, 137 2 Dermer, C. D., & Schlickeiser, R. 1994, ApJS, 90, 945 1 Evans, P. A., Beardmore, A. P., Page, K. L., et al. 2009, MNRAS, 397, 1177 7 Fan, J. H., Yang, J. H., Liu, Y., et al. 2016, ApJS, 226, 20 5 Fan, X.-L., Li, S.-K., Liao, N.-H., et al. 2018, ApJ, 856, 80 2, 7, 8 Gatuzz, E., Wilms, J., Hammerich, S., & Arcodia, R. 2024, A&A, 683, A213 10 Ghisellini, G., Celotti, A., Fossati, G., Maraschi, L., & Comastri, A. 1998, MNRAS, 301, 451 1 Ghisellini, G., Righi, C., Costamante, L., & Tavecchio, F. 2017, MNRAS, 469, 255 5 Ghisellini, G., & Tavecchio, F. 2009, MNRAS, 397, 985 1 Gu, M. F., Lee, C. U., Pak, S., Yim, H. S., & Fletcher, A. B. 2006, A&A, 450, 39 2, 8 Hayashida, M., Nalewajko, K., Madejski, G. M., et al. 2015, ApJ, 807, 79 2 Hovatta, T., & Lindfors, E. 2019, New Astron. Rev., 87, 101541 1 Hovatta, T., Valtaoja, E., Tornikoski, M., & Lahteenmaki, A. 2009, A&A, 494, 527 1 Isler, J. C., Urry, C. M., Coppi, P., et al. 2017, ApJ, 844, 107 2 Junkkarinen, V. T., Cohen, R. D., Beaver, E. A., et al. 2004, ApJ, 614, 658 10 Krauß, F., Wilms, J., Kadler, M., et al. 2016, A&A, 591, A130 2 Liodakis, I., Romani, R. W., Filippenko, A. V., et al. 2018, MNRAS, 480, 5517 2 Maraschi, L., Ghisellini, G., & Celotti, A. 1992, ApJ, 397, L5 1 Masci, F. J., Laher, R. R., Rusholme, B., et al. 2019, PASP, 131, 018003 5 McCall, C., Jermak, H., Steele, I. A., et al. 2024, arXiv e-prints, arXiv:2404.12835 2, 8, 9 Otero-Santos, J., Acosta-Pulido, J. A., Becerra Gonz'alez, J., et al. 2022, MNRAS, 511, 5611 2, 8, 9 Patel, S. R., & Chitnis, V. R. 2020, MNRAS, 492, 72 10 Raiteri, C. M., Villata, M., Ibrahimov, M. A., et al. 2005, A&A, 438, 39 10 Raiteri, C. M., Villata, M., Larionov, V. M., et al. 2008, A&A, 491, 755 2, 8 Schlafly, E. F., & Finkbeiner, D. P. 2011, ApJ, 737, 103 9 Schleicher, B., Arbet-Engels, A., Baack, D., et al. 2019, Galaxies, 7, 62 3 Sikora, M., Begelman, M. C., & Rees, M. J. 1994, ApJ, 421, 153 1 Stickel, M., Padovani, P., Urry, C. M., Fried, J. W., & Kuehr, H. 1991, ApJ, 374, 431 5 Sun, S.-S., Wang, Z.-X., & Xing, Y. 2023, Research in Astronomy and Astrophysics, 23, 075006 2 Urry, C. M., & Padovani, P. 1995, PASP, 107, 803 1, 5 Villata, M., Raiteri, C. M., Balonek, T. J., et al. 2006, A&A, 453, 817 8 Willingale, R., Starling, R. L. C., Beardmore, A. P., Tanvir, N. R., & O'Brien, P. T. 2013, MNRAS, 431, 394 7 Zhang, B.-K., Zhao, X.-Y., & Wu, Q. 2022, ApJS, 259, 49 2, 5, 7, 8, 9, 10 Zhang, H., Yan, D., & Zhang, L. 2023, ApJ, 944, 103 2\n```  \nAppendix A: COLOR MAGNITUDE DIAGRAMS OF 44 BLAZARS  \nFig. A.1: Color magnitude ( zg -zr versus zr ) diagrams for 24 blazar targets. The blue line indicates the fit to the BWB part and the red line the fit to the stable part. The magnitudes at the turning points, with the uncertainties from the fitting, are given in Table 1. Note that the panels have different scales, which are set to include all the data points for each target. No error regions from the fitting are shown for clarity.  \n<!-- image -->  \nFig. A.2: Same as Figure A.1 for 20 blazar targets.  \n<!-- image -->"}