High-energy Emission Component, Population, and Contribution to the Extragalactic Gamma-Ray Background of Gamma-Ray-emitting Radio Galaxies

Yasushi Fukazawa; Hiroto Matake; Taishu Kayanoki; Yoshiyuki Inoue; Justin Finke

doi:10.3847/1538-4357/ac6acb

1. Introduction

The formation of supermassive black holes is one of the most intriguing questions in modern astrophysics. They are believed to coevolve with their host galaxies through feedback processes (Fabian 2012). Active galactic nucleus (AGN) activity drives feedback via relativistic jets or fast outflows. The evolutionary history of these AGN activities is the key to understanding the feedback history in the universe. The luminosity function (LF) of AGNs and their contribution to the cosmic X-ray and gamma-ray background radiation allow us to unveil the evolutionary history of AGNs. These two are complementary to each other. The former gives the differential history while the latter gives the integrated history.

The evolutionary history of blazars (AGNs with jets pointed toward Earth) and their contribution to the cosmic gamma-ray background (CGB) have been well studied using data from the Neil Gehrels Swift Observatory's Burst Alert Telescope (BAT) in hard X-rays and Fermi Large Area Telescope (LAT) in gamma-rays (see e.g., Ajello et al. 2009, 2012, 2014; Toda et al. 2020). However, although blazars are the brightest and most easily detected AGNs with jets owing to beaming, they represent <1% of the AGN jet population. The most numerous AGN jet population is misaligned AGNs (MAGNs), i.e., radio galaxies. No definition of MAGNs is given, but typically MAGNs are sources with jets with viewing angle $\gt \tfrac{1}{{\rm{\Gamma }}}$ , where Γ is a jet's Lorentz factor.

Studies of the evolutionary history of radio galaxies have been mainly conducted in the radio band (see, e.g., Willott et al. 2001), where they are easiest to detect. However, due to synchrotron self-absorption, the radio emission from the relativistic jet is dominated by the downstream component (Blandford & Königl 1979), reflecting the past activity of AGNs. High spatial resolution radio observations or gamma-ray observations can allow us to investigate the current jet activity.

Fermi-LAT has been monitoring the entire GeV gamma-ray sky since its launch in 2008. Although radio galaxies are several orders of magnitude fainter than blazars in the gamma-ray band, early Fermi-LAT observations reported a small number of radio galaxies (Abdo et al. 2010d). Based on those samples, past studies reported that radio galaxies contribute 10%–50% of the CGB using radio LFs (e.g., Inoue 2011; Di Mauro et al. 2014). Ajello et al. (2015) argued that, together, blazars and star-forming galaxies make up 100% of the CGB in the 0.1–800 GeV band. However, there is considerable uncertainty in the radio galaxies' contribution owing to a limited sample size. An accurate estimate will allow improved estimates of a possible component of the CGB from dark matter annihilation.

Here, the 4th Fermi-LAT Gamma-ray source catalog (4FGL-DR2; Ballet et al. 2020) contains 61 MAGNs based on its 10 yr survey. This number is significantly increased from the past gamma-ray catalogs: 11 MAGNs in the first catalog (1FGL; Abdo et al. 2010b), 12 in the second (2FGL; Nolan et al. 2012), and 21 in the third (3FGL; Acero et al. 2015).⁸ The larger statistical sample of gamma-ray-detected radio galaxies in 4FGL-DR2 allows us to constrain their gamma-ray LF (GLF) and their CGB contribution using only gamma-ray data. Others have tried to estimate the radio galaxy contribution to the CGB using radio LFs and ratios of gamma-ray to radio flux (e.g., Inoue 2011; Di Mauro et al. 2014).

We can do further studies of other features of AGN jets with this statistical sample, such as the relationship between blazars and radio galaxies and that between gamma-ray-selected and radio-selected radio galaxies. In the radio band, the radio galaxies are often classified as Fanaroff & Riley (1974) type I (FR I) and type II (FR II), where FR I galaxies are brighter in the middle and fainter at the edges and FR II galaxies are fainter in the middle and brighter at the edges. Fanaroff & Riley (1974) and subsequent work (e.g., Ledlow & Owen 1996) indicated that FR I galaxies are overall fainter in radio than FR II galaxies (although see Mingo et al. 2019). FR I radio galaxies are widely considered the parent population of BL Lac type blazars (Padovani & Urry 1990), and FR II radio galaxies are considered the parent population of flat-spectrum radio quasar (FSRQ) type blazars. FSRQs are blazars with high luminosities and strong broad emission lines in their optical spectra, while BL Lacs are blazars with lower luminosities and weak or absent broad emission lines. Recently, radio galaxies with core-dominant radio emission have been found. These cannot be classified as FR I or FR II and have been dubbed FR 0 (Baldi & Capetti 2010; Ghisellini 2011). The first Fermi Gamma-ray Catalog contained seven FR I galaxies and four FR II galaxies (Abdo et al. 2010d), and thus FR II galaxies did not seem to be as bright in the GeV gamma-ray band. The MAGN sample of the 4FGL-DR2 can elucidate these questions by constraining the GLF and comparing it to the GLF of blazars and to the radio LF of radio galaxies.

In addition, we can combine gamma-ray data from Fermi-LAT with X-ray observations to investigate the spectral energy distribution (SED) shape of radio galaxies from X-rays to gamma-rays. In the case of blazars, the SEDs have a controversial luminosity dependence, the so-called blazar sequence (Fossati et al. 1998; Ghisellini et al. 1998, 2017; Kubo et al. 1998; Finke 2013), while the SEDs of radio galaxies have not been well studied. Constraining the SED from X-rays to gamma-rays is also essential for estimating the contribution of radio galaxies to the CGB spectrum.

We construct the GLF and the X-ray and gamma-ray SED of radio galaxies and establish their contribution to the CGB. In Section 2 we summarize general properties of our sample of gamma-ray-emitting radio galaxies, including their X-ray properties. X-ray and gamma-ray SEDs are presented in Section 3. The GLF and the contribution to the CGB are described in Sections 4 and 5, respectively. The relations to blazars and radio-selected radio galaxies are discussed in Section 6. Conclusions are presented in Section 7. Throughout this paper, we assume cosmological parameters of Ω_Λ = 0.7, Ω_M = 0.3, and H₀ = 70 km s⁻¹ Mpc⁻¹.

2. Sample of Radio Galaxies

2.1. Gamma-Ray Data

Fermi-LAT has been surveying the gamma-ray sky since 2008 (Atwood et al. 2009). The latest catalog, the 4FGL-DR2, is based on its 10 yr of survey data, from 2008 August 4 to 2018 August 2. The previous catalog, the Fermi-LAT 4th catalog (4FGL), was based on 8 yr of data and was described in detail in Abdollahi et al. (2020). The 4FGL-DR2 has 5788 gamma-ray sources detected at >4σ in the energy range from 50 MeV to 1 TeV, including 61 MAGNs. We quote the GeV gamma-ray photon index, gamma-ray flux in 0.1–100 GeV, six-band flux, redshift, and classification from the 4FGL-DR2. Note that the detection threshold in 4FGL-DR2 is the same as that in the previous catalogs from 1FGL to 4FGL.

The 61 4FGL-DR2 MAGNs are classified into four types: 43 radio galaxies (RDG), 5 compact steep-spectrum radio sources (CSS), 2 steep-spectrum radio quasars (SSRQ), and 11 AGNs, as summarized in Table 1. Hereafter we call this sample the 4FGL-DR2 MAGNs. Two new radio galaxies (RDG) were added in the 4FGL-DR2 that were not in the 4FGL. CSSs are powerful radio sources associated with an AGN, characterized by their small radio size and steep slope of the radio spectrum, with a peak around 100 MHz. SSRQs are high-luminosity AGNs with lobe-dominated radio emission and radio spectral slopes >0.5 at frequencies of several GHz. AGNs as defined in the 4FGL-DR2 are nonblazar AGNs whose existing data do not allow an unambiguous determination of their AGN types. In this paper, we treat them as radio galaxies. For many of them, we give their radio morphological FR type from the literature in Table 1. Although typically FR I galaxies are associated with lower luminosities and FR II galaxies with higher luminosities, recent work has indicated that this is not always the case and that many low- and high-luminosity objects can be found with both morphological types (Mingo et al. 2019).

Table 1. Our Sample 4FGL-DR2 Radio Galaxies

