Osaka Feedback Model III: Cosmological Simulation CROCODILE

We generate an initial condition using MUSIC2-monofonIC³³3The code’s website is https://bitbucket.org/ohahn/monofonic/src/master/ (Michaux et al., 2021; Hahn et al., 2021). The phase fixing technique (Angulo & Pontzen, 2016) is used to suppress the impact of cosmic variance, but we run only a single simulation for each model and do not include paired simulations for the analyses in this paper. We adopt the following cosmological parameters from Planck Collaboration et al. (2020), ($\Omega_{m}$, $\Omega_{b}$, $\Omega_{\Lambda}$, $h$, $n_{s}$, $10^{9}A_{s}$) = (0.3099, 0.0488911, 0.6901, 0.67742, 0.96822, 2.1064). ⁴⁴4The best-fit parameters of the baseline model 2.20 base_plikHM_TTTEEE_lowl_lowE_lensing_post_BAO_Pantheon of Planck 2018 cosmological parameter table as of May 14, 2019 (https://wiki.cosmos.esa.int/planck-legacy-archive/images/b/be/Baseline_params_table_2018_68pc.pdf). The initial condition is generated at $z_{\rm ini}=39$ using the third-order Lagrangian perturbation theory (3LPT). The transfer function at $z_{\rm ini}$ is obtained by back-scaling the transfer function at the reference redshift $z_{\rm ref}=2$ computed by CLASS ⁵⁵5The code’s website is http://class-code.net/ (Blas et al., 2011). The simulation box size is $L_{\rm box}=50\,{h^{-1}\,{\rm cMpc}}$, and the total number of particles is $N_{\rm p}=2\times 512^{3}$; the initial particle mass of dark matter and gas is $m_{\rm DM}=6.74\times 10^{7}\,{h^{-1}\,{{\rm M_{\odot}}}}$ and $m_{\rm gas}=1.26\times 10^{7}\,{h^{-1}\,{{\rm M_{\odot}}}}$, respectively.

We use the cosmological N-body/SPH code GADGET4-Osaka (Romano et al., 2022a, b), which is a privately developed branch of GADGET-4⁶⁶6The code’s website is https://wwwmpa.mpa-garching.mpg.de/gadget4/ (Springel et al., 2021). The Newtonian gravity for dark matter and baryons is solved using the third-order TreePM method (Xu, 1995; Bagla, 2002; Springel, 2005). The same gravitational softening length is used for all particle types, and it is set to $\epsilon_{\rm grav}=3.38\,{h^{-1}\,{\rm ckpc}}$ but limited to physical $0.5\,{h^{-1}\,{\rm kpc}}$ at all times.

The SPH method (for reviews, see Rosswog, 2009; Springel, 2010a; Price, 2012) is used to solve the governing equations of hydrodynamics. Appendix A describes the numerical treatment of hydrodynamics.

The heating rate, $\Gamma$, and the cooling rate, $\Lambda$, are computed using the Grackle⁷⁷7The code’s website is https://grackle.readthedocs.io/en/latest/ library (Smith et al., 2017), which solves the non-equilibrium primordial chemistry network of 12 chemical species, H, D, He, H₂, HD, and their ions, with radiative cooling by metals and photo-heating and photo-ionization by uniform ultraviolet background (UVB). The metal cooling rate is precomputed using the photoionization code Cloudy (Ferland et al., 2013), and the UVB model by Haardt & Madau (2012) is assumed. In Grackle, we set the parameters related to the onset of UVB to UVbackground_redshift_on=8 and UVbackground_redshift_fullon=6, and then the cosmic reionization occurs at $z\sim 7$. Our simulations do not include the heating by the cosmic rays and stellar radiation, and we ignore the self-shielding of the UVB by atomic hydrogen to avoid overcooling of the neutral gas.

We follow 12 metal elements, H, He, C, N, O, Ne, Mg, Si, S, Ca, Fe, and Ni, produced by core-collapse SNe, type-Ia SNe, and AGB stars. The metallicity floor is $Z_{\rm floor}=10^{-6}\,Z_{\odot}$, where $Z_{\odot}=0.0134$ is the solar metallicity (Asplund et al., 2009).

GADGET4-Osaka solves the production and destruction of dust as described in Romano et al. (2022a) following the full dust grain size distribution; however, we do not focus on dust in this work and omit the description of the dust module here. The analysis of dust will be made in our future work.

