NuSTAR Perspective on High-redshift MeV Blazars

The hard X-ray band has been crucial to studying some of the most powerful blazars (see e.g., Tavecchio et al. 2000; Massaro et al. 2004a, 2004b, 2006; Donato et al. 2005). More recently, the outstanding sensitivity of the Nuclear Spectroscopic Telescope Array (NuSTAR, 3–79 keV, Harrison et al. 2013) has enabled us to find and study some of the most distant and luminous ones (e.g., Sbarrato et al. 2013; Tagliaferri et al. 2015; Ajello et al. 2016; Paliya et al. 2016; Sbarrato et al. 2016; Marcotulli et al. 2017). Harboring highly relativistic jets pointed closely at the observer (θ_V < 1/Γ, θ_V being the viewing angle and Γ the bulk Lorentz factor, Γ ∼ 10−15, Urry & Padovani 1995), this subclass of the Active Galactic Nuclei (AGNs) is home to some of the most energetic particle acceleration and radiation processes known in astrophysics. The boost in flux ascribed to relativistic beaming, arising from the peculiar orientation of the jets, renders them visible at redshifts well beyond z = 2 (the farthest blazar detected so far is at z = 5.47, Romani et al. 2004), making them extraordinary beacons with which to explore the early universe. Their typical double-hump spectral energy distribution (SED) spans the whole electromagnetic spectrum and is shaped by the nonthermal processes occurring in the jets. Relativistic electrons, spiraling along the magnetic field lines, undergo both synchrotron and inverse Compton (IC) processes. The first produces a peak in the SED located between infrared and X-ray frequencies. The second instead results in a peak located between X- and γ-ray energies. If the electrons interact with a source of low-energy photons within the jet, this is referred to as Synchrotron Self Compton (SSC, e.g., Ghisellini & Maraschi 1989), whereas if the photons are external to the jet (i.e., the accretion disk, the torus, and/or the broad line region, BLR), it is referred to as External Compton process (EC, e.g., Sikora et al. 1994). Based on their optical spectra, blazars are usually classified either as BL Lacertae objects (BL Lacs) or flat spectrum radio quasars (FSRQs), the first showing weak or no emission lines, the second showing broad (EW > 5 Å) ones. Following Abdo et al. (2010a), these sources can also be classified according to the position of the synchrotron peak (${\nu }_{\mathrm{peak}}^{{\rm{S}}}$), with low-, intermediate-, and high-synchrotron peak (LSP, ISP, HSP) blazars having, respectively, ${\nu }_{\mathrm{peak}}^{{\rm{S}}}\lt {10}^{14}$ Hz, ${10}^{14}\lt {\nu }_{\mathrm{peak}}^{{\rm{S}}}\lt {10}^{15}$ Hz, and ${\nu }_{\mathrm{peak}}^{{\rm{S}}}\gt {10}^{15}$ Hz. FSRQs usually belong to the LSP class and at the high-luminosity end of such subclass are the so-called "MeV blazars," whose high-energy peak falls in (or close to) the MeV band. With bolometric luminosities exceeding 10⁴⁸ erg s⁻¹, these are among the most powerful objects in the universe. In fact, they host powerful relativistic jets (Ghisellini et al. 2014), are usually found at high-redshift (z > 2, e.g., Ajello et al. 2009; Ghisellini et al. 2010; Ackermann et al. 2017; Marcotulli et al. 2017), and typically host billion solar mass black holes (e.g., Ghisellini et al. 2010; Paliya et al. 2017a).

In accordance with the blazar sequence (Ghisellini et al. 1998), as both source luminosity and distance increase, the two SED peaks drift toward lower energies and the Compton Dominance (CD, the ratio between IC and synchrotron peak luminosities) increases (CD > 1). As a consequence, MeV blazars are bright in the hard X-rays (above 10 keV) and therefore represent perfect targets for NuSTAR. Moreover, high-energy emission dominates their energetics (see Madejski & Sikora 2016 for a recent review), hence good coverage of the IC hump is required in order to obtain a precise determination of the peak position. This, in turn, allows us to constrain jet parameters such as viewing angle, power, bulk Lorentz factor, as well as the shape of underlying electron distribution (see, e.g., Celotti & Ghisellini 2008; Ghisellini & Tavecchio 2009; Gao et al. 2011; Chen et al. 2015; Finke 2016; Ghisellini et al. 2019; Paliya et al. 2019b), keys to untangle the physics of these cosmic monsters. Undoubtedly, NuSTAR is the best instrument to sample the rising part of their high-energy spectrum, also enabling us to detect possible hard X-ray peculiarities of these sources (e.g., flattening of the spectrum, see Paliya et al. 2016, and variability, e.g., Sbarrato et al. 2016). Lacking an all sky MeV instrument (e.g., AMEGO, Rando 2017), complementary to hard X-ray observations are the γ-ray ones obtained with the Fermi-Large Area Telescope (Fermi-LAT; Atwood et al. 2009). Encompassing a range of energies from 20 MeV to >300 GeV, it allows the sampling of the falling part of the IC hump (e.g., Ackermann et al. 2015) and, together with NuSTAR, it is fundamental to place tight constraint on the jet power.

On another important note, the SED peaks' shifts also reveal the underlying thermal emission from the accretion disk at optical-UV energies (e.g., Ghisellini et al. 2010), which can be modeled to extract the black hole mass and the accretion disk luminosity (Shakura & Sunyaev 1973). This fact renders MeV blazars excellent probes for addressing open issues such as the accretion disk-jet connection in jetted quasars (e.g., Paliya et al. 2017a), as well as the understanding of black hole growth and evolution in the early stages of the universe (Volonteri et al. 2011). Indeed, since the detection of a single source with jets aligned to our line of sight implies the existence of 2Γ² similar sources, at the same z, with same black hole masses but misaligned jets, finding more high-redshift blazars allows us to probe their intrinsic population and constrain the space density of supermassive black holes (e.g., Sbarrato et al. 2014).

Even though blazars are the most numerous extragalactic sources in Fermi-LAT catalogs (see Abdo et al. 2010b; Nolan et al. 2012; Acero et al. 2015), few are detected at z > 3, probably since their peak high-energy emission lies below the Fermi-LAT band. Using the improved Pass 8 data set, taking advantage of the higher detection threshold of the Fermi-LAT, five new γ-ray emitting MeV blazars have been found to lie at z > 3, the farthest located at z = 4.31 (Ackermann et al. 2017). Since, on average, high-redshift sources lack good quality multifrequency observations, we started a multiwavelength campaign in order to characterize these newly found γ-ray emitting blazars. In particular, we wanted to obtain good X-ray spectral coverage, and, importantly, sensitive hard X-ray data, in order to derive their X-ray properties and tightly constrain the jet emission.