No.	4FGL Name	Galaxy Name	Class ^a	z	${\mathrm{log}}_{10}{L}_{1.4}$ ^b	RG ^c	Flag ^d
1	J0009.7−3217	IC 1531	rdg	0.025	−0.08	I^B18	001110+
2	J0013.6+4051	4C +40.01	agn	0.255	2.50	⋯	011100+
3	J0038.7−0204	3C 17	rdg	0.220	2.39	I/II^ZB95	111000+
4	J0057.7+3023	NGC 315	rdg	0.016	−0.44	II^ICS99	011110+
5	J0153.4+7114	TXS 0149+710	rdg	0.022	−0.49	I^L01	011111−
6	J0237.7+0206	PKS 0235+017	rdg	0.022	−1.07	I^CB88I	000110+
7	J0308.4+0407	NGC 1218	rdg	0.029	1.13	I^ZB95	011111+
8	J0312.9+4119	B3 0309+411B	rdg	0.136	1.52	II^ICS99	001110−
9	J0316.8+4120	IC 310	RDG	0.019	−0.86	H/T^M93	001111−
10	J0319.8+4130	NGC 1275	RDG	0.018	1.20	I^ZB95	111111−
11	J0322.6−3712e	Fornax A	RDG	0.006	−2.30	I^ZB95	011111+
12	J0334.3+3920	4C +39.12	rdg	0.021	−0.97	0^RBC20	001111−
13	J0418.2+3807	3C 111	rdg	0.049	1.64	II^OL89	111100−
14	J0433.0+0522	3C 120	RDG	0.033	−0.16	I^OL89	111100+
15	J0519.6−4544	Pictor A	rdg	0.035	1.25	II^ZB95	111110+
16	J0521.2+1637	3C 138	css	0.759	4.34	C^LRL83	001110−
17	J0522.9−3628	PKS 0521−36	AGN	0.056	1.37	I/II^B16	111111+
18	J0627.0−3529	PKS 0625−35	rdg	0.055	0.88	I^W04	111111−
19	J0708.9+4839	NGC 2329	rdg	0.019	−0.24	I^D21	001110+
20	J0758.7+3746	NGC 2484	rdg	0.043	1.05	I^G94	001100+
21	J0840.8+1317	3C 207	ssrq	0.681	3.72	II^LRL83	111100+
22	J0858.1+1405	3C 212	ssrq	1.048	4.17	II^LRL83	011000+
23	J0910.0+4257	3C 216	css	0.670	3.87	II^P06	111000+
24	J0931.9+6737	NGC 2892	rdg	0.023	−0.24	I^J82I	111110+
25	J0958.3−2656	NGC 3078	rdg	0.009	−1.29	I^WH84I	001100+
26	J1012.7+4228	B3 1009+427	agn	0.365	1.56	II^K18	001111+
27	J1116.6+2915	B2 1113+29	rdg	0.047	1.00	II^ZB95	000010+
28	J1118.2−0415	PMN J1118−0413	agn	0.000	⋯	⋯	111100+
29	J1144.9+1937	3C 264	rdg	0.022	0.81	I^LRL83	001111+
30	J1149.0+5924	NGC 3894	rdg	0.011	−0.04	I^X95I	001110+
31	J1219.6+0550	NGC 4261	rdg	0.007	0.31	I^ZB95	001110+
32	J1230.8+1223	M87	rdg	0.004	−1.70	I^ZB95	111111+
33	J1236.9−7232	PKS 1234−723	rdg	0.024	0.18	I^FJ02	011100−
34	J1306.3+1113	TXS 1303+114	rdg	0.086	0.88	I^CMB17a	001100+
35	J1306.7−2148	PKS 1304−215	rdg	0.126	1.11	⋯	111110+
36	J1325.5−4300	Cen A	RDG	0.002	−0.43	T^ZB95	111111−
37	J1331.0+3032	3C 286	css	0.850	4.70	C^LRL83	111100+
38	J1346.3−6026	Cen B	rdg	0.013	0.47	I^JLM01	111110−
39	J1356.2−1726	PKS B1353−171	agn	0.075	0.39	⋯	001100+
40	J1443.1+5201	3C 303	rdg	0.141	2.13	II^LRL83	001110+
41	J1449.5+2746	B2 1447+27	rdg	0.031	−0.88	⋯	000110+
42	J1449.7−0910	1RXS J1449−0910^e	agn	0.000	⋯	⋯	000110+
43	J1459.0+7140	3C 309.1	css	0.910	3.14	C^LRL83	111100+
44	J1516.5+0015	PKS 1514+00	rdg	0.052	0.70	II^CMB17b	111100+
45	J1518.6+0614	TXS 1516+064	rdg	0.102	1.12	I^CMB17a	000111+
46	J1521.1+0421	PKS B1518+045	rdg	0.052	0.42	I^CMB17a	001110+
47	J1543.6+0452	CGCG 050−083	agn	0.040	−0.43	⋯	111111+
48	J1630.6+8234	NGC 6251	rdg	0.024	0.02	I/II^OL89	111110+
49	J1724.2−6501	NGC 6328	rdg	0.014	−0.58	GPS/CSO^T97	111100−
50	J1824.7−3243	PKS 1821−327	agn	0.355	3.41	⋯	011100−
51	J1829.5+4845	3C 380	css	0.695	4.50	C^LRL83	111110+
52	J1843.4−4835	PKS 1839−48	rdg	0.111	2.21	I^ZB95	001010−
53	J2114.8+2026	TXS 2112+202	agn	0.000	⋯	⋯	001110−
54	J2156.0−6942	PKS 2153−69	rdg	0.028	1.75	II^F98	011000+
55	J2227.9−3031	PKS 2225−308	rdg	0.056	−0.08	I?^ZB95	001010+
56	J2302.8−1841	PKS 2300−18	rdg	0.129	1.58	II?^ZB95	001110+
57	J2326.9−0201	PKS 2324−02	rdg	0.188	1.16	⋯	111010+
58	J2329.7−2118	PKS 2327−215	rdg	0.031	0.15	⋯	111110+
59	J2334.9−2346	PKS 2331−240	agn	0.048	0.63	II^B16	111100+
60	J2338.1+0325	PKS 2335+03	agn	0.270	2.50	II^H17	011110+
61	J2341.8−2917	PKS 2338−295	rdg	0.052	−0.58	⋯	001100

Notes.

^aSource class in 4FGL-DR2. Uppercase letters represent objects identified as the same objects found in other wavelengths, and lowercase letters represent objects positionally associated with objects found in other wavelengths.^bRadio luminosity in units of 10²⁴ W Hz⁻¹ at 1.4 GHz from Angioni (2020) and NED. "0" represents that no data are available.^cRadio morphology type. I: FR I; II: FR II; C: core; J: jet; T: transition; H/T: head/tail. Superscript represents references for morphology. B16: Bassani et al. (2016); B18: Bassi et al. (2018); CMB17a: Capetti et al. (2017a); CMB17b: Capetti et al. (2017b); CB88I: Condon & Broderick (1988); D21: Das et al. (2021); F98: Fosbury et al. (1998); G94: Giovannini et al. (1994); H17: Hernández-García et al. (2017); ICS99: Ishwara-Chandra & Saikia (1999); J82I: Jenkins (1982); JLM01: Jones et al. (2001); K18: Kuźmicz et al. (2018); LRL83: Laing et al. (1983); L01: Lara et al. (2001); LJ02: Lloyd & Jones (2002); M93: Mack et al. (1993); OL89: Owen & Laing (1989); P06: Punsly (2006); RBC20: Rulten et al. (2020); T97: Tingay et al. (1997); W04: Wills et al. (2004); WH84I: Wrobel & Heeschen (1984); X95I: Xu et al. (1995); ZB95: Zirbel & Baum (1995). For CB88I, J82I, WH84I, and X95I, we determined FR I or FR II based on the radio image in the reference paper.^dDetection (1) or nondetection (0) at each of six bands (0.1–0.3 GeV, 0.3–1.0 GeV, 1.0–3.0 GeV, 3.0–10.0 GeV, 10.0–30.0 GeV, 30.0–300.0 GeV) in 4FGL-DR2. The right end flag is + for ∣b∣ ≥ 20° or 0 for ∣b∣ < 20°.^eThe exact catalog name is 1RXS J144942.2–091018.

Download table as: ASCIITypeset images: 1 2

Figure 1 shows a relationship between GeV gamma-ray and radio luminosity for our sample. We also plot radio galaxies not detected by Fermi-LAT from the radio flux-limited sample (Mingo et al. 2014), where we give upper limits on the GeV gamma-ray luminosity for an upper limit flux of 10⁻¹² erg s⁻¹ cm⁻² (0.1–300 GeV). Interestingly, radio galaxies in the radio flux-limited sample are dominated by high-luminosity objects (most are FR II), while our GeV gamma-ray flux-limited sample contains many low-luminosity ones. We will discuss this issue in Section 6.2.

Table 1 also summarizes the availability of the six spectral bands in 4FGL-DR2 and the flag for whether the source latitude b is ∣b∣ < 20° or not, since the Galactic diffuse emission can bring significant systematic uncertainties. We use these pieces of information in the LF modeling described later. Some objects have a detection flag at one or two energy bands since some energy bands (1–3 GeV, 3–10 GeV, and 10–30 GeV) have better sensitivity than others.

2.2. X-Ray Data

X-ray spectral information, combined with GeV gamma-rays, is also essential to investigating the high-energy emission from MAGNs. We searched available X-ray observational data on the 4FGL-DR2 MAGNs using the XMM-Newton/EPIC (Turner et al. 2001), Chandra/ACIS (Weisskopf et al. 2000), and Swift/XRT (Gehrels et al. 2004; Burrows et al. 2005) archival lists. When available, we chose the data having the highest photon statistics for each MAGN. Since multiple X-ray satellite data sets are available for some MAGNs, we selected the data in the order of priority: first XMM-Newton/EPIC, then Chandra/ACIS, and finally Swift/XRT. Then, we analyzed 20 data sets from XMM-Newton/EPIC, 15 data sets from Chandra/ACIS, and 9 data sets from Swift/XRT. During most of these observations, the instruments were pointed directly at the 4FGL-DR2 MAGNs; however, some were observed off-axis. For the seven gamma-ray-bright radio galaxies, we refer to the results of Suzaku data analysis in Fukazawa et al. (2015), instead of analyzing other satellite data. As a result, we obtained X-ray data on 51 galaxies. AGNs often show time variability in the X-ray band, but we confirmed that the variability amplitude is at most a factor of 2–3 by using Swift/XRT data for about 20 objects. Our results presented here are not affected significantly by such time variability.

XMM-Newton/EPIC data were analyzed with SAS version 15.0.0 in the standard way; reprocessing and elimination of background flaring time region were done before making spectra. Photons within 60'' of the object were extracted for spectral analysis, and background spectra were extracted in the annular region with radii from 500'' to 550''. Chandra/ACIS data were analyzed with CIAO version 4.11 in the standard way, using chandra_repro and specextract for reprocessing and extraction of spectra, respectively. Photons within 2'' of the object were extracted for spectral analysis. For some objects observed in the off-axis position, we set a larger extraction radius to collect enough photons. Swift/XRT data were retrieved from the UK Swift Science Data Center and analyzed with HEASoft version 6.19 in the standard way, using xselect. Photons within 30'' of the object were extracted for spectral analysis, and background spectra were extracted in the annular region with radii from 210'' to 250''.

The X-ray spectra obtained above were fitted with XSPEC version 12.9.0o. We adopt a spectral model of a single power law with the Galactic interstellar medium absorption, with hydrogen column density from Dickey & Lockman (1990). Thermal emission components also appear in the soft X-ray band in some MAGNs (IC 1531, NGC 315, NGC 1316, NGC 2484, and NGC 4261), which originate in the hot interstellar or intragroup medium. In such cases, the thermal plasma model apec (Smith et al. 2001) was added in the analysis to improve the goodness of fit. Two objects (NGC 4261 and NGC 6251) have an excess absorption column density well above the Milky Way value, so we added an additional absorption component. Table 2 summarizes the resulting X-ray photon indices and fluxes of our sample galaxies. Detailed results of X-ray spectral analysis will be presented in a separate paper (Kayanoki & Fukazawa 2022).

Table 2. X-Ray and Gamma-Ray Properties of 4FGL-DR2 Radio Galaxies

No. ^a	Data ^b	α_G ^b	α_X ^c	F_G(L_γ) ^d	F_X(L_X) ^e	ν_HE ^f	L_HE ^g
1	xmm	2.15 ± 0.13	1.97 ± 0.03	−11.75(42.42)	−12.26(41.91)	21.9	42.1
	0202190301
2	chandra	2.19 ± 0.15	1.89 ± 0.07	−11.73(44.57)	−11.62(44.68)	21.4	44.7
	5669
3	0	2.89 ± 0.10	0	−11.35(44.81)	0.00(0.00)	19.6	45.5
4	xmm	2.37 ± 0.11	1.16 ± 0.02	−11.45(42.34)	−11.91(41.88)	21.6	42.0
	0305290201
5	swift	1.89 ± 0.11	1.96 ± 0.23	−11.47(42.57)	−11.76(42.28)	24.9	41.8
	00010149002
6	rass	2.15 ± 0.17	2.0	−11.82(42.21)	−13.15(40.87)	25.0	41.5
7	f15	1.97 ± 0.05	2.32 ± 0.04	−11.05(43.23)	−11.92(42.36)	24.2	42.5
8	xmm	2.69 ± 0.34	1.91 ± 0.02	−11.39(44.31)	−11.35(44.34)	21.3	44.3
	0306680301
