text:
# Bumpy Power Spectra and Galaxy Clusters

## 1 Introduction

The distribution of matter on large scales in the Universe is supposed to have evolved by gravitational interactions from seeds originating in quantum fluctuations that were stretched to cosmological dimensions by inflation. Ordinarily, this initial spectrum of fluctuations is assumed to be a featureless power law. This requires the least number of parameters and is produced in the simplest models of inflation. However, different observations, such as cluster redshift surveys (Einasto et al. 1997) and galaxy surveys (Hamilton & Tegmark 2000, Gaztanaga & Baugh 1998, Broadhurst et al. 1990) reveal possible traces of features in the matter power spectrum. In addition, a lower than expected second peak in the CMB spectrum as measured by BOOMERANG (Lange et al. 2001) appears to be incompatible with the simplest models of cold dark matter and a scale invariant primordial spectrum.

While the results, particularly from the surveys are still not beyond statistical doubt, a genuine feature(s) in the primordial spectrum cannot be ruled out. Several mechanisms have been proposed that could generate features in the primordial spectrum during the epoch of inflation. They commonly involve an extension of the simplest one field inflation model, e.g. by coupling the inflaton to a massive particle (Chung et al. 1999) or considering two-field inflation (Lesgourgues, Polarski & Starobinsky 1998). There are also more exotic ideas (see e.g. Martin, Riazuelo & Sakellariadou 1999). Another idea proposed is to Taylor expand the primordial power spectrum to include higher order terms that account for a running spectral index (Lidsey et al. 1997, Hannestad, Hansen & Villante 2000). The latter is able to introduce a very broad negative or positive bend into the power spectrum.

While possibly providing a motivation for observed features on a specific scale (i.e. scales of about 100 Mpc), the inflationary mechanisms proposed appear in principle capable of producing features on other scales as well. Here we examine the effect of primordial features on smaller scales $`k0.40.8h\mathrm{Mpc}^1`$, where due to (the onset of) non-linear evolution and the problem of biasing, a connection with the primordial spectrum is much harder to establish. We have previously examined the effects of such bumpy power spectra on the cosmic microwave background (Griffiths, Silk and Zaroubi 2001).

In Section 2 we will present the modifications applied to an otherwise unperturbed SCDM power spectrum. These spectra were then used as input to cosmological $`N`$-body simulations as described in Section 3. Section 4 deals with the complete analysis of these simulations with respect to the evolution of the power spectrum, velocity statistics, the masses of galaxy clusters, the halo-halo correlation function, and density profiles for a selection of halos. In Section 5 we try to link our numerical results to analytical prediction mainly based on the Press-Schechter theory (Press & Schechter 1974). We close with a discussion of our main results in Section 6.

## 2 The Power Spectra

We focus our attention on an SCDM model with parameters as given in Table 1. Although in this paper we are looking at properties of clusters and particularly their abundance evolution, we decided to use only COBE-normalised spectra with spectral index $`n=1`$ for modifications. This is because the features we add to the spectra are on scales that at the present time are subject to non-linear evolution and thus are expected to lead to a deviation from the $`\sigma _8\mathrm{\Omega }`$ relation (Eke et al. 1996) that would otherwise be used to cluster-normalise the spectra. We also focused our attention mostly on negative amplitude features – or dips – in the spectrum of SCDM models. In addition to the significance of features for structure evolution on relatively small scales, dips will reduce the overall power in the spectrum, pushing the normalisation and possibly other properties towards what is expected in OCDM models.

The corresponding SCDM power spectrum was calculated using the publicly available CMBFAST code (Seljak & Zaldarriaga 1996) and hence is COBE-normalised. For Model 1b we lowered the amplitude to reach a normalisation that agrees with the cluster abundance as described in Eke et al. (1996). All other models are based on Model 1a with no further modification to the amplitude than introduced by the artificially added feature. The modified power spectra follow the equation:

$$P_{\mathrm{mod}}(k)=P(k)(1\pm A\mathrm{exp}[0.4(\frac{\mathrm{log}k\mathrm{log}k_0}{\sigma _{\mathrm{mod}}})^2])$$
(1)
where P(k) is the unmodified spectrum.

These features are completely Gaussian in log-space and can be described by their width $`\sigma _{\mathrm{mod}}`$, height $`A`$, and the location $`k_0`$ of the bump (+) and dip (–), respectively. Their parameters are also summarised in Table 1. Each modification (characterised by the Model no. 2–6) was applied to the SCDM spectrum corresponding to Model 1a keeping the COBE normalisation fixed. A visual impression of the resulting power spectra can be found in Fig. 1 where all spectra are plotted (linearly extrapolated to redshift $`z=0.0`$). These spectra are now used as an input to our initial conditions generator for the cosmological $`N`$-body simulations.

Again, we decided to use a COBE-normalised SCDM power spectrum as reference model because our intention was to investigate the influence of artificially introduced features on scales corresponding to galaxy clusters (i.e. 8$`h^1\mathrm{Mpc}`$). For this reason we chose the COBE normalisation even though this SCDM model seems to be rather unattractive or even ruled out nowadays. However, by taking away power on cluster scales (as done for the majority of our ’feature’ models) the value of $`\sigma _8`$ drops as we do not apply any re-normalisation of the power spectra. The purpose was not to find a new standard model which fits the observational data better but rather to analyse the influence of such features on an interesting range of cluster quantities; and this is more easily followed in a more distinctive structure formation scenario.

## 3 The $`N`$-body Simulations

The simulations were carried out using Couchman’s AP3M code (Couchman 1991). All simulations were performed with $`128^3`$ particles in a box of side length 256$`h^1\mathrm{Mpc}`$; the (comoving) force resolution was fixed at 100$`h^1\mathrm{kpc}`$ for all runs. We evolved the particle distribution from redshift $`z=30.0`$ until $`z=0.0`$ in 5000 steps. The box size was chosen such that the scales on which the Gaussian features had been added lie comfortably within the range covered by the simulations, e.g. our simulations cover the k-range from $`k_{\mathrm{min}}=0.0245`$ (limit set by box size) to $`k_{\mathrm{max}}=1.571`$ (limit set by particle number for representing the waves initially present), and the modifications lie clearly within that range (cf. Fig. 1). The particle-particle summation part of the AP3M code guarantees that we properly follow the evolution of all initially present waves especially for models 2–4 (cf. Couchman 1991).

For identifying particle groups within our numerical simulations we used the standard friends-of-friends algorithm (Davis et al. 1985) with the linking-lengths $`ll=0.2`$ for the SCDM, and $`ll=0.17`$ for the OCDM model, respectively (cf. Knebe & Müller 1999).

## 4 Analysis

### 4.1 Power Spectrum Evolution

In Fig. 2 we plot the evolution of the dark matter power spectrum for models 1 through 6 (Model 1b and the OCDM model are left out for clarity as there is nothing unexpected to observe). We see that the spectrum appears to be converging towards the unperturbed spectrum at late times. The feature nearly vanishes completely, leaving no further imprint in the power spectrum. The thin solid line crossing the whole plot corresponds to the unmodified, COBE-normalised SCDM spectrum linearly extrapolated to $`z=0.0`$. Comparing it with the spectra derived from the numerical simulations, we notice a significant boost in power on small scales due to non-linear evolution. However, Model 6 with the very prominent dip at 16 $`h^1\mathrm{Mpc}`$ nearly matches the linearly extrapolated SCDM power spectrum: at late times the non-linear boost of power almost exactly compensates the lacking power at the location of the dip.

We may conclude that even when starting with a prominent feature on small scales, there might only be little if any evidence for it left in the present universe when looking at the dark matter power spectrum P(k). In any case, the evolution of the models has to be different, and one needs to think of other ways to investigate the influence of such features in the present day universe or to detect them in the evolution of P(k) and related quantities. As far as the spectrum itself is concerned current galaxy clustering surveys offer little prospect of detecting these features in the evolution of P(k). Observations of the Ly$`\alpha `$ forest, however, appear to be a promising tool for constraining P(k) at high redshift (see e.g. Croft et al. 2000, Weinberg et al. 1998), particularly on the scales under investigation here.

### 4.2 Velocity Statistics

Next we consider the distribution of velocities for a) the dark matter particles and b) friends-of-friends particle groups. The peculiar velocity field results from the gravitational acceleration that develops from initial density fluctuations in the early universe. Clusters of galaxies can therefore be used as tracers of the large-scale peculiar velocity field (Bahcall, Cen & Gramann 1994).

In Fig. 3 we show the probability distribution $`P(v)`$ of dark matter particles with peculiar velocities in the range $`v\pm dv`$ along with the integrated distribution $`P(>v)`$. The same is plotted for friends-of-friends groups identified in all runs in Fig. 4. All curves are normalised by the total number of particles and total number of particle groups in the respective model. Here we also show the data for SCDM Model 1b as well as the fiducial OCDM model. Even though there are differences of a factor of two at the high velocity end of the distribution, it might be difficult to observationally distinguish these differences between modified and unmodified models.

The deviations between the cluster-normalised SCDM model and the COBE-normalised one are bigger than any differences between our ’feature’ models. This deviation in velocities tending to lower values in the cluster-normalised model can be explained by the lower value for $`\sigma _8`$: the amplitudes of the initial density fluctuations are smaller and hence it takes longer to accelerate particles (and clusters) to high velocities; in the course of a simulation the velocity distribution function is always similar to a Maxwellian distribution whose peak gradually moves from low to high velocities. Apart from this we get increasingly more high velocity particles leading to a bigger ’tail’ in the distribution. However, the most prominent difference is in the high velocity tail of the distribution for clusters. Only here it is possible to discriminate between models as our dip models show a significant drop in the integrated probability. We do not observe any galaxy clusters with peculiar velocities higher than 3000 $`\mathrm{km}\mathrm{s}^1`$ whereas there is a distinctive number of these objects in the unmodified Model 1a (and the bump Model 2). A detailed check showed that this tail is mainly due to objects with masses $`M<10^{14}`$$`h^1\mathrm{M}_{}`$. When only taking into account particle groups with heavier than $`10^{14}`$$`h^1\mathrm{M}_{}`$all curves for the ’feature’ models fall on top of each other. However, the location and the width of the feature does not seem to have any significant influence.

### 4.3 The Masses of Galaxy Clusters

The most basic property of a galaxy cluster is its mass $`M`$. Nevertheless, this quantity can provide a lot of information especially when using the (cumulative) distribution of objects with a certain mass $`M`$; and the evolution of the abundance of massive clusters within a given mass range is indeed one of the corner stones of the currently favoured $`\mathrm{\Lambda }`$CDM model (Bahcall et al. 1999).

#### 4.3.1 Mass Function and Press-Schechter Prediction

In Fig. 5 we show the cumulative mass function of particle groups for all our models at a redshift of $`z=0.0`$. Even though we could not find well-pronounced imprints of the modifications in the power spectra at redshift $`z=0.0`$ (cf. Fig. 2), we clearly see differences in the amplitude and slope of the cumulative mass function $`n(>M)`$. These deviations are mainly at the low mass end of the resolvable mass range, where a positive feature leads to an excessive number of objects and a negative feature to a lack of groups. This is in general agreement with an excess/lack of objects on mass scales corresponding to where the features are located. However, we always end up with the same number of galaxy clusters for masses $`M>10^{15}`$$`h^1\mathrm{M}_{}`$(besides for Model 1b which shows too few massive groups due to the low normalisation). Unfortunately we are not able to resolve particle groups lighter than about $`210^{13}`$$`h^1\mathrm{M}_{}`$ and moreover, the simulation volume is still too small to get a statistically significant number of objects heavier than about $`310^{15}`$$`h^1\mathrm{M}_{}`$.

