text:
# BINARY AGGREGATIONS IN HIERARCHICAL GALAXY FORMATION: THE EVOLUTION OF THE GALAXY LUMINOSITY FUNCTION

## 1. INTRODUCTION

Understanding of galaxy formation has undergone impressive developments in recent years. On the observational side, a key step in the study of the statistical properties of high-redshift galaxies has been taken with the discovery of the ”Lyman-break” galaxies at $`z34`$ (Steidel et al. 1996; Adelberger et al. 1998); other important steps have been taken with the estimate of the cosmic star formation rate (see Madau, Pozzetti & Dickinson 1998; Fontana et al. 1999), and the measurement of the galaxy luminosity functions (LFs) at redshifts up to $`z4`$ for UV luminosities extending over the range $`24M18`$ (Steidel et al. 1999 from spectroscopic data, Pozzetti et al. 1998 and Poli et al. 2001 from photometric redshifts).

These results, with their unprecedented span in both cosmic time and magnitude, boosted the developement of an “ab initio” theory of galaxy formation. The theory has been firmly rooted into the cosmological framework with the development of semi-analytic models (SAMs, Kauffmann et al. 1993; Cole et al. 1994; Somerville & Primack 1999; Poli et al. 1999; Wu, Fabian & Nulsen 2000; Cole et al. 2000) which link in a single computational structure two main blocks: the dynamical history of galaxies as they emerge from primordial dark matter (DM) density perturbations and grow through hierarchical merging events; the baryonic processes in the galactic structures, namely, gas cooling, star formation, rise and fall of stellar populations, energy feedback from Supernovae.

Such advances in both theory and observations are pushing the comparison between models and data to higher and higher degrees of accuracy. A key role is played by the $`z`$-resolved LFs. In fact, the earliest SAMs were calibrated to fit the locally observed LF (still with considerable uncertainties in the normalization and in the slope at the faint end), and were tested at higher redshift through integrated counts and $`z`$-distributions. The recent observations are producing a progressive convergence of the local LFs found by different groups (see Zucca et al. 1997, Cross et al. 2001) down to $`M_B16`$, with the additional indication that at the faint end the shape is more complex than represented by a Schechter form (see Marzke, Huchra & Geller 1994; Loveday 1997). In addition, the resolved LFs now begin to describe the evolution of the galaxy population out to $`z4`$, thus providing a differential test for the model predictions.

The comparison of the models with such data supports the grand design of hierarchical galaxy formation, but indicates that some of the processes included in the SAMs require an improved treatement. For example, the LFs predicted at $`z3`$ overestimate the number of faint ($`M_{1700}17`$) galaxies by a factor $`58`$ when compared with observations based on photometric reshifts (see Somerville, Primack & Faber 2001; Poli et al. 2001), and the excess goes beyond the estimated incompleteness of the data faintwards of $`m_B26`$. A similar trend, though less evident, is found in the redshift range around $`z1`$ (Poli et al. 2001), while at $`z0`$ the observed upturn in the LF at the faint is not accounted for in the current SAM predictions. Addressing the above critical points is crucial to fully understand the physical processes driving the galaxy evolution.

The shape of the predicted LF at faint and intermediate luminosities is affected by two main processes. The first concerns the effect of Supernovae (SN hereafter) which heat and partially expell the galactic gas. A stronger SN feedback would suppress star formation in smaller haloes thus decreasing their B and UV luminosity and flattening the LF at the faint end. However, in the framework of the simple parametrizations currently adopted in SAMs it is extremely difficult to increase such feedback without destroying the agreement of the model with other observables; in particular, the predicted luminosities for small spiral galaxies would be too faint when compared with existing data concerning the Tully-Fisher relation (see Cole et al. 2000). Attempts at non- parametric treatements of the feedback have been made by Mac Low & Ferrara (1999); Goodwin, Pearce & Thomas (2000); Monaco (2001). Different sources of feedback, like that arising from the photoionization of the inter-galactic medium by stars and quasars, have been recently included in the SAM framework (Benson et al. 2001).

The other component affecting the shape of the LF is the mass distribution of the galaxies, which is determined by the detailed dynamical processes taking place inside the host DM haloes. Among these, tidal stripping and binary aggregations of satellite galaxies play a relevant role. However, the former affects mainly the very low-mass end of the mass distribution at circular velocities $`v2050`$ km/s, and has a minor impact on the faint end of the LF (see Benson et al. 2001); so we shall focus here on binary aggregations. Such a process is treated in current SAMs only under the approximation of orbital decay toward a central dominant galaxy. Here we extend the treatement to include aggregations between all galaxies in common DM haloes. On the other hand, we shall adopt the standard SAM prescriptions concerning star formation to derive the galaxy LF from the mass distribution; no additional starburst recipes are associated with the aggregation events between satellite galaxies. This allows us to single out the dynamical effects of our description without introducing new (and uncertain) free parameters in the model, as we discuss in §6.

The paper is organized as follows: the derivation of the mass function in the framework of SAMs, and the motivation for our improved treatment are discussed in §2. A technical description of our approach is presented in §3, where we show how we fit our treatment of binary aggregations into the canonical SAM framework. In §4 we briefly recall the basic prescriptions that we share with other SAMs to correlate the gas cooling, the star formation and evolution, and the SN feedback with the dynamical history of the galaxies. The LFs we derive are compared with the data in §5. In §6 we discuss our results and present our conclusions.

## 2. THE DYNAMICAL EVOLUTION OF GALAXIES: OVERVIEW AND SPECIFIC MOTIVATIONS

In SAMs the galaxy mass distribution is derived from the merging histories of the host DM haloes, under the assumption that the galaxies contained in each halo coalesce into a central dominant galaxy if their dynamical friction timescale is shorter than the halo survival time; the surviving galaxies (commonly referred to as satellite galaxies) retain their identity and continue to orbit within the halo. While the histories of the DM condensations rely on a well established framework (the extended Press & Schechter theory, EPST, see Bower 1991; Bond et al. 1991; Lacey & Cole 1993), the recipe concerning the galaxy fate inside the DM haloes is guided by a posteriori consistency with the outputs of high-resolution N-body simulations.

In reality, additional dynamical processes complement the dynamical friction in driving the evolution of the mass distribution. Among these, binary aggregations between satellite galaxies in common haloes have previously been considered by Cavaliere, Colafrancesco & Menci (1991, 1992); Cavaliere, & Menci (1993), as a process that would flatten the shape of the mass function (and hence of the LF) at small/intermediate masses. This is because such masses aggregate into larger units, while the large masses are so few that their binary encounters are unlikely. In the above papers, the aggregation-driven evolution of the mass distribution was computed in terms of the non-linear Smoluchowski kinetic equation with a mass-dependent aggregation rate. The results showed that the effectiveness of the process depends critically on the environment, which in those models was a given input. To describe at the same time the galaxy dynamics and the evolution of the host haloes by hierarchical clustering, high-resolution N-body simulations or SAMs are required.

On the N-body side, recent works (see Klypin et al. 1999 for pure DM, and Murali et al. 2001 for hydrodynamical simulations) indicate a complex galaxy growth. While at very low circular velocities ($`v2050`$ km/s) tidal stripping may affect the mass distribution (Gnedin & Ostriker 1997; see also Taylor & Babul 2000; Taffoni et al. 2001), at larger masses ($`v100`$ km/s) binary aggregations of satellite galaxies do play a relevant role which complements the coalescence into a central galaxy through dynamical friction. Indeed, Murali et al. (2001) analyze the competing role of the two growth modes of the simulated galaxies in terms of an evolutionary equation which includes also a binary aggregation term of the Smoluchowski type.

On the SAM side, efforts to insert the aggregations between satellite galaxies have recently been started by Somerville & Primack (1999). However, they adopt a cross section (derived from the N-body simulations of galaxy encounters by Makino & Hut 1997) valid for equal galaxies in clusters with velocity dispersion much higher than the internal galaxy dispersion, a condition that allowed them to adopt a one- body treatement for the aggregations.

Here we take a step forward, and consider in closer detail the actual two-body dynamics of aggregations; we shall also adopt a cross section valid down to small groups with velocity dispersions close to those internal to galaxies, as is the case at early times in the hierarchical clustering picture.

In fact, we develop a SAM including both dynamical friction and binary aggregations. Instead of employing a Monte Carlo simulation, as usual for SAMs, we follow the evolution of the galaxy mass distribution by solving numerically a set of evolutionary equations, as in Poli et al. (1999). The subset describing the two-body dynamics is constituted by the (non-linear) kinetic Smoluchowski equation. This modifies the mass function resulting from dynamical friction by an amount which depends on the properties of the host DM haloes, which in turn evolve according to the EPST theory. The basic description of gas cooling, star formation and evolution, and SN feedback is kept unchanged with respect to the standard SAM prescriptions given, e.g., by Poli et al. (1999) and Cole et al. (2000).

We first compare our results with those from existing N-body simulations, and then compare our LFs with those obtained from spectroscopic and deep imaging surveys down to faint magnitutes ($`I_{AB}27.2`$), at redshifts up to $`z4`$. This allows us to test the effects of binary aggregations on the SAM predictions over a wide range of cosmic times and masses.

## 3. THE GROWTH OF THE GALAXY MASSES

We consider the number $`N(m,M,t)dmdM`$ per Mpc3 of galaxies with mass in the range $`dm`$ about $`m`$, residing in haloes with masses in the range $`M`$ to $`M+dM`$ at cosmic time $`t`$. At the initial time $`t_0`$ we assign one galaxy to each halo, a condition which formally translates into $`N(m,M,t_0)=N_H(M,t_0)\delta (mM)`$ where $`\delta `$ is the Dirac delta function. Our default choice for the halo mass distribution $`N_H(M,t)`$ is the standard Press & Schechter (1974) expression, whose dependence on the cosmological parameters and on the spectrum of primordial density perturbation is recalled in Appendix A; but we shall also explore the effects of adopting the Sheth & Tormen (1999) mass distribution, also recalled in Appendix A.

Given the merging history of DM haloes described by EPST, two main processes affect the evolution of $`N(m,M,t)`$: the orbital decay of satellite galaxies onto a central dominant galaxy due to dynamical friction, and the binary aggregations between satellite galaxies. In the following two subsections we investigate the effects of the above processes on the evolution of $`N(m,M,t)`$. A schematic view of the two dynamical processes driving the growth of galaxies in our model is given in fig. 1. The corresponding evolution of the mass distribution of galaxies is presented in the last subsection 3.4.

Fig. 1. - A schematic representation of the various terms contributing to the evolution of the galaxy velocity function $`N(v,V)`$ inside DM haloes of circular velocity $`V`$. The first three processes marked with DF correspond to the construction and destruction terms (eq. 1) from dynamical friction following the merging of the host haloes (see text). The last two processes, marked with BA, represent binary aggregations of satellite galaxies (see eq. 3) inside the DM host halo. Their combined effects drive the evolution of $`N(v,V,t)`$.

### 3.1. Coalescence driven by dynamical friction

