# X-ray and GeV-$\gamma$-ray emission property of TeV Compact Symmetric Object
PKS 1413+135 and Implication for Episodic Jet Activity

Ying-Ying Gan School of Physics, Beijing Institute of Technology, Beijing
100081, People’s Republic of China<EMAIL_ADDRESS>Jin Zhang† School of
Physics, Beijing Institute of Technology, Beijing 100081, People’s Republic of
China<EMAIL_ADDRESS>Su Yao Max-Planck-Institute für Radioastronomie,
Auf dem Hügel 69, 53121 Bonn, Germany  Hai-Ming Zhang School of Astronomy and
Space Science, Nanjing University, Nanjing 210023, People’s Republic of China
Yun-Feng Liang Guangxi Key Laboratory for Relativistic Astrophysics, School of
Physical Science and Technology, Guangxi University, Nanning 530004, People’s
Republic of China  En-Wei Liang Guangxi Key Laboratory for Relativistic
Astrophysics, School of Physical Science and Technology, Guangxi University,
Nanning 530004, People’s Republic of China

###### Abstract

PKS 1413+135, a compact symmetric object (CSO) with a two-side pc-scale
structure in its miniature radio morphology, is spatially associated with 4FGL
J1416.1+1320 and recently detected with MAGIC telescopes. We comprehensively
analyze its X-ray and GeV $\gamma$-ray observation data for revealing its high
energy radiation physics. It is found that the source is in a low-flux stage
before MJD 58500 and experiences violent outbursts after MJD 58500 in the GeV
band. The flux at 10 GeV varies $\sim 3$ orders of magnitude, and the GeV-flux
variation is accompanied by the clear spectral variation, which is
characterized as a soft log-parabola spectrum in the low-flux state and a hard
power-law spectrum in the bright flares. The variability amplitude of X-rays
is lower than that of $\gamma$-rays, and no correlation of variability between
$\gamma$-rays and X-rays is observed. Fitting the broadband spectral energy
distribution during a GeV outburst with a multi-zone leptonic model, we show
that the GeV $\gamma$-rays are attributed to the external Compton process
while the X-rays are hybrid of several components. The predicted TeV
$\gamma$-ray flux during the GeV outburst is consistent with the detection of
MAGIC telescopes. These results, together with its CSO radio morphology, imply
that PKS 1413+135 has the episodic nuclear jet activities. The weak
$\gamma$-ray emission before MJD 58500 may be from its sub-pc-/pc-scale jet
component powered by previous activities, and the violent outbursts with
short-timescale variability after MJD 58500 could be attributed to the
recently restarted jet activity.

galaxies: active—galaxies: jets—radio continuum: galaxies—gamma rays: galaxies

##  1 Introduction

Compact symmetric objects (CSOs), a sub-class of active galactic nuclei
(AGNs), are defined as those with symmetric twin radio jet structure on both
sides of their nuclei (e.g., Readhead et al., 1996). They are thought to be a
class of misaligned AGNs (Wilkinson et al., 1994; Readhead et al., 1996).
Generally, low polarization, low radio variability, low core luminosity,
hosted in elliptical galaxies, and lack of the optically thick component at
short wavelengths are presented in CSOs (Wilkinson et al., 1994; Readhead et
al., 1996). CSOs are characterized by an overall size being less than 1 kpc.
The small size of CSOs may result from the frustrated jet due to the dense
interstellar medium (O’Dea et al., 1991; Carvalho, 1994, 1998). The ages
estimated by the advance speed, the lobe supply timescale, and the
synchrotron-loss timescale for CSOs are consistent with each other within the
uncertainties, indicating that CSOs are young with typical ages of 3$\times
10^{3}-10^{4}$ yr, almost certainly $\ll 10^{5}$ yr (Readhead et al., 1996).
The non-difference between these compact radio sources and larger radio-loud
objects at the mid- and far-infrared emission also supports that CSOs are
young (Heckman et al., 1994). Thus, CSOs are an important fraction of radio-
loud AGNs for understanding the formation and evolution of powerful jets in
extragalactic radio sources.

Gamma-ray emission is a critical probe to study the AGN jets. As a new
population of $\gamma$-ray sources, six CSOs have been detected with the
_Fermi_ -LAT in the GeV band, i.e., PMN J1603–4904 (Müller et al., 2014,
2015), PKS 1718–649 (Migliori et al., 2016), NGC 3894 (Principe et al., 2020),
TXS 0128+554 (Lister et al., 2020), CTD 135 (Gan et al., 2021), and PKS
1413+135 (Principe et al., 2021). It is debated that whether the $\gamma$-rays
of CSOs are from the pc-scale lobes or from the core jet region.
Interestingly, the $\gamma$-rays in some CSOs may indicate a characteristic of
recently restarted jet activity (Lister et al., 2020; Gan et al., 2021).
Recently, PKS 1413+135 was announced to be detected by the Major Atmospheric
Gamma-ray Imaging Cherenkov (MAGIC) telescopes on January 12, 2022 (MJD 59591,
Blanch et al., 2022). It is the first CSO detected at the very-high-energy
(VHE; $E>100\rm\,GeV$) band. It would be a good sample to study the
$\gamma$-ray emission of young radio sources.

The high-resolution observations with the very large array (VLA), US very long
baseline interferometry (VLBI) Network, and very long baseline array (VLBA)
show no kpc-scale extended structure is presented for PKS 1413+135, and there
is a mini-triple structure with a counter-jet on pc-scale (Perlman et al.,
1994). Using the high-dynamic range VLBA maps at 3.6, 6, 13, and 18 cm,
Perlman et al. (1996) revealed a complex, two-side pc-scale structure for PKS
1413+135 and no evidence of superluminal motion was obtained, and thus they
suggested that PKS 1413+135 likely is a young radio source with an age $\leq
10^{4}$ yr. It is intriguing that the multi-band optical images of PKS
1413+135 show a surface brightness profile that can be well fitted by an
exponential disk, suggesting that PKS 1413+135 is hosted in a spiral galaxy
(McHardy et al., 1991). The Hubble Space Telescope (HST) observation further
reveals that the galaxy contains a previously unresolved dust lane and very
likely is an early-type spiral galaxy viewed edge-on (McHardy et al., 1994).
In addition, an extremely large column density is needed in the analysis of
_Einstein_ X-ray data (Stocke et al., 1992). The enormous extinction, but the
lack of the thermal IR emission and strong narrow emission lines, suggests
that PKS 1413+135 would not be hosted in this spiral galaxy and be a
background source behind the spiral galaxy (Stocke et al., 1992). The redshift
of z=0.247 for PKS 1413+135 is derived from the redshift HI absorption
(Carilli et al., 1992), which should be associated with the spiral galaxy. The
U-shaped symmetric achromatic variability from 15 to 234 GHz of PKS 1413+135
was proposed to be caused by gravitational lensing of the foreground spiral
galaxy (Vedantham et al., 2017). However, no sign of image multiplicity or
distortion due to the gravitational lensing is found, in either optical or
radio bands (Stocke et al., 1992; Perlman et al., 1996; Vedantham et al.,
2017; Readhead et al., 2021). It was thus argued that PKS 1413+135 should be
located at a redshift range of $0.247<z<0.5$ (Vedantham et al., 2017; Readhead
et al., 2021). Due to the strong obscuring and contamination by the
foreground, the high energy photons, especially $\gamma$-rays in the GeV–TeV
bands, would be powerful probes for studying the nature of the source.

For revealing the radiation properties of PKS 1413+135, we analyze its data in
the X-ray and $\gamma$-ray bands observed with the _Chandra_ , _XMM-Newton_ ,
_Swift_ -XRT, and _Fermi_ -LAT over the past $\sim$ 13.5 yr. Description of
the data reduction is presented in Section 2. We investigate its temporal and
spectral variations in both X-ray and $\gamma$-ray bands (Section 3 and
Section 4). Together with the data in the low-energy bands, we compile its
broadband spectral energy distributions (SEDs) and study its $\gamma$-ray
emission via the SED modeling in Section 5. Discussion on the $\gamma$-ray
origin of PKS 1413+135 and its AGN type is given in Section 6. A summary is
presented in Section 7. Throughout, if not otherwise specified, the results
are derived on the basis of $z=0.247$. $H_{0}=71$ km s-1 Mpc-1, $\Omega_{\rm
m}=0.27$, and $\Omega_{\Lambda}=0.73$ are adopted in this paper.

##  2 Observations and Data Reduction

###  2.1 GeV $\gamma$-Ray Observations with Fermi-LAT

