# X-Ray Observations of 1ES 1959+650 in Its High-activity State in 2016–2017 with AstroSat and Swift

Sunil Chandra<fnref n="1" /> ID, Markus Boettcher<fnref n="1" /> ID, Pranjupriya Goswami<fnref n="2" /> ID, Kulinder Pal Singh<fnref n="3" />,<fnref n="4" /> ID, Michael Zacharias<fnref n="1" />,<fnref n="5" /> ID, Navpreet Kaur<fnref n="6" /> ID, Sudip Bhattacharyya<fnref n="4" /> ID, Shashikiran Ganesh<fnref n="7" /> ID, and Daniela Dorner<fnref n="8" /> ID

*Received 2021 March 27; revised 2021 May 8; accepted 2021 May 11; published 2021 September 8*

## Abstract

We present a comprehensive multifrequency study of the high-frequency-peaked BL Lac object 1ES 1959+650 using data from various facilities during the period 2016–2017, including X-ray data from AstroSat and Swift during the historically high X-ray flux state of the source observed until 2021 February. The unprecedented quality of X-ray data from high-cadence monitoring with AstroSat during 2016–2017 enables us to establish a detailed description of X-ray flares in 1ES 1959+650. The synchrotron peak shifts significantly between different flux states, in a manner consistent with a geometric (changing Doppler factor) interpretation. A time-dependent leptonic diffusive shock acceleration and radiation transfer model is used to reproduce the spectral energy distributions and X-ray light curves, to provide insight into the particle acceleration during the major activity periods observed in 2016 and 2017. The extensive data of Swift-XRT from 2015 December to 2021 February (exp. = 411.3 ks) reveals a positive correlation between flux and peak position.

*Unified Astronomy Thesaurus concepts: Blazars (164); Active galactic nuclei (16); Black hole physics (159); Jets (870); High energy astrophysics (739)*

*Supporting material: machine-readable tables*

## 1. Introduction

Blazars are a subclass of active galactic nuclei with a jet of relativistic plasma streaming along or very close to the line of sight (Blandford & Rees 1978). The observed electromagnetic (EM) emission, which is predominantly nonthermal in nature, is considered to be emanating from the relativistic jet. The observed features of blazars include superluminal motion of radio jet components, high optical polarization, and strong continuum emission variable at timescales ranging from a few minutes to years, across the entire EM spectrum.

The broadband spectral energy distribution (SED; $\nu$ vs. $\nu F_\nu$ plot) of blazars consists of two distinct broad continuum hump-like structures, with the first one peaking somewhere in submillimeter to soft X-rays, whereas the second one peaks at MeV to TeV energies (Urry & Padovani 1995). The low-energy component of the SED is mostly due to the synchrotron emission from relativistic electrons/positrons gyrating around the magnetic field in the relativistic jet. This emission component, in some cases, is superimposed by significant thermal contributions from, e.g., the accretion disk, a hot corona accompanying the accretion disk, and/or an obscuring dusty torus. On the other hand, the physical mechanisms behind high-energy emission (MeV to TeV) are not well established, and two families of models, namely, (a) leptonic models and (b) hadronic models, both appear to be viable mechanisms to explain the X-ray through $\gamma$-ray emission. In hadronic models the high-energy emission is produced by relativistic protons through proton synchrotron radiation and photo-pion production, followed by pion decay

and electromagnetic cascades (e.g., Mannheim & Biermann 1992; Mannheim 1993; Nellen et al. 1993; Protheroe & Mücke 2001; Mücke et al. 2003; Böttcher et al. 2013). Leptonic models assume the jet protons to be cold enough not to contribute to the radiative output, and high-energy emission is produced by inverse Compton scattering of low-energy seed photons by the ultrarelativistic leptons ($e^-/e^+$). The seed photons may come from the synchrotron radiation field in the emission region, which are up-scattered by the same leptons that produced the synchrotron radiation (synchrotron self-Compton [SSC]; Ghisellini & Maraschi 1989; Bloom & Marscher 1996; Mastichiadis & Kirk 1997). Alternatively, if the seed photons originate external to the emission region (e.g., from the accretion disk, the dusty torus, or the broad-line region), then the process is termed as external Compton (Dermer et al. 1992; Ghisellini et al. 1998).

1ES 1959+650 ($z$ = 0.048; Perlman et al. 1996) is a prominent high-synchrotron-peaked blazar. It was first detected in X-rays during the Slew Survey with the Einstein Imaging Proportional Counter (Elvis et al. 1992), followed by BeppoSAX (Beckmann et al. 2002), RXTE, Swift, and XMM-Newton (Tagliaferri et al. 2003; Massaro et al. 2008) in later years. This source is also a prominent TeV $\gamma$-ray emitter, with the first detection at TeV energies, reported by the Utah Seven-Telescope Array collaboration in 1998 (Nishiyama 1999).

The historical observations establish 1ES 1959+650 to be a high-frequency-peaked BL Lac object (HBL) in which the synchrotron peak of the broadband SED appears in the UV to X-ray band (Krawczynski et al. 2004; Abdo et al. 2010;

<footnote n="1">¹Center for Space Research, North-West University, Potchefstroom, 2520, South Africa; sunil.chandra355@gmail.com, sunil.chandra@nwu.ac.in</footnote>

<footnote n="2">²Department of Physics, Tezpur University, Assam, India</footnote>

<footnote n="3">³Indian Institute of Science Education and Research, Mohali, Punjab, India</footnote>

<footnote n="4">⁴Department of Astronomy and Astrophysics, Tata Institute of Fundamental Research, 1 Homi Bhabha Road, Mumbai 400005, India</footnote>

<footnote n="5">⁵Laboratoire Univers et Théories, Observatoire de Paris, Université PSL, CNRS, Université de Paris, F-92190 Meudon, France</footnote>

<footnote n="6">⁶KTH Royal Institute of Technology, Sweden</footnote>

<footnote n="7">⁷Physical Research Laboratory Ahmedabad, India</footnote>

<footnote n="8">⁸Department of Physics and Astronomy, Julius Maximilian University Würzburg, Am Hubland, D-97074 Würzburg, Germany</footnote>
Kapanadze et al. 2016a). This source exhibits strong flux variability across almost the entire EM spectrum. The flux increase of up to 3–4 orders of magnitude in the optical, X-ray, and TeV energy bands during the short/erratic flares has been witnessed for 1ES 1959+650 (Krawczynski et al. 2004; Perlman et al. 2005; Kapanadze et al. 2016b). The rapid variability and its frequency dependence provide crucial insight into the physical processes of particle acceleration and radiation mechanisms, as well as the geometry and size of the emission region (e.g., Dondi & Ghisellini 1995).

Recent high-sensitivity X-ray observations have found several high flux states and strong X-ray outbursts of this source. XMM-Newton and RXTE-PCA observations in 2002–2003 revealed strong X-ray flares with flux variations by a factor of up to ~4.2 (Krawczynski et al. 2004; Perlman et al. 2005). Many of such frequently occurring strong X-ray flares were reported by Kapanadze et al. (2016b) during 2006–2014 using Swift-XRT observations. The source underwent a number of active states and an unprecedented X-ray flaring activity during 2015 August–2016 January that was observed by Swift-XRT. The observed count rate was reported to vary by a factor of ~5.7 with a maximum value above 20 counts s⁻¹, with simultaneous high flux activity in the TeV energy band (Kapanadze et al. 2016a; Kaur et al. 2017; Patel et al. 2018; MAGIC Collaboration et al. 2020). During this large flare, the synchrotron peak position of the SED showed a tendency to shift toward the higher X-ray energies accompanied by a hard X-ray spectral index. The detailed X-ray spectral studies further confirmed the harder-when-brighter trend (Tagliaferri et al. 2003; MAGIC Collaboration et al. 2020). However, during most of these epochs with an X-ray flare, the TeV counterpart was found to be in low flux states. On the other hand, in several multiwavelength campaigns, “orphan” flares in very high energy (VHE; used for TeV) γ-rays (not accompanied by a simultaneous X-ray flare) were reported in 2002 June (Krawczynski et al. 2004) and 2012 April—June (Aliu et al. 2014).

The uncorrelated variability is inconsistent with the simplest one-zone SSC models, which are often successful in reproducing the broadband emission of HBLs but have proven to be inadequate to explain several aspects of emission in many studies (Krawczynski et al. 2004; Patel et al. 2018; MAGIC Collaboration et al. 2020). “Orphan” flares hint at a more complex geometry and/or underlying particle distribution, such as those invoked in multiple-component SSC models and/or external Compton models (e.g., the synchrotron mirror model, Böttcher & Dermer 1999) within the leptonic schemes. Recently, Shah et al. (2021) have reported an anticorrelation between the photon index and X-ray flux using a broken power law for analyzing only a segment of X-ray data presented here.

In this work, we present a detailed investigation of the X-ray spectral and light-curve features of 1ES 1959+650 observed by AstroSat in 2016–2017. Our main focus here is to understand the distinct, irregular X-ray outbursts observed during this period by both the Soft X-Ray Telescope (SXT) and Large Area Proportional Counter (LAXPC) instruments aboard AstroSat.

These data are supplemented with simultaneous/quasi-simultaneous XRT data extending before and after 2016–2017, spanning over 6 yr from 2015 January 30 to 2021 February 09, and also other multiwavelength data to probe the evolution of the underlying nonthermal particle distribution. In order to consistently

fit the SEDs and the light curves obtained during these erratic flares observed with AstroSat, we adopt the time-dependent multizone shock acceleration and radiation transfer model, as described by Böttcher & Barknmh (2019), to investigate the nature of shocks responsible for the observed spectral variability. We further provide a detailed analysis of the time-resolved spectra and light curves and their correlation over the span of ~6 yr.

The paper is structured as follows. Section 2 describes the multiwavelength observation details and the data analysis procedures. In Section 3, we provide the results of the timing and spectral analysis and the detailed correlation study. Section 3 contains the results and the interpretation through the modeling of snapshot SEDs and light curves. In Section 4 we summarize our work, followed by a comprehensive discussion.

## 2. Observations and Data Analysis

1ES 1959+650 was observed in campaign mode during 2016 and 2017 at various epochs, representing different flux states, using a number of observing facilities, including AstroSat, Swift, and Mt. Abu Infrared Observatory, India (MIRO). The details of the observing epochs and the respective total exposure times are mentioned in the Table 4. PASS8 photon data from Fermi-LAT are also analyzed to study the high-energy (GeV) emission. The following subsections provide the details of the observations and analysis procedures.

### 2.1. Fermi-LAT

The PASS8 (P8R3) Fermi-LAT photon data and corresponding spacecraft data from the beginning of 2015 November to the end of 2017 December are downloaded from the LAT data center<fnref n="9" /> with a search radius of 30° and in an energy range of 30 MeV to 500 GeV. The Fermitools package (version 1.2.1 conda-release) distributed by the Fermi Science Support Center, installed with the most recent release of point-source (4FGL) and extended-source catalogs, is used to analyze the data. The python package, *fermipy*<fnref n="10" /> (Wood et al. 2017) is used, which facilitates handy wrappers for various procedures of LAT data analysis, as described by the instrument teams, including model optimization, the localization, sanity checks, product extractions, etc. The initial selection of parameters includes a bin size of 0.1 pixels for map creation, a zenith angle of accepted events of 90° to exclude or eliminate most of the contamination from secondary γ-rays contributed by Earth's limb, an energy range 100 MeV–500 GeV, event type 3, and event class 128. The P8R3_SOURCE_V2 instrument response functions (IRFs) are used. The initial source model (XML) is created by including all the point-like and extended sources located within 25° radius of the location of 1ES 1959+650 as listed in the Fourth Fermi-LAT Source catalog (4FGL; Abdollahi et al. 2020), as well as the Galactic diffuse (gll_iem_v07.fits) and isotropic background emission (iso_-P8R3_SOURCE_V2_v1.txt). The source model for 1ES 1959+650 imported from the 4FGL catalog is “Log Parabola”; however, due to poor photon statistics for the duration of interest to this work, a re-optimization of the source model is performed after forcing the spectral shape of 1ES 1959+650 to “PowerLaw”. The spectral parameters of sources within 5° of 1ES 1959+650 are kept variable, while others are kept frozen

<footnote n="9">⁹ https://fermi.gsfc.nasa.gov/cgi-bin/ssc/LAT/LATDataQuery.cgi</footnote>

<footnote n="10">¹⁰ http://fermipy-readthedocs.io/en/latest/</footnote>
<table><caption>**Table 1**
Details of the Data Used for the Present Study</caption><thead><tr><th>S.N.</th><th>Instrument</th><th>Total Exposure</th><th>Epoch of Observations</th></tr></thead><tbody><tr><td>1</td><td>Fermi-LAT</td><td>...</td><td>57037.0 to 58072.0</td></tr><tr><td>2</td><td>AstroSat</td><td>143.9 ks</td><td>57666.2-57666.7<br>57708.4-57709.2<br>57695.9-57699.4<br>58051.1-58052.6</td></tr><tr><td>4</td><td>Swift</td><td>411.2 ks</td><td>57052.1-59254.4</td></tr><tr><td>5</td><td>MIRO</td><td>...</td><td>57690.81-57696.85</td></tr><tr><td></td><td></td><td>R band: 8.4 hr</td><td>”</td></tr><tr><td></td><td></td><td>B band: 0.65 hr</td><td>”</td></tr><tr><td></td><td></td><td>V band: 0.65 hr</td><td>”</td></tr><tr><td>6</td><td>FACT</td><td>12.4 days</td><td>57632.88-57719.83<br>57997.88-58098.83</td></tr></tbody></table>

Note. MIRO: Mt. Abu Infrared Observatory, Rajasthan, India; FACT: First G-APD Cerenkov Telescope, La Palma, Spain https://fact-project.org/.