We use the Smagorinsky-Lilly type turbulent metal and dust diffusion model (Smagorinsky, 1963; Shen et al., 2010; Romano et al., 2022a) with a diffusion parameter $C_{\rm diff}=2\times 10^{-4}$.

The nonthermal Jeans pressure floor (Hopkins et al., 2011; Kim et al., 2016) is applied to avoid artificial fragmentation as described in the Appendix A.2. The lower limit of the temperature is set to the higher value of the CMB temperature and $T_{\rm min}=15$ K, and the upper limit is set to $T_{\rm max}=10^{9}$ K.

The simulations presented in this article are performed at the SQUID supercomputer at Cybermedia Center, Osaka University, equipped with dual Intel Xeon Platinum 8368 processors (2.4 GHz) per compute node. Running the Fiducial model took $1.12\times 10^{5}$ CPU hours.

In the following subsections, we describe our treatment of subgrid star formation, stellar feedback, and black hole physics. In this work, SNe and AGB stars are considered as sources of stellar feedback, while early stellar feedback, e.g., stellar radiation and stellar wind, is neglected. The adopted values of parameters introduced in the following subsections are summarized in Table 1.

We assume star formation to occur when the hydrogen number density is higher than the threshold density, $n_{\rm thres}=0.1\,{\rm cm^{-3}}$, and the temperature is lower than the threshold temperature, $T_{\rm thres}=10^{4}\,{\rm K}$. Each gas particle is allowed to spawn $n_{\rm spawn}$ star particles at most, and the mass of the star particle is $m_{*}=m_{\rm gas}/n_{\rm spawn}$, where $m_{\rm gas}$ is the mass of the gas particle, and $n_{\rm spawn}=2$.

We use the simple stellar population (SSP) approximation; a star particle represents a cluster of stars with the same metallicity, and their mass function follows an adopted initial mass function (IMF). For the IMF, its dependence on metallicity and redshift is considered based on the star cluster formation simulation by Chon et al. (2022). They have investigated IMF at a metallicity range of $10^{-4}\leq Z\leq 10^{-1}$ and redshift range of $0\leq z\leq 20$ and quantified the mass fraction of excess component from a Salpeter-like component,

with $x=\min\{1+\log_{10}(Z/Z_{\odot}),0\}$, where $z$ is the redshift. The metallicity range investigated by Chon et al. (2022) is limited to $Z\leq 0.1\,Z_{\odot}$, and we use the value of $f_{\rm massive}$ at $Z=0.1\,Z_{\odot}$ for the higher metallicity because IMF at $Z=0.1\,Z_{\odot}$ and $z=0$ is already similar to present day Salpeter-like IMF. Figure 1 shows the variation of IMF depending on $f_{\rm massive}$. We first assume the Chabrier IMF (Chabrier, 2003) with stellar mass range $0.1\,{\rm M_{\odot}}<M<100\,{\rm M_{\odot}}$, and then add a log-flat component at $5\,{\rm M_{\odot}}<M<100\,{\rm M_{\odot}}$ with mass fraction $f_{\rm massive}$ to allow for a treatment of top-heavy IMF.

When a star particle is spawned from a gas particle, $f_{\rm massive}$ is computed from the metallicity of the gas particle and the redshift at that time, which is used later for computing energy and metal yield by SNe and AGB stars. In computing the metal and energy yields, we use the metallicity of the star particle with a floor of $Z_{\rm yield,floor}=0.03\,Z_{\odot}$. This metallicity floor is introduced to consider the metal enrichment by unresolved early star formation, and the $Z_{\rm yield,floor}$ is similar to the observed metallicity of galaxies with $M_{*}\sim 10^{7}\,{\rm M_{\odot}}$.

The star formation rate for a gas particle follows the local Schmidt law (Schmidt, 1959),

where $t_{\rm ff}=\sqrt{3\pi/32G\rho}$ is the local free fall time and $\epsilon_{*}=0.01$ is the star formation efficiency. Star particles are spawned stochastically (Katz, 1992; Springel & Hernquist, 2003; Shimizu et al., 2019) following the star formation rate.

When massive stars ($\gtrsim 8\,{\rm M_{\odot}}$) end their life, their major path is to explode as core-collapse supernovae. Typically, there is one massive star per $100\,{\rm M_{\odot}}$ stars, and the kinetic energy of a supernova is $10^{51}\,{\rm erg}$; the specific SN energy is $\zeta_{\rm SN}=10^{49}\,{\rm erg}\,{\rm M_{\odot}}^{-1}$. However, the specific energy can vary due to multiple factors, e.g., stellar IMF, lower and upper limits of SN progenitor mass, and energy per SN. The specific SN energy adopted in previous works varies by more than a factor of two (Keller & Kruijssen, 2022).