In the following we shall frequently adopt the circular velocities as the proper quantity to mark the depth of the potential wells of galaxies and of their host haloes. As for the latter, this is related to the mass through the relation $`V=\sqrt{GM/R}`$ where the limiting radius $`R`$ is the radius within which the mean mass density is $`200\rho _c`$, where $`\rho _c`$ is the critical density at the redhsift $`z`$ where the halo is identified. The relation between $`R`$ and $`M`$ is given, e.g., in Navarro, Frenk & White (1997) and Mo, Mao & White (1998) and takes the form $`R=\mathrm{1.63\hspace{0.17em}10}^2(M/h^1M_{})^{1/3}[\mathrm{\Omega }_0/\mathrm{\Omega }(z)]^{1/3}(1+z)^1h^1kpc`$ for Einstein-de-Sitter ($`\mathrm{\Omega }=1`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0`$), open ($`\mathrm{\Omega }_0<1`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0`$), and flat ($`\mathrm{\Omega }_0+\mathrm{\Omega }_\mathrm{\Lambda }=1`$) universes.

As to a galaxy inside a host DM halo, the circular velocity $`v`$ is related to $`m`$ through $`m=v^2r_{tid}/G`$; here $`r_{tid}`$ is the tidal radius, within which the mean density of the galaxy exceeds the average density of the host halo interior the pericentre of the galactic orbit. This ensures that the galactic subhalo survives the tidal stripping due to the host halo potential wells and retains its own identity. Our computation of $`r_{tid}`$ is given in Appendix B.

The above relations allow us to relate the mass and the circular velocity distributions by applying the proper Jacobian.

The evolution of $`N(v,V,t)`$ driven by dynamical friction is computed as follows. At each time step, we first compute the conditional probability $`d^2P_H(V^{}V,t)/dV^{}dt`$ that a given halo with circular velocity $`V`$ at time $`t`$ has a progenitor with circular velocity $`V^{}`$ at $`t^{}=t\mathrm{\Delta }t`$; this is calculated in the framework of the EPST, and its expression (depending on cosmology and on the perturbation spectrum as given in Bower 1991; Bond et al. 1991; White & Frenk 1991, Lacey & Cole 1993), is recalled in Appendix A for a generic $`t^{}`$ (see eq. A-3). Similarly, from the EPST we compute the inverse conditional probability $`d^2P_H(VV^{},t)/dV^{}dt`$ that a given halo with circular velocity $`V`$ at time $`t`$ ends up in a halo with circular velocity $`V^{}`$ at $`t+\mathrm{\Delta }t`$ (see eq. A-2).

According to the canonical prescriptions usually adopted by SAMs, during halo merging a galaxy contained in one of the parent halos contributes to enrich the dominant galaxy of the common halo of circular velocity $`V`$ if its coalescence time $`\tau _{df}(v)`$ is shorter than the halo survival time $`\tau _l(V)`$ predicted in the EPST. The circular velocity of the merger is then equated to that of the host halo. Or else, for $`\tau _{df}>\tau _l`$, the galaxy retains its identity. Thus the increment in the number $`N(v,V,t)`$ of galaxies of given $`v`$ in halos of given $`V`$ is linked to the increment of DM halos with $`V`$, whose contribution from progenitors $`V^{}`$ is $`N_H(V^{},t)d^2P_H(V^{}V,t)/dV^{}dt`$; the link occurs through the probabilities: i) of forming one dominant galaxy of circular velocity $`v=V`$ by coalescing galaxies contained in the parent halos with lower velocities $`v^{}<v`$; ii) of finding galaxies with velocities $`v`$ which have not coalesced despite the merging of their host haloes. The corresponding decrement is due to the inclusion of galaxies with current velocity $`v`$ into a larger halo $`V^{}>V`$. The construction and destruction terms described above are schematically represented in the upper part of fig. 1. The corresponding evolution of $`N(v,V,t)`$ in a single timestep $`\mathrm{\Delta }t`$ is expressed by

$`N(v,V,t+\mathrm{\Delta }t)N(v,V,t)=`$ (1)
$`=`$ $`\mathrm{\Delta }t\delta (vV){\displaystyle _0^v}𝑑v^{}{\displaystyle _v^{}^V}𝑑V^{}N_H(V^{},t){\displaystyle \frac{dP_H(V^{}V,t)}{dV^{}dt}}{\displaystyle \frac{N(v^{},V^{})}{N_T(V^{})}}prob\left[\tau _{df}(v^{})<\tau _l(V)\right]`$
$`+`$ $`\mathrm{\Delta }t{\displaystyle _v^V}𝑑V^{}N_H(V^{},t){\displaystyle \frac{d^2P_H(V^{}V,t)}{dV^{}dt}}{\displaystyle \frac{N(v,V^{})}{N_T(V^{})}}\left\{1prob\left[\tau _{df}(v)<\tau _l(V)\right]\right\}`$
$``$ $`\mathrm{\Delta }t{\displaystyle _V^{\mathrm{}}}𝑑V^{}{\displaystyle \frac{d^2P_H(VV^{},t)}{dV^{}dt}}N(v,V),`$
where $`N_T(V)=𝑑v^{}N(v^{},V,t)`$ is the total number of galaxies per $`\mathrm{Mpc}^3`$, per unit halo rotational velocity. The first term on the right-hand side implies that if any galaxy with circular velocity $`v^{}<V`$ is included in a halo with circular velocity $`V`$ and its orbit decays (due to dynamical friction) to the center of such a halo within the halo survival time, then a central galaxy with velocity $`v=V`$ is formed inside the halo.

Coalescence is caused by loss of galaxy energy and orbital angular momentum due to dynamical friction to the halo material. The timescale of such process is usually determined through orbit averaging of the Chandrasekhar formula (see Lacey & Cole 1993); the result is given in Cole et al. (2000) as follows:

$$\tau _{df}=\mathrm{\Theta }\tau _{dyn}\frac{0.3722}{ln(M/m)}\frac{M}{m}.$$
(2)
Here $`\tau _{dyn}\pi R/V`$ is the dynamical time of the halo, and $`\mathrm{\Theta }=[J/J_c(E)]^{0.78}[r_c(R)/R]^2`$ contains the dependence on the initial energy $`E`$ and angular momentum $`J`$ of the galaxies, in terms of the angular momentum $`J_c`$ and the radius $`r_c`$ of the circular orbits corresponding to $`E`$. Although the values of $`\mathrm{\Theta }`$ are statistically distributed (approximately following a lognormal function with $`log_{10}\mathrm{\Theta }=0.14`$ and $`(log_{10}\mathrm{\Theta }log_{10}\mathrm{\Theta })^2^{1/2}0.26`$, see Tormen 1997), in the following we shall use for $`\mathrm{\Theta }`$ its average value. Note that the orbit averaging procedure (adopted to derive eq. 2) is not appropriate for very elongated orbits for which the local density varies on a much shorter timescale than that given by eq. (2); however, such orbits are expected to be a minority ($`15\%`$, see Ghigna et al. 1998).

The probability $`prob\left[\tau <\tau _l(V)\right]`$ for a DM halo of velocity $`V`$ to have a survival time $`\tau _l`$ larger than a given value $`\tau `$ has been computed by Lacey & Cole (1993, see their eq. 2.21) in the framework of the EPST, and is recalled in Appendix A.

### 3.2. Binary aggregations of satellite galaxies

In addition to the above coalescence process, we include binary aggregations. So for each DM halo with circular velocity $`V`$, we compute the further evolution of the galaxy mass distribution in a timestep due to binary aggregations. This is described by the Smoluchowski equation, which we write in terms of masses for the sake of simplicity:

$`N(m,M,t+\mathrm{\Delta }t)`$ $``$ $`N(m,M,t)=`$ (3)
$`{\displaystyle \frac{1}{2}}\mathrm{\Delta }t{\displaystyle _0^m}𝑑m^{}N(m^{},M,t)N(mm^{},M,t)\tau _{agg}^1(m^{},mm^{},V)`$
$``$ $`\mathrm{\Delta }tN(m,M,t){\displaystyle _0^{\mathrm{}}}𝑑m^{}N(m^{},M,t)\tau _{agg}^1(m,m^{},V),`$
where $`\tau _{agg}^1(m,m^{},V)`$ is the aggregation rate depending on the DM halo where the galaxies $`m`$ and $`m^{}`$ reside. The first term describes the construction of galaxies with mass $`m`$ from smaller ones with mass $`m^{}`$ and $`mm^{}`$, while the second represents the destruction of galaxies $`m`$ due to their aggregation with others. A schematic representation of the two terms governing the binary aggregations is given in the lower part of fig. 1.

The aggregation rate is governed by $`\tau _{agg}^1=\mathrm{\Sigma }V_{rel}/(4\pi R^3/3)`$ (see Cavaliere, Colafrancesco & Menci 1992). Here $`V_{rel}`$ is the average relative velocity of galaxies in the DM halo whose rms value is equal to twice the halo 1-D velocity dispersion $`\sigma _VV/\sqrt{2}`$, and $`\mathrm{\Sigma }`$ is the cross section. The latter is given for nearly grazing, weakly hyperbolic encounters by Saslaw (1985), and by Cavaliere, Colafrancesco, & Menci (1992). It includes a geometrical term (proportional to the area of the galaxies), and a focussing factor $`1+(v/V_{rel})^2`$ that accounts for the enhancement of $`\mathrm{\Sigma }`$ in slow encounters with resonance between the internal and the orbital degrees of freedom (see also Binney & Tremaine 1986). Thus, the average rate for binary aggregations is

$$\tau _{agg}^1=\pi (r^2+r^2)\left(1+G\frac{m+m^{}}{r+r^{}}\frac{1}{V_{rel}^2}\right)\frac{V_{rel}}{4\pi R^3/3},$$
(4)
where the average is over the relative velocities $`V_{rel}`$. The distribution of the encounter velocities is assumed to be Maxwellian, namely

$$g(V_{rel})=\sqrt{\frac{2}{\pi }}\frac{V_{rel}^2}{(\sqrt{2}\sigma _V)^3}e^{V_{rel}^2/4\sigma _V^2}.$$
(5)
Note that for encounters between equal galaxies with $`r^{}=r`$ and $`N(M^{})=N(M)=N`$ the scaling of the aggregation rate (4) reduces to $`\tau _{agg}^1r^2(1+v^2/V_{rel}^2)V_{rel}`$. Performing the average over the distribution $`g(V_{rel})`$ in terms of the rescaled variable $`yV_{rel}/v`$ yields $`\tau _{agg}^1r^2v^4\sigma _V^3R^3I(x)`$, where the function $`I(x)𝑑yy^3exp(y^2/x^2)+𝑑yyexp(y^2/x^2)`$ tends to a constant when the ratio $`x\sigma _V/v\mathrm{}`$. Thus the aggregation rate (4) reduces to the expression given by Makino & Hut (1997) – originally derived analytically by Mamon (1992) – and the r.h.s. of eq. (3) becomes proportional to $`N^2r^2v^4\sigma _V^3R^3`$ (the effective rate adopted by Somerville & Primack 1999) in the proper limits of large encounter velocities relative to the internal galaxy velocity dispersion, and of encounters between equal galaxies.

