# A gamma-ray flare from TXS 1508+572: characterizing the jet of a $z=4.31$
blazar in the early Universe

Andrea Gokus Department of Physics & McDonnell Center for the Space Sciences,
Washington University in St. Louis, One Brookings Drive, St. Louis, MO 63130,
USA  Markus Böttcher Centre for Space Research, North-West University,
Potchefstroom 2520, South Africa  Manel Errando Department of Physics &
McDonnell Center for the Space Sciences, Washington University in St. Louis,
One Brookings Drive, St. Louis, MO 63130, USA  Michael Kreter Centre for Space
Research, North-West University, Potchefstroom 2520, South Africa  Jonas
Heßdörfer Julius-Maximilians-Universität Würzburg, Fakultät für Physik und
Astronomie, Institut für Theoretische Physik und Astrophysik, Lehrstuhl für
Astronomie, Emil-Fischer-Str. 31, D-97074 Würzburg, Germany  Florian Eppel
Julius-Maximilians-Universität Würzburg, Fakultät für Physik und Astronomie,
Institut für Theoretische Physik und Astrophysik, Lehrstuhl für Astronomie,
Emil-Fischer-Str. 31, D-97074 Würzburg, Germany  Matthias Kadler Julius-
Maximilians-Universität Würzburg, Fakultät für Physik und Astronomie, Institut
für Theoretische Physik und Astrophysik, Lehrstuhl für Astronomie, Emil-
Fischer-Str. 31, D-97074 Würzburg, Germany  Paul S. Smith Steward Observatory,
The University of Arizona, Tucson, AZ 85721, USA  Petra Benke Max-Planck-
Institut für Radioastronomie, Auf dem Hügel 69, D-53121 Bonn, Germany  Julius-
Maximilians-Universität Würzburg, Fakultät für Physik und Astronomie, Institut
für Theoretische Physik und Astrophysik, Lehrstuhl für Astronomie, Emil-
Fischer-Str. 31, D-97074 Würzburg, Germany  Leonid I. Gurvits Joint Institute
for VLBI ERIC (JIVE), Oude Hoogeveensedijk 4, 7991 PD, Dwingeloo, The
Netherlands  Faculty of Aerospace Engineering, Delft University of Technology,
Kluyverweg 1, 2629 HS, Delft, The Netherlands  Alex Kraus Max-Planck-Institut
für Radioastronomie, Auf dem Hügel 69, D-53121 Bonn, Germany  Mikhail Lisakov
Instituto de Física, Pontificia Universidad Católica de Valparaíso, Casilla
4059, Valparaíso, Chile  Max-Planck-Institut für Radioastronomie, Auf dem
Hügel 69, D-53121 Bonn, Germany  Felicia McBride Department of Physics and
Astronomy, Bowdoin College, Brunswick, ME 04011, USA  Eduardo Ros Max-Planck-
Institut für Radioastronomie, Auf dem Hügel 69, D-53121 Bonn, Germany  Florian
Rösch Julius-Maximilians-Universität Würzburg, Fakultät für Physik und
Astronomie, Institut für Theoretische Physik und Astrophysik, Lehrstuhl für
Astronomie, Emil-Fischer-Str. 31, D-97074 Würzburg, Germany  Max-Planck-
Institut für Radioastronomie, Auf dem Hügel 69, D-53121 Bonn, Germany  Jörn
Wilms Remeis Observatory & Erlangen Centre for Astroparticle Physics,
Universität Erlangen-Nürnberg, Sternwartstr. 7, 96049 Bamberg, Germany

###### Abstract

Blazars can be detected from very large distances due to their high
luminosity. However, the detection of $\gamma$-ray emission of blazars beyond
$z=3$ has only been confirmed for a small number of sources. Such observations
probe the growth of supermassive black holes close to the peak of star
formation in the history of galaxy evolution. As a result from a continuous
monitoring of a sample of 80 $z>3$ blazars with Fermi-LAT, we present the
first detection of a $\gamma$-ray flare from the $z=4.31$ blazar TXS 1508+572.
This source showed high $\gamma$-ray activity from February to August 2022,
reaching a peak luminosity comparable to the most luminous flares ever
detected with Fermi-LAT. We conducted a mutliwavelength observing campaign
involving XMM-Newton, Swift, the Effelsberg 100-m radio telescope and the Very
Long Baseline Array. In addition, we make use of the monitoring programs by
the Zwicky Transient Facility and NEOWISE at optical and infrared wavelengths,
respectively. We find that the source is particularly variable in the infrared
band on daily time scales. The spectral energy distribution collected during
our campaign is well described by a one-zone leptonic model, with the
$\gamma$-ray flare originating from an increase of external Compton emission
as a result of a fresh injection of accelerated electrons.

Blazars (164) — Gamma-ray astronomy (628) — High-redshift galaxies (734) —
High energy astrophysics (739) — Relativistic jets (1390) — Radiative
processes (2055) — Flat-spectrum radio quasars (2163)

††journal: ApJ††facilities: Fermi-LAT, NuSTAR, Swift/XRT, XMM-Newton, Zwicky
Transient facility, NEOWISE, 100-m Effelsberg radio telescope††software:
fermipy (Wood et al., 2017), ISIS (Houck & Denicola, 2000)

##  1 Introduction

Among jetted active galactic nuclei (AGN), those with their jet pointing
towards the Earth are called blazars and appear particularly variable and
luminous due to relativistic beaming (Urry & Padovani, 1995). Blazars can be
broadly classified as either BL Lacertae (BL Lac) objects or flat-spectrum
radio quasars (FSRQs), for which the distinction is based on the existence of
optical emission lines with an equivalent width of at least 5Å(FSRQ) or the
lack thereof (BL Lacs).

Due to their extreme luminosities, we are able to detect quasars and blazars
at large distances, currently up to $z=7.6$ (Wang et al., 2021), which enables
glimpses into the early Universe. Among those found at high redshift, some
seem to be powered by the heaviest specimens of black holes, exceeding a
billion solar masses (e.g., Lobanov et al., 2001; Ghisellini et al., 2010;
Ackermann et al., 2017; Belladitta et al., 2022; Burke et al., 2024). Given
that these extremely massive black holes appear to exist only about one
billion years after the Big Bang, the circumstances under which black holes
can grow so fast are not yet understood (see, e.g., Inayoshi et al., 2020, for
a recent review). While our observations are biased towards finding the most
luminous and therefore extreme sources, current estimations for their number
densities present a hurdle for the application of our existing models that
describe black hole formation and growth (e.g., Johnson & Haardt, 2016). Some
of the suggestions for black hole formation for heavy seeds are remnants of
the supernovae of massive first-generation stars (Madau & Rees, 2001), or
matter collapsing into SMBHs right away (e.g., Begelman et al., 2006).
However, even with black hole seeds of $\sim$ 100 M⊙ the time scales needed
for SMBHs to grow to $>10^{9}$ M⊙ are too long even when assuming accretion at
the Eddington limit throughout the entire time. Alexander & Natarajan (2014)
have suggested that the first black holes might have been able to grow fast
through supra-exponential accretion. A recent study by Lai et al. (2024)
argues that SMBHs at $z\sim 5$ are only growing slowly, thereby requiring
either initial seed masses $>10^{8}$ M⊙ or extremely rapid growth at $z>5$.
Simulations of rapidly spinning black holes accreting above the Eddington
limit have revealed that their growth is accompanied by the formation of
powerful jets (McKinney et al., 2014; Sądowski et al., 2014). Hence, high-
redshift ($z>3$) blazars are ideal targets in order to learn more about the
growth of SMBHs, while also taking into account that their powerful large-
scale jets influence their host galaxies and galaxy clusters, which affects
galaxy evolution.

In order to properly assess the physical properties of these blazars, it is
essential to obtain a multiwavelength data set for these objects, in
particular, information about the high-energy emission that is needed to
measure the full power of the jet. While at X-ray energies it is possible to
probe distances up to $z>6$ (e.g., Belladitta et al., 2020; Medvedev et al.,
2020; Migliori et al., 2023), the faint $\gamma$-ray emission from such
distant sources is significantly affected by attenuation due to extragalactic
background light. Studies on high-z blazars are currently limited to a small
sample of objects with $z\geq 3$. The current version of the AGN catalog based
on data from the Fermi Large Area Telescope (LAT) lists 11 (Ajello et al.,
2020, 2022), and other studies have identified three more high-$z$ blazars
(Liao et al., 2018; Kreter et al., 2020).

A blazar’s broadband spectral energy distribution (SED) consists of two broad,
non-thermal components, where the low-energy component can be generally
attributed to leptonic synchrotron emission, while the high-energy component
can be explained either by leptonic Inverse Compton processes (e.g., Maraschi
et al., 1992; Dermer & Schlickeiser, 1993; Sikora et al., 1994; Blandford &
Levinson, 1995; Bloom & Marscher, 1996; Boettcher & Schlickeiser, 1997;
Błażejowski et al., 2000), but also additional hadronic interactions such as
photon-pion interactions and proton synchrotron emission (e.g., Mannheim &
Biermann, 1992; Mannheim, 1993; Aharonian, 2000; Mücke & Protheroe, 2001;
Aharonian, 2002; Mücke et al., 2003; Böttcher et al., 2013). In addition, some
blazars feature a thermal component from their accretion disks, which are
particularly present in high-z sources for which the emission from the disk is
redshifted to optical and even infrared wavelengths. As the entire emission
from high-$z$ blazars becomes redshifted, the peak of the high-energy
component drifts to MeV energies, justifying the term ‘MeV blazar’ (Bloemen et
al., 1995). As a result, the $\gamma$-ray spectra of these objects are steep
(Ackermann et al., 2017).

In Section 2, we outline our strategy of detecting $\gamma$-ray flares from
high-$z$ blazars and introduce TXS 1508+572. Our data reduction is described
in Section 3. We compute the $\gamma$-ray luminosity displayed by TXS 1508+572
during its peak in Section 4.1 and assess the changes in the $\gamma$-ray,
X-ray, and radio spectra in Section 4.2. In addition, our results of the
multiwavelength data analysis include a variability study of TXS 1508+572
during the flaring state as well as some long-term variability in Section 4.3
and a modeling of the broadband SED in Section 4.4. Finally, we discuss all
our findings and conclude in Section 5. Throughout this paper, we assume a
flat cosmology following Planck Collaboration et al. (2016), which reported
H0=67.8 km/s, $\Omega_{\lambda}=0.692$, and $\Omega_{M}=0.308$. With those
parameters, TXS 1508+572 has a luminosity distance of $\sim 40$ Gpc.

