text:
# Implications of Cosmological Gamma-Ray Absorption

## 1 Introduction

Understanding the evolution of large scale structure in the Universe is a major goal of modern observational cosmology. Numerical simulations of hierarchical structure formation in a globally homogeneous universe are now tractable (e.g., Nagamine, Cen, $`\&`$ Ostriker 2000; Kauffmann et al. 1999), but connecting the evolving structures to observable fluxes of electromagnetic radiation involves uncertain empirical descriptions of star formation, supernova feedback, and the dust-gas interplay. The necessary input comes from extensive observational campaigns, such as deep galaxy surveys, which measure the number of galaxies, their morphological types, colors, fluxes, and distances in presumably representative solid angles out to redshifts of z $`6`$. The wealth of detailed information derived from these observations can significantly complicate the effort to link theories of galaxy evolution and large scale structure formation. It is helpful to single out global quantities for which predictions can be compared with observations. One such quantity is the cosmic star formation rate (SFR) and its associated metagalactic radiation field (MRF) (the MRF at $`z=0`$ is commonly referred to as Extragalactic Background Light, EBL). The contribution of galaxies to the MRF is most significant between the far-infrared and the ultraviolet, while at longer wavelengths cosmic 2.7 K microwave background (CMB) radiation from the big bang dominates. At shorter wavelengths, accretion-powered active galactic nuclei provide much of the high-energy background (e.g., Mushotzky et al. 2000; Sreekumar 2000). In the GeV regime the MRF seems to originate from unresolved blazars and galaxies. Gamma rays emitted by novae, supernovae, and $`\gamma `$-ray bursts contribute the bulk of the observed background in the window around 1 MeV (e.g., Watanabe et al. 1999; Weidenspointner 1999, Ruiz-Lapuente et al. 2000). In principle, the evolution of the MRF should be predictable with structure formation models (e.g., Sommerville & Primack 1999), so that the observed MRF could be used to infer either the role of AGNs, low surface brightness objects, decays of relic particles, or to single out cosmological parameters. However, these models still rely on a wealth of uncertain parameters, and we are far from the ultimate goal of a first principles theory of the MRF. For an overview sie Hauser & Dwek (2001) and references therein.

Three basic methods are commonly employed for computations of the MRF from luminosity functions undergoing (i) forward evolution from a theoretically determined initial state, (ii) backward evolution from an observationally given final (present-day) state, or (iii) evolution that is directly observed over some range in redshift.

Method (i) starts from the theoretical framework of structure formation and evolution and predicts how luminosity functions evolve forward in cosmic time. The semi-analytic models of galaxy formation (e.g., White $`\&`$ Frenk 1991; Baugh et al. 1998; Kauffmann, White, $`\&`$ Guiderdoni 1993; Sommerville $`\&`$ Primack 1999, Granato et al. 2000) are based on structure formation studies with dissipationless N-body simulations. These studies yield luminosity functions for different morphological galaxy types that are in reasonable agreement with the observations. However, regarding the MRF, predictions from these models generally fail to satisfactorily reproduce observed cosmic emissivities. To better match these emissivities the models often require significant adjustments in the prescriptions of star formation, supernova feedback, and the inclusion of further astrophysical effects which presently can not be calculated from first principles. A simplified model for the infrared and sub-mm range was developed by Guiderdoni et al. (1998), while Malkan & Stecker (1998, 2001), en route of method (ii), determine the infrared MRF from local luminosity functions obtained with IRAS. In a very substantial paper, Franceschini et al. (2001) employ recent ISO data and more detailed models for the IR emission.

Method (iii) computes the MRF directly from the global SFR inferred from tracers of cosmic chemical evolution, such as the various Lyman $`\alpha `$ systems (Salamon & Stecker 1998, Pei, Fall & Hauser 1998), or from deep galaxy surveys (e.g., Madau et al. 1998, Rowan-Robinson 2001, Franceschini 2000). The spectral energy distribution (SED) for the globally averaged stellar population residing in galaxies can be estimated with population synthesis models (e.g. Bruzual & Charlot 1993) available for various input parameters, of which the initial mass function (IMF) and metallicity are the most important ones. Reprocessing by gas and dust can be taken into account explicitly via some model of the evolution of the dust and gas content in galaxies, in combination with assumed dust properties derived from local observations in the Milky Way.

Observational attempts to determine or constrain the present-day background face severe problems due to emissions from the Galaxy, which can introduce large systematic errors. Nevertheless, a number of studies with COBE FIRAS (Fixsen et al. 1998) and COBE DIRBE (Hauser et al. 1998) have resulted in highly significant detections of a residual diffuse IR background, providing an upper bound on the MRF in the IR regime. Similarly, the cumulative flux from galaxies detected in deep HST or ISO exposures provide lower limits to the present-day MRF. In the UV, measurements of the proximity effect provide an estimate of the MRF at high redshifts (e.g., Giallongo et al. 1996). To constrain cosmic evolution of the MRF one can also utilize the fact that high-energy gamma rays (from blazars or gamma-ray bursts) originating at large redshifts are attenuated by pair creation from interactions with low-energy MRF photons (e.g. Stanev& Franceschini 1998, Renault et al. 2001), which is the subject of a subsequent publication in this series.

