X-Ray Studies of Blazar 1ES 1959+650 Using Swift and XMM-Newton Satellite

High synchrotron energy peaked blazar 1ES 1959+650 was studied with the Swift and XMM-Newton satellites in a total of 127 observations during the period 2018 June–2020 December. We extensively studied its flux and spectral variability on intraday and long-term timescales. Discrete correlation function analysis between soft and hard X-ray bands indicates soft as well as hard lags. The results were used to constrain the magnetic field of the emitting region, which was found to be 0.64 ± 0.05 G. On long-term timescales, the distribution of fluxes shows lognormal behavior, which could be attributed to the minijets-in-a-jet model or might be due to the propagation of relativistic shocks down the jet. The spectral energy distribution around the synchrotron peak is well described by the log-parabola model. Spectral parameters like the peak energy E p , curvature β, and peak luminosity L p were derived from spectral analysis. Their correlations were studied to constrain the acceleration processes of the emitting particles. E p shows strong correlation with L p during the high state of the source, which indicates spectral changes might be caused by the variations of the average electron energy. Low values of the curvature parameter β and a weak correlation between E p and β indicate the coexistence of stochastic/statistical acceleration of electrons in the emitting region. Implications of other results are also discussed.


Introduction
Blazars, a subclass of active galactic nuclei, comprise BL Lacertae objects and flat-spectrum radio quasars (FSRQs). Blazars are characterized by rapid flux variability at timescales ranging from a few minutes to years across the entire electromagnetic spectrum as well as by high optical polarization, and the emission is predominantly nonthermal in nature, which emanates from a relativistic jet streaming along, or aligned very close to, our line of sight (Blandford & Rees 1978).
The spectral energy distribution (SED) of blazars consists of a double-peaked hump with the first peak in the submillimeter to soft X-ray bands and the second peak 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. However, the physical mechanisms behind the high-energy emission are not well understood and are thought to originate from the Compton upscattering of synchrotron photons by the same population of relativistic electrons (synchrotron self-Compton (SSC); Ghisellini & Maraschi 1989;Mastichiadis & Kirk 1997) or by external seed photons originating from the accretion disk, broad-line region, and torus components of the blazar (external Compton; Dermer et al. 1992;Ghisellini et al. 1998). The other models that appear to be viable mechanisms to explain the X-ray through γ-ray emission are hadronic models where the highenergy emission is produced by relativistic protons through proton synchrotron radiation and photopion production, followed by pion decay and electromagnetic cascades (e.g., Mannheim & Biermann 1992;Böttcher et al. 2013 and references therein).
Blazars are divided into three classes based on the position of the synchrotron peak in their SED (Padovani & Giommi 1995). Low-frequency BL Lac objects exhibit a synchrotron peak in the IR-optical band, intermediate-frequency BL Lac objects have their synchrotron peaks at optical-UV frequencies, and highfrequency BL Lac objects show a synchrotron peak in the UV to X-ray band.
The blazar 1ES 1959+650 is at a redshift of z = 0.048 (Perlman et al. 1996) and is a prominent high synchrotron peaked (HSP) blazar in which the synchrotron peak of the SED appears in the UV-X-ray band (Krawczynski et al. 2004;Abdo et al. 2010;Kapanadze et al. 2016b). In X-rays, it was first detected by the Slew Survey with the Einstein Imaging Proportional Counter (Elvis et al. 1992), followed by BeppoSAX (Beckmann et al. 2002), RXTE, Swift, XMM-Newton (Tagliaferri et al. 2003;Massaro et al. 2008), and the Nuclear Spectroscopic Telescope Array (NuSTAR) (Pandey et al. 2017) in later years.
It has shown strong flux variability in the optical, X-ray, and TeV energy bands (Krawczynski et al. 2004;Kapanadze et al. 2016bKapanadze et al. , 2018aKapanadze et al. , 2018bKaur et al. 2017;Patel et al. 2018;Wang et al. 2018Wang et al. , 2019MAGIC Collaboration et al. 2020). High X-ray flaring activity of the source has been reported by Perlman et al. (2005) and Krawczynski et al. (2004) using XMM-Newton and RXTE-PCA observations in 2002. Kapanadze et al. (2016b reported frequent X-ray flares of this source during 2005-2014 using Swift-XRT observations. This blazar, 1ES 1959+650, underwent an unprecedented X-ray flaring activity during 2015 August-2016 January. During this period, it varied by a factor of ∼5.7 with a maximum value above 20 counts s −1 , along with high flux activity in the TeV energy band (Kapanadze et al. 2016b;Kaur et al. 2017;Patel et al. 2018). However, in several multiwavelength campaigns, orphan flares in γ-rays (which are not simultaneous with X-ray flares) were found in 2002 June (Krawczynski et al. 2004) and these flares cannot be explained with conventional one-zone SSC models. They invoke multiple-component SSC models; external Compton models where the variations of the external photon intensity in the jet frame can cause orphan γ-ray flares; or a magnetic field aligned along the jet axis, in which case the observer would not see the synchrotron flare, but the electrons would scatter SSC γ-rays in the observer's direction and thus the observer would be able to see the inverse Compton flare. Bottcher (2005) proposed that orphan flares are difficult to reconcile with the standard leptonic SSC model and suggested that they may originate from relativistic protons interacting with an external photon field supplied by electron synchrotron radiation reflected off a dilute reflector. Chandra et al. (2021) modeled the SED of this blazar using a single-zone time-dependent SSC model, which could explain the flares successfully. Patel et al. (2018) studied this source during observation period 2016 June-July and explained its broadband SED using a two-zone SSC model where the inner zone is mainly responsible for producing the synchrotron peak and the high-energy γ-ray part, whereas the second zone explains the less variable optical-UV and lowenergy γ-ray emission.
The X-ray spectral index hardens with increasing flux level in the long-term duration (Patel et al. 2018;MAGIC Collaboration et al. 2020) and also during a number of flares (Wang et al. 2018). Kapanadze et al. (2016bKapanadze et al. ( , 2018aKapanadze et al. ( , 2018b also showed the "harder-when-brighter" trend in the blazar 1ES 1959+650. Recently, Chandra et al. (2021) presented an extensive analysis of 1ES 1959+650 from the period 2016-2017 to 2021 February using the X-ray data from AstroSAT and Swift and found that the synchrotron peak shifts significantly with different flux states. Shah et al. (2021) also used AstroSAT observations to study the anticorrelation between the photon index and the X-ray flux using a broken power law.
Around the synchrotron peak, the SED is curved and can be well described by a log-parabolic model (e.g., Landau et al. 1986;Massaro et al. 2004;Tramacere et al. 2007;Chen 2014;Gaur et al. 2018;Pandey et al. 2018). It is characterized by the peak energy (E p ), peak luminosity (L p ), and curvature parameter β. The log-parabolic spectral distribution arises when the acceleration probability is a decreasing function of the electron energy. After analyzing large X-ray observations, Massaro et al. (2004) and Tramacere et al. (2007Tramacere et al. ( , 2011 suggested that the observed anticorrelation expected between E p and β could be used as a signature of a stochastic component in the acceleration process. For the blazar population, an apparent anticorrelation is expected between E p and L p (e.g., Tramacere et al. 2007;Massaro et al. 2008;Kapanadze et al. 2016bKapanadze et al. , 2017Kapanadze et al. , 2018a, which might be associated with the change in average electron energy, beaming factor, or magnetic field (e.g., Tramacere et al. 2009).
In the present work, we analyzed Swift observations of 1ES 1959+650 during the period 2018 June-2020 December on 125 nights. We also studied two XMM-Newton satellite observations available during this period. Our aim was to study the temporal/spectral variability of this source on timescales of minutes to years covering different flux states.
Studying flux and spectral variability during different flux states is important as distinct physical processes may play a dominant role in different flux states. Flux and spectral variations of blazars on diverse timescales arise either from a purely intrinsic phenomenon, such as the interaction of relativistic shocks with particle density or magnetic field irregularities in the jet (e.g., Marscher 2014) or the production of minijets-in-a-jet (e.g., Giannios et al. 2009), or from extrinsic mechanisms. Extrinsic mechanisms involve the geometrical effects that result from the bending of the jets, either through instabilities (e.g., Pollack et al. 2016) or through orbital motion (e.g., Valtonen & Wiik 2012;Fromm et al. 2013;Larionov et al. 2020). Long-term variability is generally attributed to a mixture of intrinsic and extrinsic mechanisms, which include shocks propagating down twisted jets or relativistic plasma blobs moving downstream a helical structure in magnetized jets (e.g., Marscher et al. 2008). The SEDs of blazars are explained by leptonic and hadronic models (Böttcher et al. 2013) and flux and spectral variability can be explained by merely changing the SED parameters adopting a common set of physical mechanisms (e.g., Patel et al. 2018;Prince et al. 2019;Ghosal et al. 2022). We fit the spectra using a log-parabolic model and derived the spectral parameters (i.e., E p , L p , and β) that varied with different flux states of the source. In this work, we analyzed correlations between these spectral parameters during different flux states of the source, which could provide tight observational constraints on the acceleration and injection processes of the emitting electrons.
The paper is structured as follows. Section 2 describes the observation and data analysis procedures of the Swift and XMM-Newton satellites. Section 3 provides the analysis techniques used to quantify variability, variability timescales, and power spectral density (PSD) as well as the spectral analysis and various models used in the study. In Section 4, we describe the results and their interpretation. The results are discussed in Section 5 and summarized in Section 6.

Swift
The Neil Gehrels Swift Observatory is a multiwavelength facility equipped with the X-Ray Telescope (XRT), the Burst Alert Telescope, and the Ultraviolet/Optical Telescope (UVOT) (Gehrels et al. 2004). We retrieved the Swift-XRT (Burrows et al. 2005) data from the publicly available HEASARC data archive. Data reduction was performed with the XRT Data Analysis Software (v.3.6.0), which is a part of the HEASoft package (v.6.28). Level 1 unscreened event files were reduced, calibrated, and cleaned with the use of the XRTPIPELINE script (version 0.13.5). The latest calibration files of Swift CALDB were used with remote access. 3 All the observations were in Windowed Timing mode. Events were selected with standard filtering criteria of 0-2 grades for the Windowed Timing observations. Sources with center pixels lying inside a two-pixel radius of bad pixels were not used in the analysis. The source region was extracted with a circular region of a 20-pixel radius. For the Windowed Timing mode, the source should appear in the middle of the 1D image; therefore the background can be selected in regions on either side of the source.
XRTPRODUCTS was used to obtain the light curve and spectrum of the source and background. The obtained source light curve was corrected for resultant loss of the effective area, bad/hot pixels, and vignetting with the use of the XRTLCCORR task. Ancillary response files were created using the XRTMKARF task. Source spectra were binned to ensure at least 20 counts per bin in order to use the χ 2 fitting method.

XMM-Newton
Blazar 1ES 1959+650 was observed on 2019 July 5 and 2020 July 16 by the XMM-Newton satellite (Jansen et al. 2001). We used the European Photon Imaging Camera (EPIC)pn instrument data. EPIC-pn is the most sensitive photonimaging camera and the least affected by photon pileup effects (Strüder et al. 2001). Data reduction was performed with the use of the XMM-Newton Science Analysis System (SAS) for the light-curve extraction. We extracted the high-energy (10 keV < E < 12 keV) light curve for the full frame of the exposed CCD in order to identify the flaring particle background. We restricted our analysis to the 0.3-10 keV energy range, as data below 0.3 keV are markedly contaminated by noise events and data above 10 keV are usually dominated by background flares. The source region was extracted using a circle of 40″ radius centered on the source. The background light curve was obtained from the region that corresponds to a circular annulus centered on the source with an inner and outer radius of 50″ and 60″, respectively. Pileup effects were examined using the SAS task EPATPLOT. We found a pileup in our observations. In order to remove the pileup, the central 7.5″ radius region was removed during the extraction of the light curve. Source light curves were obtained for the 0.3-10 keV energy band (corrected for background flux and given in units of counts s −1 ), sampled with a fixed bin size of 0.5 ks.

Excess Variance
Excess variance (s XS 2 ) is a measure of source intrinsic variance (Edelson et al. 2002;Vaughan et al. 2003a) evaluated by subtracting the variance that arises from measurement errors (s err 2 ) from the total variance of the observed light curve (S 2 ). If a light curve consists of N measured flux values x i with corresponding finite uncertainties σ err,i arising from measurement errors, then the normalized excess variance (s NXS 2 ) is calculated as follows:¯( is the mean square error, and S 2 is the sample variance of the light curve, given by The fractional rms variability amplitude, F var (Edelson et al. 1990;Rodríguez-Pascual et al. 1997), which is the square root of s NXS 2 , is thus¯( The uncertainty on F var is given by Vaughan et al. (2003a):

Hardness Ratio
The hardness ratio (HR) is useful for characterizing spectral changes over a broad X-ray energy range (e.g., Sivakoff et al. 2004;Park et al. 2006). The energy ranges of (0.3-2) keV and (2-10) keV are used here as soft and hard bands, respectively. HR is then defined as and the error in HR (σ HR ) is calculated as follows: where S and H are the net count rates in the soft (0.3-2 keV) and hard (2-10 keV) bands, respectively, while σ S and σ H are their respective errors (e.g., Pandey et al. 2017).

Doubling/Halving Timescales
The characteristic halving/doubling timescale τ, depending on the increase or decrease of the flux, is the shortest flux variability time, which is calculated from where F(t 1 ) and F(t 2 ) are the fluxes of the light curve at times t 1 and t 2 , respectively. We consider the timescales when the difference in flux is significant at the 3σ level (e.g., Foschini et al. 2011;Dhiman et al. 2021).

Discrete Correlation Function
The discrete correlation function (DCF) is calculated for light curves of two different energy bands. We consider here 0.3-2 keV and 2-10 keV as the two energy bands. The function is used to investigate a correlation between two unevenly sampled time series data. For two such discrete data sets a i and b j , the unbinned discrete correlation is defined as (Edelson & Krolik 1988 for all measured pairs (a i , b j ) with the pairwise lag ΔT ij = t j − t i ; e a andē a , b and s b , and a and σ b are the measurement error, mean, and standard deviation associated with data set a i and b j , respectively. DCF(τ) is obtained by averaging the number N of pairs lying in the range ij DCF evaluates the cross-correlation and possible time lags between the soft and hard energy bands. The obtained DCF is fitted with the Gaussian function (Edelson & Krolik 1988) of the form where τ lag is the time lag at which DCF peaks and σ is the width of the Gaussian function (as used in Gaur et al. 2015).

PSD
PSD is the measure of variability power as a function of temporal frequency. It is used to characterize temporal variations in flux and noise processes in general. A periodogram is a tool that is used to find hidden periodicities including any quasiperiodic oscillations. It is defined as the modulus-squared of the discrete Fourier transform (DFT) of the data (Vaughan 2005). The obtained power spectrum consists of red noise, which dominates over the measurement error, i.e., the Poisson noise at lower frequencies, and white noise, which dominates red noise at higher frequencies and becomes equivalent to the Poisson noise level of the data. The red noise part of the power spectrum is then fitted with a model of the form P( f ) = N f −m , where N is the normalization and m is the power-law spectral index (m > 0) (van der Klis 1989; González-Martín & Vaughan 2012). The equations used in determining PSDs are provided in the Appendix.

Spectral Analysis
The Xspec software package version 12.11.1 is used for spectral fitting. The Galactic absorption n H is fixed to 1.0 × 10 21 cm −2 (Lockman & Savage 1995) and the Xspec routine cflux is used to obtain unabsorbed flux and its error. Massaro et al. (2004Massaro et al. ( , 2008 found that blazars' spectra are curved, which arises from log-parabolic electron distributions. Therefore, they are well described by the log-parabola model (e.g., Tramacere et al. 2007Tramacere et al. , 2009. We fit each spectrum using models defined as follows: 1. Power-law model, which is defined by k E − α . It is characterized by the photon index α, redshift z, and normalization k. 2. The log-parabolic model, logpar. It is characterized by the photon index α, curvature β, and normalization k. 3. Another form of the log-parabolic model, i.e., the eplogpar model, is used to calculate the synchrotron peak E p . Details of the equations are provided in the Appendix.

Results
We studied 125 archived observations of the Swift satellite of the TeV blazar 1ES 1959+650 during the period 2018 June-2020 December. We also analyzed two publicly archived XMM-Newton observations of this blazar, which were made on 2019 July 5 and 2020 July 16. These observations were used to study the flux and spectral variability of this blazar on intraday as well as long-term timescales.

Intraday Flux and Spectral Variability
Blazar 1ES 1959+650 was studied with the XMM-Newton observations made on 2019 July 5 and 2020 July 16. The obtained light curves are shown in Figure 1. Flux variability was calculated using excess variance and the amplitude was calculated to be 1.95% and 3.12%, respectively. We also analyzed Swift-XRT observations of this source during the period 2018 June-2020 December on a total of 125 nights. The sample of light curves are shown in Figure 2. We found significant flux variability (i.e., F var > 3σ) in five of these observations with amplitude varying between 4.11% and 7.34%, which are presented in Table 1.
In order to calculate the spectral variability of these observations, HR analysis was performed, the results of which are presented in Table 2. HR analysis yielded no spectral variability for any individual light curve, which can be seen from the plots of HR versus counts per second (0.3-10 keV) shown in Figure 2.

Variability Timescale
The doubling/halving timescale τ var defined in Equation (7) was used to calculate the shortest variability timescales. The shortest variability timescale we found was t var = 15.28 ks. The emission region size was constrained with the equation

Cross-correlated Variability and PSD
Cross-correlation studies were performed using the XMM-Newton data as they are long observations with high cadence. The cross-correlation between the soft band (0.3-2 keV) and the hard band (2-10 keV) was obtained using DCF and the corresponding DCF plots are shown in the rightmost panels of Figure 1. The DCF plots were fitted with a Gaussian function (as described in Section 3.4) and we obtained a time lag of  −0.94 and 0.36 ks for the observations performed on 2019 July 5 and 2020 July 16, respectively. The results of DCF are provided in Table 3. Significant correlation at positive/negative lags means that soft/hard variations are leading the variations in the hard/soft bands, respectively. During the observation performed on 2019 July 5, the lags are negative, i.e., the variations at (2-10) keV are leading those at (0.3-2) keV. Therefore, in this case, variations at lower energies are slower than variations at higher energies. The reverse situation is seen during the second observation. We found that there is a hard lag, i.e., the (2-10) keV band is lagging behind the soft one (0.3-2 keV). The X-ray emission of HSP blazars lies at the top of the synchrotron hump and is characterized by the energy-dependent acceleration and cooling mechanisms.
Following Zhang et al. (2002), the acceleration timescale t acc and cooling timescale t cool of relativistic electrons in the observed frame can be expressed as a function of the observed photon energy E (in keV): where z is the source's redshift, B is the magnetic field in gauss, δ is the Doppler factor of the emitting region, and ξ is the parameter describing how fast the electrons can be accelerated. One can note that t acc and t cool both depend on the photon energy in inverse fashion. Higher-energy electrons cool faster compared to lower-energy electrons but accelerate slower than lower-energy electrons. Therefore, if t cool is greater than t acc , the cooling process dominates (Kirk et al. 1998). In such a case, higher-energy photons will lead to lower-energy photons and a soft lag is expected: If t acc is comparable to t cool in the observed energy range, acceleration processes dominate in the emitting region and a hard lag is expected. The time lag in an acceleration-dominated system is expressed as Using the observed lags, one can discern the physical parameters of the emitting region as follows: where τ hard and τ soft refer to the observed hard and soft lags (in seconds) between the low-energy E l and high-energy E h bands (in keV), respectively (E l and E h are the logarithmic-averaged energies of the given energy bands). Equation (17) was used to calculate the magnetic field with redshift z = 0.048, τ soft = 940 s, and the Doppler factor δ = 15 (Patel et al. 2018). The magnetic field B of the emitting region was found to be 0.64 ± 0.05 G, which is consistent with values provided in the literature. Tagliaferri et al. (2003) andMAGIC Collaboration et al. (2020) performed SED modeling and obtained magnetic field values close to the value we obtained.

Long-term Flux and Spectral Variability
We studied the long-term flux and spectral variability of blazar 1ES 1959+650, observed within the period 2018 June-2020 December. Based on the flux state of the blazar, we divided the total epoch into three periods. The flux state of the source was defined with respect to the average counts during our observing run.  Figure 3. The light curves for different periods are shown in Figure 4. During long-term timescales, the source shows significant spectral variability. The HR, α, β, and E p show long-term variations, which can be seen in Figure 4 for different periods.
In Period 1, the source exhibits two flares with flux that reaches up to 128.82 × 10 −11 erg cm −2 s −1 on 2018 September 16. The fractional variability amplitude F var is found to be ∼29.96% and significant spectral variability is also found during this period. The photon index hardens with an increase in flux of the source, varying between 2.17 and 1.58. Lower values of the curvature parameter are found, which range between 0.84 and 0.31. The peak energy shows a positive correlation with the flux of the source and reaches up to 2.45 keV. During Period 2, the source exhibits small flares with the mean flux level of this period reaching up to 52.93 × 10 −11 erg cm −2 s −1 . F var is found to be 23.12% and significant spectral variability is also found during this period, which can be seen in Figure 4. The photon index hardens with an increase in flux of the source, varying between 2.25 and 1.56. As in Period 1, β has lower values, which range between 0.82 and 0.23. During Period 3, the source exhibits a relatively low flux state with the highest flux reaching up to 61.66 × 10 −11 erg cm −2 s −1 . F var is found to be 21.43% and significant spectral variability is found during this period. The photon index hardens with an increase in flux of the source, varying between 2.14 and 1.34. β has lower values ranging between 0.99 and 0.24. E p is positively correlated with flux and reaches up to 2.88 keV.

Lognormality of Flux Distributions
The nature of variability in a long-term period can be quantified or explained with the flux distribution of the variable source (Uttley et al. 2005). Blazars exhibit a lognormal flux  Narayan & Piran 2012). As blazars have strongly magnetized jets pointing toward us, flux distributions using the minijets-ina-jet model (Giannios et al. 2009) were studied by Biteau & Giebels (2012) and they found that the flux from a single randomly oriented minijet will follow a Pareto distribution. The flux integrated from many isotropically oriented minijets could lead to an α-stable distribution, which could converge to a lognormal distribution when subjected to experimental uncertainties. Therefore, nonlinear flux distributions can arise from small Gaussian perturbations. This could provide an explanation for the lognormal flux distributions in blazars during flaring states. In this scenario, the flux distribution has been found to follow the rms-flux relation (Biteau & Giebels 2012).
An alternative interpretation for the non-Gaussian distribution of blazar variability light curves is provided by Sinha et al. (2018). They explained it using a small perturbation in the acceleration timescale, which can result in the variability of the particle number density, which is a linear combination of Gaussian and lognormal processes. The dominant shape of the resultant flux distribution is determined by the relative weight of these processes Bhatta & Dhital 2020). They also demonstrated that a perturbation in the acceleration timescale leads to a Gaussian distribution of its photon index, whereas a perturbation in the particle cooling rate produces neither of the two distributions (Shah et al. 2018Sinha et al. 2018;Khatoon et al. 2020Khatoon et al. , 2022, and references therein).
The flux distribution of 1ES 1959+650 was studied by Patel et al. (2018) using radio to γ-ray data. It was also studied by Bhatta & Dhital (2020) using decade-long Fermi/LAT observations. Duda & Bhatta (2021) used maximum likelihood estimation methods to study flux distributions.
In order to investigate the lognormality of flux distributions in our observations, we fit histograms of the X-ray data observed by Swift-XRT within the period 2018 June to 2020 December with Gaussian and lognormal distributions. Figure 5 shows the lognormal and normal flux distribution of the blazar for the total epoch used in this analysis. We used the Anderson-Darling test (e.g., Anderson & Darling 1952, 1954Stephens 1977;D'Agostino & Stephens 1986;Jäntschi & Bolboacă 2018;Shah et al. 2018) to quantify the nature of the flux distribution of blazar 1ES 1959+650. We obtained that the blazar 1ES 1959+650 follows a lognormal behavior in the total epoch with a p-value = 0.20. We also found a significant linear correlation between the rms and flux for the total epoch, which indicates that the variability might arise from the minijets-in-ajet model. However, as recently shown by Scargle (2020), the linear rms-flux relationship can be obtained from intrinsically  additive processes; therefore this result might be used with caution. We also calculated the photon index distribution for the total epoch and found it to be well fitted with lognormal as well as normal distributions. As discussed in Sinha et al. (2018), small temporal fluctuations in the intrinsic timescales in the acceleration region are capable of producing particle distributions with non-Gaussian signatures and significant flux-rms correlations. Therefore, we cannot rule out the possibility of an acceleration-due-to-shock scenario.

Relation between Spectral Parameters
Blazar 1ES 1959+650's spectra were fit by log-parabolic and power-law models. The results of both models are presented in Table 4. We used the F-test to compare the fitting results of these two models. We found that all the spectra are well fitted by a log-parabolic model. Then, we derived the spectral fitting parameters, i.e., the location of the synchrotron peak (E p ) and the peak luminosity (L p ), with the log-parabolic model. The results show the variation of E p in the range of 0.46-2.88 keV whereas L p varies in the range of 0.51-2.50 erg cm −2 s −1 . Spectral variations of the photon index appear between 1.34 and 2.25 and the curvature (β) varies between 0.21 and 0.99. Correlations between the various spectral parameters are shown in Figure 6 and are studied with Spearman's rank correlation coefficient ρ s and their corresponding p-values. The results are presented in Table 5.
A positive correlation is expected between the photon index α and curvature β in the first-order Fermi acceleration scenario. This correlation is predicted for the energy-dependent acceleration probability (EDAP) process (Massaro et al. 2004), where the probability p i that a particle undergoes an acceleration step i, with the corresponding energy g i q and energy gain ε, is given by p i = g g i q , where g and q are positive constants. Therefore, as the energy of the particle increases, the probability of the particle's acceleration decreases. According to Massaro et al. (2004), a linear relationship is expected between the spectral index s and curvature r, . The synchrotron emission produced by a differential energy spectrum of the form ( )  Massaro et al. 2004). In our observations, we found a weak negative correlation between α and β, which is expected when g > γ 0 (Kapanadze et al. 2020a). It implies that there exists an electron population with a very low initial energy γ 0 in the emission zone. Kapanadze et al. (2020a) found a negative correlation between these quantities for Mrk 421 for some of its observational period during 2015 December-2018 April. The coexistence of second-order Fermi acceleration/ stochastic acceleration could also weaken the α-β correlation. Katarzyński et al. (2006) have shown via simulations that electrons can be accelerated at the shock front via EDAP but can gain energy via the stochastic mechanism after escaping the shock front. The combined effect of both processes could result in a weak or merely modest correlation between α and β (Kapanadze et al. 2016b(Kapanadze et al. , 2020aKapanadze 2018). The synchrotron peak energy (E p ) and the luminosity (L p ) of a source follow a power-law relation of the form µ L E p p a (Rybicki & Lightman 1979). If the electron distribution in the emitting region follows a log-parabolic distribution, the peak luminosity is given by p p 2 (e.g., Tramacere et al. 2009). The parameter γ p represents the peak of n(γ)γ, where γ is the electron Lorentz factor; N ∼ n(γ p )γ p is the total electron number; B represents the magnetic field; and δ is the Doppler beaming factor. If a = 1, the spectral changes are mainly caused by the variations of the average electron energy, but the total electron number remains constant; if a = 1.5, they are mainly caused by the variations of the average electron energy, but the total electron number also changes; if a = 2, they are correlated with the changes of the magnetic field, B; and if a = 4, they might be dominated by the variations of the beaming factor, δ. We found significant correlations between E p and L p during our observations. We fit our observational data using the equation and found a to be <1 during the whole observational period. During Period 1, a ∼ 1, which confirms that the spectral changes might be caused by the variations of the average electron energy.
The correlation between E p and β provides us clues about the acceleration mechanism (e.g., Massaro et al. 2004;Tramacere et al. 2011), such as whether it is a statistical/stochastic mechanism. These two mechanisms are produced by a logparabolic electron distribution, resulting in a log-parabolic SED.
In the statistical acceleration process, the electron energy distribution follows the log-parabolic law and the acceleration efficiency of the emitting electrons is inversely proportional to their energy (Massaro et al. 2004). Therefore, in such a process, E p and β follow a correlation of the form Const. 2 5 p , with the assumption of β = r/4, where r is the curvature of the electron energy distribution (Chen 2014). For fluctuations of the fractional acceleration gain process, electron energies follow the lognormal law, and energy gain fluctuations are a random variable around the systematic energy gain (Tramacere et al. 2011). In such cases, E p and β follow the correlation Const. 3 10 p given β = r/4 (Chen 2014). The second scenario is described by the stochastic acceleration process. Here, a momentum diffusion term is included in the kinetic equation, which leads to energy gain fluctuations in the diffusive shock acceleration process (Tramacere et al. 2011). In this process, E p and β follow the relation Const. 1 2 p , where β = r/4 (Chen 2014). Therefore, the theoretically expected values of C are 10/3, 5/2, and 2 for the fractional acceleration gain fluctuation, EDAP, and stochastic acceleration processes, respectively. We found a significant negative correlation between E p and 1/β in our observations. We fit our observational data using the equation 1/β = + C E D log p but did not get a value of C close to the above values; hence we cannot explain our observational data using any of the above acceleration mechanisms. However, the coexistence of the above acceleration mechanisms might be possible (Wang et al. 2019), which could lead to an overall weakening of the correlation.
A positive correlation between flux and E p was found during our observations, which indicates a shift of the synchrotron peak to higher energies. Also, we found significant correlation between flux and S p , which is expected as peak height increases as flux increases (Holder et al. 2003;Kapanadze et al. 2016aKapanadze et al. , 2016bKapanadze et al. , 2018aKapanadze et al. , 2018bWang et al. 2019;Chandra et al. 2021). Near the peak energy of the emission, the cooling timescale shortens and can complete with the acceleration timescales (Tramacere et al. 2009), which leads to an anticorrelation between E p and β if the cooling timescale is shorter than that of EDAP or stochastic acceleration.
We found significant correlation between flux and α, which is an indication of hardening of the spectra as flux increases, which is very common among HSP blazars. It indicates that hard X-rays vary more rapidly than soft X-rays (Giebels et al. 2002;Holder et al. 2003;Tagliaferri et al. 2003;Cui 2004) or  that there is an injection of fresh electrons with an energy distribution harder than that of the previously cooled electrons (Mastichiadis & Kirk 2002).

Discussion
We observed the blazar 1ES 1959+650 during the period 2018 June-2020 December when the source showed different flux states (high/low) to study its flux and spectral variability on intraday and long-term timescales using the Swift satellite. The source was studied for intraday flux variability on a total of 125 nights and it exhibited significant variability on only five nights. We did not find any spectral variability during these observations. The source was studied using two observations of XMM-Newton held on 2019 July 5 and 2020 July 16 and it exhibited significant flux variability in both of these observations with low-amplitude variation of 1.95% and 3.12%, respectively. The source favors the shock-in-jet model, where intraday variation is triggered by the interaction of the shock front with jet inhomogeneities, by turbulence behind a shock front, and by smallest-size jet turbulent structures (Marscher & Gear 1985;Wagner & Witzel 1995;Sokolov et al. 2004). Smallest-size jet structures are considered to produce rapidly variable emission according to light travel arguments. We found a flux doubling timescale of 15.27 ks (between MJD 58,670.07 and MJD 58,670.13), which leads to a size of the emitting region of 6.56 × 10 15 cm and a black hole mass estimate of 2.95 × 10 8 M e , which was obtained as follows: where G is the gravitational constant (e.g., Zhang et al. 2021). Yuan et al. (2015) reported optical timescales in the range of 23 minutes to 3.72 hr with a Kerr black hole mass in the range Two XMM-Newton observations of 1ES 1959+650 suggest the presence of soft as well as hard lags, which implies that the t acc and t cool of the emitting region change from epoch to epoch. This is consistent with previous studies. Hard lags suggest that the flux evolution is dominated by acceleration processes while soft lags suggest that the emission region is dominated by cooling mechanisms. The first-order Fermi acceleration process (e.g., Kirk et al. 1998, and references therein) and statistical/ stochastic processes involving second-order Fermi acceleration (Becker et al. 2006;Katarzyński et al. 2006) are the most acceptable models for particle acceleration and therefore, hard delays are modulated by variations of the acceleration parameter ξ. We estimated the value of the magnetic field to be 0.64 ± 0.05 G using the soft lags and the values were found to be close to those reported in literature for our source.
PSD analysis was done using two observations of XMM-Newton, which were used to characterize the variability on intraday timescales. Accretion disk models typically produce PSD slopes between −1.30 and 2.10 (Zhang & Bao 1991;Mangalam & Wiita 1993;Kelly et al. 2011) while jet-based models yield steeper slopes between −1.70 and 2.90 (Calafut & Wiita 2015;Pollack et al. 2016;Wehrle et al. 2019). The observed PSD slopes of −2.41 and −2.15 for intraday light curves are steeper compared to those predicted by accretion disk models and are more consistent with jet-based models. However, the small number of PSDs used here provides a tentative hint favoring fluctuations originating in jets.
On long-term timescales, the source exhibits a high state during Period 1 with two flares with flux reaching up to 34.82 counts s −1 . During Period 2, the flux of the source decreases to 10.35 counts s −1 and during Period 3, the flux of the source is lowest, reaching up to 7.43 counts s −1 . The source shows significant spectral variability throughout our observations. We studied the flux distributions of our source during different observational periods and found that the source follows a lognormal behavior during the total epoch. Kapanadze et al. (2020a) studied the lognormality of Mrk 421 using multiwavelength data and found that lognormal fits were preferable to normal fits for most of their data set. The flux variability of the source is attributed to the propagation of shocks downstream the relativistic jets (Sokolov et al. 2004). The formation of these shocks might be related to turbulence/inhomogeneities occurring in the accretion disk (Kushwaha et al. 2016;Sinha et al. 2017;Kapanadze et al. 2020b;Kushwaha & Pal 2020). However, this is not always the case. Many observations including flaring events show deviations from a lognormal distribution, which indicate that these flares might be triggered by the interaction of shocks with jet inhomogeneities, which could be related to jet instabilities (Marscher 2014). All the spectra are well described by a log-parabolic model yielding spectral curvature ranging between 0.23 and 0.99 and a photon index that varies between 1.34 and 2.25. During our observations, the peak energy E p varies in such a way that it shifts to higher energies as flux increases. We studied the correlation between the spectral parameters derived from the log-parabolic model. We found a weak negative correlation between α and β, which is due to the coexistence of stochastic and statistical acceleration processes (Kapanadze 2018;Kapanadze et al. 2018b). The spectral hysteresis analysis of 1ES 1959+650 shows an interplay between the acceleration and cooling timescales of the emitting particles and the flux variability timescale (Kapanadze et al. 2018a(Kapanadze et al. , 2018b. The correlation between E p and L p also follows a power law as described in Section 4.6. During Period 1, we found spectral changes are caused by the variations of the average electron energy. The anticorrelation between E p and β is expected for the efficient stochastic acceleration of electrons by magnetic turbulence, which is not seen in our observations. The weak correlation between these parameters implies the coexistence of stochastic and statistical acceleration processes in the emitting region.

Conclusion
The main findings of this work are summarized as follows: 1. Swift-XRT and XMM-Newton EPIC-pn observations were used to study the HSP 1ES 1959+650 during the period 2018 June-2020 December on a total of 127 nights. Significant variability was detected on a total of seven nights with the flux variability amplitude varying between 1.95% and 3.12%. HR analysis shows no significant spectral variability in any of the nights. The flux doubling timescale was found to be 15.27 ks and the black hole mass was calculated to be 2.95 × 10 8 M e . 2. Using XMM-Newton observations, cross-correlation between the soft band (0.3-2 keV) and hard band (2-10 keV) was performed using the DCF method. Both DCF plots are correlated and hard and soft lags of 360 and 940 s, respectively, were found, providing a magnetic field strength of ∼0.64 ± 0.05 G in the jet. 3. PSD analysis was performed using XMM-Newton observations and the power-law slopes were found to be −2.41 and −2.15, which favor the jet-based model (as discussed in Section 5).
4. Intraday light curves were checked for lognormal behavior and it was found that they are well modeled by normal as well as lognormal distributions. As suggested by Gaidos et al. (1996) and Narayan & Piran (2012), variations on timescales of minutes/hours are independent of accretion disk fluctuations and could be attributed to some linear/nonlinear perturbations in the physical parameters used to model the relativistic jets in blazars. The most plausible model to explain short-term variability is turbulence in the jet behind the reconfinement shock that contains multiple synchrotron-emitting cells of different sizes within a single emitting region, which yields remarkable light curves, PSDs, and polarization variations (Marscher 2014;Pollack et al. 2016). 5. The source exhibits lognormal behavior on long-term timescales. Long-term variations are explained by superposition of small flares on long-term trends. As blazar jets are highly magnetized, variability may be incorporated by the minijets-in-a-jet model as we found a linear correlation between the rms and flux. However, as the photon index distribution is well fitted by Gaussian as well as normal distributions, we cannot rule out the possibility of propagation and evolution of relativistic shocks through the jet leading to variability in X-ray bands. These shocks could be related to an abrupt increase of the plasma injection rate at the jet base owing to instabilities in the accretion disk. 6. On long timescales, the source shows high-flux as well as low-flux states. A log-parabolic model is required to describe the X-ray spectra of this source yielding spectral curvature values that range between 0.23 and 0.99 and a photon index that ranges between 1.34 and 2.25. The position of the synchrotron SED peak E p varies between 0.46 and 2.88 keV. The synchrotron peak E p strongly correlates with flux, which implies that E p increases as flux increases. HR analysis on long-term timescales indicates that the source follows the "harder-whenbrighter" trend. 7. The source shows a weak correlation between the photon index α and curvature β, which could be due to the combined effect of statistical/stochastic processes. Electrons can be accelerated at the shock front via EDAP but they can gain energy via the stochastic mechanism after escaping the shock front (Katarzyński et al. 2006). 8. The blazar 1ES 1959+650 shows a low spectral curvature (β ∼ 0.23-0.99) and an anticorrelation between E p and 1/β was found, which might be due to the coexistence of stochastic and statistical acceleration processes.