A Wide and Deep Exploration of Radio Galaxies with Subaru HSC (WERGS). IV. Rapidly Growing (Super)Massive Black Holes in Extremely Radio-loud Galaxies

2.1.1. VLA/FIRST

2.1.2. Subaru/HSC-SSP Data

The HSC-SSP is an ongoing wide and deep imaging survey covering five broadband filters (g, r, i, z, and y band; Aihara et al. 2018a; Bosch et al. 2018; Furusawa et al. 2018; Kawanomoto et al. 2018; Komiyama et al. 2018; Huang et al. 2018), consisting of three layers (wide, deep, and ultra deep). Yamashita et al. (2018) utilized the wide-layer S16A data release, which contains observations from 2014 March to 2016 January with a field coverage of 154 deg² (based on six fields; XMM-LSS, GAMA09H, WIDE12H, GAMA15H, HECTOMAP, and VVDS), and forced photometry of 5σ limiting magnitude down to 26.8, 26.4, 26.4, 25.5, and 24.7 for the g, r, i, z, and y bands, respectively (Aihara et al. 2018b). The average seeing in the i band is 06, and the astrometric root-mean-squared uncertainty is about 40 mas. After removing spurious sources flagged by the pipeline, Yamashita et al. (2018) cross matched the FIRST sources with a search radius of 1''. Yamashita et al. (2018) also required detections with S/N > 5 in the r, i, and z bands to qualify as an optical counterpart. This initial sample of Yamashita et al. (2018) comprised 3579 sources.

We then applied additional cuts to this sample. To remove radio galaxies from the initial WERGS sample whose radio emission might be dominated by the star formation of host galaxies and not by the AGNs, we set a lower-limit to the radio luminosity of L_{1.4 GHz} > 10²⁴ WHz⁻¹, which is equivalent to the radio emission from the host galaxies with a star formation rate of SFR ≈ 250 M_⊙ yr⁻¹ (Condon 1992; Condon et al. 2013). This cut is also supported by a steep decline in the radio luminosity function of starburst galaxies above L_{1.4 GHz} =10²³ WHz⁻¹ (e.g., Kimball et al. 2011; Padovani 2016; Tadhunter 2016). This criterion reduces the sample size to 3147.

Next, we limit our sample to compact sources in the radio band. This is important in order to reduce false optical identification by mismatching the optical sources with the locations of spatially extended radio-lobe emission. Yamashita et al. (2018) discussed this and categorized radio compact sources using the ratio of the total integrated radio flux density to the peak radio flux density f_int/f_peak. They treated the source as compact if the ratio fulfilled either of the two following equations,

$$\begin{eqnarray}&&{f}_{\mathrm{int}}/{f}_{\mathrm{peak}}\lt 1+6.5\times {({f}_{\mathrm{peak}}/\mathrm{rms})}^{-1}\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\mathrm{log}({f}_{\mathrm{int}}/{f}_{\mathrm{peak}})\lt 0.1,\end{eqnarray} \tag{ 2 }$$

where f_peak/rms is the S/N and rms is a local rms noise in the FIRST catalog, and the above two equations are obtained from the study by Schinnerer et al. (2007; see also the original criterion in Ivezić et al. 2002). We set an additional conservative requirement that there shall be no FIRST source in the surrounding 1' (corresponding to ≲500 kpc at z ∼ 1) to select the isolated radio core emission and avoid coincident matching with radio lobes of other nearby FIRST sources. This reduces the sample to 2583 sources.

A further 10 sources were removed in crowded regions of the HSC footprint due to multiple HSC detections within < 1'' of the HSC counterpart, bringing the sample down to 2573 sources.

Finally, we restricted the sample to sources with spectroscopic redshift (spec-z; z_spec) or reliable photometric redshift (photo-z; z_phot) values. The spec-z were obtained from SDSS DR12 (Alam et al. 2015), the GAMA project DR2 (Driver et al. 2011; Liske et al. 2015), and WiggleZ Dark Energy Survey project DR1 (Drinkwater et al. 2010). In this study, we limit the sample to the spec-z range 0.3 < z_spec <1.6 to cover the same redshift range as the photo-z sample. For the sources without spec-z, we utilized the photo-z estimated using the Mizuki SED-fitting code which is one of the standard photo-z packages for the HSC-SSP survey. The method utilizes the photometries of the five HSC-SSP bands (see Tanaka 2015 and Tanaka et al. 2018, for more details) for the SED fitting. Yamashita et al. (2018) discussed that the HSC-SSP photo-z derived by MIZUKI is reliable for z_phot < 1.6 based on comparison with spectroscopic redshift in the COSMOS field. In addition, the redshifts of sources with z_phot < 0.3 are sometimes erroneous because they lack the Balmer-break tracer in the HSC optical bands. Therefore, we limit the sample to 0.3 < z_phot < 1.6. We further imposed requirements based on the reliable photo-z fitting quality (reduced χ² < 3, σ/z_phot < 0.2, and σ/(1 + z_phot) < 0.1), following Toba et al. (2019). The resulting WERGS sample contains 1770 sources. Given that our sample largely relies on photo-z results, we further discussed how possible erroneous photo-z might affect our main results in Appendix A.

To study the AGN and host-galaxy properties in WERGS, Toba et al. (2019) compiled optical, near-IR, mid-IR, far-IR, and radio data for the WERGS sample. Toba et al. (2019) performed spectral energy distribution (SED) fitting with CIGALE (Boquien et al. 2019) and inferred physical properties, including the IR luminosity contributed from AGNs (IR AGN luminosity; L_AGN,IR), dust-extinction-corrected stellar mass (M_⋆), and SFR estimated from the decomposed host-galaxy IR emission. ¹⁴ For the estimation of L_AGN,IR, Toba et al. (2019) considered the contamination of synchrotron radiation in the IR bands by extrapolating the power law from the radio bands. Thus, L_AGN,IR is obtained purely from the dust emission heated by AGNs. The bolometric AGN luminosity is estimated by using the conversion L_AGN,bol ≃ 3 × L_AGN,IR (Delvecchio et al. 2014; Inayoshi et al. 2018). Because Toba et al. (2019) restricted themselves to objects in an area of ∼95 deg², in which multi-wavelength information from u band to far-IR was available, our sample is reduced to 1055 sources.

