format: text

text:
# SNe heating and the chemical evolution of the intra-cluster medium

## 1 Introduction

Clusters of galaxies are the largest virialized systems in the Universe and are used as probes for the properties of the cosmic large-scale structure. They arise from high peaks of the primordial density field and during gravitational collapse they collect material over a scale of about 10 Mpc. Their typical size is of about 1 Mpc and their mass ranges from $`10^{13}M_{}`$ (the group edge) to $`10^{15}M_{}`$ (rich clusters). About 10-15$`\%`$ of their total mass is in the form of diffuse hot gas in hydrostatic equilibrium within the cluster potential well at a temperature of about 1-10 keV. Observations show that the intra-cluster gas has a mean iron abundance of one third of the solar value (Renzini 1997; White 2000), with no evidence for evolution out to $`z0.4`$ (Matsumoto et al. 2000). This fact suggests that the cluster galaxies must have lost their nuclear processed gas into the intra-cluster medium (ICM). In recent years, measurements of the abundances of other heavy elements such as Si and O in the ICM have also become available. In particular, ASCA results (Mushotzky et al. 1996) suggested \[$`\alpha `$/Fe\]ICM $``$ +0.2 dex. Later Ishimaru and Arimoto (1997) reanalyzed the same data adopting the meteoritic solar value for the Fe abundance, as opposed to the photospheric one used by the previous authors, and concluded that \[$`\alpha `$/Fe\]ICM $``$ 0. The photospheric solar Fe value used to be less certain than the meteoritic one, although more recent measurements (Grevesse et al. 1996) of the photospheric Fe abundance are in agreement with the meteoritic one. Very recently, Tamura et al. (2001) derived \[$`\alpha `$/Fe\] ratios from XMM-Newton observation of the cluster Abell 496 which are again roughly solar. Therefore, we can conclude that the data show a roughly solar value for the $`\alpha `$-elements/Fe ratio in the ICM. This ratio represents a very strong constraint to model the chemical enrichment in galaxies and in the cluster gas. Abundance gradients in the ICM of some clusters have been revealed by studies (Finoguenov et al. 2000; White 2000; De Grandi & Molendi 2001) based on data from ASCA, ROSAT and Beppo-SAX. However, firm conclusions on abundance gradients in clusters are still premature since the abundance determinations are known to be sensitive to the assumed temperature distribution.

From the theoretical point of view the first attempts to model the chemical enrichment of the ICM were by Larson and Dinerstein (1975), Vigroux (1977) and Hinnes and Biermann (1980). Their models were quite simple and they did not consider the element Fe, but only the global metal content Z and did not integrate the galactic contributions over the cluster mass function. Matteucci & Vettolani (1988) started a more detailed approach by computing the evolution of Fe and $`\alpha `$-elements under the assumption that type Ia SNe are the major producers of Fe and by integrating the galactic contributions over the Schechter (1976) luminosity function of cluster galaxies. Their approach was followed by David et al. (1991), Arnaud et al. (1992), Renzini et al. (1993), Elbaz et al. (1995), Matteucci and Gibson (1995), Gibson and Matteucci (1997), Loewenstein and Mushotzky (1996), Martinelli et al. (2000), Chiosi (2000) among others. The majority of these papers assumed that galactic winds (mainly from ellipticals) are responsible for the ICM chemical enrichment. Alternatively, the abundances in the ICM could be due to ram pressure stripping (Hinnes and Biermann 1980) or to pre-galactic Population III stars (White & Rees 1978; Loewenstein 2001).

Early attempts to model the ICM thermodynamics were originally based on the assumption that the gas is only affected by gravitational processes (Kaiser 1986; Evrard 1991), such as adiabatic compression and accretion shocks. Such models, which assume the ICM to behave in a self-similar fashion, were soon recognized to fail in accounting for several observational facts, such as the shape and evolution of the $`L_XT`$ relation, and the gas density profiles. For example, gas heating only from gravitation predicts $`L_XT^2\left(1+z\right)^{3/2}`$, whereas observations indicate $`L_XT^3`$ for $`T3`$ keV (e.g. Arnaud & Evrard 1999), with negligible evolution out to $`z1`$ (Borgani et al. 2001a), and an even steeper slope for colder systems (Helsdon & Ponman, 2000). At the same time, the shallower gas profile for poor clusters and groups turns into an excess of ICM entropy in central regions with respect to the expectation from self-similar scaling (Ponman et al. 1999). This discrepancy points towards the need of heating the gas before the cluster collapses by some non-gravitational source. This heating would bring the gas on a higher adiabat preventing it from reaching high densities and, therefore, would suppress its X-ray emissivity. Both semi-analytical (Tozzi & Norman 2001; Loewenstein 2001) and numerical (Bialek et al. 2001; Borgani et al. 2001b) simulations including gas-preheating suggest that about 1 keV per particle of extra-energy is required to break self-similarity. The possibility that SNe can provide this extra-energy has been explored already by several authors (Bower et al. 2001; Wu et al. 2000; Kravtsov & Yepes 2000; Valageas & Silk 1999; cf also Voit & Bryan 2001). All of these studied concluded that SNe alone can hardly provide the necessary 1 keV per particle and that the energy input from Active Galactic Nuclei (AGN) is required. However, in all of these studies no detailed chemical evolution and stellar lifetimes nor the contribution from SNIa were taken into account.

In this paper we model in a self-consistent way both the chemical and thermodynamical history of the ICM. To do that we adopt detailed chemical evolution models for elliptical galaxies where SN-driven galactic winds are triggered by the energy that SNe of both types (II, Ia) inject into the interstellar medium. To this purpose we take into account in great detail both nucleosynthesis prescriptions and energy feed-back between SNe and ISM. In particular, we explore several recipes for SN feedback, in order to check whether SNe can supply the necessary extra energy to break the self-similarity. We then integrate the chemical and energetical contributions from different SNe over the cluster luminosity function and predict their evolution as a function of the cosmic time. In section 2 the chemical evolution model for the ellipticals is described. In section 3 we model the enrichment of the ICM as a function of the cosmic time. In particular, we use a phenomenological evolution for the cluster luminosity function instead of adopting semi-analytical modelling for the formation of galaxy clusters. In section 4 we compare model results with observations and we draw our main conclusions in section 5. In the following, conversion from cosmic time to redshift is done for a cosmological model with $`\mathrm{\Omega }_m=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$ and $`H_o=70\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$.

## 2 The chemical evolution model for ellipticals

The adopted chemical evolution models (both one-zone and multi-zone) for elliptical galaxies are those of Matteucci & Gibson (1995) and Martinelli et al. (1998). We recall here the main assumptions of these models. First of all, they belong to the category of the so-called monolithic models since they assume that the collapse of gas into the potential well of a dark matter halo involves a large gas cloud. This collapse occurs at high redshifts on free-fall timescales, and the star formation rate is quite high since the beginning, thus converting gas into stars before the gas has time to substantially cool and form a disc. This is particularly true for the one-zone model whereas in the multi-zone model the process of galaxy formation is more complex in the sense that the internal regions continue to form stars for a longer period than the external ones. In the multi-zone model of Martinelli et al. (1998) the galactic wind occurs first in the external galactic regions because of the shallower potential well as compared to the most internal ones. Then after the onset of galactic wind the star formation stops and the region evolves passively. The time of the onset of the galactic wind is therefore a crucial parameter in these models since it fixes the abundances in the galactic stars and those in the ISM. The time for the occurrence of the galactic wind is dictated by the condition that the thermal energy of the gas be equal to or larger than its potential energy. The basic equations of chemical evolution following the temporal behaviour of the abundances of several species (H, He, C, N, O, Ne, Mg, Si, S, Fe) can be found in the above mentioned papers and we do not recall them here.

### 2.1 The energetics

The condition for the onset of a galactic wind in galaxies can be written as:

$$\left(E_{th}\right)_{ISM}E_{Bgas}$$
(1)
where $`E_{Bgas}`$ is the binding energy of the gas (see section 2.4) and $`\left(E_{th}\right)_{ISM}=\left(E_{th}\right)_{SN}+\left(E_{th}\right)_W`$ indicates the thermal energy of the gas due to injection of energy from SNe and stellar winds, respectively.

The thermal energy of the gas due to both SN types (II, Ia) is defined as:

$$\left(E_{th}\right)_{SN}=_0^tϵ_{SN}\left(tt^{}\right)R_{SN}\left(t^{}\right)𝑑t^{}$$
(2)
where $`R_{SN}`$ represents the SN rate of either type II or type Ia objects (see later). Note that the computation of $`\left(E_{th}\right)_{SN}`$ after the onset of the galactic wind is performed by integrating from the time of occurrence of the wind $`t_{GW}`$ and not from zero. For computing the thermal energy of the gas due to SNe, Matteucci and Vettolani (1988) and Matteucci and Gibson (1995) assumed the formulation of Cox (1972) (hereafter Cox72) where the efficiency of energy transfer from SNe into the ISM is taken to be constant at a value:

$$ϵ_{SN}=0.72E_o\mathrm{erg}$$
(3)
for $`t_{SN}t_c`$ years, where $`t_c`$ is the cooling time, $`E_o=10^{51}`$ erg is the explosion energy and $`t_{SN}=tt^{^{}}`$ is the time elapsed from the SN explosion, whereas it evolves as:

$$ϵ_{SN}=0.22E_o\left(t_{SN}/t_c\right)^{0.62}\mathrm{erg}$$
(4)
for $`t_{SN}>t_c`$.

