Active Galactic Nucleus Feedback in an Elliptical Galaxy with the Most Updated AGN Physics: Parameter Explorations

In a previous work, we have proposed a sub-grid model of active galactic nucleus (AGN) feedback by taking into account the state-of-the-art AGN physics, and used that model to study the effect of AGN feedback on the evolution of an isolated elliptical galaxy by performing two dimensional high-resolution (i.e., the Bondi radius is well resolved) simulations. In that work, typical values of model parameters were adopted. In the present work, we extend that study by exploring the effects of uncertainties of parameter values. Such a study is also useful for us to understand the respective roles of various components of the model. These parameters include the mass flux and velocity of AGN wind and radiative efficiency in both the hot and cold feedback modes, and the initial black hole (BH) mass. We find that the velocity of AGN wind in the hot mode is the most important quantity to control the typical accretion rate and luminosity of AGN, and the mass growth of the BH. The effect of the wind on star formation is less sensitive. Within the limited parameter range explored in the current work, a stronger AGN wind suppresses star formation within ~100 pc but enhances star formation beyond this radius, while the star formation integrated over the evolution time and the whole galaxy roughly remain unchanged. AGN radiation suppresses the BH accretion in a mild way, but dust is not considered here. Finally, a smaller initial BH mass results in a more violent evolution of the BH accretion rate. The corresponding AGN spends more time in the high-luminosity state and the percentage of BH mass growth is higher. Our results indicate the robustness of AGN feedback in keeping the galaxy quenched.


INTRODUCTION
In the centre of every massive galaxy with a bulge there exists a supermassive BH (see, e.g., Kormendy & Ho 2013 for a review). Observations have found tight correlations between the mass of the BH and the properties of the classical bulge, including its stellar mass (Magorrian et al. 1998;Häring & Rix 2004;Kormendy & Ho 2013), luminosity (Marconi & Hunt 2003;Gültekin et al. 2009), and stellar velocity dispersion (Gebhardt et al. 2000;Ferrarese & Merritt 2000;Tremaine et al. 2002). These correlations suggest the coevolution of the central BH and its host galaxy, and AGN feedback might play a major role in this coevolution.
The literatures regarding AGN feedback has increased greatly in the past twenty years (Harrison et al. 2018). Due to the complexities of this relatively young topic, AGN feedback is studied mainly ★ E-mail: yzy1116@shao.ac.cn † E-mail: fyuan@shao.ac.cn by numerical simulations. The gap of the physical scales between the BH and the galaxy can be as large as nine orders of magnitude, indicating all the current simulations require sub-grid assumptions and approximations. Di  and Springel et al. (2005) first studied AGN feedback by using hydrodynamical cosmological simulations. Although the AGN feedback was simply employed through a free parameter (the feedback efficiency f , which determines the fraction of the AGN bolometric luminosity that couples to the gas near the BH), they found a tight correlation between the BH mass and the stellar velocity dispersion. On the other hand, by including AGN feedback, both semi-analytic modeling (Croton et al. 2006;Bower et al. 2006; Monaco et al. 2007) and hydrodynamical cosmological simulations (Vogelsberger et al. 2014;Dubois et al. 2014;Khandai et al. 2015;Schaye et al. 2015;Pillepich et al. 2018) obtained much better fittings to observations, e.g., the galaxy luminosity function at the massive end of the galaxy mass, the "downsizing problem", and the "cooling flow problem". Therefore, it is generally believed that AGN feedback does play an important role in the galaxy evolution. One of the most important quantities for the study of AGN feedback is the mass accretion rate of the BH, because it determines the power of AGN. However, in cosmological simulations, the scales relevant to the BH accretion and the ejection of AGN winds are not directly resolved due to resolution limitation. The BH accretion rate usually is calculated by the gas ∼ 1 kpc from the central BH, which could overestimate or underestimate the real BH accretion rate with the uncertainty as large as ∼ 300 (Negri & Volonteri 2017). Moreover, the AGN outputs such as radiation and matter outflows are often not speficied and the interactions between these AGN outputs and the gas in the galaxy are usually treated in a phenomenological parameterized approach.
Since AGN feedback occurs in a single galaxy rather than cosmological scale, to investigate the details of how AGN feedback works, a perhaps better approach is zooming in on a galaxy with high resolution which can resolve the Bondi radius B = 2 BH / 2 s,∞ , which is roughly ten pc for the BH mass of 10 9 and gas temperature of 10 8 K (Ciotti & Ostriker 1997Ciotti et al. 2009;Shin et al. 2010;Ciotti et al. 2010;Novak et al. 2011;Gan et al. 2014;Ciotti et al. 2017;Yuan et al. 2018, hereafter Paper I;Zeilig-Hess et al. 2019). Within the Bondi radius, the gravity is dominated by the BH, so this is the regime of accretion flow.
Various types of accretion disks have been well investigated by the accretion disk community, including cold accretion mode at high accretion rate regime (e.g. Frank et al. 2002) and hot accretion mode in the low accretion rate regime (e.g. Yuan & Narayan 2014). The cold accretion mode correspond to radiative or quasar feedback mode, while the hot accretion mode corresponds to maintenance or radio or kinetic feedback mode. Some of these names are confusing or even misleading so in Paper I we suggest to call them cold and hot feedback modes respectively. Since the Bondi radius is commonly regarded as the outer boundary of accretion flow, once it is resolved as it is in our simulations, we can precisely calculate rather than estimate the accretion rate. Given the BH accretion rate, we can then adopt the accretion knowledge to calculate the outputs, including radiation and matter outflow, and further calculate their interaction with the gas in the host galaxy.
Taking an isolated elliptical galaxy as an example, Paper I investigated the effects of AGN feedback by establishing a sub-grid model of AGN feedback with state-of-the-art accretion physics. The inner radius of the simulation domain is several times smaller than the Bondi radius. The mass accretion rate is thus precisely calculated and the accretion (and feedback) mode can be determined. Both radiation and momentum-driven wind were considered in the two modes. The properties of the wind in the cold mode were adopted from observations, while in the hot mode they were described based on the 3D GRMHD simulation results of Yuan et al. (2015) due to the rarity of observational data. The jet in the hot mode was temporarily omitted.
In Paper I, they focus on the case of low angular momentum of the galaxy. They examined the respective roles of radiation and wind feedback, and found that both can suppress star formation and cause the variability of the AGN, but the wind was by momentum interaction while radiation was by radiative heating. Wind was believed to play a more important role in suppressing star formation and the BH accretion rate. In the second paper of this series, Yoon et al. (2018) (Paper II) extended this model to the high angular momentum case of the elliptical galaxy. They found that while some results were qualitatively similar to those in Paper I, other results, such as star formation and black hole growth, showed a significant difference due to the mass concentration in the galactic disk as a consequence of galactic rotation. More recently, Yoon et al. (2019) specifically focused on the role of hot mode feedback. They find that although the AGN power in the hot mode is much lower compared to the cold mode, the hot mode feedback still plays an important role, because most time of the AGN stays in this mode.
One remaining question in Paper I is that, since the parameters of the sub-grid model have some uncertainties, to what extent do these uncertainties effect the galaxy evolution. We will accomplish this goal in this paper. The paper is structured as follows. In Section 2, we briefly introduce the framework of our models, including the galaxy setups, stellar feedback, star formation, and AGN feedback. In Section 3 we introduce the parameters explored in this paper. Then in Section 4, we show our results. Finally, in Section 5, we summarize our results and compare with the previous works.