Some of the selected sources have intensive star formation rates, reaching SFR ≈ 250M_⊙ yr⁻¹, so we applied an additional cut to remove possible star-formation-dominated radio galaxies using the ratio of IR and radio luminosity (q_IR; Helou et al. 1985; Ivison et al. 2010) defined by

$$\begin{eqnarray}&&{q}_{\mathrm{IR}}=\mathrm{log}\left(\displaystyle \frac{{L}_{\mathrm{IR}}/3.75\times {10}^{12}}{{L}_{1.4\,\mathrm{GHz}}}\right),\end{eqnarray} \tag{ 3 }$$

where L_IR is the total IR luminosity in units of W derived from CIGALE and 3.7 × 10¹² is the frequency in Hz corresponding to 80 μm, which makes q_IR a dimensionless quantity. We set q_IR < 1.68 (Del Moro et al. 2013) to select radio-excess (meaning jet-emission-dominated) sources. As a result, our final sample reduces to 1026 sources, and the resulting final number of ERGs fulfilling $\mathrm{log}{{ \mathcal R }}_{\mathrm{rest}}\gt 4$ is 39 sources; they are shown as orange stars in Figure 1. We refer to the remaining 987 WERGS sources with $1\lt \mathrm{log}{{ \mathcal R }}_{\mathrm{rest}}\lt 4$ as "normal radio galaxies (NRGs)" (shown with blue circles in Figure 1) hereafter.

Table 1 shows the number of sources for the full WERGS sample and the ERGs after each selection cut. ERGs suffer from the photo-z quality cut more than the full WERGS sample because ERGs contain a relatively large fraction of the optically faintest sources with i_AB > 25, which sometimes makes the reliable photo-z estimation difficult. The IR selection cut is another main factor in reducing the number of sources. This is mainly because the region available for the IR data is limited in the area of ∼95 deg² out of the 154 deg² (Toba et al. 2019). Note that this IR selection reduces the similar fraction of sources for the full WERGS sample (60%) and ERGs (55%).

3.1.1.  $\mathrm{log}{{ \mathcal R }}_{\mathrm{rest}}$ versus g-band Magnitude

3.1.2.  L–z Plane

Figure 2 shows the radio and IR AGN luminosities of our sample as a function of redshift. The 1.4 GHz radio luminosity (L_{1.4 GHz} W Hz⁻¹) is taken from Yamashita et al. (2018) or Toba et al. (2019), in which the integrated flux density (f_int) obtained by the VLA/FIRST final catalog (Helfand et al. 2015) is used, and k-correction is applied by assuming a power-law radio spectrum with f_ν ∝ ν^α , where α is estimated with $\alpha \,=\mathrm{log}({f}_{1.4\mathrm{GHz}}/{f}_{150\mathrm{MHz}})/\mathrm{log}({\nu }_{1.4\mathrm{GHz}}/{\nu }_{150\mathrm{MHz}})$ for the objects having TGSS 150 MHz data (Intema et al. 2017) and α = −0.7 for all others (e.g., Condon 1992).

Zoom In Zoom Out Reset image size

As shown in the left panel of Figure 2, ERGs show relatively high radio luminosities with a median of $\langle \mathrm{log}({L}_{1.4\mathrm{GHz}}/{\rm{W}}\,{\mathrm{Hz}}^{-1})\rangle =26.3$, which is one order of magnitude higher than that of the NRGs ($\langle \mathrm{log}({L}_{1.4\mathrm{GHz}}/{\rm{W}}\,{\mathrm{Hz}}^{-1})\rangle =25.1$). However, in terms of the IR AGN luminosity shown in the right panel, there is no clear difference between ERGs and the NRGs $\langle \mathrm{log}({L}_{\mathrm{AGN},\mathrm{IR}}/\mathrm{erg}\,{{\rm{s}}}^{-1})\rangle =44.9$ for ERGs and $\langle \mathrm{log}({L}_{\mathrm{AGN},\mathrm{IR}}/\mathrm{erg}\,{{\rm{s}}}^{-1})\rangle =44.6$ for the NRGs. We therefore assume ERGs to be intrinsically radio-louder sources than NRGs and we will discuss this in Section 4.1.

The left panel of Figure 3 shows the distribution of our RGs in the SFR and M_⋆ plane, where both values are obtained from the optical+IR SED fitting with CIGALE (Toba et al. 2019). It is known that most star-forming galaxies follow the main-sequence (MS) whose normalization moves upward with redshift, as galaxies are more gas-rich at higher z (e.g., Brinchmann et al. 2004; Noeske et al. 2007; Elbaz et al. 2011; Whitaker et al. 2012; Schreiber et al. 2015). Above the MS, galaxies are referred to as starburst galaxies, producing stars more efficiently than regular star-forming galaxies. Pearson et al. (2018) measured the SFR and M_⋆ using multi-wavelength data from UV to far-IR employing the CIGALE code, and their MS is delineated with the shaded area in Figure 3, whose bottom and top boundaries correspond to redshifts of z = 0.3 and z = 1.6 with the scatter of σ = 0.35 at z = 0.3 and σ = 0.24 at z = 1.6, respectively. In this study, we define the sources as starburst galaxies if they are located above the yellow shaded MS area, which is roughly consistent with the criteria in previous studies (e.g., Elbaz et al. 2011).