Some supernovae, called hypernovae (HNe) or superluminous supernovae, have an order of higher explosion energy of $10^{52}\,{\rm erg}$. The number fraction of HNe in a solar metallicity environment is about one percent, but it is observationally suggested that the HN fraction increases at subsolar metallicity (Moriya et al., 2018; Gal-Yam, 2019). From a cosmological hydrodynamic simulation of galaxies, Kobayashi et al. (2006) suggested that the HN fraction of 50% is necessary to explain the zinc abundance.

In this work, we use the yield model by Nomoto et al. (2013), which covers the progenitor mass range of $13-40\,{\rm M_{\odot}}$. We assume that some core-collapse SNe with progenitor masses of $20-40\,{\rm M_{\odot}}$ explode as HNe, and its number fraction in that mass range is 50% at $Z<10^{-3}$ and 1% otherwise. We generate a yield table as a function of stellar age using the chemical evolution library CELib (Saitoh, 2017). The bottom panel of Figure 1 shows the specific energy of core-collapse SNe as a function of metallicity in our model. The specific energy is higher for a low-metallicity star due to the top-heavy IMF and the higher HN fraction.

The integration of the metallicity- and redshift-dependent top-heavy IMF with the metallicity-dependent HN fraction is the key innovation in CROCODILE, leading to a significantly higher $\zeta_{\rm SN}$ (by an order of magnitude) at low metallicities and high redshifts compared to previous simulations such as EAGLE and IllustrisTNG. In those simulations, $\zeta_{\rm SN}$ is also introduced as a decreasing function of metallicity to fine-tune the stellar mass function at $z=0$. However, in their cases, the SN energy is boosted only by a factor of three from the canonical value at low metallicity, based on the idea that the effective SN energy would be higher due to reduced radiative energy loss in a lower metallicity environment (Crain et al., 2015; Pillepich et al., 2018b). In contrast, our SN feedback model introduces a new consideration, resulting in a significant difference in one of the major feedback parameters. Additionally, we include the metallicity dependence on momentum feedback in our SN feedback model, as detailed later in Section 2.4.1.

We divide the Type-II SN feedback from a star particle into $n_{\rm event,snii}$ events so that the SN energy and metal yield are deposited gradually rather than instantaneously (Shimizu et al., 2019), and in this work, we adopt $n_{\rm event,snii}=2$. The SN feedback event occurs following the yield table considering the stellar age. The energy and metal output of SN is distributed to surrounding $N_{\rm ngb,fb}=8\pm 2$ gas particles using the feedback model developed in Oku et al. (2022) with some updates, as we briefly summarize in the following. The model consists of two components: local mechanical feedback and galactic wind feedback.

For Type-Ia SNe, the element yield table by Seitenzahl et al. (2013) is used in combination with a power-law type delay-time distribution function of $\Psi(t)=4\times 10^{-13}\,{\rm SN}\,{\rm yr}^{-1}\,{\rm M_{\odot}}^{-1}\,(t/1\,{\rm Gyr})^{-1}$ with a minimum delay of 40 Myr (Maoz & Mannucci, 2012). We divide the Type-Ia SN feedback from a star particle into $n_{\rm event,snia}=8$ events. The energy per Type-Ia SN is fixed to $10^{51}\,{\rm erg}$, and we use Eq. (3) to convert the energy to feedback momentum. We did not consider thermal galactic wind feedback for Type-Ia SNe because their degree of clustering is low, and they are not considered to be the major driver of the galactic wind.

For AGB stars, two yield tables are combined to cover the mass range of 1–8 ${\rm M_{\odot}}$; Karakas (2010) table is used in $M_{*}/M_{\odot}\in[1,6]$, Doherty et al. (2014) table in $M_{*}/M_{\odot}\in[7,8]$, and we linearly interpolate between them. We do not consider the mechanical feedback due to the stellar wind of AGB stars. We divide the AGB feedback from a star particle into $n_{\rm event,agb}=8$ events.

For the supermassive black holes (BHs), we largely follow the model used in the EAGLE project (Schaye et al., 2015; Crain et al., 2015). The model consists of seeding using the Friend-of-Friends (FoF) halo finder, torque-limited Bondi-Hoyle accretion, and thermal quasar feedback, as summarized below.

We perform simulations with varied feedback settings, black hole accretion parameters, and resolution as summarized in Table 2.