To check the validity of our results and get an impression of how they might generalise for larger/smaller mass objects we performed Press-Schechter (PS) calculations (Press & Schechter 1974) of the abundance of gravitationally bound objects. The (differential) number of objects for a given mass $`M`$ can be calculated using their formula:

$$\frac{dn}{dM}dM=\sqrt{\frac{2}{\pi }}\frac{\overline{\rho }}{M}\frac{\delta _c}{\sigma _M}\left|\frac{d\mathrm{ln}\sigma _M}{d\mathrm{ln}M}\right|\mathrm{exp}\left(\frac{\delta _c^2}{2\sigma _M^2}\right)\frac{dM}{M}$$
(2)
with the variance $`\sigma _M`$ defined as follows:

$$\sigma _M=P(k)W^2(kR)k^2𝑑k,$$
(3)
where $`W(kR)`$ is the window function (top-hat in our case) for filtering fluctuations in the power spectrum on scales characterised by $`R`$ and hence mass $`M=4\pi R^3/3`$. Using our model power spectra from Table 1 (plotted in Fig. 1) together with Eq. (3) and Eq. (2) we are able to compare our numerically achieved mass functions with the PS prediction. The results are shown in Fig. 5 as thin solid lines.

We observe a similar phenomenon as already seen in other comparisons of PS-predicted and $`N`$-body mass functions (Efstathiou et al. 1988, White, Efstathiou & Frenk 1993, Gross et al. 1998, Governato et al. 1999, Jenkins et al. 2001): the PS theory tends to show too many low mass objects (about a factor of 1.2) and too few high mass objects (again a factor of about 1.2). Apart from that the simulations agree fairly well with the Press-Schechter prediction, even in the cases where we modified the power spectra.

Anyway, to allow for better comparison between the individual models and the effect of the features on $`n(>M)`$ we also plot all mass functions derived from the numerical simulations (as already presented in Fig. 5) in one single Figure 6.

#### 4.3.2 Evolution of the Cluster Abundance

The evolution of the cluster abundance has proven to be (potentially) one of the key constraints on the density parameter $`\mathrm{\Omega }_0`$ (e.g. Eke et al. 1996, Bahcall et al. 1997, Eke et al. 1998, Bahcall et al. 1999). This leads immediately to the question of how our modifications to the dark matter power spectrum influence this important issue. We have already seen that while we might not find hints in the observational power spectrum for the features under investigation, there are differences in the mass functions $`n(>M)`$. In Fig. 7 we show the evolution of galaxy clusters with mass $`M`$ greater than $`10^{14}`$$`h^1\mathrm{M}_{}`$.

The first thing that catches the eye is the difference between the COBE and cluster-normalised SCDM model when compared to the OCDM model. Only the cluster-normalised SCDM model shows a very steep evolution of the cluster abundance whereas the evolution in the COBE-normalised one agrees more or less with the OCDM model (cf. Bahcall, Fan & Cen 1997, Fig. 2 in their paper). It is now interesting to check whether our features affect only the normalisation of these evolutionary curves or also the slope. As we can see, the variations in the amplitude of the dip (or bump) mainly shift the cluster evolution up and down, and the location has an influence primarily on the slope of the curve.

This can be understood in the following way. Features on the 8$`h^1\mathrm{Mpc}`$ scale are ’dynamically’ important, in that this scale has evolved from the linear regime ($`z1`$) into quasi-linear or non-linear regime by now. Since the relative abundance of objects increases dramatically in the non-linear regime, we consequently expect still more (less) objects at late times if power had been added (subtracted) on the corresponding scale, i.e. features on this scale are expected to affect the shape of the abundance evolution of corresponding mass objects. Conversely, features on scales that up to present have remained in the linear regime only add (subtract) power on dynamically unimportant scales, thus only affecting the normalisation. 111Though, if the difference in power is very large (c.f. abundance evolution for $`\sigma _8=1.18`$ and $`\sigma _8=0.52`$) there is a significant change in the slope of the abundance evolution, too.

Moving the dip from the non-linear scale of 8$`h^1\mathrm{Mpc}`$ to the semi-linear 16$`h^1\mathrm{Mpc}`$ mainly results in an overall lack of power and therefore we expect the cluster evolution to be closer to the cluster-normalised ($`\sigma _8=0.52`$) Model 1b. However, the evolution for Model 4 agrees fairly well with the OCDM Model. This can be seen even better in Fig. 8 where we normalised the cluster abundance at redshift $`z=0`$ to unity for all models to allow for better comparison of the slopes of the curves.

A similar effect was already observed by Barriga et al. (2000) where the influence of a step-like feature in the primordial power spectrum was studied in the context of phase transitions during inflation. However, their features are on larger scales and thus would leave traces in current epoch measurements of the CMB and large scale structure.

### 4.4 $`\sigma _8`$ Evolution

As the number of objects formed depends on the normalisation of the input power spectrum measured via $`\sigma _8`$ (cf. Eq. 3), it is interesting to check the evolution of this quantity with redshift, too. Fig. 9 shows this evolution for all models. Again, all curves are normalised to the value at redshift $`z=0`$ to allow for better comparison of the slopes. We observe no overlap of any feature model with the OCDM model: regarding the evolution of the dark matter fluctuations measured in 8$`h^1\mathrm{Mpc}`$ spheres (rather than the cluster evolution) the OCDM model still differs significantly from all other models. Moreover, it is difficult again to distinguish between our feature models as well as the differently normalised SCDM models itself. Only Model 6 shows a steeper evolution from redshift $`z1`$ to $`z0.5`$. Hence $`\sigma _8`$-evolution does not provide a discriminant of features in the power spectrum. However, Robinson, Gawiser and Silk (2000) have shown that $`\sigma _8`$ and cluster abundance combine to probe non-gaussianity.

### 4.5 Halo-Halo Correlation

In Fig. 10 we present the halo-halo correlation function for our models. We fixed the number density to $`n=510^4h^3\mathrm{Mpc}^3`$ which actually means to only use the $`N=nV`$ most massive halos with $`V=256`$$`h^1\mathrm{Mpc}`$ being our simulation volume.

We now observe more obvious differences between these models, which can be fully ascribed to our artificial modifications of the initial power spectrum. The amplitude of the cluster correlation function is actually insensitive to the amplitude of fluctuations in the density field (Croft & Efstathiou 1994), which is clearly reflected when comparing Model 1a and 1b. Moreover, there are also only moderate changes in the cluster correlation function when varying the cosmological parameters from our SCDM Model 1a/b to the OCDM model (cf. Martel & Matzner 2000). But as we move on to our feature models not only the amplitude of the correlation function increases (for the dip models), but also the slope changes (i.e. Model 6). This suggests a strongly (and non-linearly) biased formation of galaxy clusters in those models compared to a ’normal’ SCDM model.

As there was no obvious biasing observed in the dark matter power spectrum (which is nothing more than the Fourier transform of the dark matter correlation function), the effect seen in Fig. 10 is entirely due to a different cluster formation scenario in our dip models. This is caused by the fact that by adding (subtracting) power only on a certain scale the relative strength of different waves becomes more and more important. If we put less power into high-$`k`$ waves (i.e. Model 3 and 4) but leave the largest waves unmodified we eventually bias the formation process of galaxy clusters with respect to an unmodified SCDM model as can be observed in Fig. 10. Massive clusters indeed are highly biased, and the extent to which one may have difficulty in accounting for this in standard low density models may constitute the strongest signature of a possible feature at redshift $`z=0`$ in a SCDM model.

### 4.6 Density Profiles

We have already seen that we introduced a strong bias in the cluster formation process by altering the ratio of high- to low-$`k`$ waves amplitudes. This immediately raises the question of whether this will subsequently lead to deviations in the shapes of the clusters themselves. We have therefore calculated the density profiles for massive FOF groups ($`M>10^{15}`$$`h^1\mathrm{M}_{}`$) and show the results for a representative high-mass halo $`M3.510^{15}`$$`h^1\mathrm{M}_{}`$ ($``$ 1700 particles) in Fig. 11. The curves for models 2 through 6 have been shifted upwards by successive factors of 2, whereas the OCDM profile was lowered by a factor of 2.

We can see that modifying the power spectrum does not affect the slope and amplitude of the density profiles. Only the SCDM model 1b deviates from the other curves as the corresponding halo consists of only about 500 particles. Due to the resolution limits of our simulations, we are not able to quantify the substructure content of these galaxy clusters in more detail as the most massive objects contain ’only’ 2000-2500 particles. But a visual comparison of the particle distributions in those clusters provided no obvious differences.

## 5 Analytical Estimates

We have seen above that SCDM models with features at small scales may display an evolution of their cluster abundance that is similar to that in the OCDM model to within some constant offset in amplitude. The question remains of how far we can push this agreement and to what extent this also holds for cluster masses larger than the ones we were able to investigate numerically. To this end we performed Press-Schechter (Press & Schechter 1974) calculations to investigate a larger set of models. However, since we are considering the evolution of cluster abundances, this also requires a relation for the redshift dependence of the critical overdensity $`\delta _c`$, which is left unspecified in the original PS approach. We determined approximate relations for the unmodified SCDM and OCDM models by matching their cluster abundance evolution to that of the corresponding N-body simulations:

$$\delta _c^{SCDM}(a)=1.2+0.6a,\delta _c^{OCDM}(a)=1.45+0.3a$$
(4)
where $`a=1/(1+z)`$ is the cosmic scale factor at redshift z.

For the OCDM model the $`\delta _c`$ relation almost exactly reproduces the abundance evolution. The corresponding relation for the SCDM model was rather chosen to fit both unmodified SCDM and featured SCDM models reasonably well, matching the N-body evolution to within 10 percent for the unmodified SCDM model, typically to within 20 percent for models 2 to 5 and 30 percent for model 6. This is shown in figures 12 and 13. An accurate fit to the unmodified SCDM model only would have required $`\delta _c(a)=1.25+0.5a`$.

However, the agreement between N-body and PS data should be viewed with care, particularly when large sharp features are included: At the location of the features, we essentially introduce a more rapidly changing spectral index. This might contribute to a larger discrepancy between N-body and PS data. To our knowledge agreement between PS and numerical simulations has so far only been established for power spectra with constant (power law) or slowly varying (e.g. SCDM) spectral indices - it is not obvious why this should also hold when a sharply varying index is being introduced, particularly on scales that only now turn non-linear. In the following, however, we will assume that this is not a problem and use the $`\delta _c^{SCDM}`$ relation from equation (4) as a best estimate for PS predictions of the abundance evolution for other small scale features in the SCDM model. We also assume that the $`\delta _c`$ relation for OCDM describes the abundance evolution in higher mass cuts as well as it does for $`M>10^{14}`$ $`h^1\mathrm{M}_{}`$in our simulations.