It is worth noting that with these prescriptions each SN deposits effectively into the ISM at the end of its evolution only few percents of $`E_o`$. The advantage of adopting the above formulation relative to a simple parametrization of $`ϵ_{SN}`$ is that we can take into account the time dependence of the SN feed-back since different supernova remnants (SNR) contribute a different amount of energy according to their evolutionary stage. However, this formulation is derived for isolated SNe: if SNe explode in associations, there is the possibility, at least for small galaxies, that the efficiency of energy transfer into the ISM can be much larger than estimated above, especially if there is complete overlapping of the SNR with the consequent formation of a super-bubble. The crucial parameter, in the above formulation, is the the cooling time of a SNR, $`t_c`$. In fact, this timescale regulates the efficiency of gas cooling as a function of time and influences the time for the onset of the galactic wind. Unlike the previous papers (Matteucci & Gibson 1995 and Martinelli et al. 1998) which adopted the old cooling time by Cox72, we adopt here the results about the SNR evolution in the ISM of Cioffi, McKee & Bertschinger (1988) (hereafter CMB) which suggest the following cooling time depending on the metallicity:

$$t_{cool}=1.4910^4ϵ_0^{3/14}n_0^{4/7}\zeta ^{5/14}yr,$$
(5)
where $`\zeta =Z/Z_{}`$, $`n_0`$ is the hydrogen number density, $`ϵ_0`$ is the energy released during a SN explosion in units of $`10^{51}\mathrm{erg}`$ and we take always $`ϵ_0=1`$.

The old (Cox72) and new (CMB) cooling times are compared in Fig. 1, where metallicity and density of the ISM evolve in a self-consistent way as functions of time. The new cooling time is about 3 times shorter than the older one after 0.1 Gyr from the beginning of galactic evolution, and soon after 0.2 Gyr the metallicity becomes oversolar and $`t_c`$ decreases by a factor of $`10`$ with the consequence of having a much faster cooling process.

The expression for the energy injected from stellar winds is:

$$\left(E_{th}\right)_W=_0^t_{12}^{m_{up}}\phi \left(m\right)\psi \left(t^{^{}}\right)ϵ_W𝑑m𝑑t^{^{}}$$
(6)
where

$$ϵ_W=\eta _WE_W$$
(7)
where $`m_{up}100M_{}`$, $`E_W310^{47}`$ erg is the energy released into the ISM from a typical massive star ($`20M_{}`$) through stellar winds during its lifetime and $`\eta _W`$ is the efficiency of energy transfer, which we assume to be $`3\%`$ in according to Bradamante et al. (1998). This energy from stellar wind is only important before the onset of the first SN explosions after which it becomes negligible as already shown by Gibson (1994). He showed that for objects with luminous masses larger than $`10^9M_{}`$ the energetic input of stellar winds can be neglected, a result confirmed also by our calculations.

The initial mass function (IMF) is expressed by $`\phi \left(m\right)m^{\left(1+x\right)}`$, defined in the mass range $`0.1100M_{}`$, and we adopt two different prescriptions for it, as shown in Tables 1 and 2. The quantity $`\psi \left(t\right)`$ is the star formation rate (SFR) which is assumed to be:

$$\psi \left(t\right)=\nu M_{gas}\left(t\right)$$
(8)
with $`\nu `$ being the efficiency of star formation, namely the inverse of the typical timescale for star formation, and therefore expressed in units of $`\mathrm{Gyr}^1`$, and $`M_{gas}`$ the gas mass in the galaxy. The efficiency $`\nu `$ is chosen to reproduce the present time features of ellipticals. It is worth nothing, as shown by Matteucci (1994), that to reproduce the observed correlation between the magnesium index $`Mg_2`$ and the iron index $`<Fe>`$ in ellipticals one possibility is to assume an efficiency of star formation higher in more massive galaxies. This leads to the so-called “inverse wind” scenario, in the sense that large galaxies develop a wind before the small ones. From the point of view of the calculation of the metals and gas ejected into the ICM there is no difference between “classic” models (those where the galactic wind occurs first in small galaxies) and “inverse wind” models. Here we adopt the classic scenario.

The second novelty of our approach to compute the energetics is that we allow for two different efficiencies of energy transfer for different SN types. In particular, we explored the case in which SNe Ia are allowed to transfer all of their initial blast wave energy, namely:

$$ϵ_{SNIa}=E_0=10^{51}\mathrm{erg}.$$
(9)
The reason for this extremely efficient energy transfer resides in the fact that radiative losses from SNIa are likely to be negligible, since their explosions occur in a medium already heated by SNII (Bradamante et al. 1998; Recchi et al. 2001). In fact, type Ia SNe are thought to originate from long living systems, namely white dwarfs in binary systems. On the other hand, for SNe II, which explode first in a cold and dense medium, we allow for the cooling to be quite efficient. In particular, we adopt the efficiency of energy transfer deduced from the prescriptions of CMB, as described before. In the reality, the energy injected by type II SNe is perhaps higher, on average, than in CMB, especially for the last stellar generations which explode, like type Ia SNe, in an already heated and rarefied ISM. In this sense, the assumption on the energy transfer from SNe II should be regarded as a lower limit. However, this particular choice for the energetic prescriptions produces realistic galaxy models and larger amounts of energy from SNe should be rejected, as we discuss in the conclusions.

### 2.2 The stellar nucleosynthesis

Supernovae are also responsible for the chemical enrichment in heavy elements. In particular, SNe II, which should be the outcome of the explosion of single massive stars ($`M>8M_{}`$), are mostly responsible for producing the so-called $`\alpha `$-elements and part of Fe. On the other hand, the SN Ia, possibly originating from white dwarfs in binary systems exploding by C-deflagration, are mostly responsible for the production of Fe and iron-peak elements.

We adopted the following nucleosynthesis prescriptions:

* For low and intermediate mass stars ($`0.8M/M_{}8`$) the yields by Renzini and Voli (1981)
* For massive stars, type II SNe ($`M>10M_{}`$) the yields by Woosley and Weaver (1995), their case B.
* For SN Ia the yields by Nomoto et al. (1997)

### 2.3 The SN rates

As previously stated, the quantity $`R_{SN}`$ is defined as the rate of supernovae of either type Ia and II. The type Ia SN rate is computed as follows:

$$R_{SNIa}=A_{M_{Bm}}^{M_{BM}}\phi \left(M_B\right)_{\mu _m}^{0.5}f\left(\mu \right)\psi \left(t\tau _{M_2}\right)𝑑\mu 𝑑M_B$$
(10)
This formulation for the SN Ia rate was first adopted by Greggio & Renzini (1983) and Matteucci & Greggio (1986) and is based on the progenitor model made of a C-O white dwarf plus a red giant; $`M_\mathrm{B}`$ is the total mass of the binary system, $`M_{Bm}`$ and $`M_{BM}`$ are the minimum and maximum masses allowed for the adopted progenitor systems, respectively. We assume $`M_{Bm}=3M_{}`$ and $`M_{BM}=16M_{}`$. Also $`\mu =M_2/M_\mathrm{B}`$ is the mass fraction of the secondary, $`\mu _m`$ is the minimum value and $`f\left(\mu \right)`$ is the distribution function. Statistical studies (e.g. Tutukov & Yungelson 1979) indicates that mass ratios close to one are preferred, so the formula:

$$f\left(\mu \right)=2^{1+\gamma }\left(1+\gamma \right)\mu ^\gamma ,$$
(11)
is commonly adopted, where $`\gamma `$ is a parameter which we take equal to 2. This formulation for the SN Ia rate still seems to be the best to describe the chemical evolution of galaxies, as recently discussed by Matteucci & Recchi (2001). Finally, the parameter $`A`$ in eq. (10) represents the fraction of binary systems in the IMF which are able to give rise to SN Ia explosions and is a free parameter, fixed by reproducing the present time observed type Ia SN rate in elliptical galaxies.

As for the type II SN rate we can write:

$`R_{SNII}`$ $`=`$ $`\left(1A\right){\displaystyle _8^{16}}\psi \left(t\tau _m\right)\phi \left(m\right)𝑑m`$ (12)
$`+`$ $`{\displaystyle _{16}^{M_U}}\psi \left(t\tau _m\right)\phi \left(m\right)𝑑m`$
where the first integral accounts for the single stars in the range 8-16$`M_{}`$, and $`M_U`$ is the upper mass limit in the IMF.

The function $`\tau _m`$ represents the stellar lifetime of a star of mass $`m`$ and is expressed as:

$$\tau _m=10^{\left[1.338\sqrt{1.790.2232\left(7.764\mathrm{log}\left(m\right)\right)}\right]/0.11169}\mathrm{Gyr},$$
(13)
for the lifetimes of stars in the range 0.8-6.6 $`M_{}`$ (Padovani & Matteucci 1993), whereas we adopt:

$$\tau _m=1.2m^{1.85}+0.003\mathrm{Gyr}.$$
(14)
for stars with masses $`m>6.6\mathrm{M}_{}`$.

### 2.4 The potential well

The binding energy of the gas for the one-zone model is computed as in Matteucci & Gibson (1995):

$$E_{Bgas}\left(t\right)=W_L\left(t\right)+W_{LD}\left(t\right)$$
(15)
where :

$$W_L\left(t\right)=\frac{0.5GM_{gas}\left(t\right)M_{lum}\left(t\right)}{r_L}$$
(16)
is the gravitational energy of the gas due to the luminous matter, and:

$$W_{LD}\left(t\right)=\frac{GM_{gas}\left(t\right)M_{dark}}{r_L}\stackrel{~}{W}_{LD}$$
(17)
is the gravitational energy due to the dark matter, where $`r_L`$ is the effective radius and $`\stackrel{~}{W}_{LD}`$ is a term taking into account the distribution of the dark matter relative to the luminous one, in particular:

$$\stackrel{~}{W}_{LD}\frac{1}{2\pi }\frac{r_L}{r_D}\left[1+1.37\left(\frac{r_L}{r_D}\right)\right]$$
(18)
where $`\frac{r_L}{r_D}`$ is the ratio between the effective radius and the core radius of the dark matter. Realistic values for this ratio are $`0.5`$ (see Bertin et al. 1992). For all the ellipticals we set $`r_L/r_D=0.1`$ and $`\frac{M_{dark}}{M_{lum}}=10`$, since with these prescriptions one obtains realistic models for elliptical galaxies (Matteucci 1992).

For the multi-zone model we adopt the same formulation of Martinelli et al. (1998) where the galaxy is divided in several concentric shells and the potential well is calculated in each of them. In this way the external regions are less bound than the innermost ones thus creating a situation of “biased-galactic-wind” in the sense that the galactic wind starts first in the external regions and gradually develops in the more internal ones. This mechanism is in good agreement with the observed correlation between escape velocity and metallicity inside ellipticals which shows that where the escape velocity is larger then the metallicity is also larger (Carollo & Danziger 1994). The distribution of the dark matter is assumed to be the same as in the one-zone model, whereas the distribution of luminous matter follows that of Jaffe (1983). This kind of multizone model well reproduces both abundance (Martinelli et al. 1998) and color (Menanteau et al. 2001) gradients inside ellipticals.

## 3 Galactic models

We have computed several galactic models by varying some parameters such as the IMF and the SN feed-back. We will describe here only the best ones of each kind. We assume galaxies to have formed 13 Gyr ago corresponding to a formation redshift $`z_f=8.0`$ for the adopted cosmology. It is worth noting that all of the models described here reproduce the main features of elliptical galaxies, such as the present time type Ia SN rate, the metallicity of the dominant stellar population and the color-magnitude diagram, as already discussed in previous papers where we address the reader for details (Matteucci and Gibson, 1995; Matteucci, Ponzone & Gibson 1998; Martinelli, Matteucci & Colafrancesco, 1998; Menanteau et al. 2001). Here we discuss only the evolution of the abundances and energetic content of the ICM.

We discuss three models:

Model MG : the same as the best model of Matteucci & Gibson (1995). It is a one- zone model with Arimoto & Yoshii (1987, AY) IMF ($`\phi \left(m\right)m^{0.95}`$), cooling time and SNR evolution by Cox72 for both SN types. We run it for comparison with our models with new energetic prescriptions.

Model I : one- zone model with Salpeter (1955) IMF, CMB cooling time for type II SNe and SNIa without cooling.

Model II : multi-zone model (for details see Martinelli et al. 1998) with the same chemical and thermal prescriptions of model I.

In Tables 1 and 2 we show the model parameters: in column 1 is reported the initial galactic luminous mass, in column 2 the effective radius, in column 3 the star formation efficiency $`\nu \left(\mathrm{Gyr}^1\right)`$, in column 4 the adopted IMF, in column 5 the parameter A for type Ia SNe, in column 6 the type of feed-back (CMB or Cox72) for the SNII efficiency ($`ϵ_{SNII}`$) and in column 7 the SNIa efficiency ($`ϵ_{SNIa}`$). It is worth noting that the global efficiency of energy transfer, namely the computed final thermal energy of the gas relative to the initial total energy from all supernovae, turns out to be no more than $`20\%`$ for models I and II where the SN Ia deposit the entire energy budget into the ISM, whereas is only $`1.7\%`$ for the MG model. Another difference is that, thanks to the new energetic prescriptions, the energy provided by SNIa makes galactic winds continue until the present time, so models I and II can release larger masses of Fe and energy into the ICM than in the MG case.

In Tables 3 and 4 we show the model results: in particular, in column 1 we give the luminous galactic masses, in column 2 the time for the onset of the galactic wind $`t_{GW}`$, in column 3 the wind duration, in column 4 the ejected masses in the form of Fe, in column 5 the ejected masses in the form of oxygen and in column 6 the ejected total masses of gas. In column 7 the total thermal energy expressed in units of $`\mathrm{M}_{}\mathrm{pc}^2\mathrm{Gyr}^2`$ which correspond to $`1.810^{37}`$ erg. In column 8 we report the total number of SNe Ia ever exploded in each galaxy, in column 9 the number of SNe II and in column 10 the predicted type Ia SN rate at the present time. Note that the predicted SN Ia rate, expressed in units of $`SN\left(100\mathrm{y}\mathrm{r}\right)^110^{10}\mathrm{L}_\mathrm{B}`$, is in very good agreement with the average rate measured for ellipticals by Cappellaro et al. (1999) ($`R_{SNIa}=0.18\pm 0.06`$ SNu, for h=0.75). In Table 5 are shown the results of the multi-zone model, where each galaxy is divided into ten shells of width 0.1$`r_L`$. Here, for the sake of simplicity, we show only three regions: central (0-0.1$`r_L`$), middle (0.4-0.5$`r_L`$) and external (0.9-1$`r_L`$). It is worth noting that the total ejected masses and thermal energies by each galaxy in Model II are larger than the sums of the three regions shown in Table 5, which do not include the whole object, and are very similar to those of the corresponding galaxies of Model I. We do not show the total number of SNe and the average present time type Ia SN rate for Model II because they are very similar to the corresponding one-zone galaxy models (Model I). In Table 6 are shown the masses of Fe, O and metals which remain trapped into stars in cluster galaxies for Model I. The last two columns of Table 6 show the ratios $`R_O`$ and $`R_{Fe}`$ of the masses of O and Fe ejected by the galaxies relative to those locked up in stars, respectively.

As is shown in Tables 3, 4 and 5, models for large ellipticals ($`M_{lum}>10^{11}M_{}`$) with Salpeter IMF develop a galactic wind before models with a flatter IMF. This result is the consequence of the adopted cooling time of CMB which increases with metallicity, which in turn increases more rapidly in the case of a flat IMF (AY). Therefore, galaxies with a Salpeter IMF undergo galactic winds earlier and consequently eject more metals into the ICM. A similar effect was found by Gibson (1996).

## 4 Time evolution of Abundances and Energy

Firstly we make the assumption, supported by observational evidence (Arnaud et al. 1992), that only elliptical and perhaps S0 galaxies contribute to the chemical and thermal enrichment of the ICM. Then, to compute the total masses of the chemical elements, gas and total thermal energy ejected into the ICM by the cluster galaxies at any cosmic time (redshift) we find the parameters linking the luminous mass $`M_{lum}\left(z\right)`$ of the galaxy with the above quantities (shown in Tables 3, 4 and 5) via least-square fits of this kind (Matteucci & Vettolani 1988):

$$M_i^{ej}\left(z\right)=E_i\left(z\right)M_{lum}^{\beta _i\left(z\right)}\left(z\right)$$
(19)
where the subscript $`i`$ refers to a specific chemical element or to the total gas mass and $`E_i\left(z\right)`$ and $`\beta _i\left(z\right)`$ are parameters defined at any given redshift. The same procedure is applied to the total thermal energy of the gas ejected by cluster galaxies into the ICM via galactic winds. In particular, we derive the following expression for each set of models:

$$E_{th}\left(z\right)=A\left(z\right)M_{lum}^{\delta \left(z\right)}\left(z\right),$$
(20)
where $`z`$ is the redshift and $`\delta \left(z\right)`$ and $`A\left(z\right)`$ are parameters defined at any given redshift.

The total masses of metals and total gas as well as the total thermal energy injected at any time into the ICM are then given by the integrals of eq. (19) and (20) over the mass function of cluster galaxies. The mass function of cluster galaxies is obtained by means of the cluster luminosity function (LF hereafter) through the mass to light ratio. In particular, at each given cosmic time, the total mass ejected by cluster galaxies in the ICM is computed as:

$`M_{i,ICM}^{ej}\left(>M_{lum}\right)`$ $`=`$ $`E_in^{}\left(h^2k\right)^{\beta _i}10^{0.4\beta _i\left(M_K^{}2.63\right)}`$ (21)
$`\times `$ $`\left[\mathrm{\Gamma }(a,b)_{M_l}\mathrm{\Gamma }(a,b)_{M_u}\right]`$
where:

$$\mathrm{\Gamma }(a,b)_M=\mathrm{\Gamma }[\left(\alpha +1+\beta _i\right),\left(M_{lum}/h^2k\right)10^{0.4\left(M_K^{}2.63\right)}]$$
(22)
is the incomplete Eulerian $`\mathrm{\Gamma }`$ function. The quantities $`M_l`$ and $`M_u`$ represent the lowest and the largest luminous galactic masses used for the integration. We note that the baryonic masses (gas plus stars) of galaxies are changing with cosmic time either because gas is lost continuously for most of the models or because the mass of living stars decreases. Therefore, $`M_l`$ and $`M_u`$ are functions of time. The meaning of the various parameters is: $`h=H_0/100\mathrm{k}\mathrm{m}\mathrm{s}^1\mathrm{Mpc}^1`$ is the Hubble constant, $`\alpha `$ is the slope of the luminosity function, $`M_K^{}`$ is the magnitude in the K band at the break of the LF, $`n^{}`$ is the cluster richness and $`k=\frac{M}{L_K}`$ is the mass to light ratio. This quantity and in particular the $`L_K`$ luminosity is computed by means of the photometric model of Jimenez et al. (1998). We obtain a mass to light ratio M/$`L_K1`$ at z=0, in agreement with observations (e.g. Mobasher et al. 1999). The choice of the K luminosity was due to the fact that it does not vary dramatically with galaxy evolution as it is the case for B luminosity which is very sensitive to young and massive stars. To compute the total thermal energy we simply use eq. (21) by substituting $`M_{i,ICM}^{ej}`$ with $`E_{th}^{ej}`$ and $`\beta _i\left(z\right)`$ with $`\delta \left(z\right)`$. We integrate eq. (21) over the K–band LF, taking into account its evolution with redshift.