Zoom In Zoom Out Reset image size

Considering that our sample spans a wide redshift range in 0.3 < z < 1.6 and their rapid redshift evolution of MS, redshift dependence is a key factor to understanding the relation in the left panel of Figure 3. Therefore, we also calculate the ratios between the observed SFRs and those predicted from the MS relations at a given stellar-mass, that is, SFR/SFR(MS), and its redshift dependence are shown in the right panel of Figure 3. We also compile how the SFR and M_⋆ relation evolves with redshift in Appendix B and Figure 10 (top panels).

Figure 3 shows several sequences of sources, notably at around sSFR ∼ 10^−10.5, ∼ 10⁻⁹, and ∼ 10⁻⁸ yr⁻¹, and the sources locate scarcely between the sequences. This does not originate from the nature properties, but more mainly from the complications of the limited bin numbers of the SED fitting parameters of star-forming history (SFH). In our case, the most notable contribution is considered to be f_burst, the mass fraction of the burst population. Toba et al. (2019) set a limited parameter binning of f_burst = 0.001, 0.1, and 0.3, which roughly corresponds to the bindings at sSFR ∼ 10^−10.5, 10⁻⁹, and 10⁻⁸ yr⁻¹. This binding might disappear once we follow the same SFH parameter binnings of Pearson et al. (2018), which was impossible for our studies because of the additional parameters of AGN and radio components which increase the computational time exponentially. Therefore, instead we assume that SFR in this study has a large uncertainty of ∼0.5 dex. Please see the more detailed discussion on the choice of the SFH parameters and its affects on the SFR estimations (e.g., Schreiber et al. 2015; Ciesla et al. 2017; Pearson et al. 2018).

Although our SFR estimation by Toba et al. (2019) might harbor large uncertainties as discussed above, Figure 3 shows three important suggestions. One is that most ERGs are distributed above the MS by a factor of ∼10 (see the right panel), suggesting that most ERGs are in starburst mode and they contain a copious amount of gas. This is indirect evidence for the availability of a cold gas supply fueling both star formation and the central SMBHs. The specific star formation rate defined by sSFR = SFR/M_⋆ reaches sSFR ∼ 10⁻⁸ yr⁻¹ for ERGs, suggesting that ERGs are in a rapid stellar-mass assembly phase with mass doubling times of only ∼100 Myr. This value is also comparable to one starburst duration (e.g., McQuinn et al. 2010) and thus most of the stellar component is expected to be dominated by a very young population.

Second, the onset of AGN feedback on the host galaxies has apparently not (yet) happened to most ERGs considering their location above the MS. This indicates a different picture from the conventional local RGs, which show a preference for massive host galaxies (M_⋆ > 10¹¹ M_⊙) and fall below the MS because their strong AGN feedback quenches star formation (e.g., Seymour et al. 2007). Note that our sample requires IR detections for the CIGALE SED fitting to derive stellar mass and SFR, therefore, the WERGS not detected in IR do not appear in Figure 3. In addition to that, our sample requires radio-compact morphologies, which removes a large population of local, radio-extended RGs. These two selection criteria may introduce a selection bias since it is possible that these IR non-detection or radio-extended sources are clustered below the SF main sequence. It is nonetheless intriguing that ERGs are located at the top edge of the distribution within the WERGS sample with IR detections.

The third point is that, if we limit ourselves to the sources with z < 1 (red star points), most of these are clustered at the low-M_⋆ end, typically at $\mathrm{log}{M}_{\star }/{M}_{\odot }\lt 10$. This suggests that the low-z subsample of ERGs is potentially important with regard to the study of BH seeds. Considering that all of our sources fulfill (1) L_{1.4 GHz} > 10²⁴ WHz⁻¹, that is, they lie above the threshold where their radio luminosity could be produced by star formation, as discussed in Section 2.1.2, and (2) the radio-excess selection with q_IR < 1.68, their radio emission can only be achieved by the AGN jet, requiring the existence of an SMBH at the center. Therefore, ERGs contain massive BH candidates residing in low-mass galaxies.

Figure 4 shows the distribution of ERGs and NRGs in the M_⋆ and L_{1.4 GHz} plane and illustrates the third point above from a different aspect. ERGs populate regions of comparatively small stellar mass and high L_{1.4 GHz}, showing that higher ${{ \mathcal R }}_{\mathrm{rest}}$ sources tend to have smaller M_⋆ as well as higher L_{1.4 GHz}. The combined results from above and from Figure 1 suggest that smaller M_⋆ sources tend to be observed as optically fainter sources and are therefore observed as high ${{ \mathcal R }}_{\mathrm{rest}}$ sources. We will discuss this point later in Section 4.1.

Zoom In Zoom Out Reset image size

Considering that most ERGs are very faint ($\left\langle {g}_{\mathrm{AB}}\right\rangle =24.5$) in optical, as shown in Figure 1 and Section 3.1.1, obtaining the BH virial mass from spectroscopy (e.g., McLure & Dunlop 2004) is time-consuming and difficult even if they are type-1 AGNs. This means that the direct measurement of the Eddington ratio (λ_Edd ≡ L_AGN,bol/L_Edd, where L_Edd is Eddington luminosity L_Edd ≃ 1.26 × 10³⁸(M_BH/M_⊙) erg s⁻¹), is difficult at this stage. Instead, we investigate the SMBH properties through the specific black hole accretion rate (hereafter, sBHAR), which may be considered a proxy for the Eddington ratio. The sBHAR is conventionally defined as sBHAR = L_AGN,bol/M_⋆ $\mathrm{erg}\,{{\rm{s}}}^{-1}\,{M}_{\odot }^{-1}$ (e.g., Mullaney et al. 2012). The left panel of Figure 5 shows the distribution of sBHAR for ERGs and NRGs. The median sBHAR of ERGs is $\langle \mathrm{log}(\mathrm{sBHAR}/\mathrm{erg}\,{{\rm{s}}}^{-1}\,{M}_{\odot }^{-1})\rangle =35.3$, which is an order of magnitude higher than that of the NRGs with $\langle \mathrm{log}(\mathrm{sBHAR}/\mathrm{erg}\,{{\rm{s}}}^{-1}\,{M}_{\odot }^{-1})\rangle =34.0$. Based on Best & Heckman (2012), the boundary of radiatively inefficient and efficient AGNs is at an Eddington ratio of $\mathrm{log}{\lambda }_{\mathrm{Edd}}\sim -2$, which corresponds to $\mathrm{log}(\mathrm{sBHAR}/\mathrm{erg}\,{{\rm{s}}}^{-1}\,{M}_{\odot }^{-1})=32$–33. Thus, almost all of our sources are considered to be radiatively efficient AGNs.