Finally, note that after performing the average of eq. (4) over the distribution $`g(V_{rel})`$ only $`\sigma _V`$ enters the computation; in other words, we use an average description for encounter velocities typical of the environment considered. The geometrical cross section contained in eq. (4) does not properly describe single events with $`V_{rel}v`$, but these are rare for our typical values of $`\sigma _V/v`$ in the range from $`1`$ to about $`4`$. In cases with larger ratios $`\sigma _V/v`$ the effect of the aggregations is suppressed since $`\tau _{agg}\sigma _V^3`$ as shown above, but even inside rich clusters the averaged cross section that we adopt is consistent with the simulation results, see Makino & Hut (1997); for numerical computations of the role of interactions inside clusters see also Lanzoni (2000).

In sum, our cross section provides an accurate average description of aggregations in systems where the velocity dispersion is close to the galaxy circular velocity and the focussing term is relevant; on the other hand, it constitutes an effective average approximation for encounters in environments with large velocity dispersion up to the scale of rich clusters, where the aggregations are disfavoured anyway. As a final global check, we have verified that the insertion of an artificial cutoff at $`V_{rel}=4v`$ in $`\mathrm{\Sigma }`$ does not change our results.

### 3.3. Numerical solutions of the equations: test cases

The equations (1) and (3) describing the evolution of $`N(m,M,t)`$ are integrated numerically on a grid of circular velocities and cosmic times with step $`\mathrm{\Delta }t=10^2H_o^1`$.

In order to test our numerical code, we first run the computation for two relevant simple cases where analytic solutions are available. In the limit $`\tau _{df}0`$ (i.e., when merging of the host haloes is promptly followed by coalescence of the galaxies within them by dynamical friction) the solution $`N(v,V,t)`$ of eq. (1) when integrated over the circular velocity $`V`$ must yield the Press & Schechter mass distribution. The comparison between the numerical and the analytic solution in this case is performed at three different times in the top panel of fig. 2, which shows that the numerical solutions remain close to the analytic Press & Schechter form over the whole range of $`v`$; in fact, the relative deviation is always smaller than 5%.

To test the code section concerning the Smoluchowski equation, we numerically solve eq. (3) in the case of constant aggregation rate $`\tau _{agg}^1=\mathrm{const}`$ for galaxies within a host halo of given mass $`M`$. In this case the exact solution (Smoluchowski 1916; Trubnikov 1971) has the form $`N(m,t)=[A_o/m_{}^2(t)]exp[m/m_{}(t)]`$, where the constant $`A_o`$ is related to the total mass $``$ contained in the system, and $`m_{}(t)=m_0(t/t_0)`$ is a characteristic mass linearly growing with time. Trubnikov (1971) has shown that such a solution holds for very general initial conditions after a transient time. The comparison with the numerical solution is performed in the bottom panel of fig. 2, where the initial condition (dotted line) has been chosen to be $`N(m,t_0)(m/m_0)^1exp[m/m_0]`$ normalized as to yield a total mass $`=\mathrm{1.54\hspace{0.17em}10}^2m_0`$. Note how the numerical solution progressively flattens at small masses during the transient, to approach the exact solution (computed at $`t=3t_0`$) with a relative error smaller than $`3\%`$ over the whole mass range. Note also that the numerical solution conserves the total galaxy mass with high accuracy.

Fig. 2 - Comparison between the numerical and the analytic solutions of eqs. (1) and (3) for two test cases. Top panel: the numerical solutions of eq. (1) in the limit of vanishing $`\tau _{df}`$ (dashed curves) are compared with the Press & Schechter form (solid curves). The curves refer to redshifts $`z=8`$, 2.6 and 0 (from left to right at the high-mass end). Bottom panel: the numerical solution of the Smoluchowski eq. (3), computed for a constant aggregation time (dot- dashed curves) at 6 equally spaced time intervals from $`t_0`$ to $`3t_0`$, are compared with the corresponding analytic solution (given in the text) computed at $`t=3t_0`$ (heavy solid line). The dotted line represents the initial condition. The inset shows the variation with time of the total mass $``$ corresponding to the numerical solution.

### 3.4. The evolution of the galaxy circular velocity distribution

Having tested our code, we proceed to compute the complete evolution of $`N(v,V,t)`$. At each time step, we first compute the change of $`N(v,V,t)`$ due to dynamical friction (the right hand side of eq. 1), and then the further change due to binary aggregations of satellites (the right hand side of eq. 3) using the physical values of $`\tau _{df}`$ and $`\tau _{agg}`$ given in eqs. (2) and (4). A discussion of the role of the two terms is given in the final §6.

Once the complete evolution of $`N(v,V,t)`$ is found for all times, we compute the probability of finding a galaxy with circular velocity $`v`$ in a halo of circular velocity $`V`$ as $`f(v,V,t)dv=N(v,V,t)dv/N_H(V,t)`$. Then the total density distribution of galaxies with circular velocity $`v`$ is computed from $`N(v,t)=_v^{\mathrm{}}𝑑VN_H(V,t)f(v,V,t)`$, where the halo distribution $`N_H(V,t)`$ takes the canonical Press & Schechter form. This is the distribution of circular velocities irrespective of the halo to which the galaxies belong (the global velocity distribution).

The “transition probability” for galaxies is given by $`p(v^{},t^{},v,t)=_v^{}^{\mathrm{}}_v^{\mathrm{}}𝑑V𝑑V^{}\frac{dP_H(V^{},t^{},V,t)}{dV^{}}f(v^{},V^{},t^{})f(v,V,t)`$, where $`\frac{dP_H(V^{},t^{},V,t)}{dV^{}}`$ is the fraction of mass in haloes of circular velocity $`V^{}`$ at time $`t^{}`$, and later at $`t>t^{}`$ in haloes with $`V`$; this is given by EPST (see Appendix A).

The evolution of the galaxy circular velocity distribution $`N(v,t)`$ resulting from the full dynamics is shown in fig. 3 for a CDM power spectrum of primordial density perturbations and for our reference set of cosmological/cosmogonical parameters: $`\mathrm{\Omega }_0=0.3`$, $`\mathrm{\Omega }_\lambda =0.7`$, and Hubble constant $`h=0.7`$ in units of 100 km/s/Mpc. For comparison, we also show the evolution corresponding to dynamical friction alone (as usually considered in SAMs), and the evolution of the halo velocity distribution as given by the Press & Schechter formula.

Fig. 3 - The galaxy velocity function at four different redshifts. The velocity is measured in units of 100 km/s. The solid line is the velocity function resulting from the full model including dynamical friction and satellite aggregation (eqs. 1 and 3). The function $`N(v)`$ resulting from dynamical friction alone is shown as a dashed line, while the velocity function of the host haloes (computed from the Press & Schechter formula) is also shown for comparison as a dotted line. In the bottom-right panel we also show the velocity function resulting from the N-body simulations by Klypin et al. (1999); the solid squares correspond to a simulation with a box size of $`60h^1`$ Mpc, and the open circles to a box size of $`7.5h^1`$ Mpc which allows for a finer mass resolution. Both cases are computed for $`\mathrm{\Lambda }`$-CDM cosmology with parameters as in the text.

First consider the effect of dynamical friction. Compared to the halo velocity distributions, the galaxy distribution shows a delayed evolution for $`z<3`$, particulary evident at low redshift. This is due to the longer timescale on which dynamical friction occurs (eq. 2) as compared to the survival time of the halo (which is of order $`\tau _{dyn}`$). As a consequence, while at high redshift the merging of haloes is promptly followed by coalescence of the galaxies within them, at later times the number of galaxies accumulating in the haloes increases due to the longer time taken by dynamical friction to occur. At small masses this implies more galaxies than haloes, while at large $`v`$ the increase in the number of haloes is not followed by a corresponding increase of massive galaxies since these are formed on a longer timescale.

The effect of binary aggregations is to flatten the slope of the galaxy velocity distribution at the low-mass end ($`v50150`$ km/s). Such a process is already effective at $`z2.5`$ and continues down to $`z0`$, where our $`N(v)`$ is fully consistent with available results from N-body simulations as shown by fig. 3. However, in the range from $`z1`$ to $`z0`$, the aggregations deplete more efficiently objects with intermediate mass ($`v100150`$ km/s), thus producing an upturn of the local velocity distributions at low $`v`$. This is because at low $`z`$ the galaxies are hosted in haloes that are typically more massive than at higher $`z`$, according to the hierarchical clustering. The large relative velocities $`V_{rel}`$ that galaxies have in such deeper potential wells suppress the aggregation of galaxies with low $`v`$ due to the $`v/V_{rel}`$ term in the cross sections (4), as stressed by Cavaliere & Menci (1997).

To further assess the consistency of our results with N-body simulations and to guide our interpretation of the aggregation role we plot in fig. 4 the overall merging rates predicted by our model and compare them with those obtained in recent hydrodynamical simulations in a cosmological framework (Murali et al. 2001). In the top panel we show the net number increment as a function of the cosmic time for all galaxies with baryonic mass larger than $`\mathrm{5.6\hspace{0.17em}10}^{10}M_{}`$ (the baryonic mass associated with the DM mass in our model will be derived in §4). We also show as a dashed line the construction rate (i.e., the sum of the positive terms in eq.  and , also illustrated in fig. 1), to show that it dominates over the destruction terms for such galactic masses.

Fig. 4 - Contributions to the number density of galaxies with masses above $`\mathrm{5.6\hspace{0.17em}10}^{10}M_{}`$. Top panel: the net rate (creation minus destruction) predicted by the full model including dynamical friction and aggregations of satellite galaxies is shown as a solid line. This is compared to that obtained from the simulations by Murali et al. (2001, represented by the dots) with the same cosmological parameters adopted here. The dashed line shows the construction rate alone, as predicted by the model. Bottom panel: the destruction rate (the difference between the two curves in the top panel) resulting from the model (solid heavy line) is compared with the simulations (dots). The contributions to such destruction rate from dynamical friction (dotted line) and from aggregations of satellite (dashed line) are also shown.

Notice how simulations and our model agree in yielding a net rate which declines with time after a peak at $`z2.53`$. The decline of such a rate for $`z1`$ is consistent with the observed one (see LeFevre et al. 2000; Carlberg et al. 2000) within the uncertainties entering the comparison.

To assess the relative role of dynamical friction and satellite aggregations we show in the bottom panel the destruction rate for the same mass threshold, compared with the N-body results. Note the smooth decline at low redshift in both the simulations and the model (the heavy solid line) after the peak at $`z2.5`$. Fig. 4 shows in detail that the position of the peak is mainly determined by dynamical friction; satellite aggregations begin to contribute significantly to the destruction rate at $`z2.53`$. Then they sustain the destruction rate during its slow decline to lower $`z`$. Thus, first the hierarchical clustering drives the construction (top panel) and the destruction (bottom panel) of the basic galactic blocks within a short stretch of cosmic time. Then binary aggregations further deplete the number of galaxies. The evolution of the destruction rate corresponding to binary aggregations is due to the shifting balance of different terms: first, the aggregation efficiency increases due to the growth of the galaxy sizes and of the number of galaxies contained in the host halos (see the aggregation rate in eq. 4); at later times, the above terms are balanced and/or overwhelmed by the increase of the galaxy relative velocities, which grow in time following the increasing mass of the host haloes and tend to suppress the binary aggregation rate eq. (4).

Such a picture is confirmed by the fact that when a lower mass threshold is considered, the high-$`z`$ shape of the destruction rate remains similar to that in fig. 4, while the low-$`z`$ tail is appreciably suppressed. This is because for small galaxies the effect of the increase of the relative velocities with time is stronger, and prevents binary aggregations from sustaining the destruction rate at small $`z`$. So at $`z0`$ the velocity functions in fig. 3 has an upturn at small masses.

