Temporal and spectral study of PKS 0208-512 during the 2019-2020 flare

Blazars are the class of active galactic nuclei (AGN) with relativistic jets pointing towards the observer within a few degrees (Urry & Padovani 1995). The jet emission is predominantly non-thermal in nature covering the entire electromagnetic (EM) spectrum ranging from radio to very high energy $\gamma$-ray. Their broadband spectral energy distribution (SED) is characterized by two prominent peaks, with the low energy peak generally observed between infrared (IR) to soft X-ray, while the high energy one falls at Mev/GeV energies. The low energy component is well explained by the synchrotron emission from relativistic electrons spiraling along the magnetic fields in the jet. The high energy component is generally modelled as the inverse-Compton (IC) scattering of low energy photons by the relativistic electrons. The source of low energy target photons is either internal (synchrotron emission) or external (accretion disk, broad-line region, dusty torus, etc.) to the jet (Boettcher et al. 1997; Ghisellini & Madau 1996; Błażejowski et al. 2000). Such a model is commonly referred as a leptonic model (Ghisellini & Tavecchio 2009; Ghisellini et al. 2014). If the target photons for the inverse-Compton scattering is synchrotron photons itself then the process is referred to as synchrotron self-Compton (SSC; Sikora et al. 2009). On the other hand, the scattering of external photons is known as external-Compton (EC; Dermer et al. 1992; Sikora et al. 1994; Ghisellini et al. 1998) process.

One of the intriguing properties of blazars is that they show strong stochastic variability and occasional strong flares across the EM spectrum. The flux doubling time scale varies from minutes to several days (Paliya et al. 2015, Shukla & Mannheim 2020, Chatterjee et al. 2021). The rapid flux variability of the order of minutes/hours suggests the emission is produced from a region very close to the central engine (Ghisellini & Tavecchio 2008; Narayan & Piran 2012). In the past three decades, development in ${\gamma}$-ray astronomy has been an important tool to understand the spectral and temporal behavior of blazars (von Montigny et al. 1995; Abdo et al. 2010a; Abdo et al. 2010b; Madejski & Sikora 2016). Particularly, the launch of the Fermi-Large Area Telescope (LAT, Atwood et al. 2009) opened up a new window to study high energy $\gamma$-ray blazars. The recently updated Fermi-LAT catalog, 4FGL (Abdollahi et al. 2020) has more than five thousand $\gamma$-ray sources, 95% of which are blazars. A detailed study of the $\gamma$-ray blazars will help us to understand the acceleration and the dynamics of the underlying particle distribution. However, this also demands additional simultaneous broadband observation of the source in tandem with $\gamma$-rays. In the recent past, various telescopes i.e., the Neil Gehrels Swift observatory (Gehrels et al. 2004), Steward observatory (Smith et al. 2009), Owens Valley Radio Observatory (OVRO) 40-m (Richards et al. 2011) etc. operating at different energy ranges have been used to observe blazars along with Fermi-LAT. The simultaneous or quasi-simultaneous observations available from these instruments were collectively used to study the complex nature of blazars (Prince 2019; Prince et al. 2021; Shah et al. 2021).

The flat spectrum radio quasar (FSRQ) PKS 0208-512 is located at redshift z = 1.003 (Healey et al. 2008). The Parkes radio survey discovered the source (Bolton et al. 1964) and was observed regularly in ${\gamma}$-rays (Thompson et al. 1995) by EGRET onboard Compton Gamma-Ray Observatory (CGRO). The CGRO mission detected excess emission from PKS 0208-512 at soft $\gamma$-rays (1-3 MeV) during 1993 May-June. The $\gamma$-ray SED indicated a broad peak at MeV energies (Blom et al. 1995), later confirmed by INTEGRAL/SPI detection (Zhang et al. 2010). Henceforth, the source has been identified as ‘MeV blazar’.

PKS 0208-512 has been monitored regularly by the Fermi-LAT in a wide range 0.1-300 GeV and during October 2008 a moderate $\gamma$-ray flare was detected (Tosti 2008). The contemporaneous observations in optical/IR during this $\gamma$-ray flare, has been observed by SMARTS and ANDICAM (Buxton et al. 2008). Subsequently, several studies on the multiwavelength property of PKS 0208-512 have also been undertaken (Ghisellini et al. 2010; Zhang et al. 2010; Szostek 2011; Nesci et al. 2011; Blanchard et al. 2012; Chatterjee et al. 2013a; Chatterjee et al. 2013b).