In this paper, we present the broadband follow up of six sources: NVSS J013126−100931 (z = 3.51, hereafter J013126−100931), NVSS J020346+113445 (z = 3.63; hereafter J020346+113445), NVSS J064632+445116 (z = 3.39; hereafter J064632+445116), NVSS J135406−020603 (z = 3.71, hereafter J135406−020603), NVSS J151002+570243 (z = 4.31; hereafter J151002+570243), and NVSS J212912−153841 (z = 3.26, hereafter J212912−153841). Four of these objects belong to the newly detected γ-ray blazars, while two are candidate γ-ray sources lying just below the LAT detection threshold. In this work, we publish their first >10 keV detection with NuSTAR.¹⁰ Since blazars are extremely variable sources in all energy ranges (e.g., Abdo et al. 2010c; Stroh & Falcone 2013; Marchesini et al. 2016; Paliya et al. 2017b), to accurately constrain the sources' SEDs and to reliably measure the X-ray spectral parameters, we also obtain quasi-simultaneous observations from a variety of facilities. The soft X-ray and UV-optical bands are covered by the Neil Gehrels Swift Observatory X-ray Telescope (Swift-XRT; Burrows et al. 2005) and Ultraviolet and Optical Telescope (Swift-UVOT; Roming et al. 2005). With the aim of studying possible long-term variability properties of these sources in the soft X-ray band, we check if multiple Swift-XRT, as well as Chandra (Weisskopf et al. 2000) and XMM-Newton (Jansen et al. 2001) observations are available. All sources but J151002+570243 and J135406−020603 have multiple Swift-XRT observations, while Chandra and XMM-Newton data are found for J151002+570243 and J212912−153841. We complement UVOT observations with quasi-simultaneous ones from the Southeastern Association for Research in Astronomy (SARA; Keel et al. 2017) facilities, in order to detect the disk emission peak. The γ-ray spectra of the sources are taken from Ackermann et al. (2017). The data reduction pipelines are given in Section 2. In order to derive the important physical parameters of these objects, we perform the SED modeling (described in Section 3.4). Section 3 details the results of the analysis, which are then discussed in Section 4. Throughout the paper, we use cosmological parameters H₀ = 67.8 km s⁻¹ Mpc⁻¹, Ω_m = 0.308, and Ω_Λ = 0.692 (Planck Collaboration et al. 2016).

2.1. Fermi-LAT

2.2. NuSTAR

NuSTAR observed the six blazars for net exposures of ∼30–50 ks during the dates reported in Table 1. We process the data for both instrument Focal Plane Modules (FPMA and FPMB; Harrison et al. 2013) using the NuSTAR Data Analysis Software (NUSTARDAS) v1.8.0. The event files are calibrated using the task nupipeline, with the response file taken from the latest Calibration Database (CALDB, v.20180419). The generation of source and background spectra, ancillary and response matrix files, is achieved using the nuproducts script. For source regions, we select circles with radii of 20'' centered on the targets; the background events are extracted from circles of 30'' radii located in a source-free region on the same frame. The excellent sensitivity of NuSTAR provides us with a good signal-to-noise detection of all sources. Therefore, if the source has also good soft X-ray photon statistics (see Sections 2.3–2.4), NuSTAR spectra are rebinned to have at least 15 counts per bin.Visual inspection of the light curves shows no detected variability for all six sources.

Table 1.  Table of Simultaneous NuSTAR and Swift-XRT Observations

NuSTAR + Swift-XRT
Name	NHa	z	Exp. time (ks)	Exp. time (ks)	Date
			(NuSTAR)	(Swift-XRT)
J013126−100931	3.40	3.51	30	3.4	2018 Jul 21
J020346+113445	5.96	3.63	33	1.7	2016 Nov 8
J064632+445116	10.7	3.39	32	1.8	2016 Nov 28
J135406−020603	3.35	3.71	32	1.9	2017 May 5
J151002+570243	1.57	4.31	53	2.0	2017 Apr 30
J212912−153841	4.85	3.26	37	2.6	2018 Sep 26

Note.

^aGalactic column density in units of 10²⁰$\,{\mathrm{cm}}^{-2}$ (Kalberla et al. 2005).

Download table as:  ASCIITypeset image

2.3. Swift

2.4. Chandra and XMM-Newton

2.5. SARA

3.1. Soft X-Ray Long-term Variability

3.2. X-Ray Spectral Analysis

We fit the soft X-ray (0.3–15 keV) and NuSTAR FPMA and FPMB (3–79 keV) spectra using XSPEC. We include Galactic absorption (N_H), using the Galactic neutral hydrogen column densities from Kalberla et al. (2005). We use the C-statistics (Cash 1979) for sources with low soft X-ray counts ($\lt 30\,\mathrm{cts}$), otherwise we employ standard χ² statistics. Sources with good X-ray signal are: J013126−100931, J064632+445116, J151002+570243, and J212912−153841. High-redshift blazars have been found to show spectral softening at soft X-rays (Page et al. 2005; Saez et al. 2011; Eitan & Behar 2013), which has been proposed to be linked to either a property of the intrinsic electron population underlying the emission or to the presence of an intervening absorber residing in the intergalactic medium (IGM, see Arcodia et al. 2018), or both. In order to test this scenario, at first these objects are fitted with a simple power-law model along with two curved ones, a broken power law and a log-parabola, all with absorption fixed at the Galactic value. Then, if a curved model is favored over a simple power-law one (highlighting the presence of a softening of the soft X-ray part of the spectrum), we also test the possibility of an intrinsic absorber. Therefore, we further fit these objects with both a power-law and a broken power-law model including a cold absorber (ztbabsin XSPEC). All X-ray spectral fit parameters are provided in Table 3 and the fits and residuals for each fitted model can be found in the Appendix. The main results are listed below:

• 1.  
  Fits' residuals and F-test results show that both J013126−100931 and J212912−153841 favor a curved model. Fits including the absorber give a reduced χ² close to one and of the same order as the one found for the log-parabola/broken power-law case (see Table 3). Results using the F-test favor a log-parabola fit for J013126−100931, while J212912−153841 is best fitted by a broken power-law model with intrinsic absorber. The radiating electron population is therefore likely in both cases to possess an intrinsic break at low energies, with J212912−153841 also showing the signature of an intervening absorbing medium.