Zoom In Zoom Out Reset image size

This difference in sBHAR between the two subgroups is more pronounced in the right panel of Figure 5, showing sBHAR as a function of M_⋆. To illustrate the dependence of M_⋆ and λ_Edd on sBHAR more clearly, we use local scaling relations between M_BH and M_⋆. Previous studies assumed a constant ratio of M_BH/M_bulge = (1.5 − 2.0) × 10⁻³ (Marconi & Hunt 2003; Häring & Rix 2004) and M_bulge ≈ M_⋆ for simplicity (Aird et al. 2012, 2019; Mullaney et al. 2012; Delvecchio et al. 2018). In this study, we apply stellar-mass dependent values for M_BH with ${M}_{\mathrm{BH}}\propto {M}_{\star }^{\beta }$, where β = 1.4 (KH13; Kormendy & Ho 2013) and β = 1.05 (RV15; Reines & Volonteri 2015). With this relation, sBHAR can be written as a function of λ_Edd and M_⋆ by

$$\begin{eqnarray}&&\mathrm{sBHAR}\ (\mathrm{KH}13)=4.6\times {10}^{35}{\lambda }_{\mathrm{Edd}}{\left({M}_{\star }/{10}^{10}\,{M}_{\odot }\right)}^{0.4},\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&\mathrm{sBHAR}\ (\mathrm{RV}15)=3.2\times {10}^{34}{\lambda }_{\mathrm{Edd}}{\left({M}_{\star }/{10}^{10}\,{M}_{\odot }\right)}^{0.05}.\end{eqnarray} \tag{ 5 }$$

The right panel of Figure 5 shows that ERGs are more clustered in the high-Eddington-ratio regime. The expected median Eddington ratio by using Equation (4) is $\left\langle \mathrm{log}{\lambda }_{\mathrm{Edd}}\right\rangle \approx -0.4$. Some sources appear to exceed the Eddington limit. NRGs cover a broader range in sBHAR and some massive galaxies with M_⋆ > 10¹¹ M_⊙ fall into the radiatively inefficient regime with λ_Edd < −2.

We caution that the above discussion does not take into account the redshift evolution of the M_BH–M_⋆ relation, which has been suggested in some studies (e.g., Woo et al. 2008; Merloni et al. 2010; Ding et al. 2020). In this case, the estimated M_BH becomes larger at a given M_⋆, and the resulting Eddington ratio becomes smaller. Assuming the evolutional trend of Merloni et al. (2010), with M_BH/M_⋆(z) ∝ (1 + z)^0.68, and using a median redshift for the ERGs of 〈z〉 = 1.1, the normalization of sBHAR becomes slightly smaller by a factor of 1.6. As shown in Figure 5, this difference does not change our overall result. Although there are other significant uncertainties and possible systematics of the expected λ_Edd values, under the basic assumptions of Equations (4) or (5), our results suggest that higher ${{ \mathcal R }}_{\mathrm{rest}}$ sources have higher sBHAR and ERGs might contain sources achieving super-Eddington accretion. We will discuss this point in Section 4.1.

In summary, our results indicate that ERGs exhibit high sBHAR with the expected Eddington ratio of $\left\langle \mathrm{log}{\lambda }_{\mathrm{Edd}}\right\rangle \approx -0.4$ and some ERGs may be experiencing a super-Eddington phase. In addition, ERGs are starburst galaxies and their stellar mass is relatively lower than NRGs. ERGs appear to be in a phase preceding the onset of AGN feedback. If we limit the sample to z < 1, most ERGs are low-mass galaxies with M_⋆ < 10¹⁰ M_⊙. Therefore, lower-z ERGs are also candidates of massive, rapidly growing BHs.

Figures 1 and 2 indicate that ERGs or high ${{ \mathcal R }}_{\mathrm{rest}}$ sources tend to trace optically faint but radio-loud AGNs. In addition, Figure 4 shows that higher ${{ \mathcal R }}_{\mathrm{rest}}$ tends to pick up relatively lower stellar-mass populations. This means that the selection of ERGs using ${{ \mathcal R }}_{\mathrm{rest}}$ has a preference of selecting smaller M_⋆ and, therefore, ${{ \mathcal R }}_{\mathrm{rest}}$ using optical bands does not seem to trace AGN accretion disk emission in the optical band. This is not a new result, but it is a long-standing issue surrounding the use of the optical radio loudness ${{ \mathcal R }}_{\mathrm{rest}}$: it is strongly affected by host-galaxy contamination and/or the nuclear obscuration around the SMBHs (Terashima & Wilson 2003; Ho 2008). This raises the question of whether ERGs are intrinsically radio-loud or not.