9	xmm	1.85 ± 0.15	2.44 ± 0.01	−11.50(42.41)	−11.39(42.52)	25.0	41.7
	0151560101
10	f15	2.11 ± 0.01	1.73 ± 0.03	−9.48(44.39)	−10.92(42.95)	22.6	43.8
11	xmm	2.07 ± 0.06	1.69 ± 0.03	−11.21(41.68)	−12.22(40.66)	22.5	41.1
	0502070201
12	chandra	1.81 ± 0.12	2.28 ± 0.07	−11.61(42.37)	−12.04(41.95)	25.0	41.7
	857
13	f15	2.74 ± 0.06	1.65 ± 0.02	−10.84(43.91)	−10.26(44.48)	20.0	44.2
14	f15	2.74 ± 0.04	1.75 ± 0.03	−10.84(43.56)	−10.33(44.07)	20.2	44.0
15	xmm	2.43 ± 0.12	1.77 ± 0.00	−11.38(43.07)	−10.77(43.68)	20.4	43.3
	0206390101
16	chandra	2.23 ± 0.14	1.51 ± 0.14	−11.55(45.88)	−11.49(45.94)	21.5	45.8
	14996
17	xmm	2.46 ± 0.01	1.80 ± 0.00	−10.26(44.61)	−10.89(43.98)	21.6	44.3
	0065760201
18	f15	1.92 ± 0.03	2.25 ± 0.02	−10.82(44.04)	−11.09(43.77)	24.5	43.3
19	chandra	1.77 ± 0.15	2.61 ± 0.36	−11.85(42.08)	−11.89(42.03)	25.0	41.5
	5900
20	xmm	2.20 ± 0.16	1.98 ± 0.03	−11.80(42.83)	−12.42(42.21)	22.3	42.6
	0602390101
21	xmm	2.49 ± 0.11	1.59 ± 0.01	−11.54(45.76)	−11.53(45.78)	21.2	45.5
	0147670301
22	chandra	2.58 ± 0.14	1.37 ± 0.02	−11.62(46.15)	−11.64(46.13)	21.2	45.9
	434
23	chandra	2.53 ± 0.10	1.63 ± 0.13	−11.47(45.82)	−11.86(45.43)	21.7	45.4
	15002
24	chandra	2.28 ± 0.06	3.77 ± 0.52	−11.27(42.80)	−13.18(40.89)	22.8	42.1
	18038
25	chandra	2.12 ± 0.16	2.56 ± 0.20	−11.84(41.37)	−12.19(41.03)	23.0	40.8
	5902
26	swift	1.80 ± 0.09	2.11 ± 0.24	−11.47(45.19)	−11.53(45.13)	25.0	44.5
	00085352004
27	xmm	1.39 ± 0.26	3.03 ± 0.18	−12.29(42.43)	−13.31(41.40)	23.9	42.5
	0550270101
28	swift	2.64 ± 0.07	1.80 ± 0.35	−11.17(0.00)	−11.83(0.00)	0.0	0.0
	00083516001
29	xmm	2.00 ± 0.09	2.26 ± 0.02	−11.46(42.57)	−11.46(42.57)	24.3	41.8
	0602200301
30	chandra	2.18 ± 0.12	0.02 ± 0.10	−11.61(41.82)	−12.54(40.89)	23.6	41.0
	10389
31	xmm	2.08 ± 0.15	0.34 ± 0.05	−11.69(41.39)	−12.11(40.97)	22.1	41.1
	0502120101
32	chandra	2.06 ± 0.03	2.19 ± 0.03	−10.75(41.85)	−11.74(40.87)	23.7	41.1
	18838
33	rass	2.33 ± 0.12	2.0	−11.51(42.60)	−12.40(41.71)	23.0	41.8
34	0	1.98 ± 0.17	0	−11.81(43.45)	0.00(0.00)	24.5	42.7
35	swift	2.18 ± 0.08	1.37 ± 0.67	−11.27(44.35)	−12.25(43.37)	23.5	43.6
	00032810001
36	f15	2.64 ± 0.02	1.73 ± 0.03	−10.20(41.06)	−8.21(43.05)	19.7	41.5
37	chandra	2.39 ± 0.12	2.15 ± 0.17	−11.66(45.89)	−12.30(45.24)	22.7	45.2
	15006
38	xmm	2.40 ± 0.04	1.14 ± 0.01	−10.61(42.96)	−11.12(42.46)	22.3	42.5
	0092140101
39	0	2.02 ± 0.14	0	−11.97(43.17)	0.00(0.00)	20.3	44.6
40	nustar	2.05 ± 0.12	1.77 ± 0.05	−11.77(43.96)	−11.48(44.25)	21.1	44.2
	60463048002
41	swift	1.54 ± 0.17	1.70 ± 0.32	−12.08(42.25)	−12.37(41.97)	22.4	42.3
	00040619002
42	swift	2.07 ± 0.18	2.20 ± 0.19	−11.80(0.00)	−12.12(0.00)	0.0	0.0
	00087620002
43	chandra	2.50 ± 0.06	1.51 ± 0.02	−11.36(46.25)	−11.53(46.08)	21.4	45.9
	3105
44	xmm	2.54 ± 0.10	1.78 ± 0.01	−11.35(43.46)	−11.54(43.27)	21.1	43.3
	0103860601
45	xmm	1.75 ± 0.17	2.13 ± 0.01	−11.86(43.57)	−11.03(44.39)	25.0	42.9
	0018741001
46	rass	2.04 ± 0.15	2.0	−11.83(42.98)	−13.10(41.72)	23.8	42.3
47	rass	1.87 ± 0.07	2.0	−11.27(43.31)	−12.22(42.35)	24.8	42.6
48	f15	2.37 ± 0.03	1.82 ± 0.05	−10.88(43.24)	−11.59(42.53)	21.8	42.8
49	xmm	2.50 ± 0.15	1.65 ± 0.02	−11.61(42.06)	−12.12(41.55)	21.9	41.5
	0804520301
50	xmm	2.28 ± 0.11	1.73 ± 0.01	−11.59(45.04)	−11.09(45.54)	21.0	45.3
	0650591401
51	swift	2.42 ± 0.03	1.60 ± 0.05	−10.75(46.58)	−11.21(46.12)	21.8	46.1
	00081221001
52	chandra	2.03 ± 0.16	1.41 ± 0.14	−11.79(43.71)	−12.87(42.63)	23.1	43.2
	10321
53	rass	2.11 ± 0.18	2.0	−11.75(0.00)	−12.70(0.00)	0.0	0.0
54	xmm	2.87 ± 0.11	1.83 ± 0.00	−11.44(42.83)	−10.96(43.30)	20.2	43.1
	0152670101
55	chandra	1.98 ± 0.18	0	−11.91(42.96)	0.00(0.00)	24.3	42.3
	5798
56	swift	2.18 ± 0.12	1.59 ± 0.06	−11.61(44.03)	−11.05(44.59)	20.6	44.3
	00031729002
57	xmm	2.51 ± 0.14	1.75 ± 0.06	−11.46(44.54)	−11.91(44.09)	21.7	44.1
	0405860101
58	swift	2.45 ± 0.13	2.0	−11.33(43.01)	0.00(0.00)	22.6	42.3
	00090972001
59	xmm	2.50 ± 0.09	1.74 ± 0.00	−11.35(43.38)	−10.79(43.94)	20.4	43.6
	0760990101
60	0	2.37 ± 0.14	0	−11.62(44.73)	0.00(0.00)	22.9	43.9
61	xmmslew	2.24 ± 0.14	2.0	−11.71(43.09)	−12.10(42.71)	24.1	42.3

Notes.

^aObject number in Table 1.^bSatellites for X-ray data analyzed in this paper. "f15," "rass," or "xmmslwe" represents that Fukazawa et al. (2015), ROSAT All-Sky Survey Catalog, or XMM-Newton Slew Survey Catalog is referred to, respectively. The number in the second line is an observation ID.^cGamma-ray photon index from 4FGL-DR2.^dX-ray photon index.^eLogarithmic gamma-ray flux in units of erg s⁻¹ cm⁻² in 4FGL-DR2. Values in parentheses are logarithmic gamma-ray luminosity in units of erg s⁻¹ in 0.1–100 GeV from 4FGL-DR2.^fLogarithmic X-ray flux in units of erg s⁻¹ cm⁻² in 2–10 keV. Values in parentheses are logarithmic X-ray luminosity in units of erg s⁻¹.^gLogarithmic peak frequency in the SED from X-ray to gamma-ray band in units of Hz.^hLogarithmic peak luminosity in the SED from X-ray to gamma-ray band in units of erg s⁻¹. In the table, "0" represents that no data are available or a value cannot be estimated.

Download table as: ASCIITypeset images: 1 2

Pileup effects were nonnegligible for five Chandra/ACIS-selected MAGNs: NGC 2329, 3C 138, 3C 303, 3C 380, and NGC 3078. For these data, we fitted the spectra with Sherpa by including the pileup model jdpileup to obtain the spectral parameters. The pileup fraction is found to be 0.05–0.15 for these MAGNs, except for 3C 303. Since the Chandra data of 3C 303 were affected by significant pileup, instead we analyzed the NuStar (Harrison et al. 2013) data for this object.

For MAGNs, which are not significantly detected in the above X-ray data or have no X-ray pointing observations, we used the X-ray fluxes from the ROSAT All-Sky Survey (RASS; Voges et al. 1999) and XMM-Newton Slew Survey (XSS; Saxton et al. 2008). Their fluxes are also summarized in Table 2. The fluxes of five MAGNs are available in RASS, and one is available in XSS. As a result, X-ray information is available for 57 MAGNs. We set an upper limit on X-ray flux to 3 × 10⁻¹³ erg s cm⁻² for the other five galaxies, based on the RASS survey limit.

The left panel of Figure 2 shows the relation between the GeV gamma-ray luminosity and X-ray luminosity for the 4FGL-DR2 MAGNs. The X-ray luminosity L_X seems to correlate with the GeV gamma-ray luminosity L_GeV. However, this correlation could be due to the mutual dependence on redshift and a selection bias that only gamma-ray-emitting radio galaxies are plotted. We will discuss this issue in Sections 6.2 and 6.4. FR I galaxies are located at the lower-luminosity regime, while CSSs and SSRQs are located at the highest-luminosity end. FR II galaxies have luminosities located in the middle, but no clear separation with FR I galaxies is seen in this plot, regardless of a clear separation in radio luminosity. Note that the one outlier with low gamma-ray luminosity and high X-ray luminosity is Cen A.