• 2.  
  The source J151002+570243 represents an exception. Analyzing residuals and F-test results, the broken power-law model plus absorber is favored over the broken power-law (p-value = 8 × 10⁻⁵) or power-law with absorber (p-value = 2 × 10⁻³) ones. As can be seen from Figure 1, the spectrum indeed shows a break around 0.8 keV, after which it becomes softer, and a second break around 6 keV, after which it hardens again. While the second break can be explained via detailed modeling of the jet properties (see Section 4), the first one can hint to either the presence of an intrinsic absorber or a break in the electron population itself. To check whether the second option is viable, we further test a double broken power-law model¹⁸ (bkn2po in XSPEC). The fit is again preferred over the broken power-law (p-value = 6 × 10⁻⁵) or power-law with absorber (p-value = 10⁻³) models. Both the broken power-law plus intrinsic absorber and double broken power-law fits give similar reduced χ² values, and the F-test does not favor one over the other.
• 3.  
  For J020346+113445 and J135406–020603, due to low soft X-ray statistics, a simple power-law model with absorption fixed at the Galactic value is employed.

In order to cross calibrate the instruments, we include a multiplicative constant factor, fixing it equal to 1 for FPMA but leaving it free to vary for FPMB and soft X-ray instruments. For all sources with good X-ray statistics, the constants between instruments are of the order of ∼1%–10% from the fixed FPMA value. Instead, for the ones with poor soft X-ray statistics, the difference for FPMB is in the range 0.6%–10%, while a larger offset in the range 10%–30% is found for Swift-XRT.

Zoom In Zoom Out Reset image size

Table 3.  Table of X-Ray Best-fit Spectral Parameters

J013126−100931
Model	Best-fit Parameters					Fluxb	χ2/d.o.f.	F-statc	p-valuesc
PL	Γ
	1.40 ± 0.05					$(1.08\pm 0.05)\times {10}^{-11}$	192.61/179
BPL	Γ1	Γ2	Eb (keV)
	${1.13}_{-0.3}^{+0.1}$	1.52 ± 0.1	${5.3}_{-3.4}^{+2.3}$			${9.50}_{-0.6}^{+0.5}\times {10}^{-12}$	175.22/177	8.78	2 × 10−4
LPa	α	β	E0
	1.02 ± 0.17	0.23 ± 0.1	1(scale)			${\boldsymbol{(}}{\bf{9.01}}{\boldsymbol{\pm }}{\bf{0.5}}{\boldsymbol{)}}{\boldsymbol{\times }}{{\bf{10}}}^{-{\bf{12}}}$	176.29/178	16.47	7 × 10−5
PL+ABS	Γ	${N}_{{\rm{H}}}(\mathrm{intr}.)$
	1.45 ± 0.06	${4.52}_{-2.7}^{+3.6}\times {10}^{22}$				${1.04}_{-0.30}^{+0.04}\times {10}^{-11}$	183.46/178	8.77	3 × 10−3
BPL+ABS	Γ1	Γ2	Eb (keV)	${N}_{{\rm{H}}}(\mathrm{intr}.)$
	${1.18}_{-0.18}^{+0.16}$	$1.56\pm 0.10$	${6.23}_{-3.35}^{+1.91}$	${0.63}_{-0.58}^{+3.39}\times {10}^{22}$		${9.45}_{-0.72}^{+0.46}\times {10}^{-12}$	174.83/176	5.96	6 × 10−4
J064632+445116
Model	Best-fit Parameters					Fluxb	χ2/d.o.f.	F-statc	p-valuesc
PLa	Γ
	1.63 ± 0.07					${{\bf{2.27}}}_{{\boldsymbol{-}}{\bf{0.17}}}^{{\boldsymbol{+}}{\bf{0.12}}}{\boldsymbol{\times }}{{\bf{10}}}^{{\boldsymbol{-}}{\bf{12}}}$	37.82/48
BPL	Γ1	Γ2	Eb (keV)
	${1.51}_{-0.14}^{+0.12}$	${1.76}_{-0.09}^{+0.3}$	${4.11}_{-2.41}^{+14.6}$			${2.06}_{-0.45}^{+0.17}\times {10}^{-12}$	34.54/46	2.18	0.12
LP	α	β	E0
	${1.48}_{-0.2}^{+0.1}$	${0.11}_{-0.15}^{+0.17}$	1(scale)			${2.04}_{-0.31}^{+0.22}\times {10}^{-12}$	36.63/47	1.52	0.22
J151002+570243
Model	Best-fit Parameters					Fluxb	${\chi }^{2}/$ d.o.f.	F-statc	p-valuesc
PL	Γ
	1.47 ± 0.05					${1.95}_{-0.09}^{+0.15}\times {10}^{-12}$	391.61/276
BPL	Γ1	Γ2	Eb (keV)
	1.41 ± 0.02	$0.98\pm 0.20$	${6.25}_{-1.10}^{+2.15}$			$(3.36\pm 0.50)\times {10}^{-12}$	381.27/274	3.71	2 × 10−2
LP	α	β	E0
	1.44 ± 0.05	−0.02 ± 0.08	1 (scale)			$(2.05\pm 0.30)\times {10}^{-12}$	388.32/275	2.32	10−1
PL+ABSa	Γ			${N}_{{\rm{H}}}(\mathrm{intr}.)$
	1.48 ± 0.04			$(0.73\pm 0.31)\times {10}^{22}$		$(1.78\pm 0.21)\times {10}^{-12}$	376.07/275	5.38	8 × 10−4
BPL+ABSa	Γ1	Γ2	Eb (keV)	${{\bf{N}}}_{{\bf{H}}}{\boldsymbol{(}}{\bf{intr}}{\boldsymbol{.}}{\boldsymbol{)}}$
	1.52 ± 0.05	${{\bf{0.93}}}_{{\boldsymbol{-}}{\bf{0.7}}}^{{\boldsymbol{+}}{\bf{1.17}}}$	${{\bf{6.02}}}_{{\boldsymbol{-}}{\bf{1.7}}}^{{\boldsymbol{+}}{\bf{1.1}}}$	${\boldsymbol{(}}{\bf{0.86}}{\boldsymbol{\pm }}{\bf{0.32}}{\boldsymbol{)}}{\boldsymbol{\times }}{{\bf{10}}}^{{\bf{22}}}$		${\boldsymbol{(}}{\bf{3.26}}{\boldsymbol{\pm }}{\bf{0.5}}{\boldsymbol{)}}{\boldsymbol{\times }}{{\bf{10}}}^{-{\bf{12}}}$	360.27/273	7.91	4 × 10−5
