Enhanced star formation in z ∼ 6 quasar companions

ABSTRACT

Quasars powered by supermassive black holes (>10⁸ M_⊙) at z ∼ 6 are predicted to reside in cosmic overdense regions. However, observations so far could not confirm this expectation due to limited statistics. The picture is further complicated by the possible effects of quasar outflows (i.e. feedback) that could either suppress or stimulate the star formation rate (SFR) of companion galaxies, thus modifying the expected bias. Here, we quantify feedback effects on the properties and detectability of companions by comparing cosmological zoom-in simulations of a quasar in which feedback is either included or turned-off. With respect to the no-feedback case, companions (a) directly impacted by the outflow have their SFR increased by a factor of 2−3, and (b) tend to be more massive. Both effects shift the [C ii] 158 μm and UV luminosity functions towards brighter magnitudes. This leads us to conclude that quasar feedback slightly increases the effective quasar bias, boosting the number density of observable quasar companions, in agreement with what has been found around the brightest quasars of recent Atacama Large Millimeter Array (ALMA) [C ii] surveys. Deeper observations performed with James Webb Space Telescope and/or ALMA will improve the statistical significance of this result by detecting a larger number of fainter quasar companions.

methods: numerical, galaxies: evolution, galaxies: high-redshift, galaxies: star formation, (galaxies:) quasars: general, (galaxies:) quasars: supermassive black holes

1 INTRODUCTION

Massive black holes (MBHs) are thought to be ubiquitous in the centre of massive galaxies at all redshifts (Kormendy & Ho 2013, and references therein). Through gas accretion, they are the engines responsible for the huge emission of galactic centres as active galactic nuclei (AGNs). Extremely bright quasars are AGNs with enormous luminosities (L_bol > 10⁴⁶ erg s⁻¹), likely powered by the most massive MBHs, with masses M_BH ranging from 10⁸ to 10¹⁰ M_⊙.

The number of observed bright quasars at z > 6 has substantially grown in the last years (e.g. Carnall et al. 2015; Jiang et al. 2015; Reed et al. 2015; Venemans et al. 2015; Wu et al. 2015; Matsuoka et al. 2016; Reed et al. 2017; Bañados et al. 2018), proving that MBHs with M_BH ∼ 10⁹ M_⊙ were already in place when the Universe was less than 950 Myr old.

The fast growth of these compact objects represents an exciting and fully open problem in modern astrophysics and various models have been proposed to solve the issue. For a long time, it was assumed that MBH seeds were provided by (i) the massive remnants of Pop III stars (see e.g. Madau & Rees 2001; Haiman & Loeb 2001; Heger et al. 2003; Volonteri, Haardt & Madau 2003; Madau et al. 2004), formed at z ≳ 20. However, this scenario is problematic, since the accretion would have to either continue unperturbed at the Eddington rate, maintaining a very low radiative efficiency (not to excessively hinder the feeding process; see e.g. Tanaka & Haiman 2009), or (ii) exceed the Eddington limit (see e.g. Madau, Haardt & Dotti 2014; Volonteri, Silk & Dubus 2015; Lupi et al. 2016; Pezzulli et al. 2017). Other models predict the (iii) existence of very massive MBH seeds (M_BH ≳ 10⁴ M_⊙, possibly originated from direct collapse; see e.g. Loeb & Rasio 1994; Bromm & Loeb 2003; Koushiappas, Bullock & Dekel 2004; Spaans & Silk 2006; Mayer et al. 2010, 2015).

The steady and efficient gas inflow, necessary for the fast MBH growth, would be strongly favoured within highly massive DM haloes (see Efstathiou & Rees 1988). Accordingly, theoretical models including numerical simulations with MBH accretion and feedback prescriptions (see e.g. Barkana & Loeb 2001; Springel et al. 2005b; Sijacki, Springel & Haehnelt 2009; Di Matteo et al. 2012; Costa et al. 2014; Ni et al. 2018; Weinberger et al. 2018; Marshall et al. 2019; Ni et al. 2020) predict high-z MBHs to reside in the most massive DM-haloes with M_vir > 10¹²–10¹³ M_⊙, corresponding to fluctuations above 3−4σ in the cosmic density field (Barkana & Loeb 2001). A connection between MBHs and their host DM halo has been studied also via clustering measurements of AGN. Various large surveys, from the optical to the X-ray band, have shown that AGN reside in massive haloes with masses M_vir = 10¹²–10¹³ M_⊙ at low (z ∼ 0.1) and higher (up to z ≈ 2) redshift (e.g. Hickox et al. 2009; Ross et al. 2009; Cappelluti et al. 2010; Allevato et al. 2011, 2012). As a consequence, quasars at high redshift should be part of large-scale structures marked by a significant overdensity in galaxy number count that may extend over Mpc scales (e.g. Costa et al. 2014). However, this expectation has not been conclusively confirmed by observations.

On the one hand, various observations corroborated the high-density scenario hypothesis (see e.g. Steidel et al. 2005; Capak et al. 2011; Swinbank et al. 2012; Husband et al. 2013 for 2 < z < 5 and Kim et al. 2009; Morselli et al. 2014; Balmaverde et al. 2017; Decarli et al. 2017, 2018; Mignoli et al. 2020; Venemans et al. 2020, for z ∼ 6 and above). In particular, Venemans et al. (2020) followed up with the Atacama Large Millimeter Array (ALMA) a sample of 27 z ∼ 6 quasars, previously detected in [C ii], discovering 17 [C ii] bright galaxies and finding that some of the quasars present multiple (2–3) companions. Furthermore, whereas no dual AGN has been confirmed so far at z > 6 (Connor et al. 2019, 2020; Vito et al. 2019, 2021), several of them have been detected up to z ∼ 5 (see Koss et al. 2012; Vignali et al. 2018; Silverman et al. 2020).

On the other hand, more than a few observations did not reveal significant galaxy count overdensities around quasars (see e.g. Francis & Bland-Hawthorn 2004 and Simpson et al. 2014 at z ∼ 2, Uchiyama et al. 2018 at z = 4, Kashikawa et al. 2007 and Kikuta et al. 2017 at z ∼ 5, Bañados et al. 2013 and Mazzucchelli et al. 2017 at z > 5.7) and some theoretical studies consistently suggested that MBHs does not have to inevitably lie in the most massive haloes (see for instance Fanidakis et al. 2013; Orsi et al. 2016; Di Matteo et al. 2017; Habouzit et al. 2019).

These conflicting observational results could be explained by the difficulty in detecting the satellites in the neighbourhood of a luminous quasar, where the strong AGN feedback can play a fundamental role in shaping their baryonic component. AGN feedback could affect the MBH environment well beyond the host galaxy-scale radius and significantly modify the star formation (SF) activity of the orbiting companions (see Martín-Navarro, Burchett & Mezcua 2019; Martín-Navarro et al. 2021). Dashyan et al. (2019) explicitly investigated, in cosmological hydrodynamical simulations, the AGN-driven quenching effect within galactic satellites at low redshift (z < 3). The authors claim that AGN winds can decrease the SF process in AGN companions, by sweeping away their gas, out to five times the virial radius of the central galaxy. Conversely, Gilli et al. (2019) observed a radio AGN at z = 1.7, finding a possible positive effect of its feedback on the star formation of its companion, over a scale ≳ 450 kpc. Furthermore, Fragile et al. (2017) performed an isolated simulation in order to explain the positive feedback observed by Croft et al. (2006) in a star-forming galaxy, where the star formation is triggered by the radio-jets of the nearby NGC–541.

At higher redshift, the direct effect of AGN feedback on the satellite galaxies remains unclear. Several studies found, both through observations (Kashikawa et al. 2007), and numerically (Efstathiou 1992; Thoul & Weinberg 1996; Okamoto, Gao & Theuns 2008), that the ionizing background produced by the quasar is able to inhibit SF in satellites, or even suppress their assembling. In the most extreme cases, satellite haloes would still populate the quasar environment, without being detectable because of their low gas and stellar content, especially given the limited sensitivity of modern instruments. In this context, the upcoming James Webb Space Telescope (JWST; Gardner et al. 2006) will allow us to observe fainter AGN companions, hopefully alleviating our sensitivity bias. With a primary mirror of about 6.5 m, JWST reaches a sensitivity 3–5 times higher than that of Hubble Space Telescope at 1 μm, enabling the detection of the rest-frame UV emission that arises from distant faint star-forming galaxies with relative short exposure time.

Numerical simulations are powerful tools to study the AGN effect on the surrounding galaxies, for they provide the complete spatial and time distribution of matter, as well as fundamental self-consistent sub-grid models, such as SF, stellar feedback, cooling processes, etc. Recently, enormous steps forward have been made in this field and several works successfully portrayed the formation and evolution of bright quasar hosts at high redshift. For instance, Di Matteo et al. (2017) investigated the accretion efficiency of high-z MBHs, by employing advanced refinement techniques, focussing on the effects of AGN feedback on the host galaxy. Curtis & Sijacki (2016) and van der Vlugt & Costa (2019) studied how AGN feedback shapes the dynamical components of the galaxy. Costa et al. (2014) and Smidt et al. (2018) analysed how AGN X-ray luminosity and feedback affect negatively SF in the host system. Lupi et al. (2019) and Lupi et al. (2021) detailed the consequences of AGN thermal feedback on the host interstellar medium (ISM). Differently, Richardson et al. (2016), Barai et al. (2018), and Valentini, Gallerani & Ferrara (2021) studied the effect of AGN feedback on the formation and evolution of a proto-cluster, using either different numerical methods or different feedback prescriptions. Finally, Costa, Pakmor & Springel (2020) developed a state-of-the-art model to describe AGN-driven small-scale winds and tested it in isolated simulations.