Here we discuss a model of the evolving MRF that is based directly on observed emissivities (method iii), and is designed to use a minimal set of assumptions to clearly reveal the connections between input physics and output MRF. The method employed here (described in $`\mathrm{\S }`$2) is similar to the method discussed by Madau et al. (1998) or Malkan & Stecker (1998), but we specifically address redshift evolution of the MRF and the effects of dust-reemission in the infrared, the initial mass function (IMF) and metallicity. In $`\mathrm{\S }`$3 we discuss the use of population synthesis models to relate the SFR to the observed emissivities, and describe models of the dusty ISM in star forming regions that allow us to reproduce the far-infrared bump in the present-day MRF spectrum . Despite the complexity of the underlying physics involved in the production of the MRF, one can successfully model the MRF with simple modules. This approach allows us to investigate with clarity the various factors contributing to the MRF. In $`\mathrm{\S }`$4 we present the MRF spectrum as a function of redshift and discuss in detail the dependencies on cosmological models. We discuss the effect on the IR peak induced by varying assumptions about the IMF, the mean metallicity of the emitting stars and the effect by adding a new dusty population of galaxies, ULIGs/LIGs (ultraluminous/luminous infrared galaxies). Note that the MRF determined in this way does not depend on the parameters of the assumed cosmological model. However, when we refer to comoving emissivities or the cosmic star formation rate SFR(z), we adopt the flat Friedmann model with $`\mathrm{\Omega }_0=1`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0`$, and $`h=0.5`$ where $`h=H_0/(100\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1)`$. This choice of parameters is most commonly made in the observational literature, so we employ it here to allow direct comparisons.

## 2 Method

The method for calculating the MRF from a given SFR relies on an accurate knowledge of evolving stellar spectra and the reprocessing of star light in various dusty environments. Luminosity evolution of stellar populations is sensitive to the initial mass function (IMF), evolution of the mean cosmic metallicity, and the amount of interstellar extinction. Starting point of any model is the spectral energy density (SED) produced by a population of stars resulting from an instantaneous burst of star formation (commonly normalized to the mass of stars formed). Because star formation is an ongoing process with relatively short time scales of 105-7 yrs, the starburst spectra can be directly convolved with the global SFR, $`\dot{\rho _{}}(z)`$, to derive the evolution of the global luminosity density due to cosmic star formation. The SEDs are constructed from realistic stellar evolution tracks combined with detailed atmospheric models (e.g., Bruzual $`\&`$ Charlot 1993). The temporal evolution of the specific luminosity, $`L_\nu (t)`$ (in units of erg s-1Hz-1 per unit mass of stars formed) is then determined by the choices of IMF and the initial stellar metallicity. Fig. 1 shows the results for a Salpeter IMF between 0.1 and 100 M and solar metallicity. Note that the figure shows $`L_\lambda `$ as a function of wavelength. The luminosity drops rapidly as massive stars become supernovae (whose light is not included in these SEDs), and the wavelength of the bulk of the emission shifts to the red as the population ages. SEDs shown in Fig. 1 are unobscured by circumstellar gas and dust. The effects of absorption are discussed in $`\mathrm{\S }`$ 3.2.

From the population synthesis starburst models we obtain the comoving emissivity (or luminosity density) at cosmic epoch $`t`$ from the convolution

$$_\nu (t)=_{t_m}^tL_\nu (tt^{})\dot{\stackrel{~}{\rho }}_{}(t^{})𝑑t^{}(\mathrm{ergs}^1\mathrm{Hz}^1\mathrm{Mpc}^1)$$
(1)
where $`\dot{\stackrel{~}{\rho }}_{}(t)=\dot{\rho }_{}(z)`$ is the star formation rate per comoving unit volume. Rewriting Eq. (1) in terms of redshift, $`z=z(t)`$, yields

$$_\nu (z)=_z^{z_m}L_\nu (t(z)t(z^{}))\dot{\rho }_{}(z^{})\left|\frac{dt^{}}{dz^{}}\right|𝑑z^{},$$
(2)
where we assumed that star formation began at some finite epoch $`z_m=z(t_m)`$. For given evolution of the emissivity a second integration over redshift yields the energy density, or, after multiplication with $`c/4\pi `$, the comoving power spectrum of the MRF

$$P_\nu (z)=\nu I_\nu (z)=\nu \frac{c}{4\pi }_z^{z_m}_\nu ^{}(z^{})\left|\frac{dt^{}}{dz^{}}\right|𝑑z^{},$$
(3)
with $`\nu ^{}=\nu (1+z^{})/(1+z)`$. Cosmological parameters enter through $`dt/dz`$, which is given by (e.g., Peebles 1993 *Principles of Physical Cosmology*)

$$\left|\frac{dt}{dz}\right|=\frac{1}{H_0(1+z)E(z)}$$
(4)
with an “equation of state”

$$E(z)^2=\mathrm{\Omega }_\mathrm{r}(1+z)^4+\mathrm{\Omega }(1+z)^3+\mathrm{\Omega }_R(1+z)^2+\mathrm{\Omega }_\mathrm{\Lambda }.$$
(5)
The term proportional to $`\mathrm{\Omega }_\mathrm{r}`$ takes into account the contribution from relativistic components (such as the CMB and the star light). The density parameter of this component is defined as $`\mathrm{\Omega }_\mathrm{r}=u_\mathrm{r}/\rho _{\mathrm{crit}}c^2`$, where $`u_\mathrm{r}`$ refers to the relativistic energy density and $`\rho _{\mathrm{crit}}`$ is the critical density of the universe; $`\rho _{\mathrm{crit}}`$ = 3H$`{}_{}{}^{2}{}_{0}{}^{}`$/8$`\pi `$G = 10.54 h2 keV/cm3. The resulting dependencies on the Hubble constant are; emissivity $`H_0^1`$, and MRF power spectrum $`H_0^2`$. However, this scaling is only correct if we formally require that the star formation rate is a given function of redshift. Lacking a reliable derivation from first principles, this function is derived from observations that involve distance and luminosity estimates, which introduces additional powers of H0.

## 3 Emissivity

### 3.1 Spectral synthesis model

Conversion of gas to stars produces a stellar mass distribution that is commonly described by a “universal” initial mass function. A batch of stars produced in an instantaneous “burst” of star formation is often referred to as a Simple Stellar Population (SSP). Massive stars in the SSP have short lives ($``$ 107 yrs) and predominantly produce UV radiation, while long-lived, low-mass stars remain close to the main sequence even over cosmological times and produce the bulk of the “red light”. Depending on their mass, stars follow different evolutionary paths and evolve on different time scales, which causes the SED of a SSP to be a sensitive function of time. To follow the changing energy output in time we use the population synthesis code of Bruzual & Charlot 1999 (BC-Model), which is an updated version of the code documented in GISSEL96 (Leitherer et al. 1996). Figure 1 shows the resulting SEDs emitted at several distinct times after the burst. The stellar spectra used to construct these SEDs are based on Padova tracks (e.g. Girardi et al. 2000) and Lejeune stellar atmosphere models (Lejeune et al. 1997, 1998), and include a post-AGB evolutionary phase. For demonstration purposes we adopted a population with solar metallicity and Salpeter single power law IMF in the range 0.1$`M_{}<M<100M_{}`$, although different choices can readily be made with this code. The SEDs shown in Fig. (1) are of course not representative of any particular type of galaxy, because those involve star formation histories that are usually different from a single burst. However, the SEDs do resemble more closely the spectra of elliptical galaxies, for which the single star burst might be a reasonable approximation.

### 3.2 Interstellar Medium

In contrast to the intergalactic medium (IGM), which is extremely tenuous and has only a small effect on the transmitted SEDs from redshifts less than unity (e.g., Madau 2000), the intrinsic absorption by the galactic interstellar medium (ISM) is significant, and must be included to obtain correct SEDs. For simplicity, we assume a uniform distribution of gas and dust surrounding the stars of the SSP, and apply Osterbrook’s Case B recombination for optically thick clouds at an equilibrium temperature of $`10^4`$K, i.e. total absorption of all ionizing photons and reemission of 68% of the absorbed power in $`L_\alpha `$ line emission (which is subsequently absorbed by dust). The remaining energy is assumed to be reemitted in the optical regime via bremsstrahlung and recombination line/continuum emission. The assumption of homogeneity of the absorbing gas and dust layers is a simplification that could have a noticeable effect on the estimated transmission of UV radiation shortward of the Lyman edge. A more realistic approach should also take into account ionized superbubbles that are driven into the ISM by multiple supernovae (e.g., Dove, Shull, $`\&`$ Ferrara 2000).

The average metallicty of gas in galaxies slowly increases with cosmic time, but the present-day value is not known precisely (e.g., Pei, Fall, $`\&`$ Hauser 1999). We thus adopt an average extinction curve

$$A_\lambda =0.68E(BV)R(\lambda ^10.35)$$
(6)
with $`R=3.2`$ and where $`A_\lambda `$ with $`\lambda `$ \[$`\mu `$m\] determines the absorption coefficient according to $`g(\lambda )=10^{0.4A_\lambda }`$. Reemission is calculated as the sum of three modified Planck spectra

$$L_\lambda ^d(L_{bol})=\underset{i=1}{\overset{3}{}}c_i(L_{bol})Q_\lambda B_\lambda (T_i)$$
(7)
where

$`Q_\lambda \lambda ^1`$ and $`L_{bol}=L_{bol}(\tau )`$.

Two temperatures characterize warm and cold dust in galaxies. The third temperature is included to model a PAH component, which is assumed to emit like a Blackbody. Emission lines are not treated separately, because of the smoothing effects from integration over redshift. We normalize these three components realtive to each other (with the constants $`c_i,i=1\mathrm{}3`$) by using non-Seyfert galaxy observations by Spinoglio et al. (1995). We fitted the relation for all 4 bands.

Dust in the ISM of the Milky Way is known to coexist at several different temperatures, determined by the distances from various heat sources. Hot dust in spiral galaxies has temperatures ranging from 50 K to 150 K-200 K (Sauvage et al. 1997 and references therein) when in thermal equilibrium with HII regions, young massive stars, or compact accreting sources. Radiation from this dust component emerges in the mid-infrared and reprocesses only a small fraction of the emitted luminosity. Warm dust with temperatures between 25 K and 50 K corresponds to regions heated by the mean interstellar radiation field. Dust inside molecular clouds is somewhat shielded against high-energy radiation, and thus appears at low temperatures between 10 K and 25 K. Very cold dust at temperatures of 10 K, or even less, can be present in the densest parts of molecular clouds or in outer regions of the galaxy where the flux of the interstellar radiation field has dropped to the value of the intergalactic radiation field. Such very cold dust is difficult to detect, and requires sub-mm observations which so far have failed to provide unambiguous results. Therefore, we do not include very cold dust in our model. To keep the model simple, we consider variations in dust composition only because the shape of the spectrum is dependent of the total flux. It is noteworthy, however, that emission features around 10 $`\mu `$m due to Polycyclic Aromatic Hydrogen (PAH) molecules seem to be ubiquitous in galaxies (Desert el al. 1990). These PAH’s are undergoing temperature fluctuations and are generally not in thermal equilibrium. The broad emission lines of the PAHs are modeled with an additional (low flux) blackbody component, characterized by T $`425`$ K (Dwek et al. 1997).

As mentioned above, we use the non-Seyfert galaxy relations for all 4 bands by Spinoglio et al. (1995) to fix the 6 model parameters in Eq. 7 (three temperatures $`T_i`$, i=1…3 and three constants $`c_i`$, i=1…3). The relations depend on the total luminosity which is radiated by a galaxy. Although we are using a SSP and not a ”real” galaxy we can use the IR relations of galaxies, because these relations only depend on the total luminosity of a galaxy. We use $`10^{11}`$M for each SSP spectrum to get a galaxy-like total luminosity. Note that the dust-spectra $`L_\lambda ^d`$ depend on the age $`\tau `$ of the SSP, because they change with total luminosity $`L_{bol}`$.

As starting values for the temperature we use $`T_1`$=240K (consistent with a PAH blackbody of 425K), $`T_2`$=80K (warm dust) and $`T_3`$=30K (cold dust). After the fitting procedure we obtain for each SSP spectrum the associated IR-spectrum (see Fig. 2). The resulting temperatures are higher for higher total luminosities (for younger SSPs) and they are generally in the range 27K$`<T_1<`$33K, 125K$`<T_2<`$70K, and 180K$`<T_3<`$400K.

The resulting total spectra, including absorption and reemission, can be written for each SSP-spectrum as

$$L_\lambda (\tau )=[(L_{\lambda ,\mathrm{BC}}(\tau )+c_g(\tau )ϵ_\lambda )]g_\lambda +c_dL_\lambda ^d(\tau )$$
(8)
where $`L_{\lambda ,\mathrm{BC}}(\tau )`$ are the SEDs from the BC-Model with an age $`\tau `$. The optical emission of gas heated by absorbed photons is included through the quantity $`ϵ_\lambda \mathrm{exp}(h\nu /kT)`$ multiplied by the fractional energy $`c_g`$ available for this channel (32%). Both spectral components are modified by the absorption coefficient $`g_\lambda `$. The dust reemission spectra $`L_\lambda ^d(\tau )`$ are added according to Eq. (7). The IR-spectra shown in Fig. 2 are normalized with $`c_d`$ using energy conservation of the absorbed and re-emitted photons.

Published values for the color index, $`E(BV)`$, cover a wide range, which reflects large uncertainties associated with dust properties. Steidel et al. (1999) adopt $`E(BV)=0.15`$ at redshifts $`z=34`$. Madau et al. (1998) use a universal value of $`E(BV)=0.1`$, Guiderdoni (1999) prefers $`E(BV)=0.09`$ at $`z>2`$. Generally, dust extinction seems to play a more important role at high redshifts (e.g., Pettini et al. 1998). According to Madau et al (1998) the index varies as $`E(BV)=0.011(1+z)^{2.2}`$.

Using the relation $`\lambda L_\lambda (\tau )=\nu L_\nu (\tau )`$ we calculate $`L_\nu `$ to obtain the emissivity in Eq. 2

### 3.3 Star formation history

Recent deep galaxy surveys (e.g., Kennicutt 1983, 1999) or Lyman $`\alpha `$ absorber studies (Pei & Fall 1995; Pei, Fall, $`\&`$ Hauser 1999) suggest a functional shape of the SFR ($`\dot{\rho _{}}(z)`$) that can be approximated with a simple broken power law

$$\dot{\rho _{}}(z)(1+z)^\alpha $$
(9)
with $`\alpha =\alpha _m>0`$ for $`zz_{peak}`$ and $`\alpha =\beta _m<0`$ for $`z>z_{peak}`$. Plotted as a function of redshift the right hand panel of Fig. 3 shows the fit function in comparison to the more complex fit function given in Madau (1999). The cosmic star formation rate density SFR(z) has been determined with different methods and for large set of input data, as recently summarized by Ruiz-Lapuente et al. (2000). These studies suggest that the original Madau curve, Madau (1997(II)), should be considered a lower limit, and that realistic rates could be larger by a factor 2$``$3 at all redshifts. The compilation of Ruiz-Lapuente et al. (2000) (see their Fig. 1) clearly shows that we do not yet understand systematic effects well enough to obtain a reliable estimate for SFR(z). This is especially true at redshifts much beyond unity.

In our approach to modeling the MRF, the SFR function is considered to be a free fit function aimed at reproducing the emissivities derived from deep surveys. ”Measurements” of the SFR are generally based on luminosity densities, and thus model dependent. It is thus preferable to use the emissivities directly to obtain a self-consistent star formation history. We note, that the SFR parameterization used here does not contain any cosmological parameters. However, it is clear that choosing a different cosmology does change the observationally determined emissivities (Lilly et al. 1996, Ellis et al. 1996), hence requiring a different SFR.

For a given SFR, the emissivity is readily obtained from the convolution given by Eq. 2. The resulting $`_\nu (z)`$ is plotted in Fig. 3 for four different wavelengths in the optical band, and compared to the observations. Note that the steep increase at 0.28 $`\mu `$m and the shallower increase at 1.0 $`\mu `$m are reproduced by the model. The model slightly overproduces the present-day emissivity at 1.0 $`\mu `$m. However, the data point of Lilly et al. (1996) at 0.44 $`\mu `$m and $`z=0`$ falls much below the corresponding value obtained by Ellis et al. (1996), indicating the conservative nature of the Lilly et al. measurements. Generally, there is good agreement with the data.

### 3.4 Contribution from luminous infrared galaxies

Luminous Infrared Galaxies (LIGs) were first discovered with IRAS (Soifer et al. 1987). They represent a population of galaxies with IR luminosites above $`10^{11}`$L (L$`>10^{12}`$L are named ultraluminous infrared galaxies ULIGs) and high star formation rates. Most of these galaxies are dust enshrouded starburst galaxies or mergers, some of the ULIGs have also been identified as AGN (e.g. Kim et al. 1998). Although the LIGs are not so numerous today, a significant fraction of the infrared light could originate from them. Moreover, their star formation history could well be different from the SFR of optical galaxies (Franceschini et al. 2001). As these galaxies do hardly show up in optical surveys (E(B-V)$`{}_{\mathrm{LIG}}{}^{}>`$ E(B-V)OPT), there existence is not reflected in the SFR discussed in the previous section. Also, the dust temperatures in LIGs are higher, Hence we emphasize that adding another SFRLIG is a straightforward ramification of the model which would affect the infrared part of the MRF. In the following we will start to model the MRF without this contribution, but we will return to it later.

## 4 Results

### 4.1 Metagalactic radiation field

The final step in computing the MRF involves an integration of the emissivity over cosmic time using Eq. (3), where we neglected spectral modifications due to the intergalactic medium (IGM). The IGM consists predominantly of Lyman-$`\alpha `$ clouds with HI, HII and HeI gas (e.g., Madau 1997(I)) and mainly affects photons with wavelengths shortward of 911 Å which we assumed to be completely absorbed inside galaxies. Recent work by Pei, Fall & Hauser (1999) also shows that absorption in the IGM, if present, is minimal. The neglect of IGM effects can cause a slight underestimate of the star formation rate, but will not affect the estimate of the MRF spectrum, because our fit procedure uses the observed emissivities, which guarantees that all photons contributing to the MRF are accounted for. The evolution of the comoving MRF spectrum is shown in Fig. 4 for several redshifts, where the wavelength scale corresponds to a comoving reference frame.

### 4.2 Dependence on cosmology

Formally, the MRF intensity computed according to the procedures outlined in Sec. 2 exhibits an explicit cosmology dependence. This is the consequence of choosing an apparently cosmology-independent SFR. Such a choice is unrealistic, because if one were to change the cosmological model, while keeping the same SFR, one would immediately fail to reproduce the observed emissivities. However, the emissivities themselves are computed directly from observables, viz. the fluxes $``$ and redshifts of galaxies. We therefore have to take a look at the cosmology dependence of the emissivities implied by the relation

$$(z)=\frac{dL}{dV_c}=\frac{4\pi d_L^2(z)d}{dV_c}$$
(10)
where $`d_L(z)`$ is the luminosity distance and

$$dV_c=\frac{d_L^2d\mathrm{\Omega }dz}{H_0(1+z)E(z)}$$
(11)
the comoving volume element, and hence $`(z)H_0(1+z)E(z)`$. Inserting this in Eq. (3) strictly cancels the cosmology dependence. This argument means that by requiring the model to fit measured emissivities (evaluated for a given cosmology) the resulting MRF no longer explicitely depends on the choice of cosmological parameters. The background radiation field becomes, in a sense, a measured quantity itself, since both luminosity and volume scale $`d_l^2`$.

### 4.3 Bolometric flux

We complete this section with a comparison of the “bolometric” (IR - opt) flux obtained from our model with results available in the literature (Table 1). A strict lower limit on the present-day MRF flux of 28 nW m-2 sr-1 was derived from COBE and HST data (Dwek et al. 1998). The integrated flux from our model, 46 nW m-2 sr-1, is in agreement with all models using a similar SFR.

## 5 Discussion

We have developed a simple model and its evolution with time. This model is based on direct measurements of the global emissivities due to galaxies. There are no measurements of the MRF at high redshifts that could be used to directly verify our model, with the possible exception of the UV background based on the proximity effect (e.g., Giallongo et al. 1996).

In the second paper of this series we will introduce a method to constrain the MRF at high redshift using this model and high-energy observations of blazars. For now this leaves us with various measurements of the present-day MRF as the most relevant set of constraints. While the match is by no means perfect, the comparison shown in Fig. 5 suggests that our simple model is capable, without fine-tuning of parameters, of explaining the global level and the general spectral shape of the MRF. But the predicted flux in the IR-band falls short by roughly a factor two. As discussed above, the magnitude as well as the shape of the MRF are the result of a convolution of SEDs from an aging stellar population with continuous star formation at a given rate SFR(z) and the cosmological effects of redshift and time dilation. The fact that our model roughly matches both magnitude and shape of the MRF indicates that the relevant input physics has been properly taken into account. This provides confidence in the predictions of the MRF at higher redshifts, which is needed for a variety of astrophysical studies. However, there are still some deviations that need to be addressed, especially the shortfall in the IR band.

### 5.1 Population synthesis models

There are several model parameters that introduce significant uncertainties in the estimated MRF flux. Most important among them are the IMF, the dust extinction model, the parameterization of the dust emission by multiple blackbodies with different temperatures, the metallicity dependence of stellar evolution tracks, and the amplitude and form of the SFR (especially at redshifts less than unity). The choice of cosmology plays only a minor role for the MRF spectrum, as we emphasized in the previous section. We also investigated different population synthesis templates for the stellar light output (e.g., Leitherer et al. 1996 and references therein) and found them to be very similar to each other. The choice of the population synthesis model thus does not significantly affect the estimated MRF.

### 5.2 IMF

We proceed with a discussion of those parameters that, at least in principle, can significantly alter the estimated intensity, flux, or energy density. We start by comparing results based on the Salpeter and Scalo IMF (see Fig. 5. One of the distinguishing features of the Scalo IMF is the fact that it contains relatively few high-mass stars. These stars are responsible for most of the UV photons, which, after thermalization by the dust in the ISM, emerge in the FIR. On the other hand, the Scalo IMF contains a relatively large fraction of low mass stars which emit most of their light at optical wavelengths. Consequently the present-day MRF has a much weaker FIR bump when the Scalo IMF is employed. This perhaps provides an argument in favor of the Salpeter form for a global IMF, because the Scalo IMF would somewhat underproduce the FIR background.

### 5.3 Star formation rate

The star formation rate density obtained here (see Fig. 3) is higher than the SFR originally suggested by Madau (1997(II), 1999), and is somewhat different from the rate derived from structure-formation theory (Primack et al. 1998). As discussed above, recent determinations of the cosmic star formation history based on Hα emission and ISO data suggest that the Madau rate has been systematically underestimated by a factor 2$``$3 (see Flores et al. 1999 for a recent discussion). The measurements of the sub-mm SCUBA array (Hughes et al. (1998) support the notion that much of the star formation activity at high redshifts is hidden due to dust absorption. Ruiz-Lapuente, Casse, $`\&`$ Vangioni-Flam (2000) summarize many of these measurements and compare (their Fig. 1) the various functions with that of Madau et al. (1998). The function used in this study falls above the ”Madau curve”, but below most of the curves compiled by Ruiz-Lapuente et al. (2000).

### 5.4 Absorption and re-emission of the ISM

The dust and gas model we use is not based on first principles, but is founded on empirical results. A three temperature model for galactic dust spectra has been proposed by Dwek et al. (1998). A small change in the temperatures (say, $`\pm 10`$K) would only cause a small change in the shape of the MRF spectrum. This lack of sensitivity originates from the broadening of the employed modified blackbody spectra due to the integration over redshift.

Any modification of $`E(BV)`$ changes the spectral shape of the MRF from 911 Å to $`10^4`$ Å, especially the amplitude of the far-infrared bump. This is simply due to energy conservation. An increase in the extinction causes a larger fraction of UV absorption, and this energy re-emerges predominantly in the FIR. We selected a value of $`E(BV)`$ that provided an acceptable fit to the available data on emissivities as a function of redshift and the present-day background. We found that $`E(BV)=0.14`$ is the appropriate value for young SSP’s and 0.03 for old SSP. These values are reasonably well determined by the emissivity-fit alone. In any case, changing the $`E(BV)`$ value does not provide a solution to the ”missing IR flux”.

### 5.5 Comparison with data and other models at $`z=0`$

The model present-day MRF flux at optical wavelengths is consistent with lower limits from HST (Pozzetti et al. 1998, 2000). Observations by Bernstein et al. (2000) suggest the possibility of a somewhat higher MRF flux, but our results still fall within their estimated uncertainties.

Absolute measurements (albeit with large systematic errors due to the need for subtraction of dominant local foreground emission Lagache et al. 1999, Puget & Lagache 2000) are available from COBE/DIRBE and FIRAS in the IR band. A discrepancy between our model and the observations occurs at $`\mu `$m wavelengths, where the measured flux appears to be larger by a factor of two. A lower limit based on galaxy counts due to Elbaz et al (1999) at $`10^5`$Å lies in the range of possible PAH emission.

While the implications for the MFR at high redshifts are rarely stated in the literature, considerable effort has gone into the calculation of the present-day MFR, i.e., the extragalactic background light (EBL). In spite of great diversity in the computational approaches employed (see the discussion in the introduction), the models generally show the same IR deficit noticed in this study (e.g. Dwek et al 1998, Primack 1998). Dwek et al. suggested a possible solution by adding a new, distinct component of obscured galaxies which emit preferentially in the IR band.

### 5.6 Metallicity

The calculations for our ”standard model” employ stellar models with solar metallicity. Salamon & Stecker (1998) and Pei, Fall & Hauser (1999) discuss some of the possible effects caused by changes in the mean cosmic metallicity. The approximation of a fixed, high (solar) metallicity is motivated in part by structure formation simulations (e.g., Valageas $`\&`$ Silk 1999) which suggest that the mean metallicity in star forming regions is a slowly rising function for redshifts less than 2. The calculations of Valageas $`\&`$ Silk also indicate that the metallicity has been larger than 1/2 Z since z $``$ 2 (note that their estimate only includes enrichment due to SNII). The significantly lower metallicity values found in damped Lyman Alpha systems (DLAs; e.g., Pettini et al. 1997) probably correspond to the enrichment history of galactic halos. The trend of \[Zn/H\] with redshift suggests a present-day metallicity of 1/3 solar (Vladilo et al. 2000), which also indicates that DLAs do not trace the chemical evolution of proto-disks, but instead star-forming fragments which build up galaxies through mergers. Although the Milky Way does not represent a good template for cosmic chemical evolution (e.g., Prantzos $`\&`$ Silk 1989; Fields 1999) its age-metallicity relation and well known radial metallicity gradients suggest that much of its current star forming activity takes place in environments with Z $``$ Z, or even higher (see Boissier $`\&`$ Prantzos 1999 for a recent model of galactic chemical evolution). We thus expect the Z = Z = const assumption to provide a reasonable approximation. However, to test the dependence of our results on metallicity, we performed the MRF simulation for two cases: i) solar metallicity Z = 0.02, and ii) Z = 0.0001. The dashed lines in Fig. 6 show the “low-Z MRF” in comparison to the standard model with solar metallicity (solid line). Both simulations used the Salpeter IMF. The comparison indicates that metallicity effects could be important. Lower metallicity in the stellar atmospheres leads to a higher fraction of light primarily emitted in the UV, and subsequently redistributed towards the IR by interstellar dust grains, with extinction parameters determined newly from the fit to the emissivities (i.e. independent of the assumed low metallicity in the stellar atmospheres producing the bulk of the light). A rigorous treamtent of metallicity effects, in order to obtain self-consistent interstellar extinction curves and chemical evolution (Pei, Fall, & Hauser 1999), is beyond the scope of the paper.

### 5.7 ULIGs/LIGs

The approach in our model up to this point was to consider the emission from galaxies found in deep optical surveys, and to compute a mean galaxy spectrum for them. This is certainly a valid scheme to obtain a lower limit to the MRF, and, in fact, the model does not overproduce observational upper limits at any frequency. However, the infrared deficit revealed by the preceding analysis, seems to justify the inclusion of a population of luminous infrared galaxies (see Sect. 3.4). This would increase the strength of the infrared bump in the MRF, but leave the optical bump unchanged. Hence the overall performance of the model to reproduce the present-day extragalactic background spectrum can be improved. The star formation history used for this population SFRLIG (Dwek et al. 1998) has been infered from number counts in ISO, IRAS and SCUBA data, and the luminosity functions from Chary & Elbaz (2001) (and references therein). The extinction parameter is set to be E(B-V)=1.0, a value which is typical for LIGs. The resulting SFRLIG ($`\alpha `$ = 5.5, $`\beta `$ = 0, zpeak=1.0, $`\dot{\rho _{}}(z)=0.10\mathrm{M}_{}\mathrm{yr}^1\mathrm{Mpc}^3`$) is in agreement with the one calculated by Chary & Elbaz. At low redshifts the contribution of stars in LIGs is very small, but at high redshifts the number of stars formed in LIGs is comparable to the one in optically selected galaxies. This is in line with the interpretation of the IRAS (Kim & Sanders 1998) observations, which indicate that LIGs were more numerous in the past than they are today. Fig. 6 shows the EBL and the respective contributions from optical and infrared galaxies.

## 6 Conclusions

We have developed a model for the evolving MRF based on optical galaxy surveys as its main observational input, and found that this model shows a deficit at infrared wavelengths in the spectrum of the EBL (the MRF at $`z=0`$). Inclusion of obscured, infrared-emitting galaxies provides a viable solution of the problem, and we have determined their SFR from fitting the model to EBL data. The model in this form can serve as a reliable basis for obtaining predictions of the MRF at high redshifts. Observations of high-redshift gamma ray sources with next-generation gamma ray telescopes (GLAST, HESS, MAGIC, VERITAS) are expected to soon provide evidence for gamma ray attenuation due to collisions of gamma rays with low-energy photons from the MRF, thus allowing to test the model predictions in an independent way (paper II in this series). 666 At http://wwww.astro.uni-wuerzburg.de/theorie/ we provide a web tool for online calculations of the MRF, allowing users to select values of key parameters..

###### Acknowledgements.

We thank Stephane Charlot for providing the latest release of the population synthesis templates (GISSEL99), and Joel Primack for communicating the SFR function from structure formation theory prior to publication. We also thank Floyd Stecker for fruitful discussions. This work is part of a research project supported by the BMBF under code DESY-HS/AM09731M45. DH acknowledges support and hospitality during visits to the Universitäts-Sternwarte Göttingen.


format: text

text:
# Study of the 𝜌, 𝜔, ϕ→𝜂𝛾→7𝛾 Decays with an SND Detector on a VEPP-2M Collider

## Abstract

The $`e^+e^{}\eta \gamma 7\gamma `$ process was studied in the energy range $`2E=600÷1060`$ MeV with an SND detector on a VEPP-2M $`e^+e^{}`$ collider. The decay branching ratios $`B(\varphi \eta \gamma )=(1.343\pm 0.012\pm 0.055)10^2`$, $`B(\omega \eta \gamma )=(4.60\pm 0.72\pm 0.19)10^4`$, and $`B(\rho \eta \gamma )=(2.69\pm 0.32\pm 0.16)10^4`$ were measured.

Radiative decays of light vector mesons — $`\rho `$, $`\omega `$, and $`\varphi `$ — are important for understanding the behavior of the strong interaction at low energies. Although many measurements were carried out for the probabilities of the radiative decays, the achieved accuracy  is insufficient for reliable determination of the parameters of phenomenological models \[2-4\].

We report the results of investigations of the $`e^+e^{}\eta \gamma `$ process followed by the $`\eta 3\pi ^06\gamma `$ decay. Since the final state includes seven photons, the background can be substantially suppressed compared to that in other channels of $`\eta `$-meson decay and, there-fore, the systematic error may be reduced.

The experiment  was carried out in 1998 at the VEPP-2M $`e^+e^{}`$ collider with the SND detector . Two scans were performed in the energy range $`2E_0=984÷1060`$ MeV (PHI-98 experiment) with an integral luminosity of $`8.0pb^1`$ at 16 energy points and with about $`10^7`$ produced $`\varphi `$-mesons. In addition, a scan (OME-98 experiment) over 38 points in the energy range $`2E_0=360÷970`$ MeV was carried out with an integral luminosity of $`3.5pb^1`$ and with about $`310^6`$ produced $`\rho `$ and $`\omega `$ mesons.

The events of the process

$$e^+e^{}\eta \gamma ,\eta 3\pi ^0,\pi ^02\gamma $$
(1)
are characterized by the final state with seven photons, a few of which may not be detected. The extra photons may also appear due to the splitting of a shower in the calorimeter, the emission of the photons by the initial particles at large angles, or the superposition of the beam background. The main background process in the $`\varphi `$-resonance region is the $`\varphi K_SK_L`$ decay, where $`K_S`$ decays into two neutral pions and $`K_L`$ , interacting in the calorimeter, produces extra “photons.” An additional background is formed by the $`e^+e^{}\omega \pi ^0+X`$ process followed by the $`\omega \pi ^0\gamma `$ decay, where $`X`$ are extra photons. An analysis of the experimental data has demonstrated that the QED process $`e^+e^{}3\gamma `$, being superposed with other events, may also result in the required event configuration.

Taking the above-listed background events into account, we selected events in two steps. At the first step, among the events in which six or more photons and no charged particles were detected, we selected those satisfying the following conditions imposed on a total energy release $`E_{tot}`$ in the calorimeter and the total momentum $`P_{tot}`$ of photons:

$$E_{tot}/2E_0<1.2,P_{tot}/2E_0<0.2/c,E_{tot}/2E_0cP_{tot}/2E_0>0.7.$$

For the selected events, we performed a kinematic reconstruction using the measured angles, the energies of the photons, and energy-momentum conservation. As a result, the energies of the photons were determined more accurately and the $`\chi ^2`$ values specifying the degree of certainty of a process were determined:

* $`\chi ^2`$ for the assumption of the $`e^+e^{}n\gamma `$ process with $`n6`$ or 7 ;
* $`\chi _{3\gamma }^2`$ for the assumption of the $`e^+e^{}2(3)\gamma +X`$ process ;
* $`\chi _{\omega \pi ^0}^2`$ for the assumption of the $`e^+e^{}\omega \pi ^0+X`$ process .

Further selection was carried out with the restrictions

$$\chi ^2<30,\chi _{3\gamma }^2>20,\chi _{\omega \pi ^0}^2>20.$$

Figure 1 shows the distributions of the selected events in the recoil mass

$`M_{rec\gamma }`$ of the highest energy photon. It is seen that the desired process prevails in all scans. The events with $`M_{rec\gamma }>600`$ MeV in Fig. 1b are determined by the process

$$e^+e^{}\varphi K_SK_L.$$
(2)
Finally, we select the events satisfying the condition $`400<M_{rec\gamma }<600`$ MeV.

The number $`N(s)`$ of the observed events at a given energy is described by the formula

$$N(s)=L(s)\left[ϵ(s)\beta (s)\sigma (s)+\sigma _b(s)\right],s=4E_0^2,$$
(3)
where $`L(s)`$ is the integral luminosity, $`ϵ`$ is the detection efficiency determined by a simulation, $`\beta `$ is the factor representing the radiative corrections, $`\sigma _b`$ is the cross section for background processes, and $`\sigma `$ is the cross section for the desired process (1).

When determining the background from process (2), inaccuracy in the simulation of the interaction of the $`K_L`$ meson with a substance in the calorimeter is possible. Figure 1b demonstrates that the contribution of process (2) dominates for $`M_{rec\gamma }>600`$ MeV, while the contribution of the desired process (1) is negligible. For this reason, the number of events of process (2) in the range $`400<M_{rec\gamma }<600`$ MeV was determined from the number of experimental events in the interval $`600<M_{rec\gamma }<800`$ MeV by taking into account the simulated ratio of the numbers of the $`K_SK_L`$ events that fall into the mass ranges $`400<M_{rec\gamma }<600`$ MeV and $`600<M_{rec\gamma }<800`$ MeV.

The energy dependence of the resulting cross section (Fig. 2) was parametrized

by the vector-dominance formulas  including the contributions of the $`\rho `$, $`\omega `$, and $`\varphi `$ resonances:

$$\sigma (s)=\frac{F(s)}{s^{3/2}}\left|\underset{V=\rho ,\omega ,\varphi }{}\sqrt{12\pi \underset{B_{VP}}{}\frac{m_V^3}{F(m_V^2)}}\frac{\mathrm{\Gamma }_Ve^{i\phi _V}}{D_V(s)}\right|^2,$$
(4)
where $`F(s)=[(sm_\eta ^2)/2\sqrt{s}]^3`$, $`D(s)=m_V^2si\sqrt{s}\mathrm{\Gamma }_V(s)`$, and the products

$$\underset{B_{VP}}{}=Br_{Ve^+e^{}}Br_{V\eta \gamma }Br_{\eta 3\pi ^0}Br_{\pi ^02\gamma }^3$$

are the free parameters of the approximations. The relative phase shifts of the resonances were fixed at $`\phi _\rho =\phi _\omega =0`$, $`\phi _\varphi =\pi `$. The approximation gives the following results:

$`{\displaystyle \underset{B_{\varphi P}}{}}`$ $`=`$ $`(1.249\pm 0.011\pm 0.035)10^6,`$
$`{\displaystyle \underset{B_{\omega P}}{}}`$ $`=`$ $`(1.01\pm 0.16\pm 0.03)10^8,`$
$`{\displaystyle \underset{B_{\rho P}}{}}`$ $`=`$ $`(3.77\pm 0.45\pm 0.11)10^9.`$ (5)
where the first error is statistical and the second, systematic, error is due to the contributions of errors in the determination of the detection efficiency and the error in the measurement of the luminosity. The luminosity was measured from elastic electron-positron scattering at large angles and from the process of two-quantum annihilation. The difference in the results of the two methods does not exceed 1%. The accuracy of the theoretical formulas used for the simulation of elastic scattering and experimental conditions provides an estimate of about 2% for the accuracy of the luminosity in the case under consideration. In order to estimate the systematic errors in the detection efficiency, the stability of the results to a change in the selection conditions was examined: we added restrictions on the polar angle of the photons and on the number ($`N_\gamma =7`$) of particles and used only the completely reconstructed $`e^+e^{}\eta \gamma 7\gamma `$ events. As was discussed above, due to the emission of the initial particles at large angles and superposition of preceding events, extra spurious photons appear in the SND calorimeter. For this reason, one of the tests of the kinematic reconstruction was carried out with the exclusion of photons with energies lower than 50 MeV and with a polar angle less than $`36^{}`$. In addition, two scans of the $`\varphi `$ meson were independently processed. All tests demonstrated the stability of the results, and the total systematic error in the efficiency, with the inclusion of all effects, was estimated at 2%. This estimate is treated as independent of the systematic error in luminosity.

Using the tabular values of $`Br_{Ve^+e^{}}`$, $`Br_{\eta 3\pi ^0}`$, and $`Br_{\pi ^02\gamma }`$ from , we obtain from Eqs. (Study of the $`\rho `$, $`\omega `$, $`\varphi \eta \gamma 7\gamma `$ Decays with an SND Detector on a VEPP-2M Collider) the values

$`Br_{\varphi e^+e^{}}Br_{\varphi \eta \gamma }`$ $`=`$ $`(4.017\pm 0.035\pm 0.124)10^6,`$
$`Br_{\omega e^+e^{}}Br_{\omega \eta \gamma }`$ $`=`$ $`(3.25\pm 0.51\pm 0.10)10^8,`$
$`Br_{\rho e^+e^{}}Br_{\rho \eta \gamma }`$ $`=`$ $`(1.21\pm 0.14\pm 0.04)10^8,`$
$`Br_{\varphi \eta \gamma }`$ $`=`$ $`(1.343\pm 0.012\pm 0.055)10^2,`$
$`Br_{\omega \eta \gamma }`$ $`=`$ $`(4.60\pm 0.72\pm 0.19)10^4,`$
$`Br_{\rho \eta \gamma }`$ $`=`$ $`(2.69\pm 0.32\pm 0.16)10^4,`$
$`\sigma _{\varphi \eta \gamma }`$ $`=`$ $`(56.75\pm 0.50\pm 1.68)nb,`$
$`\sigma _{\omega \eta \gamma }`$ $`=`$ $`(0.78\pm 0.12\pm 0.02)nb,`$
$`\sigma _{\rho \eta \gamma }`$ $`=`$ $`(0.300\pm 0.036\pm 0.009)nb,`$ (6)
where $`\sigma _{V\eta \gamma }=12\pi Br_{Ve^+e^{}}Br_{V\eta \gamma }/m_V^2`$ and the errors of the tabular values are included in the systematic errors. Using data (Study of the $`\rho `$, $`\omega `$, $`\varphi \eta \gamma 7\gamma `$ Decays with an SND Detector on a VEPP-2M Collider) obtained above, the results of previous measurements with the SND detector , and the data on the widths of the $`\rho `$ and $`\omega `$ mesons , we derive the ratios

$`{\displaystyle \frac{B(\eta \pi ^+\pi ^{}\pi ^0)}{B(\eta \gamma \gamma )+B(\eta 3\pi ^0)}}`$ $`=`$ $`0.304\pm 0.012,`$
$`{\displaystyle \frac{B(\eta 3\pi ^0)}{B(\eta \gamma \gamma )}}`$ $`=`$ $`0.826\pm 0.024,`$
$`{\displaystyle \frac{\mathrm{\Gamma }_{\rho \eta \gamma }}{\mathrm{\Gamma }_{\omega \eta \gamma }}}`$ $`=`$ $`10.5\pm 2.1.`$ (7)
Our results (5)-(7) are in agreement with the data of other experiments \[8-12\]. The branching ratios for the $`\varphi `$, $`\omega \eta \gamma `$ decays are measured with an accuracy close to the tabular one , and the branching ratio for the $`\rho \eta \gamma `$ decay is determined with the doubly improved accuracy. Note that the quantity $`Br_{\varphi e^+e^{}}Br_{\varphi \eta \gamma }Br_{\eta 3\pi ^0}`$ (Study of the $`\rho `$, $`\omega `$, $`\varphi \eta \gamma 7\gamma `$ Decays with an SND Detector on a VEPP-2M Collider) was measured with noticeably higher accuracy than the branching ratio $`Br_{\varphi \eta \gamma }`$ (Study of the $`\rho `$, $`\omega `$, $`\varphi \eta \gamma 7\gamma `$ Decays with an SND Detector on a VEPP-2M Collider), because the leptonic width of the $`\varphi `$ meson is known with an accuracy of 2.7%, which is noticeably worse than the statistical accuracy of our measurement.

This work was supported by the “Russian Universities” Foundation (project no. 3N-339-00) and the Russian Foundation for Basic Research (project nos. 99-02-16813, 00-02-17478, and 00-02-17481).


text:
# 1. Background

## 1. Background

Required reading: Symbol-Crunching with the Transfer-Matrix Method in Order to Count Skinny Physical Creatures \[Z1\].

In \[Z1\], Zeilberger used finite transfer matrix method and developed two Maple packages ANIMALS and FreeANIMALS to count square lattice animals confined to a strip. Here we define a bijection between hexagonal lattice animals and a special class of square lattice animals and adopt the packages HexANIMALS and HexaFreeANIMALS accompanying \[Z1\] to enumerate hexagonal lattice animals.

Definition 1.1 A hexagonal lattice animal (hexagonal polyominoes hence forth) on the hexagonal lattice is an edge-connected set of lattice cells on the hexagonal lattice. Two animals are equivalent if they are translations of each other.

Figure 1 shows the three hexagonal polyominoes comprised of two hexagons.

The number of nonequivalent hexagonal polyominoes with $`n`$ cells, $`a(n)`$, is given by Sloane’s sequence A001207 for $`n35`$ \[NJAS\]. For example, $`a(1)=1`$, $`a(2)=3`$, $`a(3)=11`$, $`a(4)=44`$, $`a(5)=186`$, $`a(6)=814`$, $`a(7)=3652`$, $`a(8)=16689`$, $`a(9)=77,359`$, and $`a(10)=362,671`$.

A hexagonal polyominoes in which each column contains at most $`k`$ contiguous blocks of cells is referred to as a $`k`$-board polyominoes. Figure 2 shows examples of one and two-board polyominoes with $`7`$ hexagons. The number of nonequivalent one-board polyominoes with $`n`$ hexagonal cells, $`b(n)`$, is given by Sloane’s sequence A059716 for $`n24`$ \[GV\]. For example, $`b(1)=1`$, $`b(2)=3`$, $`b(3)=11`$, $`b(4)=42`$, $`b(5)=162`$, $`b(6)=626`$, $`b(7)=2419`$, $`b(8)=9346`$, $`b(9)=36,106`$, and $`b(10)=139,483`$.

The generating function that enumerates one-board hexagonal polyominoes is computed in \[K\], which can also be computed using the Maple package LEGO in \[Z0\] by taking $`p(a,b)=a+b`$ (for definition of $`p(a,b)`$ refer to \[Z0\]). To the best of the authors’ knowledge, $`k`$-board hexagonal polyominoes have not been enumerated for $`k2`$. In this article, in addition to computing generating functions (and series expansions) for hexagonal lattice animals that fit into a prescribed but arbitrary height, we also compute the first $`12`$ terms of the sequence that enumerate the number of board-pair-hexagonal polyominoes, the analog of square board-pair-polyominoes computed in \[Z2\]. Fisrt we recall some definitions and introduce the essential details of the Maple implementation of the transfer-matrix method \[S\] from \[Z1\].

## 2. The Transfer-Matrix Method

Let $`G(V,E)`$ be a vertex weighted directed graph consisting a finite set of vertices, $`V`$, and a finite set of directed edges, $`E`$. A path $`P`$ in G is a sequence $`v_1,e_1,v_2,e_2,\mathrm{},v_{m1},e_m,v_m`$ where $`v_iV`$($`1im`$) and $`e_iE`$($`1im`$) is an edge from $`v_{i1}`$ to $`v_i`$. The weight of a path $`P`$ is the sum of the weights of all vertices participating in the path $`P`$. In this case $`Wt(P)=Wt(v_1)+Wt(v_2)+\mathrm{}+Wt(v_m)`$.

For $`T`$ a subset of $`V`$, we want to compute the weight enumerator (generating function),

$$F(z)=\underset{j=0}{\overset{\mathrm{}}{}}a_jz^j$$

where $`a_j`$ is the number of paths with weight $`j`$ that starts at any vertex in $`V`$ and ends at a vertex in $`T`$. For $`vV`$, let $`F_v(z)`$ be the generating function that enumerates all paths starting at vertex $`v`$ and ends at a vertex in $`V`$. Then, clearly,

$$F_v(z)=\underset{\begin{array}{c}paths\text{P}\\ init(\text{P})=v,fin(\text{P})T\end{array}}{}t^{wt(\text{P})}$$
(1)
where $`init(P)`$ is the initial vertex of $`P`$ and $`fin(P)`$ is the terminal vertex of $`P`$. Let $`N(v)`$ be the set of vertices adjacent to $`v`$. Then, $`\left\{F_v(z)\right|vV\}`$ satisfies the following $`\left|V\right|`$ (where $`\left|V\right|`$ is the number of vertices) system of equations in $`\left|V\right|`$ unknowns, namely $`\left\{F_v(z)\right|vV\}`$,

$$F_v(z)=\mathrm{𝟏}_T(v)z^{wt(v)}+z^{wt(v)}\underset{uN(v)}{}n(v,u)F_u(z)$$
(2)
where $`n(v,u)`$ is the number of paths from $`v`$ to $`u`$ and

$$\mathrm{𝟏}_𝐓(𝐳)=\{\begin{array}{cc}\mathrm{𝟏}& \mathrm{𝐢𝐟}𝐳𝐓\\ \mathrm{𝟎}& \mathrm{𝐨𝐭𝐡𝐞𝐫𝐰𝐢𝐬𝐞}\end{array}$$
(3)
Since each $`F_u(z)`$ exists as a power series in $`z`$, our system has a solution. Once we solve for this sytem, the answer to our original problem is then

$$F(z)=\underset{vT}{}F_v(z)$$

As an example consider the following vertex-weighted directed graph on four vertices with weights as shown next to the vertices.

In this example, the generating function that enumerates all paths according to their weight is found by first finding the four generating functions $`\{F_u(z),F_v(z),F_w(z),F_t(z)\}`$ for all paths starting at vertex $`u,v,w,`$ and $`t`$ respectively. These four generating functions are related by

$$f_t(z)=z,f_w(z)=f_v(z)=zf_u(z)=\frac{z(1+z+z^2+z^3+z^4)}{1z^5}.$$

Solving these system and adding together we get the required generating function that enumerates all paths according to their weight as

$$F(z)=\frac{2z+2z^2+z^3+3z^4+3z^5+z^6}{1z^5}.$$

For the hexagonal polyominoes case, we follow the structure of \[Z1\] and automate the problem and use Maple do the hard part. First, we recall the following definitions from \[Z1\]:

Definition 2.1 A Combinatorial Markov Process is a six tuple (V, E, init, fin, Start, Finish), where $`V`$ is a finite set of vertices, $`E`$ is the set of directed edges, and $`init,fin:EV`$ are functions that assign the initial and terminal vertex to an edge, respectively.

Definition 2.2 A vertex Weighted Combinatorial Markov Process is a seven tuple (V, E, init, fin, Start, Finish, wt), where (V, E, init, fin, Start, Finish) is a Combinatorial Markov Process and $`wt:V(E)^+`$ is a function that assigns weights to vertices.

In Maple, we represent a Combinatorial Markov Process with vertices $`V=\{1,2,\mathrm{},n\}`$ and no multiple edges as a four tuple \[S,T,ListOfOutGoingNeighbors,ListOfWeights\] where S,T $`V`$, and ListOfOutGoingNeighbors and ListOfWeights are lists of length $`n`$ whose $`i^{th}`$ element is the set of vertices adjacent to vertex $`i`$ and the weight of vertex $`i`$, respectively.

Similarly, the Combinatorial Markov Process for a multiple edge graph is a four tuple \[S,T,ListOfOutGoingNeibors,ListOfWeights\] where S, T, and ListOfWeights are as previously defined and ListOfOutGoingNeighbors is a list of length $`n`$ whose $`i^{th}`$ element is a multi-set $`j_1^{m_1}j_2^{m_2}\mathrm{}j_k^{m_k}`$ represented in the form $`\{[j_1,m_1],[j_2,m_2],\mathrm{},[j_k,m_k]\}`$, which means that from vertex $`i`$ there are $`m_1`$ edges to vertex $`j_1`$, $`m_2`$ edges to vertex $`j_2`$, and so on.

If we successfully model hexagonal polyominoes as a weighted Combinatorial Markov Process as it is done for the squre lattice animals in \[Z1\] then we can employ the Maple package MARKOV given in \[Z1\] to automatically find the generating function and series expansion of the generating function. The next section describes how weighted Combinatorial Process for hexagonal polyominoes is created.

## 3. Maple Representation of hexagonal polyominoes

We could develop a grammar that describes height-restricted hexanimals from first principles and apply the same methods as Zeilberger \[Z1\] to encode the grammar as a Combinatorial Markov Process; however, this would result in rewriting much of the code that was already created for the square lattice animals. Our approach here is to define a bijection between hexagonal lattice animals and a special class of square lattice animals and adopt the Maple codes ANIMALS and FreeANIMALS in \[Z1\]. In the remaining sections, we describe this processes.

We overlay a hexagonal polyominoe on the square lattice so that its leftmost vertex (vertices of the square contained inside the hexagonal) lie on the line $`x=0`$ and its bottom most cell lie on the line $`y=0`$ as shown in Fig.4(a). We then construct a polyominoes comprised of the square cells which have diagonally opposing corners falling on the same hexagonal cell as shown in Fig. 4(b).

As described in detail by Zeilberger \[Z1\], we represent a polyomino by the set of coordinates of the bottomleft corner of each of its cells. We further encode a polyomino as a word in the alphabet consisting of the non-empty subsets of $`\{0,1,\mathrm{},k\}`$ where $`k`$ is the maximum $`y`$ coordinate in the set representation of the polyomino. For example, the animal in Fig. 4 (b) is represented by the set

$$\{(0,2),(0,3),(1,1),(1,2),(1,5),(1,6),(2,0),(2,1),(2,2),(2,3),(2,4),(2,5),(3,1),(3,2)\}$$

and encoded as the word (in interval notation)

$$\{[2,3]\},\{[1,2],[5,6]\},\{[0,5]\},\{[1,2]\}.$$

There are two important properties of the polyomino words that encode hexanimals. First, every interval has even size (i.e., it represents a contiguous block of an even number of cells); this is obvious from the mapping which associates each hexagonal cell with a pair of square cells. Second, within a letter, the starting position of every interval has the same parity (odd or even), and adjacent letters have different parity; this is again obvious from the mapping of each hexagonal cell to a pair of square cells and the vertical shift of adjacent hexagonal columns. Note, however, that the parity of the leftmost column may be either odd or even depending on the hexanimal.

We will call the class of square polyominoes to which hexagonal polyominoes are mapped parity polyominoes. It is clear that the mapping between hexagonal animals and parity polyominoes is unique and reversible, and thus defines a bijection. As a result of this bijection, we can count hexagonal animals by counting parity polyominoes. Zeilberger’s code counts general polyominoes (a much larger class than parity polyominoes), so in order to use it to count parity polyominoes, we need to restrict the alphabet. In the following sections of this paper, we describe how to modify Zeilberger’s code to count parity polyominoes while minimizing the amount of changes required; we assume that the reader is familiar with \[Z1\].

## 4. Counting Globally Skinny Hexagonal Animals

Zeilberger’s ANIMALS package uses the transfer-matrix method \[S\] to count polyominoes whose set representation has $`y`$ coordinates in the range $`0yk`$ for some arbitrary $`k`$. Before the transfer-matrix method can be applied, a grammar needs to be constructed to describe how to combine letters into words that encode polyominoes. A key routine in the generation of this grammar is PreLeftLetters(a,b) which returns the set of all possible subsets of the integers in the range $`[a,b]`$ written in interval notation. For example, PreLeftLetters(0,2) returns the following subsets:

$$\{\left\{\right\},\left\{[0,0]\right\},\left\{[0,1]\right\},\left\{[1,1]\right\},\left\{[0,2]\right\},\left\{[1,2]\right\},\left\{[2,2]\right\},\{[0,0],[2,2]\}\}$$

For parity polyominoes, we are only interested in those subsets whose intervals have even size and the same parity for the starting position. Rather than rewriting PreLeftLetters from scratch, we can use a trick in Maple to subclass PreLeftLetters and alter its parameters without having to modify the original PreLeftLetters code:

```
read ANIMALS:

origPreLeftLetters := subs(PreLeftLetters=origPreLeftLetters,
                           eval(PreLeftLetters)):

PreLeftLetters := proc(a, b)
    local halfS, S, i, PreLet,  PreLet1:

    halfS := origPreLeftLetters(0, floor((b-a+1)/2)-1) minus {{}}:
    S := {{}}:
    for i from 1 to nops(halfS) do
        PreLet  := {seq([a+2*halfS[i][j][1], a+2*halfS[i][j][2]+1],
                    j=1..nops(halfS[i]))}:
        PreLet1 := {seq([PreLet[j][1]+1, PreLet[j][2]+1],
                    j=1..nops(PreLet))}:
        S := S union {PreLet}:
        if max(seq(PreLet1[j][2], j=1..nops(PreLet1))) <= b then
            S := S union {PreLet1}:
        fi:
    od:
    S:
end:

```

The first line loads Zeilberger’s original ANIMALS package. The second line makes a copy of PreLeftLetters; the subs function is designed to modify PreLeftLetters so that any recursive calls are made to the copy rather than the original. The remaining lines redefine PreLeftLetters so that it produces the type of intervals that are required by parity polyominoes. The new PreLeftLetters uses the original PreLeftLetters to generate a set of subsets (in interval notation) for half the required range (stored in halfS). The for loop goes through each interval and doubles its size; this results in intervals that have even size and even parity (stored in PreLet). To create intervals with even size and odd parity, each interval in PreLet is shifted down by one, and the result is stored in PreLet2. The final condition checks that the shifted intervals do not go out of the range $`[a,b]`$.

For example, PreLeftLetters(0,5) would call origPreLeftLetters(0,2) and get the subsets shown above. The for loop would then modify these subsets to generate the following result which represents the parity polyomino preletters in the range $`[0,5]`$:

$$\{\left\{\right\},\left\{[0,1]\right\},\left\{[1,2]\right\},\left\{[0,3]\right\},\left\{[1,4]\right\},\left\{[2,3]\right\},\left\{[3,4]\right\},\left\{[0,5]\right\},\left\{[2,5]\right\},\left\{[4,5]\right\},\{[0,1],[4,5]\}\}$$

Once the preletters are generated, the grammar is extended by determining which of the preletters are valid extensions of existing polyomino words; the ANIMALS routine that performs this check is PreLetToLet. For parity polyominoes, not only must the preletter satisfy the same connectivity constraints as for general polyominoes, but the parity of the preletter’s interval starting positions must be different than the letter it follows. As we did with PreLeftLetters, a simple extension to the original code is that is required:

```
origPreLetToLet := subs(PreLetToLet=origPreLetToLet,
                        eval(PreLetToLet)):