The right panel of Figure 2 shows the photon index relation between the GeV gamma-ray band and the X-ray band. MAGNs can be divided into two classes: one class has a soft GeV gamma-ray photon index of ≳2 with a hard X-ray photon index of ≲2, and another has the opposite. The former group tends to contain higher-luminosity galaxies such as CSSs and SSRQs, while the latter tends to contain lower-luminosity galaxies. This relation between X-ray and GeV gamma-ray photon indices is similar to that of blazars (Sambruna et al. 2010). There are two outliers. One is B2 1447+27, which has a hard X-ray photon index (1.70) and a hard gamma-ray photon index (1.54). The quality of the Swift/XRT X-ray spectrum is not good, and we cannot rule out that this spectrum is affected by excess absorption or other emission components. The other one is NGC 2894, which has a very soft X-ray photon index (3.77) and a soft gamma-ray photon index (2.28). This Chandra/ACIS X-ray spectrum is also not of good quality, and some absorption or additional components might affect the power-law model parameters.

3. Spectral Energy Distribution from X-Ray to GeV Gamma-Ray

We can derive the SED from the X-ray to gamma-ray band by simultaneously fitting the X-ray spectral data and the six flux band measurements from Fermi-LAT. We adopt the following polynomial function so as to represent a flatter slope in the X-ray band and a steep drop in the highest-energy regime:

$\begin{eqnarray}\begin{array}{rcl}\mathrm{log}\nu F(\nu ) & = & A{\left(\mathrm{log}\nu -\mathrm{log}{\nu }_{0}\right)}^{2}\left(2\mathrm{log}\nu \right.\\ & & \left.+\,\mathrm{log}{\nu }_{0}-3\mathrm{log}{\nu }_{1}\right)+B,\end{array}\end{eqnarray} \tag{ 1 }$

where F and ν are the energy flux and frequency, respectively, and ν₀, ν₁, A, and B are fitting parameters. The parameter A is negative, and ν₀ and ν₁ are limited to be in the range of 10¹⁶–10¹⁷ Hz and > 10¹⁷ Hz, respectively. This function has a flux peak at ν₁, and the differential coefficient becomes zero at ν₀. Four free parameters were obtained by maximizing a log-likelihood $\mathrm{ln}L=-{\sum }_{i}{\left(F({\nu }_{i})-{F}_{i}\right)}^{2}/\left(2\delta {F}_{i}^{2}\right)$ with the Markov Chain Monte Carlo (MCMC) technique, where ν_i, F_i, and δ F_i are the measured frequency, flux, and flux error, respectively.

When the X-ray photon index is larger than 2, the X-ray emission could be synchrotron radiation (Fukazawa et al. 2015), and thus it is not appropriate to describe both X-ray and GeV gamma-ray with this function. In such a case, −1 or 10⁹ is added to the log-likelihood when the model is smaller or larger than the data flux, respectively, so that the model does not exceed the data; in other words, we use X-ray data as an upper limit for the high-energy emission component. When there are no available X-ray data, we use an upper limit flux of 10⁻¹¹ erg s⁻¹ cm⁻² at 1 keV and fit in the same way as above so that the model does not exceed this flux. If only RASS or XSS data are available, we treat that flux as an upper flux at 1 keV, as above.

The left panel of Figure 3 shows the SED curves obtained from X-rays to gamma-rays for our sample of MAGNs. SED curves seem to depend on the luminosity, which is clearly seen as a dependence on a gamma-ray-to-X-ray luminosity ratio on the SED peak luminosity shown in the right panel of Figure 3. Table 2 summarizes the estimated peak frequencies and luminosities of the high-energy components of MAGNs. Note that the high-energy component in the X-ray band has contributions not only from jet emission but also from coronal emission. Synchrotron peak frequencies and luminosities from jet emission are also interesting and have been studied in the past (Fukazawa et al. 2015; Keenan et al. 2021), and thus we quote these values from 4LAC (Lott et al. 2020) for our MAGNs. For comparison, we also quote the values of blazars. Synchrotron peak values are taken from 4LAC, while IC peak values are from Abdo et al. (2010a), where the peak was estimated by fitting the available SED data with polynomial functions. Note that these peak values are obtained by just fitting SEDs by polynomial functions without modeling non-jet components in our study. For the 4LAC paper (Ajello et al. 2020), if non-jet components are clearly seen in the SED, the data affected by such components are ignored in the fitting. However, for MAGNs, non-jet components are not clearly seen in the SED. Therefore, these peak frequencies and luminosities could be affected by non-jet emission components.

Figure 4 shows the relation between peak frequency and peak luminosity for the synchrotron and high-energy components. It can be seen from the top left panel that the synchrotron peak frequency of radio galaxies has a wide distribution, as for blazars, from infrared to X-rays bands but tends to concentrate toward a lower frequency. However, we note that the infrared, optical, and X-ray bands could be significantly contaminated by the emission from host galaxies and AGN disks/coronae, and thus the SED peak position could be affected. There is a weak correlation between the peak frequencies of the synchrotron and high-energy components, but a scatter could again be attributed to the contamination from other emission components. The highest synchrotron peak frequency around 10^17.5 Hz is TXS 1516+064, which has a flat GeV gamma-ray photon index of 1.75.

The top right panel shows that the peak luminosities of the synchrotron and high-energy components have a tight correlation in radio galaxies and blazars. Note that the peak luminosity of the high-energy component of radio galaxies is somewhat lower than that of the synchrotron component; this is in contrast to the Compton dominance of high-luminosity blazars, especially FSRQs. Contamination of other emission components to the synchrotron peak luminosity makes this trend less certain. These properties will be discussed again in Section 6.2.

The bottom left panel shows that radio galaxies do not seem to follow the anticorrelation between synchrotron peak frequency and luminosity, sometimes called the blazar sequence (Fossati et al. 1998; Kubo et al. 1998). This could also be caused by contamination of non-jet emission components.

The bottom right panel shows that the high-energy component has an anticorrelation between peak frequency and luminosity for radio galaxies that is similar to blazars. The wide range of peak frequencies is also similar to that of blazars. Note that the outlier object with a low peak frequency of 10^19.7 Hz and a low luminosity of 10^41.5 erg s⁻¹ is Cen A. This trend is consistent with the right panel of Figure 2. This feature can be understood as follows. When the GeV gamma-ray spectrum is soft, implying that the GeV band is from the high-energy tail of the inverse Compton (IC) scattering component, the X-ray emission corresponds to the low-energy tail of the IC component, which will be hard, or the X-rays could be from coronal emission, with typical X-ray photon index around 1.7–1.9 in Seyfert galaxies (Nandra & Pounds 1994). When the GeV gamma-ray spectrum is hard, implying that the GeV emission comes from the low-energy tail of the IC component, the X-ray emission corresponds to the high-energy tail of the synchrotron radiation, producing a soft spectrum.

4. GeV Gamma-Ray Luminosity Function

4.1. Trial Luminosity Function

An LF represents a source number density as a function of luminosity and redshift. The luminosity-dependent density evolution (LDDE) model is often used for radio-quiet AGNs and blazars (Ueda et al. 2003, 2014; Narumoto & Totani 2006; Ajello et al. 2015). It is given by

$\begin{eqnarray}\begin{array}{rcl}{\rm{\Phi }}({L}_{\gamma },z) & = & \displaystyle \frac{{d}^{2}N}{{{dzdL}}_{{\rm{X}}}}\\ & = & \displaystyle \frac{A}{\mathrm{ln}(10){L}_{\gamma }}{\left[{\left(\displaystyle \frac{{L}_{\gamma }}{{L}_{* }}\right)}^{{\gamma }_{1}}+{\left(\displaystyle \frac{{L}_{\gamma }}{{L}_{* }}\right)}^{{\gamma }_{2}}\right]}^{-1}\end{array}\end{eqnarray} \tag{ 2 }$

$\begin{eqnarray}&&\ \times \,{\left[{\left(\displaystyle \frac{1+z}{1+{z}_{c}({L}_{\gamma })}\right)}^{{p}_{1}^{* }}+{\left(\displaystyle \frac{1+z}{1+{z}_{c}({L}_{\gamma })}\right)}^{{p}_{2}^{* }}\right]}^{-1},\end{eqnarray} \tag{ 3 }$

where ${z}_{c}{({L}_{\gamma })={z}_{c}({L}_{\gamma }/{L}_{z})}^{\alpha }$ . Here z and L_γ are the redshift and gamma-ray luminosity in a certain energy band, respectively. We assume L_z = 10^42.5 erg s⁻¹. The model parameters are a normalization constant A, a characteristic gamma-ray luminosity L_*, a characteristic redshift z_c, two luminosity indices γ₁ and γ₂, two redshift indices ${p}_{1}^{* }$ and ${p}_{2}^{* }$ , and α. We adopt

$\begin{eqnarray}&&{p}_{1}^{* }\,=\,{p}_{1}+\tau \left({\mathrm{log}}_{10}{L}_{\gamma }-{\mathrm{log}}_{10}{L}_{p}\right)\end{eqnarray} \tag{ 4 }$

and

$\begin{eqnarray}&&{p}_{2}^{* }\,=\,{p}_{2}+\delta \left({\mathrm{log}}_{10}{L}_{\gamma }-{\mathrm{log}}_{10}{L}_{p}\right),\end{eqnarray} \tag{ 5 }$

which were used for the GLF of blazars (Ajello et al. 2015). We set L_p = 10^42.5 erg s⁻¹.

Since our MAGN sample lacks high-redshift objects, we cannot determine the redshift dependence precisely. Therefore, we apply eight previously established redshift evolution models of AGNs from the literature. Models BLLz and FSRQz assume the indices p₁ and p₂ to be the same as those of GLFs of BL Lacs (LDDE₁ in Ajello et al. 2014) and FSRQs (Ajello et al. 2012), respectively, with a free z_c and a single power law with a cutoff at higher luminosity for luminosity dependence. We fix L_* = 10⁴⁶ erg s⁻¹, higher than the luminosity in each band of our sample galaxies, and γ₂ is fixed to 5.0. These two LFs assume that ${p}_{1}^{* }$ and ${p}_{2}^{* }$ do not depend on the luminosity (τ = δ = 0). Models BLLp2 and FSRQp2 assume the redshift dependence of p₁ and z_c to be the same as those of GLFs of BL Lacs (LDDE₁; Ajello et al. 2014) and FSRQs (Ajello et al. 2012), respectively, with a low-redshift index p₂ allowed to be free and the same luminosity dependence as BLLz and FSRQz. Models BLL0, FSRQ0, and BLAZAR0 assume the same LF shape as that of GLFs of BL Lacs (LDDE₂; Ajello et al. 2014), FSRQs (Ajello et al. 2012), and blazars (FSRQ + BL Lac; Ajello et al. 2015), respectively; parameters other than A and L_* are fixed to those from these papers. Thus, the luminosity dependence of the LF is also the same as that of blazars by assumption. The last model RLF0 assumes a similar redshift dependence to that of the radio galaxy LF obtained in the radio band (Willott et al. 2001), which was used in the previous study of GLFs of radio galaxies (Inoue 2011; Di Mauro et al. 2014). It assumes a single power law with an exponential cutoff for luminosity dependence with the index as a free parameter. See Table 3 for details on fixed and free parameters.

