Fast Optical Variability of the TeV Blazar PKS 1725+123 Observed by SVOM-VT and Insights from Multiwavelength Follow-up Observations

We produced the SVOM-VT light curves of PKS 1725+123 in both the VT_R and VT_B bands with a time bin of 5 minutes, as shown in Figure 1. The source exhibits significant variability on daily timescales during the SVOM-VT observations. The brightness of PKS 1725+123 in both the VT_R and VT_B bands declines significantly from 2025 August 20 to August 23, then increases, and subsequently decreases again on 2025 September 5. The brightest magnitude is 13.62 ± 0.01 mag in the VT_R band and 14.02 ± 0.02 mag in the VT_B band observed on 2025 August 20. The brightness variation of the source is approximately 0.89 mag in the VT_R band and approximately 0.98 mag in the VT_B band from August 20 to August 23.

Zoom In Zoom Out Reset image size

To quantify the significance of variability, a weighted mean flux 〈F〉 is estimated using (M. A. McLaughlin et al. 1996)

where N is the number of data points and F_i and σ_i are the flux and its error for the ith data point. We assume a constant flux of 〈F〉 and subsequently calculate the χ² value by

and the associated probability p(χ²) = 1 − p(>χ²). The results indicate that the significance of the variability in both the VT_R and VT_B bands far exceeds the 5σ confidence level, with the estimated values of p(>χ²) being considerably lower than 5.7 × 10⁻⁷.

To investigate whether PKS 1725+123 exhibits significant variability on hourly timescales, we reproduce the VT light curves obtained from several consecutive orbital observations conducted within a 24 hr period. As shown in Figure 2, there are two observation epochs: consecutive four-orbit observations from August 20 to August 21 and consecutive six-orbit observations from August 22 to August 23. The data analysis uses a time bin of 5 minutes. Using Equations (1) and (2), we estimate the significance of variability for these light curves. The results indicate that the source exhibits significant variability on hourly timescales during both observation epochs in both the VT_R and VT_B bands, with a confidence level far exceeding 5σ. The weighted mean flux 〈F〉 for each observation epoch is also presented as dashed lines in Figure 2.

Zoom In Zoom Out Reset image size

We further assess the significance of variability for each orbit observation, as presented in Table 1. The data analysis also uses a time bin of 5 minutes. The results suggest that the source shows significant variability with a confidence level exceeding 5σ in the VT_B band during the second and fourth orbit observations on August 20–21, as displayed in Figure 2. In contrast, there is no variability during the second orbit observation and variability with a 3.5σ confidence level during the fourth orbit observation in the VT_R band. Evidently, the variability of PKS 1725+123 in the optical band can be as short as several minutes. Nevertheless, no significant (>5σ) variability is observed for the remaining orbit observations.

The time evolution of the color index is also presented in Figures 1 and 2, where the color index is defined as the magnitude in the VT_B band minus that in the VT_R band. It can be observed that the color index, which accompanies the magnitude variation, also exhibits evident variation on daily timescales (Figure 1), with a confidence level exceeding 5σ. Subsequently, we plot the color index versus the magnitudes in both the VT_R and VT_B bands for all the observational data, as displayed in Figures 3(a) and (b). It appears that the data are distributed into two different clusters. Considering the uncertainties associated with both the color index and the magnitude, we use the bootstrap method (B. Efron 1979) to estimate the correlation coefficient (r) between the color index and the magnitude. We obtain r = −0.78 ± 0.05 and r = 0.10 ± 0.17 when the VT_R magnitude values are higher (fainter) and lower (brighter) than 13.8 mag, respectively. Similarly, we obtain r = −0.61 ± 0.08 and r = 0.45 ± 0.11 when the VT_B magnitude values are higher (fainter) and lower (brighter) than 14.2 mag, respectively. These results likely indicate a redder-when-brighter (RWB) trend for the source if it is fainter than 13.8 mag in the VT_R band (or 14.2 mag in the VT_B band).

Zoom In Zoom Out Reset image size