As seen above, one might expect that adding (subtracting) power only at specific scales leads to a larger (smaller) abundance of objects with mass corresponding to these scales. This is not obvious for scales that are already in the non-linear regime, as we have a coupling across scales, and also objects not only get newly created but now have also been partly incorporated into larger objects, i.e. the effect of a feature on a specific scale spreads out, not only in the evolution of the power spectrum as we have seen above, but also in the abundance of objects of corresponding mass. However, for broad features we may assume that the latter is less important, since a whole range of scales and corresponding masses are affected.

To check this, we placed a broad Gaussian dip in the SCDM model ($`A=0.55`$, $`\sigma _{\mathrm{mod}}=2.5`$) at 10$`h^1\mathrm{Mpc}`$, which is expected to affect a similarly broad range of masses centred on $`M10^{14}10^{15}M_{}h^1`$. A broad feature of this kind is naturally generated by e.g. including a running spectral index term when expanding the primordial perturbation spectrum (Hannestad, Hansen & Villante 2000).

Comparing Fig. 12 with Fig. 14 confirms that the abundance of masses in these ranges ($`M>10^{14}`$, $`M>810^{14}`$$`h^1\mathrm{M}_{}`$, the mass range covered by e.g. Bahcall & Fan 1998) is significantly suppressed. What is more, the abundance evolution for the mass cuts agrees in logarithmic slope and magnitude with that in the OCDM model to within a factor of less than two and three across all redshifts up to z = 1 for the lower and upper mass cuts respectively.

Given that observations of the most distant massive clusters currently only constrain the abundance to within one or two orders of magnitude (see e.g. Bahcall & Fan 1998), SCDM with broad features towards small scales, mimicking (in shape and absolute magnitude) the abundance evolution expected in OCDM models, can therefore not be ruled out. A similar result is obtained by Barriga et al.(2001) who consider primordial step-like features in the context of phase transitions during inflation, however their features are on larger scales in the linear regime and thus would also leave traces in current epoch measurements of the CMB and large scale structure.

## 6 Conclusions

In this paper we presented a series of simulations all based on the COBE normalised SCDM model but with Gaussian features added to an otherwise unperturbed power spectrum (called ’feature’ models). Such bumps (or dips) might naturally arise from non-standard inflationary theories and our main purpose was to investigate their influence on the large-scale structure of the Universe as measured via galaxy clusters. We analysed the evolution of the dark matter power spectrum, the large-scale velocity field represented by the velocity distribution of galaxy clusters, the evolution of the cluster abundance, the halo-halo correlation function, and density profiles of clusters. We furthermore compared our numerical results to analytical predictions based on the Press-Schechter formalism (Press & Schechter 1974).

When comparing the modified SCDM models to a fiducial OCDM model, we can see that when choosing an appropriate scale for the added feature the histories of cluster abundance evolution might be indistinguishable. However, there is still a discrepancy in the overall normalisation left at the present epoch, with only the slopes of the cluster abundance evolution coinciding. But as observations of the most distant massive clusters only constrain the current abundance to within one or two orders of magnitude (cf. Bahcall & Fan 1998), an SCDM model including such a broad feature cannot be ruled out.

Finally we remark that whatever the origin or nature is of these bumps and dips, their effect on the evolutionary history of the SCDM model is much more moderate than changes in the cosmological parameters of the model itself, e.g. lowering the normalisation from $`\sigma _8=1.18`$ (COBE normalisation) to $`\sigma _8=0.52`$ (cluster normalisation). It is therefore an observational challenge to find traces of such features, and the best place to search for them might be the cluster-cluster correlation function which appears to be highly biased with respect to the ’normal’ SCDM model. This is due to an unusual ratio of power for high- and low-$`k`$ waves, which will be reflected in the aforementioned ’more than biased’ cluster formation scenarios.

## Acknowledgements

We are grateful to the referee Mirt Gramann for helpful comments and suggestions. We also acknowledge the use of Hugh Couchman’s AP3M code. RRI gratefully acknowledges a Graduate Studentship from Oxford University and support from St Cross College, Oxford.


format: text

text:
# Interpreting the Cosmic Infrared Background: Constraints on the Evolution of the Dust Enshrouded Star Formation Rate

## 1 Introduction

The Extragalactic Background Light (EBL) in the infrared, also referred to as the Cosmic Infra-Red Background (CIRB), is a record of the emission, absorption and re-radiation of photons integrated over cosmic history. It provides a valuable constraint on theories of galaxy formation and evolution. The EBL at near-infrared wavelengths is due to redshifted radiation from stars. At mid-infrared (MIR) wavelengths, the background is due to redshifted emission from dust that consists of the polycyclic aromatic hydrocarbon (PAH) features and very small grains transiently heated to $`T200K`$ in individual galaxies. At far-infrared (FIR) wavelengths, the dominant contributor is thought to be cold dust ($`T20K`$) that is heated by the ambient interstellar radiation field in galaxies. The recent detection of this background at 2.2 $`\mu `$m, 3.5 $`\mu `$m, 140 $`\mu `$m and 240 $`\mu `$m using DIRBE/COBE data, and in the 125$``$2000 $`\mu m`$ range using FIRAS/COBE measurements by various groups (Puget et al., 1996; Dwek & Arendt, 1998; Fixsen et al., 1998; Hauser et al., 1998; Schlegel et al., 1998; Lagache et al., 1999; Gorjian, Wright & Chary, 2000; Wright & Reese, 2000; Wright, 2001), has indicated that the intensity of the optical/near-infrared background is roughly equal to that of the far-infrared background. This implies that about 50% of the integrated rest-frame optical/ultraviolet (UV) emission from stars and other objects is thermally reprocessed by dust and radiated at mid- and far-infrared wavelengths. Thus, star formation rates that are derived from rest-frame optical/UV luminosities of galaxies are a lower limit to the true star formation rate (e.g. Madau, Pozzetti, & Dickinson, 1998; Meurer et al., 1999; Steidel et al., 1999; Yan et al., 1999).

The first good evidence of this came from the IRAS sky survey which revealed a new population of galaxies with $`L_{IR}=L(81000\mu m)10^{11}L_{\mathrm{}}`$ (see review by Sanders & Mirabel, 1996). Those with $`L_{IR}10^{12}L_{\mathrm{}}`$ were classified as ultraluminous infrared galaxies (ULIGs) while galaxies with $`10^{12}>L_{IR}10^{11}L_{\mathrm{}}`$ were classified as luminous infrared galaxies (LIGs) 111Previously, the term LIG was used for all objects with $`L_{IR}10^{11}L_{\mathrm{}}`$. We use “infrared luminous galaxies” when referring to both LIGs and ULIGs collectively.. These objects exhibited the largest known star formation rates of all local galaxies but had $``$90% of the bolometric luminosity being emitted in the far-infrared (40$``$500 $`\mu `$m) indicating that dust reprocessing is a significant parameter that needs to be considered in estimates of star formation in certain galaxies (e.g. Soifer et al., 1986). However, in the local universe, the integrated bolometric luminosity density of “normal” optically selected galaxies is $`L_{Bol}=4\times 10^8L_{\mathrm{}}Mpc^3`$ while that of infrared luminous galaxies is $`8\times 10^6L_{\mathrm{}}Mpc^3`$ i.e. 50 times less (Soifer et al., 1987). This seemed to indicate that the contribution from LIGs and ULIGs is sufficiently small that they need to be considered only as extreme cases.

Spectroscopic follow-up of the faint IRAS population which covered a relatively small redshift range ($`z<0.27`$) indicated that infrared luminous galaxies were more numerous in the past than they are today and may have a significantly larger contribution to the integrated luminosity density than inferred from observations of the local universe (Kim & Sanders, 1998). Deeper observations that trace the far-infrared luminosity of galaxies to high redshift are difficult since cirrus and confusion noise rapidly begins to dominate.

The ISOCAM guaranteed time extragalactic surveys in conjunction with the ELAIS survey and observations of the lensing cluster Abell 2390 covered a range of flux densities between 50 $`\mu `$Jy and 50 mJy at 15 $`\mathrm{\mu m}`$ (Altieri et al., 1999; Elbaz et al., 1999; Serjeant et al., 2000). The differential counts resulting from these surveys revealed that the counts of galaxies increase quite rapidly as $`S_\nu ^3`$ at brighter flux levels ($`S_\nu >0.4`$ mJy) and then flatten out as $`S_\nu ^{1.6}`$ at fainter levels. The observed mid-infrared counts are an order of magnitude higher than expected if the local mid-infrared luminosity function were not evolving with redshift. This rapid increase in mid-infrared luminous galaxies has been modeled as a $`(1+z)^{4.5}`$ luminosity evolution in the 15 $`\mu `$m local luminosity function by Xu (2000), and as a combination of number density and luminosity evolution by Franceschini et al. (2001). This evolution is much stronger than observed in the UV by Cowie, Songaila & Barger (1999) who find that the comoving UV luminosity density evolves as $`(1+z)^{1.5}`$, instead of the $`(1+z)^{3.9\pm 0.75}`$ initially proposed by Lilly et al. (1996). Furthermore, observations of galaxies in the local universe have shown that the mid-infrared and infrared luminosities are well correlated (Section 2). The mid-infrared luminosity of $`70`$% of the sources seen in the ISOCAM surveys, particularly in the Hubble Deep Field-North and flanking fields (HDFN+FF) translates to an infrared luminosity greater than $`10^{11}L_{\mathrm{}}`$ implying that the majority of them are LIGs and ULIGs (Elbaz et al., 2001). At $`z0.8`$, the mid-infrared luminosity density derived from the ISOCAM 15 $`\mu `$m sources is $``$7$`\times `$10$`{}_{}{}^{7}L_{\mathrm{}}^{}Mpc^3`$ while the $`81000\mu m`$ luminosity density adopting the mid to far-infrared correlation seen in the local universe is $``$5$`\times `$10$`{}_{}{}^{8}L_{\mathrm{}}^{}Mpc^3`$. In comparison, at $`z`$0, LIGs and ULIGs contribute $`10^6`$ $`L_{\mathrm{}}Mpc^3`$ to the 12 $`\mu `$m and 15 $`\mu `$m luminosity density and 7.8$`\times `$106 $`L_{\mathrm{}}Mpc^3`$ to the infrared luminosity density, as derived from the local luminosity function (LLF) of Soifer et al. (1987); Fang et al. (1998); Xu et al. (1998). This indicates an increase by a factor of about 60 between $`z`$0 and $`z`$0.8, providing further evidence for an evolution in the infrared luminosity function (IRLF) with redshift.

Similar deep surveys have been conducted at 850 $`\mu `$m using the SCUBA instrument on the JCMT (Hughes et al., 1998; Barger, Cowie, & Sanders, 1999; Blain et al., 1999a; Eales et al., 2000). The large beam size (14$`\mathrm{}`$ FWHM) and the negative $`k`$-correction in this wavelength regime makes identification of the optical counterparts and thereby the redshift distribution of the sources very difficult. High resolution radio interferometric observations and the use of 450 $`\mu `$m/850 $`\mu `$m flux ratios have helped somewhat in localizing the sources and constraining the redshifts (Hughes et al., 1998; Barger, Cowie, & Richards, 2000). These have placed the bright ($`S_\nu >6`$ mJy) submillimeter sources at $`z13`$, which has been confirmed by the more extensive survey of Chapman et al. (2001). The implication of this is that most, if not all the submillimeter sources are extreme ULIGs with star formation rates of $`10^210^3M_{\mathrm{}}yr^1`$. Furthermore, the star formation rate density due to ULIGs must have increased by about 2 orders of magnitude between $`z`$0 and $`z13`$.