Table 3. Best-fit Parameters for the LDDE Model

Model	${\mathrm{log}}_{10}A$ ^a	L_* ^b	γ₁	γ₂	p₁	p₂	z_c	α	τ	δ	$\mathrm{ln}{{ \mathcal L }}_{\max }$ ^c	K-S_zL ^d	K-S_z ^e
`BLLz`	−11.55 ± 0.21	2.00	5.00	1.40 ± 0.05	3.40	−13.00	10^{−8.25±4.11}	0.045	0.00	0.00	−3269.6	55.7	46.3
`FSRQz`	−11.57 ± 0.19	2.00	5.00	1.41 ± 0.05	3.40	−13.00	10^{−6.38±4.17}	0.210	0.00	0.00	−3269.6	63.2	57.3
`BLLp2`	−12.98 ± 0.26	2.00	5.00	1.42 ± 0.06	3.40	6.09 ± 7.07	1.36	0.045	0.00	0.00	−3269.3	69.9	61.9
`FSRQp2`	−13.12 ± 0.25	2.00	5.00	1.62 ± 0.06	3.40	9.09 ± 14.42	0.56	0.210	0.00	0.00	−3270.7	47.5	44.3
`BLL0`	−6.68 ± 0.09	−0.13 ± 0.14	1.86	0.27	7.37	−2.24	1.09	0.045	0.00	−4.92	−3276.0	0.1	0.2
`FSRQ0`	−4.61 ± 0.20	−2.06 ± 0.17	1.58	0.21	7.35	−6.51	0.56	0.210	0.00	0.00	−3270.3	61.1	93.2
`BLAZAR0`	−6.34 ± 0.15	−0.56 ± 0.18	1.83	0.50	4.96	−3.39	0.90	0.072	−0.64	−3.16	−3269.5	72.0	74.7
`RLF0`	−10.71 ± 0.19	2.00	5.00	1.45 ± 0.04	3.40	−3.40	2.00	0.000	0.00	0.00	−3290.2	0.5	0.6

Notes.

^aIn units of Mpc⁻³.^bIn units of 10⁴⁴ erg s⁻¹.^cMaximum likelihood.^dK-S probability in units of percent for luminosity distribution.^eK-S probability in units of percent for redshift distribution.

Download table as: ASCII Typeset image

4.2. Parameter Constraint

Parameters of the trial LF Φ(L_γ, z) are determined by the likelihood analysis. In order to maximize the log-likelihood $\mathrm{ln}{ \mathcal L }$ , an MCMC method using the adaptive Metropolis algorithm (Haario et al. 2001) is applied, following Yamada et al. (2020). The log-likelihood function is defined as

$\begin{eqnarray}\begin{array}{rcl}\mathrm{ln}{ \mathcal L } & = & \displaystyle \sum _{i}^{{N}_{\mathrm{obs}}}\mathrm{ln}\left({\rm{\Phi }}({L}_{\gamma ,i},{z}_{i})\right)\\ & - & {\displaystyle \int }_{{z}_{\min }}^{{z}_{\max }}{dz}{\displaystyle \int }_{{L}_{\gamma ,\min }}^{{L}_{\gamma ,\max }}{{dL}}_{\gamma }{\rm{\Phi }}({L}_{\gamma },z)S({L}_{\gamma },z),\end{array}\end{eqnarray} \tag{ 6 }$

where S(L_γ, z) is the sky coverage function described below. We set ${z}_{\min }=0$ , ${z}_{\max }=6$ , ${L}_{\gamma ,\min }=2\times {10}^{40}\,\mathrm{erg}\ {{\rm{s}}}^{-1}$ , and ${L}_{\gamma ,\max }={10}^{46}\,\mathrm{erg}\ {{\rm{s}}}^{-1}$ . The redshift and luminosity ranges are based on the observed ranges.

Sky coverage S(L_γ, z) as a function of flux is calculated from the public 4FGL-DR2 sensitivity map⁹ at ∣b∣ > 20°. This map is for the 0.1–100 GeV band and was created assuming a single-power-law source energy spectrum with a photon index of 2.2 (Figure 5). At first we scaled it to the function for the 1–3 GeV energy band by using a power-law spectrum with a photon index of 2.2. Then, we scaled the function for the 1–3 GeV band to that for the energy bands of 0.1–0.3 GeV, 0.3–1.0 GeV, 3.0–10.0 GeV, 10.00–30.0 GeV, and 30.00–300.0 GeV by the sensitivity ratio between 1 and 3 GeV and each band, where the Fermi-LAT sensitivity curve as a function of energy¹⁰ is referred to.

**Figure 5.** Sky coverage function at 0.1–100 GeV of the Fermi-LAT 4FGL-DR2 catalog, created from the sensitivity map.
Download figure:
Standard image High-resolution image

4.3. Gamma-Ray Luminosity Function Based on the 1–3 GeV Flux

We constrain the LF using the 4FGL-DR2 1–3 GeV band flux. The Fermi-LAT has the best sensitivity in this band, and we can use the largest number of detected MAGNs. We use MAGNs, which are detected in 1–3 GeV. We exclude MAGNs whose flux is lower than the 1–3 GeV detection limit, where the sky coverage introduced in the previous subsection becomes 0. We do not include CSS and SSRQ galaxies. Since they have a higher redshift and higher luminosity than others and their gamma-ray emission mechanisms are still veiled in mystery, they are a different population from other MAGNs. We also exclude MAGNs located near the Galactic plane at ∣b∣ ≤ 20°. As a result, 34 radio galaxies are used for the GLF, where seven AGNs are included. Table 3 summarizes the parameters obtained for LFs for eight trial models in this energy band.

Overall, the eight models reproduce the data with similar likelihood values, but BLL0 and RLF0 give a significantly lower likelihood than others. BLLz and FSRQz give small z_c = 0.001, and BLLsp2 and FSRQsp2 give a positive p₂, meaning that these four models give a negative evolution. A luminosity-dependent slope is obtained to be γ₂ = 1.3–1.5. BLL0, FSRQ0, and BLAZAR0 give a characteristic gamma-ray luminosity L_* of ∼10⁴⁴, ∼10⁴², and ∼10^43.4 erg s⁻¹, respectively, which are several orders of magnitude smaller than those of blazars (10⁴⁷–10⁴⁸ erg s⁻¹). A smaller L_* for FSRQ0 than the other two could be due to a flatter slope of the luminosity dependence, γ₂ = 0.21, of the LF at higher luminosity.

In order to quantitatively look at how well the obtained LF matches the data, we calculate the observed redshift, luminosity, and cumulated source count (log N–log S) distributions, respectively, as follows:

$\begin{eqnarray}&&\displaystyle \frac{{dN}}{{dz}}=4\pi {\int }_{{L}_{\gamma ,\min }}^{{L}_{\gamma ,\max }}{\rm{\Phi }}({L}_{\gamma },z)S({L}_{\gamma },z){{dL}}_{\gamma }\displaystyle \frac{{d}^{2}V}{{dzd}{\rm{\Omega }}}\end{eqnarray} \tag{ 7 }$

$\begin{eqnarray}&&\displaystyle \frac{{dN}}{{{dL}}_{\gamma }}=4\pi {\int }_{{z}_{\min }}^{{z}_{\max }}{\rm{\Phi }}({L}_{\gamma },z)S({L}_{\gamma },z){dz}\displaystyle \frac{{d}^{2}V}{{dzd}{\rm{\Omega }}}\end{eqnarray} \tag{ 8 }$

$\begin{eqnarray}&&N(\gt {S}_{0})={\int }_{{z}_{\min }}^{{z}_{\max }}{dz}\displaystyle \frac{{d}^{2}V}{{dzd}{\rm{\Omega }}}{\int }_{{L}_{0}}^{{L}_{\gamma ,\max }}{\rm{\Phi }}({L}_{\gamma },z){{dL}}_{\gamma },\end{eqnarray} \tag{ 9 }$

where d² V/dzdΩ is the comoving volume per redshift per solid angle, L₀ = 4π D_L(z)² S₀(1 + z)^Γ−2, D_L(z) is the luminosity distance, S₀ is the observed gamma-ray flux, and Γ is the power-law photon index of the gamma-ray spectra. Here we fix Γ = 2.25, the mean value for our sample (see the next subsection).

The three panels of Figure 6 show model curves together with the observed data points. We correct data points by the sky coverage factor for the log N–log S plot. Hereafter we do not show the model curves of FSRQz and BLLp2, since they give almost identical curves to that of BLLz.

All the models reproduce the data distribution of redshift and luminosity well, but BLL0 and RLF0 models give a somewhat larger deviation from the data for the redshift and luminosity distribution. We calculate the Kolmogorov–Smirnov (K-S) probability (Eadie et al. 2006) for luminosity and redshift distributions for each model, shown in Table 3. The BLL0 and RLF0 models show statistically different behavior compared with the data based on the K-S test values obtained. Therefore, we reject these models.

Figure 6 also shows the model curve based on the LF used in Inoue (2011) for comparison. This LF (RLF) was obtained in the radio band for radio galaxies (Willott et al. 2001). For this RLF, we set ${L}_{\gamma ,\min }={10}^{39}$ erg s⁻¹ and ${L}_{\gamma ,\max }={10}^{48}$ erg s⁻¹, following Inoue (2011). As the K-S test values are small, ∼0.01, this LF does not adequately represent the GLF of radio galaxies. This indicates that GeV gamma-ray-emitting radio galaxies are not the same population as radio-emitting radio galaxies, as discussed in Section 6.2.

For the log N–log S plot, BLLz and FSRQp2 reproduce the data well, while the other models tend to underestimate the number count at higher flux. Since these LFs have a strong positive evolution (as shown later in Figure 8), low-z radio galaxies are predicted not to be dominant in the number count, even at a higher flux. Therefore, we only consider the BLLz and FSRQp2 models hereafter and show the model curves of only these models, together with the model curve of BLAZAR0 as a reference.

Figure 7 shows a visualization of LFs for two redshift regimes: 0.0 ≤ z < 0.1 and 0.1 ≤ z < 1.5. We use the "N^obs/N^mdl" method (La Franca & Cristiani 1997; Miyaji et al. 2001). We deconvolve the observed data points by dividing them by ${N}_{i}^{\mathrm{obs}}/{N}_{i}^{\mathrm{mdl}}$ , where ${N}_{i}^{\mathrm{obs}}$ and ${N}_{i}^{\mathrm{mdl}}$ are the observed and predicted number of radio galaxies in that bin, respectively, and are calculated by integrating Φ(L_γ, z)S(L_γ, z) and Φ(L_γ, z), respectively, in that bin. Note that BLLz and FSRQp2 assume a single-power-law luminosity dependence, while BLAZAR0 assumes a broken power-law luminosity dependence. All models match the data well, and they would be distinguished if radio galaxies at various redshifts were sampled. Figure 8 visualizes the obtained LFs as a function of redshift for various luminosity ranges with the "N^obs/N^mdl" method. Here we plot the LF in a redshift range where the observed radio galaxies exist. Different redshift dependence among LFs is clearly seen. In these figures, we cannot say which model matches the data the best. On the other hand, BLLz and FSRQp2 seem to reproduce the log N–log S relation the best as described above (Figure 6). This indicates that the gamma-ray-emitting radio galaxies have a lower redshift peak in the LF than blazars and other AGNs.