2BPLa	Γ1	Γ2	${{\boldsymbol{\Gamma }}}_{{\bf{3}}}$	${{\bf{E}}}_{{\bf{b}}{\bf{1}}}{\boldsymbol{(}}{\bf{keV}}{\boldsymbol{)}}$	${{\bf{E}}}_{{\bf{b}}{\bf{2}}}{\boldsymbol{(}}{\bf{keV}}{\boldsymbol{)}}$
	${{\bf{1.02}}}_{{\boldsymbol{-}}{\bf{0.34}}}^{{\boldsymbol{+}}{\bf{0.27}}}$	1.49 ± 0.04	0.93 ± 0.20	${{\bf{0.85}}}_{{\boldsymbol{-}}{\bf{0.10}}}^{{\boldsymbol{+}}{\bf{0.28}}}$	${{\bf{6.09}}}_{{\boldsymbol{-}}{\bf{1.50}}}^{{\boldsymbol{+}}{\bf{1.24}}}$	${\boldsymbol{(}}{\bf{3.41}}{\boldsymbol{\pm }}{\bf{0.5}}{\boldsymbol{)}}{\boldsymbol{\times }}{{\bf{10}}}^{-{\bf{12}}}$	355.18/272	6.97	8 × 10−5
J212912−153841
Model	Best-fit Parameters					Fluxb	${\chi }^{2}/$ d.o.f.	F-statc	p-valuesc
PL	Γ
	1.38 ± 0.02					${3.53}_{-0.07}^{+0.08}\times {10}^{-11}$	813.92/631
BPL	Γ1	Γ2	Eb (keV)
	0.91 ± 0.09	1.56 ± 0.03	${1.75}_{-0.17}^{+0.30}$			${3.00}_{-0.08}^{+0.06}\times {10}^{-11}$	542.37/629	157.462	10−56
LP	α	β	E0
	1.10 ± 0.04	0.25 ± 0.20	1(scale)			${2.57}_{-0.07}^{+0.08}\times {10}^{-11}$	569.69/630	270.085	10−50
PL+ABS	Γ	${N}_{{\rm{H}}}(\mathrm{intr}.)$
	1.51 ± 0.02	$2.16\pm 0.33\times {10}^{22}$				$(3.09\pm 0.09)\times {10}^{-11}$	560.98/630	284.06	10−53
BPL+ABSa	Γ1	Γ2	Eb (keV)	${{\bf{N}}}_{{\bf{H}}}{\boldsymbol{(}}{\bf{intr}}{\boldsymbol{.}}{\boldsymbol{)}}$
	${{\bf{1.22}}}_{{\boldsymbol{-}}{\bf{0.13}}}^{{\boldsymbol{+}}{\bf{0.09}}}$	1.57 ± 0.03	2.21 ± 0.35	${{\bf{0.98}}}_{{\boldsymbol{-}}{\bf{0.46}}}^{{\boldsymbol{+}}{\bf{24}}}{\boldsymbol{\times }}{{\bf{10}}}^{{\bf{22}}}$		${{\bf{2.95}}}_{{\boldsymbol{-}}{\bf{0.10}}}^{{\boldsymbol{+}}{\bf{0.04}}}{\boldsymbol{\times }}{{\bf{10}}}^{{\boldsymbol{-}}{\bf{11}}}$	528.13/628	113.278	10−58
J020346+113445
Model	Best-fit Parameters					Fluxb	C-stat/d.o.f.
PL	Γ
	1.59 ± 0.13					${3.26}_{-0.33}^{+0.25}\times {10}^{-12}$	376.90/456
J135406−020603
Model	Best fit parameters					Fluxb	C-stat/d.o.f.
PL	Γ
	1.32 ± 0.16					${1.99}_{-0.32}^{+0.19}\times {10}^{-12}$	398.68/477

Notes.

^aThe best-fit model extracted from the analysis is highlighted in bold. ^bObserved flux in units of 10⁻¹² erg cm⁻² s⁻¹ in the 0.3–79 keV energy band. The errors are at the 90% level of confidence for one parameter of interest and the fluxes are corrected for Galactic absorption. ^cThe F-statistic and p-values are calculated for the different models with respect to the power-law case.

Download table as:  ASCIITypeset image

3.3. Central Black Hole Mass

Black hole masses are computed following two approaches. The first is through single epoch optical spectroscopy (e.g., Shen et al. 2011), assuming a virialized BLR and adopting emission line parameters such as full width at half maximum, luminosity, and continuum emission. The errors associated with this method are of the order ∼0.4 dex (e.g., Vestergaard & Peterson 2006; Shen et al. 2011). The second is by modeling the accretion disk using a standard Shakura & Sunyaev (1973) model. Blazars' accretion disk emission falls in the IR-UV band and, if measurable, it enables us to derive the black hole mass, even when optical spectra are not available. High-redshift blazars, thanks to the shift of the low-energy SED peak, unveil this emission and are therefore ideal sources for applying this method. With good IR-UV coverage, the uncertainty associated with this approach is of about a factor of 2 and the results have been shown to be in good agreement with the virial method (see, e.g., Castignani et al. 2013; Ghisellini & Tavecchio 2015; Paliya et al. 2017a). The accretion disk SED is assumed to be a multicolor blackbody. The flux density profile of this emission is given by the following (Frank et al. 2002):

$$\begin{eqnarray}&&{F}_{\nu }={\nu }^{3}\displaystyle \frac{4\pi h\cos {\theta }_{{\rm{v}}}}{{c}^{2}{D}^{2}}{\int }_{{R}_{\mathrm{in}}}^{{R}_{\mathrm{out}}}\displaystyle \frac{R\,{dR}}{{e}^{h\nu /{kT}(R)}-1},\end{eqnarray} \tag{ 1 }$$

where k, c, and h are the Boltzmann constant, the speed of light, and the Plank constant, respectively, D is the luminosity distance, ${\theta }_{{\rm{v}}}$ is the jet viewing angle and ${R}_{\mathrm{in}}$ and ${R}_{\mathrm{out}}$ are the inner and outer disk radii. We take ${R}_{\mathrm{in}}=3{R}_{\mathrm{Sch}}$ and ${R}_{\mathrm{out}}\,=500{R}_{\mathrm{Sch}}$ (R_Sch is the Schwarzschild radius). The radial temperature profile is given by:

$$\begin{eqnarray}&&T(R)=\displaystyle \frac{3{R}_{\mathrm{Sch}}{L}_{\mathrm{disk}}}{16\pi {\eta }_{{\rm{a}}}{\sigma }_{\mathrm{SB}}{R}^{3}}{\left[1-{\left(\displaystyle \frac{3{R}_{\mathrm{Sch}}}{R}\right)}^{1/2}\right]}^{1/4},\end{eqnarray} \tag{ 2 }$$