It was reported that 4FGL J1416.1+1320 is spatially coincident with PKS
1413+135 in the _Fermi_ -LAT 12-year Source Catalog (4FGL-DR3, Abdollahi et
al., 2022). We download the Pass 8 data from the Fermi Science Support Center
(FSSC). The data are selected within a $15^{\circ}$ region of interest (ROI)
centered on the radio position of PKS 1413+135 (R.A. = $213.995\arcdeg$, Decl.
= $13.340\arcdeg$; Fey et al., 2015). The temporal coverage of the data is
from 2008 August 04 to 2022 March 22 (MJD 54682–59660) of about 13.5 yr. We
perform a binned likelihood analysis for the $\gamma$-rays of 4FGL
J1416.1+1320 using the publicly available software Fermitools (ver. 2.0.8).
Only the $\gamma$-ray photons in the energy range of 0.1–300 GeV and
satisfying the standard data quality selection criteria
“$(DATA\\_QUAL>0)\&\&(LAT\\_CONFIG==1)$” are considered in our analysis. A
zenith angle cut of $90^{\circ}$ is set to avoid the $\gamma$-ray
contamination causing by the Earth limb. We bin the data with a pixel size of
$0.2^{\circ}$ in space and 25 logarithmical energy bins. The background models
include all $\gamma$-ray sources listed in the 4FGL-DR3 Catalog and the
Galactic diffuse component (gll_iem_v07.fits) as well as the isotropic
emission (iso_P8R3_SOURCE_V3_v1.txt). The P8R3_SOURCE_V3 set of instrument
response functions (IRFs) is used.

The spectrum model of 4FGL J1416.1+1320 reported in the 4FGL-DR3 is a log-
parabola function (Abdollahi et al., 2022), i.e.,

$\frac{dN}{dE}=N_{0}(\frac{E}{E_{\rm
b}})^{-(\Gamma_{\gamma}+\beta{\log}(\frac{E}{E_{\rm b}}))},$ (1)

where the decimal logarithm is used, $E_{\rm b}$ is the scale parameter of
photon energy, $\Gamma_{\gamma}$ is the photon spectral index, and $\beta$ is
the curvature parameter (Massaro et al., 2004). So, the model for fitting the
spectrum of 4FGL J1416.1+1320 is in priority selected as a log-parabola
function in our analysis. If $\beta$ is poorly constrained, we set $\beta=0$
and the log-parabola function turns into a single power-law function (i.e.,
$dN/dE=N_{0}(E/E_{\rm b})^{-\Gamma_{\gamma}}$). The spectral parameters of all
sources lying within $8^{\circ}$ are left free, whereas the parameters of
those sources lying beyond $8^{\circ}$ are fixed to their 4FGL-DR3 values.
Also, the normalization parameters of the standard Galactic and isotropic
background templates are set free in the likelihood fit. The significance of
the $\gamma$-ray detection is quantified by adopting the maximum likelihood
test statistic (TS), which is defined as
TS$=2(\ln\mathcal{L}_{1}-\ln\mathcal{L}_{0})$, where $\mathcal{L}_{1}$ and
$\mathcal{L}_{0}$ are maximum likelihood values for the models with and
without the target source, respectively. The analysis results of the _Fermi_
-LAT observations for PKS 1413+135 are given in Table 1.

###  2.2 X-ray observations with Swift-XRT, Chandra, and XMM-Newton

PKS 1413+135 was observed by the _Swift_ -XRT, _Chandra_ , and _XMM-Newton_ in
several epochs, as listed in Table 2. It was observed by the _Chandra_
Advanced CCD Imaging Spectrometer S-array (ACIS-S) as a general observer
target for an exposure of 20 ks on March 20, 2007. It was also observed
several times by _Chandra_ again during 2019 and 2020 as a Director’s
Discretionary Time target for exposures of 4 ks. As a significant foreground
absorption of soft X-rays was found in PKS 1413+135 (Stocke et al., 1992),
firstly we attempt to get the equivalent neutral hydrogen column density of
this absorption using the observation in 2007. The data are reduced using ciao
(version 4.12) and caldb (version 4.9.2.1). The level-2 event file is created
following the standard procedure. We extract source photons from a circle
centred on the radio position of PKS 1413+135 with radius of 6 arcsec. The
background is determined in an annulus region with inner and outer radius of 7
and 14 arcsec, respectively. There are $\sim 600$ net photon counts. The
spectrum is grouped to have at least 25 counts per bin and the $\chi^{2}$
minimization technique is adopted for spectral analysis. The spectrum is
fitted by a single power-law absorbed by two absorption components, one is an
absorption at $z=0$ with neutral hydrogen column density fixed at Galactic
value $N^{\rm H}_{\rm gal}=1.56\times 10^{20}$cm-2 (HI4PI Collaboration et
al., 2016), the other one is an extragalactic foreground absorption $N^{\rm
H}_{\rm int}$ at redshift of $z=0.247$ with column density set free. The model
gives a well fit to the spectrum and an extragalactic foreground absorption
column density of $N^{\rm H}_{\rm int}=3.6^{+0.6}_{-0.5}\times 10^{22}$cm-2,
which is consistent with the results in Stocke et al. (1992) and Perlman et
al. (2002). We generate the _Chandra_ spectra at other epochs following the
procedures mentioned above. The single power-law with the Galactic ($N^{\rm
H}_{\rm gal}$) and extragalactic ($N^{\rm H}_{\rm int}$) absorption is used to
fit the spectra. As there are only no more than 150 net photon counts for each
spectrum, $C$-statistic minimization is adopted for evaluating the goodness of
the fits and $N^{\rm H}_{\rm int}$ is fixed at $3.6\times 10^{22}$cm-2 during
the fitting.

The public data of pointing observations on PKS 1413+135 by the _XMM-Newton_
are reduced with the XMM–Newton science analysis system (version 18) following
standard procedures. We generate the spectra using the data from PN CCD arrays
since PN has a larger effective area. The source events are extracted from a
circle of 32 arcsec radius centred at source position while the background
events are extracted from a circle of the same radius in a source-free region
nearby. The spectrum is binned to contain at least 25 counts per bin required
for $\chi^{2}$ analysis. Again, the single power-law with a Galactic and an
extragalactic absorption is used to fit the spectra. The absorption is also
fixed as mentioned above during the fitting.

There are 21 observations by the Neil Gehrels _Swift_ observatory for PKS
1413+135 from 2007 to 2022. The X-ray telescope (XRT) onboard the _Swift_
satellite was operating in the photon counting mode with exposure times of
0.1–10 ks. We collect the XRT data from the _Swift_ archive and reproduce the
clean events using the xrtpipeline task. The source photons are extracted from
a circle with radius of 50 arcsec, while the backgrounds are determined in an
annulus with inner and outer radius of 60 and 105 arcsec, respectively. There
remain 16 observations after excluding those without significant detection for
this source. As the statistics are too low for each spectrum, $C$-statistic
minimization is adopted to evaluate the goodness of the fits. $N^{\rm H}_{\rm
gal}=1.56\times 10^{20}$cm-2 and $N^{\rm H}_{\rm int}=3.6\times 10^{22}$cm-2
are all fixed, and the photon index is also fixed at $\Gamma_{\rm X}=1.9$ (the
average value in Figure 3(b)) during the fitting.

##  3 Temporal Analysis

We perform a likelihood fit for PKS 1413+135 with the $\sim$13.5-year _Fermi_
-LAT data. PKS 1413+135 is bright in the GeV $\gamma$-ray band, with an
average flux of $F_{0.1-300\rm\,GeV}=(2.28\pm 0.13)\times 10^{-11}$ erg cm-2
s-1 in the whole interval. Variability is very common for AGNs and its
magnitude and timescale are useful to study the property of emission region.
For estimating the variability of PKS 1413+135 at the GeV band, we follow the
definition in 2FGL (Nolan et al., 2012) and derive the variability index TSvar
(see also Abdollahi et al., 2020). We split the full $\sim$13.5-year interval
into N=13 intervals of about one year each. The source is considered to be
significantly variable if TSvar exceeds 53.9, where TS${}_{\rm var}=53.9$
corresponds to the $5\sigma$ confidence level in a $\chi^{2}_{N-1}$(TSvar)
distribution with 12 degrees of freedom. The TSvar value of PKS 1413+135 is
2074 and far beyond $5\sigma$ confidence level, meaning that its $\gamma$-rays
are significantly variable. We extract the $\gamma$-ray light curve of PKS
1413+135 using an adaptive-binning method based on a criterion of TS$\geq 9$
for each time bin, where the minimum bin size is set to be seven days, as
shown in Figure 1(b). We note that before MJD 58500, the $\gamma$-ray emission
is in a relatively low state, where the $\gamma$-ray fluxes are lower than the
average flux level for almost all time bins and even lower than the average
flux by almost an order of magnitude in some time bins. The $\gamma$-ray flux
begins to increase from a very low state since $\sim$ MJD 58200 to the average
flux level around MJD 58500. After MJD 58500 the source seems to be in the
high state with fluxes higher than the average flux, and almost all the time
bins only need the minimum bin size of seven days to meet TS$>$9\. Therefore,
the whole $\gamma$-ray light curve of PKS 1413+135 seems to present two
totally different phases with a boundary of MJD 58500, i.e., a low-flux stage
and a outburst stage.