From November, 2019 to May, 2020, PKS 0208-512 went through a major $\gamma$-ray flare and the follow-up ATel initiated a major coordinated observation of the source (Angioni 2019; Lucarelli et al. 2019; Angioni 2020; Angioni et al. 2020). The preliminary analysis of Fermi-LAT observations suggested that the daily averaged $\gamma$-ray flux at energy >100MeV enhanced up to 1.1$\pm$0.2 $\times 10^{-6}$ ph cm^-2 s^-1 on November 29, 2019, and 2.0$\pm$0.3 $\times 10^{-6}$ ph cm^-2 s^-1 on March 15, 2020 (Angioni, 2019, 2020). Further, AGILE observed the source during December 14-16, 2019 and using multi-source maximum likelihood analysis it measured a flux of 2.7$\pm$0.8 $\times 10^{-6}$ ph cm^-2 s^-1 for energy >100 MeV (Lucarelli et al. 2019).

The follow-up observations in the optical/UV and X-ray bands have been performed by Swift telescope under the Target of opportunity (ToO) program during December 17-24, 2019, and March 9-21, 2020. The rich simultaneous multi-wavelength observation of PKS 0208-512 during different flaring activity encouraged us to perform a detailed broadband timing and spectral study of the source. Our primary study is based on the $\gamma$-ray observation by Fermi-LAT and we supplement this with the simultaneous observation by Swift-XRT/UVOT in X-ray and optical/UV wavebands. The paper is organised as follows: In §2, we describe the multi-wavelength observations and data analysis procedure. In §3.1, we provide the broadband lightcurve and the identification of different flux states followed by the results obtained from the temporal and spectral analysis of $\gamma$-ray observations (§3.2 – §3.5). The results from Swift observations are given in §3.6. The detailed broadband spectral modeling of the source using synchrotron and inverse Compton emission processes is described in §4. Finally, the study results are discussed and summarized in §5. Throughout this work, we adopt a cosmology with $\Omega_{m}=0.3$, $\Omega_{\Lambda}=0.7$ and $H_{0}=70$ km s^-1 Mpc^-1.

Fermi-LAT observation of PKS 0208-512 has shown enhanced flaring activity in $\gamma$-ray during the period 2019-2020 (Angioni, 2019, 2020). A daily averaged $\gamma$-ray flux (E>100MeV) of (1.1$\pm$0.2)$\times 10^{-6}$ ph cm^-2 s^-1 was reported for this source on November 29, 2019 (Angioni 2019). Subsequently, an increased flux with 2.0$\pm$0.3 $\times 10^{-6}$ ph cm^-2 s^-1 was observed on March 15, 2020, which is the highest daily binned $\gamma$-ray flux ever detected from this source (Angioni 2020). In the 2-days binned $\gamma$-ray lightcurve extending from September 4, 2017 to July 20, 2020 (MJD 58000-59050) and integrated over the energy range 100 MeV - 300 GeV, the peak flux (1.65$\pm$0.15) $\times 10^{-6}$ ph cm^-2 s^-1 is observed on March 15, 2020. This flux value is slightly lower than the daily binned flux reported by Angioni 2020. Moreover, we noted that 2-days binned $\gamma$-ray flux detected by Fermi-LAT is $\sim$ (1.43$\pm$0.14)$\times 10^{-6}$ ph cm^-2 s^-1 on December 16, 2019, which is lower than the 2-days integrated flux $\sim$ (2.7$\pm$0.8)$\times 10^{-6}$ ph cm^-2 s^-1 detected by AGILE on the same day. Such discrepancy in flux value can be possibly associated with the detection from different instruments and the analysis procedure. The maximum fluxes in X-ray, optical (B-band) and UV (W1-band), overlapping with the $\gamma$-ray lightcurve, are found to be (1.8$\pm$0.2) $\times 10^{-11}$ erg cm^-2 s^-1, (2.14$\pm$0.06) $\times 10^{-11}$ erg cm^-2 s^-1 and (1.65$\pm$0.05) $\times 10^{-11}$ erg cm^-2 s^-1 falling on March 9, 2020 as indicated by Swift-XRT/UVOT analysis. Here, we perform a detailed temporal and spectral study of the brightest $\gamma$-ray flare witnessed during MJD 58780-59000.

To perform the broadband SED modelling of the source PKS 0208-512 during different flux states, we selected three time segments including one quiescent state (Q state) MJD 58660 - 58676, and two flaring states with MJD 58813 - 58840 (F1_b state) and MJD 58888 - 58935 (F2_b state) from the lightcurve as shown in Figure 1. The selection was made such that each segment have simultaneous $\gamma$-ray, X-ray, UV and optical observations. The details of the Swift observations for states Q, F1_b, and F2_b are given in Table 4.

