GABE: Galaxy Assembly with Binary Evolution

TL;DR

GABE introduces binary star evolution into semi-analytic galaxy formation models, improving the accuracy of luminosity predictions and enabling studies of high-energy astrophysical phenomena within a cosmological framework.

Contribution

This work is the first to incorporate binary evolution into semi-analytic galaxy formation models using the Yunnan-II stellar population synthesis model.

Findings

Model slightly reduces optical/infrared luminosity in the local universe.

Model significantly increases ultraviolet luminosity, especially in the F_UV band.

Predicted galaxy colors are bluer, with minimal change in luminosity functions.

Abstract

We developed a new semi-analytic galaxy formation model: Galaxy Assembly with Binary Evolution (GABE). For the first time we introduce binary evolution into semi-analytic models of galaxy formation by using Yunnan-II stellar population synthesis model, which includes various binary interactions. When implementing our galaxy formation model onto the merger trees extracted from the Millennium simulation, it can reproduce a large body of observational results. We find that in the local universe the model including binary evolution reduces the luminosity at optical and infrared wavelengths slightly, while increases the luminosity at ultraviolet wavelength significantly, especially in $F_{\rm UV}$ band. The resulting luminosity function does not change very much over SDSS optical bands and infrared band, but the predicted colors are bluer, especially when $F_{\rm UV}$ band is under consideration. The new model allows us to explore the physics of various high energy events related to the remnants of binary stars, e.g. type Ia supernovae, short gamma-ray bursts and gravitational wave events, and their relation with host galaxies in a cosmological context.

Tables

Table 1. Table 1 : Free parameters in GABE and their fiducial values.

Parameter	Description	Fiducial Value
$\alpha$	Star formation efficiency, Equ.(9)	0.02
$ϵ$	Amplitude of SN reheating efficiency, Equ.(13)	6.5
${\beta }_1$	Slope of SN reheating efficiency, Equ.(13)	3.5
$\eta$	Amplitude of SN ejection efficiency, Equ.(16)	0.32
${\beta }_2$	Slope of SN ejection efficiency, Equ.(16)	3.5
$V_{\max ,c}$	Characteristic $V_{\max }$ in SN efficiency, Equ.(13) and Equ.(16)	90$\mathrm{km}{s}^{-1}$
$\gamma$	Outflow return efficiency, Equ.(17)	0.3
$f$	Black hole growth efficiency in mergers, Equ.(31)	0.03
$\kappa$	Quiescent Black hole accretion rate, Equ.(18)	$8.64\times {10}^{-4}$

Equations

$$f_b(M,z)=f_b^c{\left[1+\left(2^{\alpha /3}-1\right){\left(\frac{M}{M_c(z)}\right)}^{-\alpha }\right]}^{-\frac{3}{\alpha }},$$

$$t_{\mathrm{cool}}(r)=\frac{T}{\mathrm{d}T/\mathrm{d}t}=\frac{\frac{3}{2}n(r)kT}{\mathcal{L}(r)},$$

$$t_{\mathrm{dyn},\mathrm{h}}=\frac{R_{\mathrm{vir}}}{V_{\mathrm{vir}}}=0.1H(z{)}^{-1}.$$

$${\dot{M}}_{\mathrm{cool}}=\frac{M_{\mathrm{hot}}}{t_{\mathrm{dyn},\mathrm{h}}}.$$

$${\dot{M}}_{\mathrm{cool}}=\frac{M_{\mathrm{hot}}}{t_{\mathrm{dyn},\mathrm{h}}}\frac{R_{\mathrm{cool}}}{R_{\mathrm{vir}}}.$$

$$R_{\mathrm{cool}}={\left[\frac{X{M}_{\mathrm{hot}}}{4\pi m_{\mathrm{p}}{R}_{\mathrm{vir}}{n}_{\mathrm{H},\mathrm{crit}}(T_{\mathrm{hot}},Z_{\mathrm{hot}},z)}\right]}^{\frac{1}{2}},$$

$${\Sigma }_{\mathrm{crit}}(r)=12\left(\frac{V_{\max }}{200\mathrm{km}{\mathrm{s}}^{-1}}\right){\left(\frac{r}{10\mathrm{kpc}}\right)}^{-1}{\mathrm{M}}_{\odot }\mathrm{p}{\mathrm{c}}^{-2},$$

$$M_{\mathrm{crit}}=11.5\times 10^9\left(\frac{V_{\max }}{200\mathrm{km}{\mathrm{s}}^{-1}}\right)\left(\frac{R_{\mathrm{disk},\mathrm{gas},\mathrm{d}}}{10\mathrm{kpc}}\right){\mathrm{M}}_{\odot },$$

$${\dot{M}}_{\mathrm{SF},\mathrm{disk}}=\alpha \left(M_{\mathrm{disk},\mathrm{gas}}-M_{\mathrm{crit}}\right)/t_{\mathrm{dyn}},$$

$${\dot{M}}_{\mathrm{SF},\mathrm{sph}}=M_{\mathrm{sph},\mathrm{gas}}/t_{\mathrm{sf},\mathrm{sph}},$$

$$t_{\mathrm{sf},\mathrm{sph}}$$

$$\delta M_{\mathrm{reheat}}={ϵ}_{\mathrm{reheat}}\times \delta M_{\mathrm{SF}},$$

$${ϵ}_{\mathrm{reheat}}=ϵ\times \left[0.5+{\left(\frac{V_{\max }}{V_{\max ,\mathrm{c}}}\right)}^{-{\beta }_1}\right],$$

$$\delta E_{\mathrm{reheat}}=\frac{1}{2}\delta M_{\mathrm{reheat}}{V}_{\mathrm{vir}}^2\le \delta E_{\mathrm{SN}}=\frac{1}{2}\delta M_{\mathrm{SF}}{V}_{\mathrm{SN}}^2,$$

$$\delta M_{\mathrm{eject}}=\frac{{ϵ}_{\mathrm{eject}}\delta E_{\mathrm{SN}}-\delta E_{\mathrm{reheat}}}{\frac{1}{2}{V}_{\mathrm{vir}}^2},$$

$${ϵ}_{\mathrm{eject}}=\eta \times \left[0.5+{\left(\frac{V_{\max }}{V_{\max ,\mathrm{c}}}\right)}^{-{\beta }_2}\right],$$

$${\dot{M}}_{\mathrm{return}}=-\gamma \left(\frac{V_{\mathrm{vir}}}{220\mathrm{km}{\mathrm{s}}^{-1}}\right)\left(\frac{M_{\mathrm{out}}}{t_{\mathrm{dyn},\mathrm{h}}}\right),$$

$$\dot{M_{\mathrm{BH}}}$$

$$\dot{E_{\mathrm{radio}}}=0.1\dot{M_{\mathrm{BH}}}{c}^2,$$

$$\dot{M_{\mathrm{cool},\mathrm{eff}}}=\dot{M_{\mathrm{cool}}}-2\dot{E_{\mathrm{radio}}}/V_{\mathrm{vir}}^2.$$

$$M_{\mathrm{disk},\ast }>\frac{3R_{\mathrm{disk},\ast ,\mathrm{d}}{V}_{\max }^2}{G}\equiv M_{\mathrm{disk},\ast ,\mathrm{crit}}.$$

$$M_{\mathrm{bar}}$$

$$M_{\mathrm{hot}}$$

$$M_{\mathrm{hot}}\ge 0.$$

$$\delta M_{\mathrm{hot}}$$