In order to further investigate the short timescale variability of this
source, we reanalyze the GeV observation data after MJD 58500 and still
extract the light curve using the adaptive-binning method with the criterion
of TS$\geq 9$ for each time bin, where the minimum time bin is taken as one
day. The derived flux of each time bin has a large error bar, as illustrated
in Figure 1(c). However, many serial time bins only need the minimum bin size
of one day to meet TS$>$9, and the variability on a daily timescale for PKS
1413+135 can still be observed. Note that PKS 1413+135 has been announced to
be detected at the VHE band with MAGIC telescopes on January 12, 2022 (MJD
59591, Blanch et al., 2022). This source is also in a relatively high state at
the GeV band during MJD 59589–59591, as shown in Figure 1(c). Its flux is
$F_{0.1-300\rm\,GeV}\sim 2.25\times 10^{-10}$ erg cm-2 s-1 on MJD 59588,
raises to $F_{0.1-300\rm\,GeV}\sim 6.64\times 10^{-10}$ erg cm-2 s-1 on MJD
59590, declines to $F_{0.1-300\rm\,GeV}\sim 2.51\times 10^{-10}$ erg cm-2 s-1
on MJD 59592, and then drops to $F_{0.1-300\rm\,GeV}\sim 3.0\times 10^{-11}$
erg cm-2 s-1 after MJD 59593. We can also observe that there are several time
bins with $F_{0.1-300\rm\,GeV}>10^{-9}$ erg cm-2 s-1 between MJD 58700–59000,
indicating that PKS 1413+135 is in a high-flux state at the GeV band and there
may be detectable VHE $\gamma$-rays during that time.

In the X-ray band, there are 3 _XMM-Newton_ , 6 _Chandra_ , and 16 _Swift_
-XRT observations from March 2007 to January 2022. As shown in Figure 1(a),
most of them cluster at the time after January 2020. Considering the large
errors of data points observed by the _Swift_ -XRT, a weighted mean flux
$\langle{F}\rangle$ and its variance $\sigma_{\langle{F}\rangle}$ are
estimated by (McLaughlin et al., 1996)

$\langle{F}\rangle=\left[\sum\limits_{i=1}^{N}\frac{F_{i}}{\sigma^{2}_{i}}\right]\left[\sum\limits_{i=1}^{N}\frac{1}{\sigma^{2}_{i}}\right]^{-1},$
(2)
$\sigma^{2}_{\langle{F}\rangle}=\left[\sum\limits_{i=1}^{N}\frac{1}{\sigma^{2}_{i}}\right]^{-1},$
(3)

where $N$ is the number of the data points, ${F_{i}}$ and $\sigma_{i}$ are the
flux and its error for the $i$th data point. We obtain
$\langle{F}\rangle_{2-10\rm\,keV}=(6.97\pm 0.17)\times 10^{-13}$ erg cm-2 s-1
for the 25 observational data points, which is shown as the horizontal dashed
line in Figure 1(a). Interestingly, the last two points in Figure 1(a) are
obtained with the observations on January 13–14 2022, just after the ATel of
the VHE detection for PKS 1413+135. We find that the X-ray flux declines from
$1.3\times 10^{-12}$ erg cm-2 s-1 on January 13, 2022 (MJD 59592) to
$5.9\times 10^{-13}$ erg cm-2 s-1 on January 14, 2022 (MJD 59593), a factor of
$\sim 2$ within one day. Perlman et al. (2002) reported that the ASCA
observations yield a flux of $F_{2-10\rm\,keV}=9\times 10^{-13}$ erg cm-2 s-1,
which declined a factor of $\sim$5 compared to the one seen in previous X-ray
observations, but the previous X-ray flux obtained from the _ROSAT_ and
_Einstein_ observations is almost completely based on the extrapolations.
These results indicate that the X-rays of PKS 1413+135 are also variable.

To quantify the variability of X-rays, we calculate
$\chi^{2}=\sum\limits_{i=1}^{N}\frac{(F_{i}-\langle{F}\rangle)^{2}}{\sigma_{i}^{2}}$
and the associated probability $p(\chi^{2})=1-p(>\chi^{2})$ (Chen et al., 2022
and references therein) of the light curve. We find that $\chi^{2}=245.9$ (the
number of the degrees of freedom d.o.f.$=N-1=24$) with $p\ll 10^{-7}$,
indicating that the variability of X-rays is also far beyond $5\sigma$
confidence level. We also calculate the fractional variability amplitude
$F_{\rm var}$ and its uncertainty for the X-ray light curve via (Chen et al.,
2022 and references therein)

$F_{\rm
var}=\frac{\sqrt{\sum\limits_{i=1}^{N}[(F_{i}-\langle{F}\rangle)^{2}-\sigma^{2}_{i}]/N}}{\langle{F}\rangle},$
(4) $\sigma_{\rm
var}=\sqrt{\left(\sqrt{\frac{\sigma^{2}_{\langle{F}\rangle}}{N}}\cdot\frac{1}{\langle{F}\rangle}\right)^{2}+\left(\sqrt{\frac{1}{2N}}\cdot\frac{\sigma^{2}_{\langle{F}\rangle}}{\langle{F}\rangle^{2}F_{\rm
var}}\right)^{2}},$ (5)

where $\langle{F}\rangle$ and $\sigma_{\langle{F}\rangle}$ are the mean flux
and its variance derived with Equations (2) and (3). $F_{\rm var}=0.456$ and
$\sigma_{\rm var}=0.005$ are yielded for the X-ray light curve in Figure 1(a).
Further more, we estimate the fractional variability amplitude with Equations
(2)–(5) for the $\gamma$-ray light curves in Figures 1(b) and 1(c), and obtain
$F_{\rm var}=14.990\pm 0.005$ for the whole $\gamma$-ray light curve and
$F_{\rm var}=7.472\pm 0.003$ for the $\gamma$-ray light curve after MJD 58500,
respectively. Clearly, the variability amplitude of $\gamma$-rays is much
larger than that of X-rays. We also investigate the correlation of variability
between X-rays and $\gamma$-rays using the observation data in Figures 1(a)
and 1(c), as displayed in Figure 2. Considering the large errors of data
points, we use a bootstrap method to estimate the correlation coefficient
($r$) between $F_{2-10\rm\,keV}$ and $F_{0.1-300\rm\,GeV}$, and obtain $r\sim
0.1$. Hence, no correlation of flux variation between two bands is observed.

In the radio band, the significant variations have been widely reported and
studied by the community (e.g., Perlman et al., 1996; Vedantham et al., 2017;
Readhead et al., 2021; Peirson et al., 2022). In the optical-UV band, the
intrinsic optical-UV emission of the source is affected by the foreground
edge-on spiral galaxy and is strongly extinct (Stocke et al., 1992; Perlman et
al., 2002; Readhead et al., 2021). This source is highly variable in radio and
$\gamma$-ray bands, hence its intrinsic optical emission may be potentially
variable.

##  4 Spectral Variation

The time-integrated spectrum of the $\sim$13.5-year _Fermi_ -LAT observations
derived from the likelihood fit is well described with the log-parabola
function, as shown in Figure 3(a). The fitting parameters are reported in
Table 1. As mentioned above, the $\gamma$-ray emission of PKS 1413+135 can be
divided into two stages with a division of MJD 58500. We also extract the
time-integrated spectra in the two stages. It is found that the time-
integrated spectra of the two stages are also well fitted by the log-parabola
function and a smaller curvature parameter value with a harder spectral index
is presented in the high-flux stage than in the low-flux stage. The results
are also given in Figure 3(a) and Table 1. One can observe that the spectrum
in the high-flux stage is harder than that in the low-flux stage.

To further reveal the spectral variation feature, we make time-resolved
spectral analysis for the observational data in the following time slices, MJD
58724–58726, MJD 58796–58798, and MJD 59589–59591, namely, the source being in
the high-flux stage. The time slice of MJD 58724–58726 is within the time bin
of the highest-flux point in Figure 1(b). The VHE observation with MAGIC
telescopes is conducted in the time slice of MJD 59589–59591. In the time
slice of MJD 58796–58798, _Fermi_ -LAT detects the maximum energy photon
($\sim$236 GeV) from the source in its $\sim$13.5-year observations. The
spectral analysis results are given in Figure 3(a) and Table 1. As illustrated
in Figure 3(a), the obvious spectral evolution at the GeV band is presented
for PKS 1413+135. Firstly, the flux variation in the high energy end is more
prominent than that in the low energy end. Above $10^{24}$ Hz, the flux
variation is almost three orders of magnitude. Secondly, the three time-
integrated spectra need a curved log-parabola function to fit with a softer
photon spectral index, while the three time-resolved spectra of $\gamma$-ray
flares can be well modelled by a simple power-law function with a harder
photon spectral index. For the power-law function, we can regard it as a log-
parabola function with the curvature parameter $\beta=0$. As listed in Table 1
and shown in Figure 3(a), the spectral curvature decreases and the spectrum
becomes hard along with the increase of flux, indicating a “harder when
brighter” behavior in the GeV band for PKS 1413+135.