##  2 Monitoring campaign

While several studies have been collecting multiwavelength data for MeV
blazars and have modeled their SEDs using non- or semi-contemporaneous data
(e.g., Paliya et al., 2016; Ackermann et al., 2017; Marcotulli et al., 2020),
some have reported on the analysis of flares from high-$z$ blazars using
quasi-simultaneous data sets (albeit some sources with $2<z<3$; Orienti et
al., 2014; Paliya et al., 2019; Liao et al., 2019). In order to be able to
catch a high-$z$ blazar during a $\gamma$-ray flare, we use Fermi-LAT’s superb
all-sky monitoring capabilities and set up a $\gamma$-ray monitoring program
for 80 sources for which the fifth version of the ROMA-BZCAT (5BZCAT Massaro
et al., 2009, 2015) lists $z>3$. Our method relies on a significant detection
(i.e., TS $\geq 25$, where the test statistic
$\mathrm{TS}=2\Delta\mathrm{log}(\mathcal{L})$ compares the likelihood with
and without a source; see Sect. 3.1 for details) for 30-day average time bins,
for which we perform a daily check. Using 30-day average time bins has been
proven successful in the search for high-$z$ blazars using archival data
(Kreter et al., 2020).

In order to obtain a simultaneous multiwavelength data set upon flare
detection, we set up follow-up observations using multiple facilities (100-m
Effelsberg radio telescope, XMM-Newton, Swift) and also make use of exisiting
all-sky survey facilities such as ZTF and NEOWISE.

In this work, we report on the first $\gamma$-ray flare detected for a $z>4$
blazar, TXS 1508+572, which occured in February 2022 and marked the beginning
of a high-activity phase lasting roughly 6 months.

###  2.1 TXS 1508+572

At a redshift of $z=4.31$ (Schneider et al., 2007, but note that the first
report of its redshift at $z=4.30$ was done by Hook et al. (1995)), the blazar
TXS 1508+572 is the most distant $\gamma$-ray emitter as reported by the
Fermi-LAT collaboration (Ajello et al., 2022) and the third distant
$\gamma$-ray source to have been detected with Fermi-LAT data as reported by
Liao et al. (2018) and Kreter et al. (2020). As a $\gamma$-ray emitter, the
source is known as 4FGL J1510.1+5702, but historically has been referred to as
GB 1508+5714. The source was first studied at X-ray energies with Einstein
data (Mathur & Elvis, 1995), and subsequently studied by ASCA (Moran &
Helfand, 1997). Using Chandra’s high-resolution ACIS detector, Siemiginowska
et al. (2003) and Yuan et al. (2003) independently reported on the
identification of an X-ray jet seen as extended emission from the source. The
first very long baseline interferometry (VLBI) high-resolution image, taken
with the European VLBI Network (EVN) at 5 GHz with the angular resolution of
$\sim 5$ mas in 1995, did not reveal any extended radio emission (Frey et al.,
1997). However, Very Large Array (VLA) snapshot observations taken at 1.4 GHz
were able to detect a low-brightness jet extending to $\sim 2.5$ arcsec toward
south-west (Cheung, 2004). O’Sullivan et al. (2011a) reported an optically
thin jet component at 5 GHz and 8.4 GHz roughly 2 mas to the south of the core
based on global VLBI polarimetry observations. A recent observation of the
source with LOFAR revealed emission from the counter-jet seen at 144 MHz with
sub-arcsecond spatial resolution (Kappes et al., 2022). On milli-arcsecond
scales, Titov et al. (2023) monitored TXS 1508+572 over four years at 2.3 and
8.4 GHz in order to study the apparent absolute astrometric proper motion and
detected a jet proper motion of $0.117\pm 0.078$ mas/yr.

In 2017, the object was confirmed as a $\gamma$-ray emitter (Ackermann et al.,
2017), and variability on monthly time scales has been detected (Li et al.,
2018). A period of enhanced brightness in both the $\gamma$-ray and the
optical band has been reported by Liao et al. (2020).

On 4 February 2022, we detected a $\gamma$-ray flare of TXS 1508+572 through
our Fermi-LAT monitoring pipeline and found that the source had significantly
brightened compared to the flux reported in the Fermi-LAT Fourth Source
Catalog (4FGL; Abdollahi et al., 2020; Gokus et al., 2022). We coordinated a
multiwavelength campaign to obtain simultaneous data during the flaring state,
and in addition launched a VLBI monitoring campaign using the Very Long
Baseline Array (VLBA) and the 100-m Effelsberg radio telescope, whose results
we report in Benke et al. (2024).

##  3 Multiwavelength data

Figure 1: Multiwavelength light curve of TXS 1508+572 covering a time span
from the end of 2021 to the beginning of 2023. The panels from top to bottom
show the fluxes, magnitudes, or flux densities in different, decreasing energy
ranges, starting with $\gamma$-rays at the top. The two bottom panels
illustrate changes in the single-dish radio polarization and EVPA,
respectively, measured at 60 mm and 20 mm. A dotted line marks the detection
of the flaring state of TXS 1508+572, while the three dashed lines mark the
VLBI observations presented in Benke et al. (2024).

In this section, we report on the broadband data extraction and analysis, and
show the fluxes, magnitudes, and flux densities that have been compiled in
light curves for different energy bands in Fig. 1.

###  3.1 Fermi-LAT

The Large Area Telescope (LAT) onboard the Fermi satellite (Atwood et al.,
2009) has monitored the entire sky since August 2008, and has detected 3814
AGN in 12 years (Ajello et al., 2020, 2022), with the large majority being
blazars. We extract the LAT data using ScienceTools Version 1.2.23 and fermipy
0.20.0 (Wood et al., 2017), following the standard data reduction process.
Events with an energy between 100 MeV and 300 GeV are selected within a region
of interest (ROI) of 15 degrees around the $\gamma$-ray source coordinates
given in the Fermi-LAT Fourth Source Catalog Paper (4FGL; Abdollahi et al.,
2020), and filtered by selecting SOURCE class events, and those fulfilling
DATA_QUAL$>$0 and LAT_CONFIG==1. In addition, we use only events that enter
the instrument with a maximum zenith angle of $90^{\circ}$ to avoid
contamination by $\gamma$-rays from Earth-limb effects. To model the
background, we use gll_iem_v07 as the Galactic diffuse model, and
iso_P8R3_SOURCE_V2_v1 to account for the isotropic diffusion emission.
Furthermore, we use the post-launch instrument response function
P8R3_SOURCE_V2. In addition to the background, we include all known
$\gamma$-ray sources listed in the 4FGL that lie within $20^{\circ}$ of TXS
1508+572, and use the spectral type listed in the 4FGL as model for each
source. In the case of our target source, this model is represented by a power
law. Optimization of our model parameters is done via a maximum likelihood
analysis, for which the significance of the $\gamma$-ray emission for each
source is determined via the test statistic
$\mathrm{TS}=2\Delta\mathrm{log}(\mathcal{L})$. The likelihood function
$\mathcal{L}$ describes the difference between a model with and without a
source at the given coordinates (Mattox et al., 1996).

The $\gamma$-ray spectra of TXS 1508+572 are created for the quiescent and
flaring states, which cover a time range of MJD 55197 to 59215 (2010-01-01 to
2021-01-01), and MJD 59610 to 59624 (2022-01-31 to 2022-02-14), respectively.
In the spectral fits, we keep spectral and normalization parameters free for
our target source and all sources within $3^{\circ}$ of it. For sources within
$5^{\circ}$, or with $\mathrm{TS}>500$, we leave the normalization free to
vary. Galactic and isotropic diffuse emissions are also kept free during the
modeling.

Moreover, we compute a light curve from MJD 59434 to 60094 (2021–08–08 to
2023–05–30) with 30-day binning, and MJD 59582 to MJD 60051 (2022–01–03 to
2023–04–17) with 7-day binning, which are shown in the top two panels of Fig.
1 for the 30-day and 7-day binning, respectively. The 30-day binning is chosen
such that the first TS$\geq 25$ bin coincides with the detection of increased
activity with our pipeline, i.e., on 4 February 2022. We keep the
normalization free for all sources within $3^{\circ}$ of TXS 1508+572, as well
as sources with $\mathrm{TS}>500$ over the entire time range. In addition, we
keep the spectral index of our target source as a free parameter. All diffuse
emission is kept frozen to the best-fit parameters determined in the modeling
of the entire time range. We display the $2\sigma$ upper limit for bins with
$\mathrm{TS}<25$ for the 30-day binned light curve, and $\mathrm{TS}<9$ for
the 7-day binned light curve, which equates to a detection significance of
$\sim 5\sigma$ and $\sim 3\sigma$, respectively. The uncertainties of the flux
bins are given as $1\sigma$ uncertainties.

###  3.2 X-ray data

We model all X-ray data using the Interactive Spectral Interpretation System
(ISIS, Version 1.6.2-51; Houck & Denicola, 2000) in order to determine
spectral parameters and the source flux, and utilize the C-statistics (Cash,
1979). Uncertainties are given at the $1\sigma$ confidence level. The
absorption by the interstellar medium is based on the vern cross sections
(Verner et al., 1996) and the wilm abundances (Wilms et al., 2000). For each
spectrum, we use an absorbed, pegged power law (tbabs*pegpwrlw), except for a
simultaneous fit of XMM-Newton and NuSTAR spectra, which we utilize to test
for a spectral break. We freeze the Galactic hydrogen absorption to the value
from the HI 4$\pi$ survey (HI4PI; HI4PI Collaboration et al., 2016), which is
$N_{\mathrm{H,Gal}}=1.55\times 10^{20}\,\mathrm{cm}^{-2}$. The results from
our X-ray spectral fits are presented in Sect. 4.2. The fluxes measured with
XMM-Newton and Swift and obtained from the respective best fits are shown in
the third panel from the top in Fig. 1.

####  3.2.1 XMM-Newton

