Gamma-ray Blazars Within the First 2 Billion Years

TL;DR

This paper reports the discovery of five high-redshift gamma-ray blazars beyond z=3.1, including the most distant known at z=4.31, providing insights into early black hole growth and blazar evolution.

Contribution

First detection of gamma-ray emitting blazars beyond z=3.1 with detailed spectral analysis, revealing their extreme properties and implications for black hole growth in the early Universe.

Findings

Discovered five gamma-ray blazars beyond z=3.1, including the most distant at z=4.31.

These blazars have steep gamma-ray spectra and large Compton dominance.

Black holes in these blazars are estimated to be 10^8-10^10 solar masses.

Abstract

The detection of high-redshift ($z>$3) blazars enables the study of the evolution of the most luminous relativistic jets over cosmic time. More importantly, high-redshift blazars tend to host massive black holes and can be used to constrain the space density of heavy black holes in the early Universe. Here, we report the first detection with the \fermi-Large Area Telescope of five \gm-ray emitting blazars beyond $z=3.1$, more distant than any blazars previously detected in $\gamma$-rays. Among these five objects, NVSS J151002+570243 is now the most distant known \gm-ray emitting blazar at $z=4.31$. These objects have steeply falling \gm-ray spectral energy distributions (SEDs) and, those that have been observed in X-rays, a very hard X-ray spectrum, both typical of powerful blazars. Their Compton dominance (ratio of the inverse Compton to synchrotron peak luminosities) is also very large ($>20$). All of these properties place these objects among the most extreme members of the blazar population. Their optical spectra and the modeling of their optical-UV SEDs confirm that these objects harbor massive black holes ($M_{\rm BH} \sim 10^{8-10} {\rm M}_{\odot}$). We find that, at $z\approx4$, the space density of $>10^{9} {\rm M}_{\odot}$ black holes hosted in radio-loud and radio-quiet active galactic nuclei are similar, implying that radio-loudness may play a key role in rapid black hole growth in the early Universe.

Tables

Table 1. Table 1: Basic information and results of the Fermi -LAT data analysis of high-redshift blazars.

		Basic information									Fermi-LAT data analysis
Name	Radio position (J2000)		$F_{1.4\mathrm{GHz}}$	R	$z$	RL	$M_{\mathrm{BH},s}$	Optimized position (J2000)		$R_{95\%}$	$F_{0.06-300\mathrm{GeV}}$	${\Gamma }_{\gamma }$	$L_{\gamma }$	TS
(NVSS)	hh mm ss.ss	dd mm ss.s	(mJy)	(mag)			$M_{☉}$	hh mm ss.ss	dd mm ss.s	(degrees)	(10-8 $\mathrm{ph}{\mathrm{cm}}^{-2}{s}^{-1}$)		(1048 erg s-1)
J064632+445116	06 46 32.00	+44 51 17.0	452	18.5	3.4	1253	9.1$		+45 02 40	0.24	2.03$\pm$0.40	2.68$\pm$0.10	1.4$\pm$0.4	62
J135406$-$020603	13 54 06.90	$-$02 06 03.2	733	19.2	3.7	3741	8.9†		$-$02 06 43	0.22	2.32$\pm$0.52	2.88$\pm$0.14	2.5$\pm$0.8	44
J151002+570243	15 10 02.92	+57 02 43.4	202	19.9	4.3	1850	8.5†		+57 05 07	0.18	0.86$\pm$0.30	2.55$\pm$0.16	1.1$\pm$0.5	34
J163547+362930	16 35 47.24	+36 29 30.0	152	20.6	3.6	2842	8.7†		+36 29 45	0.13	5.31$\pm$0.96	3.15$\pm$0.14	7.0$\pm$0.2	151
J212912$-$153841	21 29 12.13	$-$15 38 41.0	590	16.5	3.3	263	9.8$		$-$15 34 25	0.33	2.65$\pm$0.43	2.82$\pm$0.10	1.9$\pm$0.4	66

Table 2. Table 2: Summary of the parameters used/derived from the modeling of the SED of the objects shown in Figure  3 . The viewing angle is taken as 3 ∘ for all of them.