Note that the points with magnitudes brighter than 13.8 mag in the VT_R band (or 14.2 mag in the VT_B band) in Figure 3(a) (or Figure 3(b)) correspond to the data from the first four orbit observations shown in Figure 1 (and to the left panels in Figure 2). Therefore, we further estimate the correlation coefficient between the color index and the magnitude for each of the four orbit observation datasets. The bootstrap analysis yields r = 0.78 ± 0.18 (r = 0.93 ± 0.09 for VT_B) and r = 0.74 ± 0.18 (r = 0.90 ± 0.14 for VT_B) for the VT_R observation data of the second (Figures 3(c) and (d)) and fourth (Figures 3(e) and (f)) orbits, whereas no correlation is observed for the remaining two orbit observations. Clearly, the source exhibits a bluer-when-brighter (BWB) behavior during the second and fourth orbital observations.

The source was not significantly detected in the SVOM-MXT observations. We thus assume an absorbed power-law (PL) spectrum with a Galactic column density of N_H = 8.64 ×10²⁰ cm⁻² (HI4PI Collaboration et al. 2016) and a photon index of Γ_X = 2.0 to compute the 3σ upper limits of the flux in the 0.3–10 keV energy band. To provide a representative constraint across all pointings, we calculated a weighted-average upper limit of 4.7 × 10⁻¹² erg cm⁻² s⁻¹, where the weighting accounts for the exposure time of each observation.

No significant emission was detected by SVOM-ECLAIRs during the observational campaign. Therefore, the 3σ upper limits of the flux were calculated for each individual observation using the variance maps produced by the ECLAIRs pipeline. A 3σ weighted-average upper limit of the flux in the 4–150 keV energy band, which accounts for the exposure duration of each pointing, was finally computed as 5.8 × 10⁻¹⁰ erg cm⁻² s⁻¹, assuming an E⁻² PL spectrum.

We also conduct an analysis of the X-ray observational data of the source from Swift-XRT. The data reduction procedure of the Swift-XRT observations is described in Appendix A.1. We derive the X-ray light curves of PKS 1725+123 in the 0.3–2 keV (F_0.3−2), 2–10 keV (F₂₋₁₀), and 0.3–10 keV (F_0.3−10) bands, along with the temporal curves of the photon index (Γ_X) and hardness ratio (HR; i.e., ), as displayed in Figure 4.

Zoom In Zoom Out Reset image size

The light curves of F_0.3−2, F₂₋₁₀, and F_0.3−10 for PKS 1725+123 exhibit significant variability at a confidence level exceeding 5σ. It should be noted that, based on the Swift-XRT observations, the source on 2025 August 19 and August 21 (two data points within the shaded region) is in a historically high-flux state. This state coincides with the VHE detections and the highest SVOM-VT flux observations.

It is noteworthy that the value of F_0.3−10 seems to increase from 2025 August 19 to 21 (two data points within the shaded region). In contrast, the spectrum becomes softer. Five days later, on 2025 August 25 and 26, the flux declines to approximately half of that on 2025 August 21 but with a harder spectral index. Notably, for almost all observations, the flux in the 0.3–10 keV band is predominantly determined by the value of F₂₋₁₀, with the exception of the two observations that are concurrent with the VHE detections, as shown in the second panel from the top of Figure 4. In other words, the observed decrease in flux from 2025 August 21 to August 25 is mainly attributed to the decline in flux in the soft 0.3–2 kev energy band.

The FSRQ PKS 1725+123 has been detected in the GeV γ-ray band by Fermi-LAT, and it is associated with the γ-ray source 4FGL J1728.0+1216 (J. Ballet et al. 2023). Therefore, we conduct an analysis of the 17 yr Fermi-LAT observational data for this source. The data reduction procedure is described in Appendix A.2.

The 17 yr integrated spectrum of PKS 1725+123 in the 0.1–1000 GeV band requires a log-parabola (LP) model for adequate explanation, with a photon spectral index of Γ_γ =2.16 ± 0.03 and a curvature parameter of β = 0.05 ± 0.01, as presented in Figure 5. We generate the 17 yr Fermi-LAT light curve in the 0.1–1000 GeV band with a time bin of 30 days, as shown in Figure 7(a) in Appendix A.2. The γ-ray emission of the source remains nearly in a low and stable flux state prior to MJD 60353 (2025 February 12), after which the γ-ray flux increases, particularly during 2025 August. To further explore the flux state and the variability features in the GeV band before and after the TeV detection, we produce the light curve with a time bin of 1 day from 2025 July 6 (MJD 60862) to 2025 October 1 (MJD 60949), as shown in Figure 7(b) in Appendix A.2. Nevertheless, no significant variability is detected in the 3 month light curve despite the source being in a high-flux state.