Many of the LIGs and ULIGs in the local universe show morphological signatures of interaction, and $`>`$50% of the optical counterparts of ISOCAM HDFN galaxies show evidence for interactions (Mann et al., 1997). Surveys at visible wavelengths show a redshift evolution of the merger fraction, defined as the fraction of close pairs of galaxies, as $`(1+z)^3`$ (e.g. Le Fevre et al., 2000). Thus, if mergers were indeed a tracer of LIGs and ULIGs, this would again suggest that the bright end of the IRLF is evolving strongly. However, the faint end of the IRLF is very poorly constrained at $`z1`$ since none of the long wavelength surveys are sensitive enough to detect galaxies with $`L_{IR}<10^{11}L_{\mathrm{}}`$ at $`z>0.5`$. Meurer et al. (1999) has shown that the FIR to UV flux ratio is closely related to the UV slope for normal starbursts but that the relationship breaks down for ULIGs (Meurer et al., 2000). This indicates that the visible/near-infrared counts can potentially place constraints on the evolution of the faint end of the infrared luminosity function but we postpone this discussion to the future.

In this paper, we combine data from a variety of published surveys of nearby galaxies to determine the correlation, if any, between the luminosities at various mid- and far-infrared wavelengths. We use these correlations to generate smoothly varying spectral energy distributions for galaxies as a function of luminosity class. We assess the need for luminosity and density evolution in the 15 $`\mu `$m luminosity function of Xu et al. (1998) and therefore the 60 $`\mu `$m luminosity function of Soifer et al. (1987) based on fits to the ISOCAM 15 $`\mu `$m, ISOPHOT 90 $`\mu `$m, ISOPHOT 170 $`\mu `$m and SCUBA 850 $`\mu `$m galaxy counts, as well as the spectrum of the CIRB at $`\lambda >5`$ $`\mu `$m. The evolution of the mid-infrared luminosity function is then translated to an estimate of the dust enshrouded star formation rate density (SFR) as a function of redshift and compared with SFR values derived from optical/near-infrared surveys. We adopt a $`H_0`$=75 km s-1 Mpc-1, $`\mathrm{\Omega }_M=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$ cosmology throughout this paper unless otherwise explicitly stated.

## 2 Luminosity Correlations in the Infrared and Template Spectral Energy Distributions

It can be shown that the 12 $`\mu `$m and far-infrared luminosities of galaxies in the IRAS Bright Galaxy Sample (BGS) cannot be accurately derived from their $`B`$band luminosities (Soifer et al., 1987)222The Zwicky magnitudes, $`m_z`$, in the BGS were converted to B-luminosities using $`m_B=m_z0.14`$ and a $`B`$band zero point of 4260 Jy.. The peak to peak scatter in the $`L_B/L_{IR}`$ ratio for a fixed $`L_{IR}`$ is about a factor of 20. However, as mentioned earlier, the FIR to UV flux ratio has been shown to be closely related to the UV slope for normal starbursts (Meurer et al., 1999). This relationship breaks down for the ULIGs (Meurer et al., 2000). The phenomenon can be qualitatively explained by the fact that the UV emission arises from stars which are relatively unobscured to the observer. Regions of star formation with large optical depth i.e. H II regions, could exist in the central regions of galaxies where almost all the UV light is reprocessed to the MIR and FIR. Thus, the regions of FIR and UV emission would be unrelated, especially for ULIGs. The best observational evidence for this explanation can be seen in the Antennae galaxy (Mirabel et al., 1998) where about half the 15 $`\mu `$m emission seen by ISOCAM arises from regions that are inconspicuous at visible wavelengths. The breakdown in the FIR to UV slope correlation for ULIGs is problematic for determination of the true star formation rate from optical/UV surveys, since submillimeter observations using the SCUBA instrument indicate that ULIGs might have a larger contribution to the star formation rate density at high redshift. This suggests that it will be difficult to determine the true SFR by applying an accurate extinction correction to the optical/UV determined value.

Since the short wavelength starlight and dust emission are not closely related, an estimate of the dust enshrouded star formation rate can only be derived from other tracers such as the mid-infrared and far-infrared luminosities or using the radio to far-infrared correlation shown by Condon (1992). In the mid-infrared regime, the spectra of galaxies exhibit broad emission features at 6.2, 7.7, 8.6, 11.3 and 12.7 $`\mu `$m which are probably from polycyclic aromatic hydrocarbons (see review by Puget & Leger, 1989). These features and their associated continuum dominate the emission at mid-infrared wavelengths shortward of 10 $`\mu `$m. There is, in addition, a continuum from very small grains (VSGs) of size $`<`$10 nm which dominates the emission above $``$10 $`\mu `$m (Désert et al., 1990; Laurent et al., 2000) except for quiescent star forming galaxies. The VSGs get transiently heated to temperatures of $``$200 K by the ambient optical/UV continuum which is proportional to the star formation activity. In addition, mid-infrared measurements do not need large extinction corrections since the extinction at mid-infrared wavelengths is only about 1% that at visible wavelengths (Mathis, 1990). The radio wavelengths on the other hand are dominated by free-free emission from HII regions and synchrotron emission from supernova remnants. Although radio observations are almost confusion limited at an 8.5 GHz sensitivity of 9 $`\mu `$Jy obtained over the HDF (Richards et al., 1998), we find that they are typically as sensitive as the ISOCAM 15 $`\mu `$m observations in that they can typically probe galaxies with $`L_{IR}10^{11.4}L_{\mathrm{}}`$ at $`z1`$. Since the ISOCAM 15 $`\mu `$m observations provide the primary constraint on evolution models at $`z<1.2`$, we adopt as a starting point the local 15 $`\mu `$m luminosity function described in Xu et al. (1998) and Xu (2000). For the rest of this paper, we will use the convention defined in Sanders & Mirabel (1996):

$$L_{IR}=L(81000\mu m)=1.8\times 10^{14}\times 10^{26}[13.48L_{12}+5.16L_{25}+2.58L_{60}+L_{100}]$$
(1)
$$L_{FIR}=L(40500\mu m)=1.6\times 1.26\times 10^{14}\times 10^{26}[2.58L_{60}+L_{100}]$$
(2)
In the above, the symbol $`L_\lambda `$ is defined as $`L_\nu (\lambda \mu m)`$ in units of $`L_{\mathrm{}}Hz^1`$. $`L_{IR}`$ and $`L_{FIR}`$ are in $`L_{\mathrm{}}`$.

To use the mid-infrared luminosity function as a tracer of the dust enshrouded star formation rate, we first need to define a calibration scale. Figures 1 and 2 illustrate the accuracy with which the infrared luminosity of galaxies can be derived from their mid-infrared luminosities. Figure 1 is based on $``$300 galaxies from the IRAS BGS while Figure 2 is based on published ISOCAM and ISOPHOT observations of IRAS galaxies. The data points with $`L_{IR}>10^{10}L_{\mathrm{}}`$ are fit by a first order polynomial shown in Equations 4-6. This is shown as a solid black line in the upper panels of the figure. Objects with $`L_{IR}<10^{10}L_{\mathrm{}}`$ are not used for the fits. At these low luminosities, the fraction of the bolometric luminosity emitted in the infrared is $`<`$50% i.e. visible starlight that is not obscured by dust is dominating the radiated energy of the galaxy. On the other hand a significant fraction of the more luminous objects show disturbed morphologies, suggesting interactions with other galaxies that would result in gas-rich systems with star formation in highly obscured regions. The lower plots in both figures show the scatter in the ratio of the infrared luminosity as derived from IRAS data for these galaxies to the infrared luminosity derived from the polynomial fits. Also shown is the 1$`\sigma `$ uncertainty in the derived infrared luminosity calculated as the range within which 68% of the galaxies lie. The lowest luminosity objects ($`L_{IR}<10^9L_{\mathrm{}}`$) have been rejected in the lower plots.

The 6.7 $`\mu `$m luminosities were derived from ISO observations of $``$90 nearby starburst dominated galaxies. 44 spiral and starburst galaxies had photometry from ISOCAM (Chanial et al., 2001; Roussel et al., 2001; Laurent et al., 2000), 8 ULIGs had spectra from ISOCAM circular variable filters (CVF) observations (Tran et al., 2001) while 37 ULIGs had ISOPHOT mid-infrared spectra (Rigopoulou et al., 1999). Rigopoulou et al. (1999) obtained mid-infrared spectra with ISOPHOT of about 60 ULIGs and about 15 low luminosity starbursts and normal galaxies to study the emission features from the PAHs. Of the 60 ULIGs, about 45 had the 7.7 $`\mu `$m PAH feature detected with good signal to noise. However, the calibration and performance of the instrument is not very well determined. To assess the quality of the data set, we compared the ISOPHOT observations on 5 ULIGs and 4 low luminosity starbursts/normal galaxies to the Chanial et al. (2001) ISOCAM LW2 observations of the same galaxies333The LW2 filter is broad enough to include the 7.7 $`\mu `$m PAH feature but is centered at 6.75 $`\mu `$m.. We find that for the ULIGs, the ratio of 7.7 $`\mu `$m line+continuum flux density as published in Rigopoulou et al. to the ISOCAM 6.7 $`\mu `$m flux density lies in the range 1.5$``$3.0. In comparison, for the starbursts/normal galaxies, the ratio of 7.7 $`\mu `$m line+continuum to the ISOCAM 6.7 $`\mu `$m flux density falls in the range 0.2$``$1.6, a factor of 8. So, for assessing the correlation between the mid and far-infrared luminosities, we consider the ISOCAM data on bright IRAS galaxies as well as the ULIG sample of Rigopoulou et al. (1999), dividing the line+continuum flux value published in the latter by 2.4 and assigning a peak to peak error bar of a factor of 2. This is consistent with the range of 1.5$``$2.7 that we find for the 7.7 $`\mu `$m line+continuum to the 6.7 $`\mu `$m flux density ratio in the ISOCAM CVF observations of Tran et al. (2001).

The 15 $`\mu `$m and infrared luminosities of 120 IRAS galaxies were taken from the sample of Chanial et al. (2001) and the survey performed in the North Ecliptic Pole Region (NEPR) by Aussel et al. (2000). The NEPR sample of galaxies only has IRAS 60 $`\mu `$m luminosities available and we have converted these to a far-infrared luminosity based on a 60 $`\mu `$m to far-infrared correlation derived by combining the IRAS BGS and the IRAS PSCz catalog of Saunders et al. (2000). The far-infrared luminosity is typically about 83% of the total infrared luminosity and we have applied this conversion to be consistent with the other plots.

Clearly, the 6.7 $`\mu `$m, 12 $`\mu `$m and 15 $`\mu `$m luminosities of galaxies trace the infrared luminosity much better than the $`B`$band luminosity. The 15 $`\mu `$m-to-IR and 12 $`\mu `$m-to-IR correlations were determined from a first order polynomial fit to the data with $`L_{IR}>10^{10}L_{\mathrm{}}`$. Since the 6.7 $`\mu `$m-to-IR correlation is based on different data sets from CVF and broadband photometry with ISOCAM and spectroscopy using ISOPHOT, the polynomial fit for the 6.7 $`\mu `$m-to-IR correlation was determined by applying a $`k`$correction from 6.7 $`\mu `$m to 15 $`\mu `$m based on ISOCAM observations of nearby galaxies (Figure 3) and then using the 15 $`\mu `$m-to-IR correlation. The mid-infrared to IR correlations show a similar scatter around the correlation line which is about a factor of 5 better than the optical-to-IR correlation.

These data sets tentatively illustrate the potential of using the mid-infrared as a tracer of dust enshrouded star formation and a more homogeneous and comprehensive survey of nearby galaxies, as will be undertaken by SIRTF, will be required to either strengthen or reject this correlation.

Kennicutt (1998) has transformed the infrared luminosity of young $`<10^8`$ yr starburst galaxies to a star formation rate. If we adopt the correlations shown in the previous figures, we can translate the mid-infrared luminosity of galaxies with $`L_{IR}>10^{10}L_{\mathrm{}}`$ to an approximate estimate of the dust enshrouded SFR ($`\rho ^{}`$) using the formulae:

$`\rho ^{}(M_{\mathrm{}}yr^1)`$ $`=`$ $`1.71\times 10^{10}L_{IR}(L_{\mathrm{}})`$ (3)
$`L_{IR}`$ $`=`$ $`11.1_{3.7}^{+5.5}\times L_{15\mu m}^{0.998}`$ (4)
$`L_{IR}`$ $`=`$ $`0.89_{0.27}^{+0.38}\times L_{12\mu m}^{1.094}`$ (5)
$`L_{IR}`$ $`=`$ $`4.37_{2.13}^{+2.35}\times 10^6\times L_{6.7\mu m}^{1.62}`$ (6)
where all values are in $`L_{\mathrm{}}`$. The 1$`\sigma `$ values have been estimated by calculating the range of values within which 68% of galaxies have their observed infrared luminosities.

As mentioned earlier, the main observational constraints on models that trace the redshift evolution of the IRLF are:

1. differential counts from various surveys at mid-infrared, far-infrared and submillimeter wavelengths and
2. the spectrum of the CIRB at $`\lambda >5`$ $`\mu `$m.

To use these constraints, it is necessary to know the luminosity at different wavelengths for galaxies in each luminosity bin of the infrared luminosity function. This motivates the generation of template spectra for objects of different luminosity classes. It is useful to note that many evolutionary models which have already been developed either use a mid-infrared template which is a very poor representation of the true PAH emission features (e.g. Rowan-Robinson, 2000; Pearson, 2001; Sadat, Guiderdoni & Silk, 2001) or neglect the PAH features altogether (Malkan & Stecker, 1998, 2001). It is trivial to show that this makes a critical difference in the quality of the fits to the ISOCAM mid-infrared number counts and therefore affects the evolution parameters particularly at $`z<1`$. So, it is important to have template SEDs that reproduce the observed trend in the luminosity of local galaxies at different wavelengths.

Using MIR, FIR and submillimeter data from ISOCAM, IRAS and SCUBA observations of nearby galaxies, we fit the observed trend between different mid- and far-infrared luminosities as shown in Figure 3. The top two panels in the figure show the ISOCAM observations at 6.7 $`\mu `$m and 15 $`\mu `$m for $``$50 IRAS galaxies which are described in Chanial et al. (2001). The solid lines in the panels show two first order polynomial fits, one for galaxies with 15 $`\mu `$m luminosity less than $`2\times 10^9L_{\mathrm{}}`$ and another for more luminous galaxies. This is because the ratio between the mid-infrared luminosities changes as a function of 15 $`\mu `$m luminosity, possibly due to enhanced emission from the VSG component (Laurent et al., 2000). The luminosity break corresponds to $`L_{IR}2\times 10^{10}L_{\mathrm{}}`$ which is similar to the luminosity cutoff used for deriving Equations 4-6. The panel showing the 15 $`\mu `$m to 60 $`\mu `$m trend consists of data described in Figure 2. The 60 $`\mu `$m to FIR correlation is for the IRAS BGS and PSCz galaxies while the panel showing the IR to FIR correlation is only for the IRAS BGS galaxies. The last panel shows SCUBA submillimeter data on $``$100 IRAS galaxies (Dunne et al., 2000). In addition, the correlation between IRAS 25 $`\mu `$m and 100 $`\mu `$m luminosities for galaxies in the BGS was also determined. The solid lines for the four lower panels utilize only a single first order polynomial fit to all the data points. Also shown in the panels as red triangles are the luminosities at the corresponding wavelength for the different templates that were generated as described below.

Template spectral energy distributions (SEDs) were generated between 0.1 and 1000 $`\mu `$m to reproduce the observed trend between mid-infrared and far-infrared luminosities. To generate these templates, we used the basic Silva et al. (1998) models to reproduce the ultraviolet-submillimeter SED of 4 prototypical galaxies: Arp220, NGC6090, M82 and M51. These correspond to objects of four different luminosity classes - ULIGs, LIGs, “starbursts” (SB) and “normal galaxies” respectively. ISOCAM CVF observations between 3 and 18 $`\mu `$m of these galaxies provided new data on the relative strength of the mid-infrared features and continuum (Charmandaris et al., 1999; Laurent et al., 2000; Forster-Schreiber et al., 2001; Roussel et al., 2001). The mid-infrared region of the modelled spectra were then replaced with the ISOCAM observations. In addition, corrections were made for the 17.9 $`\mu `$m silicate feature based on observations by Smith et al. (1989). The four template spectra were checked to ensure that the IRAS observed values of these four galaxies were reproduced. We then partitioned the four templates into a mid-infrared (4$``$20 $`\mu `$m) and far-infrared (20$``$1000 $`\mu `$m) component and interpolated between the four to generate a range of mid-infrared and far-infrared sample templates of intermediate luminosity. An additional set of far-infrared templates provided by Dale et al. (2000) were added to the ensemble of far-infrared templates to span a wider range of spectral shapes.

For each luminosity bin of the 15 $`\mu `$m luminosity function, the luminosity at the following wavelengths: 6.7, 12, 25, 60, 100, 850 $`\mu `$m was predicted based on the polynomial fits to the data shown in Figure 3. Of the $`100`$ mid-infrared sample templates generated as described above, the mid-infrared template that best fits the predicted 6.7, 12 and 15 $`\mu `$m luminosities was selected. Similarly, the far-infrared template that best fits the predicted 25, 60, 100 and 850 $`\mu `$m luminosities was selected. The luminosity of the templates at the corresponding wavelengths was determined by integrating over the filter curves of the instruments. Our goal was only to generate SEDs that reproduce the observed trend in luminosities at different wavelengths. Selecting a variety of sample templates provided better fits to the predicted luminosities than by just interpolating between the four SEDs generated by the Silva et al. (1998) models. The best fitting mid-infrared and far-infrared templates were then merged together to provide the final template SED for each luminosity bin. The red triangles in Figure 3 are the luminosities at the corresponding wavelengths from the final merged template SEDs. The $`B`$band luminosity of galaxies in the IRAS BGS shown in Figure 1 was also used to constrain the optical/near-infrared spectral energy distribution of galaxies but as stated before, we have not constrained the UV slope of the template SEDs. The absence of a good correlation between the $`B`$band and IR luminosities implies that the optical/near-infrared part of our SEDs is highly uncertain. This is not a major problem since we are only analyzing the dust emission in this paper. The templates for three objects with infrared luminosities of $`10^{10},10^{11}`$ and $`10^{12}L_{\mathrm{}}`$ along with the predicted luminosities at different wavelengths based on the correlations in Figure 3 are shown in Figure 4.

## 3 Evolution of the 15 $`\mu `$m and Far-Infrared Local Luminosity Function

Since our intention is to use the different mid- and far-infrared observational constraints to estimate the evolution of the dust enshrouded star formation rate with redshift, we use the 15 $`\mu `$m LLF as a tracer of dust emission in the local universe. Xu et al. (1998) and Xu (2000) derived a 15 $`\mu `$m LLF based on a correlation between ISOCAM mid-infrared and IRAS mid- and far-infrared data. In addition, estimates of the 12 $`\mu `$m LLF have been made by Rush, Malkan & Spinoglio (1993) and Fang et al. (1998). The different 12 $`\mu `$m and 15 $`\mu `$m luminosity functions are shown in Figure 5. Also shown is the predicted 60 $`\mu `$m LLF derived from the mid-infrared LLF using a mid-infrared to 60 $`\mu `$m conversion from the polynomial fit to the observations of Aussel et al. (2000) and Chanial et al. (2001) described earlier. All these are in good agreement with each other since they were essentially derived from IRAS observations of nearby galaxies.

Evolution of the luminosity function with respect to redshift can be expressed as:

$`\mathrm{\Psi }(L,z)`$ $`=`$ $`n(z)\varphi [L,z]`$ (7)
$`\varphi (L,z)`$ $`=`$ $`\varphi [\frac{L}{g(z)},0]`$ (8)
where $`\mathrm{\Psi }(L,z)`$ is the number density of galaxies as a function of luminosity $`L`$ and redshift $`z`$. The $`n(z)`$ term represents evolution in the number density of galaxies while the $`\varphi [L,z]`$ term represents luminosity evolution. $`\varphi [L,0]`$ is the local luminosity function. We consider models where $`n(z)`$ is of the form $`n(0)(1+z)^{\alpha _D}`$ up to a turnover redshift $`z_{turn}^D`$ followed by $`n(z_{turn}^D)[(1+z)/(1+z_{turn}^D)]^{\beta _D}`$ up to $`z=4.5`$. The luminosity evolution component $`g(z)=(1+z)^{\alpha _L}`$ up to $`z_{turn}^L`$ followed by $`n(z_{turn}^L)[(1+z)/(1+z_{turn}^L)]^{\beta _L}`$.

It should be emphasized that there is considerable degeneracy in the density and luminosity evolution of galaxies. While density evolution slides the luminosity function along the vertical axis, the latter slides it on the horizontal axis. However, current observations at mid- and far-infrared wavelengths detect galaxies only at the luminous end of the luminosity function, as a result of which the two are indistinguishable. This is shown in Figure 6. The figure also shows that evolving just the luminous end of the LLF i.e. $`L>5\times 10^{10}L_{\mathrm{}}`$ similar to the model proposed by Dole et al. (2000) would result in the same degeneracy with observations although it would result in a somewhat unphysical break in the luminosity function.

There is an additional degeneracy induced by the fraction of the luminosity function that is evolving. Redshift measurements of ISOCAM 15 $`\mu `$m sources in the HDFN+FF indicate that the majority of them are LIGs and ULIGs (Elbaz et al., 2001). Interestingly, when the local 60 $`\mu `$m luminosity function is compared to the Schechter function commonly used to represent the local luminosity function at visible wavelengths, an excess of galaxies is seen in the 60 $`\mu `$m LLF beyond $`L_{60\mu m}>10^{11}L_{\mathrm{}}`$ since the Schechter function drops faster at the bright end. On one hand it seems likely that just this excess of galaxies, most of which show morphological signatures of merger activity, could be evolving at high redshift. Alternatively, it is possible that a luminosity dependent fraction which approaches 100% at $`L_{60\mu m}>10^{11}L_{\mathrm{}}`$ of the LLF could be evolving. Unfortunately, the observational constraints on the faint end of the IRLF are limited since these galaxies are undetected at mid- and far-infrared wavelengths at $`z>0.5`$. The correlation between mid-infrared and visible wavelengths being poor, the counts of galaxies at visible/near-infrared wavelengths cannot be used to constrain the distribution. However, we will investigate in a future paper if the relationship between the FIR/UV flux ratio and UV slope can constrain the evolution of the faint end of the IRLF.

The principal observational constraints on the evolution of the bright end of the luminosity function then are:

* The ISOCAM differential number counts at 15 $`\mu `$m especially the “knee” in the counts slope at 0.4 mJy (Elbaz et al., 1999),
* The 15 $`\mu `$m extragalactic background light which has a lower limit from ISO counts and an upper limit from gamma-ray observations of a TeV flare in Markarian 501,
* The redshift distribution of 15 $`\mu `$m sources in the HDFN+FF which is somewhat peaked at $`z0.8`$ (Aussel et al., 2001) and indicates that 45% of galaxies with $`0.1<S_{15}<0.4`$ mJy are LIGs and 20% are ULIGs with the remaining being normal and low luminosity starburst galaxies (Elbaz et al., 2001),
* The spectrum of the cosmic far-infrared background between 100 and 850 $`\mu `$m as measured by DIRBE and FIRAS on COBE,
* ISOPHOT 90 $`\mu `$m and 170 $`\mu `$m counts and SCUBA 850 $`\mu `$m number counts (Efstathiou et al., 2000; Hughes et al., 1998; Blain et al., 1999a; Barger, Cowie, & Sanders, 1999; Eales et al., 2000; Dole et al., 2001).

The $`k`$correction of galaxies illustrated in Figure 7 clearly illustrates the range of redshifts that can be studied by observations at these wavelengths. The 850 $`\mu `$m observations are typically confusion limited at 2 mJy. However, deeper lensed surveys or using high resolution radio interferometric data can push the detection threshold down to 0.5 mJy, past the confusion limit. This allows detection of objects with $`L_{IR}>5\times 10^{11}L_{\mathrm{}}`$ out to $`z5`$ which transforms to a star-formation rate of $`>`$85 $`M_{\mathrm{}}/yr`$. On the other hand, the ISO mid- and far-infrared observations are typically dominated by galaxies at $`z<1`$ and so can constrain the low redshift turnover in the luminosity function evolution.

Figure 8 illustrates the nature of the counts if the 15 $`\mu `$m luminosity function remained equal to the local one at all redshifts, i.e. no evolution. The first plot shows the ISOCAM 15 $`\mu `$m differential counts from Elbaz et al. (1999), which include the IRAS 12 $`\mu `$m counts converted to 15 $`\mu `$m by Xu (2000) and the ELAIS 15 $`\mu `$m counts of Serjeant et al. (2000) renormalized as in Genzel & Cesarsky (2000). The remaining plots show the ISOPHOT FIRBACK 170 $`\mu `$m counts, ISOPHOT, IRAS BGS and PSCz 90 $`\mu `$m counts and SCUBA 850 $`\mu `$m integral counts. Also shown are the relative contributions from ULIGs, LIGs and $`L_{IR}<10^{11}L_{\mathrm{}}`$ galaxies to the counts at different wavelengths444 From now on we refer to the $`L_{IR}<10^{11}L_{\mathrm{}}`$ galaxies as normal/SB galaxies. Clearly, some form of redshift evolution in the luminosity function is required to fit the counts. It should be emphasized that evolutionary models should be fit to the differential counts at different wavelengths since integral counts tend to smooth over any subtle changes in the galaxy counts slope. This requires that the calibration of the data from different surveys using the same instrument be consistent and accurate.

The ISOPHOT 90 $`\mu `$m data are known to suffer from large calibration uncertainties which are not reflected in the error bars. We show in Section 3.3 that our models fit the 170 $`\mu `$m counts but consistently overpredict the 90 $`\mu `$m counts although both wavelengths probe similar populations of galaxies at $`z<1.2`$. We find that an upward correction of the 90 $`\mu `$m flux densities by 30%, which is well within the calibration uncertainties, leads to excellent agreement between our models and the data. In addition, the faint end PSCz counts are known to suffer from incompleteness (Efstathiou et al., 2000).

Any evolution of the luminosity function must have a turnover at some redshift $`z_{turn}`$ to avoid overproducing the CIRB. In theory, the turnover redshifts $`z_{turn}^D`$ ($`z_{turn}^L`$), the $`\alpha _D(\alpha _L)`$ value for the slope of the density (luminosity) evolution at $`z<z_{turn}`$ and the $`\beta _D(\beta _L)`$ value for $`z>z_{turn}`$ can be different implying that there are effectively six parameters. We however consider models with $`z_{turn}^D`$= $`z_{turn}^L`$ since it is unclear why the turnover for luminosity and density evolution, both of which are probably induced by galaxy interactions, should be different. The range of values for $`\alpha `$, $`\beta `$ and $`z_{turn}`$ selected for our models were $`1\alpha 6`$, $`3\beta 0`$ and $`0.6z_{turn}1.5`$ respectively. Evolutionary models with pure luminosity evolution and pure density evolution are also considered.

### 3.1 Constraints on Pure Density Evolution

We find that pure density evolution of the entire luminosity function cannot reproduce the counts at all the wavelengths. In this scenario, the 15 $`\mu `$m counts are dominated by the normal/SB galaxies, not by LIGs and ULIGs, which is inconsistent with observations in the HDFN+FF. Secondly, the normal/SB galaxies are unable to reproduce the break in the 15 $`\mu `$m counts seen at 0.4 mJy but instead produce a sharp break only at $`S_\nu <0.2`$ mJy.

However, density evolution models which evolve just a fraction of the 15 $`\mu `$m luminosity function with the fraction being $`<`$5% at $`L_{15\mu m}10^9L_{\mathrm{}}`$ and approaching 100% at 15 $`\mu `$m luminosities greater than 8$`\times `$10$`{}_{}{}^{10}L_{\mathrm{}}^{}`$ provide reasonable fits to the data (red broken line in Figure 14). The best fit density evolution parameters then are $`\alpha _D=12.0\pm 0.5`$ up to $`z_{turn}^D=0.7\pm 0.1`$ followed by $`0.5<\beta _D0`$ (Figure 9). The first plot shows the spectrum of the cosmic infrared background with lower limits from integrated counts of galaxies in the optical/UV from Madau & Pozzetti (2000), measurements in the near- and far-infrared using the DIRBE instrument (Hauser et al., 1998; Finkbeiner et al., 2000; Gorjian, Wright & Chary, 2000; Wright, 2001), an estimate of the far-infrared background from FIRAS (Lagache et al., 1999), and lower limits in the mid-infrared, far-infrared and submillimeter from counts of individual galaxies (Elbaz et al., 1999; Blain et al., 1999a; Matsuhara et al., 2000). Also shown is the upper limit on the CIRB from TeV observations of Mrk501 (Stanev & Franceschini, 1998) and the cosmic microwave background at $`\lambda >300`$ $`\mu `$m. The remaining plots show the counts for this evolution model.

Models with $`\alpha _D>13.0`$ result in a significant overprediction of the 170 $`\mu `$m counts while $`\alpha _D<11.0`$ underpredicts the submillimeter counts. Changing the turnover redshift $`z_{turn}^D`$ to high redshift shifts the knee in the 15 $`\mu `$m differential counts to fainter flux levels and vice-versa. The slope of the evolution at $`z_{turn}^D>0.7`$ is mainly constrained by the spectrum of the CIRB at $`\lambda >200`$ $`\mu `$m and the 850 $`\mu `$m counts.

The 170 $`\mu `$m, 90 $`\mu `$m and 15 $`\mu `$m observations all probe the population of galaxies at $`z<1.2`$. However, in this density evolution model, the 15 $`\mu `$m counts are dominated ($``$90%) by LIGs and have $``$5% contribution from ULIGs. In comparison, the 170 $`\mu `$m counts are dominated by ULIGs at redshifts between 0.5$``$1 while the 90 $`\mu `$m counts have roughly equal contributions from all three populations of galaxies. The 850 $`\mu `$m galaxies at flux densities larger than 1 mJy are mainly ULIGs at $`z>1`$ while LIGs between redshifts of $`0.62.0`$ dominate the counts at fainter flux levels.

The pure density evolution model shown appears to overpredict the contribution from LIGs ($``$90%) to the 15 $`\mu `$m number counts and underestimates the bright end 850 $`\mu `$m counts. Secondly, many LIGs and ULIGs have been morphologically associated with disturbed systems. So a steep evolution in the density of objects should reflect in an increase of the merger fraction which is defined as the fraction of galaxies in close pairs. Observationally, the merger fraction when averaged over all galaxies appears to evolve much slower with redshift, approximately as $`(1+z)^3`$ (e.g. Le Fevre et al., 2000). It is possible though that the LIGs and ULIGs have a merger fraction that increases much more rapidly than $`(1+z)^3`$ but this has not been estimated since there are no clear observational signatures of LIGs and ULIGs at visible wavelengths. Although there is no strong observational evidence in favor of pure density evolution of a fraction of the luminosity function, we show in Figure 14 that it does predict the same number density of infrared luminous galaxies at high redshift as other models and so cannot be entirely ruled out.

### 3.2 Constraints on Pure Luminosity Evolution

We have shown above that some form of luminosity evolution is required to avoid deriving large values for the slope of the density evolution. The slope of the luminosity evolution $`\alpha _L`$ is strongly constrained by the mid-infrared number counts. Values of $`\alpha _L<4.5`$ are unable to reproduce the break of the differential counts seen at a 15 $`\mu `$m flux density of 0.4 mJy while $`\alpha _L>5.5`$ results in an overproduction of the counts at bright flux densities (Figure 10). If the entire LLF is evolved (solid black line in Figure 14), $`\alpha _L=5.0`$ up to $`z_{turn}^L`$=0.8 followed by $`0.5<\beta _L0.0`$, the 90 $`\mu `$m differential counts are overproduced but all the other counts are reproduced very well.

This is only partially consistent with the results of Xu (2000) who suggested that the mid-infrared counts can be modelled by evolving the entire LLF by $`L(z)(1+z)^{4.5}`$ for $`z<1.5`$ and by $`L(z)=L(0)\times 2.5^{4.5}`$ for higher redshifts. We find that a luminosity evolution of $`(1+z)^{4.5}`$ up to $`z1`$ followed by $`L(z)=L(0)\times 2.0^{4.5}`$ overpredicts the CIRB and provides only marginal fits to the mid-infrared counts at the faint end. Extending this evolution up to $`z_{turn}^L1.5`$ severely overpredicts the CIRB as well as the observed faint end 15 $`\mu `$m and 850 $`\mu `$m counts.

The main problem with pure luminosity evolution is that the 90 $`\mu `$m counts at $`S_\nu <6`$ Jy are overproduced but as mentioned earlier, this can be resolved by rescaling the ISOPHOT 90 $`\mu `$m flux densities upwards by 30% which is within the calibration uncertainties of the instrument. In addition, the break in the 15 $`\mu `$m differential counts in this evolution model is not as sharp as observed.

As in the pure density evolution model, the counts at 15 $`\mu `$m, 90 $`\mu `$m and 170 $`\mu `$m all trace galaxies at $`z<1.2`$ but the relative contributions from LIGs, ULIGs and normal/SB galaxies differ, with the 90 $`\mu `$m counts having roughly equal contributions from all three populations, the 170 $`\mu `$m counts being dominated by ULIGs and the 15 $`\mu `$m counts being dominated by LIGs and low redshift normal/SB galaxies. Pure luminosity evolution predicts that the contribution from normal/SB galaxies to the 15 $`\mu `$m counts between 0.1 and 0.4 mJy is 35%, similar to that observed in the HDFN+FF. However, the contribution from ULIGs in the same flux density range is found to be only 6% in this model; a factor of $`34`$ smaller than that observed.

### 3.3 Combination of Luminosity and Density Evolution

We have already illustrated the degeneracy and problems with pure luminosity and pure density evolution in the fits to the number counts. As illustrated above, both of them provide reasonable fits to the spectrum of the CIRB and to the differential counts at three wavelengths. An additional degeneracy is introduced when using a combination of luminosity and density evolution with only a fraction of the LLF evolving.

The fraction of the LLF that is evolving is defined as the “dusty starburst” population. We consider these to be galaxies with $`L_B/L_{IR}<0.5`$. This corresponds to $`L_{IR}10^{10.2}L_{\mathrm{}}`$. All galaxies with $`L_{IR}10^{10.2}L_{\mathrm{}}`$ are then evolved while only about 5% of the galaxies with $`L_{IR}<10^{10.2}L_{\mathrm{}}`$ (“normal” galaxy population) is evolved, with a smooth transition between the two (purple broken line in Figure 14). We also considered models where the dusty starburst population is defined at a larger minimum luminosity ($`L_{IR}10^{11}L_{\mathrm{}}`$) but were unable to find evolution parameters that could reasonably reproduce all the data.

Our best fit model using a combination of both density and luminosity evolution for a fraction of the LLF as defined above is shown in Figure 11. Almost all the counts are reproduced quite well, with the exception of the 90 $`\mu `$m counts from ISOPHOT and PSCz which all our models consistently over estimate. We interpret this to be due to a calibration error in the ISOPHOT data as a result of which the published flux density values are low by $``$30%. Also shown in the figure is the surface density of galaxies per redshift bin contributing to the counts at different wavelengths, as derived from the model, and the observed redshift distribution of the ISOCAM 15 $`\mu `$m galaxies in the HDFN+FF (Aussel et al., 2001). The size of each redshift bin is 0.2.

In the preceding three subsections, we have shown a range of models that evolve the local 15 $`\mu `$m luminosity function and fit the observed counts at mid- and far-infrared wavelengths as well as the spectrum of the CIRB. Ultradeep SIRTF observations at 24 $`\mu `$m can potentially break the degeneracy in these models if the counts can be determined to an accuracy of 20% or better (Figure 12). The range of integral source counts that we predict based on our three evolutionary models (density, luminosity and combination of both) are 4.7, 3.8 and 3.7 per $`\mathrm{}\mathrm{}`$ for $`S_\nu >120\mu `$Jy and 19.8, 19.3 and 11.8 per $`\mathrm{}\mathrm{}`$ for $`S_\nu >22\mu `$Jy respectively. However, the integral counts can be as low as 9.1/$`\mathrm{}\mathrm{}`$ at a flux density limit of 22 $`\mu `$Jy for models which are at the extreme lower limit of the uncertainty in current observations.

## 4 The Origin of the CIRB

### 4.1 Nature of the Galaxies Contributing to the CIRB

In our evolution models, we have assumed that the contribution from AGN to the counts and the cosmic background is insignificant. Other evolutionary models, which assumed an AGN component, arrived at the same conclusion (Malkan & Stecker, 1998; Rowan-Robinson, 2000; Xu et al., 2000; Franceschini et al., 2001). Observational evidence for this assumption comes from deep Chandra observations of the HDFN proper (Brandt et al., 2001) which detected 8 of the ISOCAM 15 $`\mu `$m sources. However, only one of these is an AGN at $`z`$1 and this was already known as such from observations at visible wavelengths (see discussion in Elbaz et al., 2001; Aussel et al., 2001). It should be noted that it is insufficient to have an AGN in a galaxy to violate this assumption but that the integrated infrared light of the galaxy must be dominated by an AGN rather than by star formation. However, since the contribution of dust obscured AGN at high redshift, beyond ISOCAM detection thresholds, is unknown, our results are subject to this uncertainty. A large ($`>`$ 20%) contribution from AGN to the source counts or the CIRB will imply a weaker redshift evolution of the luminosity function.

The evolution parameters in our model are constrained strongly by the ISOCAM counts at $`z<1.2`$ and by the SCUBA counts at $`z13`$. In addition, the ISOCAM counts are dominated by LIGs while the SCUBA counts are dominated by ULIGs. Thus, barring a dramatic change in the ratio of LIGs to ULIGs between a redshift of 1 and 2, we conclude that our models have robustly determined the evolution of the luminous end ($`L_{IR}>10^{11}L_{\mathrm{}}`$) of the LLF up to $`z2`$. At $`z>>2`$, the best constraint comes from the spectrum of the CIRB. Since all our models which are almost flat beyond $`z2`$ provide values for the CIRB that are at the upper limit of the values observed by FIRAS, we conclude that these models place a strong upper limit on the estimate of dust enshrouded star formation at $`z>2`$.

Our models indicate that about 80% of the 140 $`\mu `$m CIRB is produced at $`z<1.5`$. In comparison, 90% of the 15 $`\mu `$m EBL, 65% of the 240 $`\mu `$m background and only about 30% of the 850 $`\mu `$m background, is produced within this redshift range (Table 1). We also derive that ISOCAM galaxies brighter than 0.1 mJy contribute 14.5 nW m-2 sr-1 to the 140 $`\mu `$m EBL, while galaxies brighter than 0.05 mJy produce 16.8 nW m-2 sr-1. This accounts for about 75% of the total far-infrared background (Table 2). In comparison, we find that the contribution from SCUBA detected galaxies brighter than the confusion limit of 2 mJy at 850 $`\mu `$m is only 3.4 nW m-2 sr-1 at 140 $`\mu `$m while galaxies brighter than 0.5 mJy produce 16.4 nW m-2 sr-1. This is because at 850 $`\mu `$m flux densities fainter than 1 mJy, the LIGs which dominate the ISOCAM counts and produce the majority of the 140 $`\mu `$m EBL contribute significantly to the SCUBA counts. The total EBL at 15 $`\mu `$m from our model is 3.2 nW m-2 sr-1. The EBL obtained by integrating the observed ISOCAM counts above 50 $`\mu `$Jy is 2.4$`\pm `$0.5 nW m-2 sr-1 (Elbaz et al., 2001), hence as much as 73$`\pm `$15% of the 15 $`\mu `$m background might have already been resolved by ISOCAM. The models indicate that ULIGs contribute 15% of the 15 $`\mu `$m EBL observed by ISOCAM above 0.1 mJy, LIGs about 65% and normal and low luminosity starburst galaxies the balance. In comparison, at 140 $`\mu `$m, we find that ULIGs contribute 25% of the CIRB, LIGS contribute 60% and normal/SB galaxies the balance (Figure 13). Thus, infrared luminous galaxies which appear to be indistinguishable from normal galaxies in terms of their optical/near-infrared luminosity and form a negligible part of the energy budget in the local universe, dominate the star formation and therefore the energy budget at redshifts $`z13`$.

### 4.2 Infrared Luminosity Function and the LBG Connection

The first panel of Figure 14 shows the 15 $`\mu `$m local luminosity function (LLF) along with the fraction of galaxies that are evolved in the different models. For the density evolution model, the broken red line is evolved. For the luminosity evolution model, the whole 15 $`\mu `$m LLF shown as the solid black line, is evolved. For the model with density+luminosity evolution, the broken purple line is evolved. In the pure density evolution and density+luminosity evolution model, there is a non-evolving component with a constant comoving density which corresponds to the difference between the total LLF and the evolving component.

Although the fraction of the luminosity function that is evolving and the evolutionary parameters are significantly different in our three evolutionary scenarios, we find that our models predict similar comoving number densities of infrared luminous galaxies at high redshift. This is illustrated in two panels of Figure 14 which show the derived 15 $`\mu `$m luminosity function at redshifts of 0.4 and 0.8. The models also provide strong evidence for a change in the shape of the infrared luminosity function. The comoving number density of luminous infrared galaxies has to increase by more than two orders of magnitude between redshifts of 0 and 1, to fit the ISOCAM and SCUBA counts. The faint end of the LLF cannot be enhanced by the same factor since this would lead to an over-production of the CIRB although these galaxies would be below the sensitivity limit of the long wavelength surveys. Lastly, as much as 85% of the far-infrared background can be attributed to infrared luminous galaxies. This implies that the contribution from normal and low-luminosity starburst galaxies ($`L_{IR}<10^{11}L_{\mathrm{}}`$) to the dust enshrouded star formation rate is relatively small. So, estimates of the total SFR by applying a constant extinction correction to all optical/UV selected galaxies are incorrect. We conclude that long-wavelength surveys between 15 $`\mu `$m and 850 $`\mu `$m that probe galaxies at the luminous end of the infrared luminosity function provide a very effective way of tracing the bulk of the dust obscured star formation.

The connection between infrared luminous galaxies and the Lyman-break galaxy (LBG) population is intriguing. Figure 14 shows a comparison between the LBG 60 $`\mu `$m luminosity function at $`z3`$ of Adelberger & Steidel (2000, hereafter AS00), which was derived based on an extinction correction to optical/UV data as a function of the UV slope of individual galaxies, and our equivalent 60 $`\mu `$m luminosity function at $`z3`$, which we have argued earlier is only a strong upper limit. The agreement is extremely good considering that they were estimated in completely independent ways. The AS00 luminosity function predicts almost the same luminosity function as our estimate from the pure luminosity evolution model to within 50%. It is discrepant with the luminosity functions from our other two models by as much as an order of magnitude at the faint end but only a factor of 2 at the bright end. As mentioned earlier, the long wavelength surveys which constrain our models are mainly sensitive to the evolution of galaxies at the bright end of the luminosity function at $`z<2`$. We are unable to constrain with much certainty the evolution of the faint end of the luminosity function although we do place an upper limit based on the observed intensity of the CIRB. Furthermore, our estimates at $`z3`$ are only a strong upper limit to the number density of infrared luminous galaxies since any further evolution at high redshift overproduces the CIRB at $`\lambda >200`$ $`\mu `$m while a decay in the evolution at $`z>2`$ as $`[1+(z2)]^2`$ is marginally consistent with both the submillimeter counts and the CIRB spectrum. Thus, we conclude that optical/UV surveys that trace the LBG population at $`z>3`$, after an extinction correction factor that spans the range 2$``$100, provide a good estimate of dust obscured star-formation at high redshift. They complement the results of future mid- and far-infrared surveys with SIRTF which will be able to directly observe the dust emission of LBGs with LIG-type infrared luminosities up to $`z2.5`$.

There is observational evidence for the LBG population to be distinct from the bright ($`S_\nu (850\mu \text{m})>6`$ mJy) submillimeter galaxies (Barger, Cowie, & Richards, 2000; Chapman et al., 2000). This is because most of the bright 850 $`\mu `$m galaxies are extreme ULIGs with $`L_{IR}>10^{12.6}L_{\mathrm{}}`$. The AS00 observations detect only 2 out of 831 galaxies above this detection threshold and only 27 of their LBG sample would have a submillimeter detection above the level of $``$1 mJy. Secondly, at $`z3`$, 850 $`\mu `$m observations would probe rest frame $``$200 $`\mu `$m emission. For a given far-infrared luminosity, the 60 $`\mu `$m luminosity shows a factor of 2$``$3 less scatter than the 850 $`\mu `$m luminosity among local galaxies (Figure 3) suggesting that the $`\lambda >200`$ $`\mu `$m spectral shape of galaxies might potentially have a larger scatter which would lead to uncertain flux estimates on a galaxy by galaxy basis. Lastly, we do not find any scenario where the contribution from ULIGs is greater than 30% of the comoving star formation rate at $`z3`$. Naturally, the contribution from extreme ULIGs traced by the bright submillimeter galaxies is even smaller. Thus, although the contribution to the star formation rate density from extreme ULIGs is missed in observations of the LBG population, their contribution is significant only at the level of $`<`$10% and hence they are less important to an estimate of the high-redshift dust obscured star formation.

### 4.3 The Revised Star Formation History of the Universe

Having constrained the evolution of the mid-infrared and thereby the far-infrared luminosity functions, we can derive the evolution of the dust enshrouded star formation rate with redshift. Using the equations listed in Section 2, our derived comoving SFR from all galaxies is shown in Figure 15 while the separate contribution from LIGs and ULIGs is shown in Figure 16. Figure 15 also shows the absolute minimum and maximum range of dust obscured star formation rate values. The maximum values are derived from models which marginally overproduce the CIRB and the counts. The minimum values are derived by using an evolutionary model which is marginally consistent with the observations (Section 4.2) and by only considering the contribution from LIGs and ULIGs since those are the only galaxies which are directly observed at high redshift in the 15 $`\mu `$m and 850 $`\mu `$m surveys. Also shown is the SFR as inferred from direct observations at visible/UV/near-infrared wavelengths (Lilly et al., 1996; Madau et al., 1996; Connolly et al., 1997; Cowie, Songaila & Barger, 1999; Steidel et al., 1999; Yan et al., 1999), SFR estimates obtained from extinction corrections to these observations (Madau, Pozzetti, & Dickinson, 1998; Meurer et al., 1999; Steidel et al., 1999; Thompson, Weymann & Storrie-Lombardi, 2001), and SFR derived from ISOCAM observations of the Canada France Redshift Survey (CFRS) field (Flores et al., 1999). In addition, lower limits to the SFR from radio measurements and two points representing the completeness corrected submillimeter observations are also shown (Barger, Cowie, & Richards, 2000). The 1$`\sigma `$ uncertainty in our derived rate is about 50% and is primarily dependent on the transformation from 15 $`\mu `$m to infrared luminosities which as derived earlier has a 1$`\sigma `$ of 40%, and the transformation from infrared luminosities to star formation rate which assumes a Salpeter initial mass function (IMF) and also has an uncertainty of about 30% (Kennicutt, 1998).

The dust enshrouded SFR density peaks at a redshift of 0.8$`\pm `$0.1 with a value of 0.25$`{}_{0.1}{}^{+0.12}M_{\mathrm{}}^{}yr^1Mpc^3`$. In a (0.3,0.7,75) cosmology this corresponds to 6.2 Gyr after the Big Bang. The dusty SFR then remains almost constant up to $`z2`$ which corresponds to an age of 3 Gyr beyond which this value provides a strong upper limit to the amount of dust obscuration. This is similar to the shape of the star formation history preferred by Sadat, Guiderdoni & Silk (2001) in their analysis of the CIRB. Our values are a factor of 2 larger than estimates at $`z1`$ from $`H\alpha `$ observations by Yan et al. (1999) and a factor of 3-7 larger than extinction uncorrected optical/UV observations at $`z<2`$. The models are in excellent agreement with the submillimeter data corrected for incompleteness (Barger, Cowie, & Richards, 2000) but are higher than the Steidel et al. (1999) extinction corrected points. The values we derive are systematically higher than those in Gispert, Lagache & Puget (2000) but within the uncertainties especially if the difference in the $`L_{IR}`$-to-SFR calibration coefficient is factored in. Our models also indicate a faster evolution at $`z<1`$ than the models of Blain et al. (1999b) which is not surprising since they did not use the ISOCAM data to constrain their low redshift evolution. However, our high redshift plateau is similar to their ‘Anvil-10’ model.

Recently, Xu et al. (2000) has developed a multiparameter model where the 25 $`\mu `$m luminosity function of Shupe et al. (1998) is partitioned into three components - starburst, late-type galaxies and AGN and each component is evolved independently of the other. Specifically, they evolved the starburst population in luminosity as $`(1+z)^{4.2}`$ and in density as $`(1+z)^2`$ out to $`z=1.5`$. The late-type galaxy population was evolved in luminosity as $`(1+z)^{1.5}`$ while the galaxies with AGNs evolve in luminosity as $`(1+z)^{3.5}`$. Beyond $`z=1.5`$, all the components drop-off as $`(1+z)^3`$. Using a 25 $`\mu `$m to IR luminosity conversion based on IRAS data, we have converted their evolution for starburst and late-type galaxies into a star formation rate and compared it with ours. This is shown as the black dotted line in the upper plot of Figure 15. We find that their derived rates between redshifts of 1 and 2.5 are inconsistent with our models. The motivation for this peak is not clear since only the CIRB and the SCUBA counts place constraints on the evolution at this redshift range and both can be reproduced very well by an almost flat evolutionary history at $`z>0.8`$ (See Section 3). However, their evolution at $`z<1`$ agrees reasonably well with ours since both are principally constrained by the 15 $`\mu `$m ISOCAM counts. Furthermore, their decline in the SFR at high redshift ($`z>2`$) is similar but below our lower limit.

By integrating our comoving star formation rate density over redshift and thereby cosmic time, we can derive the density of stars and stellar remnants and compare it with the total baryon density in the local Universe. Stellar lifetimes were chosen for solar metallicity stars (Bressan et al., 1993) while the mass of remnants was chosen using the recipe of Prantzos & Silk (1999, and references therein). If a Salpeter IMF is assumed, then the model predicts a local density of baryons of about 1.0$`\times 10^9M_{\mathrm{}}Mpc^3`$ which is a factor of 2 in excess of the value of 5$`\pm `$3$`\times 10^8M_{\mathrm{}}Mpc^3`$ estimated by Fukugita, Hogan & Peebles (1998). The model also predicts that 100% of the local stars and remnants would have been produced at a redshift $`z<2.0`$. If we instead use the shape of the IMF below 1 M suggested by Gould, Bahcall & Flynn (1996) which reduces the number density of low mass stars, then the density of stars and remnants resulting from the model is 7.5$`\times 10^8M_{\mathrm{}}Mpc^3`$, in agreement with the local density. Madau & Pozzetti (2000) argued for a similar IMF based on their analysis of the total extragalactic background light. Interestingly, our model predicts that the local baryon density in stars and remnants, derived by integrating the star formation in ULIGs with a redshift distribution as shown in Figure 16, is similar to that seen in local spheroids, suggesting that high redshift infrared luminous galaxies may be the progenitors of present day spheroids.

## 5 Conclusions

A variety of observational data at mid-infrared through submillimeter wavelengths trace the fraction of emission from stars that is thermally reprocessed by dust. By using the counts of galaxies at these wavelengths, it is possible to estimate the amount of star formation that is enshrouded by optically thick HII regions and thereby invisible to observations at ultraviolet and visible wavelengths. In addition, the spectrum of the cosmic infrared background at mid- and far-infrared wavelengths places an upper limit on the fraction of starlight that has undergone thermal reprocessing by dust.

We have developed a set of template spectral energy distribution (SEDs) for galaxies as a function of infrared luminosity, which reproduce existing data at 0.44, 7, 12, 15, 25, 60, 100, and 850 $`\mu `$m from ISO, IRAS and SCUBA on nearby galaxies. The 15 $`\mu `$m local luminosity function was then evolved with redshift, considering both luminosity and density evolution models into account, and using the template SEDs to fit the observed counts at 15 $`\mu `$m, 90 $`\mu `$m, 170 $`\mu `$m and 850 $`\mu `$m. A number of evolutionary models provide reasonable fits to the data and the spectrum of the CIRB. The principal reason for this is that all the long wavelength surveys are typically sensitive to only the most luminous galaxies ($`L_{IR}>10^{11}L_{\mathrm{}}`$) at $`z>0.5`$. So, evolutionary models that result in similar luminosity functions at $`L_{IR}>10^{11}L_{\mathrm{}}`$ are degenerate. However, our models accurately constrain the comoving number density of these luminous galaxies as a function of redshift. In the local universe, it is these galaxies, many of which show morphological signatures of interaction, that show an infrared-determined star formation rate (SFR) that is about an order of magnitude higher than the corresponding UV-determined SFR. By integrating the infrared luminosity of these luminous galaxies we then obtain an estimate of the dust enshrouded star formation rate. The dust enshrouded star formation rate density appears to peak at a much lower redshift than previously thought, at $`z=0.8\pm 0.1`$ with a value of 0.25$`{}_{0.1}{}^{+0.12}M_{\mathrm{}}^{}yr^1Mpc^3`$ and remains approximately constant at least till $`z2`$. Any drop off at a lower redshift would result in an underestimate of the 850 $`\mu `$m galaxy counts. Although our models do not constrain the evolution of the faint-end ($`L_{IR}<10^{10}L_{\mathrm{}}`$) of the luminosity function, their net contribution to the high-redshift dust enshrouded star-formation is negligible as can be seen in the range of evolutionary models considered. The evolution at $`z>2`$ is constrained much more weakly. Having a constant SFR between redshifts of 0.8 and 4 is consistent with the CIRB spectrum and the submillimeter counts as is a decay by a factor of 7 between redshift 2 and $``$5. However, we find that there is excellent agreement between our luminosity function and the infrared luminosity function derived from extinction correction to optical/UV observations of Lyman-break galaxies at $`z3`$. This suggests that dust obscuration is significant even at $`z>3`$ and that the dust enshrouded star formation rate is constant to within a factor of 2 between redshifts 2 and 4.

The models also provide a census of the luminosity of galaxies that contribute to the counts at different wavelengths, their redshift distribution and the relative contribution to the cosmic infrared background at $`\lambda >5`$ $`\mu `$m as a function of redshift. Furthermore, we find that ultradeep observations with SIRTF at 24 $`\mu `$m down to a sensitivity of 25 $`\mu `$Jy can potentially break the degeneracy in the evolutionary models by detecting galaxies with $`L_{IR}10^{11.5}L_{\mathrm{}}`$ out to $`z2.5`$, well beyond the turnover redshift of 0.8 that is derived from our models.

RC wishes to thank Harland Epps and Rodger Thompson for kindly funding this research through NASA grant NAG5-3042. DE wishes to thank the American Astronomical Society for its support through the Chretien International Research Grant and Joel Primack and David Koo for supporting his research through NASA grants NAG5-8218 and NAG5-3507. We wish to acknowledge Pierre Chanial for collating published data from a large number of surveys and making them available to us.