text:
# Three-body recombination in Bose gases with large scattering length

## Abstract

An effective field theory for the three-body system with large scattering length is applied to three-body recombination to a weakly-bound $`s`$-wave state in a Bose gas. Our model independent analysis demonstrates that the three-body recombination constant $`\alpha `$ is not universal, but can take any value between zero and $`67.9\mathrm{}a^4/m`$, where $`a`$ is the scattering length. Other low-energy three-body observables can be predicted in terms of $`a`$ and $`\alpha `$. Near a Feshbach resonance, $`\alpha `$ should oscillate between those limits as the magnetic field $`B`$ approaches the point where $`a\mathrm{}`$. In any interval of $`B`$ over which $`a`$ increases by a factor of 22.7, $`\alpha `$ should have a zero.

preprint: DOE/ER/40561-84-INT00
The successful achievement of Bose-Einstein condensation has triggered a large interest in interacting Bose gases. One of the main factors limiting the achievable density in these experiments is the loss of atoms through 3-body recombination. Such losses occur when three atoms scatter to form a molecular bound state (called a “dimer” for brevity) and a third atom. The kinetic energy of the final state particles allows them to escape from the trapping potential. This 3-body recombination process is interesting in its own right as it provides a unique window on 3-body dynamics.

The number of recombination events per unit time and volume can be parametrized as $`\nu _{rec}=\alpha n^3`$, where $`\alpha `$ is the recombination constant and $`n`$ the density of the gas. The calculation of $`\alpha `$ in general is a complicated problem, because it is sensitive to the detailed behavior of the interaction potential . The simplest case is 3-body recombination to a weakly-bound $`s`$-wave state. For atoms of scattering length $`a`$ and mass $`m`$, the binding energy of the dimer is $`B_d=\mathrm{}^2/(ma^2)`$ and the size of the bound state is comparable to $`a`$. Assuming that $`a`$ is the only important length scale, dimensional analysis implies $`\alpha =𝒞\mathrm{}a^4/m`$, where $`𝒞`$ is dimensionless. The problem of 3-body recombination to a weakly-bound $`s`$-wave state has been studied previously . Fedichev et al.  found that the coefficient $`𝒞`$ has the universal value $`𝒞=3.9`$. Nielsen and Macek  and Esry et al.  found that $`𝒞`$ could take any value between $`0`$ and $`𝒞_{\mathrm{max}}65`$.

An ideal means to study the dependence of $`\alpha `$ on the scattering length $`a`$ is to use Feshbach resonances , which occur when the energy of a spin-singlet molecular bound state is tuned to the energy threshold for two separated atoms by applying an external magnetic field $`B`$. Such resonances have, e.g., been observed for 23Na and 85Rb atoms . When the magnetic field is varied in the vicinity of the resonance, the scattering length varies according to

$$a(B)=a_0\left(1+\frac{\mathrm{\Delta }_0}{B_0B}\right),$$
(1)
where $`a_0`$ is the off-resonant scattering length and $`\mathrm{\Delta }_0`$ and $`B_0`$ characterize the width and position of the resonance, respectively. On the side of the resonance where $`a`$ increases towards $`+\mathrm{}`$, the spin-singlet molecule becomes a weakly-bound $`s`$-wave state.

In this Letter, we use effective field theory methods to calculate the recombination constant $`\alpha `$ for a weakly-bound $`s`$-wave state. We find that the naive scaling relation $`\alpha =𝒞\mathrm{}a^4/m`$ is modified by renormalization effects, so the coefficient $`𝒞`$ is not universal, and that its maximum value is $`𝒞_{\mathrm{max}}=67.9`$. Our results confirm those of Refs. , demonstrate that they are model-independent, and provide a new tool for precise studies of the recombination rate. On the side of a Feshbach resonance where $`a\mathrm{}`$, the $`B`$-dependence of $`\alpha `$ can be predicted in terms of one adjustable parameter. It has a remarkable behavior, with $`𝒞`$ oscillating between zero and 67.9 more and more rapidly as $`B`$ approaches the resonance.

Effective field theory (EFT) is a powerful method for describing systems composed of particles with wave number $`k`$ much smaller than the inverse of the characteristic range $`R`$ of their interaction. (For a van der Waals potential $`C_6/r^6`$, $`R=(2mC_6/\mathrm{}^2)^{1/4}`$). EFT focuses on the aspects of the problem that are universal, independent of the details of short-distance interactions, by modelling the interactions as pointlike. The separation of scales $`k1/R`$ allows a systematic expansion in powers of the small parameter $`kR`$ . Generically, the scattering length $`a`$ is comparable to $`R`$, and the expansion is effectively in powers of $`ka`$. The pointlike interactions of the EFT generate ultraviolet divergences, but they can be absorbed into the renormalized coupling constants of the effective Lagrangian. All information about the influence of short-distance physics on low-energy observables is captured by these constants. At any given order in $`kR`$, only a finite number of coupling constants enter and this gives the EFT its predictive power. The domain of validity of EFT is $`kR1`$, even in the case of large scattering length $`aR`$. Thus it should accurately describe weakly-bound states with size of order $`a`$. However, the dependence on $`ka`$ is nonperturbative for $`k1/a`$, and it is necessary to reorganize the perturbative expansion into a new expansion in $`kR`$ by resumming higher order terms to all orders in $`ka`$. There has been significant progress recently in carrying out this resummation for the 3-body system. At leading order in $`kR`$, a single 3-body parameter is necessary and sufficient to carry out the renormalization . The scattering length and this 3-body parameter are sufficient to describe all low-energy 3-body observables up to errors of order $`R/a`$. In nuclear physics, this result has allowed a successful description of low-energy neutron-deuteron scattering and the binding energy of the triton. The variation of the 3-body parameter provides a natural explanation for the Phillips line .

We will apply this EFT for systems with large scattering length to the 3-body recombination problem. For simplicity, we now set $`\mathrm{}=1`$. We start by writing down a general local Lagrangian for a non-relativistic boson field $`\psi `$:

$$=\psi ^{}\left(i\frac{}{t}+\frac{\stackrel{}{}^2}{2m}\right)\psi \frac{C_0}{2}(\psi ^{}\psi )^2\frac{D_0}{6}(\psi ^{}\psi )^3+\mathrm{}.$$
(2)
The dots denote terms with more derivatives and/or fields; those with more fields will not contribute to the 3-body amplitude, while those with more derivatives are suppressed at low momentum. In order to set up integral equations for 3-body amplitudes, it is convenient to rewrite $``$ by introducing a dummy field $`d`$ with the quantum numbers of two bosons,

$``$ $`=`$ $`\psi ^{}\left(i{\displaystyle \frac{}{t}}+{\displaystyle \frac{\stackrel{}{}^2}{2m}}\right)\psi +d^{}d{\displaystyle \frac{g}{\sqrt{2}}}(d^{}\psi \psi +\text{h.c.})`$ (3)
$`+`$ $`hd^{}d\psi ^{}\psi +\mathrm{}.`$ (4)
The original Lagrangian (2) is easily recovered by a Gaussian path integration over the $`d`$ field, which implies $`d=(g/\sqrt{2})\psi ^2/(1+h\psi ^{}\psi )`$, $`C_0=g^2`$, and $`D_0=3hg^2`$. The atom propagator has the usual non-relativistic form $`i/(\omega p^2/2m)`$. The bare propagator for $`d`$ is simply $`i`$, but the exact propagator, including atom loops to all orders, is

$$iS_d(\omega ,\stackrel{}{p})=\frac{i4\pi /(mg^2)}{1/a+\sqrt{m\omega +\stackrel{}{p}^{\mathrm{\hspace{0.17em}2}}/4iϵ}},$$
(5)
where $`a`$ is the scattering length, which is related to the bare parameter $`g`$ and the ultraviolet (UV) cutoff $`\mathrm{\Lambda }`$ by

$$a=\frac{mg^2}{4\pi }\left(1+\frac{mg^2\mathrm{\Lambda }}{2\pi ^2}\right)^1.$$
(6)
The propagator (5) has a pole at $`\omega =1/(ma^2)+\stackrel{}{p}^{\mathrm{\hspace{0.17em}2}}/(4m)`$ corresponding to a weakly-bound state. Attaching four atom lines to this propagator gives the exact two-atom scattering amplitude. Taking the incoming atoms to have momenta $`\pm \stackrel{}{p}`$, the amplitude is $`(1/aip)^1`$, confirming the identification of $`a`$ as the scattering length.

We now consider the 3-body recombination process. We take the momenta of the incoming atoms to be small compared to the momenta of the final particles, which have magnitude $`p_f`$. Using Fermi’s golden rule, the recombination coefficient can be written

$$\alpha =\frac{map_f^2}{6\sqrt{3}\pi }\left|T(p_f)\right|^2,$$
(7)
where $`T(p)`$ is the amplitude for the transition between three atoms at rest and a final state consisting of an atom and a dimer in an $`s`$-wave state with momentum $`p`$ in their center-of-momentum frame. In Eq. (7), this amplitude is evaluated on shell at the value $`p_f=2/(\sqrt{3}a)`$ prescribed by energy conservation. However, $`T(p)`$ is also defined at off-shell values of $`p`$. The first few diagrams contributing to $`T`$ are illustrated in Fig. 1. All loop diagrams are of the same order as the tree diagrams, and they therefore have to be summed to all orders. This is conveniently accomplished by solving the integral equation represented by the second equality in Fig. 1. Corrections to this equation are of order $`R/a`$. The integral equation is

$`T(p)={\displaystyle \frac{96\pi ^{3/2}\sqrt{a}}{m}}\left({\displaystyle \frac{1}{p^2}}+{\displaystyle \frac{h}{2mg^2}}\right)+{\displaystyle \frac{2}{\pi }}{\displaystyle _0^\mathrm{\Lambda }}𝑑q`$ (8)
$`\times {\displaystyle \frac{q^2T(q)}{1/a+\sqrt{3}q/2iϵ}}\left[{\displaystyle \frac{1}{pq}}\mathrm{ln}\left|{\displaystyle \frac{q^2+pq+p^2}{q^2qp+p^2}}\right|+{\displaystyle \frac{h}{mg^2}}\right],`$ (9)
where we have inserted an UV cutoff $`\mathrm{\Lambda }`$ on the integral over $`q`$. If we were allowed to set $`h=0`$ and take $`\mathrm{\Lambda }\mathrm{}`$, a rescaling of the variables in Eq. (8) would lead to $`T(p)=K(ap)a^{3/2}/(mp)`$, with $`K(x)`$ a dimensionless function. Evaluating this in Eq. (7), the scaling relation $`\alpha =𝒞a^4/m`$ would follow immediately. However, the integral equation (8) has the same properties as the one describing atom-dimer scattering  and the limit $`\mathrm{\Lambda }\mathrm{}`$ can not be taken. The individual diagrams are finite as $`\mathrm{\Lambda }\mathrm{}`$, but their sum is sensitive to the cutoff. In an EFT, the dependence on the UV cutoff is cancelled by local counterterms. In Ref. , it was shown that the dependence of the low-energy observables on the cutoff $`\mathrm{\Lambda }`$ could be precisely compensated by varying $`h`$ appropriately. Writing $`h=2mg^2H(\mathrm{\Lambda })/\mathrm{\Lambda }^2`$, it was found that $`H(\mathrm{\Lambda })`$ could be well approximated by

$$H(\mathrm{\Lambda })\mathrm{tan}\left[s_0\mathrm{ln}(\mathrm{\Lambda }/\mathrm{\Lambda }_{})\pi /4\right],$$
(10)
where $`s_01.0064`$ is determined by the asymptotic behavior of the integral equation. This expression defines a parameter $`\mathrm{\Lambda }_{}`$  that characterizes the effect of the 3-body force on physical observables. A remarkable feature of this expression is its periodicity in $`\mathrm{ln}\mathrm{\Lambda }`$. As $`\mathrm{\Lambda }`$ is increased, $`H(\mathrm{\Lambda })`$ decreases to $`\mathrm{}`$, changes discontinously to $`+\mathrm{}`$, and continues decreasing.