where ${\sigma }_{\mathrm{SB}}$ is the Stefan–Boltzmann constant and ${\eta }_{{\rm{a}}}$ is the accretion efficiency (we consider ${\eta }_{{\rm{a}}}=0.1$). With the above formulation, there are only two free parameters: the black hole mass (${M}_{\mathrm{BH}}$) and the accretion rate (${\dot{M}}_{{\rm{a}}}$). The latter is related to the intrinsic accretion luminosity via ${L}_{\mathrm{disk}}={\eta }_{{\rm{a}}}{\dot{M}}_{{\rm{a}}}{c}^{2}$. If the peak of the disk emission is visible, it is possible to obtain L_disk from the observations. The only free parameter left, therefore, is ${M}_{\mathrm{BH}}$, which scales as ${M}_{\mathrm{BH}}\,\propto {\nu }_{\mathrm{peak},\mathrm{disk}}^{-2}{{L}_{\mathrm{disk}}}^{1/2}$ (where ${\nu }_{\mathrm{peak},\mathrm{disk}}$ is the peak of the disk emission, see Calderone et al. 2013).

For our sources, we looked into the literature to find the black hole masses derived through the first method, i.e., optical spectroscopy, then we performed disk modeling to extract them through the second. Both masses are reported in Table 4. They are found to be in reasonable accordance, and match within a factor of 2, apart from J151002+570243, whose difference cannot be accounted for by considering the associated errors. In fact, we note that the SARA magnitudes detected for this source are larger (∼a factor of 2 higher in flux) than the average archival ones. This could be linked to the fact that we are observing the source in a high state (see, e.g., Kelly et al. 2009; Ruan et al. 2014; Rumbaugh et al. 2018), hence our measurement of the disk luminosity could be higher than the low/average one. However, the peak frequency of the emission is consistent with that found in the archival data, therefore the estimate on the black hole mass is not strongly affected. In the case of J013126–100931 we find that both SARA magnitudes are larger¹⁹ and the peak of the emission is at lower frequencies with respect to archival data. This could lead to an overestimate of the black hole mass, although the value would still be within the errors associated with the modeling approach. Longer term optical monitoring would be needed to study variability of the disk emission.

Table 4.  Table of Black Hole Masses, Derived both from Spectroscopic Approach and SED Modeling

	MBH,SED(M☉)	MBH,spectroscopy (M☉)
J013126−100931	6.0 × 109	⋯
J020346+113445	3.0 × 109	⋯
J064632+445116	1.5 × 109	1.3 × 109
J135406−020603	1.0 × 109	8.9 × 108
J151002+570243	6.5 × 109	3.1 × 108
J212912−153841	5.0 × 109	6.3 × 109

Download table as:  ASCIITypeset image

3.4. Broadband SED Modeling

We reproduce the broadband SED of all six sources following a simple one-zone leptonic emission model (see Dermer et al. 2009; Ghisellini & Tavecchio 2009) and here we explain it briefly. We assume that radiating electrons are confined in a spherical region canvassing the full jet cross section. This region is considered to be located at a distance R_diss from the central black hole, moving with bulk Lorentz factor Γ. The jet semi-opening angle is taken to be 0.1 rad. The energy distribution of relativistic electrons is modeled as a broken power law of the following type:

$$\begin{eqnarray}&&N(\gamma )\,\propto \,\displaystyle \frac{{\left({\gamma }_{\mathrm{break}}\right)}^{-p}}{{\left(\gamma /{\gamma }_{\mathrm{break}}\right)}^{p}+{\left(\gamma /{\gamma }_{\mathrm{break}}\right)}^{q}}\end{eqnarray} \tag{ 3 }$$

where p and q are, respectively, the slopes of the distribution before and after the break energy, γ_break. Relativistic electrons, immersed in a randomly oriented but uniform magnetic field (B), emit synchrotron radiation and, in the presence of a photon field, they also undergo IC process. We consider the sources of radiation for IC to be either the same photons emitted via synchrotron or radiation fields external to the jet. For the latter, following Ghisellini & Tavecchio (2009), we take into account (1) the accretion disk emission, (2) the X-ray corona located above the accretion disk, (3) the BLR, and (4) the dusty torus. The corona is assumed to have a cutoff power-law spectral distribution and to reprocess 30% of the disk luminosity. The BLR is considered to be a spherical shell located at a distance ${R}_{\mathrm{BLR}}={10}^{17}{L}_{\mathrm{disk},45}^{1/2}$ cm (${L}_{\mathrm{disk},45}$ being the accretion disk luminosity in units of 10⁴⁵ erg s⁻¹) from the central engine, and reemitting 10% of L_disk (Bentz et al. 2006; Kaspi et al. 2007; Tavecchio & Ghisellini 2008). The torus is also assumed to be a spherical shell, located further away at ${R}_{\mathrm{TORUS}}\,={10}^{18}{L}_{\mathrm{disk},45}^{1/2}$ cm and reprocessing 50% of L_disk (Lawrence 1991; Nenkova et al. 2002; Cleary et al. 2007). The spectral shapes of both the BLR and the torus are considered to be blackbodies peaking, respectively, at the Lyα frequency and at the characteristic torus temperature (T_TORUS = 300 K). Calculating the radiative energy densities of all four components, we are able to extract the IC spectra (both SSC and EC). In order to calculate the total jet power, we make use of the following assumptions: the kinetic power of the jet is carried only by protons considered to be cold and therefore contributing only to the inertia of the jet; and the number densities of relativistic electrons and protons are equal (e.g., Celotti & Ghisellini 2008). The total jet power is evaluated as the sum of electrons (${P}_{{\rm{e}}}$), protons (${P}_{{\rm{p}}}$), magnetic field (${P}_{{\rm{m}}}$), and radiation (${P}_{{\rm{r}}}$) powers, which are all calculated.

The simultaneous multiwavelength coverage obtained for all our sources is crucial to accurately determine the input parameters to our model (i.e., L_disk, R_diss, ${M}_{\mathrm{BH}}$, p, q, Γ, B, low- and high-energy cutoff of the electron distribution (γ_min, γ_max)). Here we list some guidelines behind our choice:

• 1.  
  Thermal component: optical data of these sources sample the accretion disk emission. It is therefore possible to gauge the peak of this radiation and hence the overall disk luminosity (${L}_{\mathrm{disk}}\sim 2{\nu }_{\mathrm{peak},\mathrm{disk}}L({\nu }_{\mathrm{peak},\mathrm{disk}})$). Following the discussion in Section 3.3, optical data also allow us to constrain ${M}_{\mathrm{BH}}$.