In this work, we analyse the environment of a powerful quasar at z ≳ 6, resulting from the simulations by Barai et al. (2018). To assess the possible effect of the AGN feedback on the surrounding companion galaxies, we take advantage of the Barai et al. (2018) suite and compare the quasar environment with a control run in which MBHs are not seeded. In detail, our goal is to (i) evaluate the impact of quasar feedback on its environment, far beyond its host galaxy radius, and (ii) provide theoretical predictions on the expected UV and rest-frame FIR luminosities of the neighbour satellites.

This paper is organized as follows. In Section 2, we describe the numerical model adopted in this work and introduce the runs with and without AGN, hereafter called AGNcone and noAGN, respectively; in Section 3, we present the sample of satellites and statistically analyse their redshift evolution and how their properties [e.g. number of satellites, star formation rate (SFR), stellar mass, gas mass, metallicity] are affected by their position in the proto-cluster; in Section 4, we focus on the effect of AGN feedback on individual satellites, and present an interpretation of our results in Section 5; in Section 6, we discuss the observational properties of the galaxy group; finally, we summarize our findings and draw our conclusions in Section 7.

2 HYDRODYNAMIC SIMULATIONS

The two runs analysed in this work belong to a suite of zoom-in cosmological simulations (Barai et al. 2018) built to follow the formation of a massive galaxy proto-cluster at z ≃ 6 through the smoothed particle hydrodynamics N-body code gadget-3 (Springel 2005; Springel et al. 2008).

Both the simulations share the same initial conditions, generated through the code music (see Hahn & Abel 2011) and assume the same recipe for the sub-grid physics, with the exception of the MBH prescription. In particular, the cosmological parameter set refers to a flat ΛCDM Universe with Ω_M,0 = 0.3089, Ω_Λ,0 = 1 − Ω_M,0 = 0.6911, Ω_b,0 = 0.0486, and H₀ = 67.74 Mpc s⁻¹ (Planck Collaboration XIII 2016). The simulated box of side 500 comoving Mpc has been evolved with only DM particles from z = 100, till z ≲ 6, with an initial mass resolution of 2 × 10¹⁰ M_⊙ per particle and a softening length of 48.72 comoving kpc. There, the Lagrangian volume of the most massive proto-cluster has been identified (through a Friends-of-Friends algorithm) with a virial mass |$M_{\rm vir} \simeq 10^{12}\, \!{\rm M}_\odot$| and a comoving virial radius r_vir ≃ 511 kpc, traced back to z = 100, refined and re-simulated along with baryons. In the most refined region – set to be, originally, a cube of side 5.21 Mpc –, the mass resolution is given by m_DM = 7.54 × 10⁶ M_⊙ for 59 1408 DM particles, and m_gas = 1.41 × 10⁶ M_⊙ for the same number of gas particles, whereas the spatial resolution is set by the gravitational softening length of all particle species (ϵ ≃ 1.476 comoving kpc). The adaptive smoothing length is computed according to the standard prescription by Springel et al. (2008) and its minimum value is set to 0.001ϵ.

The whole suite implements radiative heating and cooling using the rates provided in the tables of Wiersma, Schaye & Smith (2009) in ionization equilibrium. Metal-line cooling is also considered. Eleven element abundances (H, He, C, Ca, O, N, Ne, Mg, S, Si, Fe) are followed according to the receipt of Tornatore et al. (2007) in the presence of a redshift-dependent cosmic ionizing background (Haardt & Madau 2012).

SF is modelled following the multiphase recipes by Springel & Hernquist (2003). More specifically, gas particles denser than n_SF = 0.13 cm⁻³ are converted into collisionless star particles according to the stochastic scheme of Katz, Weinberg & Hernquist (1996). Each spawned star particle represents a stellar population described by a Chabrier (2003) initial mass function in the mass range of 0.1−100 M_⊙. Stars are allowed to explode as supernovae (SN) releasing kinetic energy via a constant-velocity outflow with v_SN = 350 km s⁻¹ (see Barai et al. 2015; Biffi et al. 2016). Metal enrichment of the interstellar medium is provided by Type Ia SN (0.8 < M/M_⊙ < 8) considering a fraction of binary of 1/10, according to Thielemann et al. (2003), Type II SN (M/M_⊙ > 8) according to Woosley & Weaver (1995), and winds from asymptotic giant branch stars following van den Hoek & Groenewegen (1997).

The two simulations differ only in terms of the prescription used to describe massive black holes. In detail, in noAGN only cooling, metal enrichment, star formation, and SN feedback are included, with no prescription for MBHs.

By contrast, in the AGNcone run MBHs are represented as sink particles that can form and grow both via accretion of gas and through mergers with other MBHs. In particular, an MBH is seeded in a halo when (i) the halo does not host any other MBH and (ii) the halo virial mass is M_h ≥ 10⁹ M_⊙ (i.e. the halo is properly resolved). When these conditions are satisfied, an M_BH = 10⁵ M_⊙ MBH is placed at the centre of mass of the halo. The simulation implements a repositioning algorithm as in Springel, Di Matteo & Hernquist (2005a) and Schaye et al. (2015).

Every MBH can accrete mass from the surrounding medium via the classical Bondi–Hoyle–Lyttleton accretion rate (Hoyle & Lyttleton 1939; Bondi & Hoyle 1944; Bondi 1952):

$$\begin{eqnarray*} \dot{M}_{\rm Bondi} = \frac{4 \pi G^2 M^2_{\rm BH} \rho }{(c^2_{\rm s}+v^2)^{3/2}}, \end{eqnarray*}$$

(1)

where G is the gravitational constant, ρ is the gas density, c_s is the sound speed, and v is the gas relative velocity with respect to the MBH.¹ The accretion rate is multiplied by a boost factor of 100, analogously to what has been done, e.g. in Springel et al. (2005a) and it is capped to the Eddington limit,

$$\begin{eqnarray*} \dot{M}_{\rm Edd} = \frac{4 \pi G M_{\rm BH} m_{\rm p}}{\epsilon _{r} \sigma _{\rm T} c}, \end{eqnarray*}$$

(2)

where m_p is the proton mass, σ_T is the Thomson cross-section, and c is the speed of light in vacuum. ϵ_r is the radiative efficiency set equal to 0.1 (see the average efficiency for an optically thick and geometrically thin accretion discs by Shakura & Sunyaev 1973). During the accretion process an MBH radiates a fraction ϵ_r of the accreted rest-mass energy

$$\begin{eqnarray*} L_{\rm rad} = \epsilon _{\rm r} \dot{M}_{\rm acc} c^2, \end{eqnarray*}$$

(3)

where |$\dot{M}_{\rm acc}$| is the rate of the inflowing gas on to the MBH, and a fraction ϵ_f of this luminosity is coupled to the interstellar medium as feedback energy:

$$\begin{eqnarray*} \dot{E}_{f} = \epsilon _{\rm f} L_{\rm rad}, \end{eqnarray*}$$

(4)

where ϵ_f = 0.05 as in, e.g., Di Matteo et al. (2008). Feedback is kinetic and is modelled through the ejection, in a bi-cone with a half-opening angle of 45 deg, of a certain mass of gas M_w with a fixed initial velocity of v_w = 10⁴ km s⁻¹, such that the kinetic luminosity |$\frac{1}{2}\dot{M_{\rm w}}v_{\rm w}^2$| is equal to |$\dot{E}_{f}$|⁠. The direction of emission is random and it is associated with the MBH when it is seeded. We note that this assumption is supported by several studies showing little or no alignment between the outflow/jet axis and the large-scale angular momentum of the host galaxy (see e.g. Hopkins et al. 2012, and references therein).

To summarize, in the control run noAGN only SN feedback is included, whereas AGNcone incorporates both SN and AGN feedback. Different feedback prescriptions result into different host galaxy properties (see Barai et al. 2018, for a detailed discussion). Here, we focus our analysis on the quasar environment. The properties of the two runs at their last snapshot are outlined in Table 1.

Table 1.

Summary table for the two analysed runs at z = 6. From left to right: (i) name of the simulation, (ii) presence of SN feedback; (iii) presence of AGN feedback, (iv) total number of MBHs, (v) accretion rate of the most accreting MBH, (vi) mass of the most massive MBH, (vii) stellar mass, and (viii) SFR of the central dominant galaxy. At z = 6 in AGNcone, the most massive MBH is also the most accreting one, even if this is not always true at higher redshift.

Run	SF/SN feedback	AGN feedback	# MBH	|$\dot{M}_{\rm acc} (\!{\rm M}_\odot\,$| yr−1)	MBH(M⊙)	|$M_{*}^{\rm cD} (\!{\rm M}_\odot )$|	SFRcD(M⊙ yr−1)
noAGN	Yes	No	0	0	0	1.5 × 1011	664
AGNcone	Yes	Yes	723	57.6	4.85 × 109	6.5 × 1010	116