to their best-fit values from the catalog. TS (Test Statistics) maps and diffuse maps are generated to look for any possible GeV source (point and/or diffuse) not included in our model, but none are found. Once the model is optimized, the best-fit spectral parameters for the GeV part of the SED are estimated using the *sed* procedure of *fermipy* with two spectral points per decade in energies. The *lightcurve* procedure of the *fermipy* package is used for generating light curves with 1-, 2-, and 3-day binning.

The SED and light-curve data points with $TS \ge 9$ (equivalent to $\ge 3\sigma$ significance) and $TS \ge 25$ ($\ge 5\sigma$) are used for spectral and temporal studies. For lower-significance points, 95% upper limits are shown.

## 2.2. X-ray Data Analysis

The data from extensive monitoring over the course of 2 yr using AstroSat are used for this study. Complementary data from the Neil Gehrels Swift Observatory are also used for various epochs. Table 2 lists the details of the data used for this work. The general FTools and several mission-specific tools distributed as part of the HEAsoft package (version 6.25) and the most recent calibration database<fnref n="11" /> are utilized as appropriate to analyze data from various facilities.

### 2.2.1. AstroSat Soft X-ray Telescope (SXT)

The SXT aboard AstroSat is a 2 m approximate Wolter-I type focusing instrument sensitive mainly in the 0.3–7.1 keV energy band (Singh et al. 2014, 2016, 2017). Its camera assembly uses an e2v CCD, identical to the one flown with XMM-Newton-MOS and Swift-XRT, at its focal plane as the main detector system. The observations were carried out in photon counting mode. The source was observed throughout all the satellite orbits when the SXT was pointed at it, taking care that the Sun avoidance angle is $\ge 45^\circ$ and the RAM angle (the angle between the payload axis and the velocity vector direction of the spacecraft) is $> 12^\circ$ to ensure the safety of the mirrors and the detector.

Level-2 data provided by the SXT payload operation center (POC) in Mumbai, India, are reduced using the most recent pipeline and calibration database (version 1.4b).

The level-2 cleaned event files are used in the `XSELECT` tool distributed with HEAsoft to extract source light curves, images, and spectra. The clean process removed events during the occultation by Earth, any contamination by the charged particles due to passage of the satellite through the South Atlantic Anomaly region, and selected events with grade 0–12 (single-quadruple events). The filter “pha_cutoff” is applied to select various energy channels (e.g., channels 30–700 for the 0.3–7.0 keV band) for the light curves.

The pileup effect is very common for the CCD-based detectors, however, due to the large point-spread function (PSF) of SXT, it becomes effective only for extremely bright sources like Crab or brighter (source count rate $> 180 s^{-1}$). 1ES 1959+650, being only moderately bright (max SXT count rate $\sim 25 s^{-1}$; see Figure 10), has no detectable pileup issue. A circular region of 16' radius around the source location, which encircles more than 95% of all photons as estimated by the standard `sxtEEFmake` module distributed through the POC website, is used to extract the source spectra and light curves. The appropriate ARF file suitable for the specific source region is generated using the command line auxiliary tool `sxtARFModule`. Because of the large PSF of SXT, we are unable to extract background products from the same frame, and hence the background spectrum provided by the POC<fnref n="12" /> is used for spectral modeling. The background correction in the light curves corresponding to various energy bands is done by subtracting a constant that is the rate of background counts for a specific energy band normalized to the area of the source region, estimated by importing the background spectrum in `XSPEC` (version 12.10.1) and applying the corresponding energy filter. This background spectrum is extracted using the data from various blank-sky observations at various locations in the sky taken over the first 3 yr of AstroSat operations. The light curves thus generated are rebinned, and hardness ratio plots are generated using the `lcurve` utility of the HEAsoft package. The spectral modeling of the SXT spectra is performed using `XSPEC`. The $n_H$ column density is fixed to $1.07 \times 10^{21} \text{ cm}^{-2}$ as estimated by the web-based tool of the Leiden/Argentine/Bonn (LAB) Survey of Galactic HI in the direction of 1ES 1959+650 (beam size of $0.266^\circ$) throughout this work to account for the Galactic photoelectric absorption in the Milky Way. Temporally resolved SXT spectra are extracted and modeled to study the temporal variability of the spectral parameters (see Section 2.2.4 for more details).

### 2.2.2. AstroSat Large Area Proportional Counters (LAXPC)

AstroSat hosts three identical units of proportional counters, filled with highly pressurized xenon gas, in a specific arrangement to provide a collective effective area of 6000 $\text{cm}^2$. This instrument is nonfocusing and has a field of view of $\sim 1^\circ \times \sim 1^\circ$. It is sensitive mainly in the 3.0–80.0 keV band (Yadav et al. 2016; Antia et al. 2017).

The LAXPC field-of-view axis is nearly coincident to the other on-axis instruments on board AstroSat, namely, CZTI, SXT, and UVIT. Thus, all sources observed with the SXT as the prime instrument are automatically observed by the hard X-ray detector LAXPC. The AstroSat-LAXPC observations of 1ES 1959+650 performed at various epochs are analyzed using the recent LAXPC pipeline package `laxpcSoft` managed

<footnote n="11">¹¹ https://heasarc.gsfc.nasa.gov/docs/heasarc/caldb/</footnote>

<footnote n="12">¹² https://www.tifr.res.in/astrosat_sxt/dataanalysis.html</footnote>
<table><caption>**Table 2**
Swift Data and the Best-fit Parameter Corresponding to Model M1</caption><thead><tr><th>ObsID</th><th>Time (MJD)</th><th>Exp. (s)</th><th>α</th><th>β</th><th>χ<sup>2</sup>/dof</th><th>F<sub>X, 0.3 - 7keV</sub> (×10<sup>-10</sup> erg cm<sup>-2</sup> s<sup>-1</sup>)</th><th>E<sub>soft</sub> (keV)</th></tr></thead><tbody><tr><td>00035025010</td><td>53884.16</td><td>3280.4</td><td>2.01 ± 0.01</td><td>0.49 ± 0.03</td><td>336.6/307</td><td>3.26 ± 0.02</td><td>0.97 ± 0.03</td></tr><tr><td>00035025118</td><td>56856.74</td><td>1276.3</td><td>1.83 ± 0.04</td><td>0.72 ± 0.09</td><td>105.1/115</td><td>1.32 ± 0.03</td><td>1.32 ± 0.09</td></tr><tr><td>00035025132</td><td>57052.09</td><td>1881.8</td><td>1.94 ± 0.03</td><td>0.52 ± 0.07</td><td>171.5/163</td><td>1.29 ± 0.02</td><td>1.15 ± 0.08</td></tr><tr><td>00035025133</td><td>57068.11</td><td>772.1</td><td>1.91 ± 0.04</td><td>0.78 ± 0.09</td><td>116.2/104</td><td>1.94 ± 0.04</td><td>1.15 ± 0.07</td></tr><tr><td>00035025134</td><td>57071.54</td><td>1322.6</td><td>1.8 ± 0.03</td><td>0.87 ± 0.08</td><td>167.5/158</td><td>1.64 ± 0.03</td><td>1.3 ± 0.06</td></tr><tr><td>00035025135</td><td>57073.91</td><td>1437.9</td><td>2.0 ± 0.03</td><td>0.63 ± 0.07</td><td>182.6/165</td><td>1.8 ± 0.03</td><td>1.0 ± 0.05</td></tr><tr><td>00035025136</td><td>57080.16</td><td>1694.1</td><td>2.08 ± 0.03</td><td>0.57 ± 0.07</td><td>180.8/164</td><td>1.47 ± 0.02</td><td>0.86 ± 0.05</td></tr><tr><td>00035025137</td><td>57087.9</td><td>1335.2</td><td>1.86 ± 0.02</td><td>0.64 ± 0.06</td><td>226.8/211</td><td>3.05 ± 0.04</td><td>1.29 ± 0.06</td></tr><tr><td>00035025139</td><td>57091.49</td><td>623.7</td><td>1.94 ± 0.03</td><td>0.62 ± 0.07</td><td>177.4/156</td><td>3.76 ± 0.07</td><td>1.12 ± 0.07</td></tr><tr><td>00035025141</td><td>57097.51</td><td>1919.8</td><td>2.06 ± 0.02</td><td>0.61 ± 0.06</td><td>215.9/194</td><td>3.36 ± 0.05</td><td>0.9 ± 0.04</td></tr><tr><td>00035025142</td><td>57099.47</td><td>1605.1</td><td>2.06 ± 0.03</td><td>0.62 ± 0.06</td><td>148.8/176</td><td>1.92 ± 0.03</td><td>0.9 ± 0.04</td></tr></tbody></table>

(This table is available in its entirety in machine-readable form.)

and distributed by the LAXPC POC<fnref n="13" /> in Mumbai, India. The background models, response functions, and gain variations are appropriately applied to generate the multiband light curves and spectra. The modeled background spectrum is properly shifted for the gain values appropriate for the time of observations using the command line utility `gainshift`, distributed with `laxpcSoft`. The resulting spectra and light curves are then used for further investigations utilizing `XSPEC` and `lcurve`. The details of observations are listed in Table 4.

### 2.2.3. Neil Gehrels Swift Observatory X-ray Telescope (XRT)

Swift has performed a number of observations covering the duration of interest to this work, sometimes even overlapping with the AstroSat pointings (see Table 2 for details of Swift observations analyzed). The X-ray data from XRT were reprocessed with the mission-specific HEAsoft tool `xrtpi-peline` (version 0.13.4) with standard input parameters as recommended by the instrument team. This step generates new cleaned event files with the most recent calibrations. The event files thus generated are used in the multimission tool `XSELECT` for extracting source and background products. We have analyzed the data taken in both operational modes, namely, the PC and WT modes. For PC mode data an annular region centered at $(\alpha = 19:59:59.929, \delta = +65:08:54.65)$, with an inner radius of 10" and outer radius of 70", is used as the source region, whereas, for the background, another annular region centered at the same location but with inner and outer radii of 150" and 350", respectively, is used. The choice of an annular source region for the PC mode data is made to mitigate the pileup effect because all the observations in this mode show count rates >0.6 counts s⁻¹. We cross-checked for the presence of another X-ray source contaminating the source or background regions. The choice of the source region for the WT mode observations is made following the recommendations by the instrument team.<fnref n="14" /> A circular region of radius 64" centered at the location of 1ES 1959+650 is used for extracting source products. An annular region centered at 1ES 1959+650 with inner and outer radii of 188".59 and 282".88, respectively, is used for the extraction of the background. This background region selection ensures symmetrical placing about 100 pixels

(the half-width of the WT window), and hence no matter where the source is in the WT window, the background region will contain $r_2 - r_1 - 1$ (where $r_1$ and $r_2$ are the inner and outer radii of the annular region) pixels in 1D (minus 1, as the end-of-window pixels are flagged as bad by the ground software processing and are therefore not available for use<fnref n="15" />). The spectra and multiband light curves thus generated are used in `XSPEC` and `lcurve` for high-end investigations. It should be noted that the `BACKSCAL` keywords of the source and background spectra are edited to proper values, as applicable to the current source and background source selections, before performing the spectral analysis.

### 2.2.4. Short-term Time-resolved Spectra from XRT + SXT

Around 400 Swift-XRT spectra taken between 2015.0 and 2021.2 are analyzed using the absorbed Log Parabola photon spectrum model to obtain the spectral changes during various flux states of the source. Additionally, the AstroSat-SXT observations are split into segments of ~3500 s duration in order to generate time-resolved spectra, which are then fitted with the aforementioned model (see Section 3.2 for the model). The time-resolved spectra from SXT are extracted by applying time filtering in `XSELECT` using different merged cleaned event files (one merged file for each individual observation). The total SXT observations thus yield 81 spectra with exposure times ≥1500 s. A similar splitting of the LAXPC observations in such small time bins results in poor spectral data and hence cannot constrain the spectral shapes beyond 8.0 keV. Therefore, the LAXPC spectra are not used for this part of the study. The best-fit XRT model parameters are obtained in the 0.3–10.0 keV band, whereas the unabsorbed fluxes are estimated for the common energy band, i.e., 0.3–7.0 keV, in order to combine the flux estimations from the two X-ray instruments. Tables 2 and 3, available only as supplementary material, provide the details of the best-fit parameters.

## 2.3. Optical/UV Observations

### 2.3.1. Neil Gehrels Swift Observatory UVOT

The UVOT observations in six optical/UV filters for all the relevant observations listed in Table 2 are analyzed using recent

<footnote n="13">¹³ https://www.tifr.res.in/astrosat_laxpc/software.html</footnote>

<footnote n="14">¹⁴ https://www.swift.ac.uk/analysis/xrt/index.php</footnote>

<footnote n="15">¹⁵ https://www.swift.ac.uk/analysis/xrt/backscal.php</footnote>
<table><caption>**Table 3**