$$\delta M_{\mathrm{disk},\mathrm{gas}}=M_{\mathrm{cool},\mathrm{eff}}-M_{\mathrm{SF},\mathrm{disk}}-M_{\mathrm{reheat},\mathrm{disk}},$$

$$\delta M_{\mathrm{disk},\ast }=M_{\mathrm{SF},\mathrm{disk}}-M_{\mathrm{bar}},$$

$$\delta M_{\mathrm{sph},\mathrm{gas}}=-M_{\mathrm{SF},\mathrm{sph}}-M_{\mathrm{reheat},\mathrm{sph}},$$

$$\delta M_{\mathrm{sph},\ast }=M_{\mathrm{SF},\mathrm{sph}}+M_{\mathrm{bar}},$$

$$\delta M_{\mathrm{out}}=M_{\mathrm{eject}}-M_{\mathrm{return}},$$

Full text

1 Introduction

In the framework of $\Lambda$CDM cosmology, semi-analytic models of galaxy formation have been very effective to simulate the formation and evolution of galaxies, and have played an indispensable role during the building of current galaxy formation theory from early 1990s to 2010s (White & Frenk 1991; Kauffmann et al. 1999; Springel et al. 2001; Croton et al. 2006; De Lucia & Blaizot 2007; Guo et al. 2011; Henriques et al. 2015). This method is very successful to explore the effect of various physical processes on galaxy formation. However, some physical models still need further investigation. Binary evolution is among one of them.

It is well known that more than $\sim 50\%$ of field stars and $\sim 70\%$ stars in massive young star clusters are in binary systems, with only a relatively small part of stars being single (e.g. Duquennoy & Mayor 1991, Sana et al. 2012). Binary stars undergo many different physical processes from single stars, such as mass transfer, mass accretion, common-envelope evolution, collisions, supernova kicks, tidal evolution, and angular momentum loss. Such binary interactions change the color and Spectral Energy Distribution (SED) of the entire stellar population (Pols & Marinus 1994; Han et al. 2002). Though binary evolution could influence photometric properties of stellar populations, which will certainly further influence the properties of galaxies, most existing semi-analytic models adopt classical single stellar population synthesis models (SPSMs), such as BC03 (Bruzual & Charlot 2003), CB07 (Bruzual 2007) and M05 (Maraston 2005).

Another motivation to use the SPSM with binary evolution is that many high energy astrophysical events are consequences of binary evolution. For example, the short gamma-ray bursts are thought to be products of double neutron star (NS-NS) or neutron star-black hole (NS-BH) mergers (Narayan et al. 1992); the supernova Ia events are due to white dwarf and non-degenerate companion evolution (Whelan & Iben 1973) or double white dwarf (WD-WD) mergers (Iben & Tutukov 1984); and the gravitational wave events are supposed to be triggered by the merger of two compact objects in the binary systems. Furthermore, the fast radio burst events are possibly connected with binary neutron star mergers (Totani 2013; Wang et al. 2016). With the rapidly increasing number of detection of such high energy events at a wide range of redshift, a large sample of such events may be investigated in a cosmological context in near future. It is timely to adopt a SPSM with binary evolution in galaxy formation model in order to explore the connection between such events, and their relationship with their host galaxies, over the history of the Universe.

In this paper, we developed a new semi-analytic model GABE (Galaxy Assembly with Binary Evolution). The model succeeds many physical recipes from existing successful semi-analytic models and includes some new ingredients, for example, adopting a SPSM with binary evolution–Yunnan-II model and a cooling table of Wiersma09 (Wiersma et al. 2009) which treats the effect of photoionization on gas cooling quite carefully. In GABE, we explicitly track the evolution of binary systems in each simple stellar population111A simple stellar population represents a set of stars which have the same age and metallicity. “Simple” is used to distinguish from the so-called complex stellar population, which is composed of multiple simple stellar populations. produced in the formation of each galaxy in the cosmic background. By recording the evolution of each binary event, especially its final fate as remnants of merging binaries, which are supposed to be associated with the high energy events mentioned above, our model could predict every such high energy event in each galaxy across cosmic time.

The structure of this paper is as follows. In Section 2, we introduce the N-body simulation we use in this work, and describe our galaxy formation model in detail. The results of our semi-analytic models are demonstrated in Section 3. In the end, we summarize our results in Section 4.

2 Galaxy Formation Models

In this section, we describe the N-body simulation and all physical models used in GABE. The physical driving models for galaxy formation have been developing gradually in the past decades (White & Frenk 1991; Kauffmann et al. 1999; Springel et al. 2001; Croton et al. 2006; De Lucia & Blaizot 2007; Guo et al. 2011; Henriques et al. 2015). In our model, we follow these successes and consider a full set of known physical models for galaxy formation. Apart from these semi-analytic models, another method of modeling galaxies’ formation and evolution is cosmological hydrodynamic simulation. Hydrodynamic simulations have the advantage to be able to study the physical processes in more detail than the semi-analytic method, while the latter requires less computational cost and therefore is more flexible in tuning and testing model parameters. Cosmological hydrodynamic simulations have achieved huge successes in recovering galaxies’ statistical properties these years (e.g. Schaye et al. 2015; Nelson et al. 2018), and there are some works which have successfully planted binary evolution into existing hydrodynamic simulations by post-processes (e.g. Mapelli et al. 2017). Whereas, there is no binary evolution in semi-analytic models yet. In this work, by combing binary SPSMs with semi-analytic models, the calculation of binary properties of galaxies can be done on-the-fly and will not be limited by the time resolution of snapshots, which is more suitable for generating mock galaxy catalogue with binaries or double compact objects for larger volume.

In Section 2.1, the N-body simulation used in GABE is introduced. The rest subsections describe in detail our galactic formation models. Note that most physical models described in this work are the same with Guo et al. (2011), which is a well-accepted semi-analytic model. Compared with Guo et al. (2011), our modifications focus on hot gas cooling (Section 2.2.2) and stellar population synthesis models (Section 2.4.1). In Section 2.2, we describe all physical processes involved in the evolution of an isolated galaxy, including reionization, gas cooling, star formation, supernova feedback, black hole growth, AGN feedback, and bar formation. A subsection on exchanges of mass, metal and angular momentum is followed as a short summary. Then we describe how to treat galaxy mergers in our model in Section 2.3, and how to calculate galaxies’ luminosity in Section 2.4. At last, we discuss our model calibration in Section 2.5.

2.1 N-body Simulation

In this study, we take advantage of the dark mater merger trees from the Millennium-I Simulation (Springel et al. 2005), which is widely used for studies of galaxy formation and evolution because of its good combination of large box size and relatively high mass resolution. Its box size is $L_{box}=685\,{\rm Mpc}$ and mass of simulation particle is $1.18\times 10^{9}\,{\rm M}_{\odot}$. Accordingly, its halo resolution is $2.36\times 10^{10}\,{\rm M}_{\odot}$, which is 20 times of the particle mass. These features make sure the following semi-analytic models can generate complete galaxy catalogue for galaxies more massive than $\sim 10^{8}\,{\rm M}_{\odot}$ and suppress cosmic variances at the same time. Millennium-I simulation assumes a WMAP1 cosmology with $\Omega_{\rm m}=0.25$, $\Omega_{\rm b}=0.045$, $\Omega_{\rm\Lambda}=0.75$, $n=1$, $\sigma_{\rm 8}=0.9$ and $H_{\rm 0}=73\,\,{\rm km}\,\,{\rm s}^{-1}\,{\rm Mpc^{-1}}$, which is derived from a combined analysis of the 2dFGRS (Colless & et al. 2001) and the first-year WMAP data (Spergel et al. 2003). The merger trees are built based on D-halo catalog (Jiang et al. 2014).