The scaling violations in $`T(p)`$ from the renormalization of the 3-body force can be expressed as a dependence on $`a\mathrm{\Lambda }_{}`$. The simple scaling relation for $`\alpha `$ is therefore replaced by $`\alpha =𝒞(\mathrm{\Lambda }_{}a)a^4/m`$. Thus the value of $`𝒞`$ is not universal. In Fig. 2, we show $`\alpha `$ as a function of $`a`$ for $`\mathrm{\Lambda }_{}a_0=1.78,\mathrm{\hspace{0.17em}4.15},\mathrm{\hspace{0.17em}7.26}`$ and 19.77, where $`a_0`$ is an arbitrary but fixed length scale. Interestingly, $`\alpha `$ appears to oscillate as a function of $`\mathrm{ln}a`$ between zero and a maximum value $`𝒞_{\mathrm{max}}`$. We find that the curves can be very well fit by the expression

$$\alpha \frac{\mathrm{}a^4}{m}𝒞_{\mathrm{max}}\mathrm{cos}^2\left[s_0\mathrm{ln}(a\mathrm{\Lambda }_{})+\delta \right],$$
(11)
with $`𝒞_{\mathrm{max}}=67.9\pm 0.7`$ and $`\delta =1.74\pm 0.02`$.

We now compare our result to those of Refs. . Using an approximate solution to the 3-body wavefunction in hyperspherical coordinates, Fedichev et al.  found that the coefficient $`𝒞`$ has the universal value $`𝒞=3.9`$, independent of the interaction potential. We find that the value of $`𝒞`$ is not universal, but can vary from zero to about 67.9, depending on the details of 3-body interactions at short distances. We suggest that the specific value in Ref.  must correspond to an implicit assumption about the short-distance behavior of the 3-body wavefunction. Our results are consistent with those of Refs. , which were obtained using the hyperspherical adiabatic approximation. Nielsen and Macek  obtained $`𝒞_{\mathrm{max}}68.4`$ by applying the hidden crossing theory. Esry et al.  used coupled channel calculations to obtain $`\alpha `$ numerically for over 120 different 2-body potentials. For $`aR`$, their empirical result has the form of Eq. (11) with $`s_0=1`$ and $`𝒞_{\mathrm{max}}=60\pm 13`$. Refs.  show that the zeroes in $`\alpha `$ arise from interference effects involving 2-body and 3-body hyperspherical adiabatic potentials. The origin of these interference effects is less obvious in our EFT approach. However, our approach has several other advantages. First, it is completely model independent. Second, it is a controlled approximation, with corrections from finite range effects suppressed by powers of $`R/a`$. Third, it has predictive power in that all other low-energy 3-body observables can be determined in terms of $`a`$ and $`\mathrm{\Lambda }_{}`$. For example, the atom-dimer scattering length $`a_d`$ can be fit by

$$a_da\left(1.41.8\mathrm{tan}[s_0\mathrm{ln}(a\mathrm{\Lambda }_{})+3.2]\right).$$
(12)
We now apply the EFT to Feshbach resonances, where the value of $`a`$ is varied by changing the external magnetic field (cf. Eq. (1)). Our formalism is only applicable close to the Feshbach resonance on the side where $`a>0`$, because only in that region is there a weakly-bound molecule with $`B_d1/(ma^2)`$. Away from the resonance, or close to the resonance but on the side where $`a<0`$, 3-body recombination must involve more deeply-bound molecules with binding energies of order $`1/(mR^2)`$. To predict $`\alpha `$ as a function of the magnetic field $`B`$, we must specify how the parameter $`\mathrm{\Lambda }_{}`$ in Eq. (10) varies as a function of $`B`$. A Feshbach resonance is characterized by nonanalytic dependence of the scattering length $`a`$ on $`B`$. In a field theory, nonanalytic dependence on external parameters arises from long-distance fluctuations . The explicit ultraviolet cutoff $`\mathrm{\Lambda }`$ in our EFT excludes long-distance effects, which could introduce nonanalytic dependence on $`B`$, from the coefficients in $``$. Thus the bare parameters $`C_0=g^2`$ and $`D_0=3hg^2`$ in Eq. (2) should be smooth functions of $`B`$ for a fixed value of $`\mathrm{\Lambda }`$. The resonant behavior of the scattering length near $`B=B_0`$ in Eq. (1) can be reproduced by approximating $`g^2`$ by a linear function of $`B`$ in the resonance region. The parameter $`h`$ should also be a smooth function of $`B`$, but we can take $`h`$ to be approximately constant over the narrow resonance region. This assumption implies via Eq. (10) that $`\mathrm{\Lambda }_{}`$ should be constant while $`a(B)`$ varies as in Eq. (1).

The behavior of the recombination coefficient near the Feshbach resonance can therefore be read off from Fig. 2. If $`\alpha `$ is measured at one value of $`B`$ for which $`a(B)|a_0|`$, it determines $`\alpha `$ as a function of $`B`$ up to a two-fold ambiguity corresponding to whether the slope of $`𝒞`$ is positive or negative at that value of $`B`$. As $`B`$ approaches $`B_0`$, $`𝒞`$ should oscillate between 0 and $`𝒞_{\mathrm{max}}`$ in the manner shown in Fig. 2. The successive zeros correspond to values of $`a(B)`$ that differ roughly by multiplicative factors of $`\mathrm{exp}(\pi /s_0)22.7`$. Thus EFT makes the remarkable prediction that there are values of the magnetic field close to a Feshbach resonance where the contribution to $`\alpha `$ from a weakly-bound state vanishes.

The loss rate of 23Na atoms from a Bose-Einstein condensate near a Feshbach resonance has been studied by Stenger et al. . Our theory applies to the low-field side of the resonance at $`B=907`$ G. Taking into account the 3 atoms lost per recombination event and a Bose factor of $`1/3!`$, the loss rate from the condensate is $`\dot{N}=\alpha Nn^2/2`$. The loss rates measured in Ref.  correspond to a coefficient $`𝒞300`$ both off and near the resonance. This value is a factor of 4 larger than our maximum value. If $`𝒞>𝒞_{\mathrm{max}}`$, the 3-body recombination rate must be dominated not by the weakly-bound Feshbach resonance but instead by molecules with much larger binding energies $`1/(mR^2)`$. Alternatively, the large loss rate in Ref.  could be due to collective effects associated with the Bose-Einstein condensate, such as a 2-body recombination process involving atoms and dimers in a molecular condensate . In Fig. 3, we show the contribution to $`\alpha `$ from the weakly-bound state as a function of the magnetic field $`B`$ for $`\mathrm{\Lambda }_{}=1.94\text{ nm}^1`$, $`1.16\text{ nm}^1`$, and $`0.65\text{ nm}^1`$. If this contribution to $`\alpha `$ could be isolated, the first zero may be wide enough to be observed by experiment. The higher zeroes, however, are increasingly narrow and very close to the resonance.

We have applied an EFT for atoms with large scattering length to the problem of 3-body recombination into a weakly-bound $`s`$-wave state. We find that the coefficient $`𝒞`$ in the scaling relation $`\alpha =𝒞\mathrm{}a^4/m`$ is not universal, but must be in the range $`0𝒞67.9`$. If the 3-body recombination rate is measured to be larger than the maximum value, there must be a large contribution from molecules that are more deeply bound. Other low-energy 3-body observables, such as the atom-dimer scattering length, can be predicted in terms of $`a`$ and $`𝒞`$. Near a Feshbach resonance as $`a\mathrm{}`$, we find that $`𝒞`$ should oscillate between zero and 67.9. In any interval of $`B`$ over which $`a`$ increases by a factor of $`\mathrm{exp}(\pi /s_0)22.7`$, $`\alpha `$ should have a zero. Assuming that it is dominated by recombination into the weakly-bound state, the 3-body loss rate should have a minimum at that value of $`B`$. If a Bose-Einstein condensate was prepared at such a value of the magnetic field, one could study its behavior with large scattering length and relatively small 3-body losses.

We thank C.H. Greene, G.P. Lepage, J.H. Macek, and U. van Kolck for useful discussions. This research was supported in part by NSF grant PHY-9800964 and by DOE grants DE-FG02-91-ER4069 and DOE-ER-40561.


format: text

text:
# Regularities in football goal distributions

## Abstract

Besides of complexities concerning to football championships, it is identified some regularities in them. These regularities refer to goal distributions by goal-players and by games. In particular, the goal distribution by goal-players it well adjusted by the Zipf-Mandelbrot law, suggesting a conection with an anomalous decay.

Regularity in some complex systems can sometimes be identified and expressed in terms of simple laws. Typical examples of such situations are found in a wide range of contexts as the frequency of words in a long text, the population distribution in big cities, forest fires, the distribution of species lifetimes for North American breeding bird populations, scientific citations, www surfing, ecology, solar flares, economic index, epidemics in isolated populations, among others. Here, universal behaviours in the most popular sport, the football, are discussed. More precisely, this work focuses on regularities in goal distribution by goal-players and by games in championships. Furthermore, the goal distribution by goal-players is connected with an anomalous decay related to the Zipf-Mandelbrot law and with Tsallis nonextensive statistical mechanics.

In many contexts, it is common that few phenomena with high intensity arise, and so do many ones with low intensity. For instance, a long text generally contains many words that are employed in few opportunities and a small number that occurs largely. The above mentioned systems are good examples too. In particular, this kind of behaviour usually occurs in football championships, because there are many players that make few goals in contrast with the topscorers.

A detailed visualization of this behaviour can be well illustrated by considering some of the most competitive and traditional championships of the world. Our particular choice of championships has been done guided by the criterion of easy accessibility of the corresponding data to anyone. Therefore, we consider, here, some of the main league football championships from Italy, England, Spain and Brazil. Each of these championships has the participation of about twenty teams, contains around three hundred games, and approximately eight hundred goals. In Fig. 1 we exhibit data of these championships. In these graphics, the abscissa presents the number of goals $`x`$ divided by an average of goals $`m`$ (total number of goals per total number of goal-players), and the ordinate indicates the quantity $`N(x)`$ of players with $`x`$ goals divided by the quantity of players with one goal, $`N(1)`$. The regular shape of the graphics presented in Fig. 1 suggests a general law to describe the distribution of goals.

In the study of the majority of the previously cited systems, the Zipf’s law, $`N(x)=a/x^b`$, arises naturally, at least in part of the analysis. In the Zipf’s law, $`a`$ and $`b`$ are constants and $`x`$ is the independent variable. In order to give a better adjustment to a large part of the data, and based on information theory, Mandelbrot proposed $`N(x)=a/(c+x)^b`$ as a generalization of the Zipf’s law, with $`a`$, $`b`$, and $`c`$ all being constants. This Zipf-Mandelbrot’s distribution also arises in the context of a generalized statistical mechanics proposed some years ago, equivalently rewritten as

$$N(x)=N_0[1(1q)\lambda x]^{\frac{1}{1q}},$$
(1)
where $`N_0`$, $`\lambda `$, a $`q`$ are real parameters. In addition, this function satisfies an anomalous decay equation,

$$\frac{d}{dx}\left(\frac{N(x)}{N_0}\right)=\lambda \left(\frac{N(x)}{N_0}\right)^q.$$
(2)
The parameter $`q`$ can be considered as a measure of how anomalous the decay is. In particular, equation (1) is reduced to the usual exponential decay, $`N(x)=N_0\mathrm{exp}(\lambda x)`$, in the limit $`q1`$.

Motivated by these physical connections, we employ the distribution (1) to adjust the goals data. Following the construction of Fig. 1, we use the number of goalplayers with one goal, $`N(1)`$, and the average goal number by goalplayer, $`m`$, to eliminate $`N_0`$ and $`\lambda `$. Furthermore, it is a good approximation to replace the discrete average with a continuous one in the present analysis, i.e.,

$$m=\frac{_0^{\mathrm{}}xN(x)}{_0^{\mathrm{}}N(x)}=\frac{1}{\lambda (32q)}(q<3/2).$$
(3)
Thus, the distribution of goals dictated by equation (1) can be rewritten as