• 2.  
  Nonthermal component: a blazar's emission region in this one-zone model is required to be compact and located within the BLR or the torus. As a consequence, synchrotron emission from this region is depleted by synchrotron self-absorption below ∼10¹² Hz. It follows that the low-energy synchrotron component of the SED is likely to be produced by another region further away along the jet, and therefore cannot be used to constrain the shape of the electron energy distribution. However, good quality X- and γ-rays allow us to exactly determine peak of the IC emission and the shape of the nonthermal continuum, which is directly linked to the shape of the electron distribution. Hence, availability of both NuSTAR and Fermi-LAT data are key to fix both p and q (see Equation (3)). Furthermore the interplay of Γ, B, and R_diss reflects on the level of the high-energy emission. Therefore, the available X-to-γ-ray data allow us to accurately estimate these parameters.

We report the parameters derived from the SED modeling in Table 5. Figure 2 shows the results for all six sources, which are discussed in the next Section.

Zoom In Zoom Out Reset image size

Table 5.  Summary of the Parameters Used/Derived from the SED Modeling of Six MeV Blazars Shown in Figure 2

Parameter	J013126−100931	J020346+113445	J064632+445116	J135406−020603	J151002+570243	J212912−153841
Slope of the particle distribution below the break energy (p)	1.5	1.8	2.05	0.9	0.5	2.1
Slope of the particle distribution above the break energy (q)	3.8	4.2	3.6	3.9	3.8	4.0
Magnetic field in Gauss (B)	0.5	0.9	1.0	1.2	1.2	0.8
Particle energy density in erg cm−3 (Ue)	0.003	0.02	0.004	0.01	0.01	0.009
Bulk Lorentz factor (Γ)	15.8	10.0	14.0	10.0	11.0	10.0
Minimum Lorentz factor (γmin)	12	1	1	1	1.7	4.0
Break Lorentz factor (γbreak)	74.07	83.85	131.22	59.80	32.42	60.32
Maximum Lorentz factor (γmax)	3e3	3e3	3e3	3e3	5e3	3e3
Dissipation distance (Rdiss,) in parsecs (RSch)	0.43 (750)	0.28 (1000)	0.57 (4000)	0.25 (2700)	0.31 (500)	0.98 (2050)
Size of the BLR in parsec (${R}_{\mathrm{BLR}}$, RSch)	0.22 (399)	0.23 (837)	0.38 (2671)	0.23 (2511)	0.27 (438)	0.72 (1514)
Accretion disk luminosity in log scale (Ldisk, erg s−1)	46.70	46.74	47.15	46.74	46.85	47.71
Accretion disk luminosity in Eddington units (${L}_{\mathrm{disk}}/{L}_{\mathrm{Edd}}$)	0.06	0.14	0.74	0.43	0.09	0.79
Jet power in electrons in log scale (${P}_{{\rm{e}}}$, erg s−1)	45.19	45.29	45.37	44.93	45.13	45.89
Jet power in magnetic field in log scale (${P}_{{\rm{B}}}$), erg s−1	45.61	45.37	46.36	45.53	45.77	46.33
Radiative jet power in log scale (${P}_{{\rm{r}}}$, erg s−1)	45.49	46.11	46.27	46.59	46.37	46.70
Jet power in protons in log scale (${P}_{{\rm{p}}}$, erg s−1)	46.90	47.78	47.98	46.89	47.12	48.11
Total jet power in log scale (${P}_{\mathrm{TOT}}$, erg s−1)	46.93	47.78	47.99	46.91	47.14	48.12

Note. A viewing angle of 3° is adopted for all of them.

Download table as:  ASCIITypeset image

High-redshift MeV blazars are among the most extreme sources of the blazar class. Harboring powerful relativistic jets and supermassive black holes at their center, they provide a unique opportunity to study jetted sources in the early universe and are ideal to address open issues such as the accretion disk-jet connection and the growth and evolution of supermassive black holes at the dawn of the universe. Since in these sources the high-energy radiation dominates the bolometric luminosity, in order to understand the physics of the jets and the properties of the emitting leptonic population, good coverage at these frequencies is required. Awaiting an all sky MeV instrument (e.g., AMEGO, Rando 2017; e-ASTROGAM, de Angelis et al. 2018) to sample the peak of the IC emission, the most effective way to study the high-energy part of the SED is via a combination of X- and γ-ray data. As a result of the drift of the nonthermal SED peaks to lower frequencies, these high-redshift sources are brighter at X-ray frequencies, displaying a hard X-ray photon index (Γ_X ≲ 1.5) and a soft γ-ray one (${{\rm{\Gamma }}}_{\gamma }\gtrsim 2.3$). Four of the objects in our study already had γ-ray data available (Ackermann et al. 2017) and, most importantly, we were able to obtain NuSTAR observations for all six of them. All our sources show X-ray photon indices <2, as can be seen in Table 3, and soft γ-ray ones (${{\rm{\Gamma }}}_{\gamma }\gt 2.5$, Ackermann et al. 2017). We point out that the sources J013126–100931 and J212912–153841 are detected in the BAT 105-months catalog (Oh et al. 2018), with a ${{\rm{\Gamma }}}_{15-150\mathrm{keV}}$ of ${1.81}_{-0.48}^{+0.53}$ and ${1.79}_{-0.43}^{+0.48}$, respectively. The unprecedented sensitivity of NuSTAR allows us to accurately measure the hard X-ray photon index (${{\rm{\Gamma }}}_{{\rm{X}}}\sim 1.55\pm 0.1$ and ${{\rm{\Gamma }}}_{{\rm{X}}}\sim 1.56\pm 0.03$, respectively, see Table 3). Furthermore, NuSTAR data are key to sample the rising part of the IC spectrum up to 60–70 keV and, in combination with Fermi ones, enable us to better constrain the peak location and consequently the bulk Lorentz factor and the shape of the electron energy distribution underlying this powerful emission. To highlight the importance of NuSTAR, in Figures 3–5 we show the zoomed SED for the three targets studied by Ackermann et al. (2017). In these plots, the orange dashed line represents the EC model taken from their paper, while our EC model is shown as the blue solid line. Note that Ackermann et al. (2017) performed the modeling using only X-ray archival observations (shown by orange data points), whereas we benefited from the availability of NuSTAR observations (red stars). SED parameters derived in this work and in Ackermann et al. (2017) are also reported in Tables 6–8 for comparison. It is evident how NuSTAR data enable us to sample the IC emission closer to the peak than previously possible, better constraining the slopes of the radiating leptonic population. For example, looking at J151002+570243 it can be seen how the value of p drastically changes from 1.8 to 0.5, almost a factor of four difference from the previous modeling. These constraints translate into a more precise determination of the IC peak location and luminosity. We point out, as extreme examples, how our modeling predicts the IC peak of J151002+570243, or the IC peak luminosity of J064632+445116, to be nearly an order of magnitude lower than previously found only using soft X-ray data. The location of the emission regions is now more accurately constrained and found for all three sources to lie outside the BLR. These factors in turn produce a more accurate measure of the jet power and particle energy distribution. The substantial difference between the two EC models reflects the necessity of having hard X-ray coverage.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