Zoom In Zoom Out Reset image size

Using the gtsrcprob tool, we estimated the maximum energy of photons detected during the entire 17 yr Fermi-LAT observation. The maximum energy of the detected photons is 103 GeV, detected on 2025 August 17 (MJD 60904), with a probability exceeding 99%. To investigate the spectral variation feature of PKS 1725+123 in the GeV band during the high-flux state in the optical band, we perform a time-resolved spectral analysis for the Fermi-LAT observational data in the following time intervals: 2025 August 16–18 (MJD 60903–60905), August 19–21 (MJD 60906–60908), and August 22–24 (MJD 60909–60911). The results of the spectral analysis are also presented in Figure 5. The three time-resolved spectra can be well modeled by a simple PL function with a flatter photon spectral index than that of the 17 yr integrated spectrum. As illustrated in Figure 5, minor flux variations are observed at the high-energy end among the three time-resolved spectra, and the spectrum appears to soften from August 16 to August 24. Specifically, the flux is (2.45 ± 1.77) × 10⁻⁹ erg cm⁻² s⁻¹ with Γ_γ = 1.72 ± 0.12 during August 16–18, (7.47 ± 2.97) × 10⁻¹⁰ erg cm⁻² s⁻¹ with Γ_γ = 1.98 ± 0.14 during August 19–21, and (6.10 ± 1.77) ×10⁻¹⁰ erg cm⁻² s⁻¹ with Γ_γ = 2.15 ± 0.50 during August 22–24. It seems that PKS 1725+123 demonstrates a harder-when-brighter trend in the GeV band.

Swift-XRT has conducted a total of 18 observations of PKS 1725+123. In this work, the X-ray spectral analysis of PKS 1725+123 is mainly based on 16 observations. The data from two observations (ObsIDs 00019638007 and 00019638009) were not utilized because the exposures were insufficient to acquire an adequate number of photons for spectral fitting. The data were processed using the XRT Data Analysis Software (XRTDAS, v.3.7.0) within the HEASoft package (v.6.34). The raw event files were calibrated using the calibration files from Swift-XRT CALDB (v.20241028) through the xrtpipeline task. Spectra of source and background photons were extracted using the xselect task. We defined the source region as a circle with a radius of 47^″ and an annulus with inner and outer radii of 71^″ and 142^″, respectively, both centered on the coordinates of PKS 1725+123 provided in M. H. Xu et al. (2019). Ancillary response files, incorporating cumulative exposure maps to account for point-spread function losses and CCD defects, were generated using the xrtmkarf task.

The spectra of PKS 1725+123 for different Swift-XRT observations were individually fitted with an absorbed PL model within Xspec (v.12.14.1; K. Arnaud et al. 1999), specifically, TBABS×POWERLAW. The PL model is described as follows:

where N₀ is the PL normalization, E₀ = 1 keV is the scale parameter of photon energy, and Γ_X is the photon spectral index. In this model, the TBABS component accounts for the Galactic photoelectric absorption. During the spectral fitting process, the column density (N_H) was fixed at the Galactic value (N_H = 8.64 × 10²⁰ cm⁻²; HI4PI Collaboration et al. 2016) since the data quality of Swift-XRT is inadequate to constrain the intrinsic absorption. To improve the quality of the spectral fitting, the grppha task was utilized to rebin the spectra, ensuring a minimum of 3 counts to guarantee the validity of the C-statistic during the spectral fitting. Moreover, a minimum of 20 counts was applied to use the χ² statistic during the spectral fitting for the observation with sufficient photon counts. The best-fit parameters are summarized in Table 2.

Table 2. Swift-XRT Spectral Analysis Results for PKS 1725+123

ObsID	Date	Exposure	N0	ΓX	Fit Statistic	F0.3−10	F0.3−2	F2−10
	(YYYY-MM-DD)	(s)	(10−4 ph cm−2 s−1 keV−1)			(10−12 erg cm−2 s−1)