We obtained a target-of-opportunity (ToO) observation with the XMM-Newton
observatory (Jansen et al., 2001) to follow up the flare as soon as possible
after the detection of the $\gamma$-ray flare. TXS 1508+572 was observed on 8
February 2022 (ObsID 0910390101) for 85 ks in the Small Window mode with the
European Photon Imaging Camera (EPIC) with the pn (Strüder et al., 2001) and
both MOS detectors (Turner et al., 2001). In addition, we extract the data
from an archival observation performed on 11 May 2002 (ObsID 0111260201),
which was taken in Full Window mode with an exposure time of 24 ks. For the
archival data, we were only able to obtain spectra from the EPIC-pn and EPIC-
MOS2 detector since the source region coincides with a CCD gap for the EPIC-
MOS1 detector. We extract spectra for both observations using the XMM-Newton
Science Analysis Software (SAS, Version 20.0.0) by using the standard methods
in order to process the observation data files and produce calibrated event
lists and images. Our source and background regions for both observations and
the respective detector are listed in the Appendix in Table 6. We extract
single and double events for the data taken with EPIC-pn, and all events for
data taken with the EPIC-MOS detectors. In both observations, pile-up is
negligible. We extract light curves of the observation in three different
energy bands using the EPIC-pn detector: for a full energy range from 0.3 keV
to 10 keV, as well as for a soft band from 0.3 keV to 2 keV and for a hard
band from 2 keV to 10 keV. We find that at the beginning and end of the
observation, the data are impacted by severe background flaring and exclude
those time ranges from any data analyses. For the remaining time range, we
compute a subtracted light curve by subtracting the background light curve
from the light curve obtained for the source region. To correct for detector
and scaling effects, we only work with light curves produced with the XMMSAS
function epiclccorr.

####  3.2.2 Neil Gehrels Swift Observatory

We requested several observations with the Neil Gehrels Swift Observatory
(Gehrels et al., 2004) over the course of a year, to gather contemporaneous
data to the VLBI monitoring with VLBA and the 100-m Effelsberg radio telescope
(Benke et al., 2024). In addition to the observations obtained in relation to
the $\gamma$-ray flare, an archival observation of TXS 1508+572 exists, which
was taken simultaneous to a NuSTAR observation in 2017. We use those data for
creating a quiescent-state SED of the source. The Swift-XRT data are extracted
using the standard procedures with the xrtpipeline (Version 0.13.7) with
HEASOFT 6.30. The source spectra are compiled from the source coordinates
(R.A. = 227.5114780∘, Decl. = 57.0447472∘) and a radius of 30 arcsec, while
the background is created using an annulus centered on the same coordinates,
with an inner radius of 40 arcsec and an outer radius of 150 arcsec.

We do not detect TXS 1508+572 with the optical/UV telescope onboard Swift.

####  3.2.3 NuSTAR

We utilize an archival NuSTAR observation taken on 30 April 2017 (ObsID
60201013002), with the exposure time of 73 ks, to obtain hard X-ray data for
the quiescent state SED. We extract the data using the standard methods with
NUSTARDAS (Version 2.1.2) and CALBD version 20230124. Using nupipeline
(version 0.4.9) we reduce data from both Focal Plane Modules A and B (FPMA and
FPMB), and create calibrated event lists and images. We choose a circular
region with a radius of 50 arcsec at R.A. = 227.5131365∘, Decl. = 57.0465406∘
and R.A. = 227.5138372∘, Decl. = 57.0473028∘ for the source region in FPMA and
FPMB, respectively. The background regions are chosen to be circular as well
with a radius of 160 arcsec, and centered on R.A. = 227.4806152∘, Decl. =
57.1443941∘, and R.A. = 227.4856466∘, Decl. = 57.1484921∘ for FPMA and FPMB,
respectively. After the data extraction, the resulting spectra have an
exposure time of $\sim 37$ ks.

###  3.3 Optical & Infrared archival data

Through Lyman-$\alpha$ ($\lambda$1215.67 Å) absorption, the intergalactic
medium affects the optical emission from distant sources (e.g., Gunn &
Peterson, 1965). With a redshift of $z=4.31$, the optical light of TXS
1508+572 is absorbed at wavelengths starting at $6443$ Å, which falls in the
middle of the range of red filters.

The optical and infrared data are displayed in two middle panels of Fig. 1.

####  3.3.1 XMM-OM

The Optical Monitor (OM) onboard XMM-Newton observed TXS 1508+572 simultaneous
to the X-ray band, and took images in the V, U, W1, W2, and M2 filters with a
net exposure of 12 ks, 12 ks, 16 ks, 20 ks, and 16 ks, respectively. We
extract the photometric information using omichain, which is part of the XMM-
Newton SAS. The source is only detected in the V, U, and W1 band. However, all
bands are affected by Lyman-$\alpha$ absorption. We include these photometric
measurements in the broadband SED for display purposes.

####  3.3.2 Zwicky Transient Facility

We obtained optical data from the Zwicky Transient Facility (ZTF; Bellm et
al., 2019; Masci et al., 2019; IRSA, 2022) survey through their public data
release (DR18). The data consist of photometric measurements in the gri filter
system. The data taken in the g-band are fully affected by absorption in the
Lyman-$\alpha$ forest, and are not taken into account in this work. The ZTF
survey data for TXS 1508+572 shown in this work cover a time range from March
2018 through February 2023, with gaps in between epochs of dense monitoring.

####  3.3.3 NEOWISE

At near-infrared wavelength, we use public data from the NEOWISE mission
(Mainzer et al., 2014; NEOWISE Team, 2020), in particular the 2024 data
release, which includes data from 25 December 2013 up to 6 June 2023. Data are
available for two wavelength bands, which are 3.4 and 4.6 $\mu\mathrm{m}$, and
we have initially selected all data available within 5 arcseconds around the
source coordinates from the online
database111https://irsa.ipac.caltech.edu/cgi-bin/Gator/nph-dd. For the data
selection, we have applied the following criteria (see description in
https://academic.oup.com/mnras/article/494/1/764/5815092, which follows
Rakshit et al. (2019)):

  1. 1.

The fit quality given as the reduced $\chi^{2}$ (w1rchi2 / w2rchi2) is less
than 5 in both photometric bands.

  2. 2.

The number of components used to perform a profile fit (nb) of the point
spread function is $<3$.

  3. 3.

The single-exposure images exhibit the best quality (qi_fact $=1$) and are
unaffected by known artifacts (cc_flags =’0000’) and not actively de-blended
(na $=0$).

Additionally, we only take into account frames with the highest quality
(qual_frame $=10$). In our analysis, we only use measured magnitudes, that is,
magnitudes for which an uncertainty is given, but no upper limits. We note
that one W2 band observation taken in December 2019 shows a magnitude of
$\sim$ 12 mag for TXS 1508+572, but for the simultaneous observation in the W1
band, no significant increase is present. Hence, even though data quality is
good and an uncertainty is available for that particular observation, we
exclude it, as it is likely that it is due to a spurious effect and not an
extremely bright and rapid flare of TXS 1508+572.

####  3.3.4 Steward Observatory

TXS 1508+572 was observed on 2023 June 19 UTC with the Steward Observatory
2.3-m Bok Telescope, located on Kitt Peak, AZ, using the SPOL
spectropolarimeter (Schmidt et al., 1992). The faintness of the quasar
dictated that the instrument be configured to provide imaging linear
polarimetry (see, e.g., Smith et al., 2007). A KPNO “nearly Mould” I filter
was used having an effective wavelength of 820 nm and a bandpass of $\sim 185$
nm, FWHM. The filter bandpass selected corresponds to an effective wavelength
of $\sim 154$ nm and FWHM of $\sim 35$ nm in the rest frame of the quasar. The
only major spectral feature expected within the bandpass is C IV $\lambda
1549$. The 9600-s observation yields $q=0.0092\pm 0.0161$ and $u=0.0087\pm
0.0163$ for the normalized linear Stokes parameters using a 6 arcsecond
circular photometric aperture centered on TXS 1508+572. As a result, only a
1-$\sigma$ upper limit of about 2.9% can be estimated for the object’s
polarization, ignoring the statistical bias associated with low-signal-to-
noise linear polarization measurements.

###  3.4 Effelsberg 100-m radio telescope

We observed TXS 1508+572 over 11 months, from February 2022 to January 2023,
with the Effelsberg 100-m telescope. These observations covered a wide
wavelength range in four bands, centered around 60 mm, 20 mm, 14 mm and 7 mm.
We performed cross-scans over the position of the point-like source in azimuth
and elevation, increasing the number of repeating scans for higher frequencies
to account for the lower flux densities. The data are then averaged, undergo a
quality-check by a semi-automatic pipeline, corrected for pointing offsets,
atmospheric opacity and elevation-dependent gain errors and lastly calibrated
using 3C 286, which is a standard calibration source. The observation and data
reduction process is described in more detail in Eppel et al. (2024). We
followed the same strategy, but we increased the number of scan repetitions to
32 for the highest frequencies and also included 60 mm observations. The radio
light curves in the different bands are shown in the third panel from the
bottom in Fig. 1. For the radio spectra, we find indices of $\alpha=-0.41\pm
0.05$, $\alpha=-0.5\pm 0.21$, and $\alpha=-0.24\pm 0.22$ for the spectra
$S\propto\nu^{\alpha}$ between 60 mm and 20 mm, 20 mm and 14 mm, and 14 mm and
7 mm, respectively.

In addition, for frequencies in the 60 mm and 20 mm band, we were able to
perform polarization measurements during 13 epochs as well. We find an average
polarization and standard deviation of $2.9\pm 1.1$ %, $2.6\pm 0.9$ % for 60
mm and 20 mm, respectively, and plot the changes of the polarization in the
second panel from the bottom in Fig. 1. In addition, we also observe some
rotations in the electric vector polarization angle (EVPA), which are shown in
the bottom panel of Fig. 1.

##  4 Results

###  4.1 Maximum $\gamma$-ray flux

Table 1: Gamma-ray detections in different binnings and during the two brightest time ranges of the long-time flaring state, including the measured fluxes and computed luminosities. Bin | TS | Photon Index | Flux [ph cm-2 s-1] | Luminosity [erg s-1]
---|---|---|---|---
MJD $59610-59617$
$59610-59617$ (7 days) | 55.54 | $2.41\pm 0.21$ | $1.2\pm 0.3\times 10^{-7}$ | $2.6\pm 1.0\times 10^{49}$
$59610-59613$ (3 days) | 37.29 | $1.4\pm 0.4\times 10^{-7}$ | $2.8\pm 1.4\times 10^{49}$
$59611-59612$ (1 day) | 24.99 | $2.6\pm 1.0\times 10^{-7}$ | $5.3\pm 2.6\times 10^{49}$
MJD $59806-59813$
$59806-59813$ (7 days) | 35.89 | $2.78\pm 0.29$ | $1.2\pm 0.3\times 10^{-7}$ | $3.2\pm 1.0\times 10^{49}$
$59809-59812$ (3 days) | 19.16 | $1.4\pm 0.5\times 10^{-7}$ | $3.6\pm 1.5\times 10^{49}$
$59811-59812$ (1 day) | 20.97 | $2.6\pm 1.0\times 10^{-7}$ | $6.8\pm 3.0\times 10^{49}$