We carried the broadband SED modelling of the three selected states, using the one-zone leptonic model (Sahayanathan et al. 2018). In this model, we assume a scenario where the non-thermal emission from the blazar jet originates from a spherical blob of radius ‘R’ moving with a bulk Lorentz factor $\Gamma$ down the jet at an angle $\theta$ with respect to the observer’s direction. The emission region is permeated with the tangled magnetic field B, and populated by a broken power-law electron distribution described as

Here, K and $\gamma$ are the normalization and the electron Lorentz factor, while the low and high energy power-law indices are denoted by p1 $\&$ p2. The $\gamma_{min}$, $\gamma_{max}$ are the minimum, and maximum Lorentz factors of the electron, and $\gamma_{b}$ is the Lorentz factor corresponding to the break energy. The observed broadband emission from the source is modelled as synchrotron and inverse Compton (IC) emission by this electron distribution. The seed photons for the IC process are the synchrotron photons themselves (SSC) and the photon field external to the jet (EC). The source of external target photons considered for the SED modelling are the IR photons from the dusty torus (EC/IR) and the dominant Lyman-alpha line emission from the broad line region (EC/BLR). The emission from the dusty torus is assumed to be a black body with temperature T $\approx$ 1000K (Schartmann et al. 2005). For numerical simplicity, the line emission from BLR is modelled as a black body at T$\approx$ 42000 K (Peterson 2006). Since the final emissivity is estimated by convolving with a relatively broad particle distribution, this simplification will have a negligible effect on the resultant IC spectrum. The numerical model estimating the synchrotron and IC emissivity and the observed flux from the source parameters is added as a local model in XSPEC (Arnaud, 1996), to perform a statistical fit of the SED. The ASCII format of the optical, UV, X-ray, and $\gamma$-ray fluxes are converted to pha format using ftflx2xsp command and the fitting is performed using $\chi^{2}$-statistics implimented in XSPEC.

The SED is mainly governed by the parameters: K, p1, p2, $\gamma_{b}$, B, R, $\Gamma$, $\theta$, temperature T and the fraction of the external target photons up scattered by the IC process. Clearly, the number of parameters exceeds the information available at four wavebands. Hence, to constrain the parameters, we enforce additional equipartition between particle and the magnetic field energy densities. The fitting is performed by first running the model on the observed SED to obtain acceptable residuals. Then the confidence intervals on parameters are obtained by freezing certain parameters to their best fit value. This is to avoid degeneracy introduced by the limited available observed information. Besides these, a 15 % error was evenly added to the data to get the reduced $\chi^{2}_{\rm red}<2$. This limit is demanded by XSPEC to obtain the confidence intervals on the fit parameters. The best fit model and the parameters corresponding to Q, F1_b, and F2_b are given in Figure 8 and Table 5. The parameters which are kept constant for all the epochs are the size of the emission region R = 2.44$\times 10^{16}$ cm (consistent with the observed variability time), viewing angle $\theta\approx$ 2^∘, the minimum and maximum values of particle Lorentz factor ($\gamma_{min}$ and $\gamma_{max}$) = 50 and $9.31\times 10^{5}$, and the energy density of the up scattered IR photon field as $\approx 4\times 10^{-12}$ erg/cm³. On the other hand, the parameters p1, p2, $\gamma_{b}$, $\Gamma$, and B are let to vary freely.

The $\gamma$-ray spectrum corresponding to the quiescent state can be readily explained by EC/IR mechanism alone; however, for the flaring states a combination of EC/IR and EC/BLR is required to fit the SED successfully. For the latter, the energy density of the up scattered BLR photons is set at $\approx 3\times 10^{-16}$ erg/cm³. If the emission region is closer to the central black hole within BLR then the $\gamma$-ray spectrum will have the contribution from both EC/BLR and EC/IR. On the other hand, if the emission region is away from BLR but within the IR torus then EC/IR will be dominant. Therefore, our broadband spectral modelling suggests during the activity state the emission region is close to the central black hole than the quiescent state. This inference; however, may banks upon the assumptions (e.g. equipartition) imposed on the theoretical SED model and the choice of fixed parameters.