Parameter	J0646+4451	J1510+5702	J2129$-$1538
Slope of the electron energy distribution before break energy ($p$)	1.8	1.8	2.2
Slope of the electron energy distribution after break energy ($q$)	4.4	4.1	4.5
Magnetic field in Gauss ($B$)	2.1	1.4	1.3
Particle energy density in erg cm-3 ($U_e^{\prime }$)	0.009	0.029	0.002
Bulk Lorentz factor ($\Gamma$)	12	11	14
Minimum Lorentz factor (${\gamma }_{\min }^{\prime }$)	1	1	1
Break Lorentz factor (${\gamma }_b^{\prime }$)	72	82	51
Maximum Lorentz factor (${\gamma }_{\max }^{\prime }$)	2e3	3e3	2e3
Size of the emission region in parsec ($R_{\mathrm{blob}}$)	0.025	0.017	0.059
Dissipation distance in parsec ($R_{\mathrm{diss}}$)	0.25	0.17	0.59
Size of the BLR in parsec ($R_{\mathrm{BLR}}$)	0.37	0.22	0.79
Black hole mass in log scale, in units of solar mass ($M_{\mathrm{BH},m}$)	9.60	9.48	9.85
Accretion disk luminosity in log scale ($L_{\mathrm{disk}}$, erg s-1)	47.11	46.65	47.78
Accretion disk luminosity in Eddington units ($L_{\mathrm{disk}}/L_{\mathrm{Edd}}$)	0.26	0.12	0.68
Fraction of the disk luminosity reprocessed by the BLR ($f_{\mathrm{BLR}}$)	0.1	0.1	0.1
Fraction of the disk luminosity reprocessed by the IR-torus ($f_{\mathrm{IR}}$)	0.5	0.5	0.5
Compton dominance (CD)	25	47	99
Jet power in electrons in log scale ($P_e$, erg s-1)	44.88	44.94	45.14
Jet power in magnetic field in log scale ($P_B$), erg s-1	46.15	45.38	46.58
Radiative jet power in log scale ($P_r$, erg s-1)	45.89	45.70	46.31
Jet power in protons in log scale ($P_p$, erg s-1)	47.38	47.45	47.89

Full text

galaxies: active — gamma rays: galaxies— galaxies: jets— galaxies: high-redshift— quasars: general

††facilities: Fermi-LAT

1 Introduction

Blazars are the most powerful Active Galactic Nuclei (AGN) with relativistic jets oriented close to the line of sight. The jet radiation dominates their broadband emission, especially at $\gamma$-rays where the inverse Compton (IC) hump of the blazar spectral energy distribution (SED) peaks around tens or hundreds of MeV. The Large Area Telescope (LAT) onboard the Fermi satellite (Atwood et al., 2009) has already detected thousands of blazars, thus confirming that they are the most numerous population in the $\gamma$-ray sky (e.g., Ackermann et al., 2015). Nonetheless, high-redshift blazars above a redshift of 3.1 are missing in the Fermi catalogs, possibly due to the shift of the IC peak to lower frequencies in which LAT is less sensitive. The newly released Pass 8 photon data set (Atwood et al., 2013), with an improved event-level analysis, substantially enhances the sensitivity of the LAT at all energies and in particular at lower energies (i.e., $<$200 MeV). This increases the capability of the LAT to detect spectrally soft, potentially high-$z$ blazars.

These objects typically have large bolometric luminosities ($L_{\rm bol.}>$1048 erg s*-1*) and host powerful relativistic jets ($P_{j}\gtrsim\dot{M}c^{2}$ for a given accretion efficiency; e.g., Ghisellini et al., 2014). In general, they harbor extremely massive black holes ($M_{\rm BH}\sim 10^{9}{\rm M}_{\odot}$; Ghisellini et al., 2010). Since blazars are highly beamed, the detection of a single blazar implies the existence of $2\Gamma^{2}$ (i.e., $\sim$400$-$600, $\Gamma$ is the bulk Lorentz factor) misaligned blazars with similar properties. Therefore, the detection of high-$z$ blazars can constrain models of supermassive black hole formation in the early Universe (see, e.g., Volonteri et al., 2011; Ghisellini et al., 2013). This suggests that the detection of new high-$z$ blazars will test the hypotheses of blazar evolution, since these high redshifts constrain the time available for such extreme objects to grow.

Motivated by the recent release of the sensitive Pass 8 dataset, we perform a systematic search for new $\gamma$-ray emitters beyond $z=3.1$ and in this Letter we report the first detection of five $z>3.1$ $\gamma$-ray emitting blazars. Throughout, we adopt a $\Lambda$CDM cosmology with the Hubble constant $H_{0}=71$ km s*-1* Mpc*-1*, $\Omega_{m}=0.27$, and $\Omega_{\Lambda}=0.73$.

2 Sample selection and Analysis