Since the global luminosity and mass of galaxies change in time, we follow the evolution of the K–band LF. To do that we adopt the observed B–band LF at $`z=0`$ by Sandage et al. (1985) (since no K-band LF is available for ellipticals in clusters), and apply the transformation from B– to K–band by Fioc & Rocca-Volmerange (1999), as well as the evolutionary corrections (Jimenez et al. 1998 and also Poggianti, 1997) for our assumed cosmological model. Finally, we also consider the possibility of morphological evolution of spiral galaxies into S0 galaxies for $`z0.4`$ (Butcher and Oemler 1978; Dressler et al. 1997; Fasano et al. 2000). To do that we simply assume that the fraction of S0 galaxies in clusters decreases by 12% every Gyr from z=0 up z=0.5 and remains constant afterwards. This assumption is based on results of Dressler et al. (1997) indicating that the fraction of S0 in clusters at z=0.5 is 2 to 3 times lower than at z=0.

It is worth reminding that we followed the evolution of the cluster galaxy population out to $`z=4.5`$, i.e. before the typical redshift of cluster formation. Therefore, our model predicts the thermal and chemical evolution of the diffuse IGM which would later ($`z<1`$) collapse to form the low–redshift ICM.

In order to derive model predictions for clusters of different $`X`$–ray temperature, $`T`$, we resort to a suitable recipe to relate the cluster richness, $`n^{}`$, to $`T`$. To this purpose, we first convert $`T`$ into the total cluster virial mass, by using the hydrostatic equilibrium relation, $`k_BTM^{2/3}`$, with normalization computed for spherical collapse and isothermal gas, as provided by eq. (2.2) of Eke, Cole & Frenk (1996), in the case of a fully thermalized gas. The total cluster mass is then converted into total cluster optical luminosity by using the relation $`M/LM^{0.2}`$, as found by Girardi et al. (2000), which gives $`ML^{1.3}\left(n^{}\right)^{1.3}`$. With this approach, results for a rich cluster, like Coma, are obtained by simply rescaling the corresponding richness:

$$n_{Coma}^{}=2.5n_{Virgo}^{}$$
(23)
## 5 Model results

In Table 7 we report the predicted total masses of Fe, O and gas plus the total thermal energy ejected by galaxies into the ICM of a poor (Virgo-like) and a rich (Coma-like) cluster. In the same Table are reported the observed total Fe and gas masses for Coma and Virgo. It is evident that our models, especially Model I and Model II, can well reproduce the total iron masses but underestimate the total gas masses. To compute the chemical abundances of the gas ejected by the cluster galaxies into the ICM we proceed as follows:

$$X_{ICM}^{el}=\frac{M_{el,ICM}^{ej}}{M_{gas,ICM}^{ej}}.$$
(24)
The abundances of Fe calculated in this way are 10-15 times larger that the solar Fe abundance. However, these abundances are not those which should be compared with the observed ones in the ICM, since the dilution due to the primordial gas (metal-free) present in the ICM is not yet taken into account. This is a well known result already discussed by previous authors (Matteucci & Vettolani, 1988; David et al. 1991; Renzini et al. 1993; Matteucci & Gibson 1995; Gibson & Matteucci 1997; Martinelli et al. 2000; Chiosi 2000). In fact, it was first pointed out by Matteucci & Vettolani (1988) that the cluster galaxies can well reproduce the total amount of Fe observed in clusters but not the total gas in the ICM. They concluded that a large fraction of the ICM has to be primordial, namely never processed inside stars, and that this fraction of primordial gas dilutes the high Fe abundances produced by the cluster galaxies down to the observed value of $`0.30.5Fe_{}`$. In fact, if we compute the ICM abundances by adopting the observed total ICM masses shown in Table 7 we obtain Fe abundances in very good agreement with observations (see Table 8). Moreover, the presence of such primordial component in the ICM is suggested by the observed ratio between the ICM mass and the mass of galaxies, which is $`5.45h^{3/2}`$ for the Coma cluster (White et al. 1993).

An important consideration deriving from eqs. (21) and (24) is that the ratio between the abundances of two chemical elements is independent of the total ICM mass and cluster richness. Another quantity which is independent of cluster richness and ICM mass is the iron mass to light ratio (IMLR), as defined by Renzini et al. (1993), and represents, together with the abundance ratios, a very good constraint to the evolution of galaxies and the ICM.

In Figs. 2 and 3 we show the iron abundances relative to the Sun, as predicted for the ICM in clusters of different richness, and therefore different temperature, at the present time. The amount of assumed dilution is exactly the difference between the predicted and the observed mass of the ICM. In each figure are shown several model predictions relative to either the one-zone or multi-zone models, with different prescriptions about IMF and morphological evolution of spirals into S0 (see captions). All models with AY IMF produce a too large iron abundance in the ICM, whereas models with Salpeter IMF give a good fit to the observations. In order to check our assumption about the high redshift of galaxy formation assumed in our models, we compare the observations with another one-zone model in which the galaxies formed in a ’hierarchical’ fashion. In particular we make our model galaxies start to form stars in the range of redshifts going from z=2.93 (z=2.92) for $`10^9M_{}`$ model, to z=2.75 (z=2.83) for $`10^{12}M_{}`$ model with Salpeter (AY) IMF. The epoch of galaxy formation has been chosen in order to satisfy the observations (e.g. Ellis et al. 1997) showing that ellipticals in dense clusters stopped most of their star formation before z$`2.53`$. In this way we have a small formation period (as suggested by the tightness of the color-magnitude relation, e.g. Bower et al. 1992) and, by imposing that all models undergo galactic wind at z=2.5, we can reproduce also the so-called ’inverse-wind model’ (Matteucci 1994), which is used to explain the overabundance of Mg relative to Fe in large ellipticals. Both Salpeter and AY models (showed in Fig. 2) predict ejected masses and energy smaller by a factor of $`2`$ than in the best case. The inverse wind scenario leads to a modest increase in \[O/Fe\] ($``$ -0.3 dex, to be compared with the values in Table 8) for Salpeter IMF, but underestimates the total amount of Fe in the ICM (see Fig. 2). Furthermore it predicts $`E_{pp}`$ lower by a factor of 2-3 than the one-zone best model. Similar conclusions can be drawn for multi-zone models. So we can strengthen our conclusion that the bulk of the ejection (in order to fit observations) has to take place at $`2.5<z<5`$ (as suggested also by Renzini 2000).

In Figs. 4 and 5 we show the total thermal energy of the gas ejected by cluster galaxies into the ICM as a function of the ICM temperature for the same models of Figs. 2 and 3. The total energy per particle is also computed by taking into account the observed mass of the ICM. For the “best models” (those with Salpeter IMF and evolution of spirals into S0), the energy per particle injected into the ICM by the cluster galaxies varies from 0.1 to a maximum of 0.2 keV for the one-zone model, and from 0.2 to 0.34 keV for the multi-zone model. It is worth noting that in Fig. 5 the predicted energy per particle for the multi-zone model with Salpeter IMF and no Spiral evolution is the highest, reaching values larger that 0.5 keV per particle. Unfortunately, this model predicts a too high Fe abundance in the ICM and should be rejected. Therefore, taking all the observational constraints into account the “best model” seems to be Model II (multi-zone) with evolution of spirals into S0.

As a final test for the energetics of the ICM we have computed the total energy budget due to SNe (II and Ia) in all galaxies in clusters by integrating over the LF the total numbers of SNe (II and Ia) exploded in each galaxy and shown in Tables 3 and 4. In particular, we take as an example Model I and obtain, for a Virgo-like cluster, a total number of SNe ever exploded of $`510^{10}`$, whereas for a Coma-like cluster a total number larger by a factor of ten. If we multiply these numbers by $`10^{51}`$ erg, which is the maximum energy contributed by a supernova, transform the result into keV and then divide by the total number of particles in the given cluster we obtain $`E_{pp}0.81\mathrm{k}\mathrm{e}\mathrm{V}`$, which is exactly what is required to break the self-similarity. Nearly the same $`E_{pp}`$ is obtained for a Virgo-like and for a Coma-like cluster, since in the Coma-like cluster the energy budget increases by a factor of ten but also the number of particles in the ICM increases by a similar factor. However, the assumption of the maximum efficiency for the SN energy transfer is not realistic since in this case the energy injected into the ISM of the galaxies is larger than their binding energy, thus inducing the complete loss of gas at very early times and preventing the formation of a realistic object. Therefore, we must conclude that, in order to obtain realistic galaxies, a substantial fraction of the total energy budget should be radiated away.

The same conclusion was reached by Renzini (2000) by taking into account simple arguments about the total Fe ejected by galaxies and the number of SNe required to produce it. However, all of these estimated energies for SNe are based on the assumption of thermal winds. If the wind velocity were $`3`$ times larger than the escape velocity more energy could be gained by the ICM (see Lloyd-Davies et al. 2000).