TXS 1508+572 showed prolonged $\gamma$-ray activity for over six months in
2022. During this period, the source exhibited flux variability on weekly time
scales (Fig. 1). In order to derive the maximum $\gamma$-ray flux and obtain a
value for the overall isotropic $\gamma$-ray luminosity this blazar was able
to produce, we compute light curves with shorter binnings that cover the two
brightest bins in the 7-day binned light curve, that is, from MJD 59610 to
59616 (31 Jan 2023 – 6 Feb 2023), and MJD 59806 to 59812 (15 Aug 2023 – 21 Aug
2023). In order to compute the luminosity, we use the spectral parameters
obtained by modeling each 7-day time span. The results are given in Table 1.
With values ranging from $(2.6-6.8)\times 10^{49}$ erg s-1, TXS 1508+572
exhibits an integrated isotropic $\gamma$-ray luminosity comparable to the
most luminous blazar flares reported by the Fermi-LAT mission since 2008: 3C
279 ($\sim 10^{49}$ erg s-1; Ackermann et al., 2016), CTA 102 ($3\times
10^{50}$ erg s-1; Gasparyan et al., 2018), 3C 454.3 ($1-4\times 10^{49}$ erg
s-1; Nalewajko, 2013, 2017), and PKS 0402-362 ($1.5\times 10^{49}$ erg s-1;
Nalewajko, 2013).

###  4.2 Spectral analysis

####  4.2.1 Fermi-LAT spectra

Table 2: Gamma-ray spectral parameters for the quiescent and flaring state for an energy range of 0.1–100 GeV, which are also used in the broadband SED. The time ranges are from 01 Jan 2010 to (incl.) 31 Dec 2019 (MJD 55197–59214) for the quiescent and 31 Jan to (incl.) 13 Feb 2022 (MJD 59610–59623) for the flaring state. Parameter | Quiescent | Flare
---|---|---
TS | 79.0 | 69.4
Flux [$10^{-9}$ ph cm-2 s-1] | $8.5\pm 0.1$ | $103\pm 22$
Photon Index | $2.99\pm 0.13$ | $2.36\pm 0.17$
Energy flux | $2.7\pm 0.4$ | $59\pm 14$
[$10^{-12}$ erg cm-2 s-1]
Lγ [$10^{48}$erg s-1] | $2.8\pm 0.2$ | $21\pm 7$

The $\gamma$-ray spectra for both the quiescent and flaring states were
modeled with a power law in the energy range from 100 MeV to 300 GeV. The best
fit results for both spectra is listed in Table 2. While for the quiescent
state an integration time of 10 years was necessary to constrain the spectral
parameters well, we were able to produce a spectrum with a similarly bound
photon index over an integration time of only 14 days during the flaring
state. The spectrum during the flare ($\Gamma\approx 2.4$) is significantly
harder than the long-term quiescent state ($\Gamma\approx 3$), and the flux is
a factor of 12 larger. Flaring blazars commonly show spectral hardening
compared to their time-averaged states (e.g., Gasparyan et al., 2018). Other
high-$z$ blazars have displayed this behavior as well (Li et al., 2018; Paliya
et al., 2019).

####  4.2.2 X-ray spectra

Figure 2: Simultaneous X-ray spectra taken in 2002 (top) and 2022 (bottom)
with XMM-Newton, and in 2017 with Swift-XRT and NuSTAR (middle), including
ratios of the spectral bins versus the best-fit model.

All available X-ray spectra from archival and dedicated observations are fit
as described in Section 3.2. Spectra taken with different detectors onboard
XMM-Newton and NuSTAR are fit simultaneously. The best fit results are listed
in Table 3, and the simultaneously fit spectra for 2002, 2017, and 2022 are
plotted in Fig. 2. The photon statistics obtained in the Swift observations
following the $\gamma$-ray flare are insufficient to assess the presence of
X-ray spectral variability. Comparing the two XMM-Newton observations, one
taken in 2002 and one directly after the detection of the $\gamma$-ray flare,
we find that the the X-ray flux is slightly elevated during the $\gamma$-ray
flare but the spectral shape has not significantly changed. The best-fit
spectral index of the NuSTAR spectrum taken in 2017 ($\Gamma=1.12\pm 0.07$) is
harder than archival and dedicated soft X-ray spectra, indicating a hardening
of the X-ray spectrum beyond 10 keV (53 keV in the source-intrinsic reference
frame).

Table 3: Best-fit results for the X-ray spectra, fitted for each instrument.
For Swift-XRT and XMM-Newton, the flux is given for an energy range of 0.5–10
keV, while the flux measured with NuSTAR is given for the range 3–79 keV.

Instrument | ObsID | Date | MJD | Net exposure | Photon Index | Flux | Fit statistic vs. | Instrument
---|---|---|---|---|---|---|---|---
|  |  |  | [ks] |  | [10-12 erg cm-2 s-1] | expect. C value & variance | constant
Swift/XRT | 00081828001 | 2017-04-30 | 57873.4 | 2.1 | $1.17^{+0.36}_{-0.28}$ | $0.78^{+0.24}_{-0.25}$ | 21.5 / $26.3\pm 5.9$ | -
00015096001 | 2022-03-28 | 59666.4 | 9.9 | $1.55\pm 0.17$ | $0.68^{+0.11}_{-0.10}$ | 84.1 / $96.7\pm 11.4$ | -
00015096002 | 2022-09-01 | 59823.5 | 4.8 | $1.54\pm 0.20$ | $0.78^{+0.15}_{-0.13}$ | 50.49 / $68.5\pm 9.3$ | -
00015096003 | 2023-03-04 | 60007.0 | 2.9 | $1.29\pm 0.15$ | $0.88^{+0.14}_{-0.12}$ | 96.4 / $108.4\pm 12.1$ | -
00015096004 | 2023-03-05 | 60008.3 | 4.9
XMM-Newton | 0111260201 | 2002-05-11 | 52405.6 | 9.3 (pn) | $1.53\pm 0.04$ | $0.553\pm 0.022$ | 829.3 / $893.9\pm 39.0$ | $0.97\pm 0.06$ (MOS2)
12.1 (MOS2)
0910390101 | 2022-02-08 | 59618.1 | 58.7 (pn) | $1.481\pm 0.012$ | $0.786\pm 0.011$ | 2729.1 / $2789.1\pm 73.6$ |  $1.003\pm 0.020$ (MOS1)
80.12 (MOS1) |  $1.010^{+0.020}_{-0.019}$ (MOS2)
81.1 (MOS2) |
NuSTAR | 60201013002 | 2017-04-30 | 57873.2 | 36.9 (FPMA) | $1.12^{+0.08}_{-0.07}$ | $3.396^{+0.356}_{-0.001}$ | 864.3 / $917.4\pm 37.9$ | $0.97\pm 0.08$ (FPMB)
36.8 (FPMB)

###  4.3 Flux variability analysis

We assess the flux variability of TXS 1508+572 in different energy bands by
computing the normalized excess variance (Nandra et al., 1997), defined as

$\sigma^{2}_{\mathrm{RMS}}=\frac{1}{N\mu^{2}}\sum^{N}_{i=1}[(X_{i}-\mu)^{2}-\sigma^{2}_{i}],$
(1)

where $N$ is the number of bins, $X_{i}$ is the flux or count rate in each
bin, $\mu$ is the mean flux or count rate, and $\sigma_{i}$ is the statistical
uncertainty associated with $X_{i}$. Depending on the cadence of observations,
we are able to examine daily to monthly time scales, partially in different
energy bands. Negative values of $\sigma^{2}_{\mathrm{RMS}}$ indicate that no
variability is present in addition to the expected noise.

####  4.3.1 Short-term variability with XMM-Newton

Table 4: Normalized excess variance in the X-rays (1 ks binning) as observed
with XMM-Newton and in several GHz radio bands with the 100-m Effelsberg
telescope. X-ray
---
Energy band |  $\sigma^{2}_{\mathrm{RMS}}$ [$10^{-2}$] | Bins
0.3 – 10 keV | $0.10\pm 0.07$ | 73
0.3 – 2 keV | $0.27\pm 0.12$ | 83
2 – 10 keV | $0.16\pm 0.2$ | 71
Radio
Wavelength |  $\sigma^{2}_{\mathrm{RMS}}$ [$10^{-2}$] | Bins
60 mm | $-0.03\pm 0.03$ | 14
20 mm | $0.06\pm 0.05$ | 16
14 mm | $-0.13\pm 0.21$ | 17
7 mm | $-1.6\pm 1.8$ | 10

The long duration of the XMM-Newton observation ($>80$ ks) enables us to probe
intraday flux variability of TXS 1508+572 at X-ray energies. We extract binned
light curves using the EPIC-pn detector with 100s-, 300s-, and 1000s-binning
for the full band (0.3–10 keV) as well as in the soft (0.3–2 keV) and hard
(2–10 keV) band (see Section 3.2.1 for details and the light curves in Fig. 8
in the Appendix). Parts of these light curves are affected by background
flaring events, which we have filtered out by excluding periods where the
background count rate lies above a certain threshold that depends on the
energy band and binning.

We compute the normalized excess variance $\sigma^{2}_{\mathrm{RMS}}$ for all
binning and energy combinations, and find that the 1000 s-binned light curve
inhibits small source-intrinsic variability, which is strongest in the 0.3 – 2
keV band (see Table 4). In the hard band, the result is consistent with the
noise level within uncertainties.

####  4.3.2 Flux variability in the Optical and Infrared

Figure 3: Optical long-term light curves in the r and i band from data taken
with ZTF: light curves in magnitudes (top), light curves in units of mJy
(middle), and normalized excess variance $\sigma^{2}_{\mathrm{RMS}}$ for each
epoch covering a time range of roughly 10 months (bottom). The dotted line
marks the time of the flare detection at $\gamma$-ray energies. The dashed
line in the bottom panel indicates $\sigma^{2}_{\mathrm{RMS}}=0$.

At optical and infrared wavelengths, we are able to probe variability at both
short- and long-term time scales. Observations were conducted from 2018 until
2023 with ZTF in the optical band (see Fig. 3), and from 2013 until 2023 with
NEOWISE in the infrared band (see Fig. 4). While NEOWISE covers the W1 and W2
every six months for a few days with cadences typically ranging from 90
minutes to a few hours, the monitoring cadence with ZTF varies and depends
more strongly on the filter and each year.