**Figure 8.** GLF for the comoving number density of our sample in various luminosity bins. Data points are deconvolved by dividing them by N^obs/N^mdl. Model curves correspond to the best-fit LDDE models at different redshift bins. Left: `BLLz`. Middle: `FSRQp2`. Right: `BLAZAR0`.
Download figure:
Standard image High-resolution image

4.4. Gamma-Ray Luminosity Function Based on the Six Energy Band Flux Measurements

We also constrain the GLF in each of six energy bands where the flux measurements are listed in 4FGL-DR2. In this case, we use the BLLz model, where free parameters are the same as in the 1.0–3.0 GeV band: A, γ₂, and z_c. The sky coverage function for each energy band is used in this fitting. Table 4 summarizes the fitting results. The selection condition for MAGNs used in this analysis is the same as that for the 1–3 GeV band. The parameters γ₂ and z_c are similar among all energy bands, but A has a dependence on the energy band.

Table 4. Best-fit Parameters for the BLLz Model Obtained in Each of Six Bands

E_band	N_gal ^a	${\mathrm{log}}_{10}A$ ^b	γ₂	z_c	μ	σ	β	$\mathrm{ln}{{ \mathcal L }}_{\max }$ ^c
0.1–0.3	13	−10.78 ± 0.26	1.33 ± 0.07	10^−7.6±3.4				−1284.9
0.3–1	19	−11.13 ± 0.22	1.36 ± 0.07	10^−4.6±2.6				−1855.2
1–3	34	−11.55 ± 0.21	1.40 ± 0.05	10^−10.0±4.1				−3269.6
3–10	33	−11.93 ± 0.23	1.47 ± 0.06	10^−4.9±2.4				−3149.7
10–30	27	−12.00 ± 0.28	1.50 ± 0.07	10^−1.6±4.5				−2571.0
30–300	8	−11.87 ± −0.42	1.48 ± 0.10	10^−1.8±4.8				−772.9

0.1–300	34	−11.50 ± 0.11	1.42 ± 0.03	10^{−9.7±−1.9}	2.26 ± 0.02	0.25 ± 0.01		−17765.1
0.1–300	34	−11.62 ± 0.23	1.52 ± 0.04	10^−0.3±0.3	2.30 ± 0.02	0.24 ± 0.01	0.06 ± 0.01	−17731.7

Notes.

^aThe number of galaxies used in the fitting.^bIn units of Mpc⁻³.^cMaximum likelihood.

Download table as: ASCII Typeset image

Next, we try a common GLF among the six energy bands, by taking into account a distribution of gamma-ray photon indices Γ, G(Γ). The observed photon index has a wide range of values, 1.4–2.5, and could cause a difference in the GLF among the six energy bands as shown above. In this case, the log-likelihood in the jth energy band is defined as

$\begin{eqnarray}&&\begin{array}{l}\mathrm{ln}{{ \mathcal L }}_{j}=\displaystyle \sum _{i}^{{N}_{\mathrm{obs}}}\mathrm{ln}\left({\rm{\Phi }}({L}_{\gamma ,i},{z}_{i})G({\rm{\Gamma }})\right)\\ \quad -\,{\displaystyle \int }_{{z}_{\min }}^{{z}_{\max }}{dz}{\displaystyle \int }_{{L}_{\gamma ,\min }}^{{L}_{\gamma ,\max }}{{dL}}_{\gamma }{\displaystyle \int }_{{{\rm{\Gamma }}}_{\min }}^{{{\rm{\Gamma }}}_{\max }}d{\rm{\Gamma }}{\rm{\Phi }}({L}_{\gamma },z)\\ \quad \times \,G({\rm{\Gamma }}){S}_{j}({L}_{\gamma },z),\end{array}\end{eqnarray} \tag{ 10 }$

where $G{({\rm{\Gamma }})={G}_{0}\exp [-({\rm{\Gamma }}-\mu )}^{2}/2{\sigma }^{2}]$ and S_j(L_γ, z) is the sky coverage function in the jth energy band. μ and σ are the mean and variance of the photon index distribution, respectively, and we set lower and upper bounds to be ${{\rm{\Gamma }}}_{\min }=1.3$ and ${{\rm{\Gamma }}}_{\max }=3.0$ , respectively. Accordingly, ${ \mathcal L }={\prod }_{j}{{ \mathcal L }}_{j}$ is maximized. Here only the BLLz model is applied for the GLF, and the obtained parameters are summarized in Table 4. The parameters obtained for Φ(L_γ, z) are almost the same as those for the 1–3 GeV band fitting. The mean and variance of the photon index are 2.26 ± 0.02 and 0.26 ± 0.04, respectively, and similar to values 2.1 and 0.26, respectively, for BL Lacs (Ajello et al. 2014). Since the L_γ-dependence of the SED is seen in the left panel of Figure 3, the gamma-ray photon index could depend on L_γ. Thus, the model, which has a mean photon index $\mu ({L}_{\gamma })={\mu }^{\star }+\beta \left({\mathrm{log}}_{10}{L}_{\gamma }-43\right)$ , was also tried, and the result is summarized in Table 4. Indeed, the L_γ-dependent coefficient β = 0.07 ± 0.01; since it is positive, it indicates a softer GeV gamma-ray spectrum for galaxies with a higher L_γ. This could not be due to the selection bias of Fermi-LAT detection, since there is no strong bias on the distribution of photon indices for newly detected faint AGNs in the latest 4th catalog (Ajello et al. 2020).

4.5. Based on Different Classes of Galaxies

In the above analysis, we include RDG and AGN classes in the 4FGL-DR2. We checked how the result changes if we use only RDG, or include CSS and/or SSRQ, using the BLLz model. The result is that the GLF parameters are not strongly affected by whether or not we include these galaxies; the normalization changed with the number of sample galaxies. When we include CSS and SSRQ, the luminosity dependence and redshift dependence became a little bit flatter because CSS and SSRQ consist of galaxies with high redshifts and luminosities. If we exclude AGN, these become steeper because AGN tends to have higher redshifts and luminosities. But the overall GLF shape does not change much.

5. Contribution to the Extragalactic Gamma-Ray Background

We calculate the contribution of gamma-ray-emitting radio galaxies to the extragalactic gamma-ray background (EGB) flux using the obtained GLF as

$\begin{eqnarray}\begin{array}{rcl}{F}_{\mathrm{EGB}}({E}_{0}) & = & {\displaystyle \int }_{{z}_{\min }}^{{z}_{\max }}\displaystyle \frac{{d}^{2}V}{{dzd}{\rm{\Omega }}}{dz}{\displaystyle \int }_{{L}_{\gamma ,\min }}^{{L}_{\gamma ,\max }}{{dL}}_{\gamma }\\ & & \times \,{\rm{\Phi }}({L}_{\gamma },z)F(z,{L}_{\gamma },{E}_{0}),\end{array}\end{eqnarray} \tag{ 11 }$

where F(z, L_γ, E₀) is a flux at an energy E₀ in the observer frame for a source with a redshift z and a gamma-ray luminosity L_γ. We apply two different methods for the calculation of F(z, L_γ, E₀). The first method treats F(z, L_γ, E₀) as a summation of a single power law with various photon indices Γ following a Gaussian distribution G(Γ), derived in Section 4.4. We add a high-energy spectral cutoff at 500 GeV, since PKS 0625–354, one of the hardest Γ radio galaxies, has a cutoff around this energy (HESS Collaboration et al. 2018b). We use the GLF obtained in the 1–3 GeV band and extrapolate the contribution of this band to other energy bands. In the second method we treat F(z, L_γ, E₀) as a single power law with an average photon index in the jth energy band and calculate the contribution in the jth energy band using the obtained GLF in each energy band. In both cases, absorption by the extragalactic background light is taken into account by using the formula of Inoue et al. (2013).

The left panel of Figure 9 shows the EGB spectrum using the former method. The contribution of radio galaxies to the EGB is predicted to ∼ a few percent from the BLLz LF model. This is smaller than that reported by previous studies, where radio LFs were used (Inoue 2011; Di Mauro et al. 2014). At higher energy, the predicted spectrum could exceed the observed EGB if there is no spectral cutoff. The FSRQp2 LF model gives a smaller contribution than BLLz because it predicts a smaller number of radio galaxies at high redshifts (see Figure 10). The BLAZAR0 LF model predicts a higher contribution, around 5%–10% of the EGB, but still smaller than previous estimates. This is caused by a higher redshift peak in these LFs than BLLz and FSRQp2; this LF predicts a larger number of radio galaxies at higher redshift.

**Figure 10.** Left: comparison of best-fit GLFs (as comoving number density). Colors of the three best-fit models are the same as in Figure 6. Curves in three luminosity ranges of 10^40.5–10⁴² erg s⁻¹, 10⁴²–10⁴⁴ erg s⁻¹, and 10⁴⁴–10⁴⁶ erg s⁻¹ are plotted for each model. As a comparison, the GLF of all blazars (Ajello et al. 2015) is also shown as blue thick dashed lines. Four curves for all blazars correspond to luminosity ranges of 10^43.8−10^46.8−10^44.8 erg s⁻¹, 10^45.8−10^46.9 erg s⁻¹, 10^46.9−10^47.9 erg s⁻¹, and 10^47.9−10^49.4 erg s⁻¹. Right: number distribution of radio galaxies against redshift, predicted by best-fit GLF models. Model curves are calculated without sky coverage function. Colors are the same as those in the left panel.
Download figure:
Standard image High-resolution image

The right panel of Figure 9 shows the contribution based on the GLF obtained in each of the six energy bands by using the BLLz model. It is almost consistent with the prediction by using a summation of a cutoff power law with a Gaussian distribution of photon indices. It predicts a smaller contribution in the lower energy band. This may be due to the result of the assumption of the Gaussian distribution for Γ in the left panel.

Next, by using the obtained SED shape of our sample radio galaxies from X-ray to GeV gamma-ray band in Section 3, we calculate a contribution to the EGB in the MeV gamma-ray band. Our sample radio galaxies are divided into the four luminosity groups; 10⁴⁰–10⁴² erg s⁻¹, 10⁴²–10⁴³ erg s⁻¹, 10⁴³–10⁴⁴ erg s⁻¹, and 10⁴⁴–10⁴⁶ erg s⁻¹ in the 1–3 GeV band, and an average SED shape is obtained in each luminosity range. These average SEDs are shown as light-green lines in the left panel of Figure 3. Extrapolating these four SED shapes to arbitrary luminosity, we calculated a contribution to the EGB in the MeV energy band. The result by using the BLLz LF model is shown in the right panel of Figure 9. The contribution of gamma-ray-emitting radio galaxies in the MeV band is ≈ a few percent.

