Chapter 2
Modeling the Universe

To simulate a population of high redshift starburst galaxies I model their spectral energy distribution from the MIR to radio frequencies. I then describe their luminosity distribution as a function of redshift, based on the well studied luminosity distribution of FIR bright galaxies in the local universe, and based on the redshift distribution of (sub)mm sources derived from their spectral index. The parameters describing the model are adjusted such that they reproduce the extragalactic background from mm to MIR wavelengths, the (sub)mm source number counts and their redshift distribution. This is done by conducting Monte Carlo simulations of a population of starburst galaxies. Finally I model the spatial distribution of the simulated sources based on numerical simulations of the dark matter distribution in order to generate artificial skymaps.

2.1 A Model SED

The SEDs of the starburst galaxies play a crucial role. The SED is used to relate flux density to bolometric luminosity for the SCUBA and MAMBO sources. Since SEDs of starburst galaxies are of such importance, they are subject to many publications by various authors. The main focus is usually on the dominating dust emission part of the SED around 100 m, as it governs the flux densities at submm wavelengths. Even more, the bolometric luminosity is mainly produced by this part of the SED.

Figure 2.1: The model fit to the SED of Arp 220. Light blue: The cold dust component with Tc = 20 K and = 2. Green: The warm dust component with Tw = 49 K and = 2. The mass ratio is R = 22. Yellow: The hot dust component with Th = 200 K and = 2 plus PAH line features. Pink: Synchrotron radiation -0.8. Dark blue: Free-free radiation -0.1. Red thick line: The model, including extinction with d = 1.15 kg m-2. The continuous black lines in the original SED belong to ISO spectral data and are not properly calibrated, so that they need not be fitted by the model. No attempt was made to reproduce the optical or NIR spectrum.

I model the SEDs of starburst galaxies from the MIR to radio frequencies by adding five components, as illustrated in Fig. 2.1. These components are three dust components (including extinction) described by greybodies, plus synchrotron and free-free radiation described by power laws. For the properties that describe these components I find relations and distributions from observations of starburst galaxies in the local universe.

The dust emission of the dust enshrouded star forming regions is best described by modified blackbodies, so called greybodies. I model the FIR-part of the SED through two greybodies. In the mid infrared (MIR) a third dust component is added. I include extinction as well as absorption and emission line features in the model. At radio frequencies, I model synchrotron and free-free radiation through power laws. The radio flux density is closely related to the FIR luminosity and can be calculated from the relation Carilli and Yun found, which is described in section 1.5.

2.1.1 Radiative Transfer

In the description of the SED model radiative transfer is important. The fundamental equations of radiative transfer are described here to get proper definitions of all quantities involved.

We consider a cylinder of length dr and a cross section dA. I is the radiation intensity (in W m^-2 Hz^-1 sr^-1) perpendicular to the bottom area going into solid angle d. Emission and absorption along dr result in a change of the intensity by an amount dI. We define the emission coefficient and the absorption coefficient ':

(2.6)

and the optical thickness

(2.7)

leads to

(2.8)

This is the basic equation of radiative transfer. It is obvious that if I < G then dI/d > 0, meaning that the intensity increases in the direction of propagation, and vice versa for I > G. In the case of thermal equilibrium the same amount of energy is absorbed as is emitted, so that dI/d = 0. This results in Kirchhoff’s law:

(2.9)

In that case the source function is given by Planck’s law:

(2.10)

where h is Planck’s constant, T is the temperature of the blackbody, and k is Boltzmann’s constant.

A formal solution of eq. (2.8) is

(2.11)

where I(0) is the background radiation intensity, coming from behind the regarded region. Assuming that this background radiation can be ignored and that the source function is constant along the radiation path, one obtains

(2.12)

On the other hand, to describe extinction by a medium with negligible emission ( = 0) in front of an object emitting I(0) one must use

(2.14)

If is independent of r one finds

(2.15)

so that ' has the dimension of an inverse length. ' is related to the cross section of the scattering particles for photons. In case of free electrons this is the cross section of Thomson scattering, _el, in case of atoms or molecules that of Rayleigh scattering, _at. The absorption coefficient is then given by

(2.16)

where n is the number of scattering particles per volume.

A quantity often used is the mass absorption coefficient

(2.17)

where is the matter density in mass per volume. is the cross section per mass and has the dimension of area per mass. To calculate from (assuming the latter being independent of r), one must use

(2.18)