00019638001	2025-03-21	1296		1.59 ± 0.25	10/13		0.92  ±  0.14
00019638002	2025-03-24	1486			10/7
00019638004	2025-03-29	438			4/1
00019638005	2025-08-19	1509		2.22 ± 0.13	5/9a
00019638006	2025-08-21	837			4/4a
00019638010	2025-08-25	1294		1.76 ± 0.21	20/22
00019638011	2025-08-26	2202			3/4a
00019638012	2025-09-07	411			8/11
00019638013	2025-09-10	870		1.52 ± 0.24	15/16
00019638014	2025-09-11	865		1.09 ± 0.23	14/16
00019638015	2025-09-12	992			4/4
00019638016	2025-09-13	1075		1.38 ± 0.18	38/27
00019638017	2025-09-14	1075			19/23		1.81 ± 0.21
00019638018	2025-09-28	2316		1.59 ± 0.16	36/32		1.10 ± 0.11	1.91 ± 0.19
00019638019	2025-09-29	1810		1.59 ± 0.20	36/21
00019638020	2025-10-01	1725		1.32 ± 0.19	18/24

Note. ^aThe χ² statistic was employed during the spectral fits of these observations.

Download table as:  ASCIITypeset image

In the Fermi-LAT 14 yr Source Catalog (4FGL-DR4; S. Abdollahi et al. 2022 ; J. Ballet et al. 2023 ), PKS 1725+123 is associated with γ-ray source 4FGL J1728.0+1216. We selected the Pass 8 data within the 0.1–1000 GeV band from a 15^∘ region of interest (ROI) centered on the position of PKS 1725+123 (R.A. = 262029, decl. = 12261). The data cover the period from 2008 August 4 to 2025 October 1 (MJD 54682–60949). The binned maximum-likelihood analysis was conducted using the publicly available software fermitools (ver. 2.2.0; Fermi Science Support Development Team 2019) and the P8R3_SOURCE_V3 instrument response function (IRF), with a bin size of 01. Photons with zenith angles exceeding 90^∘ were excluded to minimize the contamination of γ-ray emission from the Earth limb. All sources within the ROI were included in the model. The spectral parameters of the sources located within a circle with a radius of 6^∘ were left free, while the other parameters were fixed at their 4FGL-DR4 values. The background emission was modeled using the diffuse Galactic interstellar emission (“gll_iem_v07. fits”) and isotropic emission (“iso_P8R3_SOURCE_V3_V1. txt”), and only their normalization parameters were kept free.

We use the maximum test statistic (TS) to evaluate the significance of the γ-ray signals of a source, where , where and are the likelihood values for the background with and without a source. The threshold is set at TS = 25 (corresponding to 5σ). The appearance of new sources in the ROI may affect the γ-ray spectrum and light curve of the research source. To eliminate the influence of new sources, we searched for possible new γ-ray sources within a 5^∘ field of view and found no new sources beyond those listed in the 4FGL-DR4.

Through maximum-likelihood fitting, the value of TS=2495.5 was obtained for 4FGL J1728.0+1216. The 17 yr average γ-ray spectrum can be well explained in the LP form (E. Massaro et al. 2004), specifically,

where N(E) is the photon distribution as a function of energy, Γ_γ is the photon spectral index, E_b is the scale parameter of photon energy, and β is the curvature parameter. An average spectrum with Γ_γ = 2.16 ± 0.03 and β = 0.05 ± 0.01 was obtained, and the 17 yr average flux is F_{0.1−1000 GeV} =(2.12 ± 0.10) × 10⁻¹¹ erg cm⁻² s⁻¹. These results are roughly consistent with those of 4FGL-DR4.

Using the same method, we obtained the time-resolved spectra for the three time intervals: 2025 August 16–18 (MJD 60903–60905), August 19–21 (MJD 60906–60908), and August 22–24 (MJD 60909–60911). Different from the 17 yr integrated spectrum, the three time-resolved spectra can be well modeled by a simple PL function with a flatter photon spectral index.

We generated the 17 yr long-term light curve using a time bin of 30 days, as displayed in Figure 7(a). During the construction of the light curves, only the parameters of the target source and two background files were released, while the other parameters were fixed at the optimal fitting values. If TS < 9, an upper limit of flux was provided for that time bin. We fixed the target source in the form of a PL function with a photon spectral index of 2 to calculate the upper-limit value. Using the same method, we also generated the 3 month (2025 July 6 to 2025 October 1) light curve using a time bin of 1 day, as depicted in Figure 7(b).