6. Discussion

6.1. Comparison with Blazars

In this study, we adopt eight LDDE GLF models to fit the data. Among them, the BLLz and FSRQp2 models reproduced the observed redshift, luminosity, and source count distributions best. These two models have the same redshift dependence as those of GeV gamma-ray BL Lacs and FSRQs, respectively, except for the peak redshift z_c for the former and p₂ for the latter. For BLLz, z_c < 0.1, which is lower than the typical peak redshift of blazars and radio-quiet AGNs, where z_c ≈ 2. For FSRQp2, p₂ > 0. Therefore, the GeV GLF of radio galaxies is similar to that of blazars, but a lower peak redshift or a negative evolution is preferred.

The left panel of Figure 10 shows a comparison of the best-fit GLFs (BLLz and FSRQp2), together with those of blazars. Although these two models give similar number densities in the redshift ranges of our sample galaxies, their behavior at higher redshifts is different.

Radio galaxies are believed to be the parent population of blazars, so it is natural to expect that their LF shape is similar. Our sample is dominated by low-power radio galaxies with similar radio power to that of FR I galaxies, the parent population of BL Lacs. However, as seen in the left panel of Figure 10, a clear difference in LF shape exists between gamma-ray-emitting radio galaxies and blazars. Compared to the blazar GLF, both BLLz and FSRQp2 models give smaller number counts at higher redshifts and larger number counts at lower redshifts. We also show the model curve of BLAZAR0, which follows the GLF of the combined population (BL Lac and FSRQ) of gamma-ray blazars (Ajello et al. 2015). Even with this model, the GLF curves of radio galaxies appear to be different from those of blazars.

Low-luminosity high-synchrotron-peaked BL Lacs (HSPs) show a strong negative evolution (Ajello et al. 2014) like gamma-ray-emitting radio galaxies. If HSPs are the parent population of gamma-ray-emitting radio galaxies, this may be a possible reason for their GLF shapes. However, HSPs are a specific population of blazars. In fact, low-z gamma-ray-emitting radio galaxies do not necessarily have a high peak frequency or a hard GeV gamma-ray photon index as seen in the right panel of Figure 11.

**Figure 11.** Left: gamma-ray luminosity vs. gamma-ray photon index for 4FGL-DR2 MAGNs, together with blazars. Colors and symbols are the same as those of Figure 4, except for the yellow crosses, which are the BCUs. Right: redshift vs. gamma-ray photon index for 4FGL-DR2 MAGNs, together with blazars. Colors and symbols are the same as those of the left panel.
Download figure:
Standard image High-resolution image

Another possible interpretation is that emission regions of gamma-ray radio galaxies seen by Fermi may be different from those of blazars. Then, an additional component appears in low-redshift GLFs. Recent gamma-ray observations found evidence that a large population of low-z, low-luminosity radio galaxies might be due to gamma-ray emission from sources other than the core jet. Cen A and Fornax A show extended gamma-ray emission from their radio lobes (Abdo et al. 2010e; Ackermann et al. 2016). Cen A also shows the existence of kiloparsec-scale extended gamma-ray emission (Prokhorov & Colafrancesco 2019; HESS Collaboration et al. 2020), where various models are proposed to explain the origin (e.g., Tanada et al. 2019; Sudoh et al. 2020; Rieger & Duffy 2021). Hot spot emission could also contribute to gamma-ray emission in radio galaxies, as suggested for M87 HST-1 (Harris et al. 2009; Imazawa et al. 2021). This non-core emission could make up a significant fraction of GeV-detected gamma-rays from radio galaxies. Therefore, the difference in GLF shapes between radio galaxies and blazars could be due to an additional component from non-core emission, appearing in low redshifts.

To investigate the relation between blazars and gamma-ray-emitting radio galaxies further, in the right panel of Figure 10 we show the redshift distribution of radio galaxies, obtained by using Equation (9) but without sky coverage. Thus, nondetected galaxies are considered. We also show the redshift distribution of blazars for comparison. The number of radio galaxies is 10–30 times larger than that of BL Lacs.

From the viewpoint of population, the comparison with the curve of BL Lacs is natural since low-power radio galaxies dominate our samples. However, the ratio of 10–30 is unexpectedly small when we consider the beaming effect. Our sample radio galaxies have a lower gamma-ray luminosity by a factor of ∼10⁴–10⁵ than blazars. Considering the beaming correction factor δ⁴, where δ is a beaming factor, radio galaxies in 4FGL-DR2 should have a smaller δ by a factor of ∼10 than blazars. In that case, the ratio should be ∼100. Inoue (2011) proposed that only radio galaxies with a viewing angle of ≲ 24° are seen in the GeV gamma-ray band.

Alternatively, if Fermi-LAT misses faint radio galaxies at the lowest-luminosity end, the number ratio of GeV-emitting radio galaxies to BL Lacs could be underestimated. Even for nearby galaxies, the detection limit in luminosity is around 10⁴¹ erg s⁻¹, and thus radio galaxies with a luminosity of 10⁴⁰–10⁴¹ erg s⁻¹ could be a hidden population.

6.2. Comparison with Radio Galaxies Seen in the Radio Band

The RLF model, which is obtained in the radio band, cannot reproduce the observed distributions of gamma-ray-emitting radio galaxies at all (see Figure 6). This indicates that a population of GeV-emitting radio galaxies is different from or represents a peculiar population of radio galaxies seen in the radio band. This section investigates the differences between gamma-ray detected and nondetected radio galaxies.

Mingo et al. (2014) compiled X-ray properties of the 2 Jy complete sample of radio galaxies, where only 4 out of 46 galaxies are GeV gamma-ray-emitting ones. Some of our sample radio galaxies have fainter X-ray fluxes compared to the Mingo et al. (2014) sample. Thus, the smaller fraction of our sample compared with the Mingo et al. (2014) sample is due not to the X-ray flux limit but to only a small fraction of radio galaxies being gamma-ray loud. Gamma-ray-quiet radio galaxies do not follow the correlation between X-ray and gamma-ray luminosity shown in Figure 2. This ratio (4/46) is similar to the population ratio estimate in Inoue (2011), which estimated the ratio as 0.088 based on the comparison of the number counts between the radio and GeV gamma-ray bands. They claimed that a large population of radio galaxies are faint in the gamma-ray band owing to a large viewing angle, which makes the jet faint. As seen in Figure 1, the gamma-ray luminosity is different by more than two orders of magnitude among galaxies with the same radio luminosity, when considering those not detected by Fermi-LAT. Such a large difference of gamma-ray luminosity could be caused by the beaming effect; gamma-ray emission comes from the jet core, and thus the gamma-ray luminosity is strongly affected by the beaming effect, while this radio emission comes not only from the core but also from an extended outer region, and thus it is less beamed.

GeV-emitting radio galaxies of our sample comprise 22 FR I galaxies and 14 FR II galaxies. Unclassified galaxies in our sample do not tend to have a higher radio power like FR II. Therefore, observationally low radio power galaxies like FR I are the dominant class in the GeV gamma-ray band. This is contrary to the radio-band population; there are 6 FR I galaxies and 33 FR II galaxies in the 2 Jy complete sample of radio galaxies (Dicken et al. 2008). Therefore, it is clear that the population of gamma-ray-emitting radio galaxies lacks high radio power galaxies like FR II galaxies.

In contrast, there are eight CSS and SSRQ galaxies in the GeV gamma-ray band, while there are six CSS galaxies in the 2 Jy radio sample. Therefore, the fraction of CSS and SSRQ is similar in the radio and GeV gamma-ray bands. CSS and SSRQ are high-luminosity radio galaxies and thus an intermediate class between FR II and FSRQ.

As seen in Figure 11, all CSS and many FR II galaxies show steep gamma-ray spectra. Many galaxies having steep spectra are FR II. Therefore, it is likely that the gamma-ray SED of FR II galaxies tends to drop down below the Fermi-LAT energy band owing to the low peak frequency of their high-energy component. This is also seen in Figure 4, where FR II galaxies tend to have lower peak frequencies.

It has been suggested that many FSRQs are not detected by Fermi-LAT owing to their low peak frequency of the IC component (Paliya et al. 2017). If FR II galaxies are a parent population of FSRQs and the IC peak appears lower owing to less beaming, they are more likely to be missed in the GeV band. Note that 80% of radio galaxies are FR II in an X-ray flux-limited Swift/BAT sample (Rusinek et al. 2020). This is similar to the case of a radio flux-limited sample (Dicken et al. 2008); X-ray surveys do not have a strong bias for radio galaxies. In the X-ray band, FR II galaxies have high-luminosity disk/corona emission like 3C 111 and 3C 120 (Fukazawa et al. 2015). Such objects have jet emission that extends up to at least the GeV gamma-ray band. Therefore, if many FR II galaxies will be detected in the MeV gamma-ray band by future missions such as COSI (Tomsick & COSI Collaboration 2022), AMEGO-X (Fleischhack & Amego X Team 2022), and GRAMS (Aramaki et al. 2020), the picture of the population will be more established from the viewpoint of the SED shape.

The other effect to make FR II galaxies fainter in the GeV gamma-ray band might be the jet structure. Another effect that might make FR II galaxies fainter in the GeV gamma-ray band is the jet structure. FR I galaxies tend to have a structured jet, and thus emission from a sheath region with a lower Lorentz factor is less beamed, and thus GeV gamma-rays can be observed even if the jet is misaligned (Ghisellini et al. 2005). FR II emission is more beamed, and thus the misaligned situation makes it less luminous. Keenan et al. (2021) discussed a similar model to explain the synchrotron peak frequency–luminosity relation of blazars and radio galaxies.

Figure 4 shows that the peak luminosity of the X-ray to gamma-ray SED is not significantly higher than that of the synchrotron peak luminosity for any radio galaxy, unlike in FSRQs. The external Compton component has a different beaming pattern than synchrotron or synchrotron self-Compton (Dermer 1995). The Compton dominance is proportional to $\max \left({\delta }^{2}{u}_{\mathrm{ext}},{u}_{\mathrm{sy}}\right)/{u}_{{\rm{B}}}$ , where δ is the Doppler factor and u_ext, u_sy, and u_B are the energy densities of the external radiation field, synchrotron photons, and magnetic field, respectively (Finke 2013). For radio galaxies with high viewing angles, δ will be small compared to blazars. Therefore, it might be that the external Compton component cannot be observed for FR II galaxies because it is below the synchrotron self-Compton component.