In order to search for high-$z$ blazars, we start from the $\sim$1.4 million quasars included in the Million Quasar Catalog (MQC; Flesch, 2015). We select all $z>3.1$ sources and retain only radio-loud (RL) quasars with $R>10$, where $R$ is the ratio of the rest-frame 5 GHz (extrapolated from 1.4 GHz, considering a flat radio spectrum) to optical $B$ band flux density (Kellermann et al., 1989), assuming an optical spectral index $-$0.5 ($F_{\nu}\propto\nu^{\alpha}$). These 1103 objects111We exclude all the objects with a photometric redshift in MQC. represent our parent sample and we analyze LAT data for all of them according to the following procedure.

For each object we use $\sim$92 months (from 2008 August 5 to 2016 April 1) of Fermi Pass 8 source class photons and analyze them following the standard data reduction procedure222http://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/, but with a few modifications as mentioned below. We define a region of interest (ROI) of 15*∘* radius centered on each quasar and define a sky model that includes all $\gamma$-ray sources detected in the third Fermi-LAT catalog (3FGL; Acero et al., 2015), within the ROI and the isotropic and Galactic diffuse emission models (Acero et al., 2016). The parameters of all sources within the ROI and power-law normalization factors of the diffuse models are optimized so that the sky model reproduces the data as best as possible. This is done via a binned likelihood method and we measure the significance of the detection by means of the maximum likelihood test statistic TS= 2$\Delta\log\mathcal{L}$, where $\mathcal{L}$ represents the likelihood function, between models with and without a point source at the position of the quasar. The targets of interest are modeled with a simple power-law model leaving the prefactor and the photon index free to vary during the likelihood fitting. We consider a source to be significantly detected if TS$>25$ (4.2$\sigma$; Mattox et al., 1996).

Since we are dealing with faint sources, we modify the standard data reduction procedure as follows. We expand the energy range considered so it spans 60 MeV to 300 GeV. This allows the analysis to be more sensitive to spectrally soft $\gamma$-ray sources, i.e., blazars whose high-energy peak is shifted to lower energies ($\sim$1$-$10 MeV) as typical for high-$z$ blazars. Moreover, Pass 8 provides an increase in the acceptance at $<$100 MeV by up to 75 %, with respect to Pass 7, which translates into an improved sensitivity for high-$z$ blazars. Because the energy resolution becomes increasingly worse at low energies we enable the energy dispersion correction in the analysis for all sources except the empirical diffuse models. A novelty introduced by Pass 8 is the characterization of the photons in PSF (point-spread function) type events, which sub-classify the events into four quartiles by quality of angular reconstruction, with the lowest quartile (PSF0) and highest quartile (PSF3) having the worst and the best, respectively, direction reconstruction. In order to take full advantage of this improvement, we perform a component-wise data analysis for each PSF event type by considering the product of the likelihood function (sum of the logarithms) for the four PSF event types, using the SUMMED likelihood method of the Science Tools333http://fermi.gsfc.nasa.gov/ssc/data/analysis/software/.

As in the case of our target sources, there could be faint gamma-ray emitters present in the data but not in the 3FGL catalog. Thus, we rely on an iterative procedure to discover, localize, and insert such sources into the sky model. This is done by generating a residual TS map for each ROI. The spatial positions of unmodeled excesses with TS$\geq$25 are optimized and inserted into the sky model with power-law spectra. This procedure is iterated until no significant unmodeled emission remains in the ROI.

3 Results and Discussion

Our systematic search for significant $\gamma$-ray emitters among a large sample of RL quasars has led to the detection of five sources. A likelihood ratio (LR) method (see Ackermann et al., 2011, for details) associates the detected $\gamma$-ray sources and their radio counterparts (from NRAO VLA Sky Survey or NVSS; Condon et al., 1998) with high confidence (association probability $>$80%). By repeating the LAT analysis adopting 1000 random positions drawn from a randomized NVSS catalog444The NVSS catalog is randomized by mirroring the Galactic longitudes (glon) of the sources to 360${}^{\circ}-$glon., we found that the probability that any of the newly detected $\gamma$-ray sources are spurious is negligible. The basic information for these objects is presented in Table 1 where we also show the results of the LAT data analysis. As can be seen, all the objects are extremely RL and $\gamma$-ray luminous quasars. According to our analysis, NVSS J151002+570243 ($z=4.31$) is now the farthest known $\gamma$-ray emitting blazar555The blazar QSO J0906$+$6930 discovered at $z=5.48$ by Romani et al. (2004) was found to be spatially coincident with a 1.5 $\sigma$ EGRET fluctuation, but it has not been confirmed as a $\gamma$-ray source..In Figure 1, we show residual TS maps of five detected objects, along with their radio and optimized $\gamma$-ray positions.