Zoom In Zoom Out Reset image size

As described in Section 3.1, SVOM-VT has observed the optical variability of PKS 1725+123 with a timescale of less than 1 hr, which suggests a small size of the emission region. Meanwhile, the hard optical spectrum, in combination with a soft X-ray spectrum, indicates a very high synchrotron peak frequency. Further considering a very hard Fermi-LAT spectrum and a high redshift, a typical one-zone leptonic model fails to provide an adequate explanation for the broadband SED of PKS 1725+123. Nevertheless, during the optical flare state observed by SVOM-VT, the emission in the X-ray and γ-ray bands also demonstrates a high-flux state, indicating that the emission among the three bands has some kind of connection. In other words, neither the model with two independent zones nor the hadronic model is suitable (e.g., S. C. Chaudhary et al. 2026). Additionally, the detection of VHE emission require extremely high-energy electrons to produce it. It is noted that both radio observations and numerical simulations have indicated that the jets possess a spine–sheath velocity stratification structure (e.g., P. Rossi et al. 2008; A. Chhotray et al. 2017; R. C. Walker et al. 2018). When electrons pass through the shear layers associated with velocity gradients, the shear acceleration process, a Fermi-type mechanism similar to shock or stochastic acceleration, can produce a population of higher-energy electrons (F. M. Rieger & P. Duffy 2004; F. M. Rieger et al. 2007; R.-Y. Liu et al. 2017; F. M. Rieger 2019). Therefore, a two-zone leptonic model, as described in Section 4, is proposed.

The electron distribution in both emission regions is taken as a PL, characterized by an electron number density parameter N₀, a spectral index p, and a minimum () and a maximum () Lorentz factor of electrons. Both radiation regions are assumed to be spherical with a radius R, a magnetic field strength B, and a Doppler boosting factor δ, where , and Γ and θ are the bulk Lorenz factor and viewing angle of the emission region. The value of R is constrained by the variability timescale, specifically, R = Δtcδ/(1 + z). As displayed in Figure 2, a continuous decrease in flux was observed during the fourth orbit observation (a similar phenomenon was also presented in the VT_B band of the second orbit observation), with a duration of 30 minutes from the first to the last data points. Therefore, a Δt value of 30 minutes is adopted for the compact region, and it is assumed that δ = Γ. For the extended region, a Δt value of 1 day is adopted, and the same viewing angle as that of the compact region is assumed. Given that the VHE γ-rays are detected during this high-flux state, it is assumed that the emission region is not located deep within the BLRs, as this would result in strong attenuation through γγ annihilation (e.g., H. T. Liu & J. M. Bai 2006; J. Poutanen & B. Stern 2010). Therefore, the energy density of the BLR is fixed at 0.5 times the typical value of 2.73 ×10⁻² Γ² erg cm⁻³, and the spectrum of the BLR can be approximated by a blackbody with a peak at 2 × 10¹⁵ Γ Hz in the comoving frame (e.g., F. Tavecchio & G. Ghisellini 2008b; S.-J. Kang et al. 2014). The electron radiative processes have incorporated synchrotron self-absorption (G. B. Rybicki & A. P. Lightman 1979 ), the Klein–Nishina effect (F. Tavecchio et al. 1998), and EBL absorption (J. D. Finke et al. 2022).

Since it is assumed that the two regions overlap, the code also incorporates the emission resulting from the interaction of the two zones. The synchrotron radiation from the two regions provides an additional target photon field for each other for the IC scattering. Moreover, due to the relative motion between the compact region and the extended region, the energy density of the seed photons will be amplified (G. Ghisellini et al. 2005; Z.-R. Wang et al. 2024). However, the radiation contribution from the interaction of the two zones is mainly attributed to the synchrotron photons from the compact region undergoing IC scattering by the relativistic electrons within the extended region. In contrast, the flux contribution from the opposite process can be disregarded.

The fitting results are presented in Figure 6, and the fitting parameters are given in Table 3. It should be noted that the derived parameter values are based on visual assessments, and the model parameters cannot be fully constrained by the current observational data.