MODELS
In this section, we briefly introduce the main physical processes included in our models. The simulation begins with a massive elliptical galaxy at the age of 2 Gyr. As in Paper I we adopt the sub-grid physics that divides the accretion and feedback into two modes depending on the accretion rate at the inner boundary. In each mode we consider both wind and radiation, which are then injected into the simulation region through inner radial boundary. We simulate the interactions of wind and radiation with the interstellar medium (ISM) by considering simplified radiative transfer. Stellar processes such as star formation, stellar winds, supernovae (SNe) Ia and II are taken into account as well. Readers are referred to Paper I for more information.

Galaxy Model
We focus on the secular evolution of a massive isolated elliptical galaxy. The gravitational potential we adopt is dominated by the dark matter halo beyond 10 kpc, by stars from 0.05 to 10 kpc, and by the central BH within 50 pc, respectively.
Both the dark matter and stellar components are modeled through static potentials, and the potential from newly formed stars can be neglected since they are minor compared to other components. We adopt the Jaffe stellar distribution (Jaffe 1983) embedded in the dark matter halo so that the total density satisfies the isothermal sphere assumption and decreases as −2 (Ciotti et al. 2009). The Jaffe stellar distribution is described by where ★ and ★ are the total stellar mass and the scale-length of the galaxy, respectively. ★ is set to be 3 × 10 11 and ★ is the scale length that corresponds to the effective radius e = 0.7447 ★ . The total density profile is given by where 0 = 260 km s −1 is the central projected velocity dispersion. From the Faber-Jackson relation and the Fundamental Plane, we can derive the total B-band luminosity B = 5 × 10 10 B, and the effective radius e = 6.9 kpc. The initial mass of the central BH is determined according to the empirical correlation BH /(10 9 ) = 0.49 ( ★ /10 11 ) 1.17 in Kormendy & Ho (2013), thus yielding the BH mass BH = 1.8 × 10 9 for ★ = 3 × 10 11 . The initial gas fraction is negligible, and most of the gas is provided by stellar evolution in our work. The galaxy is assumed to be slowly rotating, which is supported by observations that suggest slow rotators start appearing when ★ > 2 × 10 11 (Cappellari et al. 2013;Graham et al. 2018). Since the stars are slowly rotating and the stellar wind is the only mass source in our simulations, the angular momentum of the gas is therefore low and we do not need to handle the angular momentum transfer (for the case of high angular momentum, readers are referred to Paper II; Gan et al. (2019)). To focus on the effect of model parameters, the gaseous halo and cosmological inflow are not considered in this paper, consistent with Paper I.