2.2 Evolution of Isolated Galaxies

2.2.1 Cosmic Reionization

It is now well accepted that photoheating by UV background significantly lowers gas content in dwarf sized dark matter halos. Here we use a formula put forward by Gnedin (2000):

[TABLE]

where $f_{b}^{c}$ is mean cosmic baryon fraction, which is $17\%$ in WMAP1 cosmology, and $M$ is the total halo mass, including dark matter and baryon. There are two parameters in Equ.(1): $\alpha$ and $M_{c}(z)$. With larger $\alpha$, the baryonic fraction decreases more rapidly with decreasing halo mass. $M_{c}(z)$ is the characteristic mass below which halos only contain less than 50 percent of universal baryon fraction. Here we adopt the latest results from Okamoto et al. (2008): $\alpha=2$; $M_{c}\approx 10^{7}\,h^{-1}\,{\rm M}_{\odot}$ just after reionization and $M_{c}\approx 6.49\times 10^{9}\,h^{-1}\,{\rm M}_{\odot}$ at $z=0$, as shown in the Fig. 3 of Okamoto et al. (2008).

2.2.2 Gas Cooling with Background Radiation

When gas falls into the potential well of dark matter halo, it will be shock-heated. These hot and dense gas will then cool down by radiating away the energy. If the cooling rate of the gas is fast enough, gas can cool and condense to the center of the halo at the free-fall rate. Otherwise, a shock-front will form and heat the gas to about the virial temperature of the halo. As the gas has to cool down before it can settle down to the central galaxy through cooling flows, the cooling rate of gas determines the amount of gas that can settle to the center.

The local cooling timescale of the gas is usually estimated as

[TABLE]

where $n(r)$ is the number density of the gas at radius $r$, and can be estimated by $n(r)=\rho(r)/\mu m_{\rm p}$; $\rho(r)$ is the local density of the hot gas, which is assumed to be an isothermal profile in our model: $\rho(r)\propto r^{-2}$; $\mu$ is the mean molecular mass, and $m_{\rm p}$ is the proton mass; $k$ is the Boltzmann constant; $T$ is the temperature of gas; $\mathcal{L}$ is the cooling rate per unit volume of gas with a unit of $\,{\rm erg}\,{\rm s^{-1}}\,{\rm cm^{-3}}$.

“Cooling radius”, $R_{\rm cool}$, is quite commonly used in many semi-analytic models. For a given dark matter halo, we follow Springel et al. (2001) to define cooling radius as the radius within which the gas can cool in a dynamical timescale, i.e. $t_{\rm cool}(R_{\rm cool})=t_{\rm dyn,h}$, and

[TABLE]

If $R_{\rm cool}>R_{\rm vir}$, the halo is in the rapid infall regime. Gas condenses to the halo center at the rate of free fall:

[TABLE]

If $R_{\rm cool}<R_{\rm vir}$, shock front forms at $R_{\rm cool}$ and only the gas inside $R_{\rm cool}$ can condense to the center. As the hot gas profile in halo is assumed to be isothermal, we have

[TABLE]