Table 3. Parameters of SED Fitting and Jet Powers

Parameter	Symbol	Compact Zone	Extended Zone
Electron density parameter	N0 [cm−3]	6 × 102	1 × 102
Electron spectral index	p	1.4	2
Minimum Lorentz factor for electrons		1 × 102	3 × 103
Maximum Lorentz factor for electrons		4.5 × 103	5 × 105
Radius of emission region	R [cm]	1.2 × 1015	1.9 × 1016
Bulk Lorentz factor	Γ	35	5.9
Doppler factor	δ	35	11.4
Magnetic field	B [G]	24	1
Energy density of BLR	UBLR [erg cm−3]	1.36 × 10−2	1.36 × 10−2
Equipartition ratio	UB/Ue	2 × 102	9.5 × 101
Radiative power	Pr [erg s−1]	8.58 × 1041	3.02 × 1043
Nonthermal electron power	Pe [erg s−1]	1.87 × 1043	1.77 × 1042
Magnetic field power	PB [erg s−1]	3.75 × 1045	1.69 × 1044
Jet power	Pjet [erg s−1]	3.77 × 1045	2.01 × 1044

Download table as:  ASCIITypeset image

Given the rapid variability observed in the optical band, it was assumed that there is a compact emission region during the SED modeling. The synchrotron frequency (ν_s), of an electron with a Lorentz factor of is (γ_e) is (G. B. Rybicki & A. P. Lightman 1979; F. Tavecchio et al. 1998). The cooling timescale of electrons through the synchrotron process is estimated by (G. B. Rybicki & A. P. Lightman 1979; M. Chiaberge & G. Ghisellini 1999), where m_e is the electron mass and σ_T is the Thomson cross section. For the contribution to the optical emission, the cooling timescale of electrons is about 46 minutes in the comoving frame and several minutes in the observer frame.

We also estimate the shear-acceleration timescale of electrons in the extended zone under the strong-scattering regime. The acceleration timescale is given by the following equations (R.-Y. Liu et al. 2017):

where Γ_jet = Γ = 35 denotes the bulk Lorentz factor of the jet, β_jet represents the speed of the jet in units of c, and r_g is the gyroradius of electrons. In this scenario, we adopt ξ ∼ 0.1, which represents the ratio of turbulent to ordered magnetic field energy density. We employ the Kolmogorov turbulence model (q = 5/3; A. Kolmogorov 1941) and assume an acceleration region width Δr ∼ 10¹⁵ cm with . Then, we obtain t_acc ∼ 4 × 10³ s.

It should be noted that the estimation results for the cooling and acceleration timescales of electrons are consistent with the observations.

By assuming that the the jet power is carried by relativistic electrons (P_e), magnetic fields (P_B), and radiation (P_r), i.e., P_jet = πR²Γ² c(U_e + U_B + U_r), and using the SED fitting parameters, we calculate the jet power and the powers of each component for both the compact and extended regions. U_e, U_B, and U_r are the energy densities of relativistic electrons, magnetic fields, and radiation in the comoving frame, respectively, which are given by the following equations (e.g., A. Celotti & G. Ghisellini 2008; G. Ghisellini et al. 2010):

where L_bol is the bolometric luminosity.

The derived jet power and the powers of each component for the compact and extended regions are presented in Table 3. It is found that for the compact zone, P_B/P_jet ∼ 0.995 and P_r/P_jet ∼ 0.0002, and for the extended zone, P_B/P_jet ∼ 0.841 and P_r/P_jet ∼ 0.150. These results suggest that both zones are highly magnetized; however, the compact zone exhibits low radiation efficiency, while the extended zone demonstrates high radiation efficiency. We compare the derived jet powers of PKS 1725+123 with a γ-ray emission FSRQ sample from A. Celotti & G. Ghisellini (2008) and J. Zhang et al. (2020b), and we plot P_e, P_B, and P_r as functions of P_jet in Figure 8. We combine the compact (red filled square) and extended (blue filled square) zones as a single jet emission region of PKS 1725+123, which is marked as a green filled star in Figure 8. Comparing with other γ-ray emission FSRQs, on average, PKS 1725+123 exhibits a low jet power and radiation efficiency but a highly magnetized jet.

Zoom In Zoom Out Reset image size