3 THE SATELLITE SAMPLE

We identify galaxies and their related DM haloes through the amiga halo finder code (see Knollmann & Knebe 2009), using a minimum of 20 bound particles to define a halo. The merger tree for each halo at z ≃ 6 is built by tracing back in time the constituent DM particles: their ID is matched in the progenitor structures in the previous snaphots. Baryons are assigned to their related galaxy (i) when a gas or a stellar particle is found within βr_vir of a given halo, where r_vir is the virial radius of the halo and β = 0.3 (similarly to what has been done in Rosdahl et al. 2018 and in Costa, Rosdahl & Kimm 2019) and (ii) when the velocity of the considered particles is lower than the escape velocity, determined by the DM potential well.²

In the following analysis, we select only those satellites with a contamination from low-resolution DM³ particles lower than 20 per cent in mass.

We consider only galaxies with M_vir > 10⁹ M_⊙ and with a minimum stellar mass of M_* = 10⁷ M_⊙. These thresholds are a good compromise to minimize the numerical errors, still maximizing the number of objects for statistical significance.

The selected samples at z = 10 consist of 35 galaxies for run noAGN and 36 for run AGNcone. At z ≃ 6 the difference between the number of satellites in the two runs becomes notable, being in noAGN, 30 per cent larger than in AGNcone (82 versus 56). We detail this difference in Section 3.1.

We find that at z ≃ 6, the virial radii of satellites range from about 2 to about 21 kpc and masses from 10⁹ to 10¹¹ M_⊙. The distribution of the galactic virial masses at z = 6 is shown in the left-hand panel of Fig. 1 for both the analysed runs. Even though the initial conditions are identical and the virial mass is dominated by the DM component, which is far less sensitive to the feedback prescription than baryons, small but appreciable variations are present. Feedback from MBHs seems to have a remarkable effect both in redistributing the baryonic component among the companions and in driving several satellites to coalescence, resulting in larger systems. After almost 1 Gyr from the beginning of the simulation, AGNcone galaxies are less peaked around ∼10⁹ M_⊙, than noAGN, whereas the cumulative mass of the galaxy populations (the right-hand panel of Fig. 1) shows that the total mass is conserved. This result demonstrates that no mass is lost in AGNcone through tidal disruption processes of smaller haloes.

Figure 1.

Left-hand panel: Virial mass distribution of the galaxy samples in noAGN (blue) and AGNcone (red) at z = 6. Right-hand panel: Cumulative virial mass function for the same samples. Note the highest mass bin, containing the cD galaxy of the proto-cluster.

At z ≃ 6, about half of the galaxies in the sample of AGNcone hosts at least a MBH, whose mass ranges from ∼10⁵ M_⊙, to 4.8 × 10⁹ M_⊙. There are both quiescent MBHs (or with a negligible accretion rate) and strongly accreting MBHs with   ∼60−70 M_⊙ yr⁻¹ (see Table 1). In particular, the most accreting MBH does not remain in the same object during the whole evolution of the galactic proto-cluster and it is not always the most massive one.

To summarize, the galaxies examined here exhibit a complex network of AGN, whose emitted energy varies with time and whose geometrical distribution, along with the preferred direction of emission of feedback, requires a proper model to study their influence on the satellite population.

3.1 Redshift evolution of companions

To start investigating the effect of AGN feedback on the quasar environment, we compute the redshift evolution of the number of companions, as shown in Fig. 2. The three families of lines refer to those objects enclosed within different spheres, centred on the ‘accretion-weighted centre’, c_aw (see Appendix  A), and with increasing radii: 1 central dominant (cD) galaxy virial radius (about 66 kpc at z = 6), 3 r_vir, and 5 r_vir.

Figure 2.

Number counts of satellites as a function of redshift. Solid, dotted, and dashed lines indicate the number count evaluated within 1, 3, and 5 virial radii, respectively. Blue lines mark the trend for noAGN, whereas red lines refer to AGNcone. Differences between the runs significantly increase after about z = 7, here highlighted with a vertical black line.

The number of satellites in the two simulations (red and blue lines) is very similar at any redshifts within 1 virial radius, starting from a few sources at z = 10 and reaching about 10 objects at z = 6. At larger distances from the cD, an increasing difference between the runs starts to arise: for z ≲ 7, the number of satellites in the AGNcone case is smaller than in the noAGN one. After z ∼ 7, i.e. when the outflow starts affecting the properties of the host galaxy,⁴ the number of satellites in the outer regions of AGNcone becomes smaller than the number in noAGN, whereas a substantial agreement is maintained at smaller radii. For r > r_vir, AGNcone satellite number reaches a peak at z ∼ 6.7 and is constantly reduced at lower-z. Differently, in noAGN the satellite number is always a decreasing function of redshift. The highest discrepancy of 24 satellites between the two simulations is reached at z = 6, within 5 r_vir.

3.2 Redshift evolution: satellite properties

Fig. 3 shows the redshift evolution (6 ≲ z ≲ 10; for a sub-sample of snapshots) of several satellite⁵ properties: gas mass (M_gas), star formation rate (SFR), stellar mass (M_*), and gas metallicity (Z) in solar units. The median of the gas content of our satellite population slowly decreases with time, from ∼6 × 10⁸ M_⊙ to ∼10⁸M_⊙ in about 0.5 Gyr. As a consequence, the capability to form stars of the satellites decreases: the median SFR, in fact, varies from ∼1 to 0.1 M_⊙ yr⁻¹. The median stellar mass fluctuates around M_* ∼ 3 × 10⁷ M_⊙, with a shallow increase from z = 10 to z = 6. The ISM is consequently gradually enriched in metals as stars are formed and explode as SNe, varying between ∼0.05 and 0.1 Z_⊙.

Figure 3.

Redshift evolution of satellite median properties (noAGN in blue, AGNcone in red). Clockwise from top left: gas mass, M_gas; star formation rate, SFR; total gas metallicity, Z in solar units (Z_⊙ = 0.0196, according to Vagnozzi 2019), and stellar mass, M_*. The error bars represent the 30th (error bar-lower end) and the 70th percentile (error bar-upper end) of the distribution. We note that the median stellar mass in the AGNcone run has a trend compatible with a higher cumulative SF history with respect to noAGN.

The differences between the runs are minimal and well within their error bars, thus suggesting that AGN feedback only plays a minor role in shaping the overall evolution of satellite properties. We notice however that the values of M_gas and SFR are higher in AGNcone till z = 7, whereas M_* and Z are almost always larger in AGNcone with respect to noAGN.

3.2.1 Stellar mass–metallicity relation

In Fig. 4, we compare our results with observational data, concerning the M_*−Z relation in isolated galaxies. In particular, we show the gas-phase metallicity from the oxygen versus hydrogen abundance ratio for the redshifts z = 6, 7, 8, along with a linear fit of all these satellite populations. As suggested by observations, also in our simulations, systems with increasing stellar mass are progressively more polluted in metals. The normalization of the sub-linear M_*−Z relation increases with decreasing redshift indicating that the overall metallicity floor of the galaxy group increases through cosmic times.

Figure 4.

Stellar mass–metallicity relation for noAGN (left-hand panel; blue, light blue, and cyan points for z = 6, 7, 8, respectively) and AGNcone (right-hand panel; red, light red, and pink points for z = 6, 7, 8, respectively). A linear fit of the same colour is superimposed to each distribution. The best-fitting equations for z = 6 populations are y = 0.33x + 5.6 and y = 0.25x + 6.2 for noAGN and AGNcone, respectively. For comparison, we plot over each panel some of M_*−Z relations from the gas phase in star-forming galaxies: Mannucci et al. (2010) at z = 0 (light green line with a shaded region representing 90 per cent of the SDSS galaxies), Cullen et al. (2014) at z ≳ 2 (triangles), Maiolino et al. (2008) at z ∼ 3.5 (squares), Faisst et al. (2016) at z ∼ 5 (stars), and Harikane et al. (2020) at z ∼ 6 (crosses).

Similarly to the other intrinsic properties, the general trends of the run noAGN are almost indistinguishable from the AGNcone results. If we linearly interpolate the M_*−Z distributions at z = 6, we notice that the slope of the best-fitting relation is only slightly steeper in noAGN than in AGNcone, being 0.33 ± 0.03 versus 0.25 ± 0.05, thus consistent with the errors. This indicates that the AGN feedback does not change significantly the process of gas pollution in the galaxy group as a whole. However, we also note that the difference between the zero-points is marginally larger, with 5.6 ± 0.2 versus 6.2 ± 0.3, which is compatible with a higher SF activity in the AGNcone past evolution.

The comparison between observations and our results is however not straightforward. Due to sensitivity limitations, observations are still unable to probe the low-mass end (M_* ≲ 10⁹ M_⊙) of the relation, that is instead statistically covered by our simulated data. Future deeper observational surveys are therefore necessary to reduce this bias. In general, we note that the level of metal enrichment reached by the galaxy groups of both runs seems to agree with the extrapolation from the data of z = 6 galaxies (Harikane et al. 2020) (green crosses).

3.3 Spatial distribution of companions