In the X-ray band, we generate the _Chandra_ and _XMM-Newton_ spectra at
different epochs as described in Section 2.2, and obtain the values of
$F_{2-10\rm\,keV}$ and $\Gamma_{\rm x}$ for nine observational epochs. We show
$F_{2-10\rm\,keV}$ against $\Gamma_{\rm x}$ in the
$F_{2-10\rm\,keV}-\Gamma_{\rm x}$ plane, as displayed in Figure 3(b).
Considering the small sample statistics and large errors of data points, we
also use the bootstrap method to estimate the correlation coefficient ($r$) of
the $F_{2-10\rm\,keV}-\Gamma_{\rm x}$ relation and obtain $r=0.35$, and thus
no obvious correlation between $F_{2-10\rm\,keV}$ and $\Gamma_{\rm x}$ is
found for PKS 1413+135.

##  5 $\gamma$-ray Emission Property Derived from Broadband SED Modeling

As mentioned above, the GeV $\gamma$-ray emission has a “harder when brighter”
spectral variation behavior, and the maximum energy of the detected photons is
up to $\sim$236 GeV, as listed in Table 1, where the maximum energy of the
detected photons is estimated with the _gtsrcprob_ tool. We show the
sensitivity curves of MAGIC telescopes (50 hours) and Cherenkov Telescope
Array (CTA) North array (CTA-N, 50 hours) in Figure 3(a), where the
sensitivity curves are obtained from Aleksić et al. (2016) and the CTA
webpage111https://www.cta-observatory.org/science/cta-performance/,
respectively. The high energy end of the time-resolved spectrum during MJD
59589–59591 is over the sensitivity of MAGIC telescopes and the hard spectrum
indicates that the GeV spectrum would extend to the VHE band. Interestingly,
PKS 1413+135 was announced to be detected by the MAGIC telescopes on January
12, 2022 (MJD 59591, Blanch et al. (2022)). The high flux with a very hard
spectrum ($\Gamma_{\gamma}=1.63\pm 0.14$) between MJD 58796–58798, especially
the detection of the maximum energy photon at $\sim$236 GeV, clearly implies
that the flux of PKS 1413+135 at the VHE band should be detectable by the
MAGIC telescopes during that time. As given in Table 1, the maximum energy of
the detected photons during MJD 59589–59591 and MJD 58724–58726 for PKS
1413+135 is 25.7 GeV and 55.7 GeV, respectively. Since only one photon with
energy of $\sim$55 GeV was detected during MJD 58724–58726, an upper-limit is
given for the highest energy bin in the time-resolved spectrum. Although the
average flux for the time slice of MJD 58724–58726 is highest among the six
time slices in Table 1, the slight soft spectrum, especially, the upper-limit
for the highest energy bin, may indicate that the simultaneous VHE emission of
PKS 1413+135 is below the sensitivity of MAGIC telescopes. We observe that
even the time-integrated spectrum of the whole high-flux stage (MJD
58500–59660) is also marginally over the sensitivity of MAGIC telescopes, but
the source may not be detectable by MAGIC telescopes at the VHE band when its
GeV $\gamma$-ray flux is low.

To further study the $\gamma$-ray emission property of PKS 1413+135, we
collect the data at low-energy band from the NASA/IPAC Extragalactic Database
(NED, NASA/IPAC Extragalactic Database, 2019), together with the X-ray and
$\gamma$-ray data derived in this paper, and construct its broadband SEDs in
the low and high states of $\gamma$-ray emission, respectively, as illustrated
in Figure 4. The _Fermi_ -LAT time-resolved spectrum of MJD 59589–59591,
corresponding to the first detection of the source at VHE band, is considered
as the data of high-flux state, and the highest flux point in Figure1(a)
observed by the _Swift_ -XRT on March 26, 2020 (MJD 58934) is taken as the
upper-limit of X-rays. The _Fermi_ -LAT time-integrated spectrum of MJD
54682–58500 and the _Chandra_ observation data on December 20, 2019 (MJD
58837)222It is the data point with the lowest flux in Figure1(a) and also
corresponds to a low-flux state at GeV band. are set as the low-state data. As
shown in Figure 4, it seems that there are two different components at the
radio band, similar feature being observed in some blazars (e.g., Ghisellini
et al., 2015). Using a power-law function of $F_{\nu}\propto\nu^{-\alpha}$ to
fit them, we gain $\alpha=0.66\pm 0.01$ from 80 MHz to 2.7 GHz and
$\alpha=0.09\pm 0.06$ from 4.8 GHz to 375 GHz, respectively; the steep radio
spectrum below a few GHz, $\alpha\sim 0.66$, indicates the signature of the
optically thin synchrotron emission from the large-scale jet extended regions.
the flat radio spectrum from a few GHz up to $\sim 10^{3}$ GHz, $\alpha\sim
0.09$, is produced by the superposition of the emission from several jet
structures located at sub-pc and pc scale. The $\gamma$-ray emission region is
too compact to produce the observed radio radiations below $\sim 10^{3}$ GHz
on account of the synchrotron-self-absorption (SSA) effect (e.g., Ghisellini
et al., 2015).

If the $\gamma$-ray emission region is inside of the broad line region (BLR),
there would be the absorption (Liu & Bai, 2006) and Klein–Nishina (KN,
Ghisellini & Tavecchio, 2009) effects. Especially the detection of the high-
energy photon at 236 GeV by the _Fermi_ -LAT, we thus propose that the
$\gamma$-ray emission region is outside of the BLR. And the torus would
provide the seed photons of the external Compton (EC) process. The synchrotron
radiation energy density can be estimated by $U^{{}^{\prime}}_{\rm
syn}=\frac{L_{\rm syn}}{4\pi R^{2}c\delta^{4}}$ (Tavecchio et al., 1998),
where $L_{\rm syn}=(\frac{1}{1-\alpha_{1}}+\frac{1}{\alpha_{2}-1})L_{\rm s}$,
$L_{\rm s}$ is taken as the luminosity at $\sim 5\times 10^{12}$ Hz (z=0.247),
$\alpha_{1}=0.73$ and $\alpha_{2}=1.34$ (the spectral indices at the GeV band
in Figure 4(a) and (b), respectively), the Doppler boosting factor $\delta\sim
10$ (Peirson et al., 2022). We obtain $U^{{}^{\prime}}_{\rm syn}\sim 4.9\times
10^{-7}$ erg cm-3 for $R=1$ pc and $U^{{}^{\prime}}_{\rm syn}\sim 4.9\times
10^{-3}$ erg cm-3 for $R=0.01$ pc, respectively. The energy density of torus
in the comoving frame is $U^{{}^{\prime}}_{\rm IR}\sim 3\times
10^{-4}\Gamma^{2}$ erg cm-3 (e.g., Ghisellini & Tavecchio, 2009; Kang et al.,
2014), where $\Gamma$ is the bulk Lorenz factor of the emission region. Hence,
the contribution of the EC/torus process comparing with the synchrotron-self-
Compton (SSC) component should be taken into account if the emission region is
located at sub-pc or pc scale from the core.

The electron distribution in radiation region is taken as a power-law or
broken power-law. It is characterized by an electron density parameter
($N_{0}$), a break energy $\gamma_{\rm b}$, and indices ($p_{1}$ and $p_{2}$)
in the range of $\gamma_{\rm e}$ to $[\gamma_{\min},\gamma_{\max}]$, where
$\gamma_{\rm e}$ is the Lorentz factor of electrons. The radiation region is
assumed as a sphere with radius $R$, magnetic field strength $B$, the Doppler
boosting factor $\delta$, where
$\delta=1/(\Gamma-\sqrt{\Gamma^{2}-1}\cos\theta)$, $\Gamma$ and $\theta$ are
the bulk Lorenz factor and viewing angle of the emission region. The
synchrotron (syn), SSC, and EC processes of the relativistic electrons are
considered to fit the broadband SEDs of PKS 1413+135333The same code was used
in other $\gamma$-ray emitting compact radio sources (Zhang et al., 2020; Gan
et al., 2021; Gu et al., 2022).. The spectrum from near-infrared to optical-UV
bands is seriously extinct due to an edge-on intervening Seyfert 2 galaxy
between PKS 1413+135 and the Earth (Stocke et al., 1992; Perlman et al.,
1994), and thus we would not consider these data during the SED modeling. The
KN effect and the absorption of high-energy $\gamma$-ray photons by
extragalactic background light (EBL, Franceschini et al., 2008) are taken into
account during the SED modeling.