The long-term monitoring of PKS 0208-512 by Fermi-LAT, together with the simultaneous observations in other wavebands by Swift, allow us to study the temporal and spectral behavior of the source in detail. The source showed an intense flaring activity during MJD 58780-59000. During this period, the source revealed two major outbursts, viz. F1 and F2. We obtained the fastest variability time-scale during F1 and F2 states as $\sim$ 2 days. However, using a 2-days binned lightcurve, the obtained variability time-scale of the order of 2-days may not be accurate. This may have an impact on the estimation of size and location of the emission region. Therefore, we generated the 1-day binned $\gamma$-ray lightcurve and estimated the fastest variability time of 0.95$\pm$0.35 day in F1-state, the corresponding values of R and d are obtained as $\sim$ 1.2$\times 10^{16}$ cm and $\sim$ 2.5$\times 10^{17}$ cm, respectively. These values do not deviate much from the values obtained from 2-days binned lightcurve and differ only by a factor of $\sim$2. Nevertheless, the larger error on the flux points in the 1-day binned lightcurve hinders us in locating the peaks of flares. Therefore, we considered the 2-days binned lightcurve in order to obtain the temporal characteristics of the source. Using 2-days variability time, the size and location of the $\gamma$-ray emission region are estimated as $\lesssim$ 2.6 $\times 10^{16}$ cm and $\sim$ 5.2$\times 10^{17}$ cm, respectively. Considering a black hole of mass $10^{9}$ $M_{\odot}$ (Ghisellini et al. 2011), the derived emission region size is $\sim$ 2 order of magnitude larger than the Schwarzschild radius ($\sim$ 3.0$\times 10^{14}$ cm). Relatively a longer time variability has been reported by using EGRET observations. For instance,von Montigny et al. 1995 found the variability time-scale of 8 days and this corresponds to the size of $\gamma$-ray emission as $\sim$ 3.2$\times 10^{16}$ cm for $\delta$ = 3.1. The authors considered a black hole of mass 10¹⁰ $M_{\odot}$ and the derived emission region size is $\sim$ 16 times larger than the Schwarzschild radii (1.9$\times 10^{15}$ cm). A variability time-scale of 4-days was also reported and this relates to the emission region of size $\sim$ 5$\times 10^{15}$ cm for $\delta\sim 1$ (Stacy et al. 2008). However, the constraints obtained from the pair production opacity of the $\gamma$-rays suggest a lower limit on the Doppler factor as $\delta\gtrsim 3.2$ (Stacy et al. (2008)). A variability time of 4-days would then imply the distance of the emission region as, d $\sim$ 1.06$\times 10^{17}$ cm. For a disk luminosity L_d = 18 erg/sec (Ghisellini et al. 2011), the size of the BLR can be estimated as R_BLR $\sim$ 4.24$\times 10^{17}$ cm. These results imply that the $\gamma$-rays are produced from a region close to the black hole. The multi-wavelength variability study between the fluxes in different wavebands suggests the variability is large in case of $\gamma$-rays and optical/UV compared to that of X-rays. This can be interpreted in terms of the electron energy responsible for the emission at these wavebands. The optical/UV emission corresponds to synchrotron emission from high energy electrons, whereas X-rays can be due to the IC emission from low energy electrons. Similarly, the $\gamma$-ray emission can be due to IC emission from high-energy electrons. This is also evident from the broadband SED modelling of the source.

The multiple peaks observed in the 2-days binned $\gamma$-ray flux lightcurve during states F1 and F2 are fitted with sum of exponentials. The result suggests that most of the peaks have characteristic rise and decay times of the order of 1-day. This finer sub-structures seen in $\gamma$-ray lightcurves of F1 and F2 are further divided as F1_a, F1_b, F1_c, F2_a, F2_b, F2_c, F2_d. The $\gamma$-ray spectral study of these sub-flares modelled by the PL, LP, and BPL model shows that the states F1_b, F2_b, F2_c, F2_d revealed a significant curvature with TS_curve values as 16.05, 31.23, 20.40, and 19.83, respectively. The observed curved spectrum is best fitted by the log-parabola (LP) and/or the broken power-law (BPL) functions. This also suggests significant curvature in the underlying electron distribution. Such curvatures in the electron distribution can be attributed to the energy dependence of the particle acceleration/escape time scales (Massaro et al. 2004; Jagan et al. 2018; Goswami et al. 2018). In Massaro et al. 2004, the authors showed that a curved particle distribution could originate when the acceleration probability is energy-dependent. The energy dependence of the escape timescale in the acceleration region can also be responsible for the curvature observed in the spectrum (Jagan et al. 2018; Goswami et al. 2018). Alternatively, EBL induced pair production absorption can also be a cause for curved $\gamma$-ray spectrum. However, we observed photons below 20 GeV in all the states and the EBL induced absorption optical depth for the redshift of PKS 0208-512 (z=1.003) at this $\gamma$-ray energies are very minimal (Franceschini et al. 2008). Hence it cannot account for the curvature observed in the $\gamma$-ray spectrum.