NuSTAR has also proven to be essential in high-redshift blazar studies since, combined with soft X-ray observations, it has enabled us to detect peculiar X-ray features of these sources, such as spectral variability (see Tagliaferri et al. 2015; Sbarrato et al. 2016) and spectral flattening (see Paliya et al. 2016). The first property is important for various reasons. It could relate to intrinsic changes in the electron distribution (hence injected power) or could be key to unveil acceleration processes happening in the jet, such as magnetic reconnection or shocks (e.g., Spada et al. 2001; Moraitis & Mastichiadis 2011; Rueda-Becerril et al. 2015; Christie et al. 2019). Furthermore, it has been noticed that often X-ray variability correlates with γ-ray variability, indicating that the assumption of single-electron population models is reliable. However, soft X-ray observations have sometimes not been reflective of this variation. On one hand, sensitivity limits of currently operating instruments could prohibit us from seeing such variability. On the other, the low-energy electron population emitting at these frequencies is not expected to vary significantly (see Section 3.1), while hard X-ray emission could show the signature of such a variation (see discussion in Sbarrato et al. 2016). Moreover as these sources become fainter in the γ-rays, observing variability in this high-energy regime becomes a challenge (Li et al. 2018). For all our sources, NuSTAR data do not show any variability in the observation period. In the soft X-ray regime, we find both flux and photon index variability only in the case of J212912–153841 during the different observing epochs.

The second property instead finds its cause in different scenarios. It could be related to the intrinsic behavior of the emitting electron population (see Paliya et al. 2016) or could be attributed to an intrinsic absorber along the line of sight, such as the warm-hot IGM (see Arcodia et al. 2018), or a combination of the two. The combination of soft X-ray and NuSTAR data enables us to detect spectral curvature in the X-ray spectra of three of our sources (see Section 3.2). Within the available statistics, it appears that the softening below ∼2 keV in J212912–153841 is a combination of both an intrinsic absorber probably located along the line of sight and a break in the electron distribution. The curvature below ∼5 keV for J013126–100931 instead is best explained by an intrinsic curvature of the leptonic distribution. Importantly, the shape of the X-ray spectrum constrains the behavior of the electrons underlying this emission. Indeed, if the X-rays are produced via EC rather than SSC, it is possible to constrain the electrons' low-energy cutoff, as can be seen in Figure 6. In fact, for values of ${\gamma }_{\min }\gt 1$, our model predicts a significant break which is seen in three of our sources. It follows that joint soft X-ray and NuSTAR observations are crucial in determining the minimum energy of the leptonic population underlying the emission, which in turns regulates the jet power.

Zoom In Zoom Out Reset image size

J151002+570243, the farthest blazar ever detected in γ-rays, is also the most peculiar of our sources in the X-ray regime. Indeed its X-ray spectrum not only favors a softening around ∼1 keV, but also shows a hardening >6 keV, in the energy range covered by NuSTAR, with an extreme photon index of 0.94. We note that previous studies of this source in the soft X-ray regime found it to have a hard X-ray photon index (${{\rm{\Gamma }}}_{{\rm{X}}}=1.55\pm 0.05$, Wu et al. 2013). This is another instance that reflects the necessity of NuSTAR data to constrain the true shape of the X-ray emission, and hence all jet parameters associated with it.

From a modeling point of view, all of our sources have SEDs comparable to the most powerful LSP FSRQ blazars. They show the typical nonthermal humps, as well as the one produced by the accretion disk, unveiled by the SED frequency shift. The low-energy emission is therefore well interpreted by synchrotron process and peaks in the submillimeter regime. Thanks to the optical data obtained with SARA we are able to detect the peak of the disk emission and therefore to estimate disk luminosity and black hole masses from modeling.²⁰ The high-energy part of the SED, covered by Swift-XRT, NuSTAR, and Fermi, instead peaks in (or close to) the MeV energy range and is explained by EC radiation. The X-ray spectra are very hard, which in turn suggests that the EC process is dominant with respect to the SSC (see Ajello et al. 2016 for a detailed discussion). In fact, the SSC is predicted by our modeling to be low for all of the sources, as can be seen in Figure 2. In Figure 7, the six sources here are compared to the Fermi-LAT blazars with available NuSTAR observations (L. Marcotulli et al. 2020, in preparation). Our sources fall in the same region as the most powerful FSRQs in both plots, showing high X-ray luminosities and hard X-ray indices as well as high γ-ray luminosities and soft γ-ray indices. In all six of them, the IC peak luminosity dominates the synchrotron one (CD > 1). In the adopted model, the radiative energy densities of the different AGN components are a function of the distance of the emission region from the central black hole. It follows that we can determine the location of the emission region. For our sources we find it resides outside the BLR and inside the torus (thus in agreement with Costamante et al. 2018). Therefore the torus is the predominant source of radiation for the EC. Overall, our sources resemble the most extreme MeV blazars.

Zoom In Zoom Out Reset image size

It has been found that, for γ-ray detected blazars up to z = 3, the accretion disk luminosity is lower than the jet power (Ghisellini et al. 2014). More recent results show that a positive correlation remains for γ-quiet sources (i.e., blazars not detected in γ-rays) beyond z = 3 (Paliya et al. 2017a). However, as discussed above, having γ-ray data is crucial to accurately determine the jet power. Therefore, with good γ- and X-ray data in hand, as well as optical ones, we can insert our sources in the context of the disk-jet relation. In Table 5 we see how all six have jet powers exceeding disk luminosities. This further hints to the idea that there has to be an extra reservoir of energy powering the jet rather than the accretion disk alone (e.g., energy could be supplied by the spinning black hole, Blandford & Znajek 1977). However, we need to be cautious and take into account uncertainties in our measurements. For example, in our treatment we do not consider electron-positron pairs, which, if present, would reduce the estimated jet power by a factor n_e/n_p (${n}_{e}={n}_{{e}^{+}}+{n}_{{e}^{-}}$, see Pjanka et al. 2017). Nevertheless, to avoid the Compton Rocket effect that would otherwise break the jet (see Ghisellini & Tavecchio 2010), the number of pairs could not outnumber by a large amount the number of protons. Moreover, another effect that could reduce the jet power is if we consider it to have a spine-sheath structure (Sikora et al. 2016). An underestimation of the accretion disk luminosity could also lead to a false positive relation. So these caveats have to be kept in mind when analyzing the relation between jet power and disk luminosity in blazars.