Keenan et al. (2021) proposed that FSRQs, LBLs, and high-excitation radio galaxies (HERGs; radio galaxies with high-excitation optical emission lines) have a strong jet with a wide range of jet power, while intermediate- or high-frequency-peaked BL Lacs (IBL or HBL) and low-excitation radio galaxies (LERGs; radio galaxies with low-excitation optical emission lines) have weak jets with low powers. The FSRQs, LBLs, and HERGs could correspond to FR II galaxies, and the IBLs, HBLs, and LERGs could correspond to FR I galaxies. However, some of the gamma-ray-emitting radio galaxies (PKS 0625–354, 3C 264, and TXS 1516+064) have high peak synchrotron frequencies of 10^15.5–10^17.5 Hz, and their high-energy components have peak frequencies of 10^24.5–10^25.0 Hz. Others have similar high peak frequencies of the high-energy component. These galaxies have a low gamma-ray luminosity compared to BL Lacs, and thus they do not follow the model of Keenan et al. (2021).

6.3. The `AGN` Class in 4FGL-DR2

The 4FGL-DR2 sample contains the AGN class, whose properties and counterparts are not well studied in other wavelengths. These galaxies cannot be classified as blazars with uncertain type (BCU class in 4FGL-DR2), since they do not satisfy the criteria of blazars. Their SEDs are similar to those of compact radio sources, but the uncertainty is large. As shown in Figure 11, AGNs seem to be divided into low-luminosity ones and high-luminosity ones, corresponding to FR I and FR II, respectively.

One candidate population of this class is FR 0 galaxies (Ghisellini 2011), where the radio emission is compact. In fact, NVSS radio images of most of these galaxies on the SkyView ¹¹ do not exhibit a clear jet structure. Torresi et al. (2018) systematically studied X-ray properties of FR 0 galaxies and found that their X-ray properties are almost the same as those of FR I galaxies. Therefore, some low-luminosity AGNs could be FR 0 galaxies. Paliya (2021) reported a possible gamma-ray-emitting population of FR 0 galaxies in a stacking analysis of Fermi-LAT data. Itoh et al. (2020) reported that there are a significant number of elliptical-like blazars in their blazar catalog. These objects are low-luminosity BL Lacs and apparently look like elliptical galaxies. Therefore, such objects with a misaligned jet are one candidate population of the AGN class.

D'Ammando et al. (2015) reported that PKS 0521–36, the brightest among the AGN class objects, has properties intermediate between broad-line radio galaxies (BLRGs) and SSRQs. BLRGs have a high accretion rate and thus higher luminosity. Therefore, high-luminosity AGNs could be a transition class between FR II galaxies and blazars.

6.4. X-Ray to Gamma-Ray Spectra

The study of X-ray emission from radio galaxies is important for constraining the broadband emission from jets because so far jet emission from radio galaxies has been primarily detected in the radio and gamma-ray bands. Fukazawa et al. (2015) systematically analyzed Suzaku X-ray spectra of 1FGL radio galaxies. They reported that the X-ray emission is dominated by disk/corona emission for HERGs and jet emission for LERGs.

Since HERGs have higher accretion rates, their luminosities should be higher. HERGs tend to have harder X-ray spectra, while LERGs tend to have softer spectra. This is consistent with the result of high-energy component SEDs from the X-ray to GeV gamma-ray band (Figure 4). Higher-luminosity radio galaxies show lower peak frequencies for the high-energy component, where the X-ray emission is the low-energy tail of the inverse Compton scattering component from the jet and/or disk/corona emission. Bright disk/corona X-ray emission can make peak frequencies of a high-energy component lower and thus, like the blazar sequence, could be enhanced for high-energy components of radio galaxies.

LERGs show a higher peak frequency for their high-energy components, where the X-ray emission is the high-energy tail of jet synchrotron emission. A detailed systematic X-ray study with a classification of radio galaxies will be presented by a separate paper.

The gamma-ray spectral index has a wide distribution, from 1.5 to 3.0. This range is similar to that of blazars, and it is caused by a variety of high-energy component peak frequencies. This seems to be consistent with the view that radio galaxies are the parent population of blazars. However, as we discussed in Section 6.1, the GLF shapes are different. This indicates that some of gamma-ray-emitting radio galaxies are not just misaligned blazars.

For galaxies with hard GeV gamma-ray spectra, cutoffs are not seen in the GeV band. The gamma-ray-brightest galaxy with a hard spectrum is PKS 0625–354, and its cutoff is detected around 0.5 TeV by HESS (HESS Collaboration et al. 2018b). If the cutoff energy was much higher, the contribution to the EGB could exceed the observed flux, as shown in Figure 9. This spectral cutoff is also similar to that of blazars and could be due to the Klein–Nishina effect.

The gamma-ray spectrum of Cen A has an upturn between the GeV and TeV (HESS Collaboration et al. 2018a), suggesting two emission components. In fact, TeV gamma-ray emission from the Cen A kiloparsec-scale lobe was detected by HESS (HESS Collaboration et al. 2020), suggesting particle acceleration in the kiloparsec-scale diffuse or knot regions. Rulten et al. (2020) surveyed the Fermi-LAT spectra of radio galaxies, but only Cen A shows an upturned spectrum. Considering that the peak frequency of the high-energy component of Cen A is the lowest among radio galaxies and thus is the faintest in gamma-rays, future surveys may find other gamma-ray-faint radio galaxies with upturned spectra.

6.5. Contribution to EGB

Various studies estimated the contribution of radio galaxies to the EGB (Inoue 2011; Di Mauro et al. 2014; Hooper et al. 2016; Stecker et al. 2019; Blanco & Linden 2021). Inoue (2011) estimated this contribution by using radio galaxies in the 1FGL, 3FGL, FL8Y (a precursor to 4FGL), and 4FGL. They estimated the contribution in a similar way to what Inoue (2011) did, by assuming that the GLF shape is the same as that determined by radio observation of radio galaxies (that is, RLF; Willott et al. 2001), and that the gamma-ray luminosity correlates with radio luminosity. We determined the GLF normalization by using the number of gamma-ray-emitting radio galaxies, and the contribution to the EGB was estimated to be several to several tens of percent, dependent on the gamma-ray spectral index.

Our estimation is based on only gamma-ray results and has as few assumptions as possible. Although some GLF parameters are fixed, we fitted the distributions of gamma-ray luminosities and redshifts by eight model GLFs and confirmed that they reproduce the data. In addition, we considered an energy dependence of the EGB contribution both by using a Gaussian distribution for the gamma-ray photon index and by estimating the contribution in each of six gamma-ray energy bands independently. As a result, the contribution is estimated to be 1%–10% of the 0.1–300 GeV background, dependent on the GLF shape, especially for the redshift dependence of the GLF. In other words, the radio galaxy contribution to the EGB is dominated by unresolved (nondetected) radio galaxies. In order to constrain the dark matter contribution to the EGB, it is important that the contribution of radio galaxies is estimated as accurately as possible.

Furthermore, by obtaining the X-ray to GeV gamma-ray SED, we estimated the contribution of GeV gamma-ray-emitting radio galaxies to the EGB in the MeV band to be <1%. However, as discussed in Section 6.2, there could be a large number of radio galaxies hidden in the GeV gamma-ray band. Considering that only 10% of the 2 Jy flux-limited sample is detected in the GeV band, but most of them are detected in the X-ray band, the contribution of radio galaxies to the EGB in the MeV band could be larger by a factor of up to 10.

Below 10 MeV, the EGB spectrum becomes softer with a turnover around 5 MeV, while the SEDs of GeV gamma-ray-emitting radio galaxies are assumed to connect directly to the X-ray band, and thus the contribution below 1 MeV becomes much smaller than 1%. Radio-quiet AGNs dominate below several hundred keV (Ueda et al. 2003, 2014; Gilli et al. 2007) because there are many more of them than radio galaxies and their spectra show cutoffs around several hundred keV. These cutoffs are due to thermal emission from the disk/corona, while nonthermal emission from the disk/corona is suggested and could contribute to the MeV EGB (Inoue et al. 2008, 2019; Murase et al. 2020).

Considering that radio galaxies make up 7%–10% of the objects in the Swift/BAT AGN catalog (Panessa et al. 2016; Gupta et al. 2018), the contribution of radio galaxies to the cosmic X-ray background (CXB) could be about 10% of the contribution of Seyfert galaxies (Gilli et al. 2007). This is 100 times larger than the contribution predicted from GeV-emitting radio galaxies. As discussed in Section 6.2, numerous GeV-quiet radio galaxies would have SEDs similar to that of Cen A, as shown in Figure 9. If that is the case, the contribution to the CXB could be 100 times larger than our prediction curve in the right panel of Figure 9. Note that, in such a case, the contribution of radio galaxies to the EGB becomes comparable to that of FSRQs around 1 MeV. A recent study of contributions of FSRQs to the MeV EGB suggested that FSRQ cannot explain the MeV EGB completely (Toda et al. 2020). Therefore, future MeV gamma-ray observations with such missions as COSI, AMEGO-X, and GRAMS could detect many FSRQs and radio galaxies to resolve this issue.

7. Conclusions

Utilizing the 4FGL-DR2 catalog, we systematically studied the X-ray to GeV gamma-ray spectra of 61 gamma-ray-emitting radio galaxies. We found an anticorrelation between peak frequency and peak luminosity for the high-energy component of sampled radio galaxies. This anticorrelation feature is also known to exist in blazars (see, e.g., Finke 2013; Ghisellini et al. 2017), although its origin is controversial. We note that if the disk/corona emission, rather than jet emission, dominates the X-ray fluxes of high-luminosity radio galaxies (Fukazawa et al. 2015), it would enhance this feature.

Comparing with radio- and X-ray-selected radio galaxy samples, the gamma-ray-selected sample of radio galaxies lacks high radio luminosity galaxies like FR II galaxies at the same radio flux threshold. This implies that FR II galaxies appear fainter in the gamma-ray band owing to the beaming effect and/or softer SED shape. Future MeV gamma-ray observations will be crucial to elucidating this discrepancy.

We further explored the cosmological evolution of gamma-ray-emitting radio galaxies using the same sample. For the first time, we construct the GLF of radio galaxies using gamma-ray data only. We found that gamma-ray-emitting radio galaxies favor negative evolution at all luminosity ranges, similar to HSPs. However, this trend is different from all blazars, which are on-axis radio galaxies. Therefore, gamma-ray photons may originate in different regions in gamma-ray-loud radio galaxies and blazars.

By combining the gamma-ray spectra and GLFs of radio galaxies, we also estimated their contribution to the extragalactic gamma-ray background radiation. The expected contribution is about 1%–10%. However, considering hidden GeV gamma-ray-emitting radio galaxies, the contribution to EGB could be around 10% in the MeV band.

Y.I. is supported by JSPS KAKENHI grant Nos. JP18H05458 and JP19K14772. J.F. is supported by NASA under contract S-15633Y.

High-energy Emission Component, Population, and Contribution to the Extragalactic Gamma-Ray Background of Gamma-Ray-emitting Radio Galaxies

Article metrics

Author e-mails

Author affiliations

ORCID iDs

Dates

Abstract

1. Introduction