where is the mass of one particle and

(2.19)

is the column density. Finally we define the mass column

(2.20)

2.1.2 The Absorption Coefficient of Dust

Draine & Lee (1984) examined the properties of dust grains. From laboratory data they derived models for the cross section and absorption coefficient for graphite and silicate dust grains. They constrained their models by observations of various astronomical objects (see Fig. 2.2). Bertoldi et al. (1999) focussed on the absorption features between 2 and 50 m and found a parameterized form to describe the most prominent features at 3.05 (ice), 9.66 and 19 m (silicate). The millimeter and submillimeter properties of dust were examined, e.g., by Schwartz (1982) and the absorption coefficient was found to obey a power law in frequency, .

In my model I use the parameterized form from Bertoldi et al. (1999) for wavelengths shorter than ~ 25 m:

(2.22)

Both parts are connected smoothly:

Figure 2.2: The extinction curve for the dust model Draine and Lee examined. The components that contribute to the total extinction curve are shown. The extinction has been multiplied by the wavelength for a better plotting range. Taken from Draine & Lee (1984).

Figure 2.3: The absorption coefficient used in the model ( = 2).

2.1.3 Dust Emission

To describe a dust enshrouded star forming region, the dust is assumed to be in thermodynamic equilibrium, heated by the starburst and cooled by emission of radiation, so that the source function is a blackbody. The dust is mostly optically thin ( < 1) in the FIR, so that eq. (2.12) gives

Since usually the angular size of the high redshift starburst galaxies for FIR to radio observations is smaller than the beam size, the flux density received is that of the entire galaxy. The flux density emitted by a galaxy containing i dust enshrouded starburst regions is:

(2.28)

where S_' is is the flux density observed at frequency '. is the corresponding rest frame frequency.

(2.29)

is the dust mass of the ith component and A_i the area of the ith dust region, perpendicular to the line of sight. D_L is the luminosity distance given by eq. (A.26), and

(2.30)

is the mass of dust enshrouding the starbursts. The flux density of a single dust enshrouded starburst is described by a greybody:

(2.31)

2.2 The FIR Spectrum

Table 2.1: The galaxies examined by Dunne & Eales (2001). S12 is taken from NED. cz Tw Tc R log(MD/M) S60 S12 Name [km s]-1  [K] [K] M c/Mw [Jy] [Jy] UGC 903 2518 45 21 91 7.37 7.91 0.357 NGC 958 5738 44 20 186 8.28 5.9 0.59 UGC 2369 9400 36 21 6 8.08 7.68 0.22 UGC 2403 4161 50 22 100 7.57 7.51 0.32 NGC 1614 4778 38 20 3 7.78 33.12 1.44 NGC 1667 4547 28 17 3 7.91 6.24 0.59 NGC 2856 2638 45 22 58 7.08 6.15 0.34 NGC 2990 3088 42 21 58 7.33 5.49 0.29 UGC 5376 2050 44 21 87 7.11 5.94 0.31 ARP 148 10350 37 21 12 8.28 6.95 0.12 MCG+00-29-023 7464 36 20 13 7.99 5.4 0.36 ZW 247.020 7666 39 25 6 7.50 5.91 0.15 1 ZW 107 11946 39 20 7 8.23 9.15 0.22 IR 1525+36 16009 45 19 15 8.29 7.2 - ARP 220 5452 48 18 42 8.80 103.33 0.64 NGC 5962 1963 31 18 11 7.49 8.99 0.72 NGC 6052 4712 37 20 16 7.65 6.46 0.26 NGC 6181 2379 31 19 7 7.47 9.35 0.65 NGC 7541 2665 31 17 6 7.90 20.59 0.149 NGC 520 2281 50 24 61 7.46 31.55 0.91 UGC 2982 5305 31 17 6 8.11 8.7 0.57 NGC 2623 5535 50 27 18 7.59 25.72 0.24 NGC 3110 5048 42 23 44 7.93 11.68 0.585 IR 1017+08 14390 40 19 8 8.20 6.08 - IR 1056+24 12912 40 26 6 8.14 12.53 0.96 IR 1211+03 21703 49 24 27 8.55 8.39 - NGC 4418 2179 58 21 62 7.40 42.32 - UGC 8387 7000 33 24 2 7.92 13.69 0.26 ZW 049.057 3927 44 23 28 7.74 21.06 - NGC 7592 7350 39 19 20 8.13 8.02 0.282 NGC 7679 5138 37 19 14 7.75 7.28 0.52 NGC 7714 2798 42 21 13 7.05 10.52 0.47 cz: Recession velocity; Tw: Temperature of the warm component; Tc: Temperature of the cold component; R: Mass ratio cold to warm component; MD: Total dust mass; S60, S12: Flux densities at 60 and 12 m. M = 1.989 . 1030 kg is the mass of the sun.