$$N(x)=N(1)\frac{\left[1\frac{(1q)}{(32q)m}\right]^{\frac{1}{q1}}}{\left[1\frac{(1q)}{(32q)m}x\right]^{\frac{1}{q1}}},$$
(4)
where $`q`$ becomes the unique parameter that remains to be adjusted, since $`N(1)`$ and $`m`$ are obtained directly from the data. Fig. 2 illustrates applications of equation (4) for four championships, indicating therefore the goodness of the formula (4). The same conclusion is obtained in the other championships showed in Fig. 1. Here, $`q=1.33`$ was employed as approximated value, leading to the Zipf-Mandelbrot’s exponent $`b3`$. In this way, $`q1.33`$ can be interpreted as the universal parameter for this kind of championships. Also, it is interesting to remark that $`b3`$ occurs in the distribution of scientific citations .

Other kinds of regularities in goal distributions can be identified, but with different behaviours. This fact can be verified in the distribution of goals per game. Proceeding in a similar way as done in Fig. 1, it is considered normalized scale distribution of games and goals. In this case, the abscissa is the number $`x`$ of goal divided by $`M`$, the mean goal per game of a championship (the number of goals of a championship divided by the corresponding number of games). In addition, the ordinate is given by the number of games with $`x`$ goals of a championship divided by the number of games of the corresponding championship. Fig. 3 contains this kind of graphics illustrating this regular behaviour by considering, again, the main league football championships from Italy, England, Spain and Brazil. As one can see, this figure strongly suggests a regularity in the distribution of goals for distinct championships.

Besides the numberless factors, including fluctuation due to relatively small number of teams and games that are present in a football championship, regularities arise in the goal distributions. In particular, the goal distribution for players that make goals are well adjusted by a Zipf-Mandelbrot law, suggesting a connection with ubiquitous phenomena such as anomalous diffusion.

###### Acknowledgements.

One of us, Mendes, R. S., thanks partial finantial support by CNPq (Brazilian Agency).


text:
# The ALADIN Interactive Sky Atlas

## 1 Introduction

### 1.1 The CDS

The Centre de Données astronomiques de Strasbourg (CDS) defines, develops, and maintains services to help astronomers find the information they need from the very rapidly increasing wealth of astronomical information, particularly on-line information.

In modern astronomy, cross-matching data acquired at different wavelengths is often the key to the understanding of astronomical phenomena, which means that astronomers have to use data and information produced in fields in which they are not specialists. The development of tools for cross-identification of objects is of particular importance in this context of multi-wavelength astronomy.

A detailed description of the CDS on-line services can be found, e.g., in Egret et al. (cds-amp2 (1995)) and in Genova et al. (cds-hub (1996), cds (1998), cds2000 (2000)), or at the CDS web site111*Internet address:* http://cdsweb.u-strasbg.fr/.

### 1.2 The ALADIN Project

Several sites currently provide on-line access to digitized sky surveys at different wavelengths: this is, for instance, the case of Digitized Sky Survey (DSS) at STScI (Morrison 1995adass…4..179M (1995)), and of similar implementations at other sites, providing quick access to cutouts of the compressed DSS images. SkyView at HEASARC (McGlynn et al. skyview (1997)) can generate images of any portion of the sky at wavelengths in all regimes from radio to gamma-ray. Some of these services provide simultaneous access to images and to catalogue data. The SkyCat tool, recently developed at ESO (Albrecht et al. skycat (1997)), addressed this concern in the context of the European Southern Observatory scientific environment (in view of supporting future users of the Very Large Telescope); SkyCat uses a standardized syntax to access heterogeneous astronomical data sources on the network.

Aladin has been developed independently by the CDS since 1993 as a dedicated tool for identification of astronomical sources – a tool that can fully benefit from the whole environment of CDS databases and services, and that is designed in view of being released as a multi-purpose service to the general astronomical community.

Aladin is an interactive sky atlas, allowing the user to visualize a part of the sky, extracted from a database of images from digitized surveys and observational archives, and to overlay objects from the CDS catalogues and tables, and from reference databases (Simbad and NED), upon the digitized image of the sky.

It is intended to become a major cross-identification tool, since it allows recognition of astronomical sources on the images at optical wavelength, and at other wavelengths through the catalogue data. Expected usage scenarios include multi-spectral approaches such as searching for counterparts of sources detected at various wavelengths, and applications related to careful identification of astronomical objects. Aladin is also heavily used for the CDS needs of catalogue and database quality control.

In the case of extensive undertakings (such as checking the astrometric quality for a whole catalogue), it is expected that Aladin will be useful for understanding the characteristics of the catalogue or survey, and for setting up the parameters to be adjusted while fine tuning the cross-matching or classification algorithms, by studying a sample section of objects or fields.

A discussion of the usage of such a tool for cross-identification can be found in Bartlett & Egret (xid-179 (1997)), where it is shown how *training sets* are used to build likelihood ratio tests.

The Aladin interactive atlas is available in three modes: a simple previewer, a Java interface, and an X-Window interface. We describe here mostly the Java interface which is publicly accessible on the World-Wide Web.

## 2 Access modes

After a long phase of development, (see e.g., Paillou et al. paillou (1994)), Aladin has been first distributed to a limited number of astronomy laboratories in 1997, as an X-Window client program, to be installed on a Unix machine on the user side. The client program interacts with the servers running on Unix workstations at CDS (image server, catalogue server, Simbad server) and manages image handling and plane overlays.

The strategy of having a client program on the user side is difficult to maintain on the long run. The World-Wide Web offers, with the development of Java applications (or *applets*), a way to solve this difficulty. Actually, there is still a *client* program: this is the Java applet itself, that the user receives from the WWW server. Most current Internet browsers are able to make it run properly, so that the user does not have to install anything special other than an Internet browser.

As a consequence, Aladin is currently available in the three following modes:

a pre-formatted image server provides a compressed image of fixed size ($`14.1\mathrm{}\times 14.1\mathrm{}`$ for the DSS-I) around a given object or position. When an object name is given, its position is resolved through the Simbad name resolver. Anchors pointing to the previewer are integral part of the World-Wide Web interfaces to the Simbad database222*Internet address*: http://simbad.u-strasbg.fr/Simbad and to the CDS bibliographic service333*Internet address*: http://simbad.u-strasbg.fr/biblio.html. The result page also gives access to the full resolution FITS image for download.

this is the primary public interface, supporting queries to the image database and overlays from any catalogue or table available at CDS, as well as from Simbad and NED databases. Access to personal files is not possible (due to security restrictions of the Java language). These restrictions do not apply to the *stand-alone* version, which can be installed and run on a local *Java virtual machine*.

The X-Window Aladin client provides most of the functionalities of the Aladin Java interface, plus more advanced functions, as described below (section 6).

## 3 The image database

### 3.1 Database summary

The Aladin image dataset consists of:

* The whole sky image database from the first Digitized Sky Survey (DSS-I) digitized from photographic plates and distributed by the Space Telescope Science Institute (STScI) as a set of slightly compressed FITS images (with a resolution of $`1.8\mathrm{}`$); DSS-II is also currently being integrated into the database (see below);
* Images of *crowded* fields (Galactic Plane, Magellanic Clouds) at the full resolution of $`0.67\mathrm{}`$, scanned at the Centre d’Analyse des Images (MAMA machine) in Paris;
* Global plate views ($`5\mathrm{°}\times 5\mathrm{°}`$ or $`6\mathrm{°}\times 6\mathrm{°}`$ according to the survey) are also available for all the plates contributing to the image dataset: these are built at CDS by averaging blocks of pixels from the original scans;
* Other images sets, or user-provided images, in FITS format, having suitable World Coordinate System information in the header (see e.g. Greisen & Calabretta wcs (1995)); this functionality is currently available only for the Java stand-alone version.

### 3.2 Building the database contents

The Aladin project has set up collaborations with the major groups providing digitizations of sky surveys. The original surveys are made of photographic Schmidt plates obtained at Palomar in the North, and ESO/SERC in the South, and covering the whole sky at different epochs and colours (see e.g., MacGillivray potsdam (1994)).

The database currently includes the first Digitized Sky Survey (DSS-I) produced by the Space Telescope Institute (Lasker dss (1992)), for the needs of the Hubble Space Telescope. To create these images, the STScI team scanned the first epoch (1950/1955) Palomar $`E`$ Red and United Kingdom Schmidt $`J`$ Blue plates (including the SERC J Equatorial Extension and some short V-band plates at low galactic latitude) with a pixel size of $`1.7\mathrm{}`$ ($`25\mu m`$). The low resolution and a light data compression (factor of 10) permit storage of images covering the full sky on a set of 102 CD-ROMs.

DSS-II images in the R-band (from Palomar POSS-II F and UK Schmidt SES, AAO-R, and SERC-ER), scanned with a $`1\mathrm{}`$ ($`15\mu m`$) sampling interval (see Lasker dss-ii (1994)) are gradually being included into the system, and will soon be followed by DSS-II images in the B-band (POSS-II J).

In addition, high resolution digitalization of POSS-II, SERC-J, SERC-SR, SERC-I, or ESO-R plates featuring crowded regions of the sky (Galactic Plane and Magellanic Clouds) have been provided by the MAMA facility at the Centre d’Analyse des Images (CAI), Observatoire de Paris (Guibert MAMA (1992)). Sampling is $`0.67\mathrm{}`$ per pixel ($`10\mu m`$). Currently, these high resolution images cover about 15% of the sky, and are stored in a juke-box of optical disks, with a capacity of 500 Gigabytes.

### 3.3 The image server

The image server for Aladin had to be able to deal with various survey data, in heterogeneous formats (uncompressed FITS, compressed JPEG or PMT – see Section 5, etc.). For that, an object-oriented design was chosen, allowing an easy manipulation of image calibrations and headers, through the use of object classes. Image compression or decompression, image reconstruction, and in a near future, part of the recalibration, are seen as class methods.

Images are currently divided into subimages of $`500\times 500`$ pixels (DSS-I), $`768\times 768`$ pixels (DSS-II), or $`1024\times 1024`$ pixels (MAMA).

The 1.5 million subimages are described by records stored in a relational database, encapsulated by several classes of the image management software. When an image of the sky is requested, the original subimages containing the corresponding sky area are retrieved through SQL commands, and the resulting image is built on the fly.

## 4 Usage scenarios

In this section we will focus on describing the usage of the Aladin Java interface, as it is available now (November 1999).

### 4.1 Access

The Aladin home page is available through the CDS Web server at the following address: http://aladin.u-strasbg.fr/

This site provides access to Aladin documentation, including scientific reports, recent publications, etc.

### 4.2 Query modes

The typical usage scenario starts with a request of the digitized image for an area of the sky defined by its central position or name of central object (to be resolved by Simbad). The size of the sky field is determined by the photographic survey used: it is $`14.1\mathrm{}`$ in the case of the DSS-I.

Astrometric information comes from the FITS header of the DSS image, and is generally accurate to the arcsecond (with deviations up to several arcsec. in exceptional cases, on plate edges).

In a subsequent step, the interface, illustrated by Figs. 1 and 2, allows the user to stack several information planes related to the same sky field, to superimpose the corresponding data from catalogues and databases, and to obtain interactive access to the original data.

The possible information planes are the following:

* Image pixels from the Aladin database of digitized photographic plates (DSS-I, MAMA, DSS-II); functionalities include zooming capabilities, inverse video, modification of the color table;
* Information from the Simbad database (Wenger et al. simbad (2000)); objects referenced in Simbad are visualized by color symbols overlaid on top of the black and white image; the shape and color of the symbols can be modified on request, and written labels can be added for explicit identification of the objects; these features are also available for all the other information planes;
* Records from the CDS library of catalogues or tables (VizieR444*Internet address:* http://vizier.u-strasbg.fr/, Ochsenbein et al. vizier (2000)); the user can select the desired catalogue from a preselected list including the major reference catalogues such as the Tycho Reference Catalogue (ESA tyc (1997); Høg et al. trc (1998)), GSC (Lasker et al. gsc (1990)), IRAS Point Source Catalog, or USNO A2.0 (Monet usno (1998)); the user can alternatively select the catalogues for which entries may be available in the corresponding sky field, using the VizieR query mechanism by position (see 4.3), catalogue name or keyword;
* Information from the NED database: objects referenced in the NASA/IPAC Extragalactic Database555*Internet address:* http://nedwww.ipac.caltech.edu/ (Helou et al. ned (2000)) can also be visualized through queries submitted to the NED server at IPAC;
* Archive images will gradually become available through the corresponding mission logs: Hubble Space Telescope images are currently available (see Fig. 3 for an example), and more archives will follow.
* Local, user data files can also be overlaid, but, because of current limitations of the Java applications, this feature is only available in the stand-alone version, or in Aladin X.

The stack of images and graphics is made visible to the user (under the eye icon, on the right of Fig. 2) so that each plane can be turned on and off. The status of queries is also easily visualized.

For all information planes (Simbad, VizieR, NED) links are provided to the original data. This is done in the following way: when selecting an object on the image, with mouse and cursor, it is possible to call for the corresponding information which will appear in a separate window on the Internet browser. It is also possible to select with the mouse and cursor all objects in a rectangular area: the corresponding objects are listed in a panel on the bottom of the Aladin window; this list includes basic information (name, position and, when applicable, number of bibliographical references) and anchors pointing to the original catalogue or database.

At any moment the position of the cursor is translated in terms of right ascension and declination on the sky and visualized in the top panel of the Aladin window. Additional features are available, such as a tool for computing angular distance between marked objects.

The *standalone* version gives access to additional facilities, not available through the Web, such as printing and saving the images and data.

### 4.3 The catalogue server

The ability to access all VizieR catalogues and tables directly from Aladin is a unique feature which makes it an extremely powerful tool for any cross-identification or classification work.

The “*Select around target*” request relies on a special feature – the genie of the lamp: this is the ability to decide which catalogues, among the database of (currently) over 2,600 catalogues or tables, contain data records for astronomical objects lying in the selected sky area. In order to do that, an index map of VizieR catalogues is produced (and kept up-to-date), on the basis of about ten pixels per square degree: for each such ‘pixel’ the index gives the list of all catalogues and tables which have entries in the field.

When a user hits the button “*Select around target*”, this index is queried and the list of useful catalogues is returned. It is possible, at this stage, either to list all catalogues, or to produce a subset selected on the basis of keywords. Note that, as the index “pixels” generally match an area larger than the current sky field, there is simply a good chance, but not 100%, to actually obtain entries in the field when querying one of the selected tables.

### 4.4 Cache

The images of the 30,000 most cited objects in Simbad are pre-computed and available on a cache on magnetic disk. For these objects, the image is served much faster than for other objects where the image has to be extracted from the Digitized Sky Survey.

### 4.5 Usage statistics

As the newest service developed by CDS, Aladin has not yet been widely publicized, and its usage is in a steeply growing phase. Currently about 10,000 queries are processed monthly, generating the extraction of more than 5,000 images.

## 5 Image compression

Astronomical image compression in the context of Aladin has been discussed in detail by Louys et al. (louys (1999)).

For the Aladin Java interface and for the Aladin previewer, the current choice has been to deliver to the user an image in JPEG 8-bit format, constructed from the original FITS images. JPEG is a general purpose standard which is supported by all current Internet browsers. The size of such an image does not exceed 30 kBytes, and thus the corresponding network load is very small.

In the near future, the Pyramidal Median Transform (PMT) algorithm, implemented in the MR-1 package (Starck et al. pmt (1996)), will be used within Aladin for storing or transferring new image datasets, such as additional high resolution images (see again Louys et al. louys (1999) for details). The corresponding decompression package is being written in Java code, and could be downloaded on request for use within the Java interface.

## 6 Aladin X

The Aladin X-Window interface is the testbed for further developments. It is currently only distributed for the Unix Solaris operating system. Interested potential users should contact CDS for details.

### 6.1 Source extraction

Aladin X includes a procedure for source extraction. The current mechanism will soon be replaced by SExtractor (Bertin & Arnouts 1996A&AS..117..393B (1996)).

### 6.2 Plate calibrations

While the first level astrometric calibrations are given by the digitizing machines, a second level is being developed that will allow the user to *recalibrate* the image with a new set of standards taken, for example, from the Tycho Reference Catalogue. The photometric calibrations (surface and stellar) will eventually also be performed within Aladin, by using the Guide Star Photometric Catalogs (GSPC I and II; Ferrari et al. gspc2 (1994); Lasker et al. gspc1 (1988)).

Users will thus be able to work on the details of local astrometric and photometric plate calibrations in order to extract the full information from the digitized plates.

## 7 Integration of distributed services

While the CDS databases have followed different development paths, the need to build a transparent access to the *whole set* of CDS services has become more and more obvious with the easy navigation permitted by hypertext tools. Aladin has become the prototype of such a development, by giving comprehensive simultaneous access to Simbad, the VizieR Catalogue service, and to external databases such as NED, using a client/server approach and, when possible, standardized query syntax and formats.

In order to be able to go further, the CDS has built a general data exchange model, taking into account all types of information available at the Data Center, known under the acronym of GLU for Générateur de Liens Uniformes – Uniform Link Generator (Fernique et al. glu (1998)).

More generally, with the development of the Internet, and with an increasing number of on-line astronomical services giving access to data or information, it has become critical to develop new tools providing access to distributed services. This is, for instance, the concern expressed by NASA through the AstroBrowse project (Heikkila et al. astrobrowse (1999)). A local implementation of this concept is available at CDS (AstroGlu: Egret et al. astroglu (1998)).

## 8 Future developments

An important direction of development in the near future is the possibility of providing access to images from other sky surveys or deep field observations: obvious candidates are the DENIS (Epchtein denis (1998)) and 2MASS (Skrutskie 2mass (1998)) near-infrared surveys. The first public point source catalogues resulting from these surveys are already available through Aladin, since they are included in the VizieR service. This has already proved useful for validating survey data in preliminary versions of the DENIS catalogue (Epchtein et al. denis-psc (1999)).

The CDS team will also continue to enrich the system functionality. The users play an important role in that respect, by giving feedback on the desired features and user-friendliness of the interfaces.

New developments are currently considered as additional modules which will be incorporated to the general release only when needed, possibly as optional downloads, in order to keep the default version simple and efficient enough for most of the Web applications.

On a longer term, the CDS is studying the possibility of designing *data mining* tools that will help to make a fruitful use of forthcoming very large surveys, and will be used for cross-matching several surveys obtained, for instance, at different wavelengths. A first prototype, resulting from a collaboration between ESO and CDS, in the framework of the VLT scientific environment is currently being implemented (Ortiz et al. ortiz (1999)).

###### Acknowledgements.

CDS acknowledges the support of INSU-CNRS, the Centre National d’Etudes Spatiales (CNES), and Université Louis Pasteur. We are indebted to Michel Crézé who initiated the project while being Director of the CDS, and to all the early contributors to the Aladin project: Philippe Paillou, Joseph Florsch, Houri Ziaeepour, Eric Divetain, Vincent Raclot. Collaboration with STScI, and especially with the late Barry Lasker, and with Brian McLean, is gratefully acknowledged. We thank Jean Guibert and René Chesnel from CAI/MAMA for their continuous support to the project. The Digitized Sky Survey was produced at the Space Telescope Science Institute under U.S. Government grant NAG W-2166. The images of these surveys are based on photographic data obtained using the Oschin Schmidt Telescope on Palomar Mountain and the UK Schmidt Telescope. Java is a registered trademark of Sun Microsystems.


text:
# Measuring the Nonlinear Biasing Function from a Galaxy Redshift Survey

## 1 INTRODUCTION

The fact that galaxies of different types cluster differently (e.g., Dressler 1980; Lahav, Nemiroff & Piran 1990; Santiago & Strauss 1992; Loveday et al. 1995; Hermit et al. 1996; Guzzo et al. 1997) indicates that the galaxy distribution is in general biased compared to the underlying mass distribution. Cosmological simulations confirm that halos and galaxies must be biased (e.g., Cen & Ostriker 1992; Kauffmann, Nusser & Steinmetz 1997; Blanton et al. 1999; Somerville et al. 2000). The biasing becomes even more pronounced at high redshift, as predicted by theory (e.g., Kaiser 1986; Davis et al. 1985; Bardeen et al. 1986; Dekel & Rees 1987; Mo & White 1996; Bagla 1998; Jing & Suto 1998; Wechsler et al. 1998), and confirmed by the strong clustering of galaxies observed at $`z3`$ (Steidel et al. 1996; 1998). Knowing the biasing scheme is crucial for extracting dynamical information and cosmological constants from the observed galaxy distribution, and may also be very useful for understanding the process and history of galaxy formation.

The simplest possible biasing model relating the density fluctuation fields of matter and galaxies, $`\delta `$ and $`\delta _\mathrm{g}`$, is the deterministic and linear relation, $`\delta _\mathrm{g}(xxx)=b\delta (xxx)`$, where $`b`$ is a constant linear biasing parameter. However, this is at best a crude approximation, because it is not self-consistent (e.g., it does not prevent $`\delta _\mathrm{g}`$ from becoming smaller than $`1`$ when $`b>1`$) and is not preserved in time. At any given time, scale and galaxy type, the biasing is expected in general to be nonlinear, i.e., $`b`$ should vary as a function of $`\delta `$. The nonlinearity of dark-matter halo biasing (as well as its dependence on scale, mass and time) is approximated fairly well by the model of Mo & White (1996), based on the extended Press-Schechter formalism (Bond et al. 1991). Improved approximations have been proposed by Jing (1998), Catelan et al. (1998), Sheth & Tormen (1999) and Porciani et al. (1999). It is quantified further for halos and galaxies using cosmological $`N`$-body simulations with semi-analytic galaxy formation (e.g., Somerville et al. 2000). The biasing is also expected, in general, to be stochastic, in the sense that a range of values of $`\delta _\mathrm{g}`$ is possible for any given value of $`\delta `$. For example, if the biasing is nonlinear on one scale, it should be different and non-deterministic on any other scale. The origin of the scatter is shot noise as well as the influence of physical quantities other than mass density (e.g., velocity dispersion, the dimensionality of the local deformation tensor which affects the shape of the collapsing object, etc.) on the efficiency of galaxy formation.

Dekel & Lahav (1999) have proposed a general formalism for galaxy biasing, that separates nonlinearity and stochasticity in a natural way. The density fields are treated as random fields, and the biasing is fully characterized by the conditional probability distribution function $`P(\delta _\mathrm{g}|\delta )`$. The constant linear biasing factor $`b`$ is replaced by a mean biasing function,

$$\delta _\mathrm{g}|\delta b(\delta )\delta ,$$
(1)
which can in principle take a wide range of functional forms, restricted by definition to have $`\delta _\mathrm{g}=0`$ and $`\delta _\mathrm{g}|\delta 1`$ for any $`\delta `$. The stochasticity is expressed by the higher moments about this mean, such as the conditional variance

$$\sigma _\mathrm{b}^2(\delta )ϵ^2|\delta /\sigma ^2,ϵ\delta _\mathrm{g}\delta _\mathrm{g}|\delta ,$$
(2)
scaled for convenience by the variance of mass fluctuations, $`\sigma ^2\delta ^2`$. To second order, the biasing function $`b(\delta )`$ can be characterized by two parameters: the moments $`\widehat{b}`$ and $`\stackrel{~}{b}`$,

$$\widehat{b}b(\delta )\delta ^2/\sigma ^2\mathrm{and}\stackrel{~}{b}^2b^2(\delta )\delta ^2/\sigma ^2.$$
(3)
The parameter $`\widehat{b}`$ is the natural extension of the linear biasing parameter, measuring the slope of the linear regression of $`\delta _\mathrm{g}`$ on $`\delta `$, and $`\stackrel{~}{b}/\widehat{b}`$ is a useful measure of non-linearity. The stochasticity is characterized independently by a third parameter, $`\sigma _\mathrm{b}^2ϵ^2/\sigma ^2`$. As has been partly explored by Dekel & Lahav (1999), these parameters should enter any nonlinear analysis aimed at extracting the cosmological density parameter $`\mathrm{\Omega }`$ from a galaxy redshift survey, and are therefore important to measure.

In this paper we propose a simple method to measure the biasing function $`b(\delta )`$ and the associated parameters $`\widehat{b}`$ and $`\stackrel{~}{b}`$ from observed data that are either already available, such as the PSC$`z`$ redshift survey (Saunders et al. 2000), or that will soon become available, such as the redshift surveys of 2dF (Colless 1999) and SDSS (e.g., Loveday et al. 1998) and high-redshift surveys such as DEEP (Davis & Faber 1999). Alternative methods have been proposed to measure the biasing function, using the cumulant correlators of the observed distribution of galaxies in redshift surveys (Szapudi 1998) or their bispectrum (Matarrese, Verde, Heavens 1997, Verde et al. 1998).

We first show in §2, using halos and galaxies in $`N`$-body simulations, that the difference between the cumulative distribution functions (CDFs) of galaxies and mass can be straightforwardly translated into $`\delta _\mathrm{g}|\delta `$ despite the scatter in the biasing scheme. Then, in §3, we demonstrate that for our purpose, $`C(\delta )`$ is insensitive to the cosmological model and can be approximated robustly by a cumulative log-normal distribution. This means that we do not need to observe $`C(\delta )`$, which is hard to do; we only need to measure $`C_\mathrm{g}(\delta _\mathrm{g})`$ and, independently, the rms value $`\sigma `$ of the mass fluctuations on the same scale. In §4, we slightly modify the method to account for redshift-space distortions, and use mock galaxy catalogs from N-body simulations to evaluate the associated errors. Finally, in §5, we estimate the errors due to the sparse sampling and finite volume. The method and its applications to existing and future data are discussed in §6.

## 2 BIASING FUNCTION FROM DISTRIBUTION FUNCTIONS

Let $`C_\mathrm{g}(\delta _\mathrm{g})`$ and $`C(\delta )`$ be the cumulative distribution functions of the density fluctuations of galaxies and mass respectively (at a given smoothing window). Had the biasing relation been deterministic and monotonic, it could have been determined straightforwardly from the difference between these CDFs at given percentiles,

$$\delta _\mathrm{g}(\delta )=C_\mathrm{g}^1[C(\delta )],$$
(4)
where $`C_\mathrm{g}^1`$ is the inverse function of $`C_\mathrm{g}`$. 111A similar relation has been used by Narayanan & Weinberg (1998) for “debiasing” the galaxy density field for the purpose of dynamical reconstruction. In the presence of scatter in the biasing scheme, strict monotonicity is violated, but it is possible that $`C_\mathrm{g}^1[C(\delta )]`$ is still a good approximation for $`\delta _\mathrm{g}|\delta `$, as long as the latter is monotonic.222The absence of spiral galaxies in the centers of rich clusters may result in a non-monotonic biasing function for this type of galaxies at small smoothing scales, as hinted in Blanton et al. (1999). However, using the simulations described in this section, Somerville et al. (2000) do not find non-monotonicity for late type galaxies at $`8h^1\mathrm{Mpc}`$ smoothing, as used in Figure 3 below. The validity of this approximation is addressed in the present section.

We use two cosmological $`N`$-body simulations in which both halos and galaxies were identified (Kauffmann et al. 1999). The cosmological models are $`\tau `$CDM (with $`\mathrm{\Omega }_\mathrm{m}=1`$ and $`h=0.5`$) and $`\mathrm{\Lambda }`$CDM (with $`\mathrm{\Omega }_\mathrm{m}=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$ and $`h=0.7`$). $`N=256^3`$ particles were simulated in a periodic box of comoving size $`85`$ and $`141h^1\mathrm{Mpc}`$ respectively (corresponding to a mass resolution of $`1.010^{10}h^1M_{}`$ and $`1.410^{10}h^1M_{}`$). The simulations were run using a parallel adaptive P3M code kindly made available by the Virgo Consortium (see Jenkins et al. 1998) as part of the “GIF” collaboration between the HU Jerusalem and the MPA Munich. The present epoch is defined by a linear rms density fluctuation in a top-hat sphere of radius $`8h^1\mathrm{Mpc}`$ of $`\sigma _8=0.6`$ in the $`\tau `$CDM simulation and $`\sigma _8=0.9`$ in the $`\mathrm{\Lambda }`$CDM simulation. Dark-matter halos were identified at densely sampled time steps using a friends-of-friends algorithm. Galaxies were identified inside these halos by applying in retrospect semi-analytic models (SAMs) of galaxy formation (Kauffmann et al. 1999). The SAMs simulate the important physical processes of galaxy formation such as gas cooling, star formation and supernovae feedback. At different times in the redshift range 0 to 3, we select halos by mass and galaxies by luminosity or type. We then compute density fields by applying top-hat smoothing with radii in the range $`515h^1\mathrm{Mpc}`$. We report detailed results for the case of $`8h^1\mathrm{Mpc}`$ smoothing, and refer to the scale dependence in several places.

The figures of this section illustrate the success of the approximation, equation (4), in several different cases based on the $`\tau `$CDM simulation, with top-hat smoothing of radius $`8h^1\mathrm{Mpc}`$ (hereafter TH8, or TH$`X`$ for radius $`Xh^1\mathrm{Mpc}`$), and at different redshifts. Figure 1 refers to halos of mass $`>10^{12}h^1M_{}`$ ($`>100`$ particles). On the top we show the cumulative distributions of halos and underlying mass fluctuations, $`C_\mathrm{g}(\delta _\mathrm{g})`$ and $`C(\delta )`$ (our notation does not distinguish between halos and galaxies). The errors in $`C_\mathrm{g}`$ are computed from 20 bootstrap simulations of the halo field. The errors in $`C`$, estimated in the same way, are smaller by an order of magnitude and are therefore not shown. The bottom panels show a point-by-point comparison of the TH8 fields of $`\delta _\mathrm{g}(xxx)`$ and $`\delta (xxx)`$ at points randomly chosen (1:8) from a uniform grid of spacing $`2.64h^1\mathrm{Mpc}`$ within the simulation box. The true mean biasing function $`\delta _\mathrm{g}|\delta `$ is marked by the filled circles with attached error bars. It is computed by a local linear regression of $`\delta _\mathrm{g}`$ on $`\delta `$ within each bin of $`\delta `$, adopting the value of the fitted line at the center of the bin (only every other bin is shown). Shown in comparison (solid line) is the approximation for $`\delta _\mathrm{g}|\delta `$ obtained by equation (4) from the CDFs, and the corresponding $`1\sigma `$ error range based on the bootstrap realizations (dotted lines).

As can be seen in Figure 1, the approximation is excellent over most of the $`\delta `$ range — the deviation at $`z`$=0 is within the $`1\sigma `$ errors up to $`\delta 1.4`$ (corresponding to $`97\%`$ of the volume). Systematic deviations show up at higher $`\delta `$ values, where the scatter becomes larger and the mean biasing function flatter, making the deviations from monotonicity larger. In order to quantify the quality of the approximation, we average the residuals (scaled by $`\sigma _\mathrm{g}`$):

$$\mathrm{\Delta }=\frac{1}{N_{\mathrm{bins}}\sigma _\mathrm{g}^2}\underset{\delta \mathrm{bins}}{\overset{N_{\mathrm{bins}}}{}}[\delta _\mathrm{g}(\delta )\delta _\mathrm{g}|\delta ]^2,$$
(5)
where $`\delta _\mathrm{g}(\delta )`$ is obtained via equation (4). We exclude the poorly recovered high-density tail by limiting the summation to those $`N_{\mathrm{bins}}`$ bins of $`\delta `$ for which $`C(\delta )<0.99`$ and $`C_\mathrm{g}(\delta _\mathrm{g})<0.99`$. The values of $`\mathrm{\Delta }`$ in the various cases studied, including halos and galaxies in $`\tau `$CDM and $`\mathrm{\Lambda }`$CDM at different redshifts, are listed in Table 1. For example, for the halos shown in Figure 1 at $`z=0`$ we obtain $`\mathrm{\Delta }=0.08`$, indicating that the typical error in the approximation $`\delta _\mathrm{g}(\delta )`$ is small compared to the actual scatter $`\sigma _\mathrm{g}`$ in the halo density field.

A complementary approach for quantifying the quality of the approximation is by testing how well it recovers the values of the moments of the biasing function, $`\widehat{b}`$ and $`\stackrel{~}{b}`$. In Table 1 we present the values of these moments for the different cases, as computed directly from the simulation and as approximated by $`\delta _\mathrm{g}(\delta )`$ (denoted by a subscript “a”). These biasing parameters are computed based on 99.9% of the volume, excluding the very highest density peaks, where the error is large (The only exception is at $`z`$=3, where we use only 99% of the volume because the errors are even larger). For the halos shown in Figure 1 at $`z=0`$, we see that $`\widehat{b}`$ and $`\stackrel{~}{b}`$ are recovered with errors of 1% and 3% respectively.

The middle panels of Figure 1 refer to $`z=1`$, where $`\widehat{b}2.2`$. The approximation of equation (4) holds well in this case up to $`\delta 0.7`$, which corresponds to $`98\%`$ of the volume. The approximation remains good despite the large scatter (compared to the $`z=0`$ case) because the steepness of the biasing function helps maintaining reasonable monotonicity. The goodness of the recovery of the biasing function, with $`\mathrm{\Delta }=0.07`$, is similar to the $`z=0`$ case. The parameters $`\widehat{b}`$ and $`\stackrel{~}{b}`$ are recovered with an accuracy of $`5\%`$ (Table 1). The right panels of Figure 1 demonstrate that the approximation is valid even at $`z`$=3, where the biasing is extremely strong, $`\widehat{b}6.6`$. The recovery of the biasing function is still good, $`\mathrm{\Delta }=0.20`$, and its moments are approximated to within $`2\%`$.

The halo biasing function in the $`\mathrm{\Lambda }`$CDM cosmology is recovered, in general, with similar success, as can be seen in the top part of Table 1. Note that in this case the recovery actually improves at higher redshift. This reflects the fact that in $`\mathrm{\Lambda }`$CDM the halo biasing scatter becomes smaller at higher redshift (see Somerville et al. 2000, Fig. 17). It results from the smaller shot noise due to the higher abundance of high-redshift halos in $`\mathrm{\Lambda }`$CDM compared to $`\tau `$CDM.

Figure 2 is analogous to Figure 1, but now for bright galaxies of $`M_\mathrm{B}5\mathrm{log}h<19.5`$. The recovery is again quantified in Table 1; it is quite similar to the case of halos. The typical error is $`\mathrm{\Delta }0.08`$, and the biasing parameters are recovered with an error of a couple to a few percent.

The performance of our method has been tested for smoothing scales in the range $`515h^1\mathrm{Mpc}`$. For the $`\tau `$CDM model, we find that the quality of the approximation is practically independent of scale throughout this range; the relative error in the biasing parameters is at the level of a few percent, and $`\mathrm{\Delta }`$ is in the range 0.1 to 0.2, rather similar to the values quoted in Table 1 for TH8 smoothing. On the other hand, for $`\mathrm{\Lambda }`$CDM we do find that the performance improves with increasing smoothing scale. With TH15 at $`z=0`$, for halos (or galaxies), the errors in the biasing parameters reduce to below 3% (1%), and $`\mathrm{\Delta }=0.07`$ (0.04), while for TH5 smoothing these errors are about 4 times larger. This difference between the two models can be attributed to a difference in the scale dependence of the biasing scatter (Somerville et al. 2000, Figure 16), which translates to an error in our method via increased deviations from monotonicity.

Before we proceed with the biasing relative to the underlying mass, we note that the relative  biasing function of two different galaxy types, $`\delta _{\mathrm{g}_2}|\delta _{\mathrm{g}_1}`$, can be directly observable from a redshift survey. Again, for a deterministic and monotonic biasing process one has

$$\delta _{\mathrm{g}_2}(\delta _{\mathrm{g}_1})=C_{\mathrm{g}_2}^1[C_{\mathrm{g}_1}(\delta _{\mathrm{g}_1})],$$
(6)
and when biasing scatter is present, the question is to what extent equation (6) provides a valid approximation for the true relative biasing function.

Figure 3 shows the relative biasing function of “early” and “late” type galaxies in the two cosmological models, at $`z=0`$ and with TH8 smoothing. These galaxy types are distinguished in the SAM $`N`$-body simulations according to the ratio of bulge to total luminosity in the V band being larger or smaller than 0.4 respectively. The large scatter in the relative biasing, due to errors in the two density fields, is reduced by including all the galaxies, without applying a luminosity cut.

As can be seen in the last three columns of Table 1, the quality of the recovery of the relative biasing function is not as good as in the case of the absolute biasing of galaxies or halos. The values of $`\mathrm{\Delta }`$ range from 0.2 to 0.56, compared to 0.08 to 0.16 in the former cases. This is expected, because in the case of relative biasing the two density fields contribute to the stochasticity or deviation from monotonicity (see also the important role of sampling errors in the recovery of the biasing function, §4.2). The moments of the relative biasing function are recovered with better than 15% accuracy at $`z1`$, and to $``$25% accuracy at $`z=3`$, in both cosmologies. In calculating the moments, unlike in Figure 3, a luminosity cut has been applied: $`M_\mathrm{B}5\mathrm{log}h<19.5`$, and 99% of the volume was used. The fact that the $`\mathrm{\Delta }`$ values are still significantly smaller than unity and the errors in the biasing parameters are not larger than 25% indicate that our method is capable of yielding meaningful estimates of the relative biasing function. In both cosmologies, the relative biasing is almost scale independent in the range 5 – 15 $`h^1\mathrm{Mpc}`$, as is the quality of the reconstruction.

## 3 THE MASS CDF: ROBUST AND LOGNORMAL

Large redshift surveys provide a rich body of data for mapping the galaxy density field in extended regions of space and computing its CDF with adequate accuracy. However, direct mapping of the mass density field is much harder. For example, POTENT reconstruction from peculiar velocities (Dekel, Bertschinger & Faber 1990; Dekel et al. 1999; Dekel 2000) yields the mass distribution in our local cosmological neighborhood (even out to $`100h^1\mathrm{Mpc}`$), which in principle enables direct mapping of the local biasing field. However, the sparse and noisy data limit the mass reconstruction to low resolution ($`10h^1\mathrm{Mpc}`$) compared to the volume sampled, which introduces large cosmic scatter in the mass CDF. New accurate data nearby, based on SBF distances (Tonry et al. 1997) do enable a promising resolution of a few Mpc (see Dekel 2000), but limited to inside the local sphere of radius $`30h^1\mathrm{Mpc}`$.

What makes the method proposed here feasible is the fact that the mass CDF is only weakly sensitive to variations in the cosmological scenario within the range of models that are currently considered as viable models for the formation of large-scale structure (e.g., Primack 1998, Bahcall et al. 1999). It has been proposed that the mass PDF can be well approximated by a log-normal distribution in $`\rho /\overline{\rho }=1+\delta `$ (e.g., Coles & Jones 1991; Kofman et al. 1994), and it has since been argued that this approximation becomes poor for certain power spectra and at the tails of the distribution (Bernardeau 1994; Bernardeau & Kofman 1995). In this section, we investigate the robustness of $`C(\delta )`$ for our purpose here, namely, in comparison with the typical difference between the CDFs of galaxies and mass (i.e., the mean biasing function) which we are trying to approximate.

We use for this purpose a suite of $`N`$-body simulations of six different cosmological models. In addition to the two high-resolution simulations of $`\tau `$CDM and $`\mathrm{\Lambda }`$CDM used in the previous section, we have simulated three random realizations of each of the three following models (all using a Hubble constant of $`h=0.5`$): standard CDM (SCDM; $`\mathrm{\Omega }_\mathrm{m}=1`$ with spectral index $`n=1`$), an extreme open CDM (OCDM; $`\mathrm{\Omega }_\mathrm{m}=0.2`$, $`n=1`$), and an extreme tilted CDM (TCDM; $`\mathrm{\Omega }_\mathrm{m}=1`$, $`n=0.6`$). These simulations were run by Ganon et al. (2000, in preparation) using a PM code (by Bertschinger & Gelb 1991), with $`128^3`$ particles in a $`256h^1\mathrm{Mpc}`$ box. The present epoch is defined in these simulations by a linear fluctuation amplitude of $`\sigma _8=1.0`$. A similar simulation was run using a constrained realization (CR) of the local universe based on the galaxy density in the IRAS 1.2Jy redshift survey under the assumption of no biasing (Kolatt et al. 1996), with $`\mathrm{\Omega }_\mathrm{m}=1`$ and the present defined in this case by $`\sigma _8=0.7`$.

Figure 4 (left) shows for the different models the deviations $`\mathrm{\Delta }C(\delta )`$ of the mass CDFs, smoothed TH8, from a cumulative log-normal distribution with the same $`\sigma `$. The log-normal probability density is

$$P_{\mathrm{ln}}(\delta )=\frac{1}{\stackrel{~}{\rho }}\frac{1}{\sqrt{2\pi }s}\mathrm{exp}\left[\frac{(\mathrm{ln}\stackrel{~}{\rho }m)^2}{2s^2}\right],$$
(7)
where

$$\stackrel{~}{\rho }=1+\delta ,m=0.5\mathrm{ln}(1+\sigma ^2),s^2=\mathrm{ln}(1+\sigma ^2)\mathrm{and}\sigma ^2=\delta ^2.$$
(8)
The cumulative log-normal distribution is obtained by integration,

$$C_{\mathrm{ln}}(\delta )=\mathrm{erf}\left[\frac{\mathrm{ln}\stackrel{~}{\rho }m}{s}\right],$$
(9)
where

$$\mathrm{erf}(x)\frac{1}{\sqrt{2\pi }}_{\mathrm{}}^xe^{t^2/2}𝑑t.$$
(10)
For the cases of OCDM, TCDM and SCDM, the CDF is obtained from the three simulations of each model put together. The errors are similar in the different cases; we therefore plot representative error bars only for the $`\tau `$CDM case.

In all the realizations that had random Gaussian initial conditions, the deviation from lognormality is less than 2%. The constrained realization shows somewhat larger deviations, but even in this case they never exceed 5%. These deviations are indeed smaller than the typical differences between $`C_\mathrm{g}(\delta )`$ and $`C(\delta )`$, which are on the order of 10% (see Figure 1).

In order to evaluate how important the contribution of $`\mathrm{\Delta }C`$ is to the error in the recovery of $`\delta _\mathrm{g}|\delta `$, we compare in the right panel of Figure 4 the true $`\delta _\mathrm{g}|\delta `$ in the $`\tau `$CDM simulation with two approximations $`\delta _\mathrm{g}(\delta )`$ based on equation (4), one using the true matter CDF and the other replacing it with a cumulative log-normal distribution of the same $`\sigma `$. The results of the two approximations are very similar; the differences between them seem to be much smaller than the differences between each of them and the true biasing function $`\delta _\mathrm{g}|\delta `$. We can conclude that for the purpose of recovering the biasing function, for the range of Gaussian cosmological models considered, $`C_{\mathrm{ln}}`$ is a good approximation for $`C`$.

The proximity of $`C`$ and $`C_{\mathrm{ln}}`$ could have been alternatively evaluated by the Kolmogorov-Smirnov (KS) statistic, $`D=\mathrm{max}\{|\mathrm{\Delta }C|\}`$. For computing the KS significance $`q(D)`$, we estimate the effective number of “independent” points by $`N_{\mathrm{eff}}=V_{\mathrm{box}}/V_{\mathrm{win}}`$, where $`V_{\mathrm{box}}`$ is the volume of the simulation box and $`V_{\mathrm{win}}`$ is the effective volume of the smoothing window. A value of $`q1`$ ($`D1`$) corresponds to a good match, and $`q1`$ ($`D1`$) to a poor match. For our $`\tau `$CDM simulation, with TH8 smoothing at $`z=0`$ and 1, we obtain $`D0.01`$ and $`q>0.9999`$, confirming that $`C_{\mathrm{ln}}`$ is a good fit. However, for the larger SCDM and OCDM simulations, although $`D`$ is still only $`0.015`$, the corresponding $`q`$ values are at the level of only a few percent. For TCDM and CR, where $`D`$ is 0.016 and 0.052 respectively, the values of $`q`$ drop to the level of a fraction of a percent, and the discrepancy seems large. This KS test indicates that the log-normal approximation is not always perfect for general purpose, as has been argued in the literature. However, our direct tests reported above demonstrate that the use of the log-normal approximation is adequate for the recovery of the mean biasing function in all these cases.

We comment in passing that while the mass CDF is well approximated for our purpose by a log-normal distribution, the shape of the halo (or galaxy) CDF is usually far from a log-normal shape. This is implied by equation (4), from which it follows that $`C_\mathrm{g}(\delta _\mathrm{g})=C[\delta _\mathrm{g}^1(\delta _\mathrm{g})]`$. One does not expect to recover a log-normal distribution from a general functional form for $`\delta _\mathrm{g}^1`$. In particular, the linear biasing model, which seems to be an acceptable approximation in some cases with large smoothing (e.g., IRAS 1.2Jy galaxies at 12$`h^1\mathrm{Mpc}`$ Gaussian smoothing; Sigad et al. 1998), leads to a $`C_\mathrm{g}(\delta _\mathrm{g})`$ that is far from log-normal. Trying to evaluate the difference between $`C_\mathrm{g}`$ and a log-normal distribution using the KS statistic, we obtain for the halos in the $`\tau `$CDM simulation, with TH8 smoothing, both at $`z=0`$ and 1, $`D0.08`$ and $`q0.05`$, namely a poor fit compared to the $`q1`$ of $`C`$ vs $`C_{\mathrm{ln}}`$. Similar conclusions are valid for galaxies.

Our method for measuring the nonlinear biasing function requires an assumed value of $`\sigma `$. Since $`\sigma `$ is known only to a limited accuracy (§6), we should check the robustness of our results to errors in $`\sigma `$. We repeated the reconstruction described in §2, both for halos and for galaxies, with perturbed values of $`\sigma `$ in a range $`\pm 20\%`$ about the true value of the simulation. Not surprisingly, we find that the analog of the linear biasing parameter, $`\widehat{b}`$, varies roughly in proportion to $`\sigma ^1`$. We also find that $`\stackrel{~}{b}`$ varies in a similar way, such that the ratio $`\stackrel{~}{b}/\widehat{b}`$, which is the natural measure of nonlinear biasing (Dekel & Lahav 1999), is a very weak function of $`\sigma `$, roughly $`\stackrel{~}{b}/\widehat{b}\sigma ^{0.15}`$. This test indicates that our method provides a robust measure of the nonlinearity in the biasing scheme, that is to a large extent decoupled from the uncertainty in the linear biasing parameter.

## 4 REDSHIFT DISTORTIONS

The densities as measured in redshift space (z-space) are in general different from the real-space (r-space) densities addressed so far, because the radial peculiar velocities distort the volume elements along the lines of sight. One approach to deal with redshift distortions is to start by recovering the full galaxy density field in r-space, using the linear or a mildly-nonlinear approximation to gravitational instability (e.g., Yahil et al. 1991; Strauss et al. 1992; Fisher et al. 1995; Sigad et al. 1998), and then compute the biasing function in r-space as outlined above. The accuracy of such a procedure would be limited by the approximation used for nonlinear gravity. Another difficulty with this approach is that it requires one to assume a priori a specific biasing scheme, already in the force calculation that enters the transformation from z-space to r-space, while this biasing scheme is the very unknown we are after; this would require a nontrivial iterative procedure.

The alternative we propose here is to actually use the z-space CDF, $`C_{\mathrm{g},\mathrm{z}}(\delta _{\mathrm{g},\mathrm{z}})`$, as provided directly from counts in cells of galaxies in a redshift survey. If the redshift distortions affect the densities of galaxies and mass in a similar way, then one may expect the biasing function in z-space to be similar to the one in real space,

$$\delta _{\mathrm{g},\mathrm{z}}|\delta _\mathrm{z}=\delta =\delta _\mathrm{g}|\delta .$$
(11)
If we only had a robust functional form for the mass CDF in z-space, $`C_\mathrm{z}(\delta _\mathrm{z})`$, then we could compute the desired biasing function all in z-space, using equation (4) but with the analogous z-space quantities. We thus need to test the validity of equation (11), and come up with a useful approximation for $`C_\mathrm{z}(\delta _\mathrm{z})`$.

Figure 5 illustrates the accuracy of equation (11). It compares the biasing functions in z-space and r-space, as derived via equation (4) and its z-space analog from the corresponding CDFs of halos and mass in the $`\tau `$CDM simulation with TH8 smoothing. The two curves are remarkably similar for $`\delta <0.60.8`$, roughly out to the 1-sigma rms fluctuation value. This is roughly the range where the biasing scatter is reasonably small and our basic method is applicable (§2, Figure 1). The curves deviate gradually as $`\delta `$ increases, partly due to stronger “fingers of god” effects at high densities. The deviation is somewhat weaker for larger-mass halos (perhaps due to a lower velocity dispersion for more massive objects as a result of dynamical friction).

The direction of the deviation from equation (11), as seen in Figure 5, can be obtained by applying linear theory of redshift distortions to the case of linear biasing in r-space, $`\delta _\mathrm{g}=b\delta `$. In linear theory, the density fluctuations in r-space and z-space are related via $`\delta _\mathrm{z}=\delta [1+f(\mathrm{\Omega }_\mathrm{m})\mu ^2]`$, where $`f(\mathrm{\Omega }_\mathrm{m})\mathrm{\Omega }_\mathrm{m}^{0.6}`$ (with a negligible dependence on $`\mathrm{\Omega }_\mathrm{\Lambda }`$, see Lahav et al. 1991) and $`\mu `$ is the cosine of the angle between the galaxy velocity vector and the line of sight. If the galaxies obey the continuity equation, then $`\delta _{\mathrm{g},\mathrm{z}}\delta _\mathrm{g}=\delta _\mathrm{z}\delta `$, which implies the following biasing relation in z-space:

$$\delta _{\mathrm{g},\mathrm{z}}=\frac{b+f(\mathrm{\Omega }_\mathrm{m})\mu ^2}{1+f(\mathrm{\Omega }_\mathrm{m})\mu ^2}\delta _\mathrm{z}.$$
(12)
Averaging over all possible directions and assuming $`\mathrm{\Omega }_\mathrm{m}=1`$, we find that the linear biasing parameter in z-space is predicted to be $`b_z=(3b+2)/5`$ for the case shown in Figure 5. This indicates that the linear biasing parameter tends to be closer to unity in z-space than in r-space. Based on our empirical tests of equation (11), we learn that the nonlinear effects (of biasing and gravity) conspire to make equation (11) a better approximation than implied by the linear approximation.

Note that while the results of Figure 5 based on our high-resolution $`\tau `$CDM simulation are quite accurate in the way they treat halos, they may suffer from significant cosmic variance due to the relatively small volume sampled, where the presence (or absence) of a few “fingers of god” could strongly affect the biasing function in the high-$`\delta `$ regime. To test the validity of equation (11) with reduced cosmic variance, we appeal to yet another set of N-body simulations (by Cole et al. 1997) which cover a much larger volume but with lower resolution. These simulations followed the evolution of $`N=192^3`$ particles in a periodic box of comoving side $`L=345.6h^1\mathrm{Mpc}`$ using an Adaptive P3M code. The cosmological models are $`\mathrm{\Lambda }`$CDM ($`\mathrm{\Omega }_\mathrm{m}=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, $`h=0.65`$, cluster-normalized to $`\sigma _8=1.05`$) and $`\tau `$CDM ($`\mathrm{\Omega }=1`$, $`h=0.25`$, cluster-normalized to $`\sigma _8=0.55`$). Nine mock catalogs were extracted from each of the parent simulations, each containing $`510^5`$ particles in a box of $`L=200h^1\mathrm{Mpc}`$. The partial overlap between the catalog volumes is thus about $`50\%`$. The central “observer” was chosen to mimic certain properties of the Local Group environment (see Branchini et al. 1999). Since the resolution of these large simulations is inadequate for a detailed halo identification based on many simulated particles in each halo, we identify individual particles as galaxies using a Monte-Carlo procedure in which the galaxies are chosen to make a random realization of an assumed nonlinear biasing function. Here we adopt the biasing function proposed by Dekel & Lahav (1999) to fit the simulated results of Somerville et al. (2000):

$$\delta _\mathrm{g}(\delta )=\left\{\begin{array}{cc}(1+b_0)(1+\delta )^{b_{\mathrm{neg}}}1\hfill & \delta 0\hfill \\ b_{\mathrm{pos}}\delta +b_0\hfill & \delta >0\hfill \end{array}\right\},$$
(13)
with $`b_{\mathrm{neg}}=2`$ and $`b_{\mathrm{pos}}=1`$. The mass density field is obtained with a Gaussian smoothing of radius $`5h^1\mathrm{Mpc}`$ at the points of a $`128^3`$ cubic grid inside a box of size $`200h^1\mathrm{Mpc}`$. Galaxy densities are obtained at the grid points based on equation (13), and then interpolated to the galaxy positions as defined by the selected particles. Given the appropriate probability distributions $`P(\delta )`$, the value of $`b_0`$ is determined for each choice of the parameters $`b_{\mathrm{neg}}`$ and $`b_{\mathrm{pos}}`$ such that $`\delta _\mathrm{g}=0`$ as required by definition. We obtain $`b_0=0.26`$ and $`b_0=0.19`$ for the models of $`\mathrm{\Lambda }`$CDM and $`\tau `$CDM respectively.

Figure 6 compares the CDFs and associated biasing functions in r-space and z-space, averaged over nine mock catalog from the large-box $`\mathrm{\Lambda }`$CDM simulation. The z-space biasing function is indeed almost indistinguishable from the r-space one (bottom panels); the differences are typically on the order of a couple of percents. The results for $`\tau `$CDM are similar.

In order to quantify this difference further, we define a statistic analogous to equation (5):

$$\mathrm{\Delta }=\frac{1}{N_{\mathrm{bins}}\sigma _\mathrm{g}^2}\underset{\delta \mathrm{bins}}{}[\delta _{\mathrm{g},\mathrm{z}}(\delta _\mathrm{z}=\delta )\delta _\mathrm{g}(\delta )]^2,$$
(14)
in which the first and second terms are the biasing functions as derived from the CDFs in z-space and r-space respectively. The summation is over bins with $`\delta <\delta _{\mathrm{max}}`$, such that $`99\%`$ of the volume is accounted for. We also compute the two moments of the observed biasing function $`\widehat{b}_{\mathrm{obs}}`$ and $`\stackrel{~}{b}_{\mathrm{obs}}`$. These three quantities, averaged over the mock catalogs, are listed in Table 2 (second column). Their deviation from the “true” values (Table 2, first column) is the systematic error. The quoted errors refer to the 1$`\sigma `$ scatter about the mean; they represent the random errors. The results are listed for the two models, $`\mathrm{\Lambda }`$CDM and $`\tau `$CDM. We conclude that the biasing function and its moments, as computed from the z-space CDFs, resemble those computed from the r-space CDFs to within 2%. Note that the Monte Carlo procedure we use to generate mock catalogs artificially reduces the amount of clustering and over-smoothes the density fields for dark and luminous particles. The net effect is to decrease the biasing moments by $`7\%`$, relative to the values implied by the biasing scheme, equation (13). However, this bias does not affect the present analysis for which “true” values are obtained from the mock catalogs themselves and not from equation (13).

Our next task is to come up with a robust CDF for the mass in z-space. We try the same log-normal distribution that was found robust for our purpose in r-space (§3), but with a proper rms in z-space, $`\sigma _\mathrm{z}`$. Based on the linear approximation for Gaussian fields in the small-angle limit (Kaiser 1987), we express $`\sigma _\mathrm{z}`$ in terms of $`\sigma `$ and $`\mathrm{\Omega }_\mathrm{m}`$ of the cosmological model by:

$$\sigma _\mathrm{z}=\left[1+\frac{2}{3}f(\mathrm{\Omega }_\mathrm{m})+\frac{1}{5}f^2(\mathrm{\Omega }_\mathrm{m})\right]^{{\scriptscriptstyle \frac{1}{2}}}\sigma .$$
(15)
We thus approximate the z-space biasing function by $`\delta _{\mathrm{ln},\mathrm{z}}(\delta _\mathrm{z})`$, as derived from the z-space CDFs but where the mass CDF is replaced by a cumulative log-normal distribution function $`C_{\mathrm{ln},\mathrm{z}}`$ (eq. ) with standard deviation $`\sigma _\mathrm{z}`$ (eq. ). The resultant biasing function, averaged over the mock catalogs, is displayed in the bottom panels of Figure 6. We see that for $`\mathrm{\Lambda }`$CDM the differences between $`\delta _{\mathrm{ln},\mathrm{z}}(\delta _\mathrm{z}=\delta )`$ and $`\delta _\mathrm{g}(\delta )`$ are at at the level of a few percent. For $`\tau `$CDM they are only a bit larger; they exceed 10% but only near $`\delta 2`$, at the tail of the distribution. The error in the biasing function $`\mathrm{\Delta }`$ defined in analogy to equation (14), and the biasing moments, are listed in Table 2 (third column, marked “z-space ln”). The systematic error $`\mathrm{\Delta }`$ is still well below 2%, but the biasing parameters are systematically underestimated by 4% and 7% in $`\mathrm{\Lambda }`$CDM and $`\tau `$CDM respectively.

Overall, it seems that our straightforward method deals with redshift distortions fairly well, without any a priori assumption about the biasing scheme.

## 5 SAMPLING ERRORS

The accuracy of the derivation of the galaxy PDF is limited by two observational factors: the finite volume sampled and the mean density of galaxies in the sample.333The additional edge effects can be greatly minimized by using a volume-limited sample and a proper choice of cell coverage (see Szapudi & Colombi 1996).

In principle, the limited volume is responsible for cosmic variance due to the fact that the long-wavelength fluctuations in the real universe are not fairly represented in the sampled volume. This is not of major concern for us here because (a) it is expected to introduce only a random error, and (b) as long as the biasing is local, the effects of long waves on the PDFs of galaxies and mass are expected to be correlated, making the local biasing function representative of the universal function despite the relatively small sampling volume.

More important is the shot noise introduced by the combination of volume and sampling density effects. For a given cell size (or smoothing length), the error can be divided into the error in the count within each cell and the error due to the finite number of cells in the sample volume. These shot-noise sources may introduce both random and systematic errors. We evaluate them by computing the mean and standard deviation over a suite of mock catalogs in which we vary either the volume or the sampling density for a fixed smoothing scale.

With TH8 smoothing, our mock catalogs from the large $`\mathrm{\Lambda }`$CDM simulation contain $`N_{\mathrm{eff}}3700`$ independent cells. However, the currently available redshift surveys allow an analysis in a much smaller volume. For example, a volume-limited subsample from the PSC$`z`$ catalog (Saunders et al. 2000), that is cut at a distance where the average galaxy separation is $`l=8h^1\mathrm{Mpc}`$ (i.e., on the order of our smoothing scale), contains only $`600`$ independent cells. We therefore estimate the error associated with reducing the sampled volume such that $`N_{\mathrm{eff}}600`$ in each mock catalog. We select from the simulation 9 such non-overlapping sub-volumes, while keeping the sampling density and smoothing length fixed. The results for $`\mathrm{\Lambda }`$CDM, averaged over the mock catalogs, are shown in the upper panels of Figure 7, and the results for the two cosmological models are summarized in Table 2 (column 4). We find no significant systematic errors due to the volume effect in a sample like PSC$`z`$ and with $`8h^1\mathrm{Mpc}`$ smoothing (except in the very high-$`\delta `$ tail for $`\tau `$CDM). The corresponding random errors in the biasing parameters are $`5\%`$ and $`6\%`$ for $`\mathrm{\Lambda }`$CDM and $`\tau `$CDM respectively.

The sampling density can be parameterized by the mean galaxy separation, $`l`$. In our large simulation $`l=2.5h^1\mathrm{Mpc}`$, much smaller than the smoothing length of $`8h^1\mathrm{Mpc}`$, but in real samples $`l`$ could be on the order of the smoothing length. To test the effect of sampling density, we produce 9 mock catalogs in which galaxies are sub-sampled at random from the original catalog such that the mean separation is $`l=`$ 6, 8, or $`10h^1\mathrm{Mpc}`$, while the smoothing length and large volume are kept fixed with $`N_{\mathrm{eff}}3700`$. The results for $`\mathrm{\Lambda }`$CDM are shown in the bottom panels of Figure 7, and for the two models in Table 2 (columns 5-7). We see that the sparse sampling artificially enhances both positive and negative density fluctuations, which enlarges the width of the galaxy PDF. This results in a steeper biasing function. For $`\mathrm{\Lambda }`$CDM, the effect becomes noticeable only when $`l8h^1\mathrm{Mpc}`$, where the systematic error in the biasing parameters is of order 10% and larger, and $`\mathrm{\Delta }`$ is of order a few percent. For $`\tau `$CDM the sampling-density effect is noticeable already for $`l6h^1\mathrm{Mpc}`$, with the error reaching $`3050\%`$ at $`ł10h^1\mathrm{Mpc}`$. A plausible explanation for why the sparse sampling is more damaging in the $`\tau `$CDM model is that the clustering in this model is weaker ($`\sigma _8`$ is smaller to match the cluster abundance which constrains $`\sigma _8\mathrm{\Omega }^{0.5}`$), and therefore the high-density regions are poorly sampled by galaxies.

In summary: the main source of error in our analysis is the sparse sampling. For recovering the biasing function with TH8 smoothing, the mean separation should be $`8h^1\mathrm{Mpc}`$.

## 6 CONCLUSION

We propose a simple prescription for recovering the mean nonlinear biasing function from a large redshift survey. The biasing function is defined by $`b(\delta )\delta =\delta _\mathrm{g}|\delta `$, and is characterized to second order by two parameters, $`\widehat{b}`$ and $`\stackrel{~}{b}`$, measuring the mean biasing and its nonlinearity respectively. The method is applied at a given cosmology, time, object type and smoothing scale, and involves one parameter that should be assumed a priori — the rms mass density fluctuation $`\sigma `$ on the relevant scale.

The main steps of the algorithm are as follows:

1. Obtain the observed cumulative distribution function in redshift space $`C_{\mathrm{g},\mathrm{z}}(\delta _{\mathrm{g},\mathrm{z}})`$, by counts in cells or with window smoothing at a certain smoothing length.
2. Assume a value for $`\sigma `$ on that scale and for the cosmological density parameter $`\mathrm{\Omega }_\mathrm{m}`$, and approximate the mass CDF in redshift space by $`C_{\mathrm{ln},\mathrm{z}}(\delta _\mathrm{z};\sigma _\mathrm{z})`$, the cumulative log-normal distribution (eq. ), with the width $`\sigma _\mathrm{z}`$ derived from $`\sigma `$ and $`\mathrm{\Omega }_\mathrm{m}`$ by equation (15).
3. Derive the mean biasing function by

$$\delta _\mathrm{g}(\delta =\delta _\mathrm{z})\delta _{\mathrm{g},\mathrm{z}}(\delta _\mathrm{z})=C_{\mathrm{g},\mathrm{z}}^1[C_{\mathrm{ln},\mathrm{z}}(\delta _\mathrm{z};\sigma _\mathrm{z})].$$
(16)

We first showed that the mean biasing function, at TH8 smoothing, can be derived with reasonable accuracy from the r-space CDFs of galaxies (or halos) and mass, despite the biasing scatter. We then demonstrated that for a wide range of CDM cosmologies the mass CDF can be properly approximated for this purpose by a log-normal distribution of the same width $`\sigma `$. Next we showed that the biasing functions in z-space and r-space are very similar, and that the z-space mass CDF can also be approximated by a log-normal distribution, with a width derived from $`\sigma `$ via equation (15). This allows us to apply the method directly to the observed CDF in a redshift survey. The errors in the recovered biasing function and its moments, in an ideal case of dense sampling in a large volume, are at the level of a few percent.

In any realistic galaxy survey the limited volume and discrete sampling introduce further random and systematic errors. For a survey like the PSC$`z`$ survey, the main source of error is the sampling density; the error does not exceed $`10\%`$ as long as the mean observed galaxy separation is kept smaller than the smoothing radius. We are currently in the process of applying this method to the PSC$`z`$ survey (E. Branchini, et al. 2000, in preparation), where a more specific error analysis will be carried out. The sampling errors are expected to be significantly smaller for the upcoming 2dF and SDSS redshift surveys.

In §2 we showed that our method works well both for halos and for galaxies, on scales 5 to 15$`h^1\mathrm{Mpc}`$, and in the redshift range $`0z3`$ over which the biasing is expected to change drastically. We obtain a similar accuracy when we vary the cosmological model, the mass of the halos in the comparison, or galaxy properties such as morphological type and luminosity. The approximation $`\delta _\mathrm{g}(\delta )`$ is consistent (the deviation is less than 1-$`\sigma `$) with the true average biasing function $`\delta _\mathrm{g}|\delta `$ over a wide range of $`\delta `$ values, which covers 98 – 99% of the volume, depending on redshift and the type of biased objects. This allows us to estimate the moments of the biasing function to within a few percent (see Table 1). The moments of the biasing function are derived from 99.9% of the volume (99% at $`z`$=3 and for relative biasing).

The method requires as external parameters the rms mass-density fluctuation $`\sigma `$ and the cosmological parameter $`\mathrm{\Omega }_\mathrm{m}`$. These can be obtained by joint analyses of constraints from several observational data sets, such as the cluster abundance (e.g., Eke et al. 1998), peculiar velocities (e.g., Dekel & Rees 1994; Zaroubi et al. 1997; Freudling et al. 1999), CMB anisotropies (e.g., de Bernardis et al. 1999), and type Ia supernovae (Riess et al. 1998; Perlmutter et al. 1999). Examples for such joint analyses are Bahcall et al. (1999) and Bridle et al. (1999).

The method is clearly applicable at $`z0`$ with available redshift surveys and especially with those that will become available in the near future, 2dF and SDSS. In the future, this method may become applicable at higher redshifts as well, where the biasing plays an even more important role. With the accumulation of Lyman-break galaxies at $`z3`$, it may soon become feasible to reconstruct their PDF by counts in cells, and our method will allow a recovery of the biasing function at this early epoch, with consequences on galaxy formation and on the evolution of structure.

We have concentrated here on smoothing scales relevant to galaxy biasing, but the method may also be applicable for the biasing of galaxy clusters, on scales of a few tens of Mpc. The biasing scatter may be larger for clusters because of their sparse sampling, but the larger mean biasing parameter for clusters may help in regaining the required monotonicity for equation (4) to provide a valid approximation to the mean biasing function. The mass PDF has been checked to be properly approximated by a log-normal distribution at smoothing scales in the range 20 to $`40h^1\mathrm{Mpc}`$, using simulations of the standard CDM and Cold+Hot DM models (Borgani et al. 1995). The errors due to sparse sampling would require a smoothing scale at the high end of this range.

In a large redshift survey which distinguishes between object types, one can measure the relative biasing function between two object types by applying equation (6) in redshift space, using the observed CDFs for the two types without appealing to the underlying mass distribution at all. The upcoming large redshift surveys 2dF and SDSS, and the DEEP survey at $`z1`$, are indeed expected to provide adequate samples of different galaxy types. Compared with the predictions of simulations and semi-analytic modeling of galaxy formation (e.g., Kauffmann et al. 1999; Benson et al. 1998; Baugh et al. 1999; Somerville & Primack 1999), the measured relative biasing function can provide valuable constraints on the formation of galaxies and the evolution of structure.

While implementing the method outlined above for measuring the mean nonlinear biasing function using current and future redshift surveys, the next challenge is to devise a practical method for measuring the biasing scatter about the mean.

###### Acknowledgements.

We thank S. Cole, A. Eldar, G. Ganon, T. Kolatt, R. Somerville and our collaborators in the GIF team, J.M. Colberg, A. Diaferio, G. Kauffmann, and S.D.M. White, for providing simulations and mock catalogs. We thank A. Maller, I. Szapudi and D. Weinberg for stimulating discussions, and V. Narayanan and M. Strauss for a helpful referee report. EB thanks the Hebrew University for its hospitality. This work was supported by the Israel Science Foundation grant 546/98.
Table 1: Recovery of the biasing function from the CDFs  $`\mathrm{\Lambda }`$CDM     halos vs. mass  galaxies vs. mass  early vs. late type   $`z`$=0 $`z`$=1 $`z`$=3  $`z`$=0 $`z`$=1 $`z`$=3  $`z`$=0 $`z`$=1 $`z`$=3  $`\widehat{b}`$ 0.67 1.21 2.98  0.89 1.31 2.38  1.11 1.32 1.28  $`\widehat{b}_\mathrm{a}`$ 0.58 1.25 2.86  0.80 1.32 2.25  1.20 1.38 1.49  $`\stackrel{~}{b}`$ 0.74 1.24 3.04  0.91 1.31 2.40  1.13 1.34 1.30  $`\stackrel{~}{b}_\mathrm{a}`$ 0.75 1.31 3.08  0.90 1.36 2.38  1.35 1.52 1.64  $`\mathrm{\Delta }`$ 0.16 0.14 0.11  0.08 0.08 0.08  0.55 0.38 0.56  $`\tau `$CDM   halos vs. mass  galaxies vs. mass  early vs. late type   $`z`$=0 $`z`$=1 $`z`$=3  $`z`$=0 $`z`$=1 $`z`$=3  $`z`$=0 $`z`$=1 $`z`$=3  $`\widehat{b}`$ 0.90 2.18 6.62  0.93 1.71 4.44  1.17 1.34 1.27  $`\widehat{b}_\mathrm{a}`$ 0.89 2.28 6.75  0.93 1.75 4.32  1.18 1.39 1.50  $`\stackrel{~}{b}`$ 0.93 2.20 7.85  0.95 1.71 4.62  1.18 1.35 1.31  $`\stackrel{~}{b}_\mathrm{a}`$ 0.96 2.30 8.00  0.98 1.76 4.63  1.26 1.46 1.65  $`\mathrm{\Delta }`$ 0.08 0.07 0.20  0.08 0.04 0.08  0.22 0.20 0.54

Table 2: Redshift distortions and sampling errors in the biasing function  $`\mathrm{\Lambda }`$CDM     True z-space z-space ln Volume $`l=6^a`$ $`l=8^a`$ $`l=10^a`$  $`\widehat{b}`$ 1.13 $`1.12\pm 0.006`$ $`1.09\pm 0.02`$ $`1.12\pm 0.05`$ $`1.17\pm 0.05`$ $`1.23\pm 0.04`$ $`1.31\pm 0.06`$  $`\stackrel{~}{b}`$ 1.14 $`1.12\pm 0.006`$ $`1.10\pm 0.02`$ $`1.13\pm 0.05`$ $`1.17\pm 0.05`$ $`1.24\pm 0.04`$ $`1.32\pm 0.06`$  $`\mathrm{\Delta }`$  $`0.001\pm 0.001`$ $`0.002\pm 0.001`$ $`0.005\pm 0.006`$ $`0.006\pm 0.006`$ $`0.016\pm 0.010`$ $`0.049\pm 0.028`$  $`\tau `$CDM   True z-space z-space ln Volume $`l=6^a`$ $`l=8^a`$ $`l=10^a`$  $`\widehat{b}`$ 1.188 $`1.18\pm 0.002`$ $`1.11\pm 0.02`$ $`1.21\pm 0.06`$ $`1.35\pm 0.07`$ $`1.55\pm 0.07`$ $`1.80\pm 0.07`$  $`\stackrel{~}{b}`$ 1.192 $`1.18\pm 0.002`$ $`1.11\pm 0.02`$ $`1.21\pm 0.06`$ $`1.36\pm 0.07`$ $`1.55\pm 0.07`$ $`1.81\pm 0.07`$  $`\mathrm{\Delta }`$  $`0.002\pm 0.0003`$ $`0.016\pm 0.011`$ $`0.072\pm 0.063`$ $`0.177\pm 0.178`$ $`0.563\pm 0.368`$ $`1.505\pm 0.564`$    a in units of $`h^1\mathrm{Mpc}`$