The radio spectrum below several GHz in the broadband SED of PKS 1413+135 is
thought to be produced by the emission of jet extended regions. The overall
size of PKS 1413+135 at radio band is $\sim$110 mas (Readhead et al., 2021),
which corresponds to the projection distance of $\sim$422 pc for z=0.247. As
reported by Stawarz et al. (2008), the external photon field is still
dominated by the dusty torus at this scale and its energy density ($U_{\rm
IR}$) can be estimated by Equation (21) in Stawarz et al. (2008). We obtain
$U_{\rm IR}\sim 9\times 10^{-9}$ erg cm-3, which is much higher than the
energy density of cosmic microwave background (CMB) light at z=0.247. Hence
the syn+SSC+EC/torus model under the equipartition condition ($U_{B}=U_{\rm
e}$) is taken to reproduce the emission of the large-scale jet extended
region, where $U_{B}$ and $U_{\rm e}$ are the energy densities of the magnetic
fields and electrons. $R$ is taken as $R=211$ pc, corresponding to the half of
overall size at the radio band. The relativistic effect is not considered,
i.e., $\delta=\Gamma=1$. Considering the limit observation data of the large-
scale jet extended region at radio band, the electron distribution is taken as
a power-law. $p_{1}$ is derived by the radio spectral index ($\alpha\sim
0.66$) below several GHz and is fixed as $p_{1}=2.32$. $\gamma_{\min}=1$ is
also fixed during the modeling. We adjust the values of $\gamma_{\max}$ and
$N_{0}$ to fit the radio spectrum below several GHz in the broadband SED and
obtain $\gamma_{\max}=2000$ and $N_{0}=0.036$ cm-3. And then we fix the
reproduced spectrum component of the large-scale jet extended region in the
following SED modeling for the low-flux and high-flux states.

The GeV emission of the low-flux and high-flux states may have different
origins, and we return to this point in Section 6.1. For the low-flux state of
GeV emission, we propose that the emission is dominated by the sub-pc-/pc-
scale jet structures. The radius of radiation region is taken as $R=0.6$ pc,
corresponding to the half of the projection distance between the radio core
and component-D8 (0.16 mas, Peirson et al., 2022) for z=0.247. Considering the
relativistic effect, $\delta=\Gamma$ is assumed during SED modeling, namely,
the viewing angle ($\theta$) being equal to the opening angle ($1/\Gamma$) of
jet (e.g., Chen & Zhang, 2021). The electron distribution of a broken power-
law is used here. $\gamma_{\min}=1$ and $\gamma_{\max}=10^{5}$ are fixed.
$p_{2}=3.84$ is fixed and constrained by the spectral index in GeV band. We
adjust the values of $B$, $\delta$, $p_{1}$, $\gamma_{\rm b}$, and $N_{0}$ to
fit the SED of the low-flux state in Figure 4(a), and then obtain $B=0.45$ G,
$\delta=2.2$, $p_{1}=1.8$, $\gamma_{\rm b}$=2609, and $N_{0}=2.81$ cm-3. For
the high-flux state at GeV band, we speculate that the outbursts after MJD
58500 in Figure 1(b) are due to the restarted activity of the central engine
and originate from the region more close to the black hole. The radius of
radiation region is taken as $R=\delta c\Delta t/(1+z)$, where $\Delta t=1$
day is taken since the variability on a daily timescale is presented in Figure
1(c). As shown in Figure 4(b), the peak of the second bump is above $10^{25}$
Hz, which needs an electron population with higher energy than that produces
the low-flux emission. We thus consider the syn+SSC+EC processes of another
electron population to reproduce the emission of the GeV outburst on the basis
of the fitting result of the low-flux state. The distribution of the new
electron population is also taken as a broken power-law. $p_{1}=2.46$,
$p_{2}=3.84$, and $\gamma_{\max}=10^{6}$ are fixed, where $p_{1}$ is derived
by the spectral index in GeV outburst while $p_{2}$ is taken the same value as
the low-energy electron population444Although no observational data are
available to constrain the $p_{2}$ value, a broken power-law (not a power-law)
distribution of electrons is generally needed to explain the broadband SEDs of
$\gamma$-ray emitting AGNs (Tavecchio et al., 1998; Aleksić et al., 2014;
Ghisellini et al., 2015), we thus still consider a broken power-law
distribution for the higher-energy electron population. . $\gamma_{\min}=50$
is roughly constrained by the upper-limit of X-ray flux. We adjust the values
of $B$, $\delta$, $\gamma_{\rm b}$, and $N_{0}$ to fit the spectrum during
$\gamma$-ray outburst on the basis of the fitting result of the low-flux
state. No observational data at optical-UV band are available, and the highest
flux point in Figure1(a) observed by the _Swift_ -XRT is taken as the upper-
limit to roughly constrain the fitting parameters. We obtain $B=0.8$ G,
$\delta=20$, $\gamma_{\rm b}=1.45\times 10^{4}$, and $N_{0}=631.8$ cm-3. The
fitting results are shown in Figure 4 and the modeling parameters are listed
in Table 3. Note that $B$ and $\delta$ are degenerate (e.g., Zhang et al.,
2012), and the model parameters cannot be totally constrained with the current
observational data.

As displayed in Figure 4, the $\gamma$-rays are absolutely produced by the EC
process while the X-rays are a hybrid of several components. The observed
spectral indices at the X-ray band cluster at $1.8<\Gamma_{\rm X}<2.2$, as
shown in Figure 3(b), suggesting a transition between the synchrotron
radiation and inverse Compton component (e.g., Ghisellini et al., 1998). The
complex origins of X-rays may result in the smaller amplitude of flux
variation than that of $\gamma$-rays. No significant correlation of
variability between X-rays and $\gamma$-rays, as shown in Figure 2, also
demonstrates that the X-rays have the different origin from $\gamma$-rays.

The observation time of the $\gamma$-ray spectrum in Figure 4(b) corresponds
to the time when PKS 1413+135 was announced to be detected firstly at the VHE
band by the MAGIC telescopes. The predicted flux by the SED modeling at the
VHE band exceeds the sensitivity of MAGIC telescopes, whether considering the
EBL absorption or not. The integral flux of the model prediction over the
MAGIC sensitivity is $\sim 1.4\times 10^{-10}$ erg cm-2 s-1 after considering
the EBL absorption. Recently, Malik et al. (2022) reported that the
attenuation of VHE photons due to the EBL absorption may be overestimated,
especially for the high redshift sources. Hence, there may be a higher
observational VHE flux for PKS 1413+135 than the model prediction. We can also
observe that PKS 1413+135 would not be detectable at the VHE band when its GeV
emission is in the low-flux state, as illustrated in Figure 4(a).

The host galaxy and redshift of PKS 1413+135 are also still debated. PKS
1413+135 may be a background source of the spiral galaxy and is located at a
higher redshift of $0.247<z<0.5$ (Readhead et al., 2021). As displayed in
Figure 4(b), the predicted flux by the SED modeling after considering the EBL
absorption is still over the sensitivity of MAGIC telescopes if PKS 1413+135
is located at z=0.5. Hence, we would not give a more accurate constraint on
its redshift with the SED modeling.

##  6 Discussion

###  6.1 Origin of $\gamma$-ray Emission for PKS 1413+135

As displayed in Figure 1(b), the GeV $\gamma$-rays of PKS 1413+135 are in a
low-flux stage before MJD 58500, then undergo the obvious outbursts, and
always stay in a high-flux stage after MJD 58500. We can also observe that the
GeV spectral features are quite different for the high-flux and low-flux
states, as displayed in Figure 3(a). The curvature log-parabola energy
spectrum with a softer spectral index is shown in the low-flux state while the
spectra in the outbursts of $\gamma$-rays can be well fitted by a power-law
function with a harder spectral index. The different spectral features in the
two states for PKS 1413+135 may reflect the differently dominant acceleration
mechanisms or acceleration regions of relativistic electrons. In the radio
band, the emission above 5 GHz is totally dominated by the radio core (Table 2
in Perlman et al., 1996; Figure 4 in Peirson et al., 2022). Except the radio
core and knot-D, the other jet structures show the steep spectra, which should
originate from the optical thin synchrotron radiations, as illustrated by the
red dashed line in Figure 4(a). The inverted spectrum of radio core is due to
the SSA effect while the flat spectrum ($\alpha\sim 0$) of knot-D may indicate
a reacceleration and/or recollimation region (Perlman et al., 1996).
Component-D8 is the closest component to the radio core located at 0.32 mas
away (Peirson et al., 2022) and is also the brightest component except radio
core in the Monitoring Of Jets in Active galactic nuclei with VLBA Experiment
(MOJAVE) images (Lister et al., 2019). We speculate that the reacceleration
spectral feature of knot-D is due to the emergence of a new component-D8 from
the radio core, and the component-D8 interacts with other jet structures to
reaccelerate the electrons in knot-D. Peirson et al. (2022) reported that
component-D8 displayed the rapid variability between 1995 and 1998. The long-
term light curve of PKS 1413+135 at 14.5 GHz taken from the University of
Michigan Radio Astronomical Observatory (UMRAO, Readhead et al., 2021) shows
three outbursts around in 1982, 1988, and 1992. Maybe the emerging of
component-D8 is connected with one of the three outbursts, and this is roughly
coincident with the derived result using the separation speed (9.2$\pm$2.2
$\mu$as/yr, Lister et al., 2019) and distance of component-D8 from the radio
core. We thus suggest that the $\gamma$-rays of PKS 1413+135 may have two
origins: one is from the sub-pc-/pc-scale jet components, which contributes
the low-flux $\gamma$-rays in Figure 1(b) and the flat radio spectral
component in Figure 4(a). another is connected with the outbursts of
$\gamma$-rays, which may be due to the ejection of a new component recently
from the core, similar to many blazars (e.g., Marscher et al., 2010; Jorstad
et al., 2013; Lisakov et al., 2017; Rani et al., 2018; Lee et al., 2019; Zhang
et al., 2020), namely, the restarted nuclear jet activity of PKS 1413+135.