Dunne & Eales (2001) fitted the FIR part of the SEDs of 32 local starburst galaxies (Table 2.1) by greybodies. They found that a single dust component did not fit the measured 450 m flux densities. They showed that a fit with two components works much better. Their model had the following parameters:

• T_w, the temperature of the warm dust component.
• T_c, the temperature of the cold dust component.
• M_D, the total dust mass.
• R = M_c/M_w, the ratio of dust masses cold to warm.
• , the absorption coefficient power law index for both components.

Dunne and Eales found that varies slightly between 1.5 and 2 and then modeled the SEDs through greybodies with a fixed = 2. They found that the temperature of the cold component varies only slightly between 17 and 27 K, while the Temperature of the warmer dust component showed a wider scatter between 28 and 58 K. In my SED model I take T_c to be fixed at 20 K and = 2. For the remaining parameters the following distributions and relations were found and adopted

• The logarithm of the ratio R = M_c/M_w is distributed uniformly between 0 and 2 and is not strongly correlated with the 60 m luminosity (Fig. 2.4). The definition for the 60 m luminosity used is that of Dunne and Eales:

  (2.32)

  where S_60m(1+z) is the flux density observed at 60 (1 + z) m in W m^-2 Hz^-1. This is not a proper luminosity, but I use this quantity for the correlations, following Dunne and Eales.

• I found that the warm component temperature T_w is related to R by

  (2.33)

  with a scatter of = 0.06 (see Fig. 2.5).

The FIR luminosity of the two dust components described by given values T_c, T_w,  and R is proportional to the total dustmass M_D. Its value is scaled so that the FIR SED has the required luminosity.

Figure 2.4: The mass ratio R plotted against LD (defined by eq. (2.32)). No strong correlation is evident. I shall assume that log R is uniformly distributed relative to LD.

Figure 2.5: The temperature of the warm dust component Tw plotted against the mass ratio R. Linear regression in the logarithmic plot was used to derive the relation between Tw and R.

2.2.1 Mid-Infrared

Dale et al. (2001) present an alternate method to assemble model SEDs. They model single regions of a galaxy. Their model is parameterized by the heating intensity

(2.34)

where T_w is the temperature of the warm component defined in section 2.1.3. The final SED is assembled from single dust regions according to a distribution assigned to U. Although I do not follow this approach, the MIR features of the Dale SEDs are used and modified for the model SEDs developed here.

Most of the MIR luminosity is given by a greybody with a fixed temperature of T_h = 200 K and = 2. The dust mass of this hot component is 5 ^. 10^-6 times the dust mass from the FIR SED. This ratio I derived from fitting the model to Arp 220 and M82. The dust of this greybody is assumed to consist mostly of PAHs. Since PAHs are destroyed at high temperatures, Dale et al. introduce a PAH damping factor

(2.35)

by which the flux density from the hot component is decreased.

Figure 2.6: The PAH damping factor D. The curve follows a Fermi-Dirac profile (see eq. (2.35)).

Boulanger et al. (1998) fitted the narrow MIR emission features, which are believed to arise from PAHs, seen in ISOCAM spectra of two objects (NGC 7023 and the Ophiuchus molecular cloud) by Lorentz functions. The Lorentz function is:

(2.36)

where ₀ is the center frequency, is the width parameter and A is the normalization, resulting in a peak amplitude of A/(). The full width at half maximum (FWHM) is 2. The fitting parameters Boulanger et al. found and that I adopt are shown in table 2.2. The fit itself can be viewed in Fig. 2.7.

Table 2.2: The Lorentz fit parameters derived by Boulanger et al. (1998) for the probable PAH features. Center wavelength [m] 6.2 7.7 8.6 11.3 12.7 15.0 A 4.3 13.4 2.3 4.0 4.0 1.0 FWHM [cm-1] 70 119 57 34 73 35

Figure 2.7: The Lorentz-fit Boulanger et al. (1998) produced for the features in the spectrum of NGC 7023.