AstroSat Data and the Best-fit Parameter Corresponding to Model M1</caption><thead><tr><th>ObsID</th><th>Time (MJD)</th><th>Exp. (s)</th><th>α</th><th>β</th><th>χ<sup>2</sup>/dof</th><th>F<sub>X, 0.3 - 7keV</sub> (×10<sup>-10</sup> erg cm<sup>-2</sup> s<sup>-1</sup>)</th><th>E<sub>exp</sub> (keV)</th></tr></thead><tbody><tr><td>A02_199T01_9000000708</td><td>57666.25</td><td>2886.3</td><td>1.85 ± 0.02</td><td>0.45 ± 0.06</td><td>267.3/230</td><td>3.71 ± 0.06</td><td>1.46 ± 0.12</td></tr><tr><td>A02_199T01_9000000708</td><td>57666.39</td><td>4041.7</td><td>1.77 ± 0.02</td><td>0.56 ± 0.05</td><td>348.9/267</td><td>3.81 ± 0.05</td><td>1.61 ± 0.1</td></tr><tr><td>A02_199T01_9000000708</td><td>57666.52</td><td>2938.6</td><td>1.81 ± 0.02</td><td>0.5 ± 0.06</td><td>228.5/238</td><td>3.94 ± 0.06</td><td>1.56 ± 0.12</td></tr><tr><td>A02_199T01_9000000708</td><td>57666.6</td><td>3430.7</td><td>1.84 ± 0.02</td><td>0.45 ± 0.05</td><td>280.9/248</td><td>3.76 ± 0.05</td><td>1.51 ± 0.12</td></tr><tr><td>G06_086T01_9000000774</td><td>57695.9</td><td>380.4</td><td>1.99 ± 0.07</td><td>0.29 ± 0.18</td><td>46.9/53</td><td>3.98 ± 0.19</td><td>1.03 ± 0.29</td></tr><tr><td>G06_086T01_9000000774</td><td>57695.93</td><td>1286.2</td><td>1.89 ± 0.04</td><td>0.47 ± 0.08</td><td>146.2/159</td><td>3.97 ± 0.08</td><td>1.32 ± 0.13</td></tr><tr><td>G06_086T01_9000000774</td><td>57696.0</td><td>924.8</td><td>1.87 ± 0.04</td><td>0.56 ± 0.1</td><td>112.1/114</td><td>3.77 ± 0.1</td><td>1.3 ± 0.13</td></tr><tr><td>G06_086T01_9000000774</td><td>57696.07</td><td>577.7</td><td>1.9 ± 0.06</td><td>0.38 ± 0.14</td><td>75.1/74</td><td>3.66 ± 0.13</td><td>1.34 ± 0.28</td></tr><tr><td>G06_086T01_9000000774</td><td>57696.13</td><td>233.0</td><td>1.88 ± 0.09</td><td>0.24 ± 0.21</td><td>33.4/33</td><td>4.02 ± 0.22</td><td>1.78 ± 1.17</td></tr><tr><td>G06_086T01_9000000774</td><td>57696.19</td><td>706.1</td><td>1.93 ± 0.05</td><td>0.25 ± 0.12</td><td>94.5/86</td><td>3.89 ± 0.12</td><td>1.39 ± 0.39</td></tr><tr><td>G06_086T01_9000000774</td><td>57696.26</td><td>370.9</td><td>1.78 ± 0.08</td><td>0.54 ± 0.18</td><td>47.1/53</td><td>3.95 ± 0.17</td><td>1.61 ± 0.38</td></tr></tbody></table>

(This table is available in its entirety in machine-readable form.)

<table><caption>**Table 4**

Details of AstroSat Observations</caption><thead><tr><th>S.N.</th><th>Date of Observation (UTC)</th><th>Time Start/Range (MJD)</th><th>SXT Exp. (ks)</th><th>LAXPC Exp. (ks)</th><th>Time Seg.</th></tr></thead><tbody><tr><td>1</td><td>2015-11-20T05:28:57</td><td>57346.27</td><td>2.9</td><td>...</td><td>PV</td></tr><tr><td>2</td><td>2016-10-05T03:49:57</td><td>57666.2–57666.6</td><td>13.4</td><td>25.4</td><td>T0</td></tr><tr><td>3</td><td>2016-11-03T20:47:00</td><td>57695.9</td><td>82.6</td><td>147.8</td><td>F1: T1, T2, T3<br>T4, T5, T6</td></tr><tr><td></td><td></td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td></td><td></td><td>57695.9–57696.4</td><td>8.6</td><td>30.0</td><td>T1</td></tr><tr><td></td><td></td><td>57696.5–57697.3</td><td>13.5</td><td>27.4</td><td>T2</td></tr><tr><td></td><td></td><td>57697.3–57698.3</td><td>21.6</td><td>18.7</td><td>T3</td></tr><tr><td></td><td></td><td>57698.3–57699.4</td><td>24.3</td><td>17.5</td><td>T4</td></tr><tr><td></td><td></td><td>57696.5–57697.8</td><td>26.6</td><td>46.3.3</td><td>T5</td></tr><tr><td></td><td></td><td>57697.8–57699.0</td><td>25.6</td><td>17.4</td><td>T6</td></tr><tr><td></td><td></td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>4</td><td>2016-11-16T10:30:11</td><td>57708.4–57709.1</td><td>16.0</td><td>33.3</td><td>T7</td></tr><tr><td>5</td><td>2017-10-25T02:28:50</td><td>58051.1</td><td>35.1</td><td>85.7</td><td>F2: T8, T9, T10</td></tr><tr><td></td><td></td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td></td><td></td><td>58051.1–58051.7</td><td>14.15</td><td>38.7</td><td>T8</td></tr><tr><td></td><td></td><td>58051.7–58052.3</td><td>11.0</td><td>25.7</td><td>T9</td></tr><tr><td></td><td></td><td>58052.3–58052.6</td><td>9.9</td><td>16.7</td><td>T10</td></tr><tr><td></td><td></td><td>...</td><td>...</td><td>...</td><td>...</td></tr></tbody></table>

**Note.** The small time segments T1–T6 and T8–T10 are used to generate time-resolved spectra to understand the various phases of flaring activities in 2016 and 2017. The time segment PV represents the first target of opportunity (ToO) observations made in 2015.

mission-specific tools such as *uvotimsum*, *uvotsource*, and *uvot2pha* distributed with the HEAsoft package. The sky images in a particular filter corresponding to individual observations are combined using *uvotimsum* to get a single frame per observation, whenever more than one image was taken. The combined images are then analyzed utilizing the tool *uvot-source* using a circular region of 5" radius centered at the sky location of 1ES 1959+650 as the source region. Another circular region of 35.76" located in a source-free region around 3.5' away from 1ES 1959+650 is used to extract background counts.

A correction due to reddening, $E(B-V) = 0.178$, due to the presence of the neutral hydrogen along the line of sight within our own Milky Way, is applied to the fluxes before using these values in SEDs. The reddening is estimated by the Python module *extinctions* using the two-dimensional dust map of the entire sky by Schlegel et al. (1998), which was recently updated by Schlafly & Finkbeiner (2011; hereafter SFD). The estimation of the same parameter using the two-dimensional

dust map at NASA/IPAC archive<fnref n="16" /> yields a value of 0.172. We also estimate this parameter using the recent three-dimensional dust map by Green et al. (2015), which turns out to be 0.180. This implies that we can safely use 0.178 measured using SFD. The empirical formalism by Cardelli et al. (1989) with $A_V = R_V$ * $E(B-V)$ and $R_V = 3.1$ is used to estimate the correction factor $A_\lambda$ for individual UVOT filters. Multiple UVOT observations taken during periods over which individual SEDs were collected were averaged to give one data point per filter.

### 2.3.2. Mt. Abu Infrared Observatory, India

In addition to the optical/UV observations from the Neil Gehrels Swift Observatory UVOT, optical photometry observations from the MIRO are used in this investigation.

<footnote n="16">¹⁶ https://irsa.ipac.caltech.edu/frontpage/</footnote>
<figure><img src="image_1.png" /><figcaption>**Figure 1.** X-ray light curves in the energy range of 0.3–7.0 keV obtained using XRT and SXT data from 2015 December to 2021 February. The red stars represent XRT fluxes, whereas the black open circles show SXT fluxes. The shaded region around the best-fit dotted curve represents the 3σ confidence interval.</figcaption></figure>

A number of optical photometric observations were made during several epochs between 2015 December and 2017 December using the MIRO, including several observations contemporaneous to the AstroSat monitoring in 2016. Table 1 provides further information about these observations. The data were obtained using the EMCCD-based optical camera installed at the f/13 Cassegrain focus of the 1.2 m telescope. The data reduction and the photometry procedures adopted are discussed in Kaur et al. (2017). Differential photometry, using several comparison stars in the same frame as the source, was used to minimize atmospheric seeing effects. The calibrated magnitudes thus obtained were converted to the fluxes and corrected for Galactic extinction. The nightly averaged fluxes in mJy are shown in the bottom panel of Figure 3.

## 2.4. Other Publicly Available Resources

For coverage at lower frequencies, we use the publicly available radio data from the Owens Valley Radio Observatory (OVRO)<fnref n="17" /> at 15 GHz (Richards et al. 2011). For completeness of the SED, we extract the quasi-simultaneous TeV data from MAGIC Collaboration et al. (2020), obtained by the MAGIC telescope, which have been corrected for γ-γ absorption by extragalactic background light (EBL). The TeV γ-ray quick-look light curve from FACT (First G-APD Cerenkov Telescope; Anderhub et al. 2013; Biland et al. 2014; Dorner et al. 2015) from 2016 September to 2017 November are also used to investigate the high-energy counterparts of the observed X-ray activities.

# 3. Results

This section presents the results of our comprehensive multiwavelength study of flux and spectral variability of 1ES 1959+650 using AstroSat and other facilities.

## 3.1. Light Curves and Flux Variations

The X-ray light curve in the energy band of 0.3–7.0 keV taken over 6 yr (between MJD 57000 and 59260) is shown in

Figure 1. The open stars symbolize the integrated fluxes from XRT, whereas the open circles represent the fluxes from SXT. The long SXT exposures are split into several small time intervals of 3500 s, within which spectra are extracted and the best-fit fluxes are used to construct the light curve (see Section 2.2.4 for details). However, in order to generate the light curve from XRT, the fluxes are extracted from the best-fit model spectra from the individual observations between 2015 January and 2021 February.

The overall long-term average flux variation trend is mathematically characterized by a broad Gaussian peaking at ~MJD 58233 and FWHM of ≈785 days. The long-term trend is superimposed by several flares.

The 6 yr long XRT light curve (Figure 1) is divided into seven segments of 300 days each, except the last one, which corresponds to 440 days. These segments, R1, R2, R3, R4, R5, R6, and R7, encompass a number of X-ray flares at different epochs sampling the different parts of the above-mentioned Gaussian function (See Figure 1). The AstroSat observations from 2016 and 2017 fall into R3 and R4, respectively. The following paragraphs summarize the quantitative analysis of the multiwavelength variability during various X-ray flares around the AstroSat monitoring.

The SXT light curves reveal that 1ES 1959+650 exhibits significant flux variations with different timescales at all the epochs as shown in Figure 2. Therefore, in order to understand the nature of the flux variability during and around the AstroSat observations, four time segments are created. The basic criterion behind the division of these variability profiles is to distinguish and characterize the X-ray outbursts (doubling/halving timescales and peak flux), which are probably related to the same physical processes that triggered the X-ray activities recorded by AstroSat. These segments are called variability profiles and are denoted by V1, V2, V3, and V4 (See Figure 2). Note that V1, V2, and V3 are subsets of R3, while V4 is a subset of R4. In order to investigate the time-dependent spectral behavior of the flaring activities observed with AstroSat, small portions of V1, V2, V3, and V4 are further subdivided into time segments denoted by T0, T1,...,T10 (See Figure 2 and Table 4 for details). These segments zoom in on the flux variations during the AstroSat monitoring. Further AstroSat monitoring data with a total exposure of ~2.9 ks from 2015 November, denoted by “PV” are also included in this list to compare the X-ray activities at earlier epochs (see Section 3.2 for details).

The characteristic flux doubling/halving timescales (Δt_D/Δt_H) of the flares in the different variability profiles (V1, V2, V3, V4), using the combined X-ray light curve from both SXT and XRT, are derived with Δt = t_d × ln 2/|ln(F2/F1)| (Saito et al. 2013). Here, F1 and F2 are the fluxes observed at a time interval of t_d. The methodology is also applied to the detailed X-ray light curve from SXT to obtain the parameters over shorter timescales (see Table 5 for details). The variations are further characterized by the fractional variability amplitude (F_var) defined in Vaughan et al. (2003) and given by

$$F_{\text{var}} = \sqrt{\frac{S^2 - \sigma_{\text{err}}^2}{x_m^2}} \qquad (1)$$

Here σ err is the mean square error of each observations and S² is the sample variance, where σ XS = S² - σ err² is the excess variance. x_m is the unweighted sample mean for N points. The

<footnote n="17">¹⁷ https://sites.astro.caltech.edu/ovroblazars/index.php?page=home</footnote>
<figure><img src="image_2.png" /><figcaption>Figure 2. Zoomed-in light curves focusing on the time around the flares observed by AstroSat in 2016 and 2017. Left: various time segments (V1, ..., V3) and subgroups (T0, ..., T7) during 2016, zooming into the flaring activities in the top and bottom panels. Right: time segment V4 and various subgroups (T8, T9, T10) during 2017, zooming into the bright flaring activities in the bottom panel. The LAXPC light curves are also shown by gray plus signs in the zoomed windows after scaling up to fit in the panels.</figcaption></figure>

<figure><img src="image_3.png" /></figure>

error in $F_{var}$ is given by

$$ \sigma_{F_{var}} = \sqrt{\left(\sqrt{\frac{1}{2N}} \frac{\sigma_{err}^2}{x_m^2 F_{var}}\right)^2 + \left(\sqrt{\frac{\sigma_{err}^2}{N}} \frac{1}{x_m}\right)^2} \quad (2) $$

The variability timescales ($\Delta t_D/\Delta t_H$), $F_{var}$, and peak flux ($F_{X, p}$) are derived for the flares in each time segment (V1, V2, V3, V4) of the long light curves shown in Figure 2, including both XRT and combined XRT/SXT observations. The values are reported in Table 5.

The multiwavelength light curves shown in Figure 3 illustrate the prominent X-ray activities and its counterparts in other energy bands. Some of the X-ray outbursts seem to have associated high-energy counterparts in GeV and TeV bands. There have been several communications from the FACT collaboration through “Astronomers Telegrams (ATels)” (Biland et al. 2016a, 2016b; Biland & FACT Collaboration 2016a, 2016b; Buson et al. 2016) reporting fluxes beyond 1 Crab unit (CU), and 36 private communications to the partners to trigger multiwavelength observations during moderately high flux states ($F_{\gamma} \ge 0.5$ CU) over the period of 2016–2017. The FACT quick look analysis (QLA) light curve<fnref n="18" /> of 1ES 1959+650 shows several flaring