Figure 4: Infrared light curves for the W1 and W2 bands taken with NEOWISE:
light curves showing all data points (top), and the average fluxes (middle)
and normalized excess variance $\sigma^{2}_{\mathrm{RMS}}$ (bottom) for each
epoch from 2014 until 2023. The dotted line marks the time of the flare
detection at $\gamma$-ray energies. The dashed line in the bottom panel
indicates $\sigma^{2}_{\mathrm{RMS}}=0$.

For both the optical and infrared light curves we calculate
$\sigma^{2}_{\mathrm{RMS}}$ for each epoch and plot them together with the
light curves in Fig. 3 and Fig. 4, respectively. We only calculate
$\sigma^{2}_{\mathrm{RMS}}$ for epochs that contain more than 5 bins, which is
not the case for some W2 NEOWISE epochs. The flux variability of TXS 1508+572
changes over time in both the infrared and optical band. An infrared flare is
observed to begin $\sim 5$ weeks before the $\gamma$-ray flare is detected,
peaks after the flare is over, and gradually decreases back to baseline (Fig
4). The short-term infrared variability peaks at the onset of the infrared
flare and gradually declines back to the baseline levels. A similar behavior
is observed in the optical data (Fig. 3).

####  4.3.3 Long-term variability in the radio band

The radio data taken with the 100-m Effelsberg telescope (light curve shown in
Fig. 1) have an average monitoring cadence of three to five weeks, depending
on observing frequency. For none of the radio bands do we find
$\sigma^{2}_{\mathrm{RMS}}$ that is significantly above the expected noise
level (see Table 4).

Interestingly, the polarization at 60 mm is more variable, but also higher
than at 20 mm for the time range of at least 3 months before mid-2022. After
that, the polarization in both bands is roughly the same apart from a dip in
the 60 mm band at MJD 59780. Similarly, the EVPA at 60mm is rotating
significantly more than the EVPA measured at 20mm.

###  4.4 SED modeling