PreLetToLet := proc(Let, PreLet)
        if PreLet <> {} and Let[1][1][1] mod 2 = PreLet[1][1] mod 2 then
                0:
        else
                origPreLetToLet(Let, PreLet):
        fi:
end:

```

Since we know that every interval in a parity polyomino letter has the same starting position parity, it suffices to check whether or not the starting positions of the first intervals have differing parities; if that is the case, PreLetToLet calls the original ANIMALS routine to ensure that the connectivity constraints are satisfied.

In order to weight the transfer-matrix vertices so that polyominoes with a specific number of cells can be enumerated, the ANIMALS routines call Weight(Let) to determine the number of cells in the letter Let. For enumerating hexanimals, we want to weight the transfer-matrix vertices with the number of hexagonal cells. Since each hexagonal cell is represented by two polyomino cells, a letter’s weight in hexagonal cells is half its weight in polyomino cells:

```
origWeight := subs(Weight=origWeight,
                   eval(Weight)):

Weight := proc(Let)
        floor(origWeight(Let)/2):
end:

```

One final change to ANIMALS is needed in order to enumerate hexanimals; the routine Khaya(L) that computes the number of polyominoes with $`L`$ cells relies on symmetry that does not exist in the hexagonal lattice. Since a hexanimal with $`L`$ cells must map to a parity polyomino whose set representation has a maximum $`y`$ coordinate of $`2L`$, we can implement Khaya(L) by a call to GFseries:

```
Khaya := proc(L)
        [GFseries(2*L, L)]:
end:

```

## 5. A User’s Manual for HexANIMALS

The modifications to ANIMALS that are described in the previous section are contained in the package HexANIMALS. To run HexANIMALS, download it and ANIMALS from the first authors’ web site to a local directory and start Maple. Once in Maple, type: read HexANIMALS and follow the on-line help.

Excluding Khaya, all the generating function and series expansion routines in ANIMALS remain unchanged. GF(n,s) computes the generating function for hexanimals embedded in the square lattice of height $`n`$ (i.e., the height of the square lattice is comprised of $`n`$ cells). GFseries(n,L) computes the list of length $`L`$ whose $`k^{th}`$ term is the number of hexanimals with $`k`$ cells embedded in the square lattice of height $`n`$. Khaya(L) computes the list of length $`L`$ whose $`k^{th}`$ term is the number of hexanimals with $`k`$ cells. Gf and Gfseries are the analogs of GF and GFseries for hexanimals whose polyomino embeddings have height exactly $`n`$.

As an example of HexANIMALS usage, the call GF(1,s) returns $`0`$ since no hexanimal can be embedded in a square lattice of height $`1`$. GF(2,s) returns $`s`$ since there is exactly one hexanimal (the single cell) that can be embedded in a square lattice of height $`2`$. GFseries(3,5) returns $`[1,2,2,2,2]`$ since the hexanimals in the square lattice of height $`3`$ with more than one hexagonal cell are chains, and each chain has two orientations depending on the parity of the starting cell. In Table 1 we compute the generating function for hexanimals embedded in a square lattice of height $`n7`$ along with the first few terms in the expansion by hexagonal cell count. Once $`n13`$, it takes too long to compute GF(n,s) exactly, but one can go much further with GFseries(n,L).

## 6. Counting Locally Skinny Hexagonal Animals

If we relax the condition that the entire hexagonal polyominoes must fit within a square lattice of fixed height, and only require that each column have height less than some arbitrary constant, then we get what Zeilberger calls locally skinny hexagonal polyominoes (as opposed to globally skinny hexagonal polyominoes).

The code for computing the generating function and series expansions for locally skinny polyominoes is contained in Zeilberger’s FreeANIMALS package. Since the changes to FreeANIMALS that are required in order to enumerate hexagonal lattice animals are virtually the same as those made to ANIMALS, we omit the details that were discussed in the previous sections and only discuss the major differences.

FreeANIMALS creates a Combinatorial Markov Process with multiple edges between letters. Each edge represents one of the possible vertical offsets of the letters. In contrast, ANIMALS creates a Combinatorial Markov Process with only one edge between letters; thus, each letter represents a particular arrangement of cells with an intrinsic offset. As a result, PreLeftLetters in HexANIMALS had to create the two possible parity offsets for each letter whereas PreLeftLetters in HexaFreeANIMALS needs to only create a single instance of the letter and the code that generates the Combinatorial Markov Process will create its own instances of the shifted letter. The resulting PreLeftLetters routine is similar to the PreLeftLetters routine in HexANIMALS except that the parity shifting code is omitted:

```
read FreeANIMALS:

origPreLeftLetters := subs(PreLeftLetters=origPreLeftLetters,
                           eval(PreLeftLetters)):

PreLeftLetters := proc(a, b)
    local halfS, S, i, PreLet:

    halfS := origPreLeftLetters(0, floor((b-a+1)/2)-1) minus {{}}:
    S := {{}}:
    for i from 1 to nops(halfS) do
        PreLet := {seq([a+2*halfS[i][j][1],
                        a+2*halfS[i][j][2]+1], j=1..nops(halfS[i]))}:
        S := S union {PreLet}:
    od:
    S:
end:

```

In addition, FreeANIMALS contains the routine PreLeftLettersk(a,b,k) which computes all possible subsets of the integers in the range $`[a,b]`$ that are represented by exactly $`k`$ intervals. PreLeftLettersk is used to compute $`k`$-board polyominoes. For example, PreLeftLettersk(0,4,2) returns the following subsets:

$$\{\{[0,0],[2,2]\},\{[0,0],[2,3]\},\{[0,0],[3,3]\},\{[0,1],[3,3]\}\}$$

The modifications to convert PreLeftLettersk for use with HexaFreeANIMALS are analogous to the changes made to PreLeftLetters as shown in the following code:

```
PreLeftLettersk := proc(a, b, k)
    local halfS, S, i, PreLet:

    halfS := origPreLeftLettersk(0, floor((b-a+1)/2)-1, k) minus {{}}:
    S := {{}}:
    for i from 1 to nops(halfS) do
        PreLet := {seq([a+2*halfS[i][j][1],
                   a+2*halfS[i][j][2]+1], j=1..nops(halfS[i]))}:
        S := S union {PreLet}:
    od:
    S:
end:

```

The only other significant difference between HexANIMALS and HexFreeANIMALS is in the specification of height restrictions. When considering globally skinny hexanimals, the height restrictions were defined in terms of the square lattice in which the hexanimals were embedded; this was because the shifting of adjacent columns in hexanimals makes the definition of the height of a hexanimal in terms of hexagonal cells difficult. In the case of locally skinny hexanimals, the height restrictions apply to individual columns and are well-defined in terms of hexagonal cells. Since the underlying algorithms of HexaFreeANIMALS operate on parity polyominoes, the height restrictions for hexanimal columns are simply doubled when applied to parity polyominoes as shown in the following code:

```
origgf := subs(gf=origgf, eval(gf)):

gf := proc(n, s)
    origgf(2*n, s):
end:

origgfSeries := subs(gfSeries=origgfSeries, eval(gfSeries)):

gfSeries := proc(n, L)
    origgfSeries(2*n, L):
end:

origgfList := subs(gfList=origgfList, eval(gfList)):

gfList := proc(resh, s)
    origgfList([seq(2*resh[i], i=1..nops(resh))], s):
end:

origgfSeriesList := subs(gfSeriesList=origgfSeriesList,
                         eval(gfSeriesList)):

gfSeriesList := proc(resh, L)
    origgfSeriesList([seq(2*resh[i], i=1..nops(resh))], L):