The parameter A is somewhat arbitrary, as only the relative peak sizes are important. The Lorentz functions representing the lines are multiplied by a factor of 6 ^. 10¹¹ (which I derived by fitting the features to those seen in Arp 220 and M82) and scaled by the flux density of the hot component.

2.2.2 Extinction

The strong 10 m absorption features in Arp 220 and M82 indicate extinction by dust along the line of sight. This dust is assumed to have the same characteristics as the dust heated by the starburst. I apply the absorption coefficient described in section 2.1.2 for this dust. Eq. (2.13) is used to calculate the apparent SED of a galaxy:

Figure 2.8: The number of objects plotted against column density times particle mass. The result suggests a normal distribution.

To get a distribution of _d in starburst galaxies, I collected the 12 m flux densities from the NASA extragalactic database (NED) for the galaxies that Dunne and Eales examined in their study (Table 2.1). The model parameters Dunne and Eales found for those galaxies were applied and the MIR component added. Then I modified _d to fit the 12 m flux density S₁₂. Fig. 2.8 shows the resulting distribution of _d for the 27 objects where I could find 12 m flux densities. For the model a normal distribution is assumed:

(2.38)

2.2.3 Radio Spectrum

By fitting the FIR SED through greybodies, one can calculate the 350 GHz flux density. The flux density at 1.4 GHz is given by the relation found by Carilli and Yun. In my model, the radio emission consists of free-free and synchrotron radiation. Following Condon (1992), a power law is assigned to both types of radiation, using different powers: synchrotron radiation ^-0.8, and free-free radiation ^-0.1. Further, the ratio between both at 1.4 GHz is 13:1 (synchrotron:free-free) (Condon, 1992). So the radio part of the SED model is:

(2.39)

where S_1.4GHz is the flux density obtained by the Carilli and Yun relation.

2.2.4 Examples

Fig. 2.1 shows the model adjusted to fit the SED of Arp 220. Fig. 2.9 shows some realizations of the model for a given 60 m luminosity L₆₀ = 10¹¹ L . The model parameters are drawn randomly from the distributions I derived above.

Assembling all components according to the derived distributions and relations results in a SED model that reproduces the SEDs of local starburst galaxies well between ~ 0.5 m and ~ 4 ^. 10⁴ GHz. At higher frequencies further emission components and line features must be considered, at lower frequencies free-free absorption affects the radio emission (Condon, 1992).

Figure 2.9: Some realizations of model SEDs with L60 = 1011 L . The model parameters are drawn randomly from the previously described distributions.

2.3 The Luminosity Distribution

I describe the luminosity distribution of starburst galaxies by a luminosity function. Based on the well observed 60 m luminosity distribution of the local universe I develop the luminosity function for starburst galaxies at high redshift. To adjust the parameters of the luminosity function I conduct Monte Carlo simulations and compare the results to the source number counts and redshift distribution derived from the (sub)mm observations. The free parameters are adjusted to match the observational constraints.

The luminosity function

(2.40)

describes the Number of objects dN per luminosity interval dL and per comoving volume element dV (as defined in eq. (A.31)). The number of objects brighter than L per solid angle is

(2.41)

When the distances to the objects are not well known, one may only be able to measure a flux limited distribution:

(2.42)

where L(S,z) is a function that relates the observed flux density, S, to the luminosity of the object. This relation is given by eq. (A.28), involving the SEDs of the objects.

The luminosity distribution changes in time as the universe evolves. The luminosity function should reflect this evolution, resulting in a explicit dependence on redshift. The cumulative counts of objects observed per solid angle are then given by

(2.43)

2.3.1 The Local Universe

The luminosity distribution of the SCUBA and MAMBO sources is only constrained by measurements at high luminosities (L_bol > _~ 10⁹ L ). To model the distribution of fainter starburst galaxies at high redshifts, I base my luminosity function on the observed luminosity distribution of FIR galaxies in the local universe. Saunders et al. (1990) examined this luminosity distribution. Their study is mainly based on the Infrared Astronomical Satellite (IRAS) 60 m survey of galaxies, supplemented by several smaller samples from ground based observations to minimize selection effects due to the flux limit. The redshifts of all these objects are well known. A total of 2818 galaxies were taken into account, the majority of which are at redshifts below z = 0.1. Their luminosity function is parameterized by the 60 m monochromatic luminosity

'

The 60 m luminosity function Saunders et al. found is shown in Fig. 2.10. It consists of a flat power law for low luminosities and an exponential cutoff at luminosities above L_*:

Figure 2.10: The luminosity function Saunders et al. (1990) derived from the IRAS 60 m survey of local galaxies.

2.3.2 Modification at High Luminosities

As the Saunders luminosity function is parameterized in L₆₀, the flux density limited number counts at other wavelengths than 60 m can only be computed by relating the flux density at the required wavelength to that at 60 m via the SEDs of the objects. In the simulations that I conduct to adjust the parameters of the luminosity function the SED of every object is known, as it is given by the model described in section 2.1. Having this knowledge, the source number counts that are flux density limited at 1200 m or 850 m can be calculated and can be compared with the observational results. At luminosities above L_** 10¹² L a second cutoff had to be introduced in order to match the observational results. The modified luminosity function is:

(2.47)

The second cutoff is not critical regarding the 60 m galaxy number counts. There are very few objects with 60 m luminosities above 10¹² L in the surveys Saunders et al. analyzed so that the Saunders luminosity function is not well constrained at these luminosities.

2.3.3 Evolution

To reflect the formation and evolution of galaxies a redshift dependence must be applied to the luminosity function. To find this dependence, the model is adjusted to match the redshift distribution of the (sub)mm deep field sources derived by the spectral index. This redshift distribution suggests an increase of source numbers with redshift to at least z 2.5. This effect is convolved with an evolution of the comoving volume per solid angle, which first increases with redshift, but then decreases again as the universe shrinks with redshift. The comoving volume per solid angle has a peak at z 2.4 in the lambda cosmology used here.

To get the evolution for very high redshifts we look at the mechanism how galaxies formed in the early universe. It is not understood and subject to many discussion. The most favored model is however, that first small density fluctuations in the initial matter distribution in the universe collapsed, and ever larger ones collapsed subsequently. Today we can still see the ”footprints“ of those initial fluctuations as anisotropies in the cosmic microwave background (CMB). Unlike dark matter, baryonic matter can dissipate energy through cooling and therefore it collapses into the centers of the dark matter condensations. This is the place where galaxies first form. By this mechanism first small galaxies form, and then later bigger galaxies form through merging of smaller ones. Thus, our luminosity function should predict only few bright objects at very high redshifts, but the number must increase as redshift decreases. The evolution then must have a peak, which the redshift distribution of (sub)mm sources suggests to be at z 2.5.

To implement the evolution into the luminosity function one may represent pure density evolution by making C dependent on z. This would mean that the relative distribution of objects according to their luminosity would not vary with epoch, but the absolute number per comoving volume would change with redshift. From the observations we cannot see whether there is a density evolution since we would need to count all galaxies.

The way chosen here is to make L_* dependent on z, resulting in a change in the relative luminosity distribution of the galaxies with redshift and in a change in absolute numbers with redshift as well. This reflects the change in the intrinsic properties of the galaxies with epoch. The evolution is described by

(2.48)

L_** is assumed not to evolve.

The parameters b and c are adjusted so that the model reproduces the redshift distribution of (sub)mm sources derived from the observations by the spectral index. The introduced evolution affects the luminosity function for the local universe, and thus the 60 m number counts. With current instruments we can only observe nearby objects at 60 m so that any evolution at higher redshifts does not affect these number counts. Saunders et al. found a pure strong density evolution, C(z) (1 + z)^7±2 to fit the 60 m number counts, but they were also able to fit the counts by making L_*(z) (1 + z)^3±1. (Blain et al., 1999) modeled the number counts at different wavelengths using a greybody as template SED and found the best fits using L_*(z) (1 + z)^3.9±0.2 in the local universe. The slope of the ansatz for evolution used in the simulation here is within this range at very low redshifts, L_*(z) (1 + z)⁴, for the final set of parameters, so that it affects the luminosity function in the correct way.

2.3.4 The Final Parameter Set

The parameters that describe the luminosity distribution were adjusted by Monte Carlo simulations. The source number counts, the source redshift distribution and the extragalactic FIR background were computed from these simulations and compared to the according observations that are described in chapter 1. The set of parameters that was found to give the best results is:

Figure 2.11: The luminosity function and its evolution. The thick black line is the Saunders luminosity function with the final parameter set, but without the additional cutoff. The thick blue line is the modeled luminosity function with the final parameter set. The black thin lines are the luminosity functions for rising z, but below the evolution peak. The peak function is the red line. The green lines show the evolution for redshifts beyond the peak. The crosses indicate the value of L*, the vertical dashed line that of L**