Figure 5: Broadband SED of TXS 1508+572 in the quiescent state (blue) and
during the 2022 flare (red). Archival data are taken from the online SED
Builder (https://tools.ssdc.asi.it/SED/) provided by the Space Science Data
Center, and shown via gray data points. Optical data from SARA taken during
the quiescent state are taken from Marcotulli et al. (2020). The dashed
vertical line indicates the frequency of the redshifted Ly$\alpha$ line. The
model fits are produced with the code of Böttcher et al. (2013), using the
parameters listed in Table 5.

The broadband SEDs in the quiescent and flaring states are well described by
the steady-state leptonic model of Böttcher et al. (2013). It is a one-zone
model, assuming a spherical emission region with radius $R_{b}$, located at a
distance $d$ from the central supermassive black hole, moving along the jet
with bulk Lorentz factor $\Gamma$, resulting in relativistic Doppler boosting
of electromagnetic radiation by a Doppler factor
$\delta=\left(\Gamma[1-\beta_{\Gamma}\,\cos\theta_{\rm obs}]\right)^{-1}$,
where $\beta_{\Gamma}$ is the normalized velocity corresponding to the Lorentz
factor $\Gamma$, and $\theta_{\rm obs}$ is the angle between the jet axis and
our line of sight, chosen to be $\theta_{\rm obs}=1/\Gamma$ in order to reduce
the number of free parameters. The code evaluates a self-consistent
equilibrium electron distribution, based on a rapid acceleration process,
injecting a power-law distribution of electrons with index $q$ between a
minimum and maximum electron Lorentz factor $\gamma_{\rm min/max}$. This
injection is self-consistently balanced with radiative energy losses and
escape from the emission region on an escape time scale $t_{\rm esc}=\eta_{\rm
esc}\,R/c$, parameterized by an escape timescale parameter $\eta_{\rm esc}\geq
1$. The resulting electron population with density spectrum $n_{e}(\gamma)$ in
the co-moving frame of the emission region carries a power of
$L_{e}=\pi\,R_{b}^{2}\,\Gamma^{2}\,\beta_{\Gamma}\,m_{e}c^{3}\,\int
d\gamma\gamma n_{e}(\gamma)$ along the jet. Radiation mechanisms included are
synchrotron radiation in a tangled magnetic field of strength $B$,
synchrotron-self-Compton (SSC), and external Compton scattering of both the
direct accretion-disk radiation, EC (disk), and an external blackbody
radiation field (temperature $T_{\rm DT}$ and radiation energy density $u_{\rm
DT}$), assumed to be isotropic in the AGN rest frame, representative of
infrared emission from warm dust, EC (DT). For each model simulation, the code
evaluates, in addition to the electron power mentioned above, the power
carried in the magnetic field (Poynting flux),
$L_{B}=\pi\,R_{b}^{2}\,\Gamma^{2}\,\beta_{\Gamma}\,c(B^{2}/8\pi)$, and the
energy partition ratio, $\epsilon_{Be}=L_{B}/L_{e}$. A set of model parameters
describing the quiescent and flaring SEDs shown in Fig. 5 are listed in Table
5. The chosen solutions are not necessarily unique given the significant
degeneracies in the model.

Fig. 5 illustrates that the SEDs in both states can be well represented with
this model, with the high-energy (X-ray through $\gamma$-ray) emission having
significant contributions from both EC (disk) and EC (DT), while SSC is sub-
dominant. The change from the quiescent to the flaring state is achieved
primarily through a harder injection spectral index of the electron spectrum
and a higher $\gamma_{\rm min}$, indicating an increased electron acceleration
efficiency. The optical emission, in our model fits, is dominated by the
thermal accretion-disk radiation in the quiescent state, in agreement with the
very low degree of polarization measured by the Steward Observatory (also in
the quiescent state, although not contemporaneous with our high-energy
observations). In the flaring state, the much harder electron-synchrotron
spectrum dominates the optical emission. One would therefore expect
significant optical polarization in this state, if the magnetic field is at
least partially ordered in the emission region. Future polarimetric
measurements during flaring states should be able to test this hypothesis.

We note that the jet plasma, especially in the flaring state, is out of
equipartition, dominated by Poynting flux, with the chosen model parameters.
The magnetic field values are also higher by a factor of a few compared to the
estimates of Benke et al. (2024) for the 43 GHz core (0.8 and 1.7 G, depending
on equipartition assumptions). This is consistent with the expectation that
the high-energy emission region is located closer to the central engine than
the 43 GHz core, as our model configuration for the high-energy emission
region is still optically thick at radio frequencies. This is, however,
inconsistent with the estimated distance of $0.32\pm 0.02$ pc for the 43 GHz
core. As there are significant degeneracies in our model parameters, a
similarly adequate fit with a smaller value of $d$ could remedy this
discrepancy (e.g., a smaller distance accompanied by an increased magnetic
field). On the other hand, the estimate of the radio-core distance relied on
the core-shift measurement under the assumption of a conical jet profile, and
on the uncertain viewing angle. A non-conical jet profile and/or different
viewing angle would obviously yield different values for the radio-core
distance.

Table 5: SED model fit parameters used for the modeling shown in Fig. 5. Parameter | Quiescent | Flare
---|---|---
$L_{e}$ [erg s-1] | $5.6\times 10^{45}$ | $2.7\times 10^{45}$
$\gamma_{\rm min}$ | 150 | 450
$\gamma_{\rm max}$ | $2.0\times 10^{5}$ | $3.0\times 10^{4}$
$q$ | 3.1 | 2.3
$B$ [G] | 5.0 | 2.8
$\eta_{\rm esc}$ | 10 | 10
$d$ [pc] | 0.4 | 0.35
$\Gamma$ | 20 | 20
$L_{\rm disk}$ [erg s-1] | $4.0\times 10^{47}$ | $2.5\times 10^{47}$
$R_{b}$ [cm] | $1.3\times 10^{16}$ | $4.0\times 10^{16}$
$M_{\rm BH}$ [M⊙] | $1.5\times 10^{10}$ | $1.5\times 10^{10}$
$T_{\rm DT}$ [K] | $1.0\times 10^{3}$ | $1.0\times 10^{3}$
$u_{\rm DT}$ [erg cm-3] | $5.0\times 10^{-4}$ | $5.0\times 10^{-4}$
$L_{B}$ [erg s-1] | $6.3\times 10^{45}$ | $1.9\times 10^{46}$
$\epsilon_{Be}$ | 1.1 | 7.0

##  5 Discussion

###  5.1 Gamma-ray luminosity of the flare

During the quiescent state, the $\gamma$-ray luminosity of TXS 1508+572 is
$2.8\times 10^{48}$ erg s-1, comparable to other $z>3$ Fermi-LAT detected
blazars Ackermann et al. (2017). During the flaring state, which lasted for
longer than six months in 2022 with two peak times in February and August,
respectively, the $\gamma$-ray luminosity increases to $>5\times 10^{49}$ erg
s-1, making this flare by TXS 1508+572 one of the most luminous that Fermi-LAT
has ever detected. Only two blazars, B3 1343+451 and CTA 102, exhibited flares
with an isotropic $\gamma$-ray luminosity $>10^{50}$ erg s-1 (Sahakyan et al.,
2020; Gasparyan et al., 2018). A few others have shown peak flare luminosities
$>10^{49}$ erg s-1: 3C 454.3 (Nalewajko, 2013, 2017), 3C 279 (Ackermann et
al., 2016), PKS 0402$-$362 (Nalewajko, 2017), and PKS 0537$-$286 (Sahakyan et
al., 2020). Among those, B3 1343+451 and PKS 0537-286 are also high-$z$
objects with $z=2.53$ and $z=3.01$, respectively, showing that blazars in the
early Universe are able to produce equally luminous flares; and given that
three out of seven of the most extreme flares have been observed for these
distant sources hints at those luminous flares being more common in high-$z$
blazars.

###  5.2 Comparison of different broadband SED modeling

The quiescent SED of TXS 1508+572 has been modeled by Ackermann et al. (2017)
using non-simultaneous data and by Marcotulli et al. (2020) using a data set
with simultaneous optical and X-ray data. All models use a one-zone leptonic
model considering both SSC and EC radiation for the high-energy component, and
in all cases EC dominates the X-ray and $\gamma$-ray emission. The SEDs in
Ackermann et al. (2017) and Marcotulli et al. (2020) include thermal emission
from the dust torus even though it is not directly constrained by
observational data in the infrared band. Our model does include potential
radiation from the dust torus as a seed photon field for the EC component, but
we find that EC from the accretion disk strongly dominates the high-energy
emission.

Our SED modeling approach differs from that of Ackermann et al. (2017) and
Marcotulli et al. (2020) in that we attempt to simultaneously describe the
quiescent and flaring states of TXS 1508+572 while changing only a small
subset of model parameters. Spectroscopic estimates of the black hole mass of
TXS 1508+572 are in the range $M_{\mathrm{BH}}=3$–$8\times 10^{8}$ M⊙
(Ackermann et al., 2017; Diana et al., 2022). A standard Shakura-Sunyaev
accretion disk (Shakura & Sunyaev, 1973) corresponding to the estimated
$M_{\mathrm{BH}}$ values is too hot to provide the low frequency seed photons
that our model favors for the origin of the high-energy SED component via the
EC (disk) process. Given that the temperature-dependent peak frequency of the
disk emission is $\nu_{\mathrm{peak}}\propto M^{-1/4}$, our model assuming a
Shakura-Sunyaev disk favors a larger black hole mass of $M=1.5\times
10^{10}\,M_{\odot}$. This discrepancy is in line with quasar spectra tending
to show ionization states compatible with disk temperatures lower than
predicted in the standard Shakura-Sunyaev theory (Bonning et al., 2013). In
addition, the spin of the black hole will also have an impact on the geometry
and temperature profile of the disk, which is not taken into account in our
SED model. In particular, a retrograde spin would push the innermost stable
circular orbit to larger distances compared to a Schwarzschild black hole
(Bardeen et al., 1972), with the effect of also lowering the effective disk
temperature.

The magnetic field values favored by our SED model are also higher than those
in Ackermann et al. (2017) and Marcotulli et al. (2020), and we find a higher
bulk Lorentz factor ($\Gamma=20$ instead of $\Gamma=11$) as well as higher
$\gamma_{\mathrm{min}}$ and $\gamma_{\mathrm{max}}$ values that are expected
if the electron acceleration efficiency increases during a gamma-ray flare.
The measurement of superluminal speeds $>15c$ with our VLBI monitoring
campaign (Benke et al., 2024) supports our higher value of $\Gamma$. The
distance of the dissipation region from the central engine is further away
than described by Ackermann et al. (2017), but at a similar distance as
determined by Marcotulli et al. (2020).

###  5.3 Multiwavelength variability

By computing the $\sigma^{2}_{\mathrm{RMS}}$ value for each energy band for
which we have monitoring data available, albeit covering different time
scales, we can compare the variability of the source at different wavelengths.
While at X-ray energies we see only a small amount of variability within one
day, the IR $\sigma^{2}_{\mathrm{RMS}}$ values for some epochs indicate more
significant variability for time ranges that cover a few days at most. For
blazars in general, the variability of IR emission seems to be correlated with
the presence of $\gamma$-ray emission (Mao et al., 2018). The variability of
TXS 1508+572 in the optical and IR regime seems to precede the $\gamma$-ray
activity given that already before the detection of the $\gamma$-ray flare the
IR flux and the intraday variability had increased. Considering that TXS
1508+572 would be classified as a low-synchrotron-peaked (LSP) source (in the
source-frame), the variability displayed at IR wavelengths would occur in the
R band for a source at $z=0$. Paliya et al. (2017) found that LSP blazars are
more variable on intra-night time scales, and for which the commonly accepted
scenario is that this variability originates from the non-thermal emission of
the jet. Therefore, we suggest that the observed variability is connected to
synchrotron emission from the jet rather than signature from the accretion
disk.

Figure 6: Change of the W1$-$W2 color with respect to the brightness of the W1
band.

We look at possible correlations between the brightness vs. a color change,
and find none at optical wavelengths. This behavior is not unusual though, as
the majority of FSRQs exhibit neither a redder-when-brighter or bluer-when-
brighter trend (Negi et al., 2022). In order to assess the change of the
infrared W1$-$W2 color with the infrared brightness of TXS 1508+572, we plot
the W1$-$W2 color over the W1 magnitude for all NEOWISE data taken
simultaneously in Fig. 6. The infrared flux from TXS 1508+572 shows a bluer-
when-brighter behavior at low fluxes and transitions to a state with no color
change for $m_{\mathrm{W1}}<15$, which could indicate the transition between a
disk-dominated and synchrotron-dominated IR flux. Indeed, the bluer-when-
brighter trend of IR emission has been observed with NEOWISE for the majority
of both FSRQs and BL Lacs (Anjum et al., 2020).

###  5.4 Radio VLBI campaign and polarization

Benke et al. (2024) present the results from a VLBI monitoring campaign
triggered by the $\gamma$-ray flare of TXS 1508+572 in 2022. The Very Long
Baseline Array (VLBA) and the Effelsberg 100-m radio telescope were used to
obtain milliarcsecond-scale images at 15 GHz, 22 GHz, and 43 GHz. In addition,
one observation at 86 GHz was performed with the VLBA and the Green Bank
Telescope, allowing to probe the jet of TXS 1508+572 at $456$ GHz in the rest-
frame frequency of the source. In general, the radio observations support the
underlying scenario in which the injection of fresh electrons in the
acceleration zone explains the observed broadband variability. By tracking the
evolution of the jet over a time range of roughly 10 months, Benke et al.
(2024) find morphological changes in the jet and apparent superluminal speeds
of the jet component motion of $\sim 14c$ – $32c$, depending on frequency. The
latter is in agreement with the bulk Lorentz factor of $\Gamma=20$ that is
favored by our SED model. The observed jet components can be traced back to an
ejection time between 2016 and 2019, and would not be causally connected to
the $\gamma$-ray activity seen in 2022. When computing a 30-day binned light
curve for TXS 1508+572 over the entire duration of the Fermi mission (see Fig.
7), we find several times when the signal from the blazar is TS $\geq 9$
($\gtrsim 3\sigma$), with four bins falling into the year 2020, which
indicates enhanced activity of the source during that time range. Assuming
that the $\gamma$-ray flare in the first half of 2022 has produced a new
outflowing component, we expect to detect a new radio component to be
resolvable by VLBI between 2025 and 2027.

Figure 7: Long-term Fermi-LAT light curve of TXS 1508+572 in 30-day binning.
Bins with TS$<9$ are plotted as arrows at the $2\sigma$ upper limit. The
dashed horizontal line marks the longterm average flux, while the dotted
vertical line marks the flare detection in 2022.

While the VLBI radio observations spatially resolve the jet within inner
several parsec of the source, the radio core itself likely consists of
multiple emission zones that cannot be resolved. Radio polarization
measurements at different frequencies from the Effelsberg observatory allow us
to look into the structure of the unresolved radio core. The Effelsberg data
(Fig. 1) shows that the time variability of the polarization degree is more
pronounced than the Stokes I variability, especially at 60 mm, indicating that
the core itself consists of several emission regions.

VLBI polarization measurements of TXS 1508+572 have been analyzed by
O’Sullivan et al. (2011b), who reported the detection of polarization at 60 mm
and 36 mm for the core, and at 60 mm for a single jet component.
Unfortunately, our VLBI observations do not fall into a time when the 60 mm
net polarization is higher than the net polarization at 20 mm, and without
synchronous data, interpretation of the apparent short timescale variability
of the polarization parameters from MJD 59680 to MJD 59750 (see Fig. 1) would
be too speculative.

##  6 Conclusions

The detection of a $\gamma$-ray flare by a $z>3$ blazar is rare: PKS
0537$-$286 is the only source which has shown multiple flares, which were so
bright that they could be detected on daily time scales (Sahakyan et al.,
2023), and which has been communicated in real time by the LAT collaboration
(Cheung, 2017; Angioni, 2020; Valverde & Forman, 2022). The study by Kreter et
al. (2020) finds nine TS $\geq 25$ detections for three $z>3$ blazars
(excluding PKS 0537$-$286) over a time span of 10 years and 8 months,
resulting in a high-$z$ blazar flare every 14 months on average.

The 2022 flare from TXS 1508+572 reported in this paper is among the most
luminous events that have been seen from this source class. For both the
quiescent and the flaring state, the broadband SED model for TXS 1508+572
requires a dominant contribution from EC emission to describe the high-energy
emission; similar to FSRQs in the local Universe displaying Compton dominance
at all activity states (e.g., Krauß et al., 2016). Our SED modeling suggests a
lower accretion disk temperature than expected in a Shakura-Sunyev disk, as
well as a high bulk Lorentz factor of $\Gamma=20$ that is in line with the
superluminal motion in the jet of TXS 1508+572 described in Benke et al.
(2024). At present, the data covering most of the synchrotron component are
sparse, and a significant portion of optical information that could help
constrain the emission from the accretion disk is lost due Lyman-$\alpha$
absorption.

Hence, constraining the low-energy component of the SED of TXS 1508+572, or
for any blazar with $z>3$ for that matter, presents a challenge. The low
observed optical polarization during a quiescent state allows us to conclude
that the optical emission from TXS 1508+572 is strongly dominated by the
thermal radiation from the accretion disk. However, we were not able to obtain
an optical polarization measurement during the flaring state. Similar to local
FSRQs, TXS 1508+572 also presents variability from the jet emission, which is
measured in the IR band. While only a small number of such systems have been
studied so far, current hard X-ray band and future MeV missions, such as the
Compton Spectrometer and Imager (COSI; Tomsick et al., 2019, 2023), or the
Advanced Particle-astrophysics Telescope (APT; Buckley et al., 2019, 2022),
are ideal to search for these distant objects, and their enhanced sensitivity
promises more detections of powerful sources in the early Universe in coming
years.

##  7 Acknowledgments

The authors thank Lea Marcotulli for providing SARA optical photometry
measurements that are used for modeling the quiescent state SED, and David J.
Thompson for editorial suggestions that improved the manuscript. JH, FE, MK,
and FR acknowledge support from the Deutsche Forschungs- gemeinschaft (DFG,
grants 447572188, 434448349, 465409577). This publication is part of the
M2FINDERS project which has received funding from the European Research
Council (ERC) under the European Union’s Horizon 2020 Research and Innovation
Programme (grant agreement No 101018682). This research has made use of the
XRT Data Analysis Software (XRTDAS) developed under the responsibility of the
ASI Science Center (ASDC), Italy. Furthermore, it made use of a collection of
ISIS functions (ISISscripts) provided by ECAP/Remeis observatory and MIT
(https://www.sternwarte.uni-erlangen.de/isis/). Part of this work is based on
observations obtained with XMM-Newton, an ESA science mission with instruments
and contributions directly funded by ESA Member States and NASA. This
publication makes use of data products from the Near-Earth Object Wide-field
Infrared Survey Explorer (NEOWISE), which is a joint project of the Jet
Propulsion Laboratory/California Institute of Technology and the University of
Arizona. NEOWISE is funded by the National Aeronautics and Space
Administration. This research is also partly based on observations obtained
with the Samuel Oschin Telescope 48-inch and the 60-inch Telescope at the
Palomar Observatory as part of the Zwicky Transient Facility project. ZTF is
supported by the National Science Foundation under Grant No. AST-2034437 and a
collaboration including Caltech, IPAC, the Weizmann Institute for Science, the
Oskar Klein Center at Stockholm University, the University of Maryland,
Deutsches Elektronen-Synchrotron and Humboldt University, the TANGO Consortium
of Taiwan, the University of Wisconsin at Milwaukee, Trinity College Dublin,
Lawrence Livermore National Laboratories, and IN2P3, France. Operations are
conducted by COO, IPAC, and UW. The ZTF forced-photometry service was funded
under the Heising-Simons Foundation grant #12540303 (PI: Graham). In parts,
this work is based on observations with the 100-m telescope of the MPIfR (Max-
Planck-Institut für Radioastronomie) at Effelsberg.

##  Appendix A Additional information about X-ray data

The extraction regions of both the newly obtained and the archival XMM-Newton
observation are listed in Table 6. The data extraction procedure is described
in the main text in Section 3.2.1.

Table 6: Extraction regions for XMM-Newton spectra. Source and background coordinates are given in degrees, and the radii are given in arcseconds. ObsId | Detector | Source region | Background region
---|---|---|---
R.A. | Decl. | Radius | R.A. | Decl. | Radius
0111260201 | pn | 227.5125 | 57.0455 | 35 | 227.2732 | 57.0657 | 60
MOS2 | 227.5113 | 57.0455 | 35 | 227.5839 | 56.9906 | 100
0910390101 | pn | 227.5121 | 57.0459 | 35 | 227.5514 | 57.0746 | 60
MOS1 | 227.5135 | 57.0455 | 27.5 | 227.3028 | 57.0376 | 60
MOS2 | 227.5113 | 57.0455 | 27.5 | 227.2892 | 57.0631 | 60

The XMM-Newton light curves in 1000 s-binning and for different energy ranges
are shown in Fig. 8.

Figure 8: XMM light curve in 1 ks binning and for different energy bands.
While we show the source and background bins for the entire light curve, we
only use the subtracted light curve for computing $\sigma^{2}_{\mathrm{rms}}$
values, as some of the bins are affecting by background flaring at the
beginning of the light curve as well as towards the end.

## References

  * Abdollahi et al. (2020) Abdollahi, S., Acero, F., Ackermann, M., et al. 2020, ApJS, 247, 33, doi: 10.3847/1538-4365/ab6bcb
  * Ackermann et al. (2016) Ackermann, M., Anantua, R., Asano, K., et al. 2016, ApJ, 824, L20, doi: 10.3847/2041-8205/824/2/L20
  * Ackermann et al. (2017) Ackermann, M., Ajello, M., Baldini, L., et al. 2017, ApJ, 837, L5, doi: 10.3847/2041-8213/aa5fff
  * Aharonian (2000) Aharonian, F. A. 2000, New A, 5, 377, doi: 10.1016/S1384-1076(00)00039-7
  * Aharonian (2002) —. 2002, MNRAS, 332, 215, doi: 10.1046/j.1365-8711.2002.05292.x
  * Ajello et al. (2020) Ajello, M., Angioni, R., Axelsson, M., et al. 2020, ApJ, 892, 105, doi: 10.3847/1538-4357/ab791e
  * Ajello et al. (2022) Ajello, M., Baldini, L., Ballet, J., et al. 2022, ApJS, 263, 24, doi: 10.3847/1538-4365/ac9523
  * Alexander & Natarajan (2014) Alexander, T., & Natarajan, P. 2014, Science, 345, 1330, doi: 10.1126/science.1251053
  * Angioni (2020) Angioni, R. 2020, The Astronomer’s Telegram, 14285
  * Anjum et al. (2020) Anjum, A., Stalin, C. S., Rakshit, S., Gudennavar, S. B., & Durgapal, A. 2020, MNRAS, 494, 764, doi: 10.1093/mnras/staa771
  * Atwood et al. (2009) Atwood, W. B., Abdo, A. A., Ackermann, M., et al. 2009, ApJ, 697, 1071, doi: 10.1088/0004-637X/697/2/1071
  * Bardeen et al. (1972) Bardeen, J. M., Press, W. H., & Teukolsky, S. A. 1972, ApJ, 178, 347, doi: 10.1086/151796
  * Begelman et al. (2006) Begelman, M. C., Volonteri, M., & Rees, M. J. 2006, 370, 289, doi: 10.1111/j.1365-2966.2006.10467.x
  * Belladitta et al. (2020) Belladitta, S., Moretti, A., Caccianiga, A., et al. 2020, A&A, 635, L7, doi: 10.1051/0004-6361/201937395
  * Belladitta et al. (2022) Belladitta, S., Caccianiga, A., Diana, A., et al. 2022, 660, A74, doi: 10.1051/0004-6361/202142335
  * Bellm et al. (2019) Bellm, E. C., Kulkarni, S. R., Graham, M. J., et al. 2019, PASP, 131, 018002, doi: 10.1088/1538-3873/aaecbe
  * Benke et al. (2024) Benke, P., Gokus, A., Lisakov, M., et al. 2024, arXiv e-prints, arXiv:2406.03135, doi: 10.48550/arXiv.2406.03135
  * Blandford & Levinson (1995) Blandford, R. D., & Levinson, A. 1995, ApJ, 441, 79, doi: 10.1086/175338
  * Błażejowski et al. (2000) Błażejowski, M., Sikora, M., Moderski, R., & Madejski, G. M. 2000, ApJ, 545, 107, doi: 10.1086/317791
  * Bloemen et al. (1995) Bloemen, H., Bennett, K., Blom, J. J., et al. 1995, A&A, 293, L1
  * Bloom & Marscher (1996) Bloom, S. D., & Marscher, A. P. 1996, ApJ, 461, 657, doi: 10.1086/177092
  * Boettcher & Schlickeiser (1997) Boettcher, M., & Schlickeiser, R. 1997, A&A, 325, 866, doi: 10.48550/arXiv.astro-ph/9703069
  * Bonning et al. (2013) Bonning, E. W., Shields, G. A., Stevens, A. C., & Salviander, S. 2013, ApJ, 770, 30, doi: 10.1088/0004-637X/770/1/30
  * Böttcher et al. (2013) Böttcher, M., Reimer, A., Sweeney, K., & Prakash, A. 2013, ApJ, 768, 54, doi: 10.1088/0004-637X/768/1/54
  * Buckley et al. (2019) Buckley, J., Bergstrom, L., Binns, B., et al. 2019, in Bulletin of the American Astronomical Society, Vol. 51, 78
  * Buckley et al. (2022) Buckley, J., Adapt, Alnussirat, S., et al. 2022, in 37th International Cosmic Ray Conference, 655, doi: 10.22323/1.395.0655
  * Burke et al. (2024) Burke, C. J., Liu, X., & Shen, Y. 2024, MNRAS, 527, 5356, doi: 10.1093/mnras/stad3592
  * Cash (1979) Cash, W. 1979, ApJ, 228, 939, doi: 10.1086/156922
  * Cheung (2004) Cheung, C. C. 2004, ApJ, 600, L23, doi: 10.1086/381366
  * Cheung (2017) —. 2017, The Astronomer’s Telegram, 10356
  * Dermer & Schlickeiser (1993) Dermer, C. D., & Schlickeiser, R. 1993, ApJ, 416, 458, doi: 10.1086/173251
  * Diana et al. (2022) Diana, A., Caccianiga, A., Ighina, L., et al. 2022, MNRAS, 511, 5436, doi: 10.1093/mnras/stac364
  * Eppel et al. (2024) Eppel, F., Kadler, M., Heßdörfer, J., et al. 2024, A&A, 684, A11, doi: 10.1051/0004-6361/202348262
  * Frey et al. (1997) Frey, S., Gurvits, L. I., Kellermann, K. I., Schilizzi, R. T., & Pauliny-Toth, I. I. K. 1997, A&A, 325, 511
  * Gasparyan et al. (2018) Gasparyan, S., Sahakyan, N., Baghmanyan, V., & Zargaryan, D. 2018, ApJ, 863, 114, doi: 10.3847/1538-4357/aad234
  * Gehrels et al. (2004) Gehrels, N., Chincarini, G., Giommi, P., et al. 2004, ApJ, 611, 1005, doi: 10.1086/422091
  * Ghisellini et al. (2010) Ghisellini, G., Della Ceca, R., Volonteri, M., et al. 2010, 405, 387, doi: 10.1111/j.1365-2966.2010.16449.x
  * Gokus et al. (2022) Gokus, A., Kreter, M., Kadler, M., et al. 2022, The Astronomer’s Telegram, 15202
  * Gunn & Peterson (1965) Gunn, J. E., & Peterson, B. A. 1965, ApJ, 142, 1633, doi: 10.1086/148444
  * HI4PI Collaboration et al. (2016) HI4PI Collaboration, Ben Bekhti, N., Flöer, L., et al. 2016, A&A, 594, A116, doi: 10.1051/0004-6361/201629178
  * Hook et al. (1995) Hook, I. M., McMahon, R. G., Patnaik, A. R., et al. 1995, MNRAS, 273, L63, doi: 10.1093/mnras/273.1.L63
  * Houck & Denicola (2000) Houck, J. C., & Denicola, L. A. 2000, in Astronomical Data Analysis Software and Systems IX, ed. N. Manset, C. Veillet, & D. Crabtree, ASP Conf. Ser. No. 216 (San Francisco: Astron. Soc. Pacific), 591
  * Inayoshi et al. (2020) Inayoshi, K., Visbal, E., & Haiman, Z. 2020, ARA&A, 58, 27, doi: 10.1146/annurev-astro-120419-014455
  * IRSA (2022) IRSA. 2022, Zwicky Transient Facility Image Service, IPAC, doi: 10.26131/IRSA539
  * Jansen et al. (2001) Jansen, F., Lumb, D., Altieri, B., et al. 2001, A&A, 365, L1, doi: 10.1051/0004-6361:20000036
  * Johnson & Haardt (2016) Johnson, J. L., & Haardt, F. 2016, PASA, 33, e007, doi: 10.1017/pasa.2016.4
  * Kappes et al. (2022) Kappes, A., Burd, P. R., Kadler, M., et al. 2022, A&A, 663, A44, doi: 10.1051/0004-6361/202141720
  * Krauß et al. (2016) Krauß, F., Wilms, J., Kadler, M., et al. 2016, A&A, 591, A130, doi: 10.1051/0004-6361/201628595
  * Kreter et al. (2020) Kreter, M., Gokus, A., Krauss, F., et al. 2020, ApJ, 903, 128, doi: 10.3847/1538-4357/abb8da
  * Lai et al. (2024) Lai, S., Onken, C. A., Wolf, C., Bian, F., & Fan, X. 2024, arXiv e-prints, arXiv:2405.10721, doi: 10.48550/arXiv.2405.10721
  * Li et al. (2018) Li, S., Xia, Z.-Q., Liang, Y.-F., Liao, N.-H., & Fan, Y.-Z. 2018, ApJ, 853, 159, doi: 10.3847/1538-4357/aaa3fb
  * Liao et al. (2019) Liao, N.-H., Dou, L.-M., Jiang, N., et al. 2019, ApJ, 879, L9, doi: 10.3847/2041-8213/ab2893
  * Liao et al. (2018) Liao, N.-H., Li, S., & Fan, Y.-Z. 2018, ApJ, 865, L17, doi: 10.3847/2041-8213/aae20d
  * Liao et al. (2020) Liao, N.-H., Li, S., Sheng, Z.-F., & Fan, Y.-Z. 2020, ApJ, 898, L56, doi: 10.3847/2041-8213/aba82e
  * Lobanov et al. (2001) Lobanov, A. P., Gurvits, L. I., Frey, S., et al. 2001, ApJ, 547, 714, doi: 10.1086/318391
  * Madau & Rees (2001) Madau, P., & Rees, M. J. 2001, 551, L27, doi: 10.1086/319848
  * Mainzer et al. (2014) Mainzer, A., Bauer, J., Cutri, R. M., et al. 2014, ApJ, 792, 30, doi: 10.1088/0004-637X/792/1/30
  * Mannheim (1993) Mannheim, K. 1993, A&A, 269, 67, doi: 10.48550/arXiv.astro-ph/9302006
  * Mannheim & Biermann (1992) Mannheim, K., & Biermann, P. L. 1992, A&A, 253, L21
  * Mao et al. (2018) Mao, L., Zhang, X., & Yi, T. 2018, Ap&SS, 363, 167, doi: 10.1007/s10509-018-3388-9
  * Maraschi et al. (1992) Maraschi, L., Ghisellini, G., & Celotti, A. 1992, ApJ, 397, L5, doi: 10.1086/186531
  * Marcotulli et al. (2020) Marcotulli, L., Paliya, V., Ajello, M., et al. 2020, ApJ, 889, 164, doi: 10.3847/1538-4357/ab65f5
  * Masci et al. (2019) Masci, F. J., Laher, R. R., Rusholme, B., et al. 2019, PASP, 131, 018003, doi: 10.1088/1538-3873/aae8ac
  * Massaro et al. (2009) Massaro, E., Giommi, P., Leto, C., et al. 2009, A&A, 495, 691, doi: 10.1051/0004-6361:200810161
  * Massaro et al. (2015) Massaro, E., Maselli, A., Leto, C., et al. 2015, Ap&SS, 357, 75, doi: 10.1007/s10509-015-2254-2
  * Mathur & Elvis (1995) Mathur, S., & Elvis, M. 1995, AJ, 110, 1551, doi: 10.1086/117627
  * Mattox et al. (1996) Mattox, J. R., Bertsch, D. L., Chiang, J., et al. 1996, ApJ, 461, 396, doi: 10.1086/177068
  * McKinney et al. (2014) McKinney, J. C., Tchekhovskoy, A., Sadowski, A., & Narayan, R. 2014, MNRAS, 441, 3177, doi: 10.1093/mnras/stu762
  * Medvedev et al. (2020) Medvedev, P., Sazonov, S., Gilfanov, M., et al. 2020, MNRAS, 497, 1842, doi: 10.1093/mnras/staa2051
  * Migliori et al. (2023) Migliori, G., Siemiginowska, A., Sobolewska, M., et al. 2023, MNRAS, 524, 1087, doi: 10.1093/mnras/stad1959
  * Moran & Helfand (1997) Moran, E. C., & Helfand, D. J. 1997, ApJ, 484, L95, doi: 10.1086/310787
  * Mücke & Protheroe (2001) Mücke, A., & Protheroe, R. J. 2001, Astroparticle Physics, 15, 121, doi: 10.1016/S0927-6505(00)00141-9
  * Mücke et al. (2003) Mücke, A., Protheroe, R. J., Engel, R., Rachen, J. P., & Stanev, T. 2003, Astroparticle Physics, 18, 593, doi: 10.1016/S0927-6505(02)00185-8
  * Nalewajko (2013) Nalewajko, K. 2013, MNRAS, 430, 1324, doi: 10.1093/mnras/sts711
  * Nalewajko (2017) —. 2017, Galaxies, 5, 100, doi: 10.3390/galaxies5040100
  * Nandra et al. (1997) Nandra, K., George, I. M., Mushotzky, R. F., Turner, T. J., & Yaqoob, T. 1997, ApJ, 476, 70, doi: 10.1086/303600
  * Negi et al. (2022) Negi, V., Joshi, R., Chand, K., et al. 2022, MNRAS, 510, 1791, doi: 10.1093/mnras/stab3591
  * NEOWISE Team (2020) NEOWISE Team. 2020, NEOWISE-R Single Exposure (L1b) Source Table, IPAC, doi: 10.26131/IRSA144
  * Orienti et al. (2014) Orienti, M., D’Ammando, F., Giroletti, M., et al. 2014, MNRAS, 444, 3040, doi: 10.1093/mnras/stu1644
  * O’Sullivan et al. (2011a) O’Sullivan, S. P., Gabuzda, D. C., & Gurvits, L. I. 2011a, MNRAS, 415, 3049, doi: 10.1111/j.1365-2966.2011.18915.x
  * O’Sullivan et al. (2011b) —. 2011b, MNRAS, 415, 3049, doi: 10.1111/j.1365-2966.2011.18915.x
  * Paliya et al. (2016) Paliya, V. S., Parker, M. L., Fabian, A. C., & Stalin, C. S. 2016, ApJ, 825, 74, doi: 10.3847/0004-637X/825/1/74
  * Paliya et al. (2017) Paliya, V. S., Stalin, C. S., Ajello, M., & Kaur, A. 2017, ApJ, 844, 32, doi: 10.3847/1538-4357/aa77f5
  * Paliya et al. (2019) Paliya, V. S., Ajello, M., Ojha, R., et al. 2019, ApJ, 871, 211, doi: 10.3847/1538-4357/aafa10
  * Planck Collaboration et al. (2016) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016, A&A, 594, A13, doi: 10.1051/0004-6361/201525830
  * Rakshit et al. (2019) Rakshit, S., Johnson, A., Stalin, C. S., Gandhi, P., & Hoenig, S. 2019, MNRAS, 483, 2362, doi: 10.1093/mnras/sty3261
  * Sahakyan et al. (2023) Sahakyan, N., Harutyunyan, G., & Israyelyan, D. 2023, MNRAS, 521, 1013, doi: 10.1093/mnras/stad517
  * Sahakyan et al. (2020) Sahakyan, N., Israyelyan, D., Harutyunyan, G., Khachatryan, M., & Gasparyan, S. 2020, MNRAS, 498, 2594, doi: 10.1093/mnras/staa2477
  * Schmidt et al. (1992) Schmidt, G. D., Stockman, H. S., & Smith, P. S. 1992, ApJ, 398, L57, doi: 10.1086/186576
  * Schneider et al. (2007) Schneider, D. P., Hall, P. B., Richards, G. T., et al. 2007, AJ, 134, 102, doi: 10.1086/518474
  * Shakura & Sunyaev (1973) Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337
  * Siemiginowska et al. (2003) Siemiginowska, A., Smith, R. K., Aldcroft, T. L., et al. 2003, ApJ, 598, L15, doi: 10.1086/380497
  * Sikora et al. (1994) Sikora, M., Begelman, M. C., & Rees, M. J. 1994, ApJ, 421, 153, doi: 10.1086/173633
  * Sądowski et al. (2014) Sądowski, A., Narayan, R., McKinney, J. C., & Tchekhovskoy, A. 2014, MNRAS, 439, 503, doi: 10.1093/mnras/stt2479
  * Smith et al. (2007) Smith, P. S., Williams, G. G., Schmidt, G. D., Diamond-Stanic, A. M., & Means, D. L. 2007, ApJ, 663, 118, doi: 10.1086/517992
  * Strüder et al. (2001) Strüder, L., Briel, U., Dennerl, K., et al. 2001, A&A, 365, L18, doi: 10.1051/0004-6361:20000066
  * Titov et al. (2023) Titov, O., Frey, S., Melnikov, A., et al. 2023, AJ, 165, 69, doi: 10.3847/1538-3881/aca964
  * Tomsick et al. (2019) Tomsick, J., Zoglauer, A., Sleator, C., et al. 2019, in Bulletin of the American Astronomical Society, Vol. 51, 98, doi: 10.48550/arXiv.1908.04334
  * Tomsick et al. (2023) Tomsick, J. A., Boggs, S. E., Zoglauer, A., et al. 2023, arXiv e-prints, arXiv:2308.12362, doi: 10.48550/arXiv.2308.12362
  * Turner et al. (2001) Turner, M. J. L., Abbey, A., Arnaud, M., et al. 2001, A&A, 365, L27, doi: 10.1051/0004-6361:20000087
  * Urry & Padovani (1995) Urry, C. M., & Padovani, P. 1995, PASP, 107, 803, doi: 10.1086/133630
  * Valverde & Forman (2022) Valverde, J., & Forman, J. 2022, The Astronomer’s Telegram, 15405
  * Verner et al. (1996) Verner, D. A., Ferland, G. J., Korista, K. T., & Yakovlev, D. G. 1996, ApJ, 465, 487, doi: 10.1086/177435
  * Wang et al. (2021) Wang, F., Yang, J., Fan, X., et al. 2021, ApJ, 907, L1, doi: 10.3847/2041-8213/abd8c6
  * Wilms et al. (2000) Wilms, J., Allen, A., & McCray, R. 2000, ApJ, 542, 914, doi: 10.1086/317016
  * Wood et al. (2017) Wood, M., Caputo, R., Charles, E., et al. 2017, in International Cosmic Ray Conference, Vol. 301, 35th International Cosmic Ray Conference (ICRC2017), 824, doi: 10.22323/1.301.0824
  * Yuan et al. (2003) Yuan, W., Fabian, A. C., Celotti, A., & Jonker, P. G. 2003, MNRAS, 346, L7, doi: 10.1046/j.1365-2966.2003.07234.x