To investigate this, we estimate the "intrinsic" radio loudness ${{ \mathcal R }}_{\mathrm{int}};$ the energy balance between the jet and accretion-disk defined by ${{ \mathcal R }}_{\mathrm{int}}={L}_{\mathrm{jet}}/{L}_{\mathrm{AGN},\mathrm{bol}}$, which is frequency-independent and also free from the contamination of host-galaxy components unlike ${{ \mathcal R }}_{\mathrm{rest}}$ (e.g., Terashima & Wilson 2003). For the accretion-disk luminosity, we use the bolometric AGN luminosity L_AGN,bol obtained from the IR AGN luminosity (L_AGN,IR) since L_AGN,IR is derived from IR SED decomposition, which traces the re-radiation of the dust powered by AGNs and is mostly insusceptible to absorption (Gandhi et al. 2009; Ichikawa et al. 2012, 2017, 2019; Asmus et al. 2015).

The total jet power of the AGNs (L_jet) can be estimated from the radio luminosity of the compact core or from the total radio emission (Willott et al. 1999; Cavagnolo et al. 2010; O'Sullivan et al. 2011). We adopt the relation between L_jet and L_{1.4 GHz} of Cavagnolo et al. (2010),

$$\begin{eqnarray}&&{L}_{\mathrm{jet}}=7.3\times {10}^{43}\,{({L}_{1.4\,\mathrm{GHz}}/{10}^{24}\,{\rm{W}}\,{\mathrm{Hz}}^{-1})}^{0.70}\,\mathrm{erg}\,{{\rm{s}}}^{-1}.\end{eqnarray} \tag{ 6 }$$

Note that the above equation is derived with the assumption that most radio sources have a spectral index of α = –0.8. This value is almost consistent with the median value of our WERGS sample; α = –0.65 for the sources with detections at both TGSS 150 MHz and VLA/FIRST 1.4 GHz or the fiducial value of α = –0.7 for the sources whose radio band detection is only in one band (i.e., VLA/FIRST 1.4 GHz only) (Condon 1992). Shabala & Godfrey (2013) recently investigated the effect of the radio source size in the jet-power estimation and found the dependence on the source size to be ∝ D^0.58, suggesting that the effect of the source size is small. Since our sample is selected only for radio-compact emission in VLA/FIRST with a spatial resolution of ∼5'' (or ≲ 20 kpc at z ∼ 1) this again mitigates the size dependence; otherwise, the sources are ultra-compact sources.

The left panel of Figure 6 shows the distribution of L_jet, spanning $44.0\lt \mathrm{log}({L}_{\mathrm{jet}}/\mathrm{erg}\,{{\rm{s}}}^{-1})\lt 46.4$ and the median values are $\left\langle \mathrm{log}{L}_{\mathrm{jet}}/\mathrm{erg}\,{{\rm{s}}}^{-1}\right\rangle =45.4$ for ERGs and $\left\langle \mathrm{log}{L}_{\mathrm{jet}}/\mathrm{erg}\,{{\rm{s}}}^{-1}\right\rangle =44.6$ for NRGs, which is equivalent to the radio-loud quasar level (e.g., Best et al. 2005; Inoue et al. 2017). This is naturally expected since we apply the simple conversion relation (Equation. 6). All of our sources are selected based on VLA/FIRST detection with the > 1 mJy flux limit and the radio luminosity range is similar to that of the FIRST-SDSS sample, although they are located in relatively higher-z (see Figure 2).

Zoom In Zoom Out Reset image size

This high L_jet is powerful enough to produce a radio jet with a size of ≳10 kpc that can disperse the interstellar medium (e.g., Nesvadba et al. 2017) and create cavities in the galactic interstellar medium and larger scale intergalactic environment to possibly quench star formation in host galaxies (e.g., McNamara & Nulsen 2007). Given the compact radio sizes of the sources in our sample (<5'' or ≲ 20 kpc at z ∼ 1) and given the fact that the host galaxies are still in an active star-forming phase (Figure 3), ERGs and most WERGS sources above or on the SF main sequence may still be in the phase before the onset of huge kinetic radio feedback.

The middle panel of Figure 6 shows the distribution of ${{ \mathcal R }}_{\mathrm{int}}$ and the median values of the two populations are $\left\langle \mathrm{log}{{ \mathcal R }}_{\mathrm{int}}\right\rangle =0.0$ for ERGs, and $\left\langle \mathrm{log}{{ \mathcal R }}_{\mathrm{int}}\right\rangle =-0.5$ for NRGs. For half of the ERGs, ${{ \mathcal R }}_{\mathrm{int}}\gt 1$, that is, the jet kinetic power is higher than the radiation from the accretion disk. The physical properties of higher ${{ \mathcal R }}_{\mathrm{rest}}$ sources are more clearly illustrated in Figure 7, which shows the distribution of ERGs and other WERGS sources in the ${{ \mathcal R }}_{\mathrm{int}}$–sBHAR plane. The yellow shaded area is the region corresponding to λ_Edd = 1 assuming Equation (4) for the range of $8\lt \mathrm{log}({M}_{\star }/{M}_{\odot })\lt 12$. The other WERGS sources are shown with a color gradation in $\mathrm{log}{{ \mathcal R }}_{\mathrm{rest}}$, with higher ${{ \mathcal R }}_{\mathrm{rest}}$ lying toward the top right of the plot. This means that higher ${{ \mathcal R }}_{\mathrm{rest}}$ sources tend to have both higher λ_Edd and higher ${{ \mathcal R }}_{\mathrm{int}}$. It is well known that radio galaxies are more radio-loud as the Eddington-scaled accretion rate decreases, producing the sequence from the top left region toward the bottom right one in Figure 7 (Ho 2002; Merloni & Heinz 2007; Panessa et al. 2007; Sikora et al. 2007; Sikora et al. 2013; Ho 2008), and overall NRGs also follow this trend. However, ERGs are located in the top right part of the plane, suggesting that ERGs contain rapidly growing BH in the center, and at the same time harbor powerful, compact radio jets.

Zoom In Zoom Out Reset image size

These radio galaxies are different from the conventional radio galaxies in the local universe (Seymour et al. 2007). Delvecchio et al. (2018) recently found that high-z (z > 1.5) radio galaxies discovered in the VLA/COSMOS survey tend to have high Eddington ratios with $\mathrm{log}{\lambda }_{\mathrm{Edd}}\gt -2$ and that they are also star-forming or starburst galaxies. This is consistent with the properties of our overall WERGS sample (see also Toba et al. 2019) and thus ERGs occupy the extreme end of highly accreting BHs in radio galaxies.

Figure 7 indicates that ERGs are a very interesting population, showing both radiative (high λ_Edd) and jet-efficient (high ${{ \mathcal R }}_{\mathrm{int}}$) emission, which were missed in previous surveys. To investigate the jet properties in more detail, we here discuss the jet efficiency of ERGs and compare the obtained values with other radio sources. The jet production efficiency is defined as

$$\begin{eqnarray}&&{\eta }_{\mathrm{jet}}={L}_{\mathrm{jet}}/{\dot{M}}_{\mathrm{BH}}{c}^{2}={\eta }_{\mathrm{rad}}{L}_{\mathrm{jet}}/{L}_{\mathrm{AGN},\mathrm{bol}}\end{eqnarray} \tag{ 7 }$$

where η_rad is the radiation efficiency of an AGN accretion disk, c is the speed of light, and ${\dot{M}}_{\mathrm{BH}}$ is the mass accretion rate onto the SMBH through the disk. In the following, we adopt a canonical value of η_rad = 0.1 based on the Soltan–Paczynski argument (Soltan 1982).

The distribution of η_jet is shown in the right panel of Figure 6. The median values are $\left\langle \mathrm{log}{\eta }_{\mathrm{jet}}\right\rangle =-1.0\pm 0.3$ for ERGs and $\left\langle \mathrm{log}{\eta }_{\mathrm{jet}}\right\rangle =-1.5\pm 0.4$ for NRGs. These values are slightly higher or consistent with the nearby radio AGN population ($\mathrm{log}{\eta }_{\mathrm{jet}}\sim -1.5$) whose host galaxies are massive (M_⋆ > 10¹¹ M_⊙; Nemmen & Tchekhovskoy 2015). However, the local radio AGN population generally has much smaller Eddington ratios of $\mathrm{log}{\lambda }_{\mathrm{Edd}}\lt -2$. This population would be located in the top left region in Figure 7, so they are clearly different from ERGs. The η_jet of ERGs is higher by a factor of ∼10 compared to the SDSS radio quasar population, whose radio jet efficiency lies at $\left\langle \mathrm{log}{\eta }_{\mathrm{jet}}\right\rangle =-2.0\pm 0.5$ (e.g., van Velzen & Falcke 2013; Inoue et al. 2017). Considering that the SDSS radio quasars on average show $\mathrm{log}{\lambda }_{\mathrm{Edd}}\sim -1$, those radio quasars would be located at the bottom center to bottom right of Figure 7, which is again different from the distribution of ERGs. Therefore, the disk–jet connection may be fundamentally different from the standard disk model, which conventionally describes the quasar/AGN accretion disk emission. One possible origin of this high production efficiency of radiation and jets is that the accretion disk of ERGs is actually in the "radiatively inefficient" state, but the physical origin of radiative inefficiency is different from the disks in the local radio galaxies. For example, the ERGs may be undergoing more rapid mass accretion (the so-called slim disk model; Abramowicz et al. 1988). Recent radiation hydrodynamical simulations suggest that when the mass accretion rate significantly exceeds the Eddington rate, radiation is effectively trapped within the accreting matter and advected to the central BH before escaping by radiative diffusion. As a result, the emergent radiation luminosity is saturated at the order of L_Edd (i.e., η_rad ≲ 0.1 at $\gg {\dot{M}}_{\mathrm{Edd}}$) and the accretion flow turns into a radiatively inefficient state even with a super-Eddington accretion rate (e.g., Ohsuga et al. 2005, 2009; McKinney et al. 2014; Inayoshi et al. 2016; Takeo et al. 2020). This phase is considered to be a key process of the BH seed growth in the high-z (z > 6) universe to describe the already known massive high-z SMBHs (e.g., Mortlock et al. 2011; Wu et al. 2015; Bañados et al. 2018; Onoue et al. 2019).

In the situation where the accretion rate is well above the Eddington rate, a relativistic jet can be launched by the Blandford–Znajek mechanism if large-scale magnetic fields exist in the innermost region (Blandford & Znajek 1977; Tchekhovskoy et al. 2011). The jet behavior of these extreme cases has been investigated by several authors (Sadowski et al. 2014; Sadowski & Narayan 2015) using a 3D general relativistic radiation magneto-hydrodynamical simulation (GRRMHD) code. They find that, in this (super-)Eddington accretion phase, BHs emit a significant amount of radiation and jet energy and both the radiation and jet power approach around the Eddington limit, which appears to be the case for ERGs and the sources in the top right region in Figure 7 (see also Blandford et al. 2019). Thus, the study of ERGs at z ∼ 1 will be complementary to studies utilizing the next generation instruments capable of observing seed BHs in the early universe at z ≳ 7 (Haiman et al. 2019) in the investigations of the growth of seed BHs.

The discussion related to the Eddington ratio depends on the M_BH estimation, which in turn relies on the local scaling relation of Kormendy & Ho (2013). However, it is not obvious that this relation also holds in the low-mass end in which many ERGs probably reside. Although the current data cannot give any constraints on this issue and a more detailed study is beyond the scope of this paper, we discuss two alternative scenarios in which our galaxies are above (case one) or below (case two) the correlation between M_BH and M_⋆ as shown in Figure 8.

Zoom In Zoom Out Reset image size

Case one reflects the possibility that ERGs might host overmassive BHs compared to the relation by Kormendy & Ho (2013), with the BH-to-stellar-mass ratio f_b,⋆ ≡ M_BH/M_⋆ ≳ 10⁻² (see the orange shaded area in Figure 8). In this case, concerns surrounding super-Eddington accretion discussed in Sections 3.3 and 4.1 are alleviated: the median Eddington ratio for ERGs can be parameterized by $\langle {\lambda }_{\mathrm{Edd}}\rangle \simeq 0.17\,{({f}_{{\rm{b}},\star }/0.01)}^{-1}$.

In addition, the possible future growth path in the plane of M_BH and M_⋆ can be discussed by assuming the current BH/stellar mass assembly rate. Since the radiation efficiency is η_rad = 0.1 in the sub-Eddington regime, the BH accretion rate is estimated as ${\dot{M}}_{\mathrm{BH}}={L}_{\mathrm{AGN},\mathrm{bol}}/({\eta }_{\mathrm{rad}}{c}^{2})\simeq 1.7\,\times \,{10}^{-3}{f}_{{\rm{b}},\star }^{-1}{\dot{M}}_{\mathrm{Edd}}$, where the median Eddington ratio above is substituted (note that this equation is valid at f_b,⋆ ≳ 1.7 × 10⁻³). Comparing the BH accretion rate to the sSFR of ERGs (sSFR = SFR/M_⋆ ≃ 10⁻⁸ yr⁻¹; see Figure 3), we obtain ${\dot{M}}_{\mathrm{BH}}/\mathrm{SFR}\simeq 4\times {10}^{-3}$, which is independent of the choice of f_b,⋆. Therefore, the BH-galaxy mass ratio for ERGs is expected to approach a canonical value observed in the local universe (e.g., f_b,⋆ ∼ 3 × 10⁻³; Kormendy & Ho 2013). This evolutional path might be analogous to that of known high-z luminous quasars whose SMBHs are likely overmassive compared to the local relation with 10⁻² ≲ f_b,⋆ ≲ 2 × 10⁻¹ (e.g., Wang et al. 2013, 2016; Trakhtenbrot et al. 2015; Decarli et al. 2018).

Case two considers the possibility that ERGs might host very undermassive BHs with f_b,⋆ ≲ 10⁻³–10⁻⁴ as shown in the blue shaded region in Figure 8. For instance, the overall normalization for late-type galaxies may be considerably lower than for the local relation of early-type galaxies (e.g., Reines & Volonteri 2015; Läsker et al. 2016; Greene et al. 2020). This "galaxy grows first, BHs come later" phase may occur because gas accretion onto the BHs is quite inefficient due to efficient stellar feedback, e.g., star formation-driven outflows, and the shallow halo potential until the stellar-mass reaches a critical value of M_⋆ ∼ 10^10.5 M_⊙ (e.g., Bower et al. 2017; Habouzit et al. 2017). However, once the galaxy reaches the critical mass, above which star formation-driven outflows no longer prevent gas accretion onto the galactic nuclei, they abruptly switch to a rapid growth phase in the absence of feedback, and finally merge into the local relation. ERGs with M_⋆ ≳ 10^10.5 M_⊙ might be in such a rapidly growing phase.

The scenario above looks promising but it is worth discussing the critical mass value because the stellar mass of most ERGs actually lies below the expected critical mass. Beside, Figure 5 shows that the median sBHAR (and therefore λ_Edd) tends to be higher for less massive ERGs with the median Eddington ratio of $\langle {\lambda }_{\mathrm{Edd}}\rangle \simeq 1.7\,{({f}_{{\rm{b}},\star }/{10}^{-3})}^{-1}$, suggesting that the extreme BH growth is possible (at least temporally) even in lower mass galaxies. This discrepancy may reflect that there is a missing understanding about the physics of AGN feedback and gas feeding in low-mass galaxies.

Another interesting point in the second scenario is that the Eddington ratio for several ERGs becomes as high as λ_Edd ∼ 10². Recent numerical simulations suggest that for a highly magnetized accretion disk around a rapidly spinning BH, the disk transits into a magnetically arrested disk (MAD) state and produces high radiative luminosity (e.g., McKinney et al. 2014; Sadowski & Narayan 2015). Through the emission process, the level of high luminosity could be achieved only when the BH mass accretion exceeds ∼${10}^{3}\,{\dot{M}}_{\mathrm{Edd}}$ (Inayoshi et al. 2016; see also Figure 5 in Inayoshi et al. 2020). Therefore, if this scenario is true, future follow-up observations of those ERGs with such high Eddington ratios would reveal the properties of possibly hyper-Eddington accreting BHs (Takeo et al. 2020).

In recent years, the theoretical study of AGN feedback has been extended to low-mass galaxies (Silk 2017). This was motivated mainly by the increasing number of intriguing observational results showing signs of AGN feedback in low-mass galaxies (Nyland et al. 2017; Penny et al. 2018; Dickey et al. 2019; Kaviraj et al. 2019; Mezcua et al. 2019). Our results show that the jet power of L_jet > 10⁴⁴ erg s⁻¹ in ERGs is equivalent to that produced from radio quasars and local radio AGNs residing in massive galaxies. This suggests the presence of a strong jet that can disperse or even expel the interstellar medium (e.g., Morganti et al. 2005; Holt et al. 2008; Nesvadba et al. 2017) or produce cavities in the host galaxies, which are often observed in the local massive counterparts (e.g., Rafferty et al. 2006; McNamara & Nulsen 2007; Bîrzan et al. 2008; Blandford et al. 2019). On the other hand, low stellar-mass ERGs are still in a starburst phase, indicative of a substantial gas reservoir, and lacking clear signatures of negative AGN feedback. ¹⁵

To investigate whether such a negative feedback from the jet is effective in the host galaxies, it is necessary to obtain high spatial resolution radio images that resolve the host galaxy down to scales of <10 kpc (Reines et al. 2020). Our sample of compact radio sources from VLA/FIRST with a spatial resolution of ∼5'' only provide poor upper-bounds of the radio jet size of <40 kpc at z ∼1. Future higher quality observations with subarcsecond resolution will provide more insight into the effect of the jet on the host galaxy.

Given that the BH mass of ERGs might reach the massive BH range below M_BH ∼ 10⁶ M_⊙, studying the properties of ERGs gives a window into understanding seed BH growth in the high-z universe. Some observational studies recently suggested that AGN feedback might also have a positive effect on the host galaxies and intensive in situ star formation is ignited in the outflowing gas launched by the AGN (Maiolino et al. 2017). Recent numerical simulations of jet-driven feedback also predict star formation triggered on galaxy scales by the overpressure of the jet during the first few Myrs of strong interaction between the jet and the interstellar medium (Wagner & Bicknell 2011; Gaibler et al. 2012; Wagner et al. 2016; Mukherjee et al. 2018). If the radio jets in low-mass galaxies have a positive impact on both star formation and BH growth, then radio-selected AGNs in low-mass galaxies are an intriguing population for understanding the BH and stellar mass assembly in the environment of shallow potential low-mass galaxies.

The small number density of ERGs (∼1 source deg⁻²) and their optical faintness of i_AB ∼ 25 suggest that ERGs were easily missed in previous optical surveys of radio sources. Our work shows that ERGs can be discovered with wide (>100 deg²) and deep (i_AB ≲ 26) optical follow-up observations of radio sources. The number of known ERGs will increase in the near future, once the Subaru/SSP survey is complete, covering an area of ∼10³ deg². The upcoming Subaru/Prime Focus Spectrograph (PSF; Takada et al. 2014) will conduct extensive optical and near-IR spectroscopic follow-up observations in the footprint of the HSC-SSP field, providing spec-z information of ERGs with i_AB ≲ 24. The forthcoming LSST survey (Ivezić et al. 2019) will cover half of the sky, and it will also increase the number of sources by another order of magnitude, resulting in ∼10⁴ sources. A significant fraction of these ERGs will also be expected to reside in low-mass galaxies. Therefore, the combination of wide area (either shallow or deep) radio and deep optical imaging and spectroscopic surveys will give us a statistically large number of massive BH candidates in the coming years.

Massive BH searches have been conducted with multi-wavelength data; optical spectroscopy using broad Hα emission (Greene & Ho 2004, 2007; Xiao et al. 2011) or narrow emission-line ratios (e.g., Barth et al. 2008; Reines et al. 2013), mid-IR (W1 − W2) colors in WISE (Satyapal et al. 2014; Sartori et al. 2015; Secrest et al. 2015; Kaviraj et al. 2019), X-ray observations (e.g., Schramm et al. 2013; Baldassare et al. 2015; Mezcua et al. 2016; Chen et al. 2017; Kawamuro et al. 2019), as well as through studies of the AGN flux variability based on the properties that lower luminosity AGN have stronger variability amplitude (Morokuma et al. 2016; Baldassare et al. 2018, 2020; Elmer et al. 2020; Kimura et al. 2020). However, radio surveys have rarely been used for searching massive BHs in low-mass galaxies.

Recently, Reines et al. (2020) employed the radio-survey approach to search for massive BHs by using VLA/FIRST and conducting VLA high spatial resolution follow-up observations of local dwarf galaxies at z < 0.1 and found more than 10 candidates. One remarkable result was that most radio cores were not located at the center of the host, but off-nucleus, a possible signature of a previous merger. If the detection of an off-nucleus core is indeed the consequence of a galaxy merger, it indicates that the orbital decay of a black hole through dynamical friction is inefficient in this BH mass range, resulting in a population of BHs that fail to reach the galactic center within a Hubble time since major mergers in dwarf galaxies becomes rare at z < 3 (Fitts et al. 2018). Currently there are many theoretical implications regarding wandering BHs (Comerford & Greene 2014; Blecha et al. 2016; Tremmel et al. 2018; De Rosa et al. 2019). High spatial resolution radio imaging can pin down the locations of such wandering BHs, and large-scale statistical studies may provide insights into the frequency and environment of wandering BHs through the combination of LSST and VLA/FIRST or the ongoing VLASS (Lacy et al. 2020), which will achieve a sensitivity down to 0.1 mJy at 2–4 GHz, giving an order of magnitude deeper radio observations than VLA/FIRST.

It could be argued that blazars might contaminate the sample of ERGs, and that therefore the optical emission is dominated by the jet component. Based on the known blazar sequence in the radio and optical bands (Fossati et al. 1998; Donato et al. 2001; Inoue & Totani 2009; Ghisellini et al. 2017), the observed radio loudness of blazars is in the range of $\mathrm{log}{{ \mathcal R }}_{\mathrm{obs}}\sim 1.5\mbox{--}3$. This range is obtained by using the luminosity range of L_{1.4 GHz} = 10⁴¹⁻⁴³ erg s⁻¹ and shifting redshift to the range of z ∼ 0–5. Therefore, blazars cannot reproduce extreme radio loudness reaching $\mathrm{log}{{ \mathcal R }}_{\mathrm{rest}}\gt 4$. Blazar contamination is unlikely also from the point of the observed optical magnitude. Considering that all ERGs are VLA/FIRST selected with a flux limit of >1 mJy, the expected observed magnitude of blazars based on the blazar sequence is high with a peak around i_AB = 21, and all cases are at i_AB < 23.5 even when changing the luminosity range of the blazar SEDs and also shifting the redshift range to z = 0–5. Based on these arguments, the contamination of blazars on ERGs is unlikely.