In general, high-redshift blazars are brighter at hard X-rays than in the $\gamma$-ray band (e.g., Romani, 2006; Sbarrato et al., 2013), probably due to the shift of the blazar SED to lower frequencies. This could be due to the intrinsic shift of the high-energy peak to lower energies as the bolometric non-thermal luminosity increases (Fossati et al., 1998; Donato et al., 2001). Another possible reason for the shift could be the fact that the high-energy emission of such high-redshift blazars is powered via IC scattering off photons from the torus rather than the broad-line region (BLR), which also contributes to the lowering of the frequency of their SED peak (Sikora et al., 2002). Alternatively, the SED peaks can also shift to lower energies provided the emission region is within the dense BLR photon field. In this case, an efficient cooling of the emitting electrons will cause the lowering of the SED peaks (Ghisellini et al., 1998). The shift of the SED causes their $\gamma$-ray spectra to become steeper and to move slightly outside, or at the limit, of the Fermi-LAT band. Indeed, all the blazars discovered here exhibit steep $\gamma$-ray spectra ($\Gamma_{\gamma}>2.5$, see Table 1). This suggests their IC peak lie at MeV energies.

In Figure 2, we compare these newly detected distant objects with the blazars included in the third catalog of Fermi-LAT detected AGN (3LAC; Ackermann et al., 2015). As can be seen in the left panel, these sources occupy the region of high $\gamma$-ray luminosities ($L_{\gamma}>10^{47}$ erg s*-1*) and soft photon indices ($\Gamma_{\gamma}>2.5$), typical of powerful blazars. The right panel of Figure 2 compares the redshift distributions of these newly discovered high-$z$ blazars with that of the 3LAC blazars. Though the population of distant blazars is small, this work opens a window for the study of high-$z$ blazars and may have a major impact on constraining the various physical parameters associated with the blazar population (see, e.g., Paliya et al., 2016, for a relevant discussion).

In order to understand the broadband behavior of these high-$z$ blazars we look into the literature for multi-frequency information. Though there is a paucity of such data, we found a few noteworthy observations that support the blazar nature of these objects. NVSS J064632+445116 was predicted as a candidate $\gamma$-ray emitter by Healey et al. (2008), whereas NVSS J135406$-$020603 is included in the ROMA-BZCAT (Massaro et al., 2015). The quasar NVSS J212912$-$153841 is a hard X-ray spectrum luminous blazar and included in the 70 month Swift-Burst Alert Telescope catalog (Baumgartner et al., 2013). NVSS J151002+570243 is one of the best studied among all of the objects and exhibits an intense and hard X-ray spectrum (e.g., Mathur & Elvis, 1995; Moran & Helfand, 1997; Wu et al., 2013).

We generate broadband SEDs of the three objects that have archival X-ray observations and model them using a simple one zone synchrotron-IC emission approach prescribed in Ghisellini & Tavecchio (2009). In brief, the model assumes a spherical emission region located at a distance $R_{\rm diss}$ from the central engine and filled with a population of highly energetic electrons that follow a broken power-law distribution. In the presence of a tangled magnetic field, the electrons lose energy via synchrotron, synchrotron self Compton (SSC), and external Compton (EC) processes. For the latter, the seed photons originate from several external AGN components: photons directly emitted from the accretion disk (e.g., Dermer & Schlickeiser, 1993), from BLR (Sikora et al., 1994) and from the infrared torus (e.g., Błażejowski et al., 2000). The calculated jet powers and SED parameters are given in Table 2 and the modeled SEDs are shown in Figure 3.

In each of the three objects, the IR-UV emission is found to be dominated by an extremely luminous accretion disk ($L_{\rm disk}>10^{46}$ erg s*-1*). The X-ray spectra, on the other hand, are hard and the entire X-ray to $\gamma$-ray band of the SED can be explained by the IC scattering off the photons originating from the BLR. A strong accretion disk radiation implies a dense BLR photon field surrounding the jet, which in turn is observable in the form of broad optical emission lines. According to our SED modeling analysis, a large BLR radiative energy density indicates that most of the high-energy emission originates from the interaction of the BLR photons with the jet electrons. This suggests that the cooling of the electrons will be efficient, and accordingly the synchrotron emission will peak at low frequencies, which is supported by the modeling results. This indicates the location of the $\gamma$-ray emitting region to be inside the BLR. However, it should be noted that with the sparse available observations, it is not possible to tightly constrain the location of the emission region. The Compton dominance (ratio of the IC to synchrotron peak luminosities) of each of the three sources is also very large ($>20$, Table 2), a characteristic feature exhibited by powerful blazars. Other SED parameters are similar to those generally observed in high-$z$ blazars (e.g., Ghisellini et al., 2010). Though the data used here are mostly non-simultaneous, they indicate a typical state of the blazar rather than any period of specific activity. Also, the $\gamma$-ray photon statistics are not good enough to search for temporal variability. Overall, the Fermi-LAT detection and the available data confirm the blazar nature of the 5 high-$z$ RL quasars.