The $\gamma$-ray luminosity of source during outbursts is up to $\sim 10^{47}$
erg s-1 ($z=0.247$), as listed in Table 1, together with the significant
variability at the $\gamma$-ray band, implying the strong Doppler boosting
effect of the $\gamma$-ray emission region. We assume $\delta=\Gamma$ during
the SED fitting, the viewing angle ($\theta$) is equal to the opening angle
($1/\Gamma$) of jet (e.g., Chen & Zhang, 2021), and thus it is $\theta\sim
3^{\circ}$. This is conflicted with the two-side pc-scale structure in the
radio morphology of PKS 1413+135 since the symmetric radio structure of CSOs
is thought to be due to a misaligned jet to the observers (Phillips & Mutel,
1980; Wilkinson et al., 1994; Readhead et al., 1996). Recently, the episodic
nuclear jet activity in $\gamma$-ray emitting CSOs, TXS 0128+554 (Lister et
al., 2020) and CTD 135 (Gan et al., 2021), was reported. We thus propose that
there are episodic nuclear jet activities in the center of PKS 1413+135. The
axis of the recently restarted jet is aligned within a few degree to the line
of sight and is different from the direction of the pre-existing components,
similar to the typical radio galaxy 3C 84 (Nagai et al., 2016). The restarted
jet may result in a new component ejected from the core, which would interact
with the surrounding materials (or the pre-existing components) and then
accelerate particles to produce the significant $\gamma$-ray outburst and
variability. One of the most striking evidence of episodic nuclear jet
activity is that two or more pairs of distinct radio components are observed
on opposite sides of the radio core (O’Dea & Saikia, 2021). The outer
components of jet and counter-jet for PKS 1413+135, which are far away from
the radio core, such as the components A, F and G (Perlman et al., 1996), may
be the remnants of the nuclear jet activity long-ago.

The overall size of the radio morphology for PKS 1413+135 is $\sim$422 pc
(projection size, $z=0.247$, Readhead et al., 2021), and an explanation of the
size below 1 kpc for compact radio sources is that they are transient or
episodic sources (O’Dea & Saikia, 2021 for a review), which cannot provide
enough energy for the jet to propagate to large scale. On the other hand, a
fairly dense nuclear medium (Perlman et al., 1996) may also frustrate the jet
propagation of PKS 1413+135. And the interaction between the jet and
surrounding medium may change the jet direction and result in the bent in the
extended jet structures. As given in Table 3, the derived values of $\delta$
by the SED modeling are 2.2 and 20 for low-flux and high-flux states,
respectively, implying that the jet of PKS 1413+135 is decelerated effectively
on sub-pc scale. If the scenario described above is true, the recent
$\gamma$-ray outbursts of PKS 1413+135 should be connected with the restarted
nuclear jet activity. And the ejected new component from the central engine
would be resolved out with the VLBI observations in the future.

###  6.2 PKS 1413+135: a CSO or a Blazar

PKS 1413+135 was classified as either a BL Lac or a “red quasar” on account of
the inverted and rapidly variable radio spectrum, the absence of optical
emission lines, its extremely steep near-IR slope, and an IR K-band
polarization of $(16\pm 3)\%$ (Bregman et al., 1981; Stocke et al., 1992). We
find that PKS 1413+135 also displays the significant flux variation in the GeV
band accompanied by the obvious spectral evolution. As shown in Figure 1(b),
the variation of $\gamma$-ray flux for PKS 1413+135 exceeds two orders of
magnitude, indicating the typical relativistic jet effect, similar to blazars.
Especially, the report of the VHE detection for PKS 1413+135 by the MAGIC
telescopes further supports that it may be a BL Lac since most of the
confirmed extragalactic VHE emission sources are BL
Lacs555http://tevcat.uchicago.edu/. Generally, the broadband SEDs of BL Lacs
can be reproduced with the single-zone syn+SSC model (e.g., Zhang et al.,
2012; Yan et al., 2014; La Mura et al., 2022). However, this simple model is
not suitable for PKS 1413+135 since the broad second bump in its broadband SED
during $\gamma$-ray outbursts.

Different from blazars, no significant superluminal motions on pc-scale jet of
PKS 1413+135 were reported. The maximum value is $1.72\pm 0.11\ c$, i.e.,
110$\pm$7 $\mu$as/yr (Lister et al., 2019), which corresponds to the
component-D3 (in Figure 2 of Peirson et al., 2022) located at $\sim$7 mas away
from the radio core. On account of the miniature radio morphology with two-
side pc-scale structure, PKS 1413+135 was suggested to be a CSO (Perlman et
al., 1996). And the counter-jet of PKS 1413+135 is also a strong radio source
(Perlman et al., 1996; Peirson et al., 2022). In contrast, blazars are
believed to be on-axially observed to their jets and an asymmetric one-side
aligned jet is observed. As shown in Figure 2, on average, the X-ray fluxes
are lower than the $\gamma$-ray fluxes two orders of magnitude for PKS
1413+135. This feature is also not like that of the TeV emitting BL Lacs, as
displayed by the broadband SEDs of TeV-selected BL Lacs in Zhang et al.
(2012). Hence, the AGN type of PKS 1413+135 is still debated.

##  7 Summary

The miniature radio morphology with two-side pc-scale structure makes PKS
1413+135 to be classified as a CSO. Interestingly, it has been detected by
MAGIC telescopes and is the first CSO detected at VHE band. In this paper, we
comprehensively analyzed the 13.5-year _Fermi_ -LAT observation data of PKS
1413+135, together with its archive X-ray data observed with _Swift_ -XRT,
_Chandra_ , and _XMM-Newton_ , to investigate its high energy radiation
properties. The significant variability is presented in the long-term GeV
$\gamma$-ray light curve with a confidence level far beyond 5$\sigma$. And the
whole $\gamma$-ray light curve shows two distinct stages, a low-flux stage
before MJD 58500 and a high-flux stage with violent outbursts after MJD 58500.
The clear spectral variation accompanying the variation of $\gamma$-ray flux
is observed, namely, the curvature log-parabola energy spectrum with a softer
spectral index in the low-flux state and the power-law spectrum with a harder
spectral index in the $\gamma$-ray outbursts, indicating a “harder when
brighter” behavior at GeV band. The maximum energy of the detected photons by
the _Fermi_ -LAT during its $\sim$13.5-year observations is up to $\sim$236
GeV. The X-rays of PKS 1413+135 are also obviously variable, but no
significant spectral variation is found. The variability amplitude of X-rays
is much smaller than that of $\gamma$-rays, and no correlation of flux
variation between two bands is presented. Attributing the broadband SEDs of
PKS 1413+135 to the radiations from the compact radiation region close to the
black hole, the sub-pc-/pc-scale jet component, and the large-scale jet
extended region, we represented the SED by a multi-zone leptonic jet model.
Its $\gamma$-rays are contributed by the EC process of relativistic electrons
while the X-rays are synthesis of several radiation components. The different
radiation origins may result in the small variability amplitude of X-ray flux
and no correlation of variability between X-rays and $\gamma$-rays. The
predicted flux by the SED modeling at the VHE band would be detectable by
MAGIC telescopes when PKS 1413+135 is in a high-flux state at the GeV band.
Considering the features of the temporal and spectral variations at the GeV
band, together with the observations in other bands, we proposed that the
$\gamma$-ray outbursts after MJD 58500 of PKS 1413+135 may be connected with
the recently restarted nuclear jet activity of the central engine and are
produced in a aligned core-jet component. The sub-pc and pc scale jet
components, which may be produced by the previous nuclear jet activities,
contribute the low-flux $\gamma$-ray emission. This scenario can be checked by
the future VLBI observations.

We thank the anonymous referee for valuable suggestions. This work is
supported by the National Natural Science Foundation of China (grants
12022305, 11973050, and 12133003), and Guangxi Science Foundation (grants
2018GXNSFGA281007 and 2019AC20334).