An important jet parameter that is obtained via modeling is the bulk Lorentz factor. In fact, assuming that the jet points toward the observer at an angle $\lesssim 1/{\rm{\Gamma }}$, by a simple geometrical argument, we can estimate the size of the parent population of these sources. It turns out that the detection of one jetted source implies the existence of 2Γ² similar sources, at the same redshift, with similar black hole mass, but with jets pointing away from the observer. From various studies, it has been found that high-redshift blazars host supermassive black holes. Using the 2Γ² correction, Sbarrato et al. (2015) derived the comoving number density of billion solar mass black holes and found that for radio-quiet blazars (i.e., that host weak or no jets) it peaks at $z\sim 2$, while for radio-loud blazars (i.e., that have strong relativistic jets) it peaks at $z\sim 4$. This fact challenges our understanding of black hole growth and evolution in the early universe since, as pointed out by Ghisellini (2015), simple accretion (even considering highly massive seeds at z = 20) cannot account for ${M}_{\mathrm{BH}}\gt {10}^{9}\,{M}_{\odot }$ at $z\gt 4$. On account of the fact that the number density strongly depends on the number of sources detected per redshift bin, finding more such sources is crucial to understand the physics of early black hole accretion. The apparent dichotomy in Sbarrato et al. (2015) between radio-quiet and radio-loud sources further hints to a connection between rapid black hole growth and the presence of powerful jets. From SED modeling, five of our sources have ${M}_{\mathrm{BH}}\gt {10}^{9}\,{M}_{\odot }$.²¹ We use the prescription in Ackermann et al. (2017) to estimate the number density (n) of these ${M}_{\mathrm{BH}}\gt {10}^{9}\,{M}_{\odot }$ radio-loud blazars in the redshift bin z = [3, 4].²² This number was derived by Sbarrato et al. (2015) to be 50 Gpc⁻³ and Ackermann et al. (2017) increased this estimate by 18 Gpc⁻³. From our calculation, we find $n\sim 37\,{\mathrm{Gpc}}^{-3}$ which leads the total space density value to be 87 Gpc⁻³. This seems to further confirm a connection between the radio-loud phase and rapid accretion. Only detecting more sources would improve our understanding of this rapid and early black hole accretion.

In this work, the multiwavelength analysis of six high-redshift (z > 3) blazars, four of which were recently discovered to be γ-ray loud by Ackermann et al. (2017), is presented. Employing the excellent capabilities of NuSTAR, we are able to precisely study X-ray properties of these sources. In combination with quasi-simultaneous Swift and Fermi-LAT observations, we accurately constrain their jet properties (i.e., jet power, underlying electron distribution, and location of the emission region). Simultaneous optical monitoring with SARA reveals the peak of the disk emission through which we can determine the accretion disk luminosities and black hole masses. The main results are summarized here:

• 1.  
  NuSTAR spectra are very hard (${\rm{\Gamma }}\lesssim 1.5$) and do not show variability. Combining NuSTAR and sensitive soft X-ray observations, we detect spectral flattening in the X-ray spectra of three sources. The curvature in J013126–100931 can be explained by a property of the underlying electron distribution (i.e., it depends on the minimum energy of the electron population). In the case of J212912–153841, it is attributed to a combination of both intrinsic curvature in the leptonic distribution and an intrinsic absorber in the IGM. The peculiar case of J151002+570243 (z = 4.31) reveals two breaks in its X-ray spectrum, the first possibly attributed to the minimum energy of the electrons producing this emission and/or an IGM absorber, while the second can be explained in the context of the one-zone leptonic emission model.
• 2.  
  Using a one-zone leptonic emission model, NuSTAR and Fermi-LAT spectra allow us to precisely pinpoint the position of the high-energy peak, which is well described by EC emission of the relativistic electrons in the jet. This enables us to determine jet parameters such as the jet power and the emission region location (for all sources outside the BLR), as well as constrain the shape of the underlying electron distribution.
• 3.  
  Inserting our sources in the context of accretion disk-jet connection, we find that all sources display jet powers larger than accretion disk luminosity, validating the possibility of an extra source of power to the jet, such as the spinning black hole.
• 4.  
  The optical data unveil the disk emission and enable us to estimate the disk luminosity and black hole mass. All sources have $M\gt {10}^{9}\,{M}_{\odot }$, which further raises the space density of supermassive black holes in the redshift bin z = [3, 4] to $87\,{\mathrm{Mpc}}^{-3}$.

L.M., V.P., and M.A. acknowledge funding under NASA contracts NNX17AC35G and 80NSSC20K0044.

The Fermi LAT Collaboration acknowledges generous ongoing support from a number of agencies and institutes that have supported both the development and the operation of the LAT as well as scientific data analysis. These include the National Aeronautics and Space Administration and the Department of Energy in the United States, the Commissariat à l'Energie Atomique and the Centre National de la Recherche Scientifique/Institut National de Physique Nucléaire et de Physique des Particules in France, the Agenzia Spaziale Italiana and the Istituto Nazionale di Fisica Nucleare in Italy, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), High Energy Accelerator Research Organization (KEK) and Japan Aerospace Exploration Agency (JAXA) in Japan, and the K. A. Wallenberg Foundation, the Swedish Research Council, and the Swedish National Space Board in Sweden.

Additional support for science analysis during the operations phase is gratefully acknowledged from the Istituto Nazionale di Astrofisica in Italy and the Centre National d'Études Spatiales in France. This work was performed in part under DOE Contract DE-AC02-76SF00515.

We thank the Swift team and the Swift PI (B. Cenko) for promptly scheduling and executing the observations.

Reported work is in part based on observations obtained with the SARA Observatory 1.0 m, 0.96 m, and 0.6 m telescopes respectively located at la Roque de los Muchachos Observatory, Canary Islands (SARA-ORM), at Kitt Peak National Observatory, Arizona (SARA-KP), and at Cerro Tololo, Chile (SARA-CT), which are owned and operated by the Southeastern Association for Research in Astronomy (saraobservatory.org).

Part of this work is based on archival data, software, or online services provided by the ASI Data Center (ASDC). This research has made use of the XRT Data Analysis Software (XRTDAS).

This Appendix contains the X-ray spectra of the six blazars as well as the residuals of the fits (see Figure 8) in the same order as listed in Table 3. The green dotted line is the best fit found analyzing the residuals as well as performing F-tests (see Section 3.2).

Zoom In Zoom Out Reset image size