The cooling rate $\mathcal{L}$ adopted in many existing semi-analytic galaxy formation models assumes the collisional ionization equilibrium (CIE) assumption, which means the ions and electrons of gas are in collisional equilibrium, and background radiation is ignored. With CIE assumption, the cooling rate of low density plasma is simply proportional to the gas density square: $\mathcal{L}=n^{2}\Lambda(T_{\rm hot},Z_{\rm hot})$. $\Lambda$ is the so-called “cooling function” and has the unit of $\,{\rm erg}\,{\rm s^{-1}}\,{\rm cm^{3}}$. $\Lambda$ only depends on the temperature and metallicity of the gas and can be conveniently calculated and tabulated (e.g. Sutherland & Dopita 1993). However, generally CGM is not in CIE state. Efstathiou (1992) found that a significant UV background radiation is already in place while galaxies form and can lengthen the cooling timescale of gas. Wiersma et al. (2009) (hereafter Wiersma09222The cooling tables are offered on their website: http://www.strw.leidenuniv.nl/WSS08/.) provided a cooling table by taking into account UV background radiation quite carefully. We adopt this table in our model.

When considering UV radiation background, the cooling rate $\mathcal{L}$ is no longer simply proportional to gas density square. However, we can still use density square to normalize the table for easier usage. In Wiersma09, the relation is $\mathcal{L}=n_{\rm H}^{2}\Lambda(T,Z,n_{\rm He}/n_{\rm H},n_{\rm H},z)$, where $T$ and $Z$ are the temperature and metallicity of gas respectively, $n_{\rm H}$ and $n_{\rm He}$ are the number densities of hydrogen and helium respectively, $z$ is the redshift.

The cooling table of Wiersma09 can not be easily used because of two reasons: (1) The cooling function depends on local density $n_{\rm H}$ which varies with radius. Therefore the equation $t_{\rm cool}(r)=t_{\rm dyn,h}$ can only be solved numerically. (2) We use the Equ.(5) of Wiersma et al. (2009) to calculate the cooling function, therefore we have to look up eight tables to get one cooling rate. This is extremely time consuming. To overcome these, we simplify the calculation as the following two steps:

(1) We decrease the radius dependence by assuming the gas has a uniform temperature, metallicity and composition distribution along radius, i.e. we estimate gas temperature as the virial temperature of host dark matter halo, $T=T_{\rm vir}=\frac{1}{2}\frac{\mu m_{\rm H}}{k}V_{\rm vir}^{2}=35.9(V_{\rm vir}/{\,\,{\rm km}\,\,{\rm s}^{-1}})^{2}{\rm K}$, and estimate the mean molecular mass and composition of gas as the fully ionized primordial gas, $\mu=0.59$ and $n_{\rm He}/n_{\rm H}=0.083$.

(2) With the assumptions above, the left side of the equation $t_{\rm cool}=t_{\rm dyn,h}$ is a function of $(T,Z,n_{\rm H},z)$, and its right side is a function of redshift $z$ (as $t_{\rm dyn,h}=0.1H(z)^{-1}$). We can see that this relation is not halo dependent, which means we can solve this equation before the running of semi-analytic model and save a lot of CPU time. Thus we re-assemble tables of Wiersma09 to a new one, $n_{\rm H,crit}(T,Z,z)$. Here $n_{\rm H,crit}$ is the gas density where $t_{\rm cool}(n_{\rm H,crit})=t_{\rm dyn,h}$. To determine $n_{\rm H,crit}$, we calculate $t_{\rm cool}$ at each density bins of Wiersma09, from high density to low density, compare $t_{\rm cool}$ with $t_{\rm dyn,h}$, then record the first highest density when $t_{\rm cool}\geq t_{\rm dyn,h}$ as $n_{\rm H,crit}$. This table is pre-made and suitable for different situations in the models.

Therefore we only need to look up this table for one time to get the $n_{\rm H,crit}$, then we can transfer the critical density into cooling radius easily by

[TABLE]

where $X=0.75$ is the mass fraction of hydrogen of primordial gas.

2.2.3 Star Formation

As gas condenses to the center of a dark matter halo, it will form a rotation supported cold gas disk because of angular momentum conservation. If the local density of gas exceeds a critical value, stars will eventually form. The critical surface density at radius r is suggested by Kauffmann (1996):

[TABLE]

where $V_{\rm max}$ is the maximum circular velocity of the dark matter halo. Integrating $\Sigma_{\rm crit}(r)$ from 0 to three times scale radius of gaseous disk, as suggested by Croton et al. (2006), the critical gas disk mass reads

[TABLE]

where $R_{\rm disk,gas,d}$ is the exponential scale radius of cold gas disk. Therefore for disks with $M_{\rm disk,gas}>M_{\rm crit}$, a certain amount of the gas will be converted into stars:

[TABLE]

where $t_{\rm dyn}=3R_{\rm disk,gas,d}/V_{\rm max}$ is the characteristic timescale at the edge of the star-forming disk, and $\alpha$ is the “star formation efficiency”, a free parameter. Here a fiducial value of 0.02 is adopted.

For cold gas in spheroid, we use a star-formation timescale scheme from Benson (2011) to calculate its star formation rate:

[TABLE]

[TABLE]

where $t_{\rm sf,sph}$ is the star-formation timescale in the spheroidal component, $V_{\rm sph}$ is the characteristic velocity of the spheroid, and $R_{\rm sph,gas}$ is the half mass radius of the gas in the spheroid.

2.2.4 Supernova Feedback

Once stars form, they are assumed to follow the evolutionary tracks of the standard stellar evolution model, become more and more metal rich and eventually release large amount of mass and energy into their surroundings. We adopt the instantaneous stellar evolution assumption, which assumes stars evolve to their final states (e.g. $t_{\rm age}\sim 10\,{\rm Gyr}$) right after their birth. Compared to the relatively large timestep used in semi-analytic models, such a simplification is acceptable. Non-instantaneous scheme is also explored by other authors, e.g. De Lucia et al. (2014). During this instantaneous evolution, we assume $43\%$ of stellar mass is recycled back to circumstance, and the mass of created metal is $3\%$ of the initial stellar mass, in agreement with the Chabrier (2003) initial mass function (IMF) we use.

As stars evolve into their final state, they may release a large amount of energy into interstellar medium (ISM) instantaneously through supernova explosion. During the process, the ISM temperature may increase and suppress further star formation. If the ejected energy is large enough, some gas may be ejected out of the galaxy. In our model, the energy of supernovae is firstly used to reheat cold gas back to the hot gas phase; if the temperature of the reheated hot gas is higher than the virial temperature of dark matter halo, the exceeded energy will be used to put gas into an ejecta component, called “outflow”, which will not be able to cool back within a certain timescale.

The amount of reheated cold gas is estimated as

[TABLE]

[TABLE]

where $\epsilon$ and $\beta_{1}$ are free parameters describing the heating efficiency of newly formed stars, and $V_{\rm max,c}$ is a characteristic velocity used to adjust the shape of $\epsilon_{\rm reheat}$. We also make sure the energy of reheated gas will not be larger than the energy of supernova:

[TABLE]

where $\delta E_{\rm reheat}$ is energy of the reheated gas, $\delta E_{\rm SN}$ is the available supernovae energy, $V_{\rm SN}=630{\rm km/s}$ is the characteristic wind speed of supernovae. If $\delta E_{\rm SN}>\delta E_{\rm reheat}$, we use a certain mount of exceeded energy to heat more hot gas to virial temperature and put them in an ejecta component:

[TABLE]

where

[TABLE]

describes the ejection efficiency. We require $\epsilon_{\rm eject}\leq 1$. Here $\eta$, $\beta_{2}$ and $V_{\rm max,c}$ are free parameters.

The ejected gas may return to the hot gas halo due to gravitational attraction in a dynamical timescale, and be available again for cooling and star formation. The return rate decreases with halo mass. By following Guo et al. (2011), we have

[TABLE]

here $M_{\rm out}$ is the total mass of the outflow, $\gamma$ is a free parameter and is set to be $0.3$ in our fiducial model. In the above scheme, supernovae feedback is more efficient to heat the cold gas and suppress star formation in low-mass galaxies.

2.2.5 Black Hole Growth and AGN Feedback

Since the first discovery of a black hole in M32 (Tonry 1984), 87 more black holes have been further discovered in the center of galaxies till the end of 2012 (Kormendy & Ho 2013). More interestingly, there is a tight scaling relation between black hole mass and stellar mass or velocity dispersion of the bulge of the host galaxy (e.g. Magorrian et al. 1998; Marconi & Hunt 2003; Gebhardt et al. 2000), indicating the growths of central black hole and the bulge are regulated by each other. A simple way to build this scaling relation is the mass averaging in galaxy mergers (Peng 2007; Jahnke & Macciò 2011). There naturally arises a scaling relation between them only if they both grow mainly through galaxy mergers. Here we follow the scheme of Croton et al. (2006) to describe the growth of black hole, in which the black hole growth through both mergers and quiet accretion are considered.

In our model, when a galaxy forms, a black hole seed is initialized in the center of the galaxy. The seed may be the remnant of first stars, with the mass of a few hundred solar masses. The initial mass of the seed normally is only a tiny fraction of its later mass because of the rapid growth, thus changing the initial seed mass almost has no effect on the evolution of galaxy. In this model we set the seed mass to be zero for simplicity.

Black hole growth in Croton et al. (2006) is divided into two different modes: “quasar” mode and “radio” mode. Quasar mode accounts for black hole growth during mergers, which will be described later. Radio mode accounts for black hole growth through accreting hot gas:

[TABLE]

where $\kappa$ is a free parameter to describe the strength of accretion. When baryons fall into black holes, $10\%$ of their rest mass energy is released to circumstance medium, heating surrounding gas and suppressing gas cooling:

[TABLE]

[TABLE]

$\dot{M_{\rm cool,eff}}$ is the effective cooling rate when considering suppression of AGN feedback.

As Equ.(18) shows, the strength of accretion depends strongly on the virial velocity of halos, $V_{\rm vir}$, and it is expected that the AGN feedback is more effective in more massive galaxies.

2.2.6 Bar Formation

Not all stellar disks are stable. Toomre (1964) pointed out that only stellar disks with large radial velocity dispersion can suppress all axisymmetric instabilities, i.e. the equivalent of pressure resulting from motion must be large enough to overcome self-gravity and other instabilities. If the surface density of a galaxy stellar disk is high and the self-gravity overwhelms the pressure, it could be unstable and may transform to a bar (Christodoulou et al. 1995; Mo et al. 1998). Efstathiou et al. (1982) used N-body simulations to explore disk instabilities, and derived a criterion to judge whether a exponential disk is stable. Here we use a modified form of De Lucia & Blaizot (2007): the cold gas disk is unstable if

[TABLE]

When a transient bar forms, it will transfer stellar mass from inner stellar disk to the bulge. After the transfer, the disk will be stable again, so we estimate the transferred mass as $M_{\rm bar}=M_{\rm disk,*}-M_{\rm disk,*,crit}$.

We follow Guo et al. (2011) to determine the radius of transferred part, $R_{\rm bar}$:

[TABLE]

where $\Sigma_{\rm disk,*,0}$ is the central surface density of the stellar disk, $R_{\rm bar}$ is the estimated radius of the bar and can be solved numerically. The mass of the transferred part is added into the existing bulge and the size of new bulge is calculated in the same way as a merger occurred, which will be described later.

2.2.7 Exchange of Mass, Metal and Angular Momentum

When a galaxy evolves in isolation, all physical models should obey conservation laws, including mass and angular momentum conservation. Fig. 1 present a cartoon plot illustrating the complex exchange of baryonic matter in our model. Below we summarize in detail the migration of mass, metal, and angular momentum between different phases in our model.

In order to make sure the total mass of galaxy meets a certain fraction $f_{\rm b}$ of host halo mass, we adjust hot gas mass as

[TABLE]

and

[TABLE]

Therefore, the exchange of mass can be described by the following equations:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The ISM of galaxies is enriched by metals released by stellar winds and supernova explosions. Metal enrichment is firstly included into semi-analytic models in White & Frenk (1991), in which metals are created in stellar evolution and then distributed into other components with mass transfer. The detailed metal exchange in our models is summarized as follows:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $y_{\rm z}=0.03$ is the so-called “yield”, controlling the amount of metals created during stellar evolution. In our model, for simplicity, we assume the composition of metal to be solar value, and do not consider different elements independently. Fig. 2 shows the metallicity in stars and gas as a function of the stellar mass at $z=0$ in our fiducial model. Our model agrees well with observational metallicity for stars. As on gaseous metallicity, our model is marginally in agreement with Blanc et al. (2019) for galaxies more massive than $10^{9.5}\,{\rm M}_{\odot}$ and is marginally in agreement with the lower limit of Lee et al. (2006) for galaxies smaller than $10^{9.5}\,{\rm M}_{\odot}$, but way below the observational result of Tremonti et al. (2004) for galaxies less massive than $10^{11}\,{\rm M}_{\odot}$. This discrepancy between different observational results is caused by systematic uncertainties in the measurement of gaseous metallicity. Kewley & Ellison (2008) demonstrated that depending on the choice of metallicity calibration, the absolute metallicity normalization varies up to 0.7 dex. López-Sánchez et al. (2012) and Blanc et al. (2019) also showed the discrepancy between different diagnostics, and which one is preferred has not come to a conclusion in the community. As illustrated in the above equations, the amount of gaseous metallicity directly relates to galactic cooling, star formation, supernova feedback and AGN feedback. However, we choose not to use gaseous metallicity as a constrain or calibration for semi-analytic models before we have clearer understanding of the systematic uncertainties in observations.

We follow Guo et al. (2011) and also consider the angular momentum transfer among different components. Generally it is assumed that the hot gas share the same spin with its host dark matter halo, and the angular momentum of cooled gas will be transferred into galaxy when stars begin to form. The full exchange of angular momentum could be summarized as the following:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The disk density profile is often assumed to be exponential, then its scale radius is determined by disk’s angular moment. By assuming a flat circular velocity curve, disk scale radius of a galaxy can be derived with

[TABLE]

[TABLE]

2.3 Hierarchical Growth

In the above, we focus on physical processes of the galaxy formation in isolation. In the standard cold dark matter cosmology, structure formation is hierarchical, namely small systems form firstly and then merge to form larger and larger systems. In this section, we will focus on physical processes relevant to mergers.

2.3.1 Central and Satellite Galaxies

When two halos approach closer and merge to a single halo in the end, the one with larger mass is called the primary halo, and the smaller one is called substructure of the primary halo and will be eventually tidal disrupted. Following Springel et al. (2001), the galaxy in the primary halo is called central galaxy (“type 0” galaxy), while the galaxy in the substructure is called satellite galaxy. If a satellite galaxy has host subhalo, it is called “type 1” galaxy; and if not, it is called “type 2” galaxy. Central galaxies stay at center of the primary halo, while satellite galaxies orbit around their main halo and suffer from tidal stripping as well as ram-pressure stripping.

2.3.2 Tidal and Ram-Pressure Stripping of Hot Gas

Tidal radius of a galaxy is a scale beyond which all dark matter and gas get tidally stripped, which can be estimated as

[TABLE]

where $M_{\rm vir,infall}$ and $R_{\rm vir,infall}$ are the virial mass and radius of the satellite’s host halo at the infall time333We define the moment at which the satellite last passes the virial radius of its host central halo as the “infall” time., respectively.

In addition to the tidal stripping, the hot gas component in satellites suffer ram-pressure stripping due to gas pressure from the primary halo. While ram-pressure stripping acts on both hot and cold gas of satellites, in our model, for simplicity, we only implement this process to more extended hot gas component. The radius, $R_{\rm rp}$, that the hot gas of a satellite is stripped by ram-pressure can be estimated by the balance of the pressure from the primary halo and the binding energy of the satellite:

[TABLE]

here $\rho_{\rm sat}\left(R_{\rm rp}\right)$ is gas density of the satellite at $R_{\rm rp}$, $V_{\rm vir,sat}$ is the virial velocity of the satellite, $\rho_{\rm cen}\left(R_{\rm orbit}\right)$ is the gas density of the central at the position of the satellite, and $V_{\rm orbit}$ is the relative velocity between the satellite and the central. In our model we adopt the minimum of [$R_{\rm tidal}$,$R_{\rm rp}$] to determine the radius beyond which hot gas of satellite is stripped:

[TABLE]

2.3.3 Dynamical Friction

When satellites orbit their central galaxies, they are dragged by the surrounding matter due to gravity, gradually lose energy and angular momentum, and will eventually merge with the central galaxy. This process is called dynamical friction (Chandrasekhar 1943). As our model is based on subhalo merging trees of N-body simulations, whether a subhalo survives in merger may depend on employed numeric resolution. As assumed by many models that when a subhalo loses its identify in a simulation, namely below 20 particles, it is assumed that the stellar component of the subhalo (i.e. a type 2 galaxy) will still survive for a while because the stellar component is more compact than dark matter. Here we adopt the dynamical friction timescale fitting formula of Jiang et al. (2008) (hereafter Jiang08):

[TABLE]

where $\epsilon$ is the circularity of the satellite orbit, $C=0.43$ is a constant, $M_{\rm cen}$ is the primary dark matter halo mass, $M_{\rm sat}$ is the satellite mass, including dark matter halo and stars, and $V_{\rm c}$ is the circular velocity. $R_{\rm vir,cen}/V_{\rm c}$ represents free fall timescale in Jiang08, thus we use the virial velocity of central galaxy $V_{\rm vir,cen}$ to estimate $V_{\rm c}$. The circularity $\epsilon$ is calculated through eccentricity $e$ (Wetzel 2011):

[TABLE]

[TABLE]

where $E_{\rm orbit}$ and $J_{\rm orbit}$ are the total orbital energy and angular momentum of satellite-central system, $\mu=M_{\rm cen}M_{\rm sat}/(M_{\rm cen}+M_{\rm sat})$ is the reduced mass.

The dynamical friction timescale is calculated when the satellite first crosses the virial radius of its primary halo. After $t_{\rm fric}$, this satellite will merge into the central galaxy. This definition of $t_{\rm fric}$ is consistent with Jiang08.

2.3.4 Remnants of Mergers

Numerical simulations demonstrate that the remnants of mergers depend on the orbits and mass ratio of the two merging galaxies (Kannan et al. 2015). In our model, for simplicity, we follow White & Frenk (1991) to model the remnants of mergers based on the mass ratio of merging galaxies. If the mass ratio $f_{\rm merger}=M_{\rm sat}/M_{\rm cen}$ is greater than $0.3$, it is called “major merger”, otherwise it is called ”minor merger”. For a major merger, all disk components of both satellite and central are destroyed to form spheroidal components, including cold gas and stars. For a minor merger, the smaller satellite will be destroyed completely in the end. Its gas component is added into the gaseous disk of the central, and the stellar component is added into the stellar bulge of the central. In both cases, the hot gas and outflow component will be simply added together respectively. The metals and angular momentum of each component will also be transferred along with mass.

2.3.5 Black Hole Growth in Mergers

A very efficient mechanism for the growth of black hole is through mergers, especially major mergers. The violent dynamical environment during mergers is advantaged for black holes to accrete gas. Following Kauffmann & Haehnelt (2000) and Croton et al. (2006), the growth of central black hole during mergers can be estimated as

[TABLE]

where $M_{\rm BH,sat}$ is the black hole mass of the satellite; $M_{\rm sat}$ and $M_{\rm cen}$ are the masses (cold gas and stars) of the satellite and the central galaxy, respectively; $M_{\rm cold}$ is the total mass of cold gas of two merging galaxies; and $V_{\rm vir}$ is the virial velocity of the primary halo. $f$ is a free parameter controlling the growth rate of central black hole. In our fiducial model, we use $f=0.03$ and can derive the tight relation between black hole mass $M_{\rm BH}$ and bulge mass $M_{\rm BH}$ very well.

2.3.6 Star Burst in Mergers

Violent interactions during mergers can compress gas and trigger star bursts. We adopt the scheme “collisional starburst” of Somerville et al. (2001). In the model, the fraction of the gas transferred into stars during a merger is

[TABLE]

where $M_{\rm cen}$ and $M_{\rm sat}$ are the total mass (cold gas and stars) of the central and satellite galaxy, respectively. Along with the burst, supernova feedback releases huge amount of energy. In a major merger, almost all the remaining cold gas is heated into the hot halo, strongly suppressing further star formation.

2.3.7 Bulge Formation

Except material from bar transfer, galaxy merger is another important channel for bulge formation. We use the energy conservation scheme of Guo et al. (2011) to estimate the radius of newly formed spheroidal remnant:

[TABLE]

where $C=0.5$ is a structure parameter describing the binding energy of each system; $\alpha_{\rm inter}=0.5$ is a parameter describing the interactive energy between systems; $M_{\rm 1}$, $M_{\rm 2}$, and $M_{\rm new,bulge}$ are the masses of progenitor systems and newly formed bulge, respectively; $R_{\rm 1}$, $R_{\rm 2}$, and $R_{\rm new,bulge}$ are the corresponding half-mass radii of progenitor systems and newly formed bulge, respectively. Equ.(33) is also used to estimate the radius of remnant after bar formation.

2.4 Luminosity of Galaxies

In order to compare our model with observations, we need to calculate the photometric properties of our model galaxies with a stellar population synthesis model:

[TABLE]

where $L_{\rm\lambda}$ is the luminosity of a galaxy at wavelength $\lambda$; $L_{\rm\lambda,i}(Z,t)$ is the luminosity of unit solar mass stellar population with metallicity $Z$ and age $t$ at wavelength $\lambda$, which is provided by SPSMs; $N_{\rm SP}$ is the total number of stellar populations of the galaxy, $M_{\rm i}$, $Z_{\rm i}$ and $t_{\rm form,i}$ is the initial mass, metallicity and formation time of the $i$th stellar population, all of which are recorded during the running of GABE; $A_{\rm\lambda,ISM}$ and $A_{\rm\lambda,BC}$ are the reddening due to ISM and birth clouds respectively, which are provided by dust extinction models.

2.4.1 Yunnan-II Model: Stellar Population Synthesis Model with Binary Evolution

The SPSMs adopted in existing semi-analytic galaxy formation models are mostly based on single star stellar evolution models, e.g. BC03, BC07 with TP-AGBs, and M05. However more than $\sim 50\%$ field stars in our Milky Way are in binary systems. Interactions between stars could produce some specific objects, such as blue stragglers (Pols & Marinus 1994) and subdwarf B stars (sdBs, Han et al. 2002), both of which would alter the Spectral Energy Distribution (SED) of the entire stellar population significantly. In our model, we introduce binary interactions into galaxy formation and evolution by using the Yunnan-II stellar population synthesis model (Zhang et al. 2004, 2005, 2010)444http://www1.ynao.ac.cn/$\sim$zhangfh/.

Yunnan SPSMs are developed by the Group of Binary Population Synthesis, Yunnan Observatories, Chinese Academy of Sciences. The main differences among different versions of Yunnan models include the stellar evolution models, the evolutionary population synthesis construction method and the condition of binary interactions. The Yunnan-I model (Zhang et al. 2002) was constructed for stellar populations without binary interactions by using the rapid single stellar evolution algorithms (Hurley et al. 2000) and used the “isochrone synthesis” technique. The Yunnan-II model was constructed for both stellar populations with and without binary interactions by using Monte Carlo simulation and the rapid single and binary stellar evolution algorithms (Hurley et al. 2000, 2002). In their latest version of Yunnan-II model (Zhang et al. 2010), the evolutionary population synthesis models of Han et al. (2007) which included sdBs were also combined. The Yunnan-III model (Zhang et al. 2013) was constructed by using single stars’ evolutionary tracks obtained from the detailed MESA stellar evolution models (Paxton et al. 2011). Thus Yunnan-II model has a binary and a single version. The Yunnan-II model with full binary interactions is called YII-sdb hereafter and the one which has switched off binary interactions is called YII-ss. YII-sdb has modeled various binary interactions, including mass transfer, mass accretion, common-envelope evolution, collisions, supernova kicks, tidal evolution, and angular momentum loss, etc.

Here we briefly describe the initial conditions and fiducial models adopted in Yunnan-II model and binary stellar evolution (BSE) algorithm, and more details can be found in their papers. In the building of Yunnan-II model, firstly $2.5\times 10^{7}$ sets of initial conditions for binary evolution are created using Monte Carlo method, including the masses of the primary and the secondary star, orbital separation, eccentricity, and metallicities of stars. The initial mass of the primary is generated by assuming a Chabrier (2003) IMF with cut-offs of $0.1\,{\rm M}_{\odot}$ and $100\,{\rm M}_{\odot}$, and a uniform initial secondary over primary mass ratio distribution is adopted. The initial orbital separation is taken to be constant in logarithm and initial eccentricity is distributed in uniform form. In the fiducial model, $50\%$ field stars are assumed to be in binary systems. Actually, if the orbital separation is wide enough, there would be no interactions between component stars and they are effectively single. The stellar evolutionary tracks of generated binaries can be derived though BSE algorithm. Most of binary interactions are modeled through prescription based approach (like galactic evolution models in semi-analytic models), and the properties of produced objects are comparable with many observations (e.g. cataclysmic variables, Algols, double-degenerates, etc, see Section 4.3 of Hurley et al. 2002 for details). Some model parameters have major influence on the binary evolution, such as the efficiency of common-envelope (CE) ejection $\alpha_{\rm CE}$ and Reimers wind mass-loss efficiency $\eta$. $\alpha_{\rm CE}$ describes the efficiency of the orbital energy dissipated by CE and decides whether the outcome of CE evolution is a close binary or a coalescence. The fiducial $\alpha_{\rm CE}$ in Yunnan-II model is 3.0, increased from 1.0, to compensate the lack of internal energy in Hurley et al. (2002). Mass-loss is another crucial process in stellar evolution, and its rate often determines the final state of a binary system. The empirical formula of Kudritzki & Reimers (1978) is used and fiducial efficiency $\eta$ is set to be $0.3$ according to the observations of Galactic globular clusters (Iben & Renzini 1983). The influence of input model parameters on the integrated SED of stellar population was discussed in detail in the Section 4 of Zhang et al. (2005). They found that the differences between models with different parameters are less than those produced between models with and without binary interactions.

In order to explore the influence of binary evolution on the photometric properties of galaxies, we run our semi-analytic models with and without binary interactions. These two runs share the same semi-analytic models’ parameters but adopting different SPSMs, and the comparison between these two runs will be demonstrated in Section 3.1. Fig. 3 shows the magnitudes as a function of age in YII-sdb and YII-ss for $1\,{\rm M}_{\odot}$ stellar population with solar metallicity, and three photometric bands are shown: GALEX $F_{\rm UV}$ for ultraviolet band, SDSS r for optical band, and 2MASS K for infrared band. In r and K band, when considering binary interactions, the stellar population is fainter for stars older than $0.01\,{\rm Gyr}$. This effect is more obvious in K band than in r band. YII-sdb model is on average $\sim 0.13$ magnitude fainter than YII-ss in K band and $\sim 0.07$ magnitude fainter in r band for stars older than $0.01\,{\rm Gyr}$. While in UV band, the situation is reversed, the stellar population considering binary interactions becomes brighter for stars older than $\sim 1\,{\rm Gyr}$ ($\sim 1.5$ magnitude brighter on average). Because of this difference, binary stars are supposed to be responsible for the UV-excess in elliptical galaxies (Han et al. 2007). In a short summary, when considering binary interactions in a stellar population synthesis model, luminosity is a bit fainter in the optical ($\sim 0.07$ magnitude) and infrared band ($\sim 0.13$ magnitude) but is significantly brighter in the far ultraviolet band ($\sim 1.5$ magnitude).

2.4.2 Dust Extinction

We adopt the dust extinction models of De Lucia & Blaizot (2007), including the extinction by diffuse ISM and birth clouds around young stars. The two reddening terms in Equ.(34) are modeled as the following.

The reddening of a slab geometry ISM is given in Devriendt et al. (1999):

[TABLE]

where $\omega_{\rm\lambda}$ is the albedo of dust (Mathis et al. 1983); $i$ is the inclination angle and can be estimated according to the direction of spin of a gaseous disk assuming the line-of-sight is along the third axis of N-body simulation; $\tau_{\rm\lambda}^{\rm z}$ is the mean face-on optical depth of the average disk at wavelength $\lambda$:

[TABLE]

[TABLE]

The first part of Equ.(36) describes the extinction at the given wavelength relative to the extinction at $5500\AA$ (Mathis et al. 1983); the second part describes the dependence of extinction on gas metallicity ($s=1.35$ for $\lambda<2000\AA$ and $s=1.6$ for $\lambda>2000\AA$, Guiderdoni & Rocca-Volmerange 1987); the last part describes the dependence on redshift (Guo & White 2009). $\langle N_{\rm H}\rangle$ is the mean H column density, and $a=1.68$ in Equ.(37).

The reddening for the birth clouds around young stars is given in Charlot & Fall (2000) and De Lucia & Blaizot (2007):

[TABLE]

[TABLE]

where $\mu$ is a random number drawn from a Gaussian distribution with a center 0.3 and a width 0.2, truncated at 0.1 and 1 (Kong et al. 2004).

2.5 Model Calibration by the Observations

In our whole set of physical models, a few parameters are set free and need to be calibrated according to observational data. Similar to Croton et al. (2006) and Guo et al. (2011), we have 9 free parameters in our models which are summarized in Tab.(LABEL:tab:parameters). We adjust these parameters within the plausible ranges provided by Croton et al. (2006) to fit the stellar mass function at $z=0$, and simultaneously to obtain a reasonable $M_{\rm sph,*}-M_{\rm BH}$ scaling relation and stellar disk sizes. The fiducial values of these parameters are listed in Tab.(LABEL:tab:parameters). By adopting these fiducial values and YII-sdb, we build our fiducial model.

Fig. 4 presents galaxy stellar mass function at $z=0$ in our fiducial model. We can see that the predicted stellar mass function agrees very well with the observational result of Li & White (2009). Only at the most massive end, with $M>10^{11.5}\,{\rm M}_{\odot}$, the number density is slightly higher than observation. To have a better understanding of the inconsistency, we plot the stellar mass of central galaxies against their host halo mass in Fig. 5. Compared to subhalo abundance matching results of Guo et al. (2010), Moster et al. (2013) and Behroozi et al. (2013), our model predicts too massive galaxies in halos $>10^{13.5}\,{\rm M}_{\odot}$. Tuning parameters in the reasonable range can only improve the situation. Note, such an excess at most massive end also exists in other semi-analytic models, such as Croton et al. (2006) and Guo et al. (2011). However more recent semi-analytic models seem overcoming this problem with improved AGN feedback models (Croton et al. 2016 and Henriques et al. 2015). We will explore this in our future work.

Fig. 6 shows the black hole to stellar bulge mass scaling relation in our fiducial model. Our results agree well with the observation yet with slightly smaller scatter than observations. This is acceptable given the fact that the systematic errors of the observation is large.

Fig. 7 is the distribution of the half mass radius of stellar disks for late type galaxies in our fiducial model. Our model prediction is shown as the blue solid curve. While the observational data from Shen et al. (2003) is presented as the black solid and dashed lines. Again, the radius of stellar disks in our model agree well with the observational data in most mass range. There is a slightly disagreement at small mass end, this is acceptable given the fact that the uncertainty in measurement of galaxies size is also large at these mass scales.

3 Results

3.1 The Influence of Binary Evolution

In the following, we examine the influence of adopting Yunnan-II SPSM on the galaxy luminosity function at different bands and galaxy color.

3.1.1 Luminosity Functions

In Fig. 8, we present the luminosity function of galaxies at $z=0$ in the SDSS g, r, i, and z bands. The results with YII-sdb and YII-ss are shown as blue and red lines respectively, while the observed data from SDSS are shown as diamonds with error bars. The ratio of the observed data or the YII-ss over YII-sdb model are shown in bottom panels. We see that the luminosity functions in both models agree well with the observational result over most magnitude scales. At the brightest end, the model prediction is higher than the observational data. Cosmic variance could be partly responsible for that. At the faint end with $M>-18$, our model prediction is lower than the observations because of the resolution is not enough for these faint galaxies in the Milleniuum simulation. Furthermore, in all bands, the difference between two Yunnan-II models are very tiny at faint galaxy end (less than $10\%$), but the difference can be $20\%\sim 50\%$ for galaxies above the character luminosity (the “knee” of luminosity function). Binary interactions have the similar behavior for galaxies’ luminosity and simple stellar population (as shown in Fig. 3). They make galaxies fainter in all optical bands, hence lower the luminosity function. Besides, this effect gets a little bit more obvious at longer wavelength. After considering binary interactions, the average shifts of luminosity functions along x-axis are $0.10$, $0.14$, $0.15$ and $0.16$ magnitude for g, r, i, and z bands respectively.

As shown in Fig. 3, the evolution of a simple stellar population with two Yunnan-II models only differ with each other by $\sim 0.07$ magnitude in optical bands, which is agree with what we see in Fig. 8. Also, it is expected that the difference would be larger in K band and $F_{\rm UV}$ band, especially when the population gets older. Therefore, we checked that in Fig. 9, in which we presented the luminosity function of galaxies from two SPSMs models in GALEX $F_{\rm UV}$ band and 2MASS K band with solid and dashed lines, respectively. The ratio of two models are also presented in the bottom panel. In $F_{\rm UV}$ band, the results of binary evolution model YII-sdb is always higher than YII-ss, this indicates galaxies get brighter if the binary evolution is included. The situation is reversed in $K$ band. Galaxies with binary evolution gets fainter, as what we see in the optical bands, but more obvious. After considering binary interactions, the average shifts of luminosity functions along x-axis are $-0.20$ and $0.18$ magnitude for $F_{\rm UV}$ and K bands respectively.

To further understand the difference on the luminosity function between two Yunnan-II models above, we explore the relation between galaxies’ magnitude and their stellar mass in $F_{\rm UV}$ and K band for two Yunnan-II models in the upper and bottom panel of Fig. 10 respectively. It is notable that, for a given stellar mass galaxy, the $F_{\rm UV}$ luminosity is always brighter when the binary evolution is considered. Specially, for galaxies with $\sim 10^{10.5}\,{\rm M}_{\odot}$ stellar mass at $z=0$, the difference is significant and reaches $\sim 1.5$ magnitude. This overall excess origins from binary interactions, which will bring about an excess in $F_{\rm UV}$ band in old stellar populations. We have check that most of this excess is contributed by early type galaxies, and the $F_{\rm UV}$ luminosity of late type galaxies is less influenced by binary interactions. The results at redshift $z=0,3.06,7.27$ are also shown by different type of lines. It is obvious that the difference is smaller at higher redshift. This is because the difference between two SPSMs gets larger when the populations in galaxies get older. The results is reversed for K band in the bottom panel. The K magnitude in YII-ss is always brighter than that with binary considered for all mass scales and redshifts. The suppression of the formation of red giants by the interactions in binary is supposed to be responsible for that.

In summary, the inclusion of binary evolution will suppress the galaxies’ luminosity in optical and infrared bands (SDSS bands and 2MASS K band) slightly, but increase the luminosity at the ultraviolet band significantly. The displacements of luminosity function along magnitude are $0.10\sim 0.16$, $0.18$ and $-0.20$ in SDSS bands, K band and $F_{\rm UV}$ band respectively. However, given the big scatter in the models (blue shaded region in Fig. 10) and huge uncertainties in the observation, the difference between two SPSMs models are still not distinguishable in the current observational data.

3.1.2 Color Distribution

Furthermore, we checked the color index u - i distribution of our model galaxies in different stellar mass bins in Fig. 11. The results from two SPSMs models, YII-sdb and YII-ss, are presented with blue and red solid curves respectively; and the observational data from SDSS is shown as black dashed curves. Generally, our fiducial model with YII-sdb can recover the bimodality of color distribution in SDSS data in stellar mass range $[10^{9.0},10^{10.5}]\,{\rm M}_{\odot}$. But at the massive end, our model predictions are bluer than observations, and at the faint end, our model galaxies are redder than observed result. As other semi-analytic model pointed out, to get the correct prediction on the color distribution, some physical models need to be updated or even some new physics are needed (Henriques et al. 2015). Moreover, galaxies in YII-sdb is $\sim 0.1$ bluer.

As the influences of binary evolution on galaxies’ flux in $F_{\rm UV}$ and K band is bigger, it is supposed that the color index between these two bands will change significantly when the binary evolution is included. In Fig. 12, the color index $F_{\rm UV}-$K is presented for both Yunnan-II models. The bimodality feature is apparent. With binary evolution, the blue peak is shifted towards blue end by $\sim 0.4$, and the red peak is shifted towards blue end by $\sim 1.6$, which is three times larger than blue peak. The galaxies in YII-sdb is brighter in $F_{\rm UV}$ band because of the existence of hot stars, so they have smaller $F_{\rm UV}-K$ and shifted to the left in every panels comparing to the model prediction of YII-ss. The red peak corresponds to old galaxies which are influenced substantially by binary evolution in $F_{\rm UV}$ band (as shown in Fig. 3), thus the shift of red peak is much larger than the blue peak. If the observational $F_{\rm UV}$ magnitudes could be constrained better in the future, it is very promising that we could tell the difference between these two model predictions with the observations.

3.2 Remnants of Binary Evolution

The remnants of the evolution of binary stars (e.g. NS-NS, NS-BH, BH-BH) are thought to be associated with many high energy events, like type Ia supernovae, short gamma-ray bursts, kilnovae and gravitational wave events. In our model, using Yunnan-II binary SPSM, we are able to study these remnants in a cosmological context. As an example, in Fig. 13 we show the cosmic number density of the double compact objects (hereafter DCOs), i.e. NS-NS, NS-BH and BH-BH as a function of cosmic time. With the number density and merger rate of these DCOs, we could study some high energy events associated with mergers of them, e.g. gravitational wave signals, kilonovae, short gamma-ray bursts etc. We will present these results in our forthcoming papers.

4 Conclusions

In this work, we develop a new semi-analytic galaxy formation model GABE. The model is based on merger trees provided by cosmological N-body simulations. It succeeds many physical recipes from existing successful models and also includes some new ingredients. For example, it adopts Yunnan-II SPSM which carefully considers binary evolution, and a cooling table which takes into account UV heating properly.

When calibrating our model to the observed stellar mass function, our model reproduces a large body of observations. We explore the influence of Yunnan SPSM model with binary evolution on the galaxy luminosity function and color. Our finding is that, with binary SPSM, the influence on luminosity is quite small over different SDSS bands, and is more substantial in infrared band and ultraviolet band. The shifts of luminosity function along magnitude after considering binary evolution are $0.10\sim 0.16$, $0.18$ and $-0.20$ magnitude in SDSS bands, K band and $F_{\rm UV}$ band respectively. Galaxies appear to be bluer with binary evolution, especially when $F_{\rm UV}$ band is under consideration. The shift of u - i color index is about $-0.1$. The shift of the blue peak of $F_{\rm UV}-K$ color index is about $-0.4$, and the one of the red peak is about $-1.6$, three times larger than the blue peak.

This is the first paper of our series works. We have used GABE to simulate the gravitational wave signals of supermassive black holes in Wang et al. (2019). Besides, the new feature of our model allows us to predict the population of double compact objects (e.g. NS-NS, NS-BH, BH-BH) which are expected to be associated with some high energy events in a cosmological context. For example, the short gamma-ray bursts which are expected to be products of NS-NS or NS-BH mergers; gravitational wave events which are supposed to be generated during mergers of two compact objects in binary systems. In addition, we will be able to explore the relation between these high energy events and their host galaxies. We will present these studies in our future works.

Acknowledgements.

We acknowledge Zhanwen Han for reading our draft and providing some very useful comments and suggestions. We also would like to thank the anonymous referee for an insightful report which led to improvements in the presentation and content of our paper. We acknowledge support from the National Key Program for Science and Technology Research and Development (2015CB857005, 2017YFB0203300) and the Chinese Natural Science Foundation (NSFC) grants (11390372, 11425312, 11503032, 11573031, 11851301, 11873051, 11573062, 11521303, 11390734, 11573033, 11622325 and 11573030). Besides, FZ acknowledges support from the YIPACAS Foundation (grant 2012048) and the Yunnan Foundation (grant 2011CI053). QG is supported by the Newton Advanced Fellowship. JP acknowledges support from the National Basic Research Program of China (program 973 under grant no. 2015CB857001).