Figure 1: Light curves of PKS 1413+135 in X-ray and $\gamma$-ray bands.
_Panel-(a)_ : the X-ray light curve obtained by the observations of _Swift_
-XRT (black squares), _Chandra_ (red circles), and _XMM-Newton_ (blue stars),
where the horizontal red dashed line represents the mean flux of all the data
points derived with Equation (2). _Panel-(b)_ : the $\sim$13.5-year
$\gamma$-ray light curve derived with an adaptive-binning method based on a
criterion of TS$\geq 9$ for each time bin, where the minimum time-bin step is
7 days. The horizontal blue solid line indicates the 13.5-year average flux of
$2.28\times 10^{-11}$ erg cm-2 s-1. The horizontal green short-dashed line is
the sensitivity of _Fermi_ -LAT, which is obtained and derived from the link
of https://www.slac.stanford.edu/exp/glast/groups/canda/lat_Performance.htm.
_Panel-(c)_ : the $\gamma$-ray light curve after MJD 58500 also obtained with
the adaptive-binning method and the criterion of TS$\geq 9$ for each time bin,
where the minimum time-bin step is one day. The horizontal red dashed line
indicates the average flux of $7.74\times 10^{-11}$ erg cm-2 s-1 while the
horizontal blue solid line is the same as in _Panel-(b)_. Figure 2:
Correlation of flux between $\gamma$-rays and X-rays for PKS 1413+135. The
solid line represents $F_{2-10\rm\,keV}=F_{0.1-300\rm\,GeV}/100$. The data are
taken from Figure 1.

Figure 3: _Panel-(a)_ : the spectra observed by the _Fermi_ -LAT for PKS
1413+135, including the time-integrated spectrum of the whole light curve
(magenta symbols) in Figure 1 together with the time-integrated spectra before
(MJD 54682–58500, cyan symbols) and after (MJD 58500–59660, orange symbols)
MJD 58500, the time-resolved spectra of MJD 58724-58726 (green symbols), MJD
58796–58798 (blue symbols), and MJD 59589–59591 (red symbols). The
corresponding colored solid lines represent the fitting results. The
sensitivity curves of MAGIC telescopes (pink dot-dashed line, 50 hr) and CTA-N
(blue dot-dashed line, 50 hr) are also given in the figure. _Panel-(b)_ :
$\Gamma_{\rm X}$ against $F_{2-10\rm\,keV}$, where $\Gamma_{\rm X}$ and
$F_{2-10\rm\,keV}$ are the photon spectral index and flux in the 2–10 keV
band. They are obtained with the observations of _Chandra_ (red circles) and
_XMM-Newton_ (blue stars).

Figure 4: Observed SEDs with model fitting for PKS 1413+135 in the low-flux
(_Panel-(a)_) and high-flux (_Panel-(b)_) states of $\gamma$-rays. The
$\gamma$-ray spectra in the two panels are the time-integrated spectrum before
MJD 58500 (_Panel-(a)_) and the time-resolved spectrum during MJD 59589–59591
(_Panel-(b)_), respectively. The data in radio–optical bands, the gray solid
and opened squares, are taken from the NED. The green solid and magenta opened
stars indicate the data of the radio core and knot-D, which are taken from
Perlman et al. (1996). In the X-ray band, the red inverted triangle represents
the highest flux point in Figure 1(a) observed by the _Swift_ -XRT on March
26, 2020 (MJD 58934), which is the integral flux value of $F_{2-10\rm\,keV}$.
The cyan bowtie is obtained with the lowest flux point in Figure 1(a) observed
by _Chandra_ on December 20, 2019 (MJD 58837), which indicates the flux at 6
keV with $\Gamma_{X}=1.83\pm 0.32$. The sensitivity curves of MAGIC telescopes
(pink dot-dashed line, 50 hr) and CTA-N (blue dot-dashed line, 50 hr) are also
presented. _Panel-(a)_ : The black solid line is the sum of each component
emission for the source located at z=0.247, synchrotron radiation (red lines),
SSC process (magenta lines), and EC process (cyan lines), where the solid and
dashed lines respectively represent the radiations from the sub-pc-/pc-scale
jet component and the large-scale jet extended region. An insert in the upper
left of Panel-(a) shows the fitting results for the spectra from 80 MHz to 2.7
GHz and from 4.8 GHz to 375 GHz with the power-law function of
$F_{\nu}\propto\nu^{-\alpha}$, and $\alpha\sim 0.66$ (orange solid line) and
$\alpha\sim 0.09$ (yellow solid line) are obtained, respectively. _Panel-(b)_
: The dark-yellow dashed line is the same line as the black solid line in
_Panel-(a)_. The colored solid lines represent the different radiation
components of the compact $\gamma$-ray emission region, synchrotron radiation
(red line), SSC process (magenta line), and EC process (cyan line), where the
cyan line at the VHE band indicates the intrinsic flux of the source without
considering the EBL absorption. The black solid line is the sum of the dark-
yellow dashed line and three colored solid lines after considering the EBL
absorption for the source located at z=0.247. The green dot-dashed line
represents the fitting result when the source is located at $z=0.5$.

Table 1: _Fermi_ -LAT Analysis Results for PKS 1413+135 Obs-date | Modelaa“LP” and “PL” represent the spectra are fitted by the Log-Parabola (LP) and Power-Law (PL) functions, respectively. | $\Gamma_{\gamma}$ | $\beta$ | $F_{\gamma}$ | $L_{\gamma}\tablenotemark{\scriptsize{b}}$ | TS | $E_{\rm max}$ | Det-timeccThe detection time of the maximum energy photon by the _Fermi_ -LAT.
---|---|---|---|---|---|---|---|---
(MJD) |  |  |  | (erg cm-2 s-1) | (erg s-1) |  | (GeV) | (MJD)
54682–59660 | LP | $2.02\pm 0.03$ | $0.03\pm 0.01$ | $(2.28\pm 0.13)\times 10^{-11}$ | $(4.13\pm 0.23)\times 10^{45}$ | 3355 | 235.9 | 58798
54682–58500 | LP | $2.34\pm 0.08$ | $0.08\pm 0.05$ | $(5.74\pm 0.61)\times 10^{-12}$ | $(1.04\pm 0.11)\times 10^{45}$ | 300 | 20.6 | 57959
58500–59660 | LP | $1.90\pm 0.03$ | $0.04\pm 0.01$ | $(7.74\pm 0.48)\times 10^{-11}$ | $(1.41\pm 0.09)\times 10^{46}$ | 4704 | 235.9 | 58798
58724–58726 | PL | $1.75\pm 0.11$ | 0.00 (fixed) | $(1.07\pm 0.39)\times 10^{-09}$ | $(1.94\pm 0.71)\times 10^{47}$ | 200 | 55.7 | 58726
58796–58798 | PL | $1.63\pm 0.14$ | 0.00 (fixed) | $(8.59\pm 4.02)\times 10^{-10}$ | $(1.56\pm 0.73)\times 10^{47}$ | 102 | 235.9 | 58798
59589–59591 | PL | $1.73\pm 0.13$ | 0.00 (fixed) | $(7.95\pm 3.43)\times 10^{-10}$ | $(1.44\pm 0.62)\times 10^{47}$ | 137 | 25.7 | 59590
54682–59063ddThe analysis results given in the 4FGL-DR3 (Abdollahi et al., 2022). | LP | $2.00\pm 0.04$ | $0.03\pm 0.01$ | $(1.94\pm 0.09)\times 10^{-11}$ |  |  |  |
Table 2: Public X-ray observations of PKS 1413+135 Obs-date |  | Mission |  | Exp. (s)
---|---|---|---|---
(1) |  | (2) |  | (3)
2007-03-20 |  | _Chandra_ |  | 20180
2019-12-20 |  | _Chandra_ |  | 4060
2019-12-29 |  | _XMM-Newton_ |  | 13000
2020-02-03 |  | _Chandra_ |  | 4060
2020-02-05 |  | _XMM-Newton_ |  | 17000
2020-03-22 |  | _Chandra_ |  | 4060
2020-05-04 |  | _Chandra_ |  | 4060
2020-06-18 |  | _Chandra_ |  | 4060
2020-07-04 |  | _XMM-Newton_ |  | 15000
2007 to 2022 |  | _Swift_ |  | 102 to 104
Table 3: SED Fitting Parameters of PKS 1413+135 Stateaa“E”, “L”, and “H” represent the parameters for the large-scale jet extended region, the sub-pc-/pc-scale jet component, and the compact radiation region close to the black hole, respectively. The parameter values correspond to the result for the source located at $z=0.247$. | $R$ | $R_{\rm diss}\tablenotemark{\scriptsize{b}}$ | $B$ | $\delta$ | $\Gamma$ | $\gamma_{\min}$ | $\gamma_{\rm b}$ | $\gamma_{\max}$ | $N_{0}$ | $p_{1}$ | $p_{2}$
---|---|---|---|---|---|---|---|---|---|---|---
|  |  | (mG) |  |  |  |  |  | (cm-3) |  |
E | 211 pc | 400 pc | 1.25 | 1 | 1 | 1 |  $\cdots$ | 2000 | 0.036 | 2.32 |  $\cdots$
L | 0.6 pc | 1 pc | 450 | 2.2 | 2.2 | 1 | 2609 | $10^{5}$ | 2.81 | 1.8 | 3.84
H |  $4.15\times 10^{16}$ cm |  $4.15\times 10^{17}$ cm | 800 | 20 | 20 | 50 | 14513 | $10^{6}$ | 631.8 | 2.46 | 3.84