end:

```

## 7. A User’s Manual for the HexaFreeANIMALS

The modifications to FreeANIMALS are contained in the package HexaFreeANIMALS. gf(n,s) computes the generating function for hexagonal polyominoes whose columns span $`n`$ hexagonal cells. gfSeries(n,L) computes the list of length $`L`$ whose $`k^{th}`$ term is the number of hexagonal polyominoes with $`k`$ cells whose columns span $`n`$ hexagonal cells. gfList(List,s) computes the generating function for polyominoes whose $`k`$-board columns span $`List[k]`$ hexagonal cells. gfSeriesList(List,L) computes the series expansion of gfList(List,s) up to $`L`$ terms.

As an example of the usage of HexaFreeANIMALS, the call gfList(,s) would compute the generating function for polyominoes whose $`1`$-board columns span $`7`$ hexagonal cells and whose $`2`$-board columns span $`5`$ hexagonal cells. gfSeriesList(,24) computes Sloane’s sequence A0059716 \[NJAS\] of $`1`$-board polyominoes.

The call gfSeriesList(,12) computes the first $`12`$ terms of the previously unknown sequence

$$1,3,11,44,186,814,3648,16611,76437,354112,1647344,7682237$$

that enumerates $`2`$-board hexanimals with a given hexagonal cell count, the analog of Sloane’s sequence A001170 for board-pair-pile polyominoes.

## 8. Conclusion

We have described the Maple packages HexANIMALS and HexaFreeANIMALS. HexANIMALS enumerates globally skinny hexagonal polyominoes in which each polyominoes is embedded in a horizontal strip of prescribed height. HexaFreeANIMALS enumerates locally skinny polyominoes in which the total height of each hexagonal polyominoes is unbounded but whose columns are composed of a bounded number of hexagonal cells. In addition, we also demonstrated used HexaFreeANIMALS to compute the first few terms of the number of board-pair-hexagonal animals.

## References

\[K\] David A. Klarner, Cell Growth Problems, Canad. J. Math. 19 (1967), pp. 851-863.

\[GV\] Markous Voge and Anthony J. Guttmann, On the number of hexagonal polyominoes, Theoret. Comp. Sci. 307 (2003), 433-53.

\[NJAS\] N.J. Sloane and S. Plouffe, The Encyclopedia of Integer Sequences, Academic Press, San Diego, 1995.

\[S\] Richard P. Stanley, ”Enumerative Combinatorics”, Vol. 1, Wadsworth, Monterey, California, 1986.

\[Z0\] Doron Zeilberger, Automated Counting of LEGO Towers, J. Difference Eq. 5 (1999), pp. 355-377.

\[Z1\] Doron Zeilberger, Symbol-Crunching with the Transfer-Matrix Method in Order to Count Skinny Physical Creatures, INTEGERS (http://www.integers-ejcnt.org) 0 (2000), A9.

\[Z2\] Doron Zeilberger, The Umbral Transfer-Matrix Method III, Counting Animals, New York J. of Mathematics 7 (2001), pp. 223-231.