2.4 Agreement with Observations

The simulation must reproduce the observations presented in chapter 1. In this section these observations are compared to the results from Monte Carlo simulations.

2.4.1 Source Number Counts

Figure 2.12 shows a one steradian Monte Carlo simulation and the cumulative counts from IRAS data, taken from Tan et al. (1999). At high flux densities, the statistic of the simulated counts gets poor. Figure 2.13 shows the observed (sub)mm counts introduced in section 1.3 and the result of a 1 deg² simulation. The number counts can be reproduced by the simulations.

Table 2.3: Background radiation limits and simulation results. Wavelength Fixsen et al. (1998) Simulation m W m-2 sr-1 W m-2 sr-1 1200 1.7 . 10-10 1.7 . 10-10 850 5.0 . 10-10 5.3 . 10-10 Wavelength: Observing wavelength; Fixsen et al. (1998): Background intensity computed from the fit by Fixsen et al.; Simulation: The background intensity from the simulation.

Figure 2.12: The 60 m cumulative counts from a 1 sr simulation. The data points are from IRAS counts, the solid line represents the simulation result. The dashed lines indicate the 1 error. Figure 2.13: The submm cumulative counts introduced in section 1.3 compared with a 1 deg2 simulation. The MAMBO (red) and the SCUBA (black) deep field counts are shown, while the solid lines represent the simulation results. The dashed lines indicate the 1 error.

2.4.2 Source Redshift Distribution

The redshift distribution of sources derived from the (sub)mm deep field maps from the spectral index redshift correlation is shown in Fig. 2.14, together with the result of a 1 deg² simulation. The counts in this plot are scaled to the maximum of the number counts of sources with radio detections. The absolute counts depend on the flux limit assigned to the imaginary detector.

The redshift distribution is matched by the simulation for redshifts up to ~ 2.5. Beyond this redshift the distribution is not well constrained, but it is likely that the sources with only lower redshift limits will fill up the apparent gap in the redshift distribution derived from observations. Deeper VLA-maps and new methods of redshift determination will provide for the missing constraints in the future.

2.4.3 FIR Background

Figure 2.15 shows the extragalactic FIR background radiation, both the fit to the FIRAS measurements from Fixsen et al. (1998) and the results from a 1 deg² simulation. Table 2.3 lists the values of the fit from Fixsen et al. and the simulation results for the MAMBO and SCUBA wavelengths. The (sub)mm background is matched well, but the simulation gives too low background intensities at shorter wavelengths, although they are within the errors of the fit. This deviation may be reduced by an alternate evolution of the luminosity function.

Figure 2.14: The redshift distribution of sources. Blue lines are data from the deep field maps, the red line is the simulated distribution for the wavelength SCUBA operates at, and green is that for MAMBO’s wavelength. The black line is the distribution of sources brighter than 1011 L . A 1 deg2 simulation was done for this plot. The counts are scaled to the maximum of the number counts of sources with radio detections.

Figure 2.15: The extragalactic FIR background. The black line represents the background derived from a 1 deg2 simulation. The red lines show the fit from Fixsen et al. (1998) and the 1 error of it. The blue arrow represents the lower limit computed from the IRAS surveys.

2.5 Star Formation

Now that the simulation is found to be consistent with the observations, an average star formation rate density can be derived from the simulated population of high redshift star forming galaxies. The FIR luminosity of these objects comes from massive bright stars that form continuously in the starburst regions. These stars have relatively short lifetimes of about 3-30 million years, so that on a timescale comparable to this lifetime an equilibrium between star formation and death establishes. This equilibrium results in a proportionality of the bolometric luminosity L_bol and the star formation rate:

(2.49)

where can be calculated from the stellar mass function, i.e., the mass composition of the stellar population of a starburst, depending on time. Omont et al. (2001) used several models for the mass function. The most probable average factor of proportionality they found is = 1.5. Adopting this for the simulations, and assuming that L_bol L_FIR, the star formation rate density, i.e., the star formation rate per comoving volume, can be derived.

Several authors evaluated star formation rate densities from optical and ultraviolet (UV) observations. This is done by relating the (rest frame) ultraviolet brightness of the observed galaxies to their star formation rate. Up to redshifts of about two this was done, e.g., by Lilly et al. (1996) and Connolly et al. (1997). Beyond this redshift Madau et al. (1996) and Steidel et al. (1999) derived average star formation rate densities. The results from their analysis are shown in Fig. 2.16 as dashed data points.