The broadband SED of PKS 0208-512 during different flux states have been fitted by a radiative model considering synchrotron, SSC, and EC processes. We find that the one-zone leptonic model can successfully explain the SED for all these three states. From the best fit parameters, it appears that no significant variation in the bulk Lorentz factor $\Gamma$ happened during different flux states, while a marginal decrease is noticed during the flare states. However, the observed information is not sufficient enough to constrain all the parameters and hence this conclusion cannot be treated as a robust one. Further, the $\Gamma$ values obtained are less compared to the typical value as observed in other FSRQs (Shah et al. 2017, 2019). Nevertheless, our broadband SED fitting suggests that the $\gamma$-ray spectra during flaring (F1_b $\&$ F2_b) states demand the EC scattering of photons from both BLR and dusty torus. On the other hand, for the quiescent (Q) state, EC scattering of IR photons alone is capable to explain the $\gamma$-ray spectrum. This result can be interpreted in terms of the location of the emission region during different activity states. For instance, during F1_b and F2_b states, the emission region is well within BLR and hence the photons from BLR as well as IR torus will be upscattered through EC process. While in the case of Q state, the emission region may lie ahead of the BLR. From Table 5, it appears $\Gamma$ decreases marginally during the flare states. This may appear to contradict to our inference since the emission region is closer to the central black hole during active states and it is expected the $\Gamma$ to be larger. However, within the estimated confidence intervals we do not find any significant change in $\Gamma$. Additionally, any deviation from the equipartition during the flaring state can also modify $\Gamma$ significantly. The flare state is also associated with a marginal increase in magnetic field and the particle energy density (through equipartition). This maybe associated with efficient particle acceleration during the flaring states as compared to the quiescent state. This is also consistent with the change in power-law index (p1) of the particle distribution where it hardens during the flaring states. Non-availability of sufficiently high energy spectra (which can probe the high energy tail electron distribution) do not let us to obtain the upper bound of p2. However, its quite evident the difference between p1 and p2 is larger than 1/2 and hence the origin of the broken power-law may not be strictly associated with the radiative cooling interpretation of a power-law electron distribution (Kardashev 1962).

For MeV blazars, X-rays can be due to the EC mechanism as well (in the framework of leptonic model, see e.g. Ghisellini & Tavecchio 2009; Sbarrato et al. 2015; Ajello et al. 2016; Marcotulli et al. 2020). However, Chatterjee et al. 2013b has performed the broadband study of PKS 0208-512 using one-zone leptonic model, where the optical, X-ray, and $\gamma$-ray emissions can be explained by the synchrotron, SSC, and EC processes. Further, explanation of X-ray emission of FSRQs through EC process demand a significant deviation from the equipartition (Sahayanathan & Godambe 2012; Shah et al. 2017). In the present work also, we find that the SSC and EC interpretation of the X-ray and $\gamma$-ray emission from PKS 0208-512 do not require any deviation from the equipartition condition. The optical/UV/X-ray observation period considered for the broadband SED modelling does not encompass the selected $\gamma$-ray period fully. However, we believe the average flux obtained during this period do not change the parameters obtained from the broadband fitting significantly unless extreme low/high optical/UV/X-ray fluxes are missed in the time gap. Alternatively, reducing the $\gamma$-ray observation period consistent with the optical/UV/X-ray observation will not deviate the average $\gamma$-ray flux significantly. On the other hand, this may impact in poor statistics due to less number of $\gamma$-ray photons. Hence, the physical parameters obtained through this study largely define the source within the ambit of the model assumptions. Future observations at very high energy by upcoming experiment Cherenkov Telescope Array (CTA Consortium & Ong 2019) have a better potential to relax some of the underlying assumptions and gain better constraint on the physical parameters.

We thank the anonymous referee for valuable comments and suggestions. We acknowledge the use of $\gamma$-ray data, software obtained from Fermi Science Support Center (FSSC), and the Swift-XRT/UVOT data from NASA’s High Energy Astrophysics Science Archive Research Center (HEASARC), service of Goddard Space Flight Center. This work has made use of the XRT Data Analysis Software (XRTDAS) developed under the responsibility of the ASI Science Data Center (ASDC), Italy. R.K. and R.G. would like to thank CSIR, New Delhi (03(1412)/17/EMR-II) for financial support. R.P. acknowledges the support by the Polish Funding Agency National Science Centre, project 2017/26/A/ST9/00756 (MAESTRO 9), and MNiSW grant DIR/WK/2018/12. Z.S. is supported by the Department of Science and Technology, Government of India, under the INSPIRE Faculty grant (DST/INSPIRE/04/2020/002319).

The data and software used in this research are available at Fermi Science Support Center (FSSC), NASA’s HEASARC webpages with the links given in the manuscript.