The Fiducial run considers the two-component SN feedback and AGN feedback. The NoSNGalWind run does not use the SN-driven hot galactic wind model. It is presented to show the effect of explicit modeling of the galactic wind in comparison to the Fiducial run. The AGNonly and SNonly runs only consider AGN and SN feedback, respectively, and are performed to investigate the impact of these feedbacks. In the NoZdepSN run, the stellar IMF is fixed to the Chabrier IMF, and the hypernova fraction is set to 0.01, independent of metallicity and redshift. This is different from other runs, where the specific energy of Type-II SN is higher for lower metallicity stars as shown in Fig. 1 due to metallicity- & redshift-dependent IMF and hypernova fraction. The NoFB run does not consider any feedback and is presented as a baseline to show the impact of feedback in comparison to other runs.

The LowCvisc and LowCviscLowMseed runs adopt black hole accretion parameters different from the Fiducial run. In these two runs, we adopt $C_{\rm visc}=2\pi$ (same value as the fiducial EAGLE simulation), which is 100 times smaller than the one adopted in the Fiducial run. The parameter $C_{\rm visc}$ modulates the limiting value of parameter $C$ in Eq. 6 through the accretion suppression factor of $C_{\rm visc}^{-1}(c_{s}/v_{\phi})^{3}$, which is usually less than unity. The lower $C_{\rm visc}$ results in a higher suppression factor (Rosas-Guevara et al., 2015; Schaye et al., 2015; Crain et al., 2015), which allows BHs to accrete at near Bondi rate, hence faster and earlier growth at $M_{*}\simeq 10^{9}-10^{10}\,{\rm M_{\odot}}$ as we will see in Section 3.2. The LowCviscLowMseed adopts the BH seed mass of $M_{\rm seed}=10^{4}\,{h^{-1}\,{{\rm M_{\odot}}}}$, 10 times smaller than that used in the other two runs.

The L25N512, L25N256, and L25N128 runs have different mass resolutions and are performed to investigate the dependence on resolution. The size of the simulation box of the L25N256 run is $25\,{h^{-1}\,{\rm cMpc}}$, and the number of employed particles is $2\times 256^{3}$, as the name indicates. The L25N256 run has the same resolution as the Fiducial run (i.e., L50N512). The L25N512 and L25N128 have eight times better and worse mass resolution than the L25N256 run. The same configurations and parameters are employed in the three runs except for the gravitational softening length. In the L25N256 run, we use the same gravitational softening length as the Fiducial run, and we set it to twice as small and as large in the L25N512 and L25N128 runs.

Simulation data is analyzed by on-the-fly friends-of-friends (FoF) and the subfind group and substructure finder (Springel et al., 2001, 2021). For visualization and further data analysis, we generate a uniform three-dimensional cartesian mesh covering the entire simulation box with $1024^{3}$ voxels, assign particle data, and take projections on the fly. The cloud-in-cell (CIC) assignment is used for dark matter and star particles, and the SPH kernel weight is used for gas particles. The SPH particle data is assigned to the mesh, conserving the total mass, metal mass, and internal energy. For only data analysis purposes, we set a floor in the SPH smoothing length at $(\sqrt{3}/2)L_{\rm voxel}$, where $L_{\rm voxel}$ is the size of a voxel, so that there are some voxel centers within the smoothing length from each SPH particle. We then assign the mass of $m_{i}W_{ij}(h_{i})/\sum_{k}^{N}W_{ik}(h_{i})$ to the $j$ -th voxel of the $i$-th particle, where $m_{i}$ and $h_{i}$ are the mass and smoothing length of the $i$-th particle, and $W_{ij}(h_{i})=W(r_{ij};h_{i})$ is the kernel weight, $r_{ij}$ being the distance from the $i$-th particle to the center of $j$-th voxel. The metal mass and internal energy are assigned to voxels in the same manner, and the metallicity and temperature of each voxel are calculated using the mass, metal mass, and internal energy of the voxel. This mesh generator module is useful for the post-process analysis of the IGM presented in Section 4.

Figure 2 shows the density projection image of the Fiducial run at $z=0$. The baryonic cosmic web, composed of intermediate-temperature ($T\sim 10^{4}\,{\rm K}$) filaments and high-temperature nodes, is apparent. The rest of the cosmic volume is filled with low-temperature voids. ¹⁰¹⁰10The movie is available at https://www.yurioku.com/research/