The best model predicts Fe abundances in the ICM in good agreement with the observational estimates, as shown in Table 8. In the same Table are shown the predicted \[$`\alpha `$/Fe\] ratios in the ICM both for poor and rich clusters. Our best model predicts \[O/Fe\] $`<0`$ and this is due to the Fe produced by type Ia SNe. In the MG model instead most of the Fe produced by type Ia SNe was not considered since galactic winds had a very short duration and the adopted IMF was flat (AY IMF). The reason for which we find here continuous wind in all the galaxies is due to the assumption that type Ia SNe can transfer all of their initial kinetic energy into the ICM. The predicted \[Si/Fe\] is higher than \[O/Fe\] and more in agreement with the observational estimates which seem to suggest a roughly solar value, although large uncertainties are still present in abundance determinations. Very recent results show that clusters dominated by a cD galaxy present strong element gradients and \[O/Fe\] less than zero. In particular, for Virgo cluster, Gastaldello & Molendi (2002) found \[O/Fe\] going from -0.52 (in the inner zone of about 2.5 kpc of radius) up to -0.28 (at about 70 kpc from the center, where the effect due to the presence of M87 is less important). These results agree with our predictions, confirming that the central region of cluster of galaxies is affected by continuous galactic winds (and, perhaps, by ram-pressure) leading to \[$`\alpha `$/Fe\] $`<0`$.

In Fig. 6 we show the evolution of the \[O/Fe\] ratio in the ICM for any cluster as a function of redshift (upper panel) and the \[O/Fe\] vs. \[Fe/H\] diagram relative to the ICM (lower panel), as predicted by the best model (solid line). As one can see, the \[O/Fe\] ratio in the ICM decreases strongly until z $`2.5`$ and flattens for low z values, thus predicting no evolution between z=1 and z=0, in agreement with observations (Matsumoto et al. 2000). The same trend (but with even stronger evolution at z $`2`$) is predicted by the model with variable redshift of galaxy formation and Salpeter IMF (dashed line). In this model \[O/Fe\] ratio starts from higher values than in the former case, due to the shorter duration of star formation in the most massive ellipticals. On the other hand, the trend predicted by the model with variable redshift of galaxy formation and AY IMF (dotted line) is flat and predicts over-solar values. It is worth noting that the lower panel in Fig. 6 represents the \[O/Fe\] vs. \[Fe/H\] for the ICM as predicted by our best model. This is the equivalent to the \[O/Fe\] vs. \[Fe/H\] in the solar neighborhood where the \[O/Fe\] ratio shows a plateau at low \[Fe/H\] followed by a linear decline towards the solar value of \[Fe/H\]. The \[O/Fe\] ratio in the ICM does not show any plateau, at variance with the solar neighborhood, owing to the fact that the contribution of type Ia SNe is relevant starting from the time of the onset of the galactic winds.

In Fig. 7 we show the predicted evolution with redshift of the total energy content (upper panel) and of the Fe mass (lower panel) in the ICM. In each panel are shown the different contributions of type II and Ia supernovae. It is evident from the first panel that the energy contribution to the ICM from type II SNe is negligible relative to the contribution of type Ia SNe, starting from the time when the galactic winds from all galaxies are important. In fact, type Ia SNe predominate both in producing energy and Fe, as evident from the lower panel. In Fig. 7 is shown also the expected $`E_{pp}`$ if an efficiency of 100% for all SN types is assumed.

Finally, we compute the IMLR for all the galaxies in clusters both for the iron which remains locked up in stars and for the iron ejected into the ICM. As it can be seen from the results of Model I (Table 6), the ratio $`R_O`$ (mass of O in stars divided by mass of O ejected into the ICM) varies roughly from to 2 to 0.6 from low mass to high mass galaxies whereas the ratio $`R_{Fe}`$ (mass of Fe in stars divided by mass of Fe ejected into the ICM) varies from 10 to 2. When integrated over the LF the total $`R_O1`$ and $`R_{Fe}5`$. This result is the consequence of the fact that most of the Fe is produced after the star formation has stopped, whereas for the oxygen is the contrary. This creates the “asymmetry” between galactic stars and ICM in the sense that we expect to find enhanced \[$`\alpha `$/Fe\] ratios in the stars and underabundances of $`\alpha `$-elements in the ICM. When the total mass of Fe in the ICM is then divided by the total blue luminosity of cluster galaxies, we obtain an ICM IMLR in the range 0.02-0.03 for all clusters, in very good agreement with the observational estimate (e.g. Renzini, 1997). The IMLR computed by means of the total mass of Fe locked up inside stars is obviously lower by a factor of 5. On the other hand, all the metals locked up in stars are comparable with all the metals ejected into the ICM, since the global metallicity is dominated by oxygen.

## 6 Conclusions

In this paper we have computed the evolution of the abundances of Fe, $`\alpha `$-elements, total mass and thermal energy of the gas ejected by the cluster galaxies into the ICM. The main assumptions are that only ellipticals and S0 galaxies pollute the ICM and that SNe II and Ia are the major responsible of the chemical enrichment and the non-gravitational heating of the ICM. To do that, we have adopted chemical evolution models for elliptical galaxies which reproduce the majority of the observational constraints and predicted the chemical enrichment of the ICM by integrating the contribution of cluster galaxies over the cluster luminosity function. Several prescriptions for the energy feed-back between SNe and ISM have been adopted in the galaxy models and their effects on the chemical and thermal history of the ICM have been studied. Relative to previous models we have considered the energy feed-back in more detail and included the effect of SNe Ia, often neglected in this kind of calculations. In particular, we have assumed that type Ia SNe can inject all of their initial blast wave energy into the ISM. We have allowed for the variation of the cluster luminosity function as a function of redshift and computed the evolution of the ICM abundances in time. The same procedure has been applied to the total thermal energy contributed by cluster galaxies through galactic winds.

Our main conclusions can be summarized as follows:

* Both one- and multi-zone models with Salpeter IMF and evolution of Spirals into S0 galaxies can reproduce the observed Fe abundance and the \[$`\alpha `$/Fe\] ratios in the ICM. In particular, there is good agreement between the predicted and observed \[Si/Fe\] ratios whereas the predicted \[O/Fe\] ratios are lower than indicated by some observations. However, the O measurements in the ICM are still very sparse and uncertain and more data are required to assess this point (but see Gastaldello & Molendi 2002 observations for Virgo). From the theoretical point of view we expect that \[O/Fe\] $`<`$ \[Si/Fe\], since Si is also produced by type Ia SNe whereas O is not.
* We predict an ICM iron-mass-to-light-ratio $`IMLR0.020.03`$ for all clusters, in very good agreement with the observed one (e.g. Renzini, 1997). We have also computed the amount of metals which remains locked up inside stars. We found that the Fe in the ICM should be 5 times more than the Fe in stars, whereas O and the global metal content should be roughly the same inside stars and in the ICM.
* Our best model is a multi-zone one (Model II), where the galaxy is divided in several shells developing the galactic wind first in the external regions and then in the internal ones.
* The best model, which reproduces the ICM abundances, can provide 0.20-0.35 keV of energy per ICM particle, depending on the cluster richness. This result has been obtained for a total efficiency of energy transfer from SNe II plus SN Ia of $`20\%`$.
* Little evolution is found both for the abundances and for the heating energy from z=0 up to z=1, in agreement with observations (Matsumoto et al. 2000).
* We predict \[$`\alpha `$/Fe\] ratios for the intergalactic medium at high redshift ($`z>3`$) to be near zero (solar) or slightly negative, indicating that already at those early times a considerable amount of Fe should be present. The existing uncertanties on SNII yields would affect slightly the value of the final \[$`\alpha `$/Fe\] but not the substance of our reasoning, as already stressed by Renzini (2000). In fact, the strong point is that the \[$`\alpha `$/Fe\] ratios are predicted to be negative or solar at the present time due to the Fe which is restored continuously as opposed to the $`\alpha `$-elements which stop to be produced starting from the time of the occurrence of the galactic winds.
* In both old and new models, while type II SNe dominate the chemical evolution of the ellipticals and are the main responsible for the onset of a galactic wind, type Ia SNe play a fundamental role in providing energy ($`8095\%`$) and Fe ($`4580\%`$) into the ICM. Again, uncertainties of a factor of 2 in SNII yields for Fe, would translate only in a $`20\%`$ (i.e. $`0.1`$ dex) uncertainty in the total Fe mass ejected into the ICM in our best model, without affecting our results. Therefore, SNIa cannot be neglected in computing the chemical and thermal evolution of the ICM.

In view of the relevance of the SN energy feedback for the thermal history of the ICM, it is worth discussing how robust is our determination of $`1/3`$ keV per particle.

By assuming an efficiency of energy transfer of 100$`\%`$ from both SN types, we obtain a significantly larger energy budget of $`1`$ keV per particle. However, with such an extreme efficiency the energy injected into the ISM is larger than the galaxy binding energy. Therefore, the natural conclusion is that a large fraction of the SN energy budget should be radiatively lost (see also Renzini et al. 1993). We estimate that at most $`35\%`$ of the overall SN efficiency is allowed in order not to unbind gas from the galaxy potential well. As an alternative route, we try to incorporate the effect of hypernovae, namely supernovae providing $`E_010^{52}`$ erg, which are believed to originate from stars with $`M>25M_{}`$. However, such stars are rare enough ($`510\%`$) to give only a modest increase of the ICM energy input from SNe. In summary, these results show that the overall efficiency of energy transfer from SNe into the ICM ranges between 20% and 35%, within the model uncertainties, so that $`E_{pp}0.4`$ keV.