Star formation and Stellar Feedback
We implement star formation by subtracting mass, momentum, and energy from the grid. We note that stellar populations and their dynamics are not explicitly tracked in the simulation. Different from Paper I, the star formation is triggered only if the temperature is lower than 4 × 10 4 K and the number density is higher than 1 cm −3 concurrently. The aim is to mimic the surface density threshold of ∼ 10 M pc −2 for star formation revealed by observations (Kennicutt 1989(Kennicutt , 1998Martin & Kennicutt 2001;Bigiel et al. 2008). The star formation rate per unit volume is given by the Kennicutt-Schmidt prescription ( Here the star formation efficiency is SF = 0.01, which is an order of magnitude smaller than in Paper I. This modification is driven by the observations on surface density relationships of galaxies (Kennicutt 1998), and local giant molecular clouds (Krumholz & Tan 2007; for a comprehensive discussion see the review by Krumholz et al. 2019). The star formation timescale is The cooling and the dynamical timescales are given by Here K ( ), , and are the Keplerian velocity at radius , internal energy density, and net cooling rate per unit volume, respectively. We adopt the formulae in Sazonov et al. (2005) to compute cooling, which includes the bremsstrahlung cooling, Compton cooling, line and recombination continuum cooling. We note that, as in many similar numerical simulation works of AGN feedback, our calculation of the star formation rate has big uncertainties. It is technically difficult to simulate the process of star formation from first principle in such large-scale simulations. These uncertainties and difficulties, including our neglect of the self-gravity of the ISM, can be absorbed in some degree in our simplified and parameterized treatment of our calculation of star formation rate.
The evolving stars inject mass and energy to the ISM, mainly during the asymptotic giant brach (AGB) phase. At the end of their lives, drastic outbursts called supernova feedback produce a large amount of energy and return most of their mass to the ISM. Based on the population of stars, SNe Ia and II are distinguished. All of these stellar feedback processes are considered in our simulations and the detailed descriptions can be found in Ciotti & Ostriker (2012). The chemical evolution and dust absorption are ignored for simplicity.

AGN feedback
The AGN feedback is divided into two modes depending on the accretion rate at the innermost grid radius. The boundary of the accretion rate can be inferred from the observations on the state transition of black hole X-ray binaries, which occurs under the critical luminosity C ∼ 2% Edd , or the equivalent critical accretion rate C ∼ 2% Edd . In principle, this critical accretion rate applies to the accretion rate at the BH horizon BH . Yet in practice, we simply rely on the accretion rate at the inner boundary ( in ) to judge the feedback mode, which is calculated by Both radiation and wind are taken into account in each mode. The radiative transport is considered in a approximated way by assuming the flow is optically thin Ciotti & Ostriker (2012). The heating terms include the Compton heating and photoionization heating driven by the central AGN. The radiation pressure is included by considering both the electron scattering and the absorption of photons by atomic lines. Wind is input via momentum, which was found to be more powerful than the thermally driven wind and better consistent with observations by previous works (Choi et al. 2012(Choi et al. , 2015. In the hot mode, the radiative efficiency used when calculating the radiation flux of AGN is a function of accretion rate (Xie & Yuan 2012), Here cold /0.057 accounts for the spin of the BH and Edd = Edd /( cold 2 ) is the Eddington accretion rate. The values of 0 and depend on BH and , which denotes the fraction of the viscously dissipated energy that directly heats electrons. Assuming = 0.1, 0 and are given by (0.12, 0.59), BH / Edd 9.4 × 10 −5 (0.026, 0.27), 9.4 × 10 −5 BH / Edd 5 × 10 −3 (0.5, 4.53), 5 × 10 −3 BH / Edd 6.6 × 10 −3 (0.057, 0). 6.6 × 10 −3 BH / Edd 2 × 10 −2 The Compton temperature in the hot mode, used to calculate the radiative heating by AGN to the ISM, is calculated based on the spectral energy distribution of low-luminosity AGNs combined from the literature (Xie et al. 2017) and is given by C,hot = 10 8 K, 10 −3 BH / Edd 0.02 5 × 10 7 K. BH / Edd 10 −3 In the hot mode, the accretion flow consists of an inner hot accretion flow plus an outer truncated thin disk. The truncation radius tr is described by (Yuan & Narayan 2014) where s is the Schwarzschild radius. The existence of a strong wind in hot accretion flows has been shown in Yuan et al. (2012). Based on 3D GRMHD numerical simulation of black hole accretion, using the "virtual particle trajectory" approach, Yuan et al. (2015) derived the mass flux and velocity of the wind: where K ( tr ) is the Keplerian velocity at tr . The velocity is ∼ 2000 km/s when ( in ) = 10 −3 Edd , and it increases to ∼ 0.08 as the ( in ) approaches 0.02 Edd . Given we can obtain the properties of the BH accretion rate and wind properties in the hot mode. Compared to jets, the opening angle of the wind is much larger, lying within ∼ 30 • − 70 • and 110 • − 150 • above and below the equatorial plane, respectively (Yuan et al. 2015). The mass flux of wind within the abovementioned tworanges is assumed to be independent of angles.
In the cold mode, the accretion rate is high and a standard thin disk model is applied. The radiative efficiency is commonly assumed cold = 0.1 (Yu & Tremaine 2002;Marconi et al. 2004), which means the BH is moderately spinning according to the thin disk model. The Compton temperature C,cold , which measures the average energy of the emitting photons, is calculated from the observed spectrum of quasars (Sazonov et al. 2004) and is given by 2 × 10 7 K.
The wind properties in the cold mode are obtained from Gofford et al. (2015). They analyzed a sample of 51 Suzaku-observed AGNs and independently detected Fe K absorption in 40 percent of the sample. After processing the data, they measured the mass flux and velocity of the wind: Obviously, the wind should not be isotropic, but the exact description of the distribution of the wind flux is still poorly constrained. Following the previous works (Novak et al. 2011;Gan et al. 2014;Ciotti et al. 2017), the mass flux of the wind is assumed to be proportional to cos 2 ( ).
To obtain the properties of the wind in the cold mode, we need to calculate the BH accretion rate and the AGN luminosity. Once the gas reaches the circularization radius cir , the accretion disk is formed. With the total mass of the gas in the disk, dg , the mass inflow rate at cir in the disk can be estimated as where the viscous timescale vis is described by (Kato et al. 2008) vis ≈ 1.2 × 10 6 yr 0.1 −1 cir 100 s 7/2 BH 10 9 .
Here is the viscosity parameter. Given we can obtain the BH accretion rate and the wind properties in the cold mode. We would like to point out that, when ( in ) is above 2% Edd , the above approach of calculating BH does not ensure BH is higher than 2% Edd , because some inflowing gas may be depleted via star formation or circularize and fall in at a slower rate. It is our future plan to improve the sub-grid physics adopted here.
Finally, an issue worth noting is whether the wind driven from the accretion flows is able to reach large scales and be injected to our simulation region.  and  study the large-scale dynamics of wind launched from hot accretion flow and thin disks in the cases of without and with magnetic field via analytical and numerical methods. They find that even when magnetic field is not included, wind in the hot mode can reach very large scales while wind launched by thin disk stops at a smaller distance due to its smaller Bernoulli parameter. They do not consider the radiation pressure, which might be the essential mechanism to drive the wind in the cold mode (King & Pounds 2003, 2015Costa et al. 2018), and the wind terminal velocity will be higher with the assistance of magnetic fields . Since in general the wind properties are function of radius, we have made sure that the properties adopted in the present work is suitable for the radius of our inner boundary of our simulation, ∼ in .

Simulation Setup
We use the parallel ZEUS-MP/2 code (Hayes et al. 2006) to study an isolated elliptical galaxy by adopting two-dimensional axisymmetric spherical coordinates. The simulations start at 2 Gyrs after the Big Bang to avoid the early stage of galaxy formation when the major merger plays a dominant role, since we only focus on an isolated elliptical galaxy. The initial conditions for the ISM are set by the very low density gas at the local thermalization temperature. The distribution of the gas is not important since the mass is dominated by the stellar winds as simulations begin. The mass injection of stellar winds reaches the peak at the beginning of the simulation, then decreases gently till = 0, which lasts about 12 Gyrs. The mesh in the polar direction is divided homogeneously by 30 grids. In the radial direction the simulation region is resolved by 120 grids to cover the range of 2.5 pc − 250 kpc, with the finest resolution up to ∼ 0.3 pc by adopting a logarithmic mesh.
How is this inner boundary compared to the scale of accretion flow, i.e., the Bondi radius? In our simulation, the gas is in general multi-phase, and consists of both inflowing and outflowing material. Only inflowing gas is taken into account when we estimate the Bondi radius since only these gas contribute to the accretion. For inflowing gas, each phase has its corresponding Bondi radius. We have calculated the mass flux-weighted value of Bondi radius by including all phases close to our inner boundary. The calculated Bondi radius as a function of time is shown in Figure 1. As we can see from the figure, the Bondi radius is about 5 times larger than the radius of our inner boundary. We have also calculated the Compton radius, which is defined by It is larger than the inner boundary radius for C = 10 8 K. This indicates that the Compton effect mainly plays a role of cooling inside the inner boundary. We carve a narrow range of (∼ 9 • ) at each pole to avoid the singularity there and also adopt the "outflow boundary condition" for polar boundaries. The AGN sub-grid model works as follows. After measuring the accretion rate at the inner boundary, we can calculate the BH accretion rate and wind properties by adopting the formulas in the corresponding mode. Then the BH mass is updated and the wind is injected in the inner grid within the opening angle with mass and momentum conservation. A temperature floor of 10 4 K is set in the cooling functions due to the considerations of spatial resolution in our simulations.

MODEL PARAMETER EXPLORATIONS
In this section, we introduce the explorations of our model parameters. We perform surveys on wind and radiation of the AGN in both hot and cold modes. To avoid confusion, each model of wind and radiation alters one parameter while keeping other parameters unchanged to control variables. Apart from that, we explore the case of a lower initial BH mass. Note that in this case the Eddington accretion rate decreases due to the lower BH mass. Therefore, the same mass flux of accretion flows yield higher Eddington ratio. All the surveyed parameters are summarized in Table 1. These parameters are explored within the constraints imposed by observations or theories, as detailed below.

Wind
In the hot mode due to the scarce of observational data, we model the wind properties by adopting simulation results of Yuan et al. (2015). This work only deals with non-spin black hole and normal magnetic field (SANE). On the basis of Yuan et al. (2015), Yang et al. (2020, in prep.) study the effects of magnetic field and BH spin on wind properties by performing new GRMHD simulations. Their results suggest that the magnetic field and BH spin can affect the wind properties significantly. Wind velocity in the case of strong magnetic field (MAD00, here the number denotes the spin of BH) is about three times higher than that in its weak magnetic field counterpart (SANE00) obtained in Paper I. HotWindVel is the model to explore this high velocity case. The mass flux of wind of MAD00 is given by  Gofford et al. (2015) and the XMM-Newton sample in Tombesi et al. (2012) respectively. The solid lines are the best fittings to the green data points given in Gofford et al. (2015). The dashed and dotted lines vary the best fittings by a factor of ten (top) and three (bottom), representing the ColdWindFluxHigh/Low models and the ColdWindVelHigh/Low models, respectively. The magenta and gray dotted lines represent the critical luminosity (2% Edd ) for the Fidu and BHmass models. This is several times lower than that of SANE00 when the mass accretion rate approaches 10 −2 Edd . On the other hand, the mass flux of the high BH spin model SANE98 given by is three times higher than that of SANE00. It is worth noting that when ( in ) is small, all three models yield W close to ( in ), because BH of all three models is orders of magnitude smaller than W . We perform the HotWindFluxLow and HotWindFlux-High models to simulate these two cases of MAD00 and SANE98, respectively. The mass flux of the HotWindVel model and the wind velocity of the HotWindFluxHigh/Low model remain unchanged under the same ( in ). This indicates the lower wind density of the HotWindVel model and higher wind density of the HotWindFlux-High/Low model given the definition of mass flux  where is the constant area that the wind blows out of the inner boundary.
In the cold mode we estimate wind properties through observations on broad absorption line (BAL) outflows (Arav et al. 1999;Tombesi et al. 2010Tombesi et al. , 2012Arav et al. 2013;Gofford et al. 2013Gofford et al. , 2015. Paper I adopts the relation between wind properties and bolometric luminosity from Gofford et al. (2015). Figure 2 compares this relationship with the observational data on UFOs (ultra-fast outflows). In addition to the data from Gofford et al. (2015), the data from Tombesi et al. (2012) is also shown in the figure. It is worth noting that because the wind injection radius, i.e., the inner boundary, in our sub-grid model is within the Bondi radius, we only choose observations of UFOs that are on sub-pc scales. From Figure 2, it is easy to see that the uncertainties of the wind flux and velocity are considerable. Therefore, we perform four models ColdWindFluxHigh/Low, ColdWindVelHigh/Low, trying to cover the scatters in the observations. The flux-variant models vary the mass flux by an order of magnitude and the velocity-variant models alter the wind velocity by a factor of three. Note that an upper limit of 10 5 km/s is applied in all the models. Similarly, the mass flux of the ColdWindVelHigh/Low model and wind velocity of the ColdWindFluxHigh/Low model remain identical.

Radiation
In addition to wind properties, the effect of radiation is studied as well in both hot and cold modes. The HotRad model studies the effect of radiation in the hot mode. As shown in Section 2.3, the radiative efficiency is a function of , which denotes the fraction of viscously dissipated energy that directly heats electrons. The value of is constrained to be ∼ 0.1 − 0.5 (Yuan & Narayan 2014), stronger radiation corresponding to larger = 0.5 is considered in the HotRad model. Figure 3 shows the radiative efficiencies corresponding to two different as a functions of BH accretion rate. The radiative efficiency corresponding to = 0.5 can be an order of magnitude higher than that of = 0.1 when the BH accretion rate is low. The radiative efficiency of the HotRad model in the cold mode remains same.
In the cold mode, the ColdRadHigh model explores the high spin of the BH with radiative efficiency cold = 0.3, while the ColdRadLow model investigates the zero spin of BH, of which the radiative efficiency cold = 0.057. In order to study the effects of radiation in the cold mode, we keep the the radiative efficiency in the hot mode identical to that of the Fidu model. Moreover, in order to study radiation, the wind properties remain unchanged under the same ( in ).

Initial BH mass
Other than wind and radiation properties, we explore the effect of the initial BH mass using the BHmass model. The initial BH mass adopted in Paper I is BH,i = 1.8 × 10 9 for ★ = 3 × 10 11 (Kormendy & Ho 2013), but recent works suggest the discrepancy of BH mass measurement between AGN hosts and ellipticals/classical bulges (Ho & Kim 2014;Reines & Volonteri 2015;Shankar et al. 2016Shankar et al. , 2019. Shankar et al. (2016) use Monte Carlo simulations to illustrate the selection bias in local, dynamically measured BH samples of ellipticals/classical bulges, in which only the massive BHs are measured. Therefore, the BH mass in Kormendy & Ho (2013) may be overestimated. In their model the BH mass is BH,i = 2.7 × 10 8 for ★ = 3 × 10 11 . The BHmass model studies this possibility, with the initial BH mass in this model being six times lower than that in Paper I.

RESULTS
In this section, we present the results of our parameter explorations. We essentially investigate the effects of parameters in four major domains of galaxy evolution: AGN luminosity, BH mass growth, star formation, and the AGN duty cycle. The fiducial model is essentially identical to the fullFB model in Paper I but with two improvements as discussed in Yoon et al. (2019). One is the calculation of star formation as described in Section 2.2, the other is the correction of a bug in computing the energy flux of the wind. This bug caused a lower energy flux of the wind in the hot mode while a higher energy flux in the cold mode. Figure 4 shows the evolution of the accretion rate at the inner boundary ( in ), BH accretion rate BH , and bolometric luminosity bol for all models. We note that the time interval of two adjacent data points is ∼ 1 Myr for clarity, thus some outbursts are filtered out in this case. For the fiducial model, ( in ) oscillates around ∼ 10 −3

AGN luminosity
Edd . The value of this "baseline accretion rate" is determined by the momentum balance between the AGN wind and the inflow at the inner boundary of our simulation, reflecting the characteristic of the hot mode in our AGN sub-grid model. When the accretion rate decreases, both the mass flux and the velocity of the wind become lower and yield a lower momentum of the wind. According to our sub-grid model, the decrease of wind momentum is more rapid than that of the inflow, thus making the accretion rate increase to reach momentum balance. Oppositely, when the accretion raises, the wind momentum grows faster than the momentum of inflow, thus making the accretion rate decrease. Therefore, we can see small oscillations of ( in ) at ∼ 10 −3 Edd . However, if the inflow is massive enough to overcome the momentum of wind, such as driven by the cold clumps, the accretion rate is able to enter the cold mode and trigger a strong AGN wind and radiation of the cold mode. The wind and radiation push the accreting gas outwards so the mass accretion rate at the Bondi radius drops drastically. In addition, they interact with the ISM of the nuclear region and heat the surrounding gas. The subsequent galactic wind can reach several kpc scale, but most of the wind can not break out of the halo. The compression of the kicked gas is likely to form the cold clumps once more, yet in a moderate way because of the lower density of the ISM. These cold clumps fall back to the centre and initiate another episode of the AGN outburst. Eventually, the surrounding gas is almost cleared out and the AGN returns back to be quiescent. Because the mass supply from the stellar wind decreases with time, the massive inflow characterized by the peaks of accretion rate occurs more frequently at the early stage of the evolution due to abundant supplies from the stellar wind. Actually, the accretion rate stays in the hot mode after ∼ 6 Gyr. However, the accretion rate for the momentum balance maintains the same level throughout the evolution, keeping the typical mass accretion rate unchanged.
As shown in the middle panel of Figure 4, the BH accretion rate oscillates around 10 −5 Edd , and the amplitude of the oscillation is larger than that of ( in ). Given that ( in ) ∼ 10 −3 Edd , this implies that 99% of the inflows are restored to the simulation region by the strong wind. The right panel shows the light curve of the AGN bol = BH 2 . Overall, we can see that the galaxy spends most of its time in the low-luminosity phase with the bolometric luminosity lying in the range 10 −6 − 10 −4 Edd , which is consistent with observations that a median Eddington ratio of ∼ 10 −5 is found for ellipticals (Ho 2009). Our typical luminosity of Fidu is 1-2 orders of magnitude lower than that of the fullFB model in Paper I, which is ∼ 10 −4 Edd . This is caused by the correction of the bug that under-produced the wind power in the hot mode.
As the parameters vary, Figure 4 shows that ( in ) varies the least among the three physical quantities, which reflects the nonlinear characteristics in our AGN sub-grid models. Here, variations are measured through typical values the curves fluctuate around and the frequency of strong oscillations, such as outbursts and sudden drops to very low accretion rate or luminosity. For the HotWindVel model, ( in ) oscillates around a value a factor of three smaller than that of the Fidu model ( Figure 5), resulting in orders of magnitude lower BH accretion rate and bolometric luminosity 1 . As mentioned above, the oscillations are caused by the momentum balance between the AGN wind and the inflow. When the wind velocity is three times higher, it requires ( in ) three times lower compared to the Fidu model to achieve a new momentum balance between the inflow and wind at in . This is because, a three-times lower ( in ) will result in both the wind flux and velocity three times lower, while the higher wind velocity in the HotWindVel model will compensate for the decrease of velocity. So the overall decrease of momentum of wind is three times, which balances with the decrease of inflow.
We find from Figure 4 that the typical accretion rates of the HotWindFluxHigh and HotWindFluxLow models vary little compared to the Fidu model. This is because, the mass fluxes of wind in these two models differ from the Fidu model only when the accretion rate is close to 10 −2 Edd , as we emphasize in Section 3.1. When ( in ) < ∼ 10 −3 Edd , the mass fluxes of the hot wind for all the models are similar, and are roughly equal to ( in ). Yet the HotWindFluxHigh model has fewer outbursts, while for the HotWindFluxLow model, more frequent outbursts can be seen during its violent evolution. Such a difference is caused by the different mass fluxes of hot wind when ( in ) approaches 10 −2 Edd . For the HotWindFluxHigh model, the wind is stronger than that of the Fidu model; thus with increasing accretion rate, the accretion mode is more likely to stay in the hot mode. For the HotWindFluxLow model, the wind is weaker, thus with the increase of accretion rate, it is easier for the accretion to enter the cold mode. Since the wind is stronger in the cold mode, it can cause larger oscillations characterized by the outbursts and suppressions of accretion rate and luminosity.
The HotRad model shows similar ( in ) and BH compared to the Fidu model, but the typical bol is an order of magnitude higher due to the higher radiative efficiency. The results of the four models concerning properties of cold wind exhibit less variation than those of the Fidu model. Comparing the High-suffix models (i.e., ColdWindVelHigh and ColdWindFluxHigh) to the Low-suffix counterparts (i.e., ColdWindVelLow and ColdWindFluxLow) of the cold mode, we can see from the figure that the strong suppression of the accretion rate appears more frequently and has larger amplitudes due to stronger wind feedback. But the effect of radiative efficiency in the cold mode is less obvious. Finally, the accretion rate of the BHmass model shows strong oscillations, which is similar to that of the HotWindFluxLow model. This is because the Eddington accretion rate is lower due to the lower BH mass. Therefore it becomes easier for the accretion to enter into the cold mode given the same amount of mass supply from the stellar wind, which produces a stronger wind. l, we find that the galaxy spends more time in the low-luminosity phase with increasing mass flux and velocity of the hot wind as shown by HotWindFluxLow, HotWindFluxHigh, and HotWindVel. Correspondingly the BH growth and the datum line luminosity decline as we mentioned above. Yet all of these three models spend less time in the cold mode. We suppose that the reason of HotWindFluxLow is unphysical as we have discussed already. Furthermore, all of these three models differ from the observations by Ho (2009), indicating these cases are a rarity. Although the galaxies of HotWindVel and HotWindFluxHigh stay longer time in the hot mode, the emitted energy fraction of the cold mode is higher than that of Fidu. This is because the energy/time shown here is the proportion of the total energy/time. Owing to the fact that they spend more time during the lower luminosity phase, in other words, the datum line luminosity is lower, the total energy they emit is lower than that of Fidu, which causes their cumulative energy above 2% Edd relatively higher. Overall, the variation of the cumulative time and energy above 2% Edd is within a factor of 2.

BH mass growth
The BH mass growth for various models are listed in the second column of Table 2. Generally all the models yield mass growth of less than 4% of the initial BH mass, suggesting that the AGNs are quiescent overall regardless of the parameters we explore. Moreover, the models with stronger wind or radiation have lower BH mass growth, suggesting that both AGN wind and radiation have negative effects on the BH mass growth, which is easy to understand. Specifically, as a result of increasing the wind velocity by a factor of three, the HotWindVel model has lower BH mass growth due to its lower typical BH accretion rate in the hot mode. The HotWind-FluxHigh model reduces the BH mass growth by over a factor of two, while the HotWindFluxLow model has higher BH mass growth due to the lower mass flux of the wind. In the cold mode, while the BHs have over 50% higher mass growth by reducing velocity and mass flux of the cold mode as shown by the ColdWindVelLow and ColdWindFluxLow models, respectively, the BH mass varies little when increasing the velocity and mass flux by the same magnitude, as shown by the ColdWindVelHigh and ColdWindFluxHigh models. This is caused by the small fraction of BH mass growth in the cold mode. Since the growth of BH mass is already dominated by the accretion in the hot mode, stronger feedback in the cold mode reduces little BH growth. However, weaker feedback in the cold mode can increase the accretion during the cold mode significantly, thus resulting in more substantial BH mass variations.
The models with stronger radiation either in the hot mode or the cold mode show negative effects on the BH accretion. Weaker * The fraction of time that galaxy spends above 0.02 Edd . ** The BH mass growth of the Fidu model is an order of magnitude smaller than that of the fullFB model in Paper I. It is reasonable considering an order of magnitude lower typical BH accretion rate for the Fidu model. radiation in the cold mode represented by the ColdRadLow produces higher BH mass growth. But the effect of radiation is smaller than the wind, consistent with Paper I. However, we should be cautious in drawing the conclusion that wind is more important than radiation in controlling the BH accretion, since we have not considered dust in our simulations, which might play an important role in the radiation feedback processes.
Finally, the BHmass model shows much smaller growth of the BH mass than that of the Fidu model. The reasons are twofold. One is that the lower BH mass provides a shallower gravitational potential that reduces the accretion. The other is that the Eddington accretion rate is lower for the BHmass model. Consequently, the accretion is easier to enter the cold mode given the same amount of mass supply from stellar wind, while the cold mode has stronger wind and radiation feedback. However, the percentage of the mass growth of the BHmass model is higher than that of the Fidu model, caused by the higher ratio of stellar mass to the initial BH mass. Figure 6 shows the distribution of the newly born stars for various models. We ignore the gravitational potential provided by the new formed stars and the stellar movements. The results differ significantly from Paper I because we now apply different density and temperature thresholds to star formation, i.e., only the gas with density over 1 cm −3 and temperature under 4 × 10 4 K can form stars, and star formation efficiency SF is reduced by an order of magnitude. Generally, we find that star formation is significantly reduced compared to that of the fullFB model in Paper I, and is highly concentrated inside 1 kpc. The new stars density at the end of the simulations, shown in the left panel of Figure 6, decreases with increasing radius due to the increasing star formation timescale and decreasing gas density. The total mass of new stars is around (4 − 9) × 10 6 for various models, shown in the middle column of Table 2, which accounts for less than 0.01% of the total stellar mass. In fact, the total star formation is lower than the BH mass growth for all the models.

Star formation
We do not find obvious correlations between the wind properties and the total star formation, as shown by the upper right panel of Figure 6, which presents the cumulative mass of new stars integrated over time for models with various wind parameters. In the simulations, the sites of star formation are determined by the distribution of the cold and dense gas. A stronger AGN wind is able to push the gas toward larger radii, and results in the accumulation and condensation of the cold gas there. Thus higher wind power increases star formation at large radii while decreasing it at small radii. On the other hand, wind feedback does strongly suppress the star formation, as shown by, e.g., the right plot of Fig. 8 in Paper I. The absence of the correlation found here is because the range of wind parameters we have explored is relatively too small.
In the center right panel, the HotRad model shows that higher radiative efficiency of the hot mode yields lower star formation inside 100 pc, caused perhaps by stronger radiative heating and radiation pressure; while at large radii, star formation is higher, which is perhaps because the radiation pressure pushes the gas from small radii to large radii, similar to the role played by wind. For radiation in the cold mode, however, stronger radiation suppresses the star formation at both small and large radius, which is likely because the radiative heating is very strong in the cold mode. The bottom right panel shows that the lower BH mass of the BHmass model produces higher star formation overall.
Apart from the spatial properties of star formation, another major difference from Paper I is the temporal evolution of the specific star formation rate (sSFR). Figure 7 shows the evolution of the sSFR for the Fidu model as an example. Throughout the lifetime of the galaxy, the sSFR is weak for most of the time with the time-averaged value of 10 −15 yr −1 . The quiescence of the galaxy, however, is punctuated by sporadic and short episodes of strong star formation outbursts that can reach the sSFR of 10 −12 yr −1 especially at the early stage. This star formation history is strongly correlated with the AGN luminosity light curve. We find the concurrent outbursts in the star formation and the AGN light curve, and the same quiescent time during 7-9 Gyr and 9-11 Gyr. This indicates the tight link between star formation and BH accretion. As mentioned before, the star formation takes place in the cold, dense clumps. If these clumps are not consumed by star formation before they fall into the central BH, a strong AGN outburst is expected afterwards. Both processes are driven by cooling flow infall of gas.

AGN Duty cycle
The duty cycles for different models are shown in the left panel of Figure 8. Each line represents the fractions of time that the galaxy spends above a given Eddington ratio. For the Fidu model, the AGN spends over 99% of its evolution time with the Eddington ratio below 2 × 10 −4 . Comparing to 80% of the total time in the fullFB model of Paper I, the AGN of the Fidu model spends more time in the low-luminosity phase, which is consistent with the lower typical AGN luminosity mentioned in Section 4.1.
In general, wind and radiation of the hot mode exert more remarkable influences on the cumulative time during the lowluminosity phase, while wind and radiation of the cold mode have larger effects above 10 −2 bol / Edd . Specifically, we can see from the top left panel of Figure 8 that the AGN spends more time below the Eddington ratio of 10 −4 bol / Edd as the velocity or mass flux of wind in the hot mode increases. For the Eddington ratio above 10 −2 bol / Edd , however, the AGN spends less time with increasing velocity or mass flux of wind in the cold mode. Particularly, the AGN of the ColdWindFluxHigh model spends little time in the cold mode throughout its evolution. This suggests that AGN wind of both hot and cold modes has negative effects on the AGN duty cycle, which is consistent with the findings that the AGN wind suppresses the BH accretion in Section 4.2.
The center left panel of Figure 8 shows the models with regard to radiation. The AGN in the HotRad model spends less time in the low-luminosity phase compared to the Fidu model. Yet we caution that BH accretion is slightly suppressed when the radiation is stronger as discussed in Section 4.2. The higher AGN luminosity is caused directly by the higher radiative efficiency. From the ColdRadLow model to the ColdRadHigh model, the AGN has decreasing cumulative time above a given Eddington ratio owing to the declining BH accretion rate which overwhelms the effect of increasing radiative efficiency. The bottom left panel exhibits the duty cycles of the BHmass model. The BHmass model shows more cumulative time of the AGN in the high-luminosity phase due to the lower value of Eddington luminosity. The specific numbers of the AGN duty cycles above 0.02 Edd are listed in the right column of Table 2. It is obvious to find that the duty cycles of the cold-mode models have stronger variations than those of the hot-mode models.
The right panel of Figure 8 shows the fraction of total energy emitted above a given Eddington ratio. Compared to the left panel, models reveal minor differences at low Eddington ratios, while the differences are amplified at high Eddington ratios. For the Fidu model, the AGN emits roughly 25% of its total energy in the cold mode. This number is larger than the 6% in the fullFB model of Paper I, which is caused by the correction of the bug that over-produces the energy flux of the wind in the cold mode. However, it is still not consistent with "the Soltan argument", which claims that AGNs emit most of their energy during the high-luminosity phase (Soltan 1982;Yu & Tremaine 2002;Marconi et al. 2004). As already discussed in Paper I, there are two main reasons for this discrepancy. The most important reason is that our simulations begin from 2 Gyr with a massive, mature elliptical galaxy. The central BH mass is already 10 9 with minor growth in its subsequent evolution via accreting mass supplied by the stellar wind. In other words, the essential parts of the BH mass growth and the energy release occur before our simulations begin. The other reason is that we only focus on an isolated galaxy without considering the gaseous halo and cosmological inflow. The fraction of the energy emitted is expected to increase if we include the gas supply from the cosmological inflow and the gaseous halo.
Among all models, the fraction of cumulative energy emitted in the cold mode (the Eddington ratio > 10 −2 bol / Edd ) is the largest for the ColdWindFluxLow model, which is ∼ 70% of the total energy, while for the ColdWindFluxHigh model, almost zero energy is emitted in the cold mode. This suggests that the properties of wind in the cold mode are crucial to the high-luminosity phase of the galaxy. For the models studying the wind of the hot mode, it is interesting to find that although the time spent in the cold mode is roughly the same for the HotWindFluxHigh and HotWindFluxLow models, the AGN of the HotWindFluxHigh model emits a larger fraction of energy in the cold mode than that of the HotWindFluxLow model. This is because, on one hand, we find that the totally emitted energy in the cold mode in the two models are roughly the same; on the other hand, the typical AGN luminosity of the HotWindFluxHigh model in the hot mode is lower than that of the HotWindFluxLow model thus the emitted energy in the hot mode of the HotWindFlux-High model is lower. This is also the reason for the larger fraction of the energy in the cold mode for the HotWindVel model.
In the center right panel, the AGN in the HotRad model emits a smaller fraction of energy in the cold mode than that of the Fidu model because of its higher total energy. The ColdRadLow model has a larger fraction of energy emitted in the cold mode due to its longer cumulative time spent in the cold mode. In particular, almost half of its energy is emitted in the cold mode for the ColdRadLow model, suggesting the non-negligible role of radiation. The Col-dRadHigh model also has a larger fraction of energy emitted in the cold mode, although its AGN spends less time in the cold mode. This is caused by the higher radiative efficiency for the ColdRad-High model. The bottom right panel presents the results for the BHmass model. The AGN of the BHmass model emits nearly half of the total energy in the cold mode.

SUMMARY AND DISCUSSION
In Paper I, we have proposed a sub-grid model of AGN feedback and used it to study the AGN feedback in an elliptical galaxy. In that work, all the model parameters are set to be at their typical values. But some uncertainties still exist although black hole accretion is a relatively mature field compared to AGN feedback. In this paper, we perform simulations to study the effects of parameters in the AGN sub-grid model of Paper I. Such a study is also useful for us to understand the role each model component plays in the feedback. The fiducial model is the updated version of the fiducial model in Paper I. Based on this model, we vary one parameter while keeping other parameters unchanged to study its effect. The models are listed in Table 1. By comparing models to the fiducial model, we are able to find the effect of this corresponding parameter on the AGN and galaxy evolution.
AGN Wind suppresses the BH accretion overall. Particularly, the wind velocity in the hot mode is the most important in controling the typical accretion rate and luminosity of the AGN. For example, a hot wind with three times higher velocity results in the typical AGN luminosity two orders of magnitude lower. When the accretion is in the cold mode, a powerful wind stifles the accretion dramatically and renders the accretion back to the hot mode. This feedbackregulated BH emits at the typical luminosity between 10 −6 Edd and 10 −4 Edd for most of the time, which is consistent with the observations on the nearby, early-type galaxies (Ho 2009). In the case of a stronger AGN wind, the mass growth of the BH decreases, and the AGN spends a larger fraction of time in the low-luminosity phase. Star formation, however, takes place in a more complicated way. A direct correlation between the AGN wind power and the total star formation integrated over cosmological time and the whole galaxy is not found, which is because the range of wind parameters explored here is not large enough.
AGN radiation also suppresses the BH accretion, although not as violently as the AGN wind. But it is premature to conclude that AGN wind plays a more important role than radiation in modulating BH accretion, because we have not considered dust in our simulations. Despite the negative effect on BH accretion, compared to the Fidu model, when the radiation of the hot mode becomes stronger, the AGN duty cycle shows a smaller fraction of time spent in the low-luminosity phase. This is caused by the larger radiative efficiency adopted. Similar to AGN wind, radiation of the hot mode also plays a role in reducing star formation at small radii while enhancing star formation at large radii. For radiation of the cold mode, however, stronger radiation suppresses the star formation at both small and large radii compared to the Fidu model.
Finally, we perform an extra model to investigate the effect of a lower initial BH mass, which is possibly suggested by the obser-vations (Ho & Kim 2014;Reines & Volonteri 2015;Shankar et al. 2016Shankar et al. , 2019. A direct consequence is the stronger oscillations of the BH accretion rate and AGN luminosity. The AGN spends more time and emits more energy in the cold mode, and correspondingly the percentage of the BH mass growth increases. The star formation is also enhanced.
In summary, however, given all the parameters we explore, the variations of the mass growth of BHs and the total star formation are within an order of magnitude, suggesting that our results are relatively insensitive to the parameters of the AGN sub-grid model.
In the initial state of our simulations, the galaxy already contains a very massive black hole and the star formation is already very low. We have shown in this paper (and our previous series of works, e.g., Yuan et al. 2018) that AGN feedback can keep the galaxy quenched. However, the gaseous halo and the cosmological inflow have so far been neglected. When these components are taken into account, it will be important to assess whether AGN feedback can still keep the galaxy quenched. This issue will be discussed in detail in our subsequent paper (Zhu et al. 2020, in preparation). Another issue is the effect of AGN feedback in the high redshift Universe. In that case, the black hole is small and the gas is abundant, so both the activity of the central AGN and star formation are much stronger. It is then interesting to ask what is the effect of AGN feedback; specifically, whether the AGN feedback can quench the galaxies. This question will be investigated in our future works.