## References

  * Abdollahi et al. (2020) Abdollahi, S., Acero, F., Ackermann, M., et al. 2020, ApJS, 247, 33.
  * Aleksić et al. (2016) Aleksić, J., Ansoldi, S., Antonelli, L. A., et al. 2016, Astroparticle Physics, 72, 76.
  * Aleksić et al. (2014) Aleksić, J., Antonelli, L. A., Antoranz, P., et al. 2014, A&A, 563, A90.
  * Blanch et al. (2022) Blanch, O., Sitarek, J., & Striskovic, J. 2022, The Astronomer’s Telegram, 15161
  * Bregman et al. (1981) Bregman, J. N., Lebofsky, M. J., Aller, M. F., et al. 1981, Nature, 293, 714.
  * Carilli et al. (1992) Carilli, C. L., Perlman, E. S., & Stocke, J. T. 1992, ApJ, 400, L13.
  * Carvalho (1994) Carvalho, J. C. 1994, A&A, 292, 392.
  * Carvalho (1998) Carvalho, J. C. 1998, A&A, 329, 845.
  * Chen et al. (2022) Chen, S., Laor, A., & Behar, E. 2022, MNRAS, 515, 1723.
  * Chen & Zhang (2021) Chen, L. & Zhang, B. 2021, ApJ, 906, 105.
  * Abdollahi et al. (2022) Abdollahi, S., Acero, F., Baldini, L., et al. 2022, ApJS, 260, 53.
  * Fey et al. (2015) Fey, A. L., Gordon, D., Jacobs, C. S., et al. 2015, AJ, 150, 58.
  * Franceschini et al. (2008) Franceschini, A., Rodighiero, G., & Vaccari, M. 2008, A&A, 487, 837.
  * Fu et al. (2022) Fu, W.-J., Zhang, H.-M., Zhang, J., et al. 2022, Research in Astronomy and Astrophysics, 22, 035005.
  * Gan et al. (2021) Gan, Y.-Y., Zhang, H.-M., Zhang, J., et al. 2021, Research in Astronomy and Astrophysics, 21, 201.
  * Ghisellini et al. (1998) Ghisellini, G., Celotti, A., Fossati, G., et al. 1998, MNRAS, 301, 451.
  * Ghisellini et al. (2015) Ghisellini, G., Haardt, F., Ciardi, B., et al. 2015, MNRAS, 452, 3457.
  * Ghisellini & Tavecchio (2009) Ghisellini, G. & Tavecchio, F. 2009, MNRAS, 397, 985.
  * Gu et al. (2022) Gu, Y., Zhang, H.-M., Gan, Y.-Y., et al. 2022, ApJ, 927, 221.
  * Heckman et al. (1994) Heckman, T. M., O’Dea, C. P., Baum, S. A., et al. 1994, ApJ, 428, 65.
  * HI4PI Collaboration et al. (2016) HI4PI Collaboration, Ben Bekhti, N., Fl?er, L., et al. 2016, A&A, 594, A116.
  * Jorstad et al. (2013) Jorstad, S. G., Marscher, A. P., Smith, P. S., et al. 2013, ApJ, 773, 147.
  * Kang et al. (2014) Kang, S.-J., Chen, L., & Wu, Q. 2014, ApJS, 215, 5.
  * La Mura et al. (2022) La Mura, G., Becerra Gonzalez, J., Chiaro, G., et al. 2022, MNRAS.
  * Lee et al. (2019) Lee, T., Trippe, S., Kino, M., et al. 2019, MNRAS, 486, 2412.
  * Lisakov et al. (2017) Lisakov, M. M., Kovalev, Y. Y., Savolainen, T., et al. 2017, MNRAS, 468, 4478.
  * Liu & Bai (2006) Liu, H. T. & Bai, J. M. 2006, ApJ, 653, 1089.
  * Lister et al. (2020) Lister, M. L., Homan, D. C., Kovalev, Y. Y., et al. 2020, ApJ, 899, 141.
  * Lister et al. (2019) Lister, M. L., Homan, D. C., Hovatta, T., et al. 2019, ApJ, 874, 43.
  * Malik et al. (2022) Malik, Z., Sahayanathan, S., Shah, Z., et al. 2022, MNRAS, 515, 4505.
  * Marscher et al. (2010) Marscher, A. P., Jorstad, S. G., Larionov, V. M., et al. 2010, ApJ, 710, L126.
  * Massaro et al. (2004) Massaro, E., Perri, M., Giommi, P., et al. 2004, A&A, 413, 489.
  * McHardy et al. (1994) McHardy, I. M., Merrifield, M. R., Abraham, R. G., et al. 1994, MNRAS, 268, 681.
  * McHardy et al. (1991) McHardy, I. M., Abraham, R. G., Crawford, C. S., et al. 1991, MNRAS, 249, 742.
  * McLaughlin et al. (1996) McLaughlin, M. A., Mattox, J. R., Cordes, J. M., et al. 1996, ApJ, 473, 763.
  * Migliori et al. (2016) Migliori, G., Siemiginowska, A., Sobolewska, M., et al. 2016, ApJ, 821, L31.
  * Müller et al. (2015) Müller, C., Krauß, F., Dauser, T., et al. 2015, A&A, 574, A117.
  * Müller et al. (2014) Müller, C., Kadler, M., Ojha, R., et al. 2014, A&A, 562, A4.
  * Nagai et al. (2016) Nagai, H., Chida, H., Kino, M., et al. 2016, Astronomische Nachrichten, 337, 69.
  * NASA/IPAC Extragalactic Database (2019) NASA/IPAC Extragalactic Database (NED). 2019, IPAC, doi: 10.26132/NED1
  * Nolan et al. (2012) Nolan, P. L., Abdo, A. A., Ackermann, M., et al. 2012, ApJS, 199, 31.
  * O’Dea et al. (1991) O’Dea, C. P., Baum, S. A., & Stanghellini, C. 1991, ApJ, 380, 66.
  * O’Dea & Saikia (2021) O’Dea, C. P. & Saikia, D. J. 2021, A&A Rev., 29, 3.
  * Peirson et al. (2022) Peirson, A. L., Liodakis, I., Readhead, A. C. S., et al. 2022, ApJ, 927, 24.
  * Perlman et al. (2002) Perlman, E. S., Stocke, J. T., Carilli, C. L., et al. 2002, AJ, 124, 2401.
  * Perlman et al. (1996) Perlman, E. S., Carilli, C. L., Stocke, J. T., et al. 1996, AJ, 111, 1839.
  * Perlman et al. (1994) Perlman, E. S., Stocke, J. T., Shaffer, D. B., et al. 1994, ApJ, 424, L69.
  * Phillips & Mutel (1980) Phillips, R. B. & Mutel, R. L. 1980, ApJ, 236, 89.
  * Principe et al. (2020) Principe, G., Migliori, G., Johnson, T. J., et al. 2020, A&A, 635, A185.
  * Principe et al. (2021) Principe, G., Di Venere, L., Orienti, M., et al. 2021, MNRAS, 507, 4564.
  * Rani et al. (2018) Rani, B., Jorstad, S. G., Marscher, A. P., et al. 2018, ApJ, 858, 80.
  * Readhead et al. (2021) Readhead, A. C. S., Ravi, V., Liodakis, I., et al. 2021, ApJ, 907, 61.
  * Readhead et al. (1996) Readhead, A. C. S., Taylor, G. B., Xu, W., et al. 1996, ApJ, 460, 612.
  * Stawarz et al. (2008) Stawarz, Ł., Ostorero, L., Begelman, M. C., et al. 2008, ApJ, 680, 911.
  * Stocke et al. (1992) Stocke, J. T., Wurtz, R., Wang, Q., et al. 1992, ApJ, 400, L17.
  * Tavecchio et al. (1998) Tavecchio, F., Maraschi, L., & Ghisellini, G. 1998, ApJ, 509, 608.
  * Vedantham et al. (2017) Vedantham, H. K., Readhead, A. C. S., Hovatta, T., et al. 2017, ApJ, 845, 89.
  * Wilkinson et al. (1994) Wilkinson, P. N., Polatidis, A. G., Readhead, A. C. S., et al. 1994, ApJ, 432, L87.
  * Yan et al. (2014) Yan, D., Zeng, H., & Zhang, L. 2014, MNRAS, 439, 2933.
  * Zhang et al. (2020) Zhang, J., Zhang, H.-M., Gan, Y.-Y., et al. 2020, ApJ, 899, 2.
  * Zhang et al. (2012) Zhang, J., Liang, E.-W., Zhang, S.-N., et al. 2012, ApJ, 752, 157.