In this section, we analyse how satellite properties are spatially distributed within the galaxy group, to check any possible correlation with the AGN activity, occurring in the AGNcone simulation.

The kinetic feedback, as modelled in our simulations, might remove gas from galaxies close to c_aw and transfer it to more peripheral systems. The way it finally affects the gas distribution of the galaxy group is anyway complex. In the DM distribution, high-density peaks tend to be clustered (Bardeen et al. 1986), implying that most massive haloes are closer to the centre of mass of the galaxy proto-cluster (‘mass segregation’). Since the centre of mass almost coincides with c_aw at z = 6, one should expect the effect of quasar feedback to be higher in the closest, more massive satellites. However, as a consequence of their deeper potential wells, massive systems might retain their gas content more easily with respect to less massive ones, thus being more resilient to the possible passage of outflows launched from the galaxy itself or coming from close companions.

To study the effect of quasar feedback on the spatial distribution of satellites, in Fig. 5 we report the satellite intrinsic properties (M_gas, SFR, M_*, Z) as a function of their distance from c_aw. Both runs show a highly dense environment, where the most gas- and metal-rich, star-forming, galaxies are preferentially located at small distances from the center of the galaxy group, independent of the presence of quasar feedback. This suggests that the satellite distribution within the galaxy proto-cluster is dominated by mass segregation and quasar feedback plays, at most, a second-order effect. In general, we find the distribution of satellites to be quite flat at distances larger than 100 kpc.

Figure 5.

Satellites intrinsic properties as a function of the distance from the accretion-weighted centre c_aw in the runs noAGN (left column) and AGNcone (right column) at z = 6. The main, central galaxy is not shown. From top to bottom: Gas mass, stellar mass, SFR, and metallicity. Brown stars refer to the most active AGN, here defined as those galaxies hosting at least an MBH with a total accretion rate |${\dot{M}_{\rm acc}\gt 1}$| M_⊙ yr⁻¹. The median of the distributions are superimposed in each panel as filled squares. Uncertainties are quantified as the 30th and 70th percentiles of each distribution. The values of the median points and the lower bounds of the error bars not shown in the plots are equal to zero.

Still, in AGNcone there is a larger number of star-forming (SFR ≳ 1 M_⊙ yr⁻¹) massive (M_* ≳ 10⁹ M_⊙) galaxies with respect to the noAGN run. This result, along with the one reported in Section 3.2 (larger values of M_gas, M_*, SFR, and Z in AGNcone with respect to noAGN), suggests to refine the comparison between the runs by focussing on individual satellites. In fact, the effect of quasar feedback on the environment could be washed-out by averaging the satellites properties over the entire population. Different galaxies may indeed perceive feedback effects at different times, depending on their relative position with respect to the most accreting MBH and on the variability of the MBH itself.

4 INDIVIDUAL COMPANIONS

The goal of this section is to compare the SFR and M_* redshift evolution of a satellite in the AGNcone run with the corresponding satellite in the noAGN run. Finding the same satellite, available in both runs in the same redshift interval is not trivial, since different feedback prescriptions result into different merger rates (see Section 3.1): in AGNcone, galaxies merge faster and more easily than in noAGN, possibly because of a more diffuse and massive stellar component around the galaxies (see Fig. 6). We thus describe in details the procedure we follow to set up our sample of galaxy couples.

Figure 6.

z = 6.3 stellar surface density maps of the central 270 kpc of noAGN (left-hand panel) and AGNcone (right-hand panel). Matched galaxies are highlighted in the panels with a circle of the same colour.

We select a sample of galaxies that are characterized by (i) the same position (within 5 comoving kpc) and (ii) the same virial mass (within 10 per cent) in both runs. Among the possible candidates from this first selection (iii) we select galaxies that, at z = 6, have different distances from c_aw and different masses. This condition allows us to probe different regions of the quasar environment in terms of mass segregation. Furthermore, we choose galaxies whose (iv) merger tree starts at least⁶ at z = 9 and reaches z = 6. Finally, (v) we exclude from our sample those galaxies that have hosted a powerful AGN for a significant amount of their evolutionary history.

Although only a small fraction of satellites hosts a powerful AGN, this last conditions is essential to isolate the effect of external from internal AGN feedback on satellites. AGN-driven outflows may, in fact, subtract part of the cold gas component from the interstellar medium of the host galaxies, affecting their SFRs and stellar masses. To select galaxies hosting powerful AGN, we first walk backwards the merger tree of each galaxy (selected at redshift z) and then we sum, at each redshift, up to the formation redshift⁷z_form, the accretion rate |${\dot{M}_{acc}}$| of all the MBHs located within βr_vir. After this, we compute the mean |$\langle \dot{M}\rangle _{\rm acc}$| of the cumulative accretion rate of the selected galaxy over the whole redshift range and exclude from the analysis those galaxies for which |${\langle \dot{M}\rangle _{\rm acc}\gt 0.1}$| M_⊙ yr⁻¹. Via this method, we always exclude at least the cD galaxy, hosting some of the most active and massive MBHs during the whole evolution.⁸

The final selection, composed by six galaxies with virial masses⁹ ∼10¹⁰–10¹¹ M_⊙, is highlighted in Fig. 6; in Fig. 7, we show the redshift evolution of stellar mass, SFR and distance from c_aw for each galaxy of the sample. At higher redshift, the values of M_* and SFR (first and third rows) are almost identical in the two simulations, for all the galaxies. After z ∼ 8, the trends start to diverge, always resulting in higher stellar mass and SFR in AGNcone satellites (the three galaxies on the left – i.e. A, B, and C – have larger differences with respect to the three rightmost galaxies D, E, and F, as it will be discussed in Section 4.1). The sudden fluctuations observable in the evolution of some satellites (e.g. object A at z ∼ 7) are due to minor mergers or close fly-bys which temporary increase the mass within βr_vir. The abrupt decrease of the red line of C is due to the upcoming coalescence of the object with the cD, at the very last snapshot. Hence, C system is majorly stripped of both its gaseous and stellar components. The increasing differences are more easily noticeable in the relative difference panels (second and forth rows), where positive values represent higher quantities in AGNcone run. As it is clear from the fifth row, satellite distances evolve quite differently among the selected objects, further increasing the generality of the sample: while B, C, and F approach the centre, galaxy D orbits increasingly far from the cD, and galaxies A and E do not significantly change their relative position during the entire simulation time.

Figure 7.

Redshift evolution of the six galaxies in the sample, as described in Section 4. From top to bottom: Stellar mass, stellar mass relative difference (i.e. |$[M_{*, \rm AGNcone}-M_{*, \rm noAGN}]/M_{*, \rm noAGN}$|⁠), SFR, SFR relative difference (i.e. |$[\rm {SFR}_{\rm AGNcone}-\rm {SFR}_{\rm noAGN}]/\rm {SFR}_{\rm noAGN}$|⁠), distance from c_aw, and |$\mathcal {E}_{\rm AGN}^{\rm cum}$|⁠. The same colour coding of Fig. 3 is applied for the two runs. Horizontal lines mark the zero-level for the relative differences (second and forth rows), the time-averaged position of the cD virial radius (fifth row), and the threshold 10⁵¹ erg for the cumulative integrated flux (solid lines in sixth row; see Section 4.1).

To visualise the relative position of the selected satellites with respect to the geometry of the bi-conical outflows, we show in Fig. 8 a Mollweide map, where the location of each galaxy at z = 6.3 is superimposed to the radial velocity map of the gas within a sphere of radius 10 kpc and centred on the most accreting MBH in that snapshot.

Figure 8.

Mollweide view of the gas radial velocity within 10 kpc of the most accreting MBH in AGNcone, at z = 6.3. The map shows the effect of the bi-conical feedback, resulting in two visible gaseous outflows. In the AGNcone run, the orientation of the outflow is randomly assigned to each MBH at the seeding time, and kept fixed for the whole MBH lifespan. The gas velocity distribution, however, also depends on the activity of the other MBHs present in the simulation. The direction of the emission axis for the most accreting MBHs is recovered by fitting the radial velocity map of the gas. The black stars mark the position of the opposite outflow axes, as they are derived from the velocity distribution, whereas the black lines represent the boundaries of the outflows, ideally confined within 45° from the axis. Black dots show the angle position of the satellites at z = 6.3. Larger coloured dots mark the galaxies selected for the matching study of Section 4. The colour coding is the same used in Fig. 6.

4.1 AGN feedback on individual companions

Given the non-isotropic emission of the out-flowing gas in AGNcone, the evaluation of the feedback effect cannot rely only on the distance from the source and the target galaxy. In addition, galaxies randomly experience the influence of various AGN during various and discontinuous accretion periods. It is then necessary to define a function of both distance and emitted power that can also take into account the relative position of a target galaxy with respect to the emission axes.

Let us consider a galaxy (top right in Fig. 9) at a distance d from an accreting MBH (bottom left), subtending an angle Ω_gal. The MBH launches an outflow (dashed lines) towards the galaxy in an opening angle Ω_f (identified by the dotted lines). We quantify the energy received by the target galaxy at a given redshift z, from any ith MBH which has been accreting with |$\dot{M}_{{\rm acc}, i}$| in the time interval Δt_snap, and whose outflow cone encompasses the galaxy centre.¹⁰