Figure 3 shows the galaxy stellar mass function (SMF) at $z=6$, 2.2, 1, and 0.1. The shaded regions are the observational constraints on SMF from Song et al. (2016); Stefanon et al. (2017, 2021) at $z=6$, McLeod et al. (2021) at $z=1$ and 2, Davidzon et al. (2017) at $z=2.25$, and Baldry et al. (2012) and Driver et al. (2022) around $z=0.1$.

The Fiducial, SNonly, and NoGalWind runs show suppressed star formation compared to the runs without SN mechanical feedback (AGNonly and NoFB) or the run with weaker SN feedback (NoZdepSN) at $M_{*}<10^{11}\,{\rm M_{\odot}}$, leading to a better consistency with the observation at $z=2.2$ and 1. This indicates that the SN mechanical feedback boosted by top-heavy IMF and hypernova is the dominant feedback source for the low-mass galaxies. We discuss the implication of the discrepancy at $z=6$ in Section 5.

At $z=0.1$, the Fiducial run nicely reproduces the low-mass end of the observed SMF. Both our simulations and the EAGLE simulation decrease towards a higher mass end with a slightly steeper slope than the observed exponential cutoff at $M_{*}\sim 10^{11.5}\,{\rm M_{\odot}}$, but capture the cutoff mass range relatively well.

In the Fiducial run, we used the metallicity- and redshift-dependent SN energy yield model, which has increasing specific SN energy towards lower SN progenitor metallicity as shown in Fig. 1. In the NoZdepSN run, the specific SN energy is fixed to $\zeta_{\rm SN}\simeq 6\times 10^{48}\,{\rm erg}\,{\rm M_{\odot}}^{-1}$. One can see that the NoZdepSN run completely overpredicts the number of low-mass galaxies at $z<2.2$. This indicates that the SN feedback with the above constant $\zeta_{\rm SN}$ is not strong enough to suppress star formation and reproduce the observed SMF.

Figure 4 shows the star-formation main sequence (SFMS) at $z=6$, 2.2, 1, and 0.1. The lines and error bars show the median SFR values and 16th and 84th percentile for each stellar mass bin in our simulations. The shaded regions are the observational constraints from Salmon et al. (2015) at $z=6$, Nakajima et al. (2023) at $z=$4-10, Karim et al. (2011) at $z=2.25$, 1.1, 0.9, and 0.3, Tomczak et al. (2016) at $z=2.25$, 1.125, and 0.875, and Whitaker et al. (2012) at $z=0-0.5$. We display only the star-forming galaxies, selected by a criterion sSFR $>10^{-2}\,{\rm Gyr}^{-1}$, because the observed data are also for star-forming galaxies.

At $z=6$, the simulations show similar SFMS with each other. The runs without stellar feedback or with weaker SN feedback show lower SFR than the observed SFMS for a given stellar mass at lower redshifts. This indicates that stellar feedback is necessary to sustain star formation activity by driving interstellar turbulence and galactic wind to delay the depletion of gas in the dark matter halo and produce star-forming galaxies at $z=0$. Runs with stellar feedback show a nice agreement with observations.

Figure 5 shows the mass–metallicity relations (MZR) at $z=6$, 2.2, 1, and 0.1. The metallicity of simulated galaxies is computed with the total metal mass in the gaseous phase. The lines and error bars show the median metallicity values and 16th and 84th percentile for each stellar mass bin in our simulations. The shaded regions are observational constraints from Nakajima et al. (2023) at $z=$4-10, Sanders et al. (2021) and Strom et al. (2022) at $z=2.3$, and Andrews & Martini (2013) and Tremonti et al. (2004) at $z=0.1$. The galaxies in the runs without stellar feedback have higher metallicity at $M_{*}\lesssim 10^{9}\,{\rm M_{\odot}}$, which indicates that stellar feedback expels metals from galaxies.

One can see that the NoSNGalWind run has slightly higher metallicity than the Fiducial run. This is because of the metal reduction by the galactic wind in the Fiducial run. The typical metal loading factor of our galactic wind model is 0.3, corresponding to 0.15 dex, which can explain the difference between the NoSNGalWind and Fiducial run. The difference among the models is small at $M_{*}\gtrsim 10^{10}\,{\rm M_{\odot}}$ because the feedback effect is weaker at the massive end, where the gravitational binding energy is stronger.