The first results for average star formation rate densities derived from SCUBA data indicated a much larger value at high redshifts. This result indicated that in the UV extinction due to dust can not be neglected. The optical observations and the derived average star formation rate densities have to be corrected for this extinction. Steidel et al. (1999) estimated this extinction and derived extinction corrected average star formation rate densities for high redshifts. They derived corrections from the slope of UV spectra that should redden with increasing extinction. The corrections Steidel et al. found are a factor of 2.7 for redshifts below 2, and 4.7 above. The corrected data are shown as solid data points in Fig. 2.16. These corrections are however highly uncertain, as the method to derive them from reddened UV spectra does not work well with luminous starburst galaxies.

From the simulations the average star formation rate density was derived. The luminosities were calculated from the known SEDs of the objects and related to the star formation rate by eq. 2.49. Fig. 2.16 shows the result of a 1 deg² simulation. Without a priori aiming to match those values, the simulation returned an average star formation rate density that is consistent with the results from optical studies.

The simulations are based on the extrapolation of the observed source number counts, which only deliver constraints at high luminosities, L_bol > _~ 10⁹ L . The starburst galaxies of intermediate luminosities contribute a significant part to the average star formation rate density derived from the model. Since their modeled number results from extrapolation, uncertainties are inflicted. The results from optical surveys are also uncertain due to the dust extinction. The question what kind of objects dominate the star formation in the early universe is still not answered.

Figure 2.16: The star formation rate density derived from the bolometric luminosities of the simulated objects of an 1 deg2 simulation. The single data points are derived from optical measurements (see text) and are adjusted to the lambda cosmology used throughout this work.

2.6 Spatial Distribution

In order to study the effects of galaxy clustering on (sub)mm observations I generate artificial maps from Monte Carlo realizations of a starburst galaxy population. Therefore it is necessary to model the spatial distribution of the objects. I derive this distribution based on numerical simulations of the dark matter distribution conducted within the ”Hubble Volume project“ (Evrard et al., 2001).

The spatial distribution of galaxies is usually described by correlation functions. The spatial correlations between galaxies are well observed in the local universe. At high redshifts however, the spatial distribution of galaxies is essentially unknown. A difficult technical problem would be to build a Monte Carlo realization from correlation functions. To derive a spatial distribution I instead use numerical simulations of the dark matter distribution, conducted for the ”Hubble Volume project“. These are based on N-body simulations of the gravitational collapse of fluctuations in the initial dark matter distribution that are assumed to be correlated to the observed anisotropies in the CMB.

Figure 2.17: This picture shows results of numerical computations of the dark matter distribution. It displays a cut through a lightcone of 10 degrees opening angle from z = 0 (bottom) to z ~ 3.2 (top). The color represents dark matter density (the lighter, the denser). Taken from Evrard et al. (2001).

The spatial distribution of baryonic matter is related to that of dark matter. Since baryonic matter has more ways to dissipate energy than dark matter does, it cools, which results in a concentration of baryonic matter in the center of regions of dark matter overdensity. Galaxies first form where the baryonic matter is densest. The spatial distribution of galaxies is therefore a biased indicator for the dark matter distribution.

The average spatial distribution both of dark matter and of galaxies may be described by an overdensity

(2.50)

where (x) is the density field and is the medium density. A simple way to model the correlation of the spatial distribution of galaxies and dark matter is ”density biasing“. The density contrast of galaxies is here assumed to be linearly related to that of dark matter:

(2.51)

where b is the biasing parameter. I use this simple ansatz to generate artificial maps, since this introduces only one additional free parameter. A problem of this simple description is that the relative galaxy density contrast may have values below -1, which is unphysical. I assume the density contrast to be -1 in this cases.

The biasing parameter may be varied to study its influence on the simulated maps. Dekel & Lahav (1999) studied numerical simulations that combined dark and baryonic matter. From these simulations they derived a biasing parameter of b = 2.42 at z = 1. I use this result as an orientation for realistic values of the biasing parameter. I do not use those combined simulations since they introduce additional and rather parameters to my simulations.

2.7 Simulated Maps