episodes. In the following subsections, the detailed variability profiles, the nature of the X-ray flux variations, and its multiwavelength associations are discussed.

### 3.1.1. Variability Profile 1 (V1)

The variability profile V1 refers to the light curves around T0 starting from MJD 57660.0 to 57670.0. It also includes four pointings by Swift. The X-ray light curves (see Figure 2, left panels) show that the observations performed during T0 have been part of a fast-varying flux state ($F_{var} \sim 0.25 \pm 0.005$) with initial rise ($\Delta t_D \sim 2.7$ days) and fall ($\Delta t_D \sim 8.7$ days). The X-ray light curves in 0.3–7.0 keV (SXT) and 3.0–30.0 keV (LAXPC) bands over T0 show significant flare-like variations (see Table 5 for details). The flux variation in the 3–30.0 keV band is more prominent than in the 0.3–7.0 keV band. As shown in the left panels of Figure 3, V1 comprises hints of flux variations in the TeV band, whereas no variations are observed in the optical/UV and GeV bands.

### 3.1.2. Variability Profile 2 (V2)

V2 corresponds to the light curves during MJD 57682.6 and 57700.0 (see Figure 2, left panels), which starts 12.5 days after the end of V1. The combined light curves clearly show that V2

<footnote n="18">¹⁸ https://fact-project.org/monitoring</footnote>
<table><caption>**Table 5**
Estimation of $\Delta t$ and $F_{\text{var}}$</caption><thead><tr><th>T<sub>tag, prof.</sub></th><th>Duration<br>(MJD)</th><th>Δt<sub>D</sub>/Δt<sub>H</sub><br>(days)</th><th>F<sub>var, com</sub>/F<sub>var, sxt</th><th>F<sub>X, p</sub><br>(×10<sup>-10</sup> erg cm<sup>-2</sup> s<sup>-1</sup>)</th><th>T<sub>tag, sxt</sub></th></tr></thead><tbody><tr><td rowspan="6"><b>V1</b></td><td rowspan="6"><b>57660.0–57670.0</b></td><td>↑ XRT: 2.69 ± 0.02</td><td>COM: 0.25 ± 0.005</td><td rowspan="3">4.12</td><td rowspan="6">T0</td></tr><tr><td>↓ COM: 8.68 ± 0.05</td><td rowspan="2">SXT: 0.02 ± 0.004</td></tr><tr><td>↑ SXT: 3.09 ± 0.17</td></tr><tr><td rowspan="2">↓ SXT: 1.69 ± 0.01</td><td rowspan="2">LAXPC: 0.07 ± 0.002</td><td rowspan="2">T0</td></tr><tr><td>↑ LAXPC: 1.01 ± 0.008</td></tr><tr><td>↓ LAXPC: 0.63 ± 0.01</td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="2"><b>V2</b></td><td rowspan="2"><b>57682.6–57700.0</b></td><td>...</td><td>COM: 0.19 ± 0.002</td><td rowspan="2">7.05</td><td rowspan="2"></td></tr><tr><td>↑ XRT: 8.59 ± 0.06</td><td>SF1: 0.31 ± 0.006</td></tr><tr><td rowspan="4">SF1</td><td rowspan="4">57682.70–57695.25</td><td>↓ XRT: 2.31 ± 0.01</td><td></td><td rowspan="4">5.74</td><td rowspan="2">T1+T2+(T3/2)<br>(T3/2)+T4</td></tr><tr><td>↑ COM: 2.92 ± 0.003</td><td>SF2: 0.15 ± 0.003</td></tr><tr><td>↓ COM: 1.95 ± 0.01</td><td rowspan="2">LAXPC: 0.26 ± 0.001</td><td rowspan="2">T1+T2+(T3/2)<br>(T3/2)+T4</td></tr><tr><td>↑ LAXPC: 1.49 ± 0.01</td></tr><tr><td rowspan="2">SF2</td><td rowspan="2">57696.30–57700.00</td><td>↓ LAXPC: 0.65 ± 0.009</td><td></td><td></td><td></td></tr><tr><td></td><td></td><td></td><td></td></tr><tr><td rowspan="10"><b>V3</b></td><td rowspan="10"><b>57700.36–57714.3</b></td><td>...</td><td>COM: 0.18 ± 0.004</td><td rowspan="2">6.78</td><td rowspan="10"></td></tr><tr><td>↑ XRT: 5.93 ± 0.07</td><td>SF1: 0.30 ± 0.007</td></tr><tr><td>↓ COM: 5.19 ± 0.08</td><td></td><td></td><td></td></tr><tr><td rowspan="5">SF1</td><td rowspan="5">57709.40–57713.80</td><td>↑ COM: 8.89 ± 0.03</td><td>SF2: 0.26 ± 0.009</td><td rowspan="5">6.94</td></tr><tr><td>↓ XRT: 1.10 ± 0.001</td><td></td><td></td><td></td></tr><tr><td>↑ SXT: 1.23 ± 0.001</td><td>SXT: 0.09 ± 0.003</td><td>T7</td><td></td></tr><tr><td>...</td><td>LAXPC: 0.15 ± 0.002</td><td>T7</td><td></td></tr><tr><td>↑ LAXPC: 0.36 ± 0.003</td><td>LAXPC, SF1: 0.17 ± 0.003</td><td></td><td></td></tr><tr><td rowspan="4">SF2</td><td rowspan="4">57709.40–57713.80</td><td>↓ LAXPC: 0.53 ± 0.005</td><td rowspan="2">LAXPC, SF2: 0.12 ± 0.003</td><td rowspan="4">12.45</td><td rowspan="4"></td></tr><tr><td>↑ LAXPC: 0.52 ± 0.004</td></tr><tr><td>↓ LAXPC: 0.18 ± 0.003</td><td></td><td></td><td></td></tr><tr><td>↑ XRT: 10.81 ± 0.04</td><td>COM: 0.18 ± 0.002</td><td></td><td></td></tr><tr><td rowspan="10"><b>V4</b></td><td rowspan="10"><b>58032.0–58058.0</b></td><td>↓ XRT: 5.06 ± 0.03</td><td>SF1: 0.56 ± 0.005</td><td rowspan="10">14.8</td><td rowspan="10"></td></tr><tr><td>...</td><td>SXT: 0.09 ± 0.008</td><td>T8+T9+T10</td><td></td></tr><tr><td>↑ XRT: 5.06 ± 0.03</td><td>T8: 0.07 ± 0.003</td><td>T8</td><td></td></tr><tr><td>↓ SXT: 2.06 ± 0.03</td><td>T9: 0.10 ± 0.009</td><td>T9</td><td></td></tr><tr><td>↑ SXT: 1.14 ± 0.01</td><td>T10: 0.04 ± 0.004</td><td>T10</td><td></td></tr><tr><td>↓ SXT: 1.60 ± 0.009</td><td>LAXPC: 0.15 ± 0.003</td><td>T8+T9+T10</td><td></td></tr><tr><td>...</td><td>T8: 0.17 ± 0.001</td><td>T8</td><td></td></tr><tr><td>↓ LAXPC: 0.63 ± 0.004</td><td>T9: 0.13 ± 0.001</td><td>T9</td><td></td></tr><tr><td>↑ LAXPC: 0.54 ± 0.003</td><td>T10: 0.07 ± 0.001</td><td>T10</td><td></td></tr><tr><td>↓ LAXPC: 0.46 ± 0.003</td><td></td><td></td><td></td></tr></tbody></table>

**Note.** T<sub>tag, prof.</sub> and T<sub>tag, sxt</sub> refer to the tags adopted for defining the variability profiles and time segments within the SXT light curves, respectively (see the left and right panels of Figure 2 for details). SF represents a small flare observed in each segment, where ↑ and ↓ represent its doubling (Δt<sub>D</sub>) and halving (Δt<sub>H</sub>) timescales, respectively, estimated for various flares. The subscripts XRT, SXT, LAXPC, and COM. refer to the data used from XRT, SXT, LAXPC, and combined from both XRT/SXT, respectively. F<sub>X, p</sub> represents the peak flux in a particular time segment.

encompasses two consecutive X-ray flares, with the first one barely covered by XRT (SF1) and the second one entirely observed by AstroSat (SF2). The peak-to-peak time difference of the SF1 and SF2 is ~6.9 days.

SF1 is highly asymmetric in X-rays and shows noticeable TeV and GeV activities. SF2, however, while also highly asymmetric in X-rays, hardly shows any variations at other energies. The fastest variations during V2 are observed with LAXPC during SF2 (Δt<sub>D</sub> = 1.5 days and Δt<sub>H</sub> = 0.65 days). Hence, the X-ray variations during both SF1 and SF2 can be characterized by “slow rise and fast decay” profiles (see Table 5).

The AstroSat light curves are further subdivided into time segments T1, T2, T3, T4, T5, and T6 (see Figure 2, left panels). In addition to the average flux variations over the AstroSat observations, the variability timescales and fractional variability indicate significant variations even during the above-mentioned T segments (see Table 5). Combining the optical/UV fluxes from MIRO and UVOT, we deduce that no

significant optical and UV counterparts of the X-ray flares during V2 are seen (see left panels of Figure 3).

### 3.1.3. Variability Profile 3 (V3)

The variability profile V3 extends from MJD 57704.3 to 57714.3, that is, it starts soon after the end of V2. The combined XRT and SXT X-ray light curves during V3 (see Figure 2, left panels) show that AstroSat pointing T7 is most probably the falling part of an X-ray flare that peaks at a flux (6.78 × 10<sup>-10</sup> erg cm<sup>-2</sup> s<sup>-1</sup>) comparable to that of V2-SF1. The V3 profile contains two (SF1 and SF2) X-ray flares of nearly similar peak fluxes, with the second one decreasing very fast (Δt<sub>H</sub> = 1.1 days). Both SF1 and SF2 show the “slow rise and fast decrease” already seen in V2. A similar trend is also seen in the LAXPC light curve (3.0–30 keV band); however, the fractional variability is higher than in the SXT band (see Table 5). There are no obvious counterparts in other
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Left: multiwavelength light curves of 1ES 1959+650. From top to bottom: FACT (quick look analysis, binned nightly), Fermi-LAT (binned in 1 day and 2 days), SXT and XRT (binned by orbit), LAXPC and UVOT (U, V, B, UVW1, UVM2, and UVW2 bands), and MIRO (B, V, and R bands, binned nightly) during 2016 October–November (corresponds to segments V2 and V3). Right: same as the left panels, but for the observations during 2017 November (corresponds to segment V4). Note that 95% upper limits instead of flux points are used for time bins where the test statistics (TS) in the LAT light curves are $TS \le 25$.</figcaption></figure>

energy bands, with the noteworthy exception of a TeV enhancement during SF2.

### 3.1.4. Variability Profile 4 (V4)

The profile V4 extends from MJD 58032.0 to 58058.0 (see Figure 2, right panels). This time segment signifies the X-ray flux variability around the major activity observed in November 2017.

V4 includes four pointings of Swift, with one coinciding with the AstroSat observations. The combined X-ray light curve during V4 shows that the X-ray flare observed with AstroSat is part of a prominent MWL activity lasting over ~24 days. The Swift light curve alone shows that the source has doubled its flux in merely 15 days from MJD 58033.0 to 58048.0, reaching a record flux of $12.45 \times 10^{-10} \text{ erg cm}^{-2} \text{ s}^{-1}$ in the 0.3–7.0 keV band (see Figure 2, right panels). The AstroSat pointing was made 3 days after the XRT peak, that is, MJD 58051.0. The AstroSat observations in V4 reveal the highest flux state ever observed ($14.8 \times 10^{-10} \text{ erg cm}^{-2} \text{ s}^{-1}$).

The observations with AstroSat lasting ~2 days reveal two major subflares. The doubling and halving times indicate a slow decline and fast rise in them. Figure 5 and Table 2 reveal a

significant shift in the synchrotron peak position ($E_{s,p}$) over the flare compared to the observations in 2016.

The LAXPC light curve shows a similar behavior to the SXT light curve but exhibits a more pronounced amplitude variability, which is underlined by a higher excess fractional variance ($EV = 0.15$) compared to the SXT light curve ($EV = 0.08$). As shown in the right panels of Figure 3, V4 comprises GeV and TeV flares nearly 3 days prior to the observed X-ray peak. The GeV/TeV fluxes double within ~1 day. Interestingly, the GeV flare exhibits a narrower profile than the TeV flare and leads the TeV flare by ~1 day. The lack of X-ray data prohibits us from judging whether the GeV/TeV peak corresponds to the (in this case unknown) X-ray maximum or whether the X-ray response is delayed. The subflare observed with AstroSat, however, corresponds to a low state in the GeV/TeV bands. No significant UV/optical counterparts are seen.

The correlated X-ray activity, along with the shift of the peak position mentioned above, implies that it is not a simple variation in electron density, but a more complicated spectral change throughout the flare that is responsible for the X-ray flare. This requires time-dependent changes in the magnetic field, the maximum electron Lorentz factor, the Doppler factor
<table><caption>**Table 6**