Is this energy enough to provide the ICM extra–heating implied by X–ray observations of galaxy clusters? Both analytical methods (e.g. Cavaliere, Menci & Tozzi, Balogh et al. 1999, Tozzi & Norman 2001) and numerical simulations (e.g. Bialek et al. 2000, Borgani et al. 2001b) indicate that $`1`$ keV per particle of extra energy is required in order to reproduce the observed slope of the $`L_X`$$`T`$ relation (e.g. Arnaud & Evrard 1999) and the excess entropy in central regions of poor clusters and groups (e.g. Ponman et al. 1999). However, this energy is determined within a factor two, owing to the uncertainties in the detail of the heating mechanism, the redshift at which the energy in dumped into the ICM and the local density of the targetted gas. Furthermore, the inclusion of radiative cooling has been also suggested to decrease the required energy (e.g. Voit & Bryan 2001). In fact, cooling causes high–density, low entropy gas to disappear from the diffuse ICM phase, thus leading to an effective increase of the entropy and to a decrease of density of hot gas, similar to non–gravitational heating. Therefore, within the uncertainty in our capability of treating the complex ICM physics, our results should be taken only as suggestive that SNe may not provide enough extra energy. In this case, alternative astrophysical sources for ICM heating should be devised, such as energy release from galactic nuclear activity (AGN; e.g. Valageas & Silk 1999; Wu, Fabian & Nulsen 2000).

Available and forthcoming observations from Chandra and XMM satellites are expected to provide a radical change in our view of the ICM. The availability of ICM observations with good spatial and spectral resolution is now opening the possibility of tracing the pattern of metal enrichment and, therefore, of better connecting the ICM chemistry to the star formation history in clusters. This calls for the need of developing more sophisticated modelizations, which are able to join a detailed treatment of chemical evolution with a careful description of ICM hydrodynamics within the cosmological framework of galaxy formation.

## Acknowledgements

We would like to thank the referee, Alexis Finoguenov, for his careful reading and his suggestions. We thank John Danziger, Simone Recchi, Daniel Thomas and Paolo Tozzi for useful comments.


text:
# Exact solution for random walks on the triangular lattice with absorbing boundaries

## 1 Introduction

The problem of a random walk on a two-dimensional lattice with a single interior source point and absorbing boundaries was first considered by Courant et al.  in 1928 where general properties of the solution were discussed. This problem was solved exactly in 1940  for the case of random walks on a square lattice with rectangular absorbing boundaries. The problem on a triangular lattice with finite absorbing boundaries was considered in 1963  where the exact solution was given for an approximation of the problem using straight boundaries rather than the true zig-zag boundaries of the triangular lattice. Indeed the authors of this paper remarked that “An explicit solution of the difference equation can hardly be obtained if these boundary conditions are used.” Other variants of the problem on the square lattice have been solved exactly  however the problem on the triangular lattice with true zig-zag triangular lattice boundaries has remained unsolved. In this paper we give the exact solution for this problem.

The problem of random walks on finite lattices with absorbing boundaries is fundamental to the theory of stochastic processes  and has numerous applications. These include potential theory , electrical networks , surface diffusion  and Diffusion-Limited Aggregation (DLA) . For example, the exact results for the square lattice problem  have been used to expedite the growth of large DLA clusters on the square lattice  and the formulae derived in this paper could similarly be used on the triangular lattice.

## 2 Field equations, boundary conditions and absorption probabilities

A schematic illustration showing random walk pathways (dashed lines) on an equi- angular triangular lattice is shown in Fig. 1. A natural co-ordinate system for the triangular lattice is to label the lattice vertices by the intersection points $`(p,q)`$ of horizontal straight lines $`p=0,1,2,\mathrm{}m+1`$ and slanted straight lines $`q=0,1,2,\mathrm{}n+1`$ – parallel to the $`p`$ axis shown in Fig. 1. In this co-ordinate system the expectation that a random walk on the lattice visits a site $`(p,q)`$, before exiting at a finite boundary, is given by the non-separable difference equation

$`F(p,q)`$ $`=`$ $`{\displaystyle \frac{1}{6}}[F(p,q1)+F(p,q+1)+F(p1,q)`$
$`+F(p+1,q)+F(p+1,q1)+F(p1,q+1)].`$
In this paper we adopt a different co-ordinate system, which may at first appear less natural but it has the advantage that the governing equations are separable.

In the co-ordinate system used here we label the vertices of the triangular lattice by the intersection points $`(p,q)`$ of horizontal straight lines $`p=0,1,2,\mathrm{}m+1`$ and vertical zig-zag lines $`q=0,1,2,\mathrm{}n+1`$ (see Fig. 2). In the absence of any source terms the expectation that a random walker visits a lattice site $`(p,q)`$ is given by the coupled homogeneous difference equations

$`F(p,q)`$ $`=`$ $`{\displaystyle \frac{1}{6}}[\widehat{F}(p1,q1)+F(p,q1)+\widehat{F}(p+1,q1)`$ (1)
$`+\widehat{F}(p+1,q)+F(p,q+1)+\widehat{F}(p1,q)],`$
$`\widehat{F}(p,q)`$ $`=`$ $`{\displaystyle \frac{1}{6}}[F(p1,q)+\widehat{F}(p,q1)+F(p+1,q)`$ (2)
$`+F(p+1,q+1)+\widehat{F}(p,q+1)+F(p1,q+1)],`$
where $`F(p,q)`$ denotes the field value at an even $`p`$ co-ordinate and $`\widehat{F}(p,q)`$ denotes the field value at an odd $`p`$ co-ordinate. This separation of even and odd field equations, which is not necessary for the problem on the square lattice, is central to our solution below.

Now consider a source term at $`(a,b)`$ and absorbing boundaries at $`p=0,m+1`$ and $`q=0,n+1`$. In Section 3 we give the solution for the case where $`m>1`$ is odd. The same methods can be used with slightly altered equations to obtain the solution for $`m`$ even. The special case of the strip with $`m=1`$ has been detailed elsewhere . Following McCrea and Whipple  we construct separate solutions, $`F_I(p,q),\widehat{F}_I(p,q)`$ for $`qb`$ and $`F_{II}(p,q),\widehat{F}_{II}(p,q)`$ for $`qb`$. The absorbing boundary conditions are thus given by

$`F(0,q)`$ $`=`$ $`0,`$ (3)
$`F(m+1,q)`$ $`=`$ $`0,`$ (4)
$`F_I(p,0)`$ $`=`$ $`0,`$ (5)
$`\widehat{F}_I(p,0)`$ $`=`$ $`0,`$ (6)
$`F_{II}(p,n+1)`$ $`=`$ $`0,`$ (7)
$`\widehat{F}_{II}(p,n+1)`$ $`=`$ $`0.`$ (8)
The ommision of a subscript $`I`$ or $`II`$ in the above indicates that the same equations are satisfied by both $`F_I`$ and $`F_{II}`$.

At $`q=b`$ we have the matching conditions

$`F_I(p,b)`$ $`=`$ $`F_{II}(p,b),`$ (9)
$`\widehat{F}_I(p,b)`$ $`=`$ $`\widehat{F}_{II}(p,b).`$ (10)
Finally, with the inclusion of the source term, we have the coupled inhomogeneous field equations at $`q=b`$,

$`6F_I(p,b)`$ $`=`$ $`6\delta _{p,a}((a+1)\text{mod}2)`$ (11)
$`+[\widehat{F}_I(p1,b1)+F_I(p,b1)+\widehat{F}_I(p+1,b1)`$
$`+\widehat{F}_I(p+1,b)+F_{II}(p,b+1)+\widehat{F}_I(p1,b)],`$
$`6\widehat{F}_I(p,b)`$ $`=`$ $`6\delta _{p,a}(a\text{mod}2)`$ (12)
$`+[F_I(p1,b)+\widehat{F}_I(p,b1)+F_I(p+1,b)`$
$`+F_{II}(p+1,b+1)+\widehat{F}_{II}(p,b+1)+F_{II}(p1,b+1)].`$
In Section III we present the solutions to Eqs. (1)–(12).

The probabilities for the random walking particle to be absorbed at any specified point on one of the four boundaries, $`p=0,p=m+1,q=0,q=n+1`$, are readily obtained by averaging over the nearest neighbour expectation values subject to the boundary conditions, Eqs. (3)–(8). For example,

$`P(2k1,0)`$ $`=`$ $`{\displaystyle \frac{1}{6}}\left[\widehat{F}_I(2k1,1)+F_I(2k2,1)+F_I(2k,1)\right];`$ (13)
$`k=1,\mathrm{}{\displaystyle \frac{m+1}{2}}.`$
## 3 Solution

### 3.1 Homogeneous equations – Separation of variables

We begin by solving the homogenous equations, Eqs. (1),(2), subject to the absorbing boundary conditions, Eqs. (3),(4). The homogenous field equations separate, on using

$`F(p,q)`$ $`=`$ $`P(p)Q(q),`$ (14)
$`\widehat{F}(p,q)`$ $`=`$ $`\widehat{P}(p)\widehat{Q}(q),`$ (15)
into