Having discussed the full dynamics, we now turn to recall the basic points of the stellar section of our SAM.

## 4. STAR FORMATION AND EVOLUTION

To link the stellar section of the model with the dynamics we adopt the standard procedure commonly used in SAMs (Somerville & Primack 1999; Cole et al. 2000); we shall give here a brief presentation of how such a procedure fits into our statistical treatement of the evolution of the mass distribution of galaxies inside host haloes.

The baryonic content $`(\mathrm{\Omega }_b/\mathrm{\Omega }_m)m`$ of the galaxy is divided into a hot phase with mass $`m_h`$ at the virial temperature $`T=(1/2)\mu m_Hv^2/k`$ ($`m_H`$ is the proton mass and $`\mu `$ is the mean molecular weight), a cold phase with mass $`m_c`$ able to radiatively cool within the galaxy survival time, and the stars (with total mass $`m_{}`$) forming from the cold phase on a time scale $`\tau _{}`$.

Initially, all the baryons (for which we adopt a density parameter $`\mathrm{\Omega }_b=0.02`$) are assigned to the hot phase. Then, for each galaxy mass $`m`$ (with circular velocity $`v`$ and radius $`r`$) we compute the baryons in the three phases as follows.

a) The cold phase. Its mass $`m_c`$ is increased from inside out by cooling processes: at each time step $`\mathrm{\Delta }t`$, we have then $`\mathrm{\Delta }m_c(v)=4\pi r_{cool}^2\rho _g\mathrm{\Delta }r_{cool}`$; for the gas density we adopt the form $`\rho _g1/(r^2+r_{core}^2)`$, with the value of $`r_{core}`$ determined by requiring that the density at the virial radius is the same that would have been obtained has no gas cooled (hence depending on the past history of $`m_c`$ and $`m_h`$ corresponding to each galaxy mass $`m`$ at time $`t`$). The cooling radius $`r_{cool}`$ at time $`t`$ is computed by equating the cooling time $`\tau _{cool}=(3/2)kT/\mu m_H\rho _g(r)\mathrm{\Lambda }(T)`$ to the current time $`t`$. The cooling function $`\mathrm{\Lambda }(T)`$ is taken from Sutherland & Dopita (1993) for a mixture of 77 % hydrogen and 23 % helium.

The gas which cools settles into a rotationally supported disk; following Mo, Mao & White (1998), its radius $`r_d(v)`$ and the rotation velocity $`v_d(v)`$ are related to the galaxy circular velocity assuming that the angular momentum of the baryons which settle into the disc is a fixed fractions $`j_d`$ of the total angular momentum $`J`$ of the galaxy dark halo. The latter is usually expressed in the adimensional form $`\lambda J/(|E|^{1/2}/Gm^{5/2})`$. We use the relations $`r_d(v,\lambda ,c,m_d,j_d)`$ and $`v_d(v,\lambda ,c,m_d,j_d)`$ given by the above authors for a fraction of cold gas $`m_d=m_c/m`$ in a Navarro, Frenk & White (1997) potential with concentration $`c(v)`$. In the following, such relations will be used only to compute the Tully-Fisher relation resulting from our model (see fig. 5 below) and the dust extinction to be applied to the galaxy integrated stellar emission (see eq. 7 and below). In such computations we will assume $`j_d=0.05`$ (as discussed by the above authors) and integrate over the full distribution of $`\lambda `$ (a lognorm expression with mean $`\overline{\lambda }=0.05`$, and dispersion $`\sigma _\lambda =0.5`$ in $`log(\lambda )`$, see Warren et al. 1992; Cole & Lacey 1996; Steinmetz & Bartelmann 1995) to obtain average values for $`v_d`$ and $`r_d`$. The disk sizes resulting from the above values of $`j_d`$ and $`m_d`$ are consistent with recent observations, as shown in Giallongo et al. (2000).

b) The stars: the mass of baryons turned into stars in the timestep $`\mathrm{\Delta }t`$ is assumed to be $`\mathrm{\Delta }m_{}=\mathrm{\Delta }tm_c/\tau _{}`$, with a star formation timescale $`\tau _{}=ϵ_{}^1\tau _d(v_d/200\mathrm{km}/\mathrm{s})^\alpha _{}`$ proportional to the dynamical time of the disk $`\tau _dr_d/v_d`$. The normalization $`ϵ_{}`$ and the exponent $`\alpha _{}`$ are free parameters. Our default choice for them is $`ϵ_{}=0.05`$ and $`\alpha _{}=1.5`$, as in the fiducial model of Cole et al. (2000). This allows us to directly compare ours with previous results, and to single out the effects of the dynamics on the LFs. Note that the star formation is linked only to the changes of the cold material following the dynamical history of the galaxies; as discussed in §6 we do not include possible starbursts associated with newly formed mergers.

c) The hot phase is replenished by part of the cool baryons, namely, those re-heated and ejected into the hot phase by SN feedback. This is the most uncertain among the baryonic processes included in SAMs, and is currently parametrized by the simple expression $`\mathrm{\Delta }m_h=\mathrm{\Delta }m_{}(v/v_h)^{\alpha _h}`$ relating the amount $`\mathrm{\Delta }m_h`$ of the mass re- heated in a time step $`\mathrm{\Delta }t`$ to the mass of stars $`\mathrm{\Delta }m_{}`$ formed. The exponent $`\alpha _h`$, a free parameter, models the increased feedback efficiency in galaxies with decreasing mass. The choice of $`\alpha _h`$ has a considerable impact on the faint end of the LF (larger values yielding flatter shapes), but its values are strongly constrained by the observed Tully-Fisher relation. In fact, increasing $`\alpha _h`$ produces fainter luminosities for small mass objects. This flattens the resulting LF at the faint end, but it also moves the predicted Tully-Fisher relation at small $`v`$ away from the observed range. Our default choice of the parameters is $`\alpha _h=2`$ and $`v_h=200`$ km $`\mathrm{s}^1`$, since these yield the best joint fit to both the local LF and the Tully- Fisher relation, as shown below.

After having computed $`\mathrm{\Delta }m_c`$, $`\mathrm{\Delta }m_{}`$ and $`\mathrm{\Delta }m_h`$ from the processes a), b), c), we compute the increments due to galaxy coalescence/aggregation. The probability $`p(v^{},t^{},v,t)`$ that a galaxy with circular velocity $`v^{}`$ at time $`t^{}`$ is included into a galaxy with velocity $`v`$ at time $`t`$ has been calculated in §3.4. Then, the average star content in a halo with circular velocity $`v`$ is updated over the time grid according to following equation:

$`m_{}(v,t)={\displaystyle \underset{i}{}}{\displaystyle _0^v}𝑑v^{}{\displaystyle \frac{N(v^{},t_i)}{N(v,t)}}p(v^{},t_i,v,t)\mathrm{\Delta }m_{}(v^{},t_i),`$ (6)
where the sum extends over the time steps $`t_i=t_0+i\mathrm{\Delta }t`$ ranging from the initial time $`t_0`$ to $`t`$. Analogous equations govern the increments of $`m_c`$ and $`m_h`$. Such a procedure is similar to that introduced by Frenk & White (1991); however, here the transition probability for the galaxies is computed on the basis of the dynamics described in §2.

From the latter equation, the corresponding integrated stellar emission $`S_\lambda (v,t)`$ at the wavelength $`\lambda `$ is computed by convolving with the spectral energy distribution $`\varphi _\lambda `$ obtained from population synthesis models:

$$S_\lambda (v,t)=_0^t𝑑t^{}\varphi _\lambda (tt^{})\dot{m}_{}(v,t^{}).$$
(7)
Note that the average star formation $`\dot{m}_{}(v,t^{})`$ of galaxies with circular velocity $`v`$ at $`t^{}`$ is that corresponding to the star mass computed in eq. (6). Substituting its expression into eq. (7) demonstrates that $`S_\lambda (v,t)`$ contains the integrated contributions of star formation in the smaller progenitor systems (with circular velocity $`v^{}`$ at times $`t^{}<t`$) which by the time $`t`$ have been included in the haloes with circular velocity $`v`$. The integrations over the time $`t^{}`$ and the velocity $`v^{}`$ entering eq. (7) account for the average building up of the stellar population in hierarchically growing galaxies, as described by Frenk & White (1991).

In the following we adopt $`\varphi _\lambda `$ taken from Bruzual & Charlot (1993), with a Salpeter IMF. The dust extinction affecting the above luminosities is computed assuming the dust optical depth to be proportional to the metallicity $`Z_{cold}`$ of the cold phase (computed assuming a constant effective yield and that the metals are re-ejected to the hot phase in the same proportion as the reheated gas $`\mathrm{\Delta }m_h`$) and to the disk surface density, so that for the $`V`$\- band $`\tau _Vm_cZ_{cold}/\pi r_d^2`$. The proportionality constant is taken as a free parameter chosen as to fit the bright end of the local LF (see below and fig. 5); this yields for the proportionality constant the value $`3.5M_{}^1\mathrm{pc}^2`$ when the stellar yield is such as to produce a solar metallicity for a $`v=220`$ km/s galaxy. To compute the extinction in the other bands different extinction curves will be considered, including the Milky Way (MW), the Small Magellanic Cloud (SMC), and the Calzetti extinctions, see Calzetti (1997).

Fig. 5 - The Tully-Fisher relation relating the disk rotation velocity to the luminosity of spiral galaxies. Left panel refers to the our fiducial case $`\alpha _h=2`$, while the right panel to the high-feedback case $`\alpha _h=5`$. The shaded areas represent the region of the $`M_Iv`$ plane allowed by the observations (Mathewson et al. 1992; Willick et al. 1996; Giovanelli et al. 1997). The solid curves represent the model predictions assuming the disk rotation velocity equal to the DM circular velocity $`v`$, while the dashed curves are computed adopting the disk circular velocity $`v_d(v)`$ as a measure of the disk rotation velocity.

The above equation (7) relates the observable properties of galaxies (specifically their integrated stellar emission $`S_\lambda `$) to their dynamical properties. The luminosity-circular velocity relation constitutes a first, key prediction that can be tested against the observed Tully-Fisher relation. The comparison between the outcomes of our model (with our default choices for the free parameters) and the observations is shown in the top panel of fig. 5. The agreement is satisfactory, although the predicted magnitudes at given $`v`$ are slighly fainter than the data when the disk circular velocity $`v_d(v)`$ is used as a measure of the disk rotation velocity, as appropriate.

To show that increasing $`\alpha _h`$ is not a viable way to get flatter LFs (as anticipated above at point c), we also show the Tully-Fisher relation with $`\alpha _h=5`$, the value initially adopted by Cole et al. (1994) to get a flat local LF; in this case the predictions at faint luminosities fail to match the observations. Other

## 5. THE EVOLUTION OF THE GALAXY LUMINOSITY FUNCTION