Best-fit Parameters for Time-resolved X-ray Spectra</caption><thead><tr><th rowspan="2">T<sub>Tag</sub></th><th colspan="4">Log Parabola</th><th colspan="4">Exponential Cutoff Power Law</th></tr><tr><th>C</th><th>α</th><th>β</th><th>F<sub>X, 0.3-8.0keV</sub> (×10<sup>-10</sup> erg cm<sup>-2</sup> s<sup>-1</sup>)</th><th>E<sub>s,p</sub> (keV)</th><th>χ<sup>2</sup>/dof</th><th>Γ</th><th>E<sub>cut</sub></th><th>F<sub>X, 0.3-8.0keV</sub> (×10<sup>-10</sup> erg cm<sup>-2</sup> s<sup>-1</sup>)</th><th>χ<sup>2</sup>/dof</th></tr></thead><tbody><tr><td>PV<sup>a</sup></td><td>...</td><td>2.48 ± 0.04</td><td>...</td><td>2.09 ± 0.04</td><td>...</td><td>1.21/197</td><td>2.37 ± 0.04</td><td>13.68</td><td>2.02 ± 0.02</td><td>1.31/195</td></tr><tr><td>T0</td><td>0.85 ± 0.03</td><td>1.98 ± 0.02</td><td>0.36 ± 0.02</td><td>3.17 ± 0.02</td><td>1.07 ± 0.07</td><td>1.6/142</td><td>1.97 ± 0.02</td><td>13.68<sup>+1.24</sup><sub>-1.06</sub></td><td>3.26 ± 0.02</td><td>1.80/142</td></tr><tr><td>T1</td><td>0.85 ± 0.03</td><td>2.07 ± 0.02</td><td>0.33 ± 0.02</td><td>3.20 ± 0.02</td><td>0.78 ± 0.06</td><td>1.08/311</td><td>2.11 ± 0.03</td><td>14.75<sup>+1.43</sup><sub>-1.21</sub></td><td>3.16 ± 0.02</td><td>0.91/288</td></tr><tr><td>T2</td><td>0.84 ± 0.03</td><td>2.03 ± 0.02</td><td>0.36 ± 0.02</td><td>3.57 ± 0.02</td><td>0.91 ± 0.05</td><td>1.29/360</td><td>2.06 ± 0.02</td><td>13.28<sup>+0.99</sup><sub>-0.87</sub></td><td>3.69 ± 0.02</td><td>1.21/330</td></tr><tr><td>T3</td><td>0.83 ± 0.02</td><td>1.99 ± 0.01</td><td>0.38 ± 0.02</td><td>4.42 ± 0.02</td><td>1.03 ± 0.03</td><td>1.78/397</td><td>2.04 ± 0.02</td><td>13.07<sup>+0.91</sup><sub>-0.80</sub></td><td>4.57 ± 0.02</td><td>1.09/367</td></tr><tr><td>T4</td><td>0.89 ± 0.03</td><td>2.09 ± 0.01</td><td>0.40 ± 0.02</td><td>3.14 ± 0.01</td><td>0.77 ± 0.03</td><td>1.55/416</td><td>2.13 ± 0.02</td><td>11.43<sup>+0.82</sup><sub>-0.76</sub></td><td>3.26 ± 0.01</td><td>1.14/386</td></tr><tr><td>T5</td><td>0.81 ± 0.02</td><td>2.01 ± 0.01</td><td>0.37 ± 0.02</td><td>4.00 ± 0.01</td><td>0.97 ± 0.03</td><td>1.60/456</td><td>2.05 ± 0.02</td><td>13.10<sup>+0.79</sup><sub>-0.73</sub></td><td>4.14 ± 0.01</td><td>1.27/426</td></tr><tr><td>T6</td><td>0.72 ± 0.02</td><td>2.05 ± 0.01</td><td>0.40 ± 0.02</td><td>3.64 ± 0.01</td><td>0.87 ± 0.03</td><td>1.67/447</td><td>2.08 ± 0.02</td><td>11.29<sup>+0.71</sup><sub>-0.71</sub></td><td>3.77 ± 0.01</td><td>1.23/412</td></tr><tr><td>T7</td><td>1.03 ± 0.04</td><td>1.98 ± 0.02</td><td>0.52 ± 0.03</td><td>3.72 ± 0.02</td><td>1.05 ± 0.05</td><td>1.84/132</td><td>1.96 ± 0.02</td><td>9.41<sup>+0.62</sup><sub>-0.56</sub></td><td>3.71 ± 0.02</td><td>1.70/132</td></tr><tr><td>T8</td><td>0.81 ± 0.02</td><td>1.84 ± 0.02</td><td>0.25 ± 0.02</td><td>10.7 ± 0.06</td><td>2.09 ± 0.19</td><td>1.33/133</td><td>1.86 ± 0.02</td><td>25.28<sup>+1.92</sup><sub>-1.68</sub></td><td>10.09 ± 0.06</td><td>1.41/135</td></tr><tr><td>T9</td><td>0.92 ± 0.03</td><td>1.85 ± 0.02</td><td>0.29 ± 0.003</td><td>10.56 ± 0.06</td><td>1.89 ± 0.12</td><td>1.45/139</td><td>1.86 ± 0.02</td><td>19.56<sup>+1.51</sup><sub>-1.32</sub></td><td>10.79 ± 0.06</td><td>1.84/127</td></tr><tr><td>T10</td><td>0.78 ± 0.02</td><td>1.65 ± 0.02</td><td>0.41 ± 0.02</td><td>12.47 ± 0.07</td><td>2.67 ± 0.18</td><td>1.89/129</td><td>1.75 ± 0.02</td><td>14.48<sup>+0.79</sup><sub>-0.72</sub></td><td>12.51 ± 0.07</td><td>1.79/133</td></tr></tbody></table>

**Notes.** The C parameter refers to the relative cross-normalization for LAXPC keeping the SXT normalization fixed to 1. The exponential cutoff power law seems to show low reduced $\chi^2$ values in comparison to that from the Log Parabola model because of higher systematic (4% instead of 3%). $E_{s,p}$ values are evaluated from the best-fit Log Parabola parameters as described in the Appendix.

<sup>a</sup> The data fit well with the power-law instead of the Log Parabola model. Also, due to poor statistics, the cutoff energy is fixed to 13.68 keV as observed during next observations.

of the emission region, or a combination thereof. A more detailed analysis is provided below.

## 3.2. X-Ray Spectral Modeling

The combined X-ray spectra from SXT and LAXPC covering the 0.3–30 keV band from various epochs are fitted to derive the time-dependent spectral behavior of the source. For these investigations, the total AstroSat observations are split into 12 segments (see Table 4 and Figure 2). These time segments are designated as PV, T0, ..., T10. The PV segment refers to the ~3 ks AstroSat pointing of 1ES 1959+650 performed on 2015 November 20. These segments are defined to highlight the changes in the X-ray spectral parameters sampling different parts of the flares and also different average flux states (see Table 4 and Figure 2 for the definitions of these time segments).

The combined SXT+LAXPC spectra from the 12 segments are individually modeled with two spectral models: (1) TBabs * LogParabola, hereafter M1, and (2) TBabs * CutoffPl, hereafter M2. The absorption model component TBabs with the WILM abundance model (Wilms et al. 2000) is used to fit the Galactic neutral hydrogen absorption in the source direction. The LogParabola and CutoffPl models fit a log parabola and a cutoff power-law shape to the intrinsic spectra, respectively. The input parameter $n_H$ is kept fixed to $0.107 \times 10^{-22} \text{ cm}^{-2}$ in models M1 and M2 throughout. This value is estimated using 21 cm observations in the source direction (GAL. Coordinate $L = 98.003370, B = 17.669746$) using the online tool<sup>19</sup> by the LAB Survey (Kalberla et al. 2005, and relevant references therein) with the default 0.27 beam. The PV segment does not have a usable LAXPC spectrum, and hence the spectral results correspond to the SXT observations (0.3–7.0 keV) only. The $\chi^2$ statistics is used for the spectral modeling throughout.

The best-fit parameters from the spectral fitting, along with their $2\sigma$ uncertainties, are listed in Table 6. The best-fit results show that the observed X-ray spectra can be represented equally well by M1 and M2. We prefer M1 over M2, as (1) it has been used in previous studies, thus allowing for an easy comparison, and (2) the fitting with M2 requires higher systematic uncertainties to be added to the statistical errors to converge the fit. On the other hand, the PV spectrum is fitted with the absorbed power-law model (TBabs*powerlaw with fixed $n_H$).

In order to get acceptable $\chi^2$ values, 4%–5% systematic uncertainties are needed in many cases. Normally, 3% systematic uncertainties are recommended for the SXT+LAXPC joint spectral fitting. The top left, top right, and bottom left panels of Figure 4 show the observed spectra for T0 to T10 with the respective best-fit models. The bottom right panel of Figure 4 presents the model-generated spectra in the energy band 0.2–12.0 keV. The colored regions around the model spectra are $1\sigma$ uncertainty intervals for the best-fit parameters.

Figure 4 clearly establishes prominent flux and spectral variations. The vertical dashed lines show the positions of the synchrotron peaks ($E_{s,p}$), which are calculated using Equation (A5) and reported in Table 6. The projected synchrotron peak position $E_{s,p}$, pv is below 0.2 keV. It also represents the faintest state as observed with AstroSat. The spectra observed in 2016 (T0, ..., T7) clearly show significant variations in $E_{s,p}$ throughout. Figure 5 illustrates the shift in the synchrotron peaks for various time bins; see also Table 4. Within the $1\sigma$ confidence intervals (Figure 5), the shifts in $E_{s,p}$ are correlated to the flux states for the outbursts in 2016 and 2017.

This investigation establishes that the spectra of 1ES 1959+650 in the 0.3–30 keV band become harder with increasing X-ray flux, and subsequently, the peak of the synchrotron emission shifts toward higher energies. The hardening of spectra and shift in $E_{s,p}$ is linked to the particle energization. These motivate us to investigate the relationship between fluxes

<footnote n="19"><sup>19</sup> https://www.astro.uni-bonn.de/hisurvey/profile/</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 4.** X-ray spectra fitted with a log parabola model. Top left: $E^2F_E$ vs. E for T0 (black) and T7 (red) extracted from AstroSat observations performed during 2016 October 05–06 and during November 16–17. Top right: same as the top left panel, but for periods T1, T2, T3, T4, T5, and T6. The various spectra represent different flux and spectral states segmenting different parts of the outburst observed from MJD 57695.9 to 57699.0 (2016 November 03–07). Bottom left: same as the top left panel, but for T8, T9, and T10 from the outburst from MJD 58051.1 to 58052.6, i.e., 2017 October 25–26. Bottom right: $E^2F_E$ vs. E plots, with the respective butterfly diagram representing the best-fit parameters and their uncertainties for all 16 X-ray spectra shown in the four panels.</figcaption></figure>

and $E_{s,p}$ further by utilizing (a) modeling SXT spectra sampling smaller portions of the flares and (b) fitting around 400 XRT spectra observed between 2015 January and 2021 February. The spectral modeling of the SXT spectra with such a sampling provides us with a unique independent data set representing the two major outbursts in 2016 and 2017, while the integration of the XRT data, sampling the X-ray variations over 6 yr, adds to a general understanding of the above relationship. Additional details are discussed in Section 3.3.

### 3.3. Correlation Analysis

This section describes the detailed correlation analysis between various spectral parameters and other dependent quantities for the time segments R1,...,R7 and two sets of AstroSat observations (Figure 1), separately. The $\alpha$ versus $F_{X,0.3-7\text{keV}}$ distribution is fitted with a straight line and is depicted in the left panel of Figure 6. The various pairs of colors and symbols represent the different R

segments. The lines are modeled with slope (m) and intercept (b) as free parameters. In order to extract the best-fit parameters, the least-squares regression (LS), maximum likelihood minimization (ML), and Monte Carlo Markov Chain (MCMC) techniques are used. For the MCMC calculations, the EMCEE tool in Python (Foreman-Mackey et al. 2013), which is an MIT licensed pure-Python implementation of the Goodman & Weare (2010) Affine Invariant MCMC Ensemble sampler, is used. The best-fit parameters derived from the LS and ML techniques are used as initial parameters for the MCMC technique. The thin colored lines in Figure 6 represent sample lines obtained with the MCMC method, whereas thick dashed lines correspond to the best fits derived with the MCMC technique.

A total of 81 time-resolved spectra extracted from SXT, as described in Section 2.2.4, are grouped into two parts. The first part covers the three epochs in 2016, and the second part is for the observations in 2017. These data sets are also modeled using the same methods as described above. The correlations of best-fit
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Positions of the synchrotron peaks estimated for the different time bins (T0 ...T10) defined in Table 4 and Figure 2.</figcaption></figure>

spectral parameters (α and β) versus X-ray flux (F<sub>X,0.3-7keV</sub>), E<sub>s,p</sub> versus F<sub>X,0.3-7keV</sub>, and α versus β are explored for all the R segments. The readers are referred to Table 7 for the details of the correlation indicators of various pairs from different R segments and respective best-fit parameters of the line.

The index α shows strong correlations with the integrated flux in the 0.3–7.0 keV band for all data segments. It is clear from this figure that most of the time the X-ray spectrum of 1ES 1959+650 gets harder with increasing flux. However, the slopes of the index versus F<sub>0.3-7.0keV</sub> relationships for the different time segments differ quite a lot. Especially, the slopes (m) for R1 and R7 are higher than for the other segments. Interestingly, R1 and R7 sample the increasing and decreasing tails of the Gaussian fit of the long-term X-ray light curve, respectively.

The (m, b) pairs, where m and b are slope and intercept for the fit on the α-F<sub>X,0.3-7keV</sub> plane, for the AstroSat observations in 2016 and R3 differ significantly, whereas the slope of AstroSat observations in 2017 matches the one in R4, though both showing different intercepts. Note that the AstroSat observations mainly sample spectra around particular outbursts, while the Swift observations are spread over a longer duration and hence represent a general behavior over long periods.

The right panels of Figure 6 show the flux dependence of the index (top) and the curvature parameter (bottom). The best-fit lines for the R segments are not plotted in the top panel for the sake of clarity. The two sets of AstroSat observations clearly represent two spectrally distinct states of variability of 1ES 1959+650. As clearly shown in Figure 6, the outburst in 2017 not only exhibits the highest X-ray flux but also follows an entirely different track in the spectral index—flux plane.