Figure 9.

Schematic representation of the model applied to AGNcone. From an accreting MBH, a gaseous outflow with mass M_w is ejected with velocity v_w in the cone with aperture Ω_f and intercepts a satellite at the distance d, subtending an angle Ω_gal. M_e is the gas envelope mass encountered by the outflow when it is ejected.

The total gas mass involved in the feedback emission process around the MBH is M_tot = M_w + M_e, where M_w is the mass of the ejected wind and M_e is the mass of the gas in the environment surrounding the MBH and entrained within the outflow. Depending on the geometry of the feedback prescription, M_e can be written as

$$\begin{eqnarray*} M_{e} = \frac{\Omega _{\rm f}}{4\pi }\frac{4\pi }{3}\rho d^3, \end{eqnarray*}$$

(5)

where ρ is the average density around the MBH. If v_f is the velocity of the gas when it hits the target galaxy, then momentum conservation¹¹ allows us to write

$$\begin{eqnarray*} M_{\rm w}v_{\rm w}+M_{\rm e}v_{\rm e} = v_{\rm f}M_{\rm tot}, \end{eqnarray*}$$

(6)

where we assume v_e = 0 when the outflow is launched.

Neglecting internal energy,¹² the energy deposited by the ith MBH on the target galaxy (⁠|$\mathcal {E}_{\rm AGN, i}$|⁠) is related to the final kinetic energy E_f through the relation:

$$\begin{eqnarray*} \mathcal {E}_{\rm AGN, i} \equiv E_{\rm f} \frac{\Omega _{\rm gal}}{\Omega _{\rm f}} = \frac{1}{2}M_{\rm tot}v_{\rm f}^{2} \frac{\Omega _{\rm gal}}{\Omega _{\rm f}}, \end{eqnarray*}$$

(7)

where the factor Ω_gal/Ω_f accounts for the fact that only a fraction of the total mass ejected by the MBH actually intercepts the target galaxy.

Isolating v_f from equation (6) and considering that the envelope mass in equation (5) fully dominates over the wind mass (i.e. M_tot ≃ M_e), equation (7) becomes

$$\begin{eqnarray*} \mathcal {E}_{\rm AGN, i} = \frac{3}{2}v_{\rm w}^2 \frac{M_{\rm w}^2}{\rho d^3} \frac{\Omega _{\rm gal}}{\Omega _{\rm f}^2}. \end{eqnarray*}$$

(8)

In equation (8), the solid angle Ω_gal – subtended by the target galaxy with respect to the MBH position – can be approximated as

$$\begin{eqnarray*} \Omega _{\rm gal} = \frac{A_{\rm gal}}{d^2}, \end{eqnarray*}$$

(9)

where A_gal is the projected area of the galaxy, as it is seen from the AGN.¹³ In our analysis, we use, with β = 0.3, in agreement with our choice to compute the properties of the haloes (see Section 3). At the same time, Ω_f = 2 × 2π(1 − cos α), with α = π/4 in AGNcone, for each of the two cones where the energy is distributed.

In order to estimate M_w, we multiply the outflow rate |$\dot{M}_{\rm w}$| into the snapshot time interval Δt_snap:

$$\begin{eqnarray*} M_{\rm w} = \dot{M}_{\rm w} \Delta t_{\rm snap}, \end{eqnarray*}$$

(10)

where |$\dot{M}_{\rm w}$| is obtained from the accretion rate |${\dot{M}_{acc}}$| via equations (3) and (4). Accordingly, the energy-conservation equation can be written as

$$\begin{eqnarray*} \dot{E}_{f} = \frac{1}{2}\dot{M}_{\rm w}v_{\rm w}^2 = \epsilon _{\rm f}\epsilon _{\rm r}\dot{M}_{\rm acc}c^2. \end{eqnarray*}$$

(11)

Finally, we can write down the form of the deposited energy, by substituting M_w in equation (8):

$$\begin{eqnarray*} \mathcal {E}_{\rm AGN, i} = 6 \pi \epsilon _{\rm r}^2 \epsilon _{\rm f}^2 c^4 \beta ^2 \frac{(\Delta t_{\rm snap})^2}{v_{\rm w}^2} \frac{\dot{M}_{\rm acc}^2 r_{\rm vir}^2}{\Omega _{\rm f}^2 \rho d^5}. \end{eqnarray*}$$

(12)

Operatively, we evaluate the gas density ρ by summing up all the gas particles in the cone subtended by the solid angle Ω_gal, with respect to the ith MBH and the accretion rate |$\dot{M}_{\rm acc}$| at the time t − Δt_snap, in agreement with what is done in Appendix  A.

The derived equation (12) has a quite intuitive dependence on the MBH accretion rate, |$\dot{M}_{\rm acc}^2$| (the stronger the AGN, the higher the effect) and on the size of the galaxy, proportional to |$r_{\rm vir}^2$| (the larger is the galaxy projected area with respect to the AGN, the higher is the energy harvested by the system). Also, the dependencies on the inverse of the distance galaxy-MBH, d⁵, and on the density of the circumgalactic medium (CGM) are expected. As a matter of fact, the density of the environment and the length of the path that the outflow has to travel both contribute to lower the final energy. In a purely momentum-conservation scenario, the kinetic energy of the ejected gas is continuously distributed on the surrounding medium and, in very dense environment, can be completely lost in the CGM before reaching the target system. A final consideration concerns the solid angle Ω_f, which favours more collimated outflows. This factor, on one hand is able to reduce the volume where the feedback energy is diluted but, on the other, it deeply affect the probability that a galaxy is affected at all by the outflow.

In order to quantify the total effect of the MBH population on each galaxy, we define the total |$\mathcal {E}_{\rm AGN}$| as

$$\begin{eqnarray*} \mathcal {E}_{\rm AGN} = \sum _{i=1}^8 \mathcal {E}_{{\rm AGN}, i}, \end{eqnarray*}$$

(13)

where the summation is carried over the eight most accreting MBHs as discussed in Appendix  A. In our model, for a given target galaxy, we consider the contribution of the ith accreting MBH only if the galaxy centre of mass is enclosed within the boundaries of the outflow launched by the AGN (see the black solid lines in Fig. 8). This equation provides a rough estimate of the energy deposited in a time interval Δt_snap by the most-accreting AGN on a target galaxy, at redshift z.

At any generic redshift z, the cumulative energy deposited on a target galaxy (⁠|$\mathcal {E}_{\rm AGN}^{\rm cum}$|⁠) is given by the integral of |$\mathcal {E}_{\rm AGN}$| between the galaxy formation redshift z_form, and z:

$$\begin{eqnarray*} \mathcal {E}_{\rm AGN}^{\rm cum} \equiv \sum _{\quad z=z_{\rm form}}^{z} \mathcal {E}_{\rm AGN}. \end{eqnarray*}$$

(14)

The results of this model are shown in the bottom line of Fig. 7, where the horizontal line, at |$\mathcal {E}_{\rm AGN}^{\rm cum}=10^{51}$| erg, refers to the energy injected in the medium by a single SN event and has the only purpose of helping the reader to distinguish the objects in relation to their received energy. According to the values of |$\mathcal {E}_{\rm AGN}^{\rm cum}$|⁠, we separate the galaxies more affected by the AGN feedback (A, B, C, on the left) from the less affected ones (D, E, F, on the right). The SFR of galaxies from the first group is more enhanced (up to a factor of 3) with respect to the second group (less than a factor of 2).

We note that in A, B, and C, |$\mathcal {E}_{\rm AGN}^{\rm cum} \gtrsim 10^{51}$| from z ∼ 8 up to the end of the simulation, while in D, E, and F |$\mathcal {E}_{\rm AGN}^{\rm cum} \gtrsim 10^{51}$| only for 7.5 < z < 6. In other words, it seems that galaxies receiving the energy by the AGN feedback for a longer time are more strongly affected, in terms of SFR and M_*. We discuss in the next section a possible interpretation of these findings.

5 INTERPRETATION

In this section, we first summarize the main findings obtained from the comparison between the runs noAGN and AGNcone, and we propose our interpretation for our results.

In both runs, we have identified a sample of galaxies in the redshift range 6 < z < 10 and characterized their properties (distance from the center of the galaxy groups, SFR, M_*, M_gas, Z). From the comparison between the samples extracted from the two simulations, we can conclude the following:

• In the AGNcone run, satellites are less numerous, especially in the outer regions, and they are hosted by more massive DM haloes (see Section 3).

• Although the differences between the median properties of the two samples at all redshifts are not noteworthy, individual AGNcone companions are more star forming and massive; the SFR enhancement in those satellites which are influenced the most by the surrounding AGN, reaches a factor up to 4 (see Sections 3 and 4).

We suggest that poinAstron. Soc. Pac. Cothe powerful quasars located at the centre. The disperse gas and stars can boost the effect of dynamical friction when two galaxies approach, thus lowering the dynamical time-scale for a merging to occur. DM–haloes follow a similar trend, suggesting that differences in the baryonic component are transferred to DM structures that, in the AGNcone run, become more massive, and therefore less in number. Interestingly, the effect is only relevant far from the cD and almost absent near the centre, where the dynamics is completely dominated by the cD environment. Strong feedback could intrinsically favour coalescence, decrease the number of systems, and speed up their bottom-up growth. These phenomena play a fundamental role in the number count of satellites around z = 6 quasars, as discussed in the next section.