The metallicity of simulated galaxies of $M_{*}\sim 10^{10}\,{\rm M_{\odot}}$ is 0.2 dex higher than the MZR by Sanders et al. (2021) at $z=2.2$, and 0.4 dex higher than the MZR by Andrews & Martini (2013) at $z=0.1$. This can be due to either inefficient metal ejection by the galactic wind or insufficient gas mass fraction in the simulated galaxies at low redshift. A low gas fraction can be caused by strong feedback suppressing the accretion of low-metallicity gas or active star formation converting gas into stars. The metallicity increases as Type-Ia SNe and AGB stars slowly supply metals into ISM with certain time delays ($\sim$100 Myr) without being affected by instantaneous feedback by Type-II SNe. We will investigate the galactic outflow properties in our simulations by comparing them with recent IFU data at $z=0$-2, as well as the gas fraction and the metal abundance pattern in ISM in our future work.

The NoZdepSN run has lower metallicity than the NoFB run, while they agree in SMF and SFMS. This is because the top-heavy IMF is not adopted in the NoZdepSN run.

Figure 6 shows the BH mass as a function of the galaxy’s total stellar mass. The points and error bars show the median value and 16th and 84th percentile of the most massive BH mass in galaxies in each stellar mass bin. The region indicated by the black solid line is the observational constraint at $z=0$ compiled by Habouzit et al. (2021). The gray diamonds are observational data by McConnell & Ma (2013).

The runs with AGN feedback have lower BH masses at $z\lesssim 3$ because the gas around BHs is blown away by the AGN feedback. In the runs with AGN feedback, the $M_{\rm BH}$–$M_{*}$ relation is consistent with observations, and the medians are almost at the lower edge of the observationally constrained region at $10^{9}\,{\rm M_{\odot}}<M_{*}<10^{10.5}\,{\rm M_{\odot}}$. This behavior is sensitive to the choice of black hole subgrid parameters, as demonstrated in Section 3.2.

The AGNonly run also shows higher BH mass than the Fiducial run at $z>1$, indicating that the stellar feedback inhibits the gas accretion onto BHs and suppresses BH growth. At $z=0.1$, the $M_{\rm BH}$–$M_{*}$ relation of the two runs converges. This convergence is interesting considering that the AGNonly and Fiducial run have distinct SMF (Fig. 3; AGNonly run significantly overpredicts the SMF relative to the Fiducial run), suggesting that the growth of supermassive BHs is self-regulated by AGN feedback, and the BH and galaxy co-evolve even if the $M_{*}$ increases more rapidly due to lack of SN feedback.

Figure 7 again shows the $M_{\rm BH}$–$M_{*}$ relation, SMF, and SFMS, but for the runs with different BH parameters, Fiducial, LowCvisc, and LowCviscLowMseed, at $z=2.2$ and 0.1.

Panels A and B show the $M_{\rm BH}$–$M_{*}$ relation. The $C_{\rm visc}$ controls the onset of BH growth, and the BH mass starts to increase at $M_{*}\sim 10^{9}\,{\rm M_{\odot}}$ for runs with low $C_{\rm visc}$ ($=2\pi$) and at $10^{10}\,{\rm M_{\odot}}$ for runs with our default $C_{\rm visc}$ ($=200\pi$) at $z=0.1$. The simulation data points of the Fiducial run are around the bottom of the region indicated by Habouzit et al. (2021), while those of the LowCvisc run are around the top of the region. At $M_{*}>10^{11}\,{\rm M_{\odot}}$, it is interesting to see that the $M_{\rm BH}$–$M_{*}$ relation for all three simulations agree due to the self-regulated growth of BH.

Panels C and D show the SMF. By comparing with panels A and B, one can see that the SMF rapidly declines at the same mass scale where the BH mass rapidly increases. The SMF of the LowCvisc run is lower than that of other runs at $M_{*}>10^{9.5}\,{\rm M_{\odot}}$ due to the early growth of the BHs and subsequent AGN feedback. In the LowCviscLowMseed run, the BH growth is delayed due to the smaller BH seed mass, as seen in the upper panels, and it shows a similar SMF to the Fiducial run.

Panels E and F show the impact of AGN feedback on the star-formation activity in the SFMS. The slope of the SFMS becomes shallower at the same stellar mass where the BH begins to grow rapidly.

Figure 8 shows the SMFs of the runs with different resolutions, L25N512, L25N256, and L25N128. We used the same subgrid and numerical parameters except for the gravitational softening length in these runs. They roughly converge and agree with the observed SMF. The SMF slightly increases with increasing resolution at $z=2.2$. This is attributed to the density dependence of the star formation efficiency; higher-resolution simulations can resolve higher-density star-forming clouds, which have shorter star formation timescales.