Figure 7 visualizes the confidence intervals and probability distribution for the best-fit parameters corresponding to the AstroSat observations in 2016 and R3, respectively. The plots are made using the Python package corner, which corresponds to the MCMC fitting method. These diagnostic plots are indicative of a reasonably good convergence and an acceptable quality of the parameter estimations. The correlations between flux enhancements and changes in the spectral parameters α and β may provide unique information about the energization of particles responsible for the flux enhancements during a

particular flare. Studies by other authors (Krawczynski et al. 2004; Kapanadze et al. 2018a, 2018b; MAGIC Collaboration et al. 2020) have shown that for a handful of blazars the curvature and index are inversely correlated with an increase in flux. More specifically, an increase in flux is accompanied by a hardening of the spectrum (decreasing α) and a decrease in curvature (decreasing β). The spectral constraints derived for the AstroSat-SXT and Swift-XRT observations confirm the strong positive correlations between the flux and the hardening of the spectrum. However, owing to the large uncertainties, the curvature parameter β (shown in the bottom right panel of Figure 6) does not show a significant correlation with the flux.

These spectral changes result in a shift of the peak position of the X-ray part of the SED and hence in the position of the synchrotron peak. Such a shift of the synchrotron peak toward higher energies with increasing integrated flux is seen in most HBLs and is equivalent to the “bluer-when-brighter” trend seen in the optical continuum spectra of many BL Lac objects (Böttcher et al. 2010; Dai et al. 2011; Kaur et al. 2017; Bhatta et al. 2018). The best-fit α and β values are used in Equation (A5) to derive the position of the peak of the synchrotron component.

The top panel of Figure 8 shows the relationship between the synchrotron peak energy (E<sub>s,p</sub>) and the spectral parameters α and β, respectively. The index parameter α is well correlated with the synchrotron peak position above an energy of ~0.5 keV, while such a statement cannot be made for the curvature parameter β.

The bottom panel of Figure 8 shows the relationship between E<sub>s,p</sub> and the integrated flux F<sub>X,0.3-7.0keV</sub>, confirming the trend of the synchrotron peak shifting to higher energies with increasing flux. While the relationship between E<sub>s,p</sub> and α is well represented by a power-law function, the relationship between F<sub>X, 0.3-7.0keV</sub> and E<sub>s,p</sub> is reasonably fitted with a linear function (All XRT: slope = 1.47 ± 0.11, intercept = 0.63 ± 0.04; All SXT: slope = 1.21 ± 0.07, intercept = 0.84 ± 0.05). In order to further test the correlations, we have performed a Spearman rank correlation analysis of these observable parameters for the spectra obtained during R1, ..., R6, and the results are reported in Table 7.

### 3.4. SEDs, Modeling, and Interpretation

The X-ray light curves shown in Figure 3 indicate the presence of clearly discernible, individual flares, in both 2016 and 2017. We interpret these events as the result of mildly relativistic shocks propagating through the jet of 1ES 1959+650. In order to model the light curves and SEDs during the period of our AstroSat observations in 2016 and 2017, we employ the time-dependent shock-in-jet model of Böttcher & Barknmh (2019). In this model, hybrid thermal + nonthermal electron distributions are generated via Monte Carlo simulations (Summerlin & Baring 2012) of diffusive shock acceleration by a mildly relativistic, oblique shock. As a representative choice of shock parameters, we assume a shock speed of $v_s = 0.71 c$ (in the comoving frame of the jet material), a magnetic field obliquity of $\Theta_{Bf1} = 32^\circ.3$, an upstream gas temperature of $5.45 \times 10^7$ K, and a compression ratio of $r = 3.71$ (see Böttcher & Barknmh 2019, for a motivation and discussion of these choices). The Monte Carlo simulations of diffusive shock acceleration parameterize the electrons' mean free path to pitch-angle scattering as $\lambda_{\text{pas}} = \eta_0 r_g p^{\alpha-1}$, where $r_g$ is an electron's Larmor radius and $p$ its momentum. $\eta_0$ and $\alpha$ are free parameters in the simulation. Since $r_g \propto p$, the mean free path scales as $\lambda_{\text{pas}} \propto p^{\alpha}$.
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Left: spectral slope $\alpha$ vs. integrated X-ray flux in the 0.3–7.0 keV band extracted from the time-resolved spectral modeling of SXT observations and the spectra from individual XRT observations starting from 2015-01-04 to 2021-02-15. The open black circles and open black diamonds refer to the SXT observations in 2016 and 2017, respectively. The other XRT observations are divided into seven time segments (R1,...,R7) and are shown by various symbols and colors as described in the legends. The dashed lines of different colors show the best-fit linear regressions using the EMCEE method. The various colored thin solid lines indicate the spread around the best-fit lines for a particular group of data. Right: best-fit linear correlations with data on $\alpha$ (top panel) and $\beta$ (without any fit—bottom panel) vs. flux. For clarity, the top panel does not include best-fit lines for the various time groups of XRT observations to illustrate the offset of the spectral index dependencies on flux between the flares in 2016 and 2017.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Diagnostics plots for the best-fit parameters for the $\alpha$ vs. flux correlations. Left: plot for the time-resolved SXT spectra from the 2016 flare. Right: plot for all the XRT spectra observed by Swift over the period MJD 57026–57637.</figcaption></figure>

The time-dependent radiative output from the resulting hybrid electron distributions is evaluated using the radiation transfer schemes of Böttcher & Chiang (2002) and Böttcher et al. (2013).

In addition to $\eta_0$ and $\alpha$, free parameters of the model are the shock-dissipated power transferred to relativistic electrons (termed “injection luminosity,” $L_{\text{inj}}$, in the following), the magnetic field
<table><caption>**Table 7**
Results from Correlation Studies</caption><thead><tr><th></th><th colspan="2">R1 = MJD 57025.0 - 57325.0</th><th colspan="2">R2 = MJD 57325.0 - 57625.0</th><th colspan="2">R3 = MJD 57625.0 - 57925.0</th></tr><tr><th></th><th>r<sub>s</sub></th><th>p</th><th>r<sub>s</sub></th><th>p</th><th>r<sub>s</sub></th><th>p</th></tr></thead><tbody><tr><td>α and F<sub>X</sub></td><td>-0.70</td><td>1.08 × 10<sup>-5</sup></td><td>-0.72</td><td>1.45 × 10<sup>-17</sup></td><td>-0.69</td><td>3.64 × 10<sup>-8</sup></td></tr><tr><td>E<sub>s,p</sub> and F<sub>X</sub></td><td>0.50</td><td>4.0 × 10<sup>-3</sup></td><td>0.79</td><td>2.64 × 10<sup>-23</sup></td><td>0.71</td><td>9.74 × 10<sup>-9</sup></td></tr><tr><td>E<sub>s,p</sub> and α</td><td>-0.79</td><td>1.26 × 10<sup>-7</sup></td><td>-0.85</td><td>6.02 × 10<sup>-29</sup></td><td>-0.95</td><td>1.24 × 10<sup>-26</sup></td></tr><tr><td>β and F<sub>X,0.3 - 7 keV</sub></td><td>-0.11</td><td>0.55</td><td>-0.40</td><td>3.37 × 10<sup>-5</sup></td><td>-0.25</td><td>8.19 × 10<sup>-2</sup></td></tr><tr><td>α and β</td><td>-0.10</td><td>0.59</td><td>0.04</td><td>0.69</td><td>0.15</td><td>0.28</td></tr><tr><td>m</td><td>-0.15 ± 0.03</td><td>...</td><td>-0.06 ± 0.006</td><td>...</td><td>-0.074 ± 0.01</td><td>...</td></tr><tr><td>b</td><td>2.20 ± 0.08</td><td>...</td><td>2.01 ± 0.03</td><td>...</td><td>2.19 ± 0.05</td><td>...</td></tr><tr><td>F<sub>var</sub></td><td>0.37 ± 0.009</td><td></td><td>0.41 ± 0.001</td><td></td><td>0.35 ± 0.002</td><td></td></tr><tr><td></td><td colspan="2">R4 = MJD 57925.0 - 58225.0</td><td colspan="2">R5 = MJD 58225.0 - 58525.0</td><td colspan="2">R6 = MJD 58525.0 - 58825.0</td></tr><tr><td></td><th>r<sub>s</sub></th><th>p</th><th>r<sub>s</sub></th><th>p</th><th>r<sub>s</sub></th><th>p</th></tr><tr><td>α and F<sub>X</sub></td><td>-0.87</td><td>1.80 × 10<sup>-19</sup></td><td>-0.68</td><td>1.75 × 10<sup>-9</sup></td><td>-0.37</td><td>1.43 × 10<sup>-2</sup></td></tr><tr><td>E<sub>s,p</sub> and F<sub>X</sub></td><td>0.84</td><td>7.55 × 10<sup>-17</sup></td><td>0.74</td><td>1.65 × 10<sup>-11</sup></td><td>0.47</td><td>1.71 × 10<sup>-3</sup></td></tr><tr><td>E<sub>s,p</sub> and α</td><td>-0.92</td><td>3.26 × 10<sup>-24</sup></td><td>-0.98</td><td>6.79 × 10<sup>-41</sup></td><td>-0.89</td><td>4.34 × 10<sup>-15</sup></td></tr><tr><td>β and F<sub>X,0.3 - 7 keV</sub></td><td>-0.13</td><td>0.32</td><td>-0.22</td><td>8.64 × 10<sup>-2</sup></td><td>0.02</td><td>0.92</td></tr><tr><td>α and β</td><td>0.05</td><td>0.71</td><td>-0.21</td><td>0.11</td><td>-0.56</td><td>1.27 × 10<sup>-4</sup></td></tr><tr><td>m</td><td>-0.057 ± 0.005</td><td>...</td><td>-0.05 ± 0.006</td><td>...</td><td>-0.03 ± 0.01</td><td>...</td></tr><tr><td>b</td><td>2.10 ± 0.04</td><td>...</td><td>2.15 ± 0.04</td><td>...</td><td>1.98 ± 0.06</td><td>...</td></tr><tr><td>F<sub>var</sub></td><td>0.36 ± 0.002</td><td></td><td>0.28 ± 0.001</td><td></td><td>0.30 ± 0.002</td><td></td></tr><tr><td></td><td colspan="2">R7 = MJD 58820.0 - 59260.0</td><td colspan="2">AstroSat - 2016 Oct., Nov.</td><td colspan="2">AstroSat - 2017 Nov.</td></tr><tr><td></td><th>r<sub>s</sub></th><th>p</th><th>r<sub>s</sub></th><th>p</th><th>r<sub>s</sub></th><th>p</th></tr><tr><td>α and F<sub>X</sub></td><td>-0.63</td><td>4.13 × 10<sup>-8</sup></td><td>-0.44</td><td>7.66 × 10<sup>-5</sup></td><td>-0.81</td><td>4.07 × 10<sup>-7</sup></td></tr><tr><td>E<sub>s,p</sub> and F<sub>X</sub></td><td>0.69</td><td>3.71 × 10<sup>-10</sup></td><td>0.52</td><td>2.04 × 10<sup>-6</sup></td><td>0.83</td><td>7.75 × 10<sup>-8</sup></td></tr><tr><td>E<sub>s,p</sub> and α</td><td>-0.90</td><td>1.09 × 10<sup>-23</sup></td><td>-0.77</td><td>6.79 × 10<sup>-16</sup></td><td>-0.68</td><td>9.71 × 10<sup>-5</sup></td></tr><tr><td>β and F<sub>X,0.3 - 7 keV</sub></td><td>0.068</td><td>0.59</td><td>-0.089</td><td>0.44</td><td>-0.11</td><td>0.58</td></tr><tr><td>α and β</td><td>-0.48</td><td>7.41 × 10<sup>-5</sup></td><td>-0.47</td><td>1.88 × 10<sup>-5</sup></td><td>-0.23</td><td>0.25</td></tr><tr><td>m</td><td>-0.09 ± 0.01</td><td>...</td><td>-0.04 ± 0.009</td><td>...</td><td>-0.05 ± 0.004</td><td>...</td></tr><tr><td>b</td><td>2.13 ± 0.05</td><td>...</td><td>2.01 ± 0.04</td><td>...</td><td>2.23 ± 0.05</td><td>...</td></tr><tr><td>F<sub>var</sub></td><td>0.26 ± 0.002</td><td></td><td>0.24 ± 0.003</td><td></td><td>0.087 ± 0.007</td><td></td></tr></tbody></table>

Note. Here $r_s$ is the Spearman rank coefficient, and $p$ denotes its null hypothesis probability. $m$ and $b$ are, respectively, the slope and the intercept of the best-fit EMCEE straight line between $α$ and $F_{X,0.3-7 keV}$. $F_{var}$ is the fractional variability amplitude.

strength $B$, the bulk Lorentz factor $Γ$ of the emission region's propagation along the jet, the viewing angle $θ_{obs}$, and the radius $R$ of the emission region, which is assumed spherical for the purpose of the radiation transfer simulations. $L_{inj}$, $B$, and $R$ are defined in the comoving frame of the emission region, while $θ_{obs}$ is the viewing angle in the observer's frame.

Since we find a satisfactory fit to snapshot SEDs and light curves of 1ES 1959+650 with synchrotron and synchrotron self-Compton (SSC) as the dominant radiation mechanisms (i.e., a pure leptonic SSC model), we do not consider any contribution from putative external radiation fields to the target photon field for Compton scattering to produce the γ-ray emission.

Our fitting procedure starts out with a quiescent-state configuration, reproducing the low-state SED of 1E 1959+650, shown by the black model curves in Figures 9(a) and 10(a), with model parameters listed in Table 8.