$`{\displaystyle \frac{6Q(q)Q(q+1)Q(q1)}{\widehat{Q}(q1)+\widehat{Q}(q)}}`$ $`=`$ $`{\displaystyle \frac{\widehat{P}(p1)+\widehat{P}(p+1)}{P(p)}}=\lambda ,`$ (16)
$`{\displaystyle \frac{6\widehat{Q}(q)\widehat{Q}(q+1)\widehat{Q}(q1)}{Q(q+1)+Q(q)}}`$ $`=`$ $`{\displaystyle \frac{P(p1)+P(p+1)}{\widehat{P}(p)}}=\kappa ,`$ (17)
where $`\kappa `$ and $`\lambda `$ are separation constants. The two coupled equations for $`P`$ and $`\widehat{P}`$ can be readily decoupled into separate equations

$`P(2k+2)+(2\lambda \kappa )P(2k)+P(2k2)=0,`$ (18)
$`\widehat{P}(2k+3)+(2\lambda \kappa )\widehat{P}(2k+1)+\widehat{P}(2k1)=0,`$ (19)
where $`k`$ is an integer. This second order linear difference equation for $`P,(\widehat{P})`$ on a lattice of even (odd) integers has the solution

$$P(2k)=A\mu ^k+B\mu ^k,$$

where

$$\mu =\frac{(2\lambda \kappa )+\sqrt{(2\lambda \kappa )^24}}{2},$$

and $`A`$ and $`B`$ are arbitrary. It follows that the solution of $`P(p)`$ (with $`p`$ even) corresponding to the boundary conditions, Eqs. (3),(4), is

$$P(p)=c\mathrm{sin}(\frac{\pi jp}{m+1})$$
(20)
with

$$\lambda \kappa =4\mathrm{cos}^2(\frac{\pi j}{m+1}).$$
(21)
The solution for $`\widehat{P}(p)`$ (with $`p`$ odd) now follows from Eq. (17) as

$$\widehat{P}(p)=\frac{2c}{\kappa }\mathrm{cos}(\frac{\pi j}{m+1})\mathrm{sin}(\frac{\pi jp}{m+1}).$$
(22)
Equations (16) and (17) can also be decoupled for $`Q(q)`$ and $`\widehat{Q}(q)`$ resulting in the same fourth order linear difference equation in each case;

$`Q(q+4)(12+\lambda \kappa )Q(q+3)+(382\lambda \kappa )Q(q+2)`$
$`(12+\lambda \kappa )Q(q+1)+Q(q)=0.`$ (23)
The zeroes of the corresponding characteristic quartic polynomial are in reciprocal pairs leading to the solutions

$`Q(q)`$ $`=`$ $`c_1\mathrm{e}^{\alpha q}+c_2\mathrm{e}^{\alpha q}+c_3\mathrm{e}^{\beta q}+c_4\mathrm{e}^{\beta q},`$ (24)
$`\widehat{Q}(q)`$ $`=`$ $`\widehat{c}_1\mathrm{e}^{\alpha q}+\widehat{c}_2\mathrm{e}^{\alpha q}+\widehat{c}_3\mathrm{e}^{\beta q}+\widehat{c}_4\mathrm{e}^{\beta q},`$ (25)
where

$`\mathrm{cosh}\alpha `$ $`=`$ $`3+{\displaystyle \frac{\lambda \kappa }{4}}+{\displaystyle \frac{1}{4}}\sqrt{(\lambda \kappa )^2+32\lambda \kappa },`$ (26)
$`\mathrm{cosh}\beta `$ $`=`$ $`3+{\displaystyle \frac{\lambda \kappa }{4}}{\displaystyle \frac{1}{4}}\sqrt{(\lambda \kappa )^2+32\lambda \kappa }.`$ (27)
The coefficients $`c_1,c_2,c_3,c_4`$ may be chosen arbitrarily with the coefficients $`\widehat{c}_1,\widehat{c}_2,\widehat{c}_3,\widehat{c}_4`$ then determined by substituting Eqs. (24), (25) into Eq. (16) or Eq. (17). From Eq. (17) we obtain

$$\widehat{c}_1=c_1f(\alpha ),\widehat{c}_2=c_2f(\alpha ),\widehat{c}_3=c_3f(\beta ),\widehat{c}_4=c_4f(\beta )$$
(28)
where

$$f(\omega )=\frac{\kappa (\mathrm{e}^\omega +1)}{2(3\mathrm{cosh}\omega )}.$$
(29)
By combining Eqs. (14),(15),(20),(22),(24),(25),(28),(29) we can write the solutions to the field equations, Eqs. (1), (2), that satisfy the boundary conditions, Eqs. (3),(4), in the form

$`F(p,q;A,B,C,D)={\displaystyle \underset{j=1}{\overset{m}{}}}\mathrm{sin}\left({\displaystyle \frac{\pi jp}{m+1}}\right)`$
$`\times [A_j\mathrm{e}^{\alpha q}(3\mathrm{cosh}\alpha )+B_j\mathrm{e}^{\alpha q}(3\mathrm{cosh}\alpha )`$
$`+C_j\mathrm{e}^{\beta q}(3\mathrm{cosh}\beta )+D_j\mathrm{e}^{\beta q}(3\mathrm{cosh}\beta )],`$ (30)
$`\widehat{F}(p,q;A,B,C,D)={\displaystyle \underset{j=1}{\overset{m}{}}}\mathrm{sin}\left({\displaystyle \frac{\pi jp}{m+1}}\right)\mathrm{cos}\left({\displaystyle \frac{\pi j}{m+1}}\right)`$
$`\times [A_j\mathrm{e}^{\alpha q}(\mathrm{e}^\alpha +1)+B_j\mathrm{e}^{\alpha q}(\mathrm{e}^\alpha +1)`$
$`+C_j\mathrm{e}^{\beta q}(\mathrm{e}^\beta +1)+D_j\mathrm{e}^{\beta q}(\mathrm{e}^\beta +1)],`$ (31)
where $`A_j,B_j,C_j,D_j`$ are arbitrary.

Before imposing the remaining absorbing boundary conditions, Eqs. (5)–(8), we represent the solutions in the two regions: $`I:qb`$; $`II:qb`$ by

$`F_I(p,q)`$ $`=`$ $`F(p,q;A,B,C,D),`$
$`\widehat{F}_I(p,q)`$ $`=`$ $`\widehat{F}(p,q;A,B,C,D),`$
$`F_{II}(p,q)`$ $`=`$ $`F(p,q;\stackrel{~}{A},\stackrel{~}{B},\stackrel{~}{C},\stackrel{~}{D}),`$
$`\widehat{F}_{II}(p,q)`$ $`=`$ $`\widehat{F}(p,q;\stackrel{~}{A},\stackrel{~}{B},\stackrel{~}{C},\stackrel{~}{D}).`$
The absorbing boundary conditions, Eqs. (5)–(8), can be used to eliminate $`C,D,\stackrel{~}{C},\stackrel{~}{D}`$. After an appropriate re-normalization of the coefficients $`A,B,\stackrel{~}{A},\stackrel{~}{B}`$ this yields

$`F_I(p,q)`$ $`=`$ $`{\displaystyle \underset{j=1}{\overset{m}{}}}\mathrm{sin}\left({\displaystyle \frac{\pi jp}{m+1}}\right)`$ (32)
$`\times \left[A_I(q,\alpha ,\beta )A_j+A_I(q,\alpha ,\beta )B_j\right],`$
$`\widehat{F}_I(p,q)`$ $`=`$ $`{\displaystyle \underset{j=1}{\overset{m}{}}}\mathrm{sin}\left({\displaystyle \frac{\pi jp}{m+1}}\right)\mathrm{cos}\left({\displaystyle \frac{\pi j}{m+1}}\right)`$ (33)
$`\times \left[\widehat{A}_I(q,\alpha ,\beta )A_j+\widehat{A}_I(q,\alpha ,\beta )B_j\right],`$
$`F_{II}(p,q)`$ $`=`$ $`{\displaystyle \underset{j=1}{\overset{m}{}}}\mathrm{sin}\left({\displaystyle \frac{\pi jp}{m+1}}\right)`$ (34)
$`\times \left[A_{II}(q,\alpha ,\beta )\stackrel{~}{A}_j+A_{II}(q,\alpha ,\beta )\stackrel{~}{B}_j\right],`$
$`\widehat{F}_{II}(p,q)`$ $`=`$ $`{\displaystyle \underset{j=1}{\overset{m}{}}}\mathrm{sin}\left({\displaystyle \frac{\pi jp}{m+1}}\right)\mathrm{cos}\left({\displaystyle \frac{\pi j}{m+1}}\right)`$ (35)
$`\times \left[\widehat{A}_{II}(q,\alpha ,\beta )\stackrel{~}{A}_j+\widehat{A}_{II}(q,\alpha ,\beta )\stackrel{~}{B}_j\right],`$
where