Concerning point 5, we note that AGN feedback do not significantly affect the evolution of the satellite population as a whole. The median of the satellite properties show only minor differences, which are though compatible with a scenario where, in noAGN, a higher M_gas results in a higher SFR, till the exhaustion of the gas reservoirs. We note that, the cumulative effect on the SFR and M_* is both due to the higher merger rate (more likely to occur far from the cD) and to the direct effect of the outflow on the satellite (stronger at smaller distances from the cD). We suggest that the enhanced SF in those satellites directly invested by the AGN outflows can be explained by (i) an increase of the fuel available for the SF process, (ii) a boost in its efficiency due to the induced shocks, (iii) the formation of stars within the outflow itself (Maiolino et al. 2017; Gallagher et al. 2019).

Several evidences of AGN positive feedback on to galaxy satellites have been reported in the literature. Gilli et al. (2019) analyse the deep multiband field of a type II radio galaxy at z = 1.7 and discover an overdensity of galaxies. These authors suggest that the star formation in satellites is promoted by the compression of their cold interstellar medium around the AGN-inflated bubbles. The possibility of a jet-induced star formation is also consistent with the SF measurement of a galaxy near the local radio-galaxy NGC–541 (Croft et al. 2006), and supported by dedicated numerical simulations (Fragile et al. 2017).

Recently, Martín-Navarro et al. (2021) measured M_* and SFR in a large sample of satellites in SDSS galaxies finding that their SFR is modulated by their relative position with respect to the central galaxies, being higher in those objects presumably affected by the outflows from the central galaxies. These authors concluded that AGN-driven outflows can influence the environment well beyond the host galaxy (1−2 R_vir), preserving (or even enhancing) the SF of those satellites along the direction of the injected energy.

There are some unavoidable limits associated with our model to quantify |$\mathcal {E}_{\rm AGN}$|⁠. On the one hand, being based only on momentum conservation, it would underestimate the final energy deposited into the target galaxy by a fully energy-conserving outflow. On the other hand, the outflow launched by AGN may interact with intervening ISM and CGM material, potentially losing a significant amount of energy, or even being stopped in the most extreme cases. This would instead lead us to overestimate the final impact on the satellite.

In addition, the limited spatial and temporal resolution of the simulation suite, the random explosion of SNe, and the possible accretion episodes of MBHs within the satellites themselves,¹⁴ may dilute the effect of feedback from external AGN. All these caveats prevent us from disentangling entirely the effects of the several physical processes affecting the CGM and enhancing the SF in satellite galaxies.

6 COMPARISON WITH OBSERVATIONS

In order to evaluate how AGN feedback affects our capability to detect an overdensity through the number counts of galaxies, we compute both the [C ii] and UV luminosity of all the simulated satellites as detailed in Appendix  B. Fig. 10 shows the satellite luminosity as a function of the distance from the centre of the system at z = 6. Horizontal dashed lines mark the observed limit luminosities L_lim of current observational campaigns at high redshift: |$L_{\rm lim, [C\, \small {II}]} = 10^{8} {\rm L}_{\odot }$| (Venemans et al. 2020, hereafter V20) and L_{lim, UV} = 10¹⁰L_⊙ (Marshall et al. 2020). In agreement with our previous results, AGNcone produces more luminous, and thus more easily detectable, galaxies than the control simulation. In details, above |$L_{\rm lim, [C\, \small {II}]}$| (L_{lim, UV}), there are 1 and 3 (1 and 5) satellites in noAGN and AGNcone, respectively.¹⁵ For a comparison with UV data, we refer the reader to Di Mascia et al. (2021a), where radiative transfer calculations are fully accounted for. Here, we focus on the comparison between our predictions and currently available ALMA data of z ∼ 6 quasars.

Figure 10.

[C ii] emission (top row) and UV emission (bottom row) as a function of distance from c_aw, for the same z = 6 sample of Fig. 5: noAGN on the left and AGNcone on the right. Horizontal dashed lines provide an estimate of the current instrumental sensitivities, derived from some of the deepest observations in the various bands. As in Fig. 5, where the brown stars mark the hosts of the most active AGN in AGNcone, the filled squares show the medians of the distributions and the error bars refer to the 30th and 70th percentiles. The values of the median points and the lower bounds of the error bars not shown in the plots are equal to zero.

In particular, we consider the results of a recent high-resolution ALMA survey of 27 (previously [C ii]-detected) quasars at z ∼ 6 by V20 see also (see also Decarli et al. 2017, 2018). The authors detected 17 companions¹⁶ with |$L_{\rm [C\, \small {II}]}\gtrsim 10^{8} {\rm L}_{\odot }$|⁠, corresponding to an average of 0.6 companions per field. We further notice that some of the V20 observed quasars present multiple (2–3) companion galaxies.

In Fig. 11, we show the number of detectable satellites in our simulations as a function of L_lim. We select all those satellites above the luminosity threshold L_lim and within a spherical volume of about 680³ kpc³, equivalent to the average volume¹⁷ observed by V20. Three subsets of the V20 companion list are also shown: the mean satellite number of the whole 27 quasar sample (green circles), the mean satellite number of the four most FIR luminous quasars (L_FIR ≳ 10¹³ L_⊙; green squares), and the satellites of the most FIR-luminous quasar (J0305-3150, with L_FIR = 1.2 × 10¹³ L_⊙; green crosses). Fig. 11 shows that the number of detected satellites in V20 increases with the FIR luminosity of the central quasar. A possible explanation for this trend is that more FIR-luminous quasars are likely hosted by galaxies with higher SFR, and therefore in more massive haloes. As a result, they tend to be more biased.

Figure 11.

Number of satellites in AGNcone (red) and noAGN (blue) with |$L_{\rm [C\, \small {II}]}\gt L_{\rm lim}$|⁠, with respect to L_lim, in a volume of about 680³ kpc³. Shaded areas show the related Poissonian errors (with a 68 per cent confidence level; Gehrels 1986). From left to right, the vertical dotted lines mark the [C ii] sensitivity of an ALMA observing programme of 10 and 1 h on source, respectively. Green symbols show the latest observational data from V20: filled circles refer to the whole sample of 17 satellites around 27 quasars, filled squares refer only to the quasar population with L_FIR ≳ 10¹³L_⊙, and Xs mark the densities around the most luminous quasar in FIR, i.e. J0305-3150.

From a comparison between our results and the most FIR-luminous quasar in the ALMA sample, we note that AGNcone (where the cD has L_FIR ∼ 5 × 10¹³L_⊙; Di Mascia et al. 2021a) agrees quite well with ALMA data, whereas noAGN predicts a number of satellites lower than observed. This implies that, for the most luminous source in FIR, the positive feedback of quasar outflows on galaxy companions is required to reproduce observational data, even though more simulations are required to verify this effect with better statistics.

In contrast, although the average observed number count is well within the Poissonian noise (shaded areas), it is systematically lower than our predictions from both our runs. This can be explained in two ways: (i) not all the quasars, but only the most FIR-luminous ones, live in overdense regions; (ii) the number of satellites detected by V20 only provides a lower limit to the actual number. The latter hypothesis can be related to observational artefacts.

First, of all, the finite angular resolution of real data can prevent us from correctly quantifying the actual number of satellites around the quasar. In fact, [C ii] emission in the bright, central regions of z ∼ 6 quasars has typical sizes of about 1−5 kpc and may therefore arise from multiple, unresolved sources. While counting satellites from simulations, we are instead assuming that all of them, even the closest to the central source (within ∼1 kpc), are spatially resolved. The inclusion of this observational artefact in our simulations would therefore automatically improve the agreement with observations.

We further note that, while ALMA observations probe distances up to 2−3 Mpc from the central quasar, the high-resolution volume of our simulations is only limited to the inner 200−300 kpc. If we consider ALMA observations only in a region from the centre within ±200 km s⁻¹ (corresponding to about 280 kpc at z = 6), the average number of satellites is 0.3 (i.e. 9 galaxies around 27 quasars), resulting¹⁸ into a mean density of about 4 × 10⁻⁸ satellites per kpc³. In the noAGN and AGNcone simulations, we retrieve mean densities of 2 × 10⁻⁸ satellites per kpc³ and 6 × 10⁻⁸, respectively. This suggests that ALMA observations in this restricted volume are in remarkable agreement with the mean density expected from of our simulations.

Finally, the large asymmetry of the volume probed by ALMA (≲100 kpc on the sky plane versus ≳2−3 Mpc along the line of sight) could introduce a further bias. In some quasars observed by V20, the number of satellites is larger than the average (e.g. J0305-3150 in Fig. 11). We can explain these cases as ‘fortunate sources’, where the distribution of satellites is somewhat aligned with the ALMA volume. In these sources, ALMA observations are able to detect the real and total amount of companions, which is also perfectly compatible with the number count predicted in AGNcone, where the density is computed without assuming a preferred line of sight. A more quantitative-statistical approach on this hypothesis will be the focus of a future work.

6.1 Predictions for ALMA