Powerful blazars are generally found to host massive black holes at their centers (e.g., Ghisellini et al., 2010; Paliya et al., 2016). It is, therefore, of great interest to determine the black hole mass of these $\gamma$-ray detected quasars. In Table 1, we report the black hole masses for each of the five sources using information that we could find/derive from the optical spectroscopic information available in the literature Torrealba et al. (2012); Alam et al. (2015). Furthermore, we also derive the masses by modeling the observed IR-UV emission for three objects (Figure 3) with a standard Shakura & Sunyaev (1973) accretion disk and the results are presented in Table 2. Both methods predict the existence of massive black holes ($\sim 10^{8-10}{\rm M}_{\odot}$) and match reasonably well within a factor of $\sim$three666One should keep in mind the large uncertainties ($\sim$0.4 dex) associated with the virial BH mass estimation methods (see, e.g., Vestergaard & Peterson, 2006; Shen et al., 2011)., except for NVSS J151002+570243, for which the modeling approach predicts a higher black hole mass. In particular, the object NVSS J212912$-$153841 hosts one of the most massive black holes, $\sim$7 $\times 10^{9}$  M*⊙*, ever found in $\gamma$-ray emitting blazars, confirmed both from optical spectroscopic and disk modeling approaches.

At redshifts between 3 and 4, the space density of black holes with $M_{\rm BH}>10^{9}\,{\rm M}_{\odot}$ hosted in jetted AGN is  50 Gpc*-3* (Sbarrato et al., 2015). This estimate is based on the luminosity function reported in Ajello et al. (2009), which, at those redshifts, relies only on five blazars. This work finds two more blazars hosting massive black holes ($M_{\rm BH}>10^{9}$ ${\rm M}_{\odot}$) in the same redshift range, considering black hole masses derived from the optical spectroscopic information. Adopting $\Gamma=13$, derived from our SED modeling, these two objects imply the presence of $\sim$675 (i.e., $2\times 2\Gamma^{2}$) similar systems, but with jets pointing in all directions, in the same redshift bin. The number density related to these two sources can be computed as $n=675/(V_{\rm MAX}\times f_{sky}\times f_{prob}\times f_{z})\approx 18$  Gpc*-3*, where $f_{sky}$=0.52 is the fraction of the sky covered by our parent sample, $f_{z}$=0.84 is the fraction of sources in the parent sample with a spectroscopic redshift, $f_{prob}$=0.66 is the fraction of $\gamma$-ray sources that are associated and $V_{\rm MAX}=134$ Gpc*-3* is the available volume777This has been computed taking into account the detection efficiency of the LAT. where these sources could have been detected (Schmidt, 1968). This brings the estimate of the space density of massive black holes hosted in jetted systems to 68${}^{+36}_{-24}$ Gpc*-3*. Already at redshift 4, this implies that there is a similar number of massive black holes hosted in radio-loud and radio-quiet systems and that, given their strong evolution, above that redshift most massive black holes might be hosted in radio-loud systems (Ghisellini et al., 2010; Volonteri et al., 2011). This clearly shows that the radio-loud phase may be a key ingredient for quick black hole growth in the early Universe. To this end, the detection of high-$z$ blazars becomes very important. Currently, the most promising approaches are, (1) lowering the energy threshold of the LAT, and (2) using NuSTAR. However, the optimal instrument would be a sensitive all-sky MeV telescope, e.g., e-ASTROGAM (Tatischeff et al., 2016) and AMEGO888https://pcos.gsfc.nasa.gov/physpag/probe/AMEGO_probe.pdf.

We are grateful to the referee for insightful comments. The Fermi-LAT Collaboration acknowledges support for LAT development, operation and data analysis from NASA and DOE (United States), CEA/Irfu and IN2P3/CNRS (France), ASI and INFN (Italy), MEXT, KEK, and JAXA (Japan), and the K.A. Wallenberg Foundation, the Swedish Research Council and the National Space Board (Sweden). Science analysis support in the operations phase from INAF (Italy) and CNES (France) is also gratefully acknowledged. Part of this work is based on archival data, software, or online services provided by the ASI Science Data Center (ASDC).