From the galaxy velocity distribution $`N(v,t)`$ derived in §3 (see fig. 3) and the $`S_\lambda v`$ relation discussed above we derive the galaxy LF predicted by our model. The results are compared with the observational LF obtained from deep surveys. In particular, our predictions at high $`z`$ are compared with the observational LFs derived by Poli et al. (2001) using photometric redshifts for galaxies in the ESO New Technology Telescope deep field and in the Hubble Deep Field North and South, down to the limiting magnitudes $`I_{AB}=27.2`$. The LFs have been computed in the rest-frame B band for $`0<z<1`$, and at the rest-frame wavelength of 1700 Å for higher $`z`$. We stress that in the magnitude range where the photometric data overlap with those from spectroscopic surveys (see, e.g., Steidel et al. 1999), the photometric and spectroscopic LFs agree to a remarkable degree of precision.

### 5.1. Comparison with data: the effect of satellite aggregations

The comparison between the observed and the predicted LF is given in figs. 6 and 7. We compare with the data both the LFs correponding to the standard model (where only coalescence driven by dynamical friction is present) and those from the complete dynamics including binary aggregations between satellite galaxies; hereafter we shall refer to such cases as DF and DF+BA. In both cases, we adopt our reference set of cosmological parameters (given in §3.4) for a Universe dominated by a cosmological constant.

Fig. 6 - The local LF computed with the DF process alone (left panel) and with the inclusion of satellite aggregations (DF+BA, right panel). The different curves (almost overlapping here) refer to different choices of the dust extinction curve; MW (dashed line), SMC (dot-long dash), and Calzetti (dotted). The data are taken from Zucca et al. (1997, filled squares) and from the 2dFGRS survey (Madgwick et al. 2001, open circles).

The low-$`z`$ LFs presented in fig. 6 show how satellite aggregations flatten the LF at faint/intermediate luminosities. Such an effect is purely dynamical, so that the agreement of the model with the Tully-Fisher relation shown in fig. 5 is unchanged. At the very faint end, note the upturn of the predicted local LF in the DF+BA case, as expected from our discussion in §3.4. The local trend toward a flatter LF shown by the DF+BA model persists at $`z1`$ (see fig. 7), where our model actually provides a better fit to the shape of the data distribution when compared to the DF case.

The high-$`z`$ LFs are compared with data in the bottom panels of fig. 7. Note that all models overpredict the number of faint galaxies as noticed in Poli et al. (2001). Several processes not yet properly inserted in the SAMs can be at the origin of the discrepancy. As discussed in sect. 4, simple variations of the stellar feedback prescritions adopted in the model are not a viable solution, being constrained by the observed Tully-Fisher relation. On the other hand, our treatement of binary aggregations considerably reduces the discrepancy, as shown by the right-hand panels in fig. 7. A residual excess of the predicted over the observed LF remains at the faint end.

Fig. 7 - The LF in the B-band at $`z=0.7`$ (top panels) and at 1700 $`\AA `$ at $`z=3`$ (bottom panels) for the DF process alone (left panels) and with the inclusion of satellite aggregations (DF+BA, right panels). The solid curves refer to the LF without dust extinction, while the other curves refer to different choices of the dust extinction curve: MW (dashed line), SMC (dot-long dash). and Calzetti (dotted). The data with the error bars correspond to the LF computed from photometric redshifts by Poli et al. (2001); the data corresponding to the spectroscopic survey by Steidel et al. (1999) are also shown by the circles in the bottom panels. In the right-bottom panel we also show the effect of adopting the Sheth & Tormen (1999) instead of the Press & Schechter expression for the halo mass distribution; the corresponding galaxy LFs (derived including DF+BA) are represented as thin curves (the line types correspond to the different extinction curves as above).

A possible origin for it can be found in the statistics of the DM haloes (the Press & Schechter mass distribution) hosting the galaxies. To investigate the issue we computed the LF from our model (including the binary aggregations of satellite galaxies) adopting the Sheth & Tormen (1999) form (recalled in the Appendix A) for the mass function of the DM haloes hosting the galaxies; this is widely held to provide a better description of the halo statistics at high redshifts compared to the canonical Press & Schechter distribution.

The result (see right-bottom panel in fig. 7) is to further reduce the excess to a factor $`2.53`$ for luminosities fainter than $`M_{1700}19`$. At present, two explanations may be offered to account for the residual excess.

First, incompleteness in the data. This has been estimated by Vanzella et al. (2001) in the HDF-S data; at the faintest limits used to compute the LF ($`I_{AB}=27.25`$), it has been found to be close to 1.6 for extended sources. Even applying this correction factor to the whole faintest bin of the observed LF, the excess still remains.

Second, additional sources of feedback; e.g., the feedback due to the photoionization of the intergalactic medium by photons escaping with high efficiency from stars and quasars (Benson et al. 2001) could reconcile the predictions with the observed LFs. In this context, the contribution of binary aggregations reduces the amount of feedback required to yield flat LFs as to match the observations.

## 6. DISCUSSION, CONCLUSIONS, AND PERSPECTIVES

We have developed a semi-analytic model (SAM) for galaxy evolution, which includes the effects of binary aggregations between satellite galaxies orbiting inside common DM haloes (see fig. 1). The corresponding evolution of the galaxy mass function (see the velocity functions in fig. 3) has been calculated in a $`\mathrm{\Lambda }`$-CDM cosmology by solving in each host halo a kinetic Smoluchowski equation grafted onto our SAM code.

The main effects of aggregations on the galaxy mass distribution are summarized as follows:

i) Starting at $`z3`$ (see fig. 4), aggregations of satellite galaxies gradually deplete the number of low/intermediate mass galaxies, and flatten the slope of the mass function (see fig. 3). The aggregation cross section strongly depends on the ratio $`v/V_{rel}`$ between the galaxy circular velocity and the velocity dispersion of the halo wherein it orbits, so the increase of the average $`V_{rel}`$ following from the hierarchical clustering results in a progressive shift of the range of circular velocities where the aggregations flatten the mass distribution. At $`z=0`$ the resulting circular velocity distribution is flatter for small/intermediate circular velocities ($`50v200`$ km/s, see fig. 3). The result agrees with the outcomes from N-body simulations.

ii) The total (dynamical friction plus binary aggregations) net merging rate peaks at $`z3.5`$ and declines sharply afterwards, in agreement with N-body results and consistent with observations. The contribution of satellite aggregations to the destruction rate (which mainly affects the low mass end of the mass distribution) becomes comparable to that produced by dynamical friction for $`z<1`$, thus sustaining the galaxy destruction rate at low redhsifts. We have checked that the construction and destruction rates corresponding to the dynamics in our model (see fig. 4) agree with those obtained from N-body simulations within their limited mass resolution.

The flattening of the mass distribution at small/intermediate masses velocities $`50v150`$ km/s) affects the galaxy luminosity functions as shown in figs. 6, 7. These, obtained on adopting for the stellar section the standard prescriptions of SAMs, have been compared with data from deep surveys with spectroscopic or photometric redshifts. The results concerning the luminosity functions are summarized as follows:

$``$ at low redshifts, the binary aggregations flatten the LF at small/intermediate luminosities, producing also an upturn at very faint magnitudes $`M_B>16`$. At magnitudes brighter than $`L_{}`$ the LFs are not appreciably affected.

$``$ at higher redshifts, the excess of the standard SAM predictions over the observed LFs at the faint end (by a factor $`8`$ at $`z3`$ for $`M_{1700}18`$, see Poli et al. 2001) is considerably reduced by the introduction of satellite aggregations computed in our present model. At $`z3`$, our LFs computed adopting the Sheth & Tormen (1999) mass distribution for the host DM haloes further reduces the excess by a factor $`3`$ in the faintest bins ($`19.5M_{1700}18`$).

The origin of the residual excess at high redshifts requires further investigation. Incompleteness at the faintest limits of the observed LF ($`I_{AB}=27.25`$) has been found to be around 1.6 for point and extended sources (Vanzella et al 2001). With regard to curing the residual excess, we note that changing the parameters in the canonical stellar feedback prescriptions is not a viable solution, since a larger feedback parameter (as required to flatten the galaxy LF still more) would result in a Tully-Fisher relation that is too faint when compared to the data. An alternative may be provided by sources of feedback other than SNe, e.g., the feedback due to photoionization of the intergalactic medium by stars and quasars (as advocated by Benson et al. 2001) could reconcile the predicted with the observed LFs. In this context, the contribution of binary aggregations to the flattening of the LF will alleviate the need for extreme feedback efficiencies. Further insight could be provided by the comparison of the model predictions with K-band luminosity functions and with the baryonic mass function of galaxies (as suggested by Salucci & Persic 1999), a point we will address elsewhere.

On the other hand, our conclusions about the faint end of the galaxy LFs are little affected by other processes recently implemented in the SAMs, namely, the additional starbursts associated with the mergers following satellite aggregations (Somerville & Primack 1999), and the effect of tidal stripping of galactic subhaloes inside the host DM haloes (Benson et al. 2001). The former authors associate starbursts (with tunable timescale and amplitude) with each newly formed merger resulting from binary aggregations. We did not implement them in our model in order to single out dynamical effects of our description without introducing new (and uncertain) free parameters in the model. On the other hand, the starbursts associated with the mergers resulting from the galaxy encounters are expected to have a only a moderate inpact on the LF at faint/intermediate luminosities. In fact, the dynamics described by eqs. (3) and (4) is such that the aggregation events mainly occur between intermediate and small mass galaxies (the first being favoured by their larger cross section and the latter by their large number), to form larger units; thus, the destruction term dominates at the small/intermediate masses while the construction term dominates at large masses. The starbursts associated with mergers are then expected to brighten mainly the bright end of the LF.

As for tidal stripping, we note that Benson et al. (2001) have shown that it plays a minor role in determining the shape of the global galaxy luminosity function. This process, if anything, affects mainly the mass distribution at the very low-mass end (circular velocities $`v2050`$ km/s), while satellite aggregations modelled here affect the mass distribution at larger scales ($`v50150`$ km/s).

We stress that our description of the dynamics of galaxies inside common DM haloes does not introduce new free parameters in the SAMs; so, it constitutes a step toward improving the physical treatement of the processes involved in galaxy formation and evolution, rather than merely enlarging the set of phenomenogical scaling laws and the associated list of free paremeters. A further step in such a direction will be constituted by a refined description of how central dominant galaxies form. In fact, although in the present paper the formation is treated in terms of dynamical friction to retain continuity with other works, in principle it could be included in the statistical description provided by the Smoluchowski kinetic equation. The process is due to the rare slow encounters for which the focussing term $`v/V_{rel}`$ contained in eq. (4) dominates. In this case, the $`V_{rel}`$-dependent cross section (i.e., eq. 4 without the velocity average) would become more than linear in mass in the low velocity tail of $`g(V_{rel})`$ (in particular, $`\mathrm{\Sigma }m^{4/3}`$ holds in such a case). In these conditions the aggregation dynamics described by eq. (3) are marked by a phase transition of gravitational nature which breaks the system of colliding galaxies, originally comprised of a number of comparable objects, into two phases: a prompt, dominant merger; and a number of satellites of much smaller mass, which gradually aggregate to the dominant object (see Cavaliere, Colafrancesco & Menci 1991; Menci, Colafrancesco & Biferale 1993). The transition takes place at time around $`\tau _{dyn}(m/M)^{5/6}`$ (the detailed expression is given by the above authors), which turns out to be comparable with the timescale actually adopted in the present paper, see eq. (2).