Most of the ALMA data available so far to study z ∼ 6 quasars are limited to shallow observations (≲1 h per source). It is therefore interesting to investigate what we can learn from deeper observations. Fig. 11 shows that, if quasars are hosted in overdense environments (as predicted by our simulations, by construction), we would be able to detect up to 6 satellites (10 as an upper limit) in a luminous quasar neighbour with 10 h of ALMA observing time.

In the case of much deeper observing programmes (L_lim < 10^6.5 L_⊙), the number of observable satellites in AGNcone is smaller than in noAGN, contrarily to what occurs in the case of shallower observations. This inverted trend is possibly due to the twofold effect of quasar feedback on the surrounding satellites: enhancing the SFR (and therefore the luminosity) of luminous satellites and lowering their intrinsic number.

Although such deep observations are beyond the capabilities of current (sub-)millimetre observatories, the discussed trend still implies that quasar feedback may leave signatures on the slope of the satellites number count, thus suggesting further investigations on this topic. However, we remark that the limited resolution of our simulation could play a role on this effect, since low-mass galaxies could be differently impacted by quasar feedback in simulations with higher resolution, leading to different results.

6.2 Predictions for JWST

Finally, we focus on JWST in order to understand how this mission will improve our knowledge of high-z quasar properties and their environment. Fig. 12 shows the expected number of satellites emitting in UV as a function of their luminosity L_UV > L_{lim, UV}. The UV emission of satellites is corrected for dust attenuation with a dust-to-metal ratio f_d = 0.08, according to the results of Di Mascia et al. (2021b, see Appendix  B for further details), and the shaded area takes into account both the minimum and maximum optical depth resulting from radiative transfer calculations and the Poissonian uncertainty.

Figure 12.

Analogously to Fig. 11, we show the number of observable satellite in the UV band with respect to the luminosity threshold L_{lim, UV}. Solid lines refer to the same UV luminosities, corrected for dust extinction of Fig. 10. The uncertainties, shown by the shaded areas, are estimated by summing the contribution of the Poissonian noise (evaluated as in Fig. 11) and the variation range from the minimum to the maximum optical depth considered. From left to right: The dotted vertical lines show the sensitivity thresholds that JWST can reach with 10 h and 1 h of observing time, respectively.

According to AGNcone, we would be able to detect between 3 and 8 satellites with 1 h of observing time, and between 4 and 10 satellites via a 10 h observing programme. We also report our predictions for the noAGN case. Although we demonstrated that the absence of AGN feedback would decrease our capability of detecting quasar companions, this represents a lower limit (1−2 satellites with 1 h of observing time and 2−6 objects with 10 h) for a deep JWST observation programme on a single quasar.

These calculations, along with the ones presented in Fig. 11, demonstrate that both ALMA and JWST will be essential to improve our understanding of high-z quasars environment, by increasing the number of observed satellites, probing different scales and emission processes. In particular, the synergy between these two observatories will be of outmost importance, since each instrument will compensate the limitation of the other one: e.g. on the one hand, the larger field of view of JWST with respect to ALMA (∼10 arcmin² versus 1 arcmin²) will allow us to probe broader regions around quasars, on the other hand, ALMA is able to reveal even those satellites that are dust-obscured and may elude JWST observations.

7 SUMMARY AND CONCLUSIONS

We have investigated the effects of quasar outflows (i.e. feedback) on the visibility of companion galaxies by comparing two cosmological zoom-in simulations of a z ∼ 6 quasar, in which AGN feedback is either included or turned-off. We have identified satellites in both runs and determined their key physical properties such as gas and stellar mass, SFR, and metallicity. Our findings can be summarized as follows:

• Within the virial radius of the central galaxy, the number of satellites in the two runs increases with time in a very similar way, resulting in 10 companions at z ≈ 6. At larger distances, the number of satellites in AGNcone is 30 per cent smaller than in noAGN; we suggest that this effect is due to AGN-driven outflows that, by dispersing stars and gas in the surrounding region, boost the dynamical friction and, consequently, increases the galaxy merger rate. The effect is negligible at smaller distances because the orbital dynamics and the CGM density are fully dominated by the presence of the central galaxy.

• The redshift evolution of the satellite median properties does not show striking differences between the runs.

  M_gas decreases from 10⁹ M_⊙ at z ∼ 10 to 10⁸ M_⊙ at z ∼ 6; the SFR decreases from 1 to 0.1 M_⊙ yr⁻¹; M_* does not evolve appreciably from ∼3 × 10⁷ M_⊙; Z increases from 0.05 to 0.2 |$\rm Z_{\odot }$|⁠.

  To get a deeper insight into the effect of AGN feedback on its environment, we have thus considered a sub-sample of satellites, located at different distances and positions with respect to the center of the group and followed the redshift evolution of their intrinsic properties (M_* and SFR) in both runs.

  For all these satellites, we found that both M_* and SFR grow faster in AGNcone with respect to noAGN.

  We argue that the SFR enhancement in those satellites engulfed by the quasar outflow can be due to an increase of the fuel available for the SF process and/or a boost in the SF efficiency due to the induced shocks. Both possibilities are linked to the kind of feedback and SF recipes implemented in the simulation: gas particles are moved away from the galaxies hosting the most accreting MBHs towards the surrounding satellites.

• We have developed a semi-analytical model based on momentum conservation to quantify the total energy deposited on a target galaxy by the surrounding active AGN, taking into account their duty cycles and accretion rates, the outflow orientation, the circumgalactic medium (CGM) density, and their relative distance from the target galaxy. We found positive feedback in numerous satellites, proportionally to the total energy received from the accreting MBHs in their whole evolutionary history.

• We have computed the [C ii]158 μm emission of satellites and studied the effect of quasar feedback on their observed number count. When compared with the most FIR-luminous quasar of the ALMA sample, the noAGN run predicts a number of satellites lower than observed, while the AGNcone run agrees quite well with ALMA data. This implies that, for the most FIR-luminous source, the positive feedback of quasar outflows on the SFR of galaxy companions is instrumental in reproducing observational data. However, the average number of satellites observed in the whole quasar sample is lower than our predictions for both the runs. This can be explained in two ways: (i) not all the quasars, but only the most FIR-luminous ones live in overdense regions; (ii) the number of satellites reported by Venemans et al. (2020) only provides a lower limit to the actual number, as a consequence of observational artefacts (e.g. finite angular resolution and anisotropic volume probed by ALMA).

• We predict that JWST will double the current ALMA detections after just 1 h of observing time and will allow us to directly observe a major part of the quasar neighbour, avoiding any bias introduced by the small ALMA field of view. Still, ALMA will be necessary to reveal those satellites around z ∼ 6 quasars that are dust-obscured and may elude JWST observations. We find that with 10 h of observing time, ALMA and JWST will detect up to 10 satellites per sources.

In conclusion, we did not find any evidence of AGN-driven quenching on the star formation of satellites surrounding high-z quasars. Thus, we rule out a scenario in which the detection of companions can be undermined by the external feedback. AGN feedback could be even instrumental to explain the high satellite number count observed around the most FIR luminous quasars and – after correcting data for ALMA observational biases – in the whole population observed so far.

ACKNOWLEDGEMENTS

The authors thank Bram Venemans and Roberto Decarli for helpful insights on ALMA [C ii] data. SG acknowledges support from the ASI-INAF n. 2018-31-HH.0 grant and PRIN-MIUR 2017. AF and SC acknowledge support from the ERC Advanced Grant INTERSTELLAR H2020/740120. Partial support from the Carl Friedrich von Siemens-Forschungspreis der Alexander von Humboldt-Stiftung Research Award (AF) is kindly acknowledged. We gratefully acknowledge computational resources of the Center for High Performance Computing (CHPC) at SNS. AL acknowledges funding from MIUR under the grant PRIN 2017-MB8AEZ. PB acknowledges support from the Brazilian Agency FAPESP (grants 2016/01355-5, 2016/22183-8). The authors greatly thank the anonymous referee for useful comments which improved the quality of this manuscript.

DATA AVAILABILITY STATEMENT

The data underlying this article will be shared on reasonable request to the corresponding author.

Footnotes

1

In the code, gas particles are swallowed according to a stochastic method described in Springel et al. (2005a).

2

The gravitational potential here is evaluated through the analytical integration of a Navarro–Frank–White profile (Navarro, Frenk & White 1996).

3

We cannot exclude the possibility of contamination-induced effects on the galactic satellites that orbits farther out in the refined zone of the simulation volume with respect to the main galaxy. For this reason we impose such constraint. The number of excluded galaxies is never larger than 1−2 per snapshot.

4

See fig. 8 in Barai et al. (2018): |$\dot{M}_{\rm out}$| increases substantially after z ∼ 8 near the cD virial radius, and by z ∼ 7, the outflows reach the outskirt of the cluster.

5

In each snapshot, we exclude from the computation the most massive galaxy in the two runs. Even if the proper way to identify a quasar companion would require to identify those galaxies which orbit around an accreting MBH host (see the method adopted in Section 4.1) in AGNcone and to look for their counterpart in noAGN, here our approximation is equivalent for the analysis.

6

This requirement takes out from the selection those objects which suffer from violent gravitational stripping processes or merge with a galaxy of equal/higher mass, long before its counterparts in the other run.

7

We define the formation redshift z_form of a galaxy, as the first redshift (the highest) where the first ancestor of the galaxy in the merger tree is identified.