Figure 2.18 shows a 1 deg² artificial map based on a simulation using pure density biasing with b = 2.42. The flux from the simulated sources is convolved with a Gaussian point spread function with 10.7 arcsec full width at half maximum (FWHM). The observing frequency of the imaginary detector is 1.2 mm. These properties reflect those of MAMBO at the IRAM 30 m telescope. No noise was added.

Figure 2.19 presents a similar skymap, but in addition to density biasing I apply a stringent luminosity biasing: The positions at the densest spots of the underlying dark matter distribution are occupied by the most luminous objects. The difference to the previous map is obvious.

Figure 2.18: A 1 deg2 simulated skymap. The ”observing wavelength“ is 1.2 mm, the FWHM beam width is 10.7 arcsec. Pure density biasing is used.

Figure 2.19: A 1 deg2 simulated skymap at an observing wavelength of 1.2 mm with a beam width of 10.7 arcsec FWHM. Galaxies nare placed spatially assuming density and luminosity biasing.

As no noise was added to the maps they represent an ideal observation. To obtain artificial maps that represent observations under realistic conditions, Frank Bertoldi used the actual MAMBO observations of the Lockman Hole and simulated an observation of the simulated source distribution shown in Fig. 2.18. Gaussian noise was added to a level similar to that of actual observations. The result is shown in Fig. 2.20 and 2.21. Negative double beam patterns generated by the differential observing method can be seen.

Figure 2.20: A simulated double beam observation of the modeled galaxy population. The observing pattern of the actual MAMBO Lockman Hole observation (Fig. 2.23) was used. Figure 2.21: Gaussian noise similar to that in the MAMBO Lockman Hole observed map.

By adding the noise map to the simulated double beam observation an artificial map of the modeled high redshift starburst galaxies observed under realistic conditions is constructed. It is shown in Fig. 2.22. The true observation of the Lockman Hole is shown in Fig. 2.23. Both maps may now be compared for noise, source distribution etc., and thus the influence of unresolved background sources and the effect of galaxy clustering may be studied. From these studies it might be possible to learn about the clustering of starburst galaxies at high redshifts.

Figure 2.22: An artificial map generated by adding the noise map of Fig. 2.21 to the simulated double beam map (Figure 2.20). This map mimics an observation of the simulated galaxy population under realistic conditions.

Figure 2.23: The MAMBO observation of the Lockman Hole.

2.8 Source Confusion

The flux from those high redshift starburst galaxies which account for the bulk of the FIR to mm background is below the detection limit of current bolometer cameras. Still the flux density of these fainter galaxies contributes ”noise“ to the measured maps. This ”contamination“ leads to source confusion, which inflicts errors to the flux densities and positions of those galaxies which are well detected in the (sub)mm maps. The magnitude of these errors depends on the observing beam width, as a smaller beam includes fewer background sources. The confusion error also depends on the flux density of the detected source, since for weak sources the contribution from the fainter galaxies is relatively stronger than for bright sources. Galaxy clustering is also affecting confusion if the unresolved sources are clustered around the detected source.

To estimate the confusion error from the artificial maps, sources were extracted by Gaussian fitting. We have done this with simulated maps from a simulation including density galaxy biasing, with no added noise, corresponding to an ideal observation. For a better understanding of the influence of source confusion to real observations this must be done for artificial maps simulating an observation under realistic conditions, e.g., the map in Fig. 2.22.

For the sources detected in the artificial maps the corresponding sources are traced from the simulated objects. The flux densities and source positions derived from the Gaussian fits are compared to the same properties of the corresponding simulated galaxy.

In Fig. 2.24 the deviation of the detected flux density from the true one is plotted against the detected flux density. At about 0.6 mJy the flux density confusion error is mostly bigger than the detected flux density. In Fig. 2.25 the position error in units of the half width at half maximum of the observation beam is plotted against the detected flux density. At about 0.6 mJy the position errors of most detections are larger than one beam. My study of source confusion is preliminary, so that the plots may suffer from errors in the tracing of the simulated galaxy corresponding to a detection. This could explain the obvious gap in the position error plot.

Confusion limits the flux level down to which sources can be detected. Lowering the noise level through longer integrations does not lower the confusion error. Only a smaller beam would permit the detection of fainter sources.

Figure 2.24: The error in flux determination due to the source confusion. Figure 2.25: The error in flux determination due to the source confusion.

To determine the confusion limit estimates for real observations these studies have to be extended to artificial maps simulating observations under realistic conditions. Conducting such studies for a variety of simulations with different galaxy biasing will show how the galaxy clustering could affect the confusion limit.