Thus, in principle, the entire dynamics of galaxies in common haloes, including satellite aggregations and formation of central dominant mergers, could be treated within SAMs using only the kinetic description given in §3.2. In practice, however, this would require the detailed consideration of the galaxy velocity distribution to describe properly the probability for the occurence of slow, resonant encounters leading to the phase transition; this goes beyond the average description provided by eq. (4). In the present paper the formation of a central merger is described only in terms of dynamical friction; this can be viewed as a mean field, average rendition of the critical phenomenon leading to the phase transition. In perspective, we plan to describe all galaxy interactions in terms of the kinetic theory; this will be the subject of a future paper.

We gratefully acknowledge constructive advice of our referee toward correcting our manuscript and improving our presentation. Work supported by partial grants from ASI and MIUR.

## APPENDIX A

The statistics of the virialized DM condensations in the framework of the hierarchical clustering is commonly described in terms of the Extended Press & Schechter Theory (EPST; Bower 1991; Bond et al. 1991; Lacey & Cole 1993). Here we summarize the results used in the text.

The number density of virialized structures of mass $`M`$ at the cosmic time $`t`$ is given by the Press & Shechter (1974) formula

$$N_H(M,t)=\sqrt{\frac{2}{\pi }}\frac{\rho _0}{M^2}\left|\frac{dln\sigma }{dlnM}\right|\frac{\delta _c(t)}{\sigma (M)}e^{\frac{\delta _c(t)^2}{2\sigma ^2(M)}},$$
(A-1)
where $`\rho _0`$ is the cosmic average matter density, $`\sigma (M)`$ is the rms density fluctuations of the linear perturbation field in spheres of mass $`M`$ (see, e.g., Lacey & Cole 1993), and $`\delta _c(t)`$ is the density threshold for collapse. For $`\sigma (M)`$ we adopt the form corresponding to Cold Dark Matter (CDM) density perturbations in a flat Universe dominated by a cosmological constant $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$. The threshold $`\delta _c(t)`$ corresponds to the value – extrapolated to the present using the linear growth factor $`D(t,\mathrm{\Omega }_0,\mathrm{\Omega }_\mathrm{\Lambda })`$ for the density perturbations – of the overdensity of a homogeneous sphere at the point where the exact non-linear theory predicts collapse to a singularity. Its normalization $`\delta _{c0}`$ at $`z=0`$ and its time evolution $`D^1(t)`$ depend on cosmology; in a critical $`\mathrm{\Omega }=1`$ Universe $`\delta _{c0}=1.686`$ and $`D^1(z)=(1+z)`$ hold; the corresponding expressions for a flat Universe with a nonzero cosmological constant are given by, e.g., Eke, Cole & Frenk (1996).

Since the variance of the density field is an inverse function of $`M`$ and the density threshold $`\delta _c(t)`$ lowers with time, larger and larger overdensities collapse to form structures of increasing mass (hierarchical clustering). For a given DM mass $`M`$ it is possible to compute the merging rates and the progenitor distributions.

In particular, the probability that a condensation of mass $`M`$ present at time $`t`$ form a halo of mass between $`M^{}`$ and $`M^{}+dM^{}`$ at time $`t+dt`$ is given by the following expression:

$$\frac{d^2P_H(MM^{},t)}{dM^{}dt}=\frac{1}{\sqrt{2\pi }}\left[\frac{\sigma ^2}{\sigma ^2(\sigma ^2\sigma ^2)}\right]^{3/2}e^{\frac{\delta _c^2(t)(\sigma ^2\sigma ^2)}{2\sigma ^2\sigma ^2}}\left|\frac{d\sigma ^2}{dM^{}}\right|\left|\frac{d\delta _c(t)}{dt}\right|,$$
(A-2)
where $`\sigma `$ and $`\sigma ^{}`$ are the rms density fluctuations corresponding to the masse $`M`$ and $`M^{}`$, respectively.

The history of the previous mergers experienced by a mass $`M`$ is instead described by the probability distribution that a given mass $`M`$ at time $`t`$ has a progenitor of mass $`M^{}`$ at time $`t^{}<t`$:

$$\frac{dP_H}{dM^{}}(M^{},t^{}M,t)=\frac{\delta _c(t^{})\delta _c(t)}{(2\pi )^{1/2}(\sigma ^2\sigma ^2)^{3/2}}\frac{M}{M^{}}\left|\frac{d\sigma ^2}{dM^{}}\right|exp\{\frac{[\delta _c(t^{})\delta _c(t)]^2}{2(\sigma ^2\sigma ^2)}\}.$$
(A-3)
The survival time of a halo of mass $`M`$ is defined as the cosmic time at which its mass has grown to $`qM`$, and hence depends on the choice of the “mass step” of the merging tree. For a generic step parameter $`q`$, the probability $`prob\left[\tau _l<\tau \right]`$ that a halo with mass $`M`$ at time $`t`$ has a survival time $`t+\tau _l`$ with $`\tau _l<\tau `$ reads (Lacey & Cole 1993):