8

We note that, in principle, this criterion might miss MBHs that are highly accreting at low mass, possibly neglecting relevant effects in the case the MBH is hosted in a small mass galaxy.

9

With an average contamination from low-resolution DM particles of a few per cent for 1 out of the 6 galaxies and null for the remaining 5 objects.

10

In the following derivation, we neglect the subscript i on the MBH properties to lighten the formalism.

11

We assume that the energy of the outflow can be radiated away during the gas migration, because of gas heating and shock fronts.

12

The energy of the outflow is almost completely kinetic, also according to the feedback receipt.

13

We note that, the correct formula for a curved surface subtending the solid angle Ω_gal is |$\Omega _{\rm gal}=2\pi (1-\frac{d}{\sqrt{r_{\rm gal}^2+d^2}})$|⁠. Our approximation fails only when d ∼ r_gal, which occurs a negligible amount of times, leading our method to overestimate the solid angle, in those few cases, at most by a factor of 1.7.

14

We remark though that we limit the impact of AGN hosted in satellites by excluding those objects with |${\langle \dot{M}\rangle _{\rm acc}\gt 0.1}$| M_⊙ yr⁻¹.

15

The same result holds even if |$L_{\rm lim, [C\, \small {II}]}=10^{7}~{\rm L}_{\odot }$| (L_{lim, UV} = 10⁹ L_⊙) – a luminosity threshold reasonably achievable with current facilities – where the median [C ii] (UV) luminosities of the satellites are 3.6 × 10⁶ and 7.8 × 10⁶ L_⊙ (4.2 × 10⁹ and 1.8 × 10¹⁰ L_⊙) in noAGN and AGNcone, respectively.

16

This number refers to satellites observed at Δv ≤ 1000 km s⁻¹ from the central quasar. The same number increases up to 19 for Δv ≤ 2000 km s⁻¹.

17

The companions in V20 are observed in a field of view of π95² kpc² on the sky plane, i.e., where the sensitivity of the primary beam equals 0.2 times the peak. We consider here the companions observed within ±1000 km s⁻¹ from the quasar, corresponding to about 2.8 Mpc at z ∼ 6.

18

These are the objects in a volume |$V_{\rm A}=\pi R_{\rm fov}^2 L \approx (200)^{3}$| kpc³, where R_fov = 95 kpc is the radius of the ALMA field of view (fov) and L = 280 kpc is the displacement along the line of sight.

APPENDIX A: ACCRETION-WEIGHTED CENTRE

Given the large number of active MBHs in the simulation, each one characterized by different positions, duty cycles, and accretion rates, the overall AGN activity in the galaxy proto-cluster is the result of a complex spatial distribution and time-dependent superposition of individual, partially overlapping, AGN activities. We thus define, for each snapshot of the simulation AGNcone, a centre of reference that takes into account both the position and the accretion rate of all the AGN in the simulation. We refer to this point in the simulation volume as the ‘accretion-weighted center’, c_aw, of the galaxy proto-cluster.

Despite the numerous MBHs in the simulation volume, their luminosity remains almost irrelevant for most of the time. In particular, at z < 10, beyond the four most accreting MBHs, rarely |$\dot{M}_{\rm acc}\gt 0.1~\!{\rm M}_\odot$| yr⁻¹. However, in order to be more conservative, we compute c_aw, taking into account the eight most accreting MBHs in the simulation volume, weighted for the accretion rate of each considered MBH:

$$\begin{eqnarray*} {\boldsymbol c}_{\rm aw}(t) = \frac{\sum _{i=1}^8 \dot{M}_{{\rm acc}, i}(t-\Delta t_{\rm snap}){\boldsymbol r}_{i}(t-\Delta t_{\rm snap})}{\sum _{i=1}^8 \dot{M}_{{\rm acc}, i}(t-\Delta t_{\rm snap})}, \end{eqnarray*}$$

(A1)

where Δt_snap ≃ 2 × 10⁷ yr is the minimum time spacing between two snapshots in the run; in other words, to compute c_aw at a given snapshot (i.e. at a given time t), we consider the accretion rate and the MBH positions in the previous snapshot (i.e. at the time t − Δt_snap). This method allows us to take into account both the time delay between the launch of the outflow and the arrival of the outflowing gas on the satellites, and the time necessary to develop the possible feedback-driven effects on the satellite properties (e.g. M_* and SFR). Indeed, if we consider, for a typical galaxy in this study, a scale distance from the feedback source of about 100 kpc and an initial wind speed of 10⁴ km s ⁻¹, the time scale needed by the outflow to reach the target galaxy is Δt_wind > 10⁷ yr ∼Δt_snap. Moreover, the typical time-sale needed to convert half of the dense gas mass into stars is Δt_SF ∼ 10⁷–10⁸ yr ∼Δt_snap (see e.g. Gao & Solomon 2004; Wu et al. 2005; Hartmann 2009).

The determination of c_aw occurs in two steps: at first (i) the eight most accreting MBHs are found in each snapshot and then (ii) their location is followed in the next snapshot by selecting either the position of the MBH with the same id, if the MBH is still present, or the position of the descendent MBH, if a coalescence event occurred in the time interval between the two consecutive snapshots.

To find a proper counterpart of c_aw in the noAGN run, where we clearly cannot apply the method described above, we proceed as follows: (i) we consider the positions of the MBHs that contributes to c_aw; (ii) we assign a DM halo to each MBH, only on the basis of its position; (iii) we identify the corresponding halo in the noAGN run by cross-matching the DM particle IDs in the two runs, selecting the halo which shares the largest fraction of particles; (iv) in this run, we choose the halo that maximizes the ratio N_s/N_DM, where N_s is the number of particles shared among the haloes and N_DM is the total number of DM particles in the halo of noAGN run; (v) we save the centre of mass of the selected halo in noAGN, as the equivalent position of the MBH in AGNcone; (vi) the position of c_aw is finally determined through equation (A1), by setting |$\dot{M}_{{\rm acc}, i}$| to the corresponding value for the ith MBH in AGNcone.

Fig. A1 shows the location of the c_aw at z = 6.3, overimposed on the SFR density map, and marked with a red cross; the most accreting MBHs are instead denoted by red filled circles.

Figure A1.

SFR density map of the central 100 kpc at z = 6.3 for noAGN (left-hand panel) and AGNcone (right-hand panel). The accretion-weighted mean centre is marked in each panel with a red cross, along with the position of the most accreting MBHs (⁠|$\dot{M}_{\rm acc}\gt 1~ \!{\rm M}_\odot$| yr⁻¹).

APPENDIX B: OBSERVATIONAL PROPERTIES

In this appendix, we describe the modelling adopted to compute the [C ii] and UV luminosities of z ∼ 6 quasar companions.

B1 [C ii] luminosity

The expected [C ii] luminosity for the satellite sample is evaluated through the relation

$$\begin{eqnarray*} \log _{10}{(L_{[C\, \rm \small {II}]}/{\rm L}_\odot)} = 7.0 + 1.2s + 0.021t + 0.012st - 0.74 t^2, \end{eqnarray*}$$

(B1)

where |$s = \log _{10}{(\rm SFR/\!{\rm M}_\odot yr^{-1})}$| and t = log₁₀(Z/Z_⊙). This relation, provided by Yue et al. (2015), is based on the Vallini et al. (2013, 2015) model that has been tested against data of z ∼ 6 galaxies and well reproduces the [C ii]-SFR relation at high-z (Pallottini et al. 2017, 2019). However, it has the limit of not considering the effect of X-ray on the [C ii] emission (e.g. Meijerink, Spaans & Israel 2007; Langer & Pineda 2015).

B2 UV luminosity

We compute the satellites expected UV luminosity at 1450 Å by considering star emission only in noAGN, and both stellar and AGN contributions in AGNcone, i.e. |$L_{\rm UV}=L_{\rm UV}^{\rm star}+L_{\rm UV}^{\rm AGN}$|⁠, where

$$\begin{eqnarray*} L_{\rm UV}^{\rm star} = 2.2\times 10^{43}~{\rm SFR}, \end{eqnarray*}$$

(B2)

according to Kennicutt & Evans (2012), and

$$\begin{eqnarray*} L_{\rm UV}^{\rm AGN}=f_{\rm UV}~L_{\rm bol}. \end{eqnarray*}$$

(B3)

f_UV is the bolometric corrections provided by Shen et al. (2020), |$L_{\rm bol} = \epsilon _{r}(1-\epsilon _{f})\dot{M}_{\rm acc}c^{2}$| is the bolometric luminosity, and |$\dot{M}_{\rm BH}$| is the sum of the accretion rates of every MBHs within βr_vir.

After this, we consider dust attenuation by applying an extinction coefficient |$e^{-\tau _{\rm eff}}$| to the UV emission. In particular, we consider the minimum and maximum values of the optical depth, as they are evaluated by Di Mascia et al. (2021a) on the same AGNcone and noAGN runs, through radiative transfer calculations, along six different lines of sight. We then average these values (τ = 1.8−2 for noAGN and τ = 1.5−2.7 for AGNcone, with a dust-to-metal ratio f_d = 0.08; Behrens et al. 2018) and apply the resulting τ_eff to each satellite.