The zoomed-in light curves in Figures 9(b) and 10(b) clearly suggest that the observed X-ray flaring behaviors in both 2016 and 2017 cannot be modeled by a single, impulsive particle acceleration event, but rather are indicative of a succession of several shocks throughout the emission region. Specifically, in order to find a satisfactory representation of both the SXT and the LAXPC light curves during 2016, a succession of four shocks of different strength is required. All shocks are

<table><caption>**Table 8**
Parameters of Our Model Fits to the Quiescent-state SEDs of 1ES 1959+650 in 2016 and 2017, Respectively</caption><thead><tr><th>Parameter (units)</th><th>2016</th><th>2017</th></tr></thead><tbody><tr><td>η<sub>0</sub></td><td>60</td><td>40</td></tr><tr><td>α</td><td>1.9</td><td>1.8</td></tr><tr><td>L<sub>inj</sub> (erg s<sup>-1</sup>)</td><td>2.5 × 10<sup>40</sup></td><td>2.8 × 10<sup>40</sup></td></tr><tr><td>B (G)</td><td>0.15</td><td>0.08</td></tr><tr><td>Γ</td><td>20</td><td>20</td></tr><tr><td>θ<sub>obs</sub> (deg)</td><td>2.87</td><td>2.34</td></tr><tr><td>R (cm)</td><td>6.e15</td><td>1.e16</td></tr></tbody></table>

characterized by an increased injection luminosity and a global decrease of the pitch-angle mean free path, parameterized by a slightly smaller value of $η_0$. A consequence of this change is more efficient particle acceleration to higher energies, resulting in the observed larger variability amplitude in the LAXPC (3–30 keV) band compared to the SXT (0.3–7 keV) band. The parameters adopted for the four shocks in our 2016 simulation are listed in Table 9. Representative snapshot SEDs are shown in Figure 9(a), and the LAXPC and SXT light curves and model fits are presented in Figure 9(b). Since both the UV and Fermi-LAT γ-ray light curves consist of only very few points
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Synchrotron peak position and unabsorbed X-ray flux. Top: dependence of the spectral slope ($\alpha$; upper panel) and the curvature ($\beta$; lower panel) on the position of the synchrotron peak ($E_{x,p}$). The filled red stars represent the XRT data, whereas the measurements from the SXT are shown by filled black circles. Bottom: correlations between the estimated synchrotron peak positions and the integrated X-ray flux measured in the 0.3–7.0 keV band. The dotted red and black lines represent the best-fit linear correlations for XRT and SXT observations, respectively. The shaded areas around the best-fit curves in both panels represent the 3$\sigma$ confidence interval.</figcaption></figure>

(and upper limits), they are not constraining for our fits and are not shown in the figure. We have, however, verified that our model predictions are consistent with the data. Figures 9(c) and 9(d) show the predicted spectral hysteresis in a hardness-intensity diagram at various photon energies and the predicted cross-correlations between the light curves at different photon energies, respectively.

The observed hardness-intensity correlations plotted in Figure 6 (left) indicate a systematic offset by a factor of ~2 in flux of the data from 2017 with respect to the 2016 data in flux. However, within the individual yearly data sets, they show consistent harder-when-brighter trends. Such a systematic offset may be explained by a slight change of the Doppler

<table><caption>**Table 9**

Parameter Variations for the Fits to MWL Light Curves and SEDs in 2016</caption><thead><tr><th>Parameter (units)</th><th>L<sub>inj</sub> (erg s<sup>-1</sup>)</th><th>η<sub>0</sub></th><th>α</th></tr></thead><tbody><tr><td>Quiescence</td><td>2.5 × 10<sup>40</sup></td><td>60</td><td>1.9</td></tr><tr><td>Shock 1</td><td>3.0 × 10<sup>40</sup></td><td>50</td><td>1.9</td></tr><tr><td>Shock 2</td><td>3.5 × 10<sup>40</sup></td><td>50</td><td>1.9</td></tr><tr><td>Shock 3</td><td>4.1 × 10<sup>40</sup></td><td>40</td><td>1.9</td></tr><tr><td>Shock 4</td><td>3.4 × 10<sup>40</sup></td><td>50</td><td>1.9</td></tr></tbody></table>

<table><caption>**Table 10**

Parameter Variations for the Fits to MWL Light Curves and SEDs in 2017</caption><thead><tr><th>Parameter (units)</th><th>L<sub>inj</sub> (erg s<sup>-1</sup>)</th><th>η<sub>0</sub></th><th>α</th></tr></thead><tbody><tr><td>Quiescence</td><td>2.8 × 10<sup>40</sup></td><td>40</td><td>1.8</td></tr><tr><td>Shock 1</td><td>3.5 × 10<sup>40</sup></td><td>30</td><td>1.7</td></tr><tr><td>Shock 2</td><td>3.0 × 10<sup>40</sup></td><td>30</td><td>1.8</td></tr><tr><td>Shock 3</td><td>3.6 × 10<sup>40</sup></td><td>25</td><td>1.8</td></tr><tr><td>Shock 4</td><td>4.3 × 10<sup>40</sup></td><td>25</td><td>1.8</td></tr><tr><td>Shock 5</td><td>5.1 × 10<sup>40</sup></td><td>15</td><td>1.8</td></tr></tbody></table>

factor (by a factor of ~1.2) without any changes in the underlying particle acceleration and emission physics. A slight change of the viewing angle (from $\theta_{\text{obs}} = 2.0^{\circ}87$ in 2016 to 2.0$^{\circ}34$ in 2017) is sufficient to reproduce such a change in Doppler factor (from 20 to 24). This forms the basis of our modeling of SEDs and X-ray light curves in 2017. Only slight modification of other model parameters is required for our quiescent-state fit to the 2017 SEDs (see Table 8).

The X-ray light curves (Figure 10(b)) from the 2017 AstroSat observations show two clearly distinct major flaring episodes near the beginning and the end of the observations. Clearly, again, multiple shocks are required in order to provide a satisfactory representation of these complex light curves. Specifically, we find a good match with the SXT and LAXPC adopting a succession of five shocks, all characterized by a larger injection luminosity compared to the quiescent state, as well as more efficient particle acceleration to higher energies. The latter is primarily achieved through a change in the mean free path parameter $\eta_0$. The parameters adopted for the five subsequent shocks are listed in Table 10.

Given the relative simplicity of our model setup and the complexity of the SXT and LAXPC light curves in both 2016 and 2017, the correspondence between the observed X-ray light curves and our model predictions is remarkable. For both years, we reproduce the observed general harder-when-brighter trend (see Figure 6) with only very moderate spectral hysteresis, which is smaller than the error bars on the SXT and LAXPC flux and spectral index points, especially for 2016. The model 1 keV versus 10 keV time lags may be considered as a proxy to the observed trends between SXT and LAXPC. For 2016, our model predicts a lag of ~1 hr of the SXT (1 keV) light curve behind the LAXPC (10 keV) one, while for 2017 this lag is predicted to be slightly larger, at ~1.5 hr. Predicted lags between X-ray and $\gamma$-ray bands are equally of the order of $\lesssim 1$ hr.

## 4. Discussion and Summary

Our analysis, presented above, clearly establishes that IES 1959+650 underwent significant flaring activity in X-rays over the 6 yr from 2015 to 2021 (see Figure 1). The overall long-
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Model fits to SEDs and light curves of 1ES 1959+650 during the AstroSat observations in 2016. (a) Snapshot SEDs at various times during the simulation. (b) Model fits to the AstroSat light curves. (c) Predicted hardness-intensity correlations in various energy bands. Note that for plotting on the same scale, the 100 keV curve has been shifted up in flux (i.e., right) by a factor of 10, while the 1 GeV curve has been shifted up in spectral index (i.e., down in the plot) by 1. (d) Predicted discrete correlation functions between light curves in various energy bands.</figcaption></figure>

<figure><img src="image_11.png" /></figure>

term X-ray variability over these years shows a profile that is symmetric and well represented by a broad (FWHM ~ 785 days) Gaussian function peaking around 2018 April. Such a long-term symmetric flux variation is suggestive of a change in the viewing angle, and hence in the Doppler boosting, as the main driver. Most probably the inner jet is curved. Unfortunately, the data coverage is not sufficient to test the jet precession theory, which predicts periodic changes in the flux. The X-ray light curves show that there are many short-timescale variations superimposed on the long-term symmetric variations throughout (Figure 2).

The detailed light curves observed with AstroSat cover different prominent X-ray flares and exhibit shorter-timescale variations (see Table 5 for timescales and the variability amplitudes). During the X-ray outburst in 2016, the source also exhibits noticeable activity at γ-ray energies, even though no clear correlation pattern emerges. Most importantly, V2-SF1 is accompanied by simultaneous GeV/TeV activity, whereas V2-SF2 (total span of ~2 days) seems to be a pure orphan X-ray flare.

In 2017 1ES 1959+650 underwent several prominent X-ray outbursts, including the one observed with AstroSat, which seems to have a twin-flare-like profile within ~2 days (see Figure 2, right panels). During this period, 1ES 1959+650 broke its historical

X-ray flux record (in the 0.3–7.0 keV band), reaching a new maximum, which has not been exceeded until the time of writing (see Figure 1). The observed X-ray flux showed two high-amplitude, short-timescale flares without any counterpart at other wavelengths. Therefore, the variations observed with AstroSat are again orphan. Interestingly, a few days prior to the X-ray observations prominent variability took place at γ-ray energies. Unfortunately, the lack of multiwavelength observations during that time precludes any correlation analyses.

The AstroSat observations of the X-ray activity periods in 2016 and 2017 are unique, as these provide such a detailed variability profile for 1ES 1959+650 for the first time. The light curves in the 0.3–7.0 keV band are highly correlated with the ones in the 3.0–30 keV band for all the epochs in 2016 and 2017.

However, during all epochs, the variability amplitude is larger in the hard X-ray band compared to the soft X-rays. The shortest variability timescale ($\Delta t$; here characterized by the doubling or halving timescales) may be used to calculate an upper limit on the size of the emission region using the light-travel-time argument, $R \le \frac{c\delta\Delta t}{(1+z)}$. The smallest inferred limits on the size of the emission region for the activity periods in 2016 and 2017 are 2.9 and 7.4 Mpc, respectively.
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Same as Figure 9, but for the AstroSat observations in 2017. Note that in panel (c) the 1 GeV hardness-intensity curve has been shifted up (i.e., right) in flux by a factor of 10 and up in spectral index by 1 (i.e., down in the plot) for plotting on the same scale as the other curves.</figcaption></figure>

The X-ray spectral investigations reveal significant changes in the spectral shapes for different flares and also for the different segments of the AstroSat observations (T0, T1,..., T10). The spectral changes follow a harder-when-brighter trend as is typical for many blazars. The broadband X-ray spectra (0.3–30 keV) are best represented by a log parabola model. A similar investigation using time-resolved spectra from SXT and long-term observations from XRT provides similar trends. However, different sets of flares (R1,...,R7) show a slightly different relation between flux and spectral index: different R's populate significantly different tracks in the $\alpha-F_{X, 0.3-7\,\text{keV}}$ plane. The R segments representing the beginning (R1) and end (R7) of our data set, on the other hand, show extreme slopes (see Table 5 for details).

The time-resolved spectroscopy of the AstroSat observations reveals a strong correlation between the slope ($\alpha$) and the flux ($F_{X,0.3-7.0\,\text{keV}}$) for both observation periods in 2016 and 2017, similar to the enveloping R segments (R3 and R4, respectively). However, the tracks in the $\alpha$ versus $F_{X,0.3-7.0\,\text{keV}}$ plane of the AstroSat data sets are significantly different from the corresponding R segments. The best-fit lines of $\alpha$ versus $F_{X, 0.3-7.0\,\text{keV}}$ for the two sets comprise significantly different tracks, indicating that the spectrum exhibits stronger hardening with increasing flux in 2017 compared to 2016 (see Table 5). This indicates that the “bluer-when-brighter” trend was stronger in 2017 than in 2016, which is also supported by

the strong, positive correlation of the synchrotron peak energy ($E_{s,p}$) with the X-ray flux.

The various X-ray SEDs, namely, T1 to T6 from the flare in 2016 and T8 to T10 from the flare in 2017, are combined to generate two broadband SEDs. A model based on diffusive shock acceleration by mildly relativistic shocks in the jet of 1ES 1959+650 is able to simultaneously reproduce those snapshot SEDs and the AstroSat light curves of both flaring episodes in 2016 and 2017. In this model, multiwavelength flares are caused by shock-generated turbulence, leading to a reduction of the electrons’ mean free path to pitch-angle scattering ($\lambda_{\text{pas}}$) and, thus, more efficient particle acceleration. The instantaneous interplay between shock acceleration and self-consistent radiative (synchrotron, SSC, and external Compton) cooling of the particles in the emission region results in characteristic flux and spectral variability patterns, consistent with the observed ones. The different flux states between the two flares are well reproduced by a change in the Doppler factor, mediated by a slight change of the viewing angle ($\Delta\theta \sim 0.5^\circ$) within $\sim 1$ yr and a reduction of the magnetic field.

This manuscript highlights the flux evolution of various X-ray flares observed between 2015 December and 2021 February, their spectral properties, and their correlations. Additionally, a time-dependent leptonic model indicates the particle acceleration responsible for the two X-ray flares in 2016 and 2017 covered in detail with AstroSat. This work additionally avails a huge data set to the community with the spectral parameters spanning 6 yr
of X-ray monitoring of 1ES 1959+650. While a detailed study of the spectral parameters of the γ-ray bands, as well as their correlation with the X-ray bands, is not possible with the currently available data, such a study may shed light on the emission mechanisms during flaring activities. Future MeV/GeV/TeV facilities like AMIGO and CTA shall provide key data for such investigations.