$`A_I(q,\alpha ,\beta )`$ $`=`$ $`(3\mathrm{cosh}\alpha )(3\mathrm{cosh}\beta )2\mathrm{sinh}\beta \mathrm{e}^{\alpha q}`$ (36)
$`(3\mathrm{cosh}\beta )\gamma (\alpha ,\beta )\mathrm{e}^{\beta q}`$
$`+(3\mathrm{cosh}\beta )\gamma (\alpha ,\beta )\mathrm{e}^{\beta q},`$
$`\widehat{A}_I(q,\alpha ,\beta )`$ $`=`$ $`2\mathrm{e}^{\alpha q}(\mathrm{e}^\alpha +1)(3\mathrm{cosh}\beta )\mathrm{sinh}\beta `$ (37)
$`\gamma (\alpha ,\beta )(\mathrm{e}^\beta +1)\mathrm{e}^{\beta q}+\gamma (\alpha ,\beta )(\mathrm{e}^\beta +1)\mathrm{e}^{\beta q},`$
$`A_{II}(q,\alpha ,\beta )`$ $`=`$ $`(3\mathrm{cosh}\alpha )(3\mathrm{cosh}\beta )2\mathrm{sinh}\beta \mathrm{e}^{\alpha q}`$ (38)
$`\mathrm{e}^{(\alpha \beta )(n+1)}\gamma (\alpha ,\beta )(3\mathrm{cosh}\beta )\mathrm{e}^{\beta q}`$
$`+\mathrm{e}^{(\alpha +\beta )(n+1)}\gamma (\alpha ,\beta )(3\mathrm{cosh}\beta )\mathrm{e}^{\beta q},`$
$`\widehat{A}_{II}(q,\alpha ,\beta )`$ $`=`$ $`(\mathrm{e}^\alpha +1)(3\mathrm{cosh}\beta )2\mathrm{sinh}\beta \mathrm{e}^{\alpha q}`$ (39)
$`\mathrm{e}^{(\alpha \beta )(n+1)}\gamma (\alpha ,\beta )(\mathrm{e}^\beta +1)\mathrm{e}^{\beta q}`$
$`+\mathrm{e}^{(\alpha +\beta )(n+1)}\gamma (\alpha ,\beta )(\mathrm{e}^\beta +1)\mathrm{e}^{\beta q},`$
and the function

$`\gamma (\alpha ,\beta )`$ $`=`$ $`4\mathrm{cosh}\alpha 4\mathrm{cosh}\beta (3\mathrm{cosh}\beta )\mathrm{sinh}\alpha `$ (40)
$`(3\mathrm{cosh}\alpha )\mathrm{sinh}\beta .`$
The coefficients defined in Eq. (36) which are shown as implicit functions of $`\alpha `$ and $`\beta `$ are also functions of $`j`$ via Eqs. (21),(26),(27). Equations (32)–(35), are general solutions to the coupled homogeneous field equations, Eqs. (1),(2), that satisfy all of the absorbing boundary conditions, Eqs. (3)–(8).

### 3.2 Inhomogeneous equations – Matching conditions

By using the two matching conditions, Eqs. (9),(10), the four arbitrary constants, $`A,B,\stackrel{~}{A},\stackrel{~}{B}`$, can be reduced to two arbitrary constants, $`A,B`$ say. The solutions in Region II can then be written as

$`F_{II}(p,q)`$ $`=`$ $`{\displaystyle \underset{j=1}{\overset{m}{}}}\mathrm{sin}\left({\displaystyle \frac{\pi jp}{m+1}}\right)\left[A_{II}^{}(q,\alpha ,\beta )A_j+A_{II}^{}(q,\alpha ,\beta )B_j\right],`$ (41)
$`\widehat{F}_{II}(p,q)`$ $`=`$ $`{\displaystyle \underset{j=1}{\overset{m}{}}}\mathrm{sin}\left({\displaystyle \frac{\pi jp}{m+1}}\right)\mathrm{cos}\left({\displaystyle \frac{\pi j}{m+1}}\right)`$ (42)
$`\times \left[\widehat{A}_{II}^{}(q,\alpha ,\beta )A_j+\widehat{A}_{II}^{}(q,\alpha ,\beta )B_j\right],`$
where

$`A_{II}^{}(q,\alpha ,\beta )=A_{II}(q,\alpha ,\beta )\mathrm{\Gamma }_1+A_{II}(q,\alpha ,\beta )\mathrm{\Gamma }_2,`$ (43)
$`\widehat{A}_{II}^{}(q,\alpha ,\beta )=\widehat{A}_{II}(q,\alpha ,\beta )\mathrm{\Gamma }_1+\widehat{A}_{II}(q,\alpha ,\beta )\mathrm{\Gamma }_2,`$ (44)
and

$`\mathrm{\Gamma }_1={\displaystyle \frac{\widehat{A}_{II}(b,\alpha ,\beta )A_I(b,\alpha ,\beta )A_{II}(b,\alpha ,\beta )\widehat{A}_I(b,\alpha ,\beta )}{A_{II}(b,\alpha ,\beta )\widehat{A}_{II}(b,\alpha ,\beta )\widehat{A}_{II}(b,\alpha ,\beta )A_{II}(b,\alpha ,\beta )}},`$ (45)
$`\mathrm{\Gamma }_2={\displaystyle \frac{\widehat{A}_{II}(b,\alpha ,\beta )A_I(b,\alpha ,\beta )A_{II}(b,\alpha ,\beta )\widehat{A}_I(b,\alpha ,\beta )}{A_{II}(b,\alpha ,\beta )\widehat{A}_{II}(b,\alpha ,\beta )A_{II}(b,\alpha ,\beta )\widehat{A}_{II}(b,\alpha ,\beta )}},`$ (46)
Finally the remaining arbitrary constants $`A,B`$ are determined from the requirement that the solutions satisfy the coupled inhomogeneous field equations, Eqs. (11),(12). This step is facilitated using the identity

$$\delta _{p,a}=\frac{2}{m+1}\underset{j=1}{\overset{m}{}}\mathrm{sin}\left(\frac{\pi ja}{m+1}\right)\mathrm{sin}\left(\frac{\pi jp}{m+1}\right).$$
(47)
Explicitly we find:
i) $`a`$ even,

$`A_j={\displaystyle \frac{T_1(j)}{T_2(j)}}B_j,`$ (48)
$`B_j={\displaystyle \frac{12}{m+1}}\mathrm{sin}\left({\displaystyle \frac{\pi ja}{m+1}}\right){\displaystyle \frac{T_2(j)}{T_1(j)T_3(j)T_4(j)T_2(j)}};`$ (49)
ii) $`a`$ odd,

$`A_j={\displaystyle \frac{T_4(j)}{T_3(j)}}B_j,`$ (50)
$`B_j={\displaystyle \frac{12}{m+1}}\mathrm{sin}\left({\displaystyle \frac{\pi ja}{m+1}}\right){\displaystyle \frac{T_3(j)}{T_2(j)T_4(j)T_1(j)T_3(j)}},`$ (51)
and

$`T_1(j)`$ $`=`$ $`2A_I(b,\alpha ,\beta )+2A_{II}^{}(b+1,\alpha ,\beta )+\widehat{A}_I(b1,\alpha ,\beta )`$ (52)
$`+\widehat{A}_{II}^{}(b+1,\alpha ,\beta )6\widehat{A}_I(b,\alpha ,\beta ),`$
$`T_2(j)`$ $`=`$ $`6\widehat{A}_I(b,\alpha ,\beta )2A_I(b,\alpha ,\beta )2A_{II}^{}(b+1,\alpha ,\beta )`$ (53)
$`\widehat{A}_I(b1,\alpha ,\beta )\widehat{A}_{II}^{}(b+1,\alpha ,\beta ),`$
$`T_3(j)`$ $`=`$ $`6A_I(j,b)2\mathrm{cos}^2({\displaystyle \frac{\pi j}{m+1}})\widehat{A}_I(b1,\alpha ,\beta )`$ (54)
$`2\mathrm{cos}^2({\displaystyle \frac{\pi j}{m+1}})\widehat{A}_I(b,\alpha ,\beta )A_I(b1,\alpha ,\beta )`$
$`A_{II}^{}(b+1,\alpha ,\beta ),`$
$`T_4(j)`$ $`=`$ $`2\mathrm{cos}^2({\displaystyle \frac{\pi j}{m+1}})\widehat{A}_I(b1,\alpha ,\beta )`$ (55)
$`+2\mathrm{cos}^2({\displaystyle \frac{\pi j}{m+1}})\widehat{A}_I(b,\alpha ,\beta )+A_I(b1,\alpha ,\beta )`$
$`+A_{II}^{}(b+1,\alpha ,\beta )6A_I(b,\alpha ,\beta ).`$
Our solution for the triangular lattice site expectation values with absorbing boundary conditions and a source point is finally given by Eqs. (32),(33) in Region I and Eqs. (41),(42) in Region II. The relevant quantities appearing in these equations are defined through the series of equations, Eqs. (21),(26),(27),(36)–(40),(43)–(46),(48)–(50).

## 4 Example

Consider the case of a triangular lattice with $`m=7,n=7`$ and a source at $`a=4,b=4`$ (Fig. 2). In the approximation to this problem using straight edge boundaries and the $`(p,q)`$ co-ordinate system of Keberle and Montet  this problem corresponds to; $`m=16,n=7,a=8,b=4`$.

We have calculated expectation values at the nearest neighbour lattice sites around the source and absorption probabilities at the boundaries; (a) using our exact results above, and (b) using the results of Keberle and Montet  for the approximation to the problem – their equations (7a), (7b), (7c). In this small model lattice system we found that the straight line boundary approximation provides reasonable results for expectation values near the source (accurate to within a few percent) but provides poor results for the absorption probabilities (Table 1). Note that the ‘absorption probabilities’ for the zig-zag boundary co-ordinates which are inside the straight line boundary of Keberle and Montet  are not true probabilities in their solution.

## 5 Conclusion

Although our solution is somewhat unwieldly, we have nevertheless solved the underlying field equations for random walks on a finite triangular lattice with a single interior source point and zig-zag absorbing boundaries and have thus calculated the associated absorption probabilities. This problem was previously considered to be intractable . We hope that our result will inspire further work in this area.

We thank Michael Barber for first drawing our attention to  and Wolfgang Schief for his warm oral translation of . MTB has been supported by The Australian Research Council.

## References