$`prob\left[\tau _l<\tau \right]`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \frac{[\delta _c(t)2\delta _c(t+\tau )]}{\delta _c(t)}}e^{\frac{2\delta _c(t+\tau )[\delta _c(t)\delta _c(t+\tau )]}{\sigma ^2(M)}}`$ (A-4a)
$`\times `$ $`[1erf(X)]+{\displaystyle \frac{1}{2}}[1erf(Y)]`$
$`X`$ $``$ $`{\displaystyle \frac{\sigma ^2(qM)[\delta _c(t)2\delta _c(t+\tau )]+\sigma ^2(M)\delta _c(t+\tau )}{\{2\sigma ^2(M)\sigma ^2(qM)[\sigma ^2(M)\sigma ^2(qM)]\}^{1/2}}}`$ (A-4b)
$`Y`$ $``$ $`{\displaystyle \frac{\sigma ^2(M)\delta _c(t+\tau )\sigma ^2(qM)\delta _c(t)}{\{2\sigma ^2(M)\sigma ^2(qM)[\sigma ^2(M)\sigma ^2(qM)]\}^{1/2}}}`$ (A-4c)
To make contact with previous works (Cole et al. 1994; Cole et al. 2000), our default choice for the mass step is $`q=2`$. We have verified that the SAM results are robust to smaller values down to $`q=1.2`$, as discussed by Cole et al. (2000).

Finally, we recall the expression for the mass function recently proposed by Sheth & Tormen (1999; see also Sheth, Mo & Tormen 2001) to provide a better fit to the N-body results. This has the form

$$N_H(M,t)=A\sqrt{\frac{2a}{\pi }}\frac{\rho _0}{M^2}\left|\frac{dln\sigma }{dlnM}\right|\left[1+\left(\frac{\sigma ^2}{a\delta _c^2(t)}\right)^p\right]\frac{\delta _c(t)}{\sigma (M)}e^{\frac{a\delta _c^2(t)}{2\sigma ^2(M)}}.$$
(A-5)
The values of the fitting coefficients $`A=0.3222`$, $`a=0.707`$ and $`p=0.3`$ are obtained by the above authors from comparison with the N-body results. The agreement of the above expression for $`N_H(M)`$ with the simulations has been confirmed by an independent analysis (Jenkins et al. 2001).

## APPENDIX B

We derive the average value for the tidal radius $`r_{tid}`$ of galaxies orbiting inside a host halo with radius $`R`$ and circular velocity $`V`$. Such a radius is that appropriate for a galactic subhalo which survives the tidal stripping of the host halo, a condition which requires the density of the galactic subhalo within $`r_{tid}`$ to exceed the density of the host halo interior to the pericentre $`r_p`$ of its orbit.

To this end we adopt the approach of Ghigna et al. (1998) who showed how the condition for survival against tidal stripping translates approximately into $`r_{tid}r_p(v/V)`$. The above expression for $`r_{tid}`$ has been tested by the above authors against high resolution N-body simulations, and has been proved to agree with the values of $`r_{tid}`$ measured in the simulations for all subhaloes except the few on very eccentric orbits (the latter have measured radii larger than expected partly due to the formation of tidal tails).

We average the above relation for $`r_{tid}`$ over the distribution of pericentres obtained by Ghigna et al. (1998) from N-body simulations, which we fit with a modified lognormal expression (in the variable $`(r_p/R0.08)`$ with logarithmic mean -1.3 and logarithmic variance 0.6). Performing the average yields for $`r_{tid}`$ the following expression:

$$r_{tid}=R(v/V)_{r_{cut}}^1p(r_p/R)d(r_p/R)(r_p/R)\alpha Rv/V,$$
(B-1)
where the lower limit $`r_{cut}`$ corresponds to the minimum tidal radius that a galactic subhalo can have without being severely distorted or disrupted. Following Bullock, Kravstov and Weinberg (2000), we adopt for $`r_{cut}`$ the radius of the peak of the galaxy velocity profile. For a Navarro et al. (1997) circular velocity profile $`r_{cut}=2.16r_{200}/c`$ holds, where $`c`$ is the concentration parameter of the subhalo, and $`r_{200}`$ is the radius where the average density of the subhalo would equal $`200\rho _c`$. The equations relating $`c`$ and $`r_{200}`$ with the circular velocity $`v`$ of the subhalo are given by the above authors. The validity of the above value for $`r_{cut}`$ has been tested against N-body results by Bullock et al. (2000).

It is easy to recast the relation $`r_{tid}=\alpha Rv/V`$ in terms of density. Substituting the relation $`v/V=(m/M)^{1/2}(R/r_{tid})^{1/2}`$ one obtains $`(r_{tid}/R)^3=\alpha ^2m/M`$ so that the ratio of the subhalo to the host halo average densities is given by $`1/\alpha ^2`$. The average over the distribution of pericentres defining the value of $`\alpha `$ yields typical values for such a ratio ranging between 3 and 5, depending on the concentration parameter $`c`$ of the subhaloes.

Note that taking the tidal radius $`r_{tid}=r_pv/V`$ as the limiting radius of the subhaloes is fully appropriate only for subhaloes much smaller than their host halo. This is the case for the bulk of the population of satellite galaxies (the one mainly contributing to binary aggregations) for which our treatment constitutes a valid approximation.


format: text

text:
# 1 Introduction

## 1 Introduction

The Crab Nebula is a supernova remnant. The supernova exploded in 1054 AD, and was monitored by contemporary Chinese astronomers (see e.g., Sollerman, Kozma, & Lundqvist 2001). Today the nebula offers a spectacular view, with a tangled web of line emitting filaments confining an amorphous part ghostly shining in synchrotron light. The beauty of the Crab makes it repeatedly appear in PR pictures, even from a southern observatory like the VLT (Fig. 1).

At the heart of the nebula resides the energetic 33-ms Crab pulsar. This m$`{}_{V}{}^{}16`$ object actually powers the whole visible nebula. The Crab nebula and its pulsar are among the most studied objects in the sky. This astrophysical laboratory still holds many secrets about how supernovae explode and about how pulsars radiate and energize their surrounding nebulae.

A main theme for pulsar research has been to understand the emission mechanism for the non-thermal pulsar radiation. This is still to be accomplished. No comprehensive model exists that can explain all the observed features of the radiation. Observationally, only recently was a broad range UV-optical spectrum of the pulsar published (Sollerman et al. 2000). We have now extended the study into the infrared (IR).

But even if most of the research on the Crab pulsar has concerned the radiation mechanism, almost all of the spin-down energy actually comes out in the particle wind. This is the power source of the Crab nebula. The stunning image of the pulsar environment obtained with CHANDRA (Fig. 2) captures a glimpse of the energetic processes at work. Direct evidence of the pulsar activity have long been seen in the system of moving synchrotron wisps close to the pulsar itself.

The detailed study of the wisps was started by Scargle (1969) using observations obtained before the authors of this article were born. Hester et al. (1995) used the HST to study the wisps at higher resolution, and presented their observations as the spectacular ’The Crab Movie’111http://oposite.stsci.edu/pubinfo/pr/1996/22.html, where the constant activity in the region around the pulsar is highlighted.

The most stunning discovery in these HST images was the knot sitting just 0.6 arcseconds from the pulsar. At 2 kpc this amounts to a projected distance of only 1000 AU. Hester et al. interpreted this feature as a shock in the pulsar polar wind. Our IR observations also allowed us to have a look at these manifestations of the magnetic relativistic wind from the pulsar.

## 2 IR photometry, reductions and results

IR imaging in the short wavelength (SW) mode of ISAAC was obtained in service mode on VLT on October 13, 2000. The exposures were taken in $`J`$s, $`H`$ and $`K`$s with a total exposure time of 156 seconds per band. The main goal of these short exposures was to properly calibrate our IR spectroscopy. However, the image quality provided by Paranal also allowed a detailed view of the central region of the Crab nebula. The near-IR images are displayed next to each other in Figure 3.

Photometry was obtained of the Crab pulsar and some of the stars in the field using PSF-fitting (DAOPHOT). We estimate that our magnitude measurements of the pulsar are correct to about 0.05 magnitudes.

In Figure 4 we plot our measurements as de-reddened fluxes together with the optical-UV data from Sollerman et al. (2000). We note that our IR fluxes deviate significantly from the most recent results published by Eikenberry et al. (1997). Our measurements give a fainter pulsar by some 0.3 magnitudes.

The pulsar magnitudes of Eikenberry et al. agree with those in the 2MASS point source catalogue. As our relatively isolated standards in the field show good agreement with 2MASS, we do not think the difference in pulsar magnitudes is due to an offset in the zero-point. Instead, our measurements of the Crab pulsar agree well with previous time-resolved photometry of the Crab (Penny 1982; Ransom et al. 1994). Such measurements generally use a large aperture and simply assume any non-varying contribution to be due to the background. By integrating under the pulsar light curve, they measure only the pulsating contribution of the flux. Our ISAAC photometry has excellent signal and image quality. In the complex region around the pulsar, this significantly improves the background subtraction. In particular, PSF subtraction excludes contributions from the wisps and the nearby knot. We believe that the difference with 2MASS is simply a matter of resolution, and that we are now able to subtract virtually all background from the pulsar.

The 3-color image of the central parts of the Crab (Fig. 5) is color coded with $`J`$s=blue, $`H`$=green and $`K`$s=red. In an attempt to keep some physical information in this image, the individual frames were scaled to make an object with a flat de-reddened Fν spectrum appear white.

This image shows the well known features of the inner Crab. The wisps are clearly visible in higher detail than previously obtained in the IR. Some filaments are also seen, most strongly in the $`J`$s band. This is most likely due to the \[Fe II\] $`1.26\mu `$m emission line. The $`K`$s band is instead dominated by amorphous synchrotron emission (Fig. 3). With suitable cuts the pulsar image appears slightly elongated. This is due to the presence of the knot first identified by Hester et al. (1995) on a HST image. To reveal this structure in our images we constructed and subtracted a PSF from the stellar images. After subtraction, the knot is clearly revealed in all three bands. A color image made out of the PSF subtracted frames is shown in Figure 6. The image directly gives the impression that the knot is redder than for example the wisps.

Quantifying this we have estimated the magnitudes of the knot as well as of the nearby wisp 1. The knot was measured within an aperture of 0.9 arcseconds and the wisp was simply measured with an aperture of 1.2 arcseconds. The spectral energy distribution of the knot is shown in Figure 4. It is clearly red. In the $`K`$s-band the flux from the knot amounts to about $`8\%`$ of the flux of the pulsar. The stationary wisp appears to have a flatter spectrum in this regime.

## 3 Optical data from VLT and HST

To extend the wavelength region over which to derive the spectral characteristics of the inner Crab components, several optical images are available in the ESO and HST data archives.

The VLT PR image (Fig. 1) was taken with FORS2 on 10 November 1999, just two weeks after first light. On the five minute $`B`$-band and the one minute $`R`$-band exposures, we could PSF-subtract the pulsar to reveal the knot. The knot appears clearly in both frames, and amounts to a few percent of the pulsar light at these wavelengths.

Data on the Crab pulsar are also available in the HST archive. Most of the observations are from Jeff Hester’s comprehensive monitoring programme of the inner parts of the nebula, which has shown just how active the Crab nebula really is. These frames allow a detailed study of both the spectral and temporal properties of the knot. In August 1995 the region was observed in three filters (F300W, F574M, F814W). We measured the de-reddened spectral index for the knot to be $`\alpha _\nu 0.8`$, which agrees well with our IR data (Fig. 4).

Furthermore, a wealth of data in the F574M filter allows a study of the temporal behaviour (Figure 7). First of all we note that the knot is indeed present in all frames. It thus appears quasistationary for more than six years, although the position appears to vary at the 0.1 arcsecond level. The de-reddened flux of the knot (Fig. 4) is measured to be $`9\times 10^{28}`$ ergs s-1 cm-2 Hz-1 but variations of the flux by at least $`50\%`$ are observed.

## 4 Discussion and Implications

For the pulsar itself we have added new information in the IR. Together with the optical-UV data in Sollerman et al. (2000) this significantly revises the observational basis for the pulsar emission mechanism. In fact, most of the theoretical efforts have been based on the old optical data from Oke (1969) and the IR continuation of Middleditch, Pennypacker, & Burns (1983). Our new results call for a fresh look on the emission mechanism scenarios for young pulsars.

For the knot, we have shown that the structure is indeed quasi-stationary, and that the emission has a red spectrum. Few models are available for the knot. Lou (1998) presented a formation scenario in terms of MHD theory, while Shapakidze & Machabeli (1999) argue for a plasma mechanism. None of these scenarios predict a very red spectral distribution.

Another area where caution may be required is in the recent claims of weak and red off-pulse emission from the Crab pulsar in the visible (Golden, Shearer, & Beskin 1999). It is clear that the knot close to the pulsar has to be seriously considered in these kinds of studies.

## 5 Future plans

The Crab pulsar and its environment continue to be the prime astrophysical laboratory for the study of the pulsar emission mechanism and the spin-down powering of pulsar nebulae. Although much observational effort has been put into this object, a modern re-investigation is likely to clean up the many contradictory measurements. Optical imaging in good seeing would require only a few minutes with VLT, and would directly determine the knot-subtracted spectral energy distribution. ISAAC in the LW range can in less than 3 hours clarify if the knot is significantly contributing to the emission at these frequencies, and would establish if the IR drop of the pulsar is real.

Most interesting is the possibility to monitor the very central parts of the pulsar environment with NAOS/CONICA. With a resolution superseding HST we will be able to monitor the structures close to the pulsar, with 2 pixels corresponding to merely 50 AU. This would provide an unique opportunity to study the structure and dynamics of the inner pulsar wind and its interaction with the surroundings.

## References

* Eikenberry, S. S., Fazio, G. G., Ransom, S. M., Middleditch, J., Kristian, J., & Pennypacker, C. R. 1997, ApJ, 477, 465
* Golden, A., Shearer, A., & Beskin, G. M. 2000, ApJ, 535, 373
* Hester, J.J, Scowen, P. A., Sankrit, R., et al. 1995, ApJ, 448, 240
* Lou, Y.-Q. 1998, MNRAS, 294, 443
* Middleditch, J., Pennypacker, C., & Burns, M. S. 1983, ApJ, 273, 261
* Oke., J. B. 1969, ApJ, 158, 90
* Penny, A. J. 1982, MNRAS, 198, 773
* Ransom, S. M., Fazio, G. G., Eikenberry, S. S., Middleditch, J., Kristian, J., Hays, K., & Pennypacker, C. R. 1994, ApJ, 431, 43
* Scargle, J. D. 1969, ApJ, 156, 401
* Shapakidze, D. & Machabeli, M. 1999, ptep.proc, 371
* Sollerman, J., Lundqvist, P., Lindler, D., et al. 2000, ApJ, 537, 861
* Sollerman, J., Kozma, C., & Lundqvist, P. 2001, A&A, 366, 197


text:
# Black Holes, Mergers, and the Entropy Budget of the Universe

## Abstract

Vast amounts of entropy are produced in black hole formation, and the amount of entropy stored in supermassive black holes at the centers of galaxies is now much greater than the entropy free in the rest of the universe. Either mergers involved in forming supermassive black holes are rare, or the holes must be very efficient at capturing nearly all the entropy generated in the process. We argue that this information can be used to constrain supermassive black hole production, and may eventually provide a check on numerical results for mergers involving black holes.

In an ideal gas of $`N`$ identical particles, the entropy $`S`$ is an extensive quantity proportional to $`N`$. Combining two ideal gases in equilibrium at equal temperature and pressure with $`N_1`$ and $`N_2`$ particles respectively, we find

$$S_{tot}N_1+N_2,$$
(1)
and since $`\mathrm{\Delta }S=S_{tot}S_1S_2=0`$, the process is reversible. In the case of black hole (BH) thermodynamics, a black hole of Schwarzschild radius $`R=2M`$ has entropy $`S=\frac{1}{4}A`$, where the area of the horizon is $`A=4\pi R^2=16\pi M^2`$. Now combining two black holes, we find (assuming no energy is lost in the process) $`R=2(M_1+M_2)`$, and $`S_{tot}=4\pi (M_1+M_2)^2,`$ so now $`\mathrm{\Delta }S=S_{tot}S_1S_2=8\pi M_1M_2`$ and this process is therefore irreversible. $`I.e`$., entropy is produced.

On the other hand, if the entropy of the new black hole is to be the same as the sum of the two initial black holes, energy must necessarily be shaken off as they coalesce. The lower bound for this is given by $`S_{tot}=S_1+S_2`$, or $`M_{tot}^2=M_1^2+M_2{}_{}{}^{2},`$ and for $`M=M_1=M_2`$ we have $`M_{tot}=\sqrt{2}M`$ and $`\mathrm{\Delta }E=2M\sqrt{2}M0.59M.`$

Considerable effort has been expended in studying the interactions of black holes with black holes, neutron stars, and other objects. While results have been obtained for the ideal case of a head-on collision between two Schwarzschild black holes, the more typical astrophysical mergers of a BH-BH binary, BH-neutron star pair or the capture of a stellar mass object in a galactic core, will take longer to unravel. By a typical merger, we mean the coalescence of two objects with different masses, angular momenta, and magnetic fields spiralling together, perhaps in the environment of accretion disks or their distorted remnants. A fuller understanding of such complicated but realistic situations will require extensive numerical analysis. In this Letter we consider the release of entropy into the inter-galactic medium in the production of large black holes by mergers involving black holes. The result can be a useful check on such numerical analyses. It can also have bearings on structure formation processes for galaxies, quasars, as well as for clusters.

The mass-energy that is thrown off in merger processes will contain some amount of entropy, and the entropy to mass ratio $`\rho _s^{out}`$ of this material will depend on the parameters of the merging system, including time ($`i.e`$., on the progress of the merger). Although the entropy of the final stable black hole will have an initial entropy to mass ratio $`\rho _s^{in}`$ linearly proportional to the black hole size as determined by black hole thermodynamics (we expect, for a Schwarzschild black hole, $`\rho _s^{in}=\frac{\frac{1}{4}Area}{M}=2\pi R`$), no such result holds for $`\rho _s^{out}(t)`$. However, the entropy to mass ratio of the total system $`\rho _s^{tot}=\rho _s^{in}+\rho _s^{out}`$ must be increasing in order to reach the Bekenstein value (and also because these processes are irreversible). Therefore we expect $`\rho _s^{out}`$ to be increasing with time if the region from which it is to be expelled contains an incoherent mixture of this material with the matter that will eventually fall into the black hole. While this perspective is clearly naive, we take it as our working hypothesis, and postpone its refinement to later work.

Let us first summarize the current state of knowledge concerning mergers. The mergers of black hole-neutron star, neutron star-neutron star, and black hole-black hole have all been studied in recent years, using various methods (numerical simulations, analytic estimates, and hydrodynamic simulations for neutron stars), in the Newtonian or post-Newtonian approximations, or treated fully relativistically when feasible. Compact binaries, containing black holes or neutron stars, lose energy in emitting gravitational radiation; slowly they spiral towards each other until they reach the innermost stable circular orbit when they start to coalesce and merge. In the merger of stellar mass black hole-neutron star binaries, a gas mass of a few $`10^1M_{}`$ is left in an accretion torus around the black hole and neutrinos are radiated. Subsequent annihilations of neutrinos and antineutrinos purportedly lead to gamma-ray bursts. But while neutrino emission from the remnant of an inspiraling binary neutron star following coalescence may be important for the cooling of the remnant, it is negligible for the emission of angular momentum of the merged objects; the consequent evolution of the remnant is dominated by the emission of gravitational waves. Estimates for the Galactic merger rates range from $`10^5yr^1`$ for neutron star mergers to an order of magnitude lower for black hole and neutron star mergers. On the stability of coalescing binary stars, one group claims that massive neutron stars, stable in isolation, individually collapse to black holes prior to merger, while another group claims that the tidal field from a binary companion stabilizes a star against gravitational collapse. Cosmological gamma-ray bursts are often thought to be associated with gravitational collapses of massive stars, but it has been suggested that the binary neutron star merger scenario is actually more favored than single stellar collapses. Head-on collisions of two equal mass, nonrotationing black holes for various initial configurations have been studied using several independent methods, with excellent agreement for total gravitational radiation between numerical results and the analytic estimates. The Binary Black Hole Grand Challenge Alliance has used a three dimensional numerical relativity code to study coalescing black hole binaries. We hope this Letter will be useful in providing a check on some of the future works related to these results.

Since we are dealing with general relativity, we cannot think of $`M`$ as proportional to the number of particles in the system. In fact, since all the particles in a black hole are hidden inside the horizon, this type of counting loses meaning, as, for instance, the number of baryons is the global quantum number, and gravity violates all global symmetries. For a discussion of the implications see .

It is interesting to compare black hole entropy with the entropy of the universe. Currently the entropy density is $`s_0=2970\left(\frac{T}{2.75K}\right)^3\frac{1}{cm^3}`$, so the entropy within our horizon today is approximately

$$S_0=s_0(cH_0^1)^32.35\times 10^{87}\left(\frac{Th^1}{2.75K}\right)^3$$
(2)
On the other hand, the entropy of a Schwarzschild black hole is

$$S_{BH}=4\pi \left(\frac{M_{}}{m_{pl}}\right)^2\left(\frac{M}{M_{}}\right)^21.05\times 10^{77}\left(\frac{M}{M_{}}\right)^2.$$
(3)
Hence, a couple of hundred thousand solar mass black holes can contain as much entropy as is free in the entire universe. There is increasing evidence that supermassive black holes ($`SMBH`$s) exist at the center of many galaxies, and that they are the sources which power active galactic nuclei and quasars. By now we know that a large fraction (at least about 30%) of galaxies contain such $`SMBH`$s with masses $`10^6M_{}M_{BH}10^9M_{}`$. (Observations by the Hubble space telescope suggest that, e.g., our own galaxy has a $`SMBH`$ with mass $`10^6M_{}`$ and M87 has one with mass $`10^9M_{}`$.) There are roughly $`10^{11}`$ galaxies in the universe. Thus the black hole entropy dominates all other sources of entropy. $`A`$ $`priori`$ this is not a problem, since as mentioned above black hole entropies are hidden behind horizons. However, the formation of these supermassive black holes does raise the question, ‘How is it possible for them to form without considerable loss of entropy to the environment as one would expect from stellar mass black holes merging at the cores of galaxies to form the supermassive black holes seen today?’

Let us give a very rough estimate of the total entropy stored in $`SMBH`$s, where we assume they are nonrotating. This gives,

$`S_{BH}^{tot}`$ $`=`$ $`{\displaystyle \underset{gal.=1}{\overset{N}{}}}4\pi {\displaystyle \frac{M_{BH}^2(gal.)}{m_{pl}^2}}`$ (4)
$``$ $`3.2\times 10^{101}\times \left({\displaystyle \frac{M}{10^7M_{}}}\right)^2\left({\displaystyle \frac{N}{10^{11}}}\right)`$ (5)
where N is the number of galaxies within our horizion, and $`M_{BH}^2(gal.)`$ is the distribution of masses of the $`SMBH`$s at the galactic cores. (For some galaxies this may be zero.) To arrive at the second line we have set all values of $`M_{BH}^2(gal.)`$ to $`M`$. Comparing with $`S_0`$ and choosing $`M=10^7M_{}`$ and $`N=3\times 10^{10}`$ implies (very conservatively) that less than one part in $``$1014 of the black hole entropy could have escaped during formation. Because there is actually a distribution of super heavy black hole masses, a more accurate calculation of their entropy would take that into account (and would be dominated by the heaviest $`BH`$s), though our order of magnitude estimate will be more than sufficient for the arguments given here.

Since big bang nucleosynthesis is highly dependent on the baryon-to-entropy ratio, and since the present ratio is in agreement with big bang nucleosynthesis ($`BBN)`$ calculations at the $`0.1MeV`$ scale, it appears that any significant entropy production between the time of $`BBN`$ and today can be ruled out. This means that in the formation of $`SMBH`$s, which presumably did not start until the structure formation era at $`z100`$, very little entropy could have escaped from the $`BH`$ forming regions. We require that in order for the change in visible entropy $`\mathrm{\Delta }S10S_0`$ from $`BBN`$ until today, the entropy released from $`BH`$ formation satisfy $`\mathrm{\Delta }S_{BH}^{tot}10^{88}`$ and so $`\mathrm{\Delta }S10^{13}S_{BH}^{tot}`$. Hence, the $`BH`$ formation must be extremely efficient at keeping entropy hidden behind the horizon. We believe the easiest way to do this is through nearly spherical collapse of large amounts of matter to directly form $`SMBH`$s, and not through mergers that are potentially much less efficient. There is indirect support for this belief: in the gravitational collapse of a rotating supermassive star, it was found recently that a supermassive black hole is formed coherently, with almost all of the matter falling into the hole, leaving very little ejected matter to form a disk, and strongly emitting gravitational waves. This is also in line with arguments given by Silk and Rees that $`SMBH`$s form in protogalactic cores as quasarlike objects before the epoch of peak galaxy formation, and argue against mergers as a primary component of $`SMBH`$ formation. The quasars then had profound effects on star formation. We conclude it is most likely that $`SMBH`$s formed in gravitational collapse events in protogalactic cores at or before galaxy formation. Alternatively, if $`SMBH`$s were to have formed via mergers, as seems less likely, then the mergers must have been extremely efficient at hiding all the entropy produced, and thus would certainly be more probable to have taken place in regions of low ambient density.

We can use the numerical results of Anninos $`etal`$.  to get an estimate of energy radiated and entropy produced in $`BH`$ mergers. For two equal-mass $`M`$, uncharged and non-rotating $`BH`$s that merge via a head-on collision, the total energy radiated in gravitons is $`E0.002M`$. This is in agreement with earlier analytic results which gave

$$E0.0104\frac{\mu ^2}{M}$$
(6)
where $`\mu =\frac{mM}{m+M}`$is the reduced mass for two $`BH`$s of mass $`m`$ and $`M`$, and so for $`m=M`$, gives $`E0.0026M`$. More mass is expected to be radiated in collisions with non-zero impact parameter and for rotating and for charged $`BH`$s, or if the $`BH`$s are surrounded by a medium or accretion disks. Hence we take this as an approximate lower bound.

For wavelength $`\lambda R=2M`$ the number of gravitons produced is $`N\frac{E}{\epsilon _\gamma }5\times 10^{72}\left(\frac{M}{M_{}}\right)^2`$. But, although an enormous number of gravitons are produced in $`BH`$ mergers, they would have a high degree of coherence, and so only a small amount of entropy could be released in the process. Where we expect a large entropy production is in mergers of $`BH`$s with neutron stars, and mergers that take place within dense media, $`e.g`$. involving dense protogalactic cores or overlapping and/or distorted accretion disks. Then much of the energy released could be in photons, with very little coherence due to local dissipation. The two merging $`BH`$s or neutron star and $`BH`$ will coalesce into a single but very perturbed black hole. The evolution into one black hole with a single horizon is still not well understood, and is the subject of ongoing numerical work. Hence the potential for entropy production from mergers is enormous and can, given sufficient numerical work, eventually be used to bound the number and types of mergers allowed since the time of the beginning of the growth of density perturbations $`z100`$. Indeed, as quasar formation peaked at around $`z3`$, and star formation peaked around $`z1.5`$, this also suggests fewer mergers and larger intitial $`BH`$s ($`QSO`$s). (In fact, $`SMBHSMBH`$ mergers may also be constrained by the fact that life exists on this planet, as accretion of dwarf galaxies and subsequent mergers would act as super gamma ray bursts. The resulting luminosity would be greater than $`L_{}`$ within $`3Mpc`$, and lasts for minutes to hours, so within the galaxy this could be fatal. However, while this thought is intriguing, it may be academic since it is known that even if galaxies merge, their $`SMBH`$s have a low probability of merging too, and orbits of several $`pc`$ are stable for times longer than the age of the Universe.)

Conclusions: Under some modest assumptions we have concluded that either mergers involving black holes that generate larger black holes are rare, or they must capture nearly all the entropy generated in the process. This has implications for galaxy formation and cosmological models. It can also serve as a guide and check on further numerical work on mergers. While our results are necessarily imprecise, the fact that less than one part in $`10^{13}`$ of the entropy stored in black holes has been released into the intergalactic environment is provocative.

We thank Chuck Evans for useful discussions and for pointing us to a number of references. This work was supported in part by the US Department of Energy under Grants No. DE-FG02-97ER-41036 and DE-FG05-85ER-40226, and by the Bahnson Fund of the University of North Carolina.