The authors acknowledge the anonymous referee for the constructive suggestions. The work of M.B. and S.C. is supported by the South African Research Chairs Initiative (grant No. 64789) of the Department of Science and Innovation and the National Research Foundation<fnref n="20" /> of South Africa. P.G. acknowledges the financial support of Indian Space Research Organisation (ISRO) under AstroSat archival data utilization program. The Fermi Science Support Center (FSSC) team is acknowledged for providing the data and analysis tools. The authors acknowledge NASA's primary data archive “HEASARC” to avail the data from Swift. This research has used the data of the AstroSat mission of the Indian Space Research Organisation (ISRO), archived at the Indian Space Science DataCentre (ISSDC). The authors acknowledge the proposers and PIs of the dedicated follow-up programs and ToO observations of 1ES 1959+650 using Swift, without which such a detailed data set would not have been possible. The authors would like to acknowledge the support from the LAXPC Payload Operation Center (POC) and SXT POC at the TIFR, Mumbai, for providing supporting data reduction. LaxpcSoft software is used for analysis. This work has been performed utilizing the calibration databases and auxillary analysis tools developed, maintained, and distributed by the AstroSat-SXT team with members from various institutions in India and abroad.

**Facilities:** AstroSat (SXT and LAXPC), Swift (XRT and UVOT), FACT, MIRO, Fermi-LAT.

**Software:** astropy (Astropy Collaboration et al. 2013), emcee<fnref n="21" /> (Foreman-Mackey et al. 2013, 2019), extinction<fnref n="22" /> (Barbary 2016), extinctions<fnref n="23" />, fermipy<fnref n="24" /> (Wood et al. 2017), HEAsoft-v 6.25²⁵ (NASA High Energy Astrophysics Science Archive Research Center (Heasarc), 2014), laxpcsoft,<fnref n="26" /> lmfit<fnref n="27" /> (Newville et al. 2014), SEtractor (Bertin & Arnouts 1996), sxtARFModule<fnref n="28" />.

## Appendix

### Synchrotron Peak Frequency

The functional form of the mathematical model used to fit the X-ray spectra of 1ES 1959+650 is given by

$$ \frac{dN}{dE} = N_0 \left( \frac{E}{E_0} \right)^{-\left[ \alpha + \beta \times \log_{10} \left( \frac{E}{E_0} \right) \right]} \quad (A1) $$

where $E_0$ is known as pivot energy, $\alpha$ is the photon index, and $\beta$ is the curvature parameter of the model. The equivalent equation valid for fitting the $\nu F_\nu$ plot, i.e., the X-ray part of the SEDs, is

$$ \nu F_{\nu} = N_{0} \left( \frac{\nu}{\nu_{0}} \right)^{\left[ 2 - \alpha - \beta \times \log_{10} \left( \frac{\nu}{\nu_{0}} \right) \right]} \quad (A2) $$

It is well established that in most of the flux states the synchrotron peak of the broadband SED of 1ES 1959+650 lies between 0.3 and 10.0 keV. Hence, it is most likely that the maxima of Equation (A2) refer to the synchrotron peak energy. That is, $\nu_{syn,peak}$ or $\nu_{s,p}$ [$E_{s,p}$] corresponds to the solution of equation $\frac{d(\nu F_\nu)}{d\nu} = 0$,

$$ \frac{d(\nu F_\nu)}{d\nu} = N_0(2 - \alpha) \left(\frac{\nu}{\nu_0}\right)^{(1-\alpha)} \left(\frac{\nu}{\nu_0}\right)^{-\beta \times \log_{10}\left(\frac{\nu}{\nu_0}\right)} \\ - N_0 \left(\frac{\nu}{\nu_0}\right)^{(2-\alpha)} \times \left(\frac{2\beta \times \log_{10}\left(\frac{\nu}{\nu_0}\right)}{\nu}\right) \\ \times \left(\frac{\nu}{\nu_0}\right)^{-\beta \times \log_{10}\left(\frac{\nu}{\nu_0}\right)}, \quad (A3) $$

i.e.,

$$ \frac{d(\nu F_\nu)}{d\nu} = N_0 \left[ (2 - \alpha) - 2\beta \times \log_{10}\left(\frac{\nu}{\nu_0}\right) \right] \\ \times \left( \frac{\nu}{\nu_0} \right)^{\left[ 1, -\alpha, -\beta, \times \log_{10}\left(\frac{\nu}{\nu_0}\right) \right]} \quad (A4) $$

This means $\frac{d(\nu F_\nu)}{d\nu} = 0 = >$

$$ \nu_{s,p} = \nu_0 \times 10^{\left(\frac{2-\alpha}{2\beta}\right)} \quad \text{or; } \quad E_{s,p} = E_0 \times 10^{\left(\frac{2-\alpha}{2\beta}\right)}, \quad (A5) $$

where $\nu_0$ corresponds to the frequency of the photons at 1 keV, i.e., $\nu_0 = 2.4180 \times 10^{17}$ Hz. Therefore, for a given index $\alpha$ and curvature $\beta$ we can estimate the position of the synchrotron peak using Equation (A5). The same formalism is derived as in Equation (3) of Massaro et al. (2004).

## ORCID iDs

Sunil Chandra https://orcid.org/0000-0002-8776-1835

Markus Boettcher https://orcid.org/0000-0002-8434-5692

Pranjupriya Goswami https://orcid.org/0000-0001-5430-4374

Kulinder Pal Singh https://orcid.org/0000-0001-6952-3887

Michael Zacharias https://orcid.org/0000-0001-5801-3945

Navpreet Kaur https://orcid.org/0000-0002-7862-1056

Sudip Bhattacharyya https://orcid.org/0000-0002-6351-5808

Shashikiran Ganesh https://orcid.org/0000-0002-7721-3827

Daniela Dorner https://orcid.org/0000-0001-8823-479X

## References

Abdo, A. A., Ackermann, M., Agudo, I., et al. 2010, ApJ, 716, 30

Abdollahi, S., Acero, F., Ackermann, M., et al. 2020, ApJS, 247, 33

Aliu, E., Archambault, S., Arlen, T., et al. 2014, ApJ, 797, 89

Anderhub, H., Backes, M., Biland, A., et al. 2013, JInst, 8, P06008

Antia, H. M., Yadav, J. S., Agrawal, P. C., et al. 2017, ApJS, 231, 10

Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33

Barbary, K. 2016, Extinction, v0.3.0, Zenodo, doi:10.5281/zenodo.804967

Beckmann, V., Wolter, A., Celotti, A., et al. 2002, A&A, 383, 410

Bertin, E., & Arnouts, S. 1996, A&As, 117, 393

<footnote n="20">²⁰ Any opinion, finding, and conclusion or recommendation in this material is that of the authors, and the NRF does not accept any liability in this regard.</footnote>

<footnote n="21">²¹ https://emcee.readthedocs.io/en/stable/</footnote>

<footnote n="22">²² https://github.com/kbarbary/extinction</footnote>

<footnote n="23">²³ https://pypi.org/project/extinctions/</footnote>

<footnote n="24">²⁴ http://fermipy-readthedocs.io/en/latest/</footnote>

<footnote n="25">²⁵ https://heasarc.gsfc.nasa.gov/lheasoft/download.html</footnote>

<footnote n="26">²⁶ (http://astrosat-ssc.iucaa.in/?q=laxpcData)</footnote>

<footnote n="27">²⁷ https://lmfit.github.io/lmfit-py/</footnote>

<footnote n="28">²⁸ https://www.tifr.res.in/astrosat_sxt/dataanalysis.html</footnote>
Bhatta, G., Mohorian, M., & Bilinsky, I. 2018, A&A, 619, A93

Biland, A., Bretz, T., Buß, J., et al. 2014, JInst, 9, P10012

Biland, A., Dorner, D., Mirzoyan, R., et al. 2016a, ATel, 9148, 1

Biland, A. & FACT Collaboration 2016a, ATel, 9139, 1

Biland, A. & FACT Collaboration 2016b, ATel, 9239, 1

Biland, A., Mirzoyan, R. & FACT Collaboration 2016b, ATel, 9203, 1

Blandford, R. D., & Rees, M. J. 1978, PhyS, 17, 265

Bloom, S. D., & Marscher, A. P. 1996, ApJ, 461, 657

Böttcher, M., & Baring, M. G. 2019, ApJ, 887, 133

Böttcher, M., & Chiang, J. 2002, ApJ, 581, 127

Böttcher, M., & Dermer, C. D. 1999, ApL&C, 39, 129

Böttcher, M., Hivick, B., Dashti, J., et al. 2010, ApJ, 725, 2344

Böttcher, M., Reimer, A., Sweeney, K., & Prakash, A. 2013, ApJ, 768, 54

Buson, S., Magill, J. D., Dorner, D., et al. 2016, ATel, 9010, 1

Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, ApJ, 345, 245

Dai, Y., Wu, J., Zhu, Z.-H., Zhou, X., & Ma, J. 2011, AJ, 141, 65

Dermer, C. D., Schlickeiser, R., & Mastichiadis, A. 1992, A&A, 256, L27

Dondi, L., & Ghisellini, G. 1995, MNRAS, 273, 583

Dorner, D., Ahnen, M. L., Bergmann, M., et al. 2015, arXiv:1502.02582

Elvis, M., Plummer, D., Schachter, J., & Fabbiano, G. 1992, ApJS, 80, 257

Foreman-Mackey, D., Farr, W., Sinha, M., et al. 2019, JOSS, 4, 1864

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Ghisellini, G., Celotti, A., Fossati, G., Maraschi, L., & Comastri, A. 1998, MNRAS, 301, 451

Ghisellini, G., & Maraschi, L. 1989, ApJ, 340, 181

Goodman, J., & Weare, J. 2010, CAMCS, 5, 65

Green, G. M., Schlafly, E. F., Finkbeiner, D. P., et al. 2015, ApJ, 810, 25

Kalberla, P. M. W., Burton, W. B., Hartmann, D., et al. 2005, A&A, 440, 775

Kapanadze, B., Dorner, D., Vercellone, S., et al. 2016a, MNRAS, 461, L26

Kapanadze, B., Dorner, D., Vercellone, S., et al. 2018a, MNRAS, 473, 2542

Kapanadze, B., Romano, P., Vercellone, S., et al. 2016b, MNRAS, 457, 704

Kapanadze, B., Vercellone, S., Romano, P., et al. 2018b, ApJ, 858, 68

Kaur, N., Chandra, S., Baliyan, K. S., Sameer, & Ganesh, S. 2017, ApJ, 846, 158

Krawczynski, H., Hughes, S. B., Horan, D., et al. 2004, ApJ, 601, 151

MAGIC Collaboration, Acciari, V. A., Ansoldi, S., et al. 2020, A&A, 638, A14

Mannheim, K. 1993, PhRvD, 48, 2408

Mannheim, K., & Biermann, P. L. 1992, A&A, 253, L21

Massaro, E., Perri, M., Giommi, P., Nesci, R., & Verrecchia, F. 2004, A&A, 422, 103

Massaro, F., Tramacere, A., Cavaliere, A., Perri, M., & Giommi, P. 2008, A&A, 478, 395

Mastichiadis, A., & Kirk, J. G. 1997, A&A, 320, 19

Mücke, A., Protheroe, R. J., Engel, R., Rachen, J. P., & Stanev, T. 2003, APh, 18, 593

NASA High Energy Astrophysics Science Archive Research Center (Heasarc) 2014, HEAsoft: Unified Release of FTOOLS and XANADU, ascl:1408.004

Nellen, L., Mannheim, K., & Biermann, P. L. 1993, PhRvD, 47, 5270

Newville, M., Stensitzki, T., Allen, D. B., & Ingargiola, A. 2014, LMFIT: Non-Linear Least-Square Minimization and Curve-Fitting for Python, v0.8.0, Zenodo, doi:10.5281/zenodo.11813

Nishiyama, T. 1999, ICRC (Salt Lake City, UT), 3, 370

Patel, S. R., Shukla, A., Chitnis, V. R., et al. 2018, A&A, 611, A44

Perlman, E. S., Madejski, G., Georganopoulos, M., et al. 2005, ApJ, 625, 727

Perlman, E. S., Stocke, J. T., Schachter, J. F., et al. 1996, ApJS, 104, 251

Protheroe, R. J., & Mücke, A. 2001, in AIP Conf. Ser. 558, High Energy Gamma-Ray Astronomy: Int. Symp., ed. F. A. Aharonian & H. J. Völk (New York: AIP), 700

Richards, J. L., Max-Moerbeck, W., Pavlidou, V., et al. 2011, ApJS, 194, 29

Saito, S., Stawarz, Ł., Tanaka, Y. T., et al. 2013, ApJL, 766, L11

Schlafly, E. F., & Finkbeiner, D. P. 2011, ApJ, 737, 103

Schlegel, D. J., Finkbeiner, D. P., & Davis, M. 1998, ApJ, 500, 525

Shah, Z., Ezhikode, S. H., Misra, R., & R, R. T. 2021, MNRAS, 504, 5485

Singh, K. P., Stewart, G. C., Chandra, S., et al. 2016, Proc. SPIE, 9905, 99051E

Singh, K. P., Stewart, G. C., Westergaard, N. J., et al. 2017, JApA, 38, 29

Singh, K. P., Tandon, S. N., Agrawal, P. C., et al. 2014, Proc. SPIE, 9144, 91441S

Summerlin, E. J., & Baring, M. G. 2012, ApJ, 745, 63

Tagliaferri, G., Ravasio, M., Ghisellini, G., et al. 2003, A&A, 412, 711

Urry, C. M., & Padovani, P. 1995, PASP, 107, 803

Vaughan, S., Edelson, R., Warwick, R. S., & Uttley, P. 2003, MNRAS, 345, 1271

Wilms, J., Allen, A., & McCray, R. 2000, ApJ, 542, 914

Wood, M., Caputo, R., Charles, E., et al. 2017, ICRC (Busan), 301, 824

Yadav, J. S., Agrawal, P. C., Antia, H. M., et al. 2016, Proc. SPIE, 9905, 99051D