text:
# A universality class in Markovian persistence

## 1 Introduction

Let us consider a stationary random process $`x(t)`$ and suppose that at $`t=0`$ this process is above a given threshold $`X`$. What is the probability $`Q(T)`$ that the process has never crossed the threshold $`X`$ up to time $`t=T`$ ? In many examples of statistical mechanics we observe

$$Q(T)\text{e}^{\theta _XT},T\mathrm{}$$
(1)
where $`\theta _X`$ is a function of the threshold. The question of determining $`\theta _X`$ is called the ”persistence problem” and has recently been the main interest of many papers in the study of out of equilibrium systems. For a short review the reader can see Ref.  and references therein. In many persistence problems, only the case where $`X=x`$ is considered and usually one writes $`\theta `$ for $`\theta _x`$.

It often happens that stationary processes can be obtained from physical nonstationary systems by rescaling the space variable and taking the new time variable $`t=\mathrm{log}\tau `$ , where $`\tau `$ is the physical time. In terms of the latter the persistence probability decays as a power law, thus defining $`\theta _X`$ as a new independant exponent  called the persistence exponent.

In this paper we present a class of stationary Markov processes $`x(t)`$ that share a universal persistence exponent when $`X=x`$. This class of processes is defined by the linear equation of motion

$$\frac{\text{d}x}{\text{d}t}=\beta x+\eta (t)$$
(2)
where $`\eta (t)`$ is a noise term consisting of upward jumps uniformly distributed along the time axis with density $`\rho `$,

$$\eta (t)=\underset{k=\mathrm{}}{\overset{\mathrm{}}{}}z_k\delta (tt_k)$$
(3)
The jump heights $`z_k`$ are independent identically distributed random variables of law $`R(z)`$ such that $`z`$ exists.

If the jump heights are all equal to $`a`$, the limit $`a0`$, $`\rho \mathrm{}`$ with $`\rho a^2`$ fixed transforms the shifted noise $`\eta (t)\rho a`$ into Gaussian white noise. In this case Eq. (2) is the usual Langevin equation, which is therefore included as a limit in our class of processes.

The interest of this work is to solve the persistence problem for the process defined by Eq. (2). Simulations have shown that for $`T`$ large $`Q(T)`$ decays exponentially with an inverse decay rate $`\theta _X`$ depending on the threshold and on the shape of $`R(z)`$, such as a Dirac peak, an exponential, or a power law. However if $`X=x=\rho z/\beta `$ the result $`\theta =\beta `$ appears for any law $`R(z)`$. This equality was known, until now, only for the Langevin equation, with Gaussian $`\eta (t)`$ .

## 2 A simple method for the calculation of $`\theta `$

In this section we show that $`\theta _X`$ is the lowest eigenvalue of a differential operator. Then when the threshold is set to $`X=x`$, the persistence problem is solved straightforwardly.

### 2.1 Transformation into an eigenvalue problem

Let $`P^{}(x,t)`$ be the probability distribution of $`x`$ at time $`t`$. It is the solution of the master equation

$$_tP(x,t)=\beta _x(xP(x,t))+\rho _0^{\mathrm{}}P(xz,t)R(z)\text{d}z\rho P(x,t)$$
(4)
For the persistence problem only the subclass of functions $`x(t)`$ that are greater than $`X`$ during the time interval $`T`$ is relevant. Let $`P_+^{}(x,t)\text{d}x`$ be the probability that $`x(t^{})`$ has not crossed the level $`X`$ for $`t^{}[0,t]`$ (given that $`x(0)>X`$) and that $`x(t)`$ is in the interval $`[x,x+\text{d}x]`$. $`P_+^{}(x,t)`$ is then solution of the same master equation but with the absorbing boundary condition $`P_+^{}(x,t)=0`$ when $`x<X`$. Because of this condition, $`P_+^{}(x,t)`$ will decay to zero as $`t\mathrm{}`$. We will write it as a superposition of decay modes $`P_+^\lambda (x)\text{e}^{\lambda _Xt}`$, where $`P_+^\lambda (x)`$ is an eigenfunction satisfying the eigenvalue equation corresponding to Eq. (4)

$$\lambda _XP_+^\lambda (x)=\beta \frac{\text{d}}{\text{d}x}(xP_+^\lambda (x))+\rho _0^{\mathrm{}}P_+^\lambda (xz)R(z)\text{d}z\rho P_+^\lambda (x)$$
(5)
with the condition $`P_+^\lambda (x)=0`$ when $`x<X`$. Let $`\lambda _X^{}`$ be the eigenvalue of the slowest decaying mode. Then we have

$$Q(T)\text{e}^{\lambda _X^{}T}_X^{\mathrm{}}P_+^\lambda ^{}(x)\text{d}x\text{as}T\mathrm{}$$
(6)
so that

$$\theta _X=\lambda _X^{}$$
(7)
if $`_X^{\mathrm{}}P_+^\lambda ^{}(x)\text{d}x`$ exists.

Finally let us call $`P_+^\theta (x)P_+^\lambda ^{}(x)`$ the persistence mode. It is necessary to have

$$\begin{array}{c}\hfill P_+^\theta (x)0\text{for}xX\text{and}_X^{\mathrm{}}P_+^\theta (x)\text{d}x<\mathrm{}\end{array}$$
(8)
in order to define $`P_+^\theta (x,t)`$ as a probability distribution. Moreover the definition of the persistence problem ensures that $`x(0)`$ is chosen within the stationary distribution of the process. This distribution is such that $`x`$ exists. Therefore we will only focus on solutions $`P_+^\theta (x)`$ with

$$_X^{\mathrm{}}xP_+^\theta (x)\text{d}x<\mathrm{}$$
(9)
As a result solving the persistence problem is now equivalent to calculating the lowest eigenvalue $`\lambda _X^{}`$ of Eq. (5) that satisfies Eqs. (8) and (9).

### 2.2 Calculation of $`\theta `$

We will now show that $`\lambda _x^{}=\beta `$ irrespective of the jump size distribution. To this end we define

$$f_{\lambda X}^{(n)}=_X^{\mathrm{}}x^nP_+^\lambda (x)\text{d}x,n=0,1,\mathrm{}$$
(10)
and we integrate Eq. (5) on $`x`$ from $`X`$ to $`\mathrm{}`$. As $`P_+^\theta (x)=0`$ when $`x<X`$, we obtain the first relation

$$\lambda _Xf_{\lambda X}^{(0)}=\beta XP_+^\lambda (X)$$
(11)
We now multiply Eq. (5) by $`x`$ and integrate as before. After one integration by parts we get

$$\lambda _Xf_{\lambda X}^{(1)}=\beta X^2P_+^\lambda (X)+(\beta +\rho )f_{\lambda X}^{(1)}\rho _X^{\mathrm{}}\text{d}x_0^{\mathrm{}}\text{d}zxP_+^\lambda (xz)R(z)$$
(12)
Using the absorbing boundary condition, we can write

$$_X^{\mathrm{}}\text{d}x_0^{\mathrm{}}\text{d}zxP_+^\lambda (xz)R(z)=_X^{\mathrm{}}\text{d}yP_+^\lambda (y)_0^{\mathrm{}}\text{d}z(y+z)R(z)$$
(13)
and the RHS of this equation simply reduces to $`f_{\lambda X}^{(1)}+zf_{\lambda X}^{(0)}`$. Then Eqs. (12) and (13) give the second relation

$$\lambda _Xf_{\lambda X}^{(1)}=\beta X^2P_+^\lambda (X)+\beta f_{\lambda X}^{(1)}\rho zf_{\lambda X}^{(0)}$$
(14)
We then eliminate $`P_+^\lambda (X)`$ from Eqs. (11) and (14) and obtain finally using $`\rho z=\beta x`$

$$\left(f_{\lambda X}^{(1)}Xf_{\lambda X}^{(0)}\right)\lambda _X=\beta \left(f_{\lambda X}^{(1)}xf_{\lambda X}^{(0)}\right)$$
(15)
which is valid for any $`\lambda _X`$ for which $`f_{\lambda X}^{(1)}<\mathrm{}`$. One can remark that Eqs. (11) and (14), and therefore Eq. (15), could not be obtained if the jump heights were not all positive. So this property is essential in the present work.

Let us now take $`X`$ equal to the average of the process: $`X=x`$. Eq. (15) shows that we must have

$$\lambda _x=\beta \text{or}f_{\lambda x}^{(1)}/f_{\lambda x}^{(0)}=x$$
(16)
If in Eq. (16) we set $`\lambda _x=\lambda _x^{}=\theta `$, we see from Eqs (8) and (10) that $`f_{\theta x}^{(1)}/f_{\theta x}^{(0)}`$ is the average of the process $`x`$ in the persistence mode, *i.e.*, given that it remains always above $`X=x`$. This quantity is therefore larger than $`x`$ and the second one of Eqs. (16) cannot be satisfied. Finally the only possibility is $`\lambda _x^{}=\beta `$ which proves the main result of this note.

For $`\lambda _x\lambda _x^{}`$ Eq. (16) shows that

$$_x^{\mathrm{}}\text{d}x(xx)P_+^\lambda (x)=0\lambda _x>\lambda _x^{}$$
(17)
It follows that for $`\lambda _x\lambda _x^{}`$ any eigenfunction $`P_+^\lambda (x)`$ with $`f_{\lambda X}^{(1)}<\mathrm{}`$ must have at least one zero for $`x<x<\mathrm{}`$. One can also see Eq. (17) as an orthogonality relation between two eigenvectors of Eq. (5) that are $`xx`$ (with the eigenvalue $`2\beta `$) and $`P_+^\lambda (x)`$.

One final remark is necessary. We have studied the persistence problem in the most natural situation, that is when the condition $`f_{\lambda x}^{(1)}<\mathrm{}`$ is fulfilled. However one could be interested in other cases where the initial distributions $`P(x,0)`$ has a large tail for $`x\mathrm{}`$ that do not ensure this condition. If one starts with such a distribution, one can expect the persistence problem to be governed by the shape of this initial distribution, with different answers depending on the choice of the latter. Therefore large initial distributions where $`f_{\lambda X}^{(1)}=\mathrm{}`$ are simply cases of no interest in our study.

### 2.3 Perturbation for $`X`$ close to $`x`$

Very few cases in statistical physics offer the possibility of calculating the quantity

$$x_X=f_{\theta X}^{(1)}/f_{\theta X}^{(0)}$$
(18)
which is the average of the process in the persistence mode. However as $`f_{\theta x}^{(n)}=f_{\beta x}^{(n)}`$ one can assume that

$$f_{\theta X}^{(n)}=f_{\beta x}^{(n)}+\stackrel{~}{f}_{\theta X}^{(n)}$$
(19)
where $`\stackrel{~}{f}_{\theta X}^{(n)}0`$ if $`|Xx|0`$. Rewriting Eq. (15) for $`\lambda _X=\lambda _X^{}`$ one gets

$$\theta _X=\beta \frac{f_{\theta X}^{(1)}xf_{\theta X}^{(0)}}{f_{\theta X}^{(1)}Xf_{\theta X}^{(0)}}$$
(20)
which can also be written

$$\theta _X=\beta \frac{1}{1\frac{Xx}{x_Xx}}$$
(21)
Finally one can suppose that $`X`$ is close to $`x`$ and perturbs Eq. (21) to find

$$\theta _X=\beta (1+\frac{Xx}{x_xx}+\mathrm{})|Xx|0$$
(22)
Unfortunately we cannot go further here. We do not know how to calculate $`x_x`$ except in the Gaussian limit which will be the focus of next section.

## 3 The special case of Gaussian noise

Our goal here is to illustrate the method of Sec. 2 in the special known case of Gaussian process. If one considers constant jumps with $`R(z)=\delta (za)`$ and sets $`\rho a^2=\mathrm{\Gamma }`$, one can take the limit $`a0`$ with $`\mathrm{\Gamma }`$ fixed in Eq. (4) and find

$$_tP(x,t)=\beta _x((xx)P(x,t))+\frac{\mathrm{\Gamma }}{2}_{xx}P(x,t)+𝒪(a)$$
(23)
In this limit Eq. (23) leads to the Fokker-Planck equation for the zero average process $`x(t)x`$ which hence becomes the stationary Gaussian Markovian process also called the Ornstein-Uhlenbeck process. This shows that the Gaussian case, which is of high interest in statistical physics, is included as a special limit in the class of processes studied in the present work.

### 3.1 Calculation of $`\theta `$

For simplicity we will now write $`x`$ instead of $`xx`$ and $`X`$ instead of $`Xx`$. We thus consider the zero average Gaussian process $`x(t)`$ and we want to determine the persistence exponent $`\theta _X`$ related to the probability of having never crossed the level $`X`$ upto time $`T`$. As in Sec. 2 it is equivalent to find the lowest eigenvalue $`\lambda _X^{}`$ of

$$\lambda _XP_+^\lambda (x)=\beta \frac{\text{d}}{\text{d}x}(xP_+^\lambda (x))+\frac{\mathrm{\Gamma }}{2}\frac{\text{d}^2}{\text{d}x^2}P_+^\lambda (x),P_+^\lambda (X)=0$$
(24)
The solution of this equation is

$$P_+^\lambda (x)=\text{e}^{y^2/4}(A_1D_{\lambda /\beta }(y)+A_2D_{\lambda /\beta }(y))$$
(25)
where $`y=\sqrt{2\beta /\mathrm{\Gamma }}x`$ is dimensionless and $`D_\nu `$ is the parabolic cylinder function of index $`\nu `$ . The integration constants $`A_{1,2}`$ are such that

$$A_1D_{\lambda /\beta }(Y)+A_2D_{\lambda /\beta }(Y)=0$$
(26)
with $`Y=\sqrt{2\beta /\mathrm{\Gamma }}X`$. After one performs the same integrations as in Sec. 2, Eq. (24) gives the two relations

$$\lambda _Xf_{\lambda X}^{(0)}=\frac{\mathrm{\Gamma }}{2}\frac{\text{d}P_+^\lambda }{\text{d}x}|_X$$
(27)
$$\lambda _Xf_{\lambda X}^{(1)}=\beta f_{\lambda X}^{(1)}+\frac{\mathrm{\Gamma }}{2}X\frac{\text{d}P_+^\lambda }{\text{d}x}|_X$$
(28)
Eliminating $`\frac{\text{d}P_+^\lambda }{\text{d}x}|_X`$ from Eqs. (27) and (28) one finally gets

$$(f_{\lambda X}^{(1)}Xf_{\lambda X}^{(0)})\lambda _X=\beta f_{\lambda X}^{(1)}$$
(29)
which is strictly equivalent to Eq. (15). In the same way as in Sec. 2 it now follows that if $`X=0`$ then $`\theta =\lambda _0^{}=\beta `$, which is the known result of Gaussian persistence .

Another remark can be made here. The solution of Eq (24) that vanishes at the origin is proportional to (with $`A_1=A_2=A`$)

$$\stackrel{~}{P}(y)=A\text{e}^{y^2/4}(D_{\lambda /\beta }(y)D_{\lambda /\beta }(y))$$
(30)
which can be reduced to

$$\stackrel{~}{P}(y)=A^{}y\text{e}_1^{y^2/2}F_1((1\lambda /\beta )/2,3/2,y^2/2)$$
(31)
where $`A^{}`$ is simply a numerical constant . Using the result

$$_0^{\mathrm{}}\text{e}^uu_1^{1/2}F_1((1\nu )/2,3/2,u)\text{d}u=0\nu >1$$
(32)
one has immediately $`f_{\lambda 0}^{(1)}=0`$ for any $`\lambda >\beta `$, which is a direct proof of Eq. (17) in the Gaussian case.

### 3.2 Perturbation for small $`X`$

Our method offers a faster way than the one presented in , in the Gaussian case, to calculate $`\theta _X`$ when $`X`$ is close to $`x`$ (which is $`0`$ here).

As $`P_+^\beta (x)=A^{}y\text{e}^{y^2/2}`$ , one finds $`x_x=\sqrt{\pi \mathrm{\Gamma }/4\beta }`$ and Eq. (22) gives immediately

$$\theta _X=\beta (1+2\sqrt{\frac{\beta }{\pi \mathrm{\Gamma }}}X+\mathrm{}),X0$$
(33)
which is the well known perturbative result for the Gaussian persistence exponent .

## 4 Application

We give here some examples of Markovian but non-Gaussian processes defined by Eq. (2).

### 4.1 1D Percolation

Let us consider an infinite 1D lattice. Each site can be occupied or empty. At an initial time, say $`t=0`$, we start with a totally empty lattice. Then, at each time step, one site is chosen randomly among non occupied ones and is then set occupied forever. The density of occupied sites $`p`$ is therefore an increasing function of time. It goes from $`0`$ to $`1`$ as time goes from $`0`$ to $`\mathrm{}`$ (in the limit of an infinite lattice, which is the case we consider here). We define clusters in a standard way, that is as sets of adjacent occupied sites. It is then possible to study percolation on this lattice as a dynamical process. de Freitas and dos Santos Lucena  have already performed this study numerically.

The interest here is to consider a semi-infinite lattice and the cluster that contains its origin. When $`p`$ goes from 0 to 1, the size $`s(p)`$ of this cluster will grow from 0 to $`\mathrm{}`$. It is very easy to see that

$$s=\frac{p}{1p}\text{and}\sqrt{\mathrm{\Delta }s^2}=\frac{\sqrt{p}}{1p}$$
(34)
The evolution equation for $`s(p)`$ is also quite simple. This random process is constant except when new sites are connected to the cluster which creates upward jumps of height $`j`$. Hence

$$\frac{\text{d}s}{\text{d}p}=\underset{k=0}{\overset{\mathrm{}}{}}j_k\delta (pp_k)$$
(35)
where $`k`$ is the jump index. If $`\stackrel{~}{R}(j)\mathrm{\Delta }p`$ is the probability to have a jump $`j`$ when the density of occupied sites goes from $`p`$ to $`p+\mathrm{\Delta }p`$ then

$$\stackrel{~}{R}(j)\mathrm{\Delta }p=p^{j1}(1p)\frac{\mathrm{\Delta }p}{1p}$$
(36)
The quantity $`\mathrm{\Delta }p/(1p)`$ is simply the probability to make a jump during $`\mathrm{\Delta }p`$. The normalised distribution law for the jumps is then

$$\stackrel{~}{R}(j)=p^{j1}(1p)$$
(37)
and the density of jumps is $`\rho (p)=1/(1p)`$. As $`p`$ tends to 1, the typical time scale of the dynamics of the process $`s(p)`$, which is $`1/\rho `$, goes to zero. The aim is now to deal with a stationary process. If the new time scale is

$$t=\mathrm{log}\frac{1}{1p}$$
(38)
the jump times $`t_k`$ are now uniformly distributed with a density 1 along the $`t`$ axis. Then we rescale the process $`s(p)`$, setting

$$x(t)=\frac{s(p)}{\sqrt{\mathrm{\Delta }s^2}}\text{and}z_k=\frac{j_k}{\sqrt{\mathrm{\Delta }s^2}}$$
(39)
and we take the limit $`p1`$ or equivalently $`t\mathrm{}`$. The jumps $`z`$ have now the probability density $`R(z)`$ so that $`\stackrel{~}{R}(j)\mathrm{\Delta }j=R(z)\mathrm{\Delta }z`$. In the limit $`p1`$

$$R(z)=\text{e}^z$$
(40)
If we write Eq. (35) in terms of $`x`$, $`t`$, and $`z`$, the evolution equation is finally

$$\frac{\text{d}x}{\text{d}t}=x(t)+\underset{k=0}{\overset{\mathrm{}}{}}z_k\delta (tt_k)$$
(41)
in the limit $`t\mathrm{}`$. This gives a second example of noise of Eq. (3) with an exponential distribution law for the jump height. The value of the persistence exponent is then $`\theta =1`$.

### 4.2 Other examples

A second example of noise defined by Eq (3) is given by the study of random walks. The time $`t(S)`$ needed for the walker to visit $`S`$ distinct sites on a 1D lattice is a random process which, under the same kind of scaling as in the last subsection, satisfies Eq. (2) . The jump height distribution $`R(z)`$ is then a power law with a cutoff. In this case the value of $`\theta `$ is $`2`$. We will not go further about this result in the present paper and the reader can see references  and  for more details.

Finally the process $`x(t)`$ defined by Eq. (2) is simply the electric signal that one would observe in an RC circuit receiving random charges due to, for instance, the presence of radioactive elements. The persistence exponent is then $`\theta =1/RC`$, that is the inverse relaxation time of the circuit.

## 5 Conclusion and discussion

With a very simple calculation we are able to find the persistence exponent of a large class of stationary Markov processes which includes the Gaussian case. As this quantity is universal within this entire class, it demonstrates the existence of universality classes related to non-Gaussian persistence exponents.

To go further one should study the properties of the master equation (4) because the level-dependent exponent $`\theta _X`$ is entirely determined by the knowledge of $`x_X`$. However the nontrivial asymptotic behaviour $`\theta _XX\mathrm{log}X+X^{1/3}(\mathrm{log}X)^{2/3}+\mathrm{}`$, $`X\mathrm{}`$, found in Ref.  when the jumps are constant shows that this investigation is likely to be very difficult.

## Acknowledgments

The author thanks Professor H. J. Hilhorst for numerous fruitful discussions about the present work and also Professor L. dos Santos Lucena of the Universidade Federal do Rio Grande do Norte, Natal, Brazil for pointing out the relation between 1D percolation and noise with discrete upward jumps.


format: text

text:
# The Bispectrum of IRAS Redshift Catalogs

## 1 Introduction

Observations suggest that the distribution of matter on cosmological scales is nearly homogeneous, with structure that in our present understanding has grown by gravitational instability from small, nearly Gaussian, initial fluctuations. On large scales, departures from homogeneity are small, and their evolution can be studied in perturbation theory. To linear order in perturbation theory, an initially Gaussian fluctuation distribution remains Gaussian, but at higher orders this is no longer true. Details of the departures depend on the nonlinear couplings; thus, non-Gaussianity can probe the dynamics of structure formation. One useful non-Gaussian measure is the three-point correlation function, or the bispectrum in the Fourier domain, a quantity that vanishes identically for a Gaussian distribution but is nonzero and has a characteristic behavior as a function of the shape of the triangle defined by its three wavevector arguments if large-scale structure is formed as a result of gravitational instability (Fry, 1984). When applied to the galaxies, the bispectrum also depends on galaxy biasing, i.e. on the relation between the galaxy and mass distributions (Fry & Gaztañaga, 1993; Fry, 1994) and on any non-Gaussianity of the primordial density fluctuations (Fry & Scherrer, 1994).

Early measurements of the bispectrum in the Lick galaxy catalog did not find the expected gravitational instability signature (Fry & Seldner, 1982). In angular catalogs such as Lick, the density field is integrated over the radial selection function, which in principle can change the dependence on shape, although it appears that in practice the effects of projection on the shape dependence of the bispectrum are not large (Fry & Thomas, 1999; Buchalter, Kamionkowski, & Jaffe, 2000). Results on the moderate scales probed by this catalog are also affected by nonlinear evolution beyond leading order in perturbation theory, which weakens the configuration dependence of the bispectrum in a way that mimics the effects of bias (Scoccimarro et al., 1998). A more recent determination of the angular three-point correlation function in the APM survey (Frieman & Gaztañaga, 1999), which probes significantly larger scales than Lick, shows an amplitude and configuration dependence consistent with gravitational instability from Gaussian initial conditions and no bias, though with large errors.

Recent galaxy redshift surveys have mapped the distribution of galaxies out to large distances. The largest distances are obtained for deep pencil beam and slice surveys; however, surveys with highly anisotropic one- and two-dimensional geometries have window functions which extend to large wavenumber in one or two directions, corresponding to the smallest survey dimension(s). As a consequence, nonlinearity on this scale tends to wash out the configuration dependence of three-point statistics. Thus, although the three-point correlation function has been measured in the Las Campanas Redshift Survey (Jing & Börner, 1998), the survey geometry is effectively two-dimensional, with redshift slices of thickness about $`10h^1\mathrm{Mpc}`$; this scale is not large enough to allow a direct comparison with perturbation theory. Thus, it is very difficult to extract quantitative information from this survey about biasing and the nature of the initial conditions.

To compare with predictions of the configuration dependence of the bispectrum from perturbation theory, we wish to observe the bispectrum directly in an analysis that to the extent possible is not affected by the complications of projection and nonlinear evolution. In this paper, we report measurement of the bispectrum in several galaxy redshift catalogs derived from the IRAS Point Source Catalog (Chester, Beichmann, & Conrow, 1987). Working with three-dimensional (redshift) data avoids the effects of radial projection. Moreover, the large volume and nearly full sky coverage of these surveys allow us to make measurements at small wavenumber (large scales), where clustering is well-described by leading order perturbation theory. Although working in redshift space introduces distortions arising from the modulation of position induced by peculiar velocities, these can be computed (Hivon et al., 1995; Verde et al., 1998; Scoccimarro, Couchman, & Frieman, 1999).

In Section 2 we review the leading order nonlinear perturbation theory predictions for the bispectrum, including its dependence on configuration shape, the distortions introduced by peculiar velocities in redshift space, and the effects of a bias in the galaxy distribution relative to matter. In Section 3 we summarize properties of the IRAS-based redshift catalogs, and we present our estimators in Fourier space, using the optimal weighting of Feldman, Kaiser, & Peacock (1994, FKP). In Section 4 we summarize the data analysis method based on comparison with mock catalogs developed in the companion paper (Scoccimarro, 2000b, Paper I), and in Section 5 we present the results of our analysis. Section 6 contains a final discussion.

## 2 The Bispectrum in Perturbation Theory

We compute statistics in the Fourier domain, since it naturally separates clustering on different physical scales. The Fourier amplitude $`\stackrel{~}{\delta }(\text{k})`$, the transform of the fractional density contrast $`\delta (\text{x})`$, is defined by

$$\stackrel{~}{\delta }(\text{k})=d^3x\delta (\text{x})e^{i\text{k}\text{x}},\delta (\text{x})=\frac{d^3k}{(2\pi )^3}\stackrel{~}{\delta }(\text{k})e^{i\text{k}\text{x}},$$
(1)
where the mass density is $`\rho (\text{x})=\overline{\rho }[1+\delta (\text{x})]`$. The first moments of $`\stackrel{~}{\delta }(\text{k})`$ are the power spectrum and the bispectrum,

$`\stackrel{~}{\delta }(\text{k}_1)\stackrel{~}{\delta }(\text{k}_2)`$ $`=`$ $`(2\pi )^3\delta _\mathrm{D}(\text{k}_1+\text{k}_2)P(k_1),`$ (2)
$`\stackrel{~}{\delta }(\text{k}_1)\stackrel{~}{\delta }(\text{k}_2)\stackrel{~}{\delta }(\text{k}_3)`$ $`=`$ $`(2\pi )^3\delta _\mathrm{D}(\text{k}_1+\text{k}_2+\text{k}_3)B(\text{k}_1,\text{k}_2,\text{k}_3).`$ (3)
The “momentum-conserving” Dirac $`\delta `$-function reflects statistically the underlying homogeneity and implies that the bispectrum is defined for configurations of wavenumbers that form closed triangles. If the initial fluctuations are Gaussian, then the bispectrum vanishes in linear perturbation theory; however, mode-mode couplings when nonlinearities are included generate a non-zero $`B`$ even for an initially Gaussian distribution.

In gravitational instability with Gaussian initial conditions, the bispectrum at second (leading nonzero) order in perturbation theory is (Fry, 1984)

$$B_{123}=Q_{12}P_1P_2+Q_{13}P_1P_3+Q_{23}P_2P_3,$$
(4)
where $`P_1`$ denotes $`P(k_1)`$, etc., and $`B_{123}=B(k_1,k_2,k_3)`$. Here $`Q_{ij}`$ is determined entirely by the configuration geometry,

$$Q_{ij}=(1+\kappa )+\mathrm{cos}\theta _{ij}\left(\frac{k_i}{k_j}+\frac{k_j}{k_i}\right)+(1\kappa )\mathrm{cos}^2\theta _{ij},$$
(5)
where $`\mathrm{cos}\theta _{ij}=\widehat{\text{k}}_i\widehat{\text{k}}_j`$. The factor $`\kappa `$ is $`\kappa =\frac{3}{7}`$ in an Einstein-de Sitter universe, with critical matter density $`\mathrm{\Omega }_m=1`$, and it never varies much from this value. For an open universe with $`\mathrm{\Omega }_m`$ near 1, $`\kappa \frac{3}{7}\mathrm{\Omega }_m^{2/63}`$, and for a flat universe with cosmological constant, $`\kappa \frac{3}{7}\mathrm{\Omega }_m^{1/143}`$ (Bouchet et al., 1992, 1995); while in both cases $`\kappa \frac{1}{2}`$ as $`\mathrm{\Omega }_m0`$.

It is useful to define the reduced bispectrum,

$$Q\frac{B_{123}}{P_1P_2+P_1P_3+P_2P_3}.$$
(6)
For a power-law spectrum, $`Pk^n`$, at second order in perturbation theory, the reduced bispectrum $`Q`$ is independent of time and of scale leaving only the dependence on triangle shape. This behavior contrasts sharply with that of $`P`$ and $`B`$, which vary strongly with scale. Moreover, the expressions above imply that $`Q`$ is very insensitive to the value of $`\mathrm{\Omega }_m`$. The solid curve in Figure The Bispectrum of IRAS Redshift Catalogs shows the characteristic dependence of $`Q`$ on configuration shape, here parameterized by the angle $`\theta _{12}`$ between two sides of ratio $`k_2/k_1=\frac{1}{2}`$, for a power-law spectrum with spectral index $`n=1.4`$ and for $`\mathrm{\Omega }_m=1`$. This ratio of sides minimizes the total range of scales covered as $`\theta `$ varies from 0 to $`\pi `$, but in fact $`Q`$ is insensitive to the ratio $`k_2/k_1`$ for $`\theta <\pi /2`$. The characteristic dependence on shape, with an enhancement for colinear configurations ($`\theta =0`$ or $`\theta =\pi `$), reflects the expected anisotropy of gravitational collapse. The dotted curve shows the behavior for $`\mathrm{\Omega }_m0`$. Other curves show redshift space and bias results, which we discuss next.

In redshift space, where radial distances are assigned from recession velocities by Hubble’s law, peculiar velocities induce distortions in clustering statistics (Kaiser 1987; for a review, see Hamilton 1997). In the plane parallel approximation, where the line of sight is taken as a fixed direction $`\widehat{\text{z}}`$, the density in space $`\stackrel{~}{\delta }(\text{k})`$ and redshift space $`\stackrel{~}{\delta }_s(\text{k})`$ are related to linear order by $`\stackrel{~}{\delta }_s(\text{k})=[1+f_1(\mathrm{\Omega })\mu ^2]\stackrel{~}{\delta }(\text{k})`$, where $`f_1=\mathrm{\Omega }^{4/7}`$ and $`\mu k_z/k`$. Averaging over orientations, the power spectrum monopole in redshift space then becomes

$$P_s(k)=\left(1+\frac{2}{3}f_1+\frac{1}{5}f_1^2\right)P(k),$$
(7)
while the bispectrum monopole is of the same form as as in equation (4) above, with (Hivon et al., 1995)

$`Q_{ij,s}`$ $`=`$ $`a_0+a_2\mathrm{cos}^2\theta _{ij}+a_4\mathrm{cos}^4\theta _{ij}+(a_1\mathrm{cos}\theta _{ij}+a_3\mathrm{cos}^3\theta _{ij})\left({\displaystyle \frac{k_i}{k_j}}+{\displaystyle \frac{k_j}{k_i}}\right)`$ (8)
$`+(1\mathrm{cos}^2\theta _{ij})\left[{\displaystyle \frac{x_1(k_i^2+k_j^2)+x_2k_ik_j\mathrm{cos}\theta _{ij}}{k_i^2+k_j^2+2k_ik_j\mathrm{cos}\theta _{ij}}}\right]`$
where the coefficients are

$`a_0`$ $`=`$ $`(1+\kappa )+{\displaystyle \frac{2}{3}}(2+\kappa )f_1+{\displaystyle \frac{1}{15}}(11+\kappa )f_1^2+{\displaystyle \frac{6}{35}}f_1^3+{\displaystyle \frac{6}{315}}f_1^4+\kappa f_2\left({\displaystyle \frac{1}{3}}+{\displaystyle \frac{6}{15}}f_1+{\displaystyle \frac{1}{35}}f_1^2\right)`$ (9)
$`a_1`$ $`=`$ $`1+{\displaystyle \frac{4}{3}}f_1+{\displaystyle \frac{14}{15}}f_1^2+{\displaystyle \frac{12}{35}}f_1^3+{\displaystyle \frac{15}{315}}f_1^4`$ (10)
$`a_2`$ $`=`$ $`(1\kappa )+{\displaystyle \frac{2}{3}}(2\kappa )f_1+{\displaystyle \frac{1}{15}}(23+\kappa )f_1^2+{\displaystyle \frac{30}{35}}f_1^3+{\displaystyle \frac{48}{315}}f_1^4\kappa f_2\left({\displaystyle \frac{1}{3}}+{\displaystyle \frac{6}{15}}f_1{\displaystyle \frac{3}{35}}f_1^2\right)`$ (11)
$`a_3`$ $`=`$ $`{\displaystyle \frac{4}{15}}f_1^2+{\displaystyle \frac{8}{35}}f_1^3+{\displaystyle \frac{20}{315}}f_1^4`$ (12)
$`a_4`$ $`=`$ $`{\displaystyle \frac{2}{15}}(1\kappa )f_1^2+{\displaystyle \frac{4}{35}}f_1^3+{\displaystyle \frac{16}{315}}f_1^4{\displaystyle \frac{4}{35}}\kappa f_1^2f_2`$ (13)
$`x_1`$ $`=`$ $`{\displaystyle \frac{2}{15}}\kappa f_1f_2,x_2={\displaystyle \frac{4}{35}}\kappa f_1^2f_2;`$ (14)
$`\kappa `$ is as given before, and $`f_1`$ and $`f_2`$ are logarithmic derivatives of the first and second order perturbation theory growth factors $`D_1`$ and $`D_2`$ with respect to the cosmological scale factor, $`f_1=d\mathrm{log}D_1/d\mathrm{log}a`$, etc. For $`\mathrm{\Omega }_m`$ near 1, $`f_1\mathrm{\Omega }^{4/7}`$, and $`f_22f_1`$. Setting $`f_1=f_2=0`$ gives $`a_0=1+\kappa `$, $`a_1=1`$, $`a_2=1\kappa `$, and $`a_3=a_4=x_1=x_2=0`$, reproducing equation (5), the result in real space.

Figure The Bispectrum of IRAS Redshift Catalogs shows the predictions for $`Pk^n`$ with $`n=1.4`$, a good approximation for the IRAS samples over the range of scales we use. The long-dashed curve shows the perturbation theory result for $`Q_s(\theta )`$ in redshift space for $`\mathrm{\Omega }_m=1`$; the short-dashed curve $`\mathrm{\Omega }_m=0.3`$. The behavior depends on $`\mathrm{\Omega }_m`$, but more weakly than the factor $`f(\mathrm{\Omega })=\mathrm{\Omega }^{4/7}`$ that appears in all peculiar velocities. Between $`\mathrm{\Omega }=1`$ and $`\mathrm{\Omega }=0.3`$, $`f`$ varies by a factor of 2, but $`Q`$ by at most about 15%, while the redshift space result for $`\mathrm{\Omega }_m0`$ is the same as in real space. Over the entire range of $`\mathrm{\Omega }`$ the characteristic behavior is very much the same.

To translate these results for the mass distribution into predictions for the galaxy distribution, we must take into account the effects of bias. In the simplest case, a local deterministic bias model, in which the galaxy density is a function of the local mass density, the galaxy number contrast can be expressed as a series in the mass density contrast,

$$\delta _g(\text{x})=\underset{k=0}{\overset{\mathrm{}}{}}\frac{b_k}{k!}\delta ^k(\text{x}).$$
(15)
To leading order in perturbation theory, the resulting galaxy power spectrum is simply related to the mass power spectrum, $`P_g(k)=b_1^2P(k)`$, where $`b_1`$ is the linear bias factor. In this model, the galaxy bispectrum amplitude is (Fry & Gaztañaga, 1993; Fry, 1994)

$$Q_g=\frac{Q}{b_1}+\frac{b_2}{b_1^2}.$$
(16)
Thus, the effect of $`b_1>1`$ is to reduce the dependence of $`Q`$ on triangle shape, while $`b_2`$ introduces a constant offset. Since biasing operates in real space, instead of redshift space, equation (16) does not hold for redshift-space quantities; however, it turns out that the deviations from equation (16) are small (Scoccimarro, Couchman, & Frieman, 1999). The dot-dashed curve in Figure The Bispectrum of IRAS Redshift Catalogs shows the effect on $`Q`$ expected for $`1/b_1=1.25`$ and $`b_2/b_1^2=0.5`$.

We also compare our results below with previous determinations of the skewness $`S_3`$. One-point moments of the galaxy overdensity in a cell of volume $`V`$, $`\overline{\xi }_p=\delta _g^p_c`$, correspond to integrals over the $`p`$point correlation functions $`\xi _p`$,

$$\overline{\xi }_p=\frac{1}{V^p}_Vd^3x_1\mathrm{}d^3x_p\xi _p.$$
(17)
These moments also display hierarchical scaling to leading order in nonlinear perturbation theory,

$$\overline{\xi }_p=S_p\overline{\xi }^{p1},$$
(18)
where the parameters $`S_p`$ depend on the shape of the volume and on the spectral index. For a spherical tophat filter and power-law spectrum $`Pk^n`$, the skewness $`S_3`$ is (Fry & Gaztañaga, 1993)

$$S_3=\frac{1}{b_1}\left[\frac{34}{7}(3+n)\right]+\frac{3b_2}{b_1^2},$$
(19)
where the result in brackets corresponds to the skewness in the mass induced by gravity from Gaussian initial conditions (Juszkiewicz et al., 1993; Bernardeau, 1994).

The relation of equation (16) assumes the bias is local and deterministic. Locality here means that the scales we consider below, $`k0.2h\mathrm{Mpc}^1`$, are assumed to be much larger than the characteristic scale of galaxy formation, where nonlocal processes other than gravity become important. Deterministic bias assumes that the scatter in the relation between $`\delta _g`$ and $`\delta `$ can be neglected, with biasing described by the mean relation in equation (15). Deviations from deterministic bias can be introduced by adding a random field $`ϵ`$, whose variance describes the scatter about the mean, to the right side of equation (15). However, at least for models in which the scatter is local, that is, uncorrelated at large scales, $`ϵ(\text{x}_1)ϵ(\text{x}_2)\delta _\mathrm{D}(\text{x}_1\text{x}_2)`$, the use of equation (15) together with bispectrum measurements to constrain bias parameters results in a good recovery of the mean biasing relation (Paper I).

We note that the use of the bispectrum to recover information about galaxy bias does not depend exclusively on the validity of equation (16); other biasing schemes may predict other relationships between galaxy and dark matter bispectra that can similarly be tested. For example, in non-Gaussian models, bias and non-Gaussianity interact in a non-trivial way (Fry & Scherrer, 1994). We will consider an example of this class of models below when we discuss $`\chi ^2`$ initial conditions. Alternatively, one can consider models in which the bias is local in Lagrangian space (e.g. Mo & White (1996)) rather than Eulerian space as assumed in equation (15). In this case, the conditions for galaxy formation depend on the linearly extrapolated density field \[and thus equation (15) holds for the initial field configurations which are then evolved by gravity\], and the prediction for the reduced bispectrum in equation (16) changes to (Catelan et al., 1998)

$$Q_g=\frac{1}{b_1}Q_{123}+\frac{b_2}{b_1^2}+\left(\frac{b_11}{b_1^2}\right)\frac{\mathrm{\Delta }Q_{12}P_1P_2+\mathrm{\Delta }Q_{23}P_2P_3+\mathrm{\Delta }Q_{31}P_3P_1}{P_1P_2+P_2P_3+P_3P_1},$$
(20)
where $`\mathrm{\Delta }Q_{ij}=(k_i/k_j+k_j/k_i)\mathrm{cos}\theta _{ij}`$. The additional source of shape dependence in the third term can then be tested against observations. It is interesting to note that the difference between Eulerian and Lagrangian bias does not show up in large-scale measurements of the power spectrum; to leading (linear) order, both schemes agree. It is the sensitivity of the bispectrum to non-linear effects that makes possible to distinguish between these schemes.

## 3 Catalogs and Optimally Weighted Fourier Transform

We apply our analysis to three IRAS redshift subsamples: the updated QDOT survey (Efstathiou et al., 1990) and the 2 Jy (Strauss et al., 1992) and 1.2 Jy (Fisher et al., 1995) redshift surveys. The QDOT survey chooses at random one in six galaxies from the IRAS point source catalog with a $`60\mu \mathrm{m}`$ flux $`f_{60}>0.6\mathrm{Jy}`$. In this sample there are 1824 galaxies with galactic latitude $`|b|>10\mathrm{°}`$, with redshifts that correspond to distances $`20h^1\mathrm{Mpc}<R<500h^1\mathrm{Mpc}`$. We have used a revised version of the QDOT database in which a redshift error that afflicted approximately 200 southern galaxies has been corrected, and we have converted all redshifts to the local group frame. The other two samples are shallower but denser than the QDOT catalog. The 2 Jy catalog, complete to a flux limit $`f_{60}>2\mathrm{Jy}`$, contains 2072 galaxies; the 1.2 Jy catalog, with $`f_{60}>1.2\mathrm{Jy}`$, contains 4545 galaxies. Although all three catalogs are large-angle surveys, they are not full-sky: because of obscuration due to the Galactic plane, the existence of some bright sources, and other technical details, these surveys cover a solid angle of $`9.282\mathrm{sr}`$, about three quarters of the sky. In a subsequent paper (Paper III) we will analyze the IRAS PSCz survey (Saunders et.al.,, 2000).

To deal with the angular mask and the redshift distribution, we follow the formalism developed by FKP, computing moments of a weighted transform of the difference between the data and a synthetic random Poisson catalog with the same radial distribution and angular mask,

$$\stackrel{~}{F}(\text{k})=d^3xw(\text{x})[n_g(\text{x})\alpha n_s(\text{x})]e^{i\text{k}\text{x}},$$
(21)
where $`n_g(\text{x})=_i\delta _\mathrm{D}(\text{x}\text{x}_i)`$ is the point distribution for the galaxy data and $`n_s(\text{x})`$ similarly for the synthetic catalog. The factor $`\alpha =N_g/N_s`$ scales the synthetic catalog to the same density as the data, and $`w(\text{x})=1/[1+\overline{n}(\text{x})P_0]`$ is the “optimal” weight function for computing the power spectrum (and bispectrum as well, see Paper I); for detailed analysis of the formalism we follow, see FKP. The subtraction removes a contribution to the power from the shape of the radial selection function. The Fourier transform is performed precisely, summed over objects at their exact positions, rather than using interpolation and a fast Fourier transform, for frequencies that are uniformly spaced on a grid, $`\text{k}=(n_x,n_y,n_z)k_f`$, where $`k_f=0.005h\mathrm{Mpc}^1`$. The results in this paper are obtained by setting the fiducial power in the weight function to $`P_0=8000`$, 2000, and $`2000(h^1\mathrm{Mpc})^3`$ for QDOT, 2 Jy, and 1.2 Jy respectively. These values of $`P_0`$ reflect the different effective depths of the samples; in practice, the inferred value of $`Q`$ is not very sensitive to this parameter.

Including contributions from discreteness, the second moment of $`\stackrel{~}{F}`$ is related to the power spectrum of the galaxy distribution as

$$|\stackrel{~}{F}(\text{k})|^2=\frac{d^3k^{}}{(2\pi )^3}|G(\text{k}^{})|^2P(|\text{k}\text{k}^{}|)+(1+\alpha )d^3xw^2\overline{n},$$
(22)
where the window function $`G(\text{k})`$ is

$$G(\text{k})=d^3xw(\text{x})\overline{n}(\text{x})e^{i\text{k}\text{x}},$$
(23)
and $`\overline{n}(\text{x})=n_g(\text{x})`$ is the mean space density of galaxies, given the angular and luminosity selection criteria. As noted above, for anisotropic survey geometries the window function $`G(\text{k})`$ is broad in at least one dimension; as a result, a range of scales contribute to $`|\stackrel{~}{F}(\text{k})|^2`$. For the IRAS catalogs, however, $`G(\text{k})`$ is a rapidly falling function of $`|\text{k}|`$; in this case, for $`k`$ not too small, equation (22) becomes

$$|\stackrel{~}{F}(\text{k})|^2=P(k)d^3xw^2\overline{n}^2+(1+\alpha )d^3xw^2\overline{n}.$$
(24)
Similarly, the three-point function of $`\stackrel{~}{F}`$ becomes

$$\stackrel{~}{F}_1\stackrel{~}{F}_2\stackrel{~}{F}_3=B_{123}d^3xw^3\overline{n}^3+(P_1+P_2+P_3)d^3xw^3\overline{n}^2+(1\alpha ^2)d^3xw^3\overline{n}.$$
(25)
Subtracting the discreteness contributions, our estimators for the power spectrum and bispectrum are

$$\widehat{P}(k)=\frac{|\stackrel{~}{F}(\text{k})|^2}{d^3xw^2\overline{n}^2}(1+\alpha )\frac{d^3xw^2\overline{n}}{d^3xw^2\overline{n}^2},$$
(26)
$$\widehat{B}_{123}=\frac{\stackrel{~}{F}_1\stackrel{~}{F}_2\stackrel{~}{F}_3}{d^3xw^3\overline{n}^3}(\widehat{P}_1+\widehat{P}_2+\widehat{P}_3)\frac{d^3xw^3\overline{n}^2}{d^3xw^3\overline{n}^3}(1\alpha ^2)\frac{d^3xw^3\overline{n}}{d^3xw^3\overline{n}^3},$$
(27)
and the reduced bispectrum estimator is

$$\widehat{Q}=\frac{\widehat{B}_{123}}{\widehat{P}_1\widehat{P}_2+\widehat{P}_1\widehat{P}_3+\widehat{P}_2\widehat{P}_3}.$$
(28)
## 4 Data Analysis: Mock Catalogs

Having obtained the bispectrum amplitude $`Q`$ following the steps in the previous section, we are faced with interpreting the results and comparing to theoretical predictions such as those in Section 2. To do this, we must take into account a number of effects present in real redshift surveys that can potentially modify these predictions. In this section we summarize results from Paper I, where these effects are discussed in detail using mock catalogs constructed from perturbation theory and $`N`$-body simulations. We refer the reader to Paper I for a more in-depth analysis.

### 4.1 Nonperturbative Redshift Distortions

There are two ways in which the perturbation theory calculations of redshift space distortions described in Section 2 potentially need to be improved, namely, by going beyond both the plane-parallel and perturbative redshift mapping approximations. In the plane-parallel approximation assumed above, the structures being observed are considered to be sufficiently far away that the line of sight can be considered a fixed direction. Since the IRAS surveys cover most of the sky, this is not necessarily a good approximation. In this paper, we only consider the angle-averaged (monopole) power spectrum and bispectrum, which are very insensitive to the radial nature of the distortions (see Figs. 2 and 3 in Paper I). In our likelihood analysis, we nevertheless use mock catalogs in which the redshift-space mapping is correctly performed radially.

By working at large scales, $`k0.2h\mathrm{Mpc}^1`$, we can apply leading order nonlinear perturbation theory to describe the shape depence of the bispectrum in the absence of redshift distortions (Scoccimarro et al., 1998). However, in redshift space, another effect beyond nonlinear dynamics enters: the mapping from real to redshift space is itself is nonlinear. The discussion of redshift distortions in Section 2 assumed a perturbative expansion for this mapping as well as the dynamics. The validity of this expansion is in fact rather restricted (Scoccimarro, Couchman, & Frieman, 1999); over the range of scales we consider, nonlinear corrections to the redshift-space mapping are non-negligible (see Fig. 3 in Paper I). We take these effects into account by using second order Lagrangian perturbation theory (2LPT) numerical realizations in which the redshift-space mapping is done exactly.

### 4.2 Finite Volume Effects

The finite volume of a survey can alter the amplitude and scale or shape dependence of the bispectrum, potential systematic effects that we must consider before we can directly compare theory with observations. \[For a general discussion of the same issues in one-point statistics, see Szapudi, Colombi, & Bernardeau (1999) and references therein.\] Since we work at scales small compared to the extent of the survey, the integral constraint bias is negligible, thus there are basically two manifestations of finite volume effects: estimator bias and statistical bias.

Estimator bias, the fact that the estimator of the reduced bispectrum in equation (28) does not obey $`\widehat{Q}=Q`$, has two significant sources. Although the power spectrum and bispectrum estimators are themselves unbiased, $`\widehat{Q}`$ is a non-linear combination of them, thus an estimation bias arises (Hui & Gaztañaga, 1999; Szapudi, Colombi, & Bernardeau, 1999). In Paper I, we found that taking into account the estimator bias is important to meaningfully compare the bispectrum obtained in different surveys. However, this is not the full story. One can consider these estimators as random fields, each of which can be described by a probability distribution function (PDF) and its moments. Normally, one requires an estimator to be unbiased (a condition on its first moment) and of minimum variance (a requirement on the second moment). However, the higher order moments of estimators can also play a role. For example, if the skewness of the estimator PDF is positive (negative), even if the estimator is unbiased, in a given realization the most likely value of the estimator will be systematically lower (higher) than the mean. We call this statistical bias. Since in observations we deal with a single sample, this effect must be taken into account (Szapudi, & Colombi, 1996).

We must also account for correlations between estimators. In an infinite universe, Fourier modes are statistically independent. For a survey of finite volume, translation invariance is broken, leading to correlations between different Fourier modes at large scales. We confine our analysis to sufficiently small length scales, $`k0.05h\mathrm{Mpc}^1`$, that the window function of the survey can be approximated by a delta function (see Paper I for a detailed discussion). However, Fourier modes will be correlated with their neighbors in $`k`$-space over the inverse width of the selection function of the survey. In addition, correlations arising from shot noise dominate at small scales.

A proper analysis of these surveys reqires we take these effects into consideration. In Paper I we have studied this problem by using mock catalogs drawn from numerical realizations of perturbation theory as well as $`N`$-body simulations. In order to measure the correlation matrix of the bispectrum and the PDF of its estimator, many ($`400`$) realizations of the survey under consideration are necessary. This is efficiently done using a numerical implementation of second order Lagrangian perturbation theory (2LPT), which is orders of magnitude faster than $`N`$-body simulations and reproduces the correct bispectrum at large scales, including nonlinear aspects of redshift distortions mentioned above. In order to investigate the domain of validity of 2LPT, we have run N-body simulations which confirm that, over the range of scales we consider ($`0.05k0.2h\mathrm{Mpc}^1`$), 2LPT is a very good approximation to the full nonlinear theory.

### 4.3 Likelihood Analysis

The survey geometry and radial selection function introduce biases and correlations that are best dealt with in a likelihood analysis. We proceed in the following way: Given a set of measured reduced bispectrum amplitudes $`\{Q_m\}`$, $`m=1,\mathrm{},N_T`$, where $`N_T`$ is the number of closed triangles in the survey, we diagonalize their covariance matrix so that the $`Q`$-eigenmodes $`\widehat{q}_n`$,

$$\widehat{q}_n=\underset{m=1}{\overset{N_T}{}}\gamma _{mn}\frac{Q_m\overline{Q}_m}{\mathrm{\Delta }Q_m},$$
(29)
satisfy

$$\widehat{q}_n\widehat{q}_m=\lambda _n^2\delta _{nm},$$
(30)
where $`\overline{Q}Q`$ and $`(\mathrm{\Delta }Q)^2Q^2\overline{Q}^2`$. These $`Q`$-eigenmodes have “signal to noise” ratio

$$\left(\frac{S}{N}\right)_n\frac{1}{\lambda _n}\left|\underset{m=1}{\overset{N_T}{}}\gamma _{mn}\frac{\overline{Q}_m}{\mathrm{\Delta }Q_m}\right|.$$
(31)
The physical interpretation of the $`Q`$-eigenmodes becomes clear when ordered in terms of their signal to noise. The ($`n=1`$) eigenmode with highest signal to noise, $`S/N3`$, corresponds to all weights $`\gamma _{m1}>0`$; it represents the overall amplitude of the bispectrum. The $`n=2`$ $`Q`$-eigenmode, with $`S/N1`$, has $`\gamma _{m2}>0`$ for nearly flat triangles and $`\gamma _{m2}<0`$ for nearly equilateral triangles—it represents the configuration dependence of the bispectrum. Higher order ($`n>2`$) eigenmodes contain further information, such as variations of the amplitude and shape with scale; although they generally have $`S/N<1`$, there are many of them, so they cannot be neglected.

If the joint PDF of the $`Q_m`$ were Gaussian, then diagonalizing the covariance matrix would guarantee that the new eigenmodes are statistically independent—the full PDF would be the product of their individual PDF’s. Although this is not necessarily so in the general non-Gaussian case, we will assume that this independence is a good approximation here (supported by the results of Paper I). We can then write down the likelihood as a function of the parameters,

$$(\alpha _1,\alpha _2)\underset{i=1}{\overset{N_T}{}}\text{P}_i[\nu _i(\alpha _1,\alpha _2)],$$
(32)
where $`\alpha _11/b_1`$, $`\alpha _2b_2/b_1^2`$, and

$$\nu _i(\alpha _1,\alpha _2)\frac{1}{\lambda _i}\underset{j=1}{\overset{N_T}{}}\gamma _{ji}\frac{Q_j(\alpha _1\overline{Q}_j+\alpha _2)}{\mathrm{\Delta }Q_j}.$$
(33)
Here the $`Q_j`$ ($`j=1,\mathrm{},N_T`$) are the data, and the mean $`\overline{Q}_j`$, standard deviation $`\mathrm{\Delta }Q_j`$, and the non-Gaussian PDF’s $`\text{P}_i(\nu _i)`$ are extracted from the mock catalogs. As described in Section 2, the reduced bispectrum is mostly sensitive to the power spectrum shape, with a very small dependence on $`\mathrm{\Omega }_m`$ through redshift distortions and negligible dependence on the normalization $`\sigma _8`$ at large scales. Our mock catalogs are based on the $`\mathrm{\Lambda }`$CDM model, with $`\mathrm{\Omega }_m=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, and $`\sigma _8=0.7`$; in redshift space, for this model $`\sigma _8=0.84`$, in agreement with the observed normalization for IRAS galaxies (Fisher et al., 1994). The power spectrum shape of this model is also in good agreement with the observed IRAS galaxy power spectrum over the range of scales we consider (Feldman, Kaiser, & Peacock, 1994; Fisher et al., 1993). If we change the power spectrum shape, the resulting bias parameters do change, as the prediction for $`Q`$ is sensitive to spectral index; for example, for a shape parameter $`\mathrm{\Gamma }=0.5`$, as in the SCDM model, $`Q`$ is related to that in the $`\mathrm{\Lambda }`$CDM model ($`\mathrm{\Gamma }=0.21`$) by approximately $`Q(\mathrm{\Gamma }=0.21)1.3Q(\mathrm{\Gamma }=0.5)0.05`$ (Fig. 5, Paper I), so that $`b_1^{SCDM}b_1^{\mathrm{\Lambda }CDM}/1.3`$ and $`(b_2/b_1^2)^{SCDM}(b_2/b_1^2)^{\mathrm{\Lambda }CDM}0.05/b_1^{SCDM}`$; thus, the inferred bias parameters for IRAS galaxies would be smaller if we assumed SCDM rather than the $`\mathrm{\Lambda }`$CDM model. To assess the sensitivity of our results to a change in $`\mathrm{\Omega }_m`$, we also consider mock catalogs for a model with $`\mathrm{\Omega }_m=1`$ but with the same power spectrum normalization and shape as the $`\mathrm{\Lambda }`$CDM model above. We do not consider mock catalogs with different values of $`\sigma _8`$, since the dependence of $`Q`$ on this parameter is very weak (see Paper I). The main effect of a change in $`\sigma _8`$ would be to decrease (increase) the error bars we infer if IRAS galaxies have linear bias parameter $`b_1<1`$ ($`b_1>1`$), due to shot noise reduction (increase). In this regard, the error bars we derive on the bias parameters can be considered conservative.

Important results obtained from results of mock catalogs as described in Paper I include

* Estimator bias can be as large as a factor of two. However, equation (16) continues to hold for biased estimates.
* The PDF of the bispectrum estimator is generally skewed and has non-Gaussian exponential tails or even power-law tails in the case of sparse sampling (QDOT) and $`\chi ^2`$ initial conditions.
* The cross-correlation coefficient induced by the window function of the survey and shot-noise between different power spectrum bins or different triangles is not significantly less than unity.

A maximum likelihood method to constrain bias parameters from redshift surveys has been proposed before by Matarrese, Verde, & Heavens (1997). However, their method does not take into account these effects and thus cannot be applied as it stands to IRAS surveys.

## 5 Results

The bispectrum is defined for closed triangles in Fourier space. After averaging over orientations, we can characterize a given triangle by the lengths of its three sides, $`k_1`$, $`k_2`$, $`k_3`$, where we order $`k_1k_2k_3`$, or by any combination of parameters that completely defines a triangle. Figures The Bispectrum of IRAS Redshift Catalogs, The Bispectrum of IRAS Redshift Catalogs, and The Bispectrum of IRAS Redshift Catalogs show the scale dependence of $`Q(k_1)`$ for the QDOT, 2 Jy, and 1.2 Jy surveys, for triangles in which the two shorter sides $`k_2`$ and $`k_3`$ are separated by angle $`\theta _{23}`$, i.e., $`\widehat{\text{k}}_2\widehat{\text{k}}_3=\mathrm{cos}\theta _{23}`$. The four panels show results for four ranges of the angle $`\theta _{23}`$; small $`\theta _{23}`$ corresponds to nearly flat triangles, i.e., all three sides almost colinear, while triangles with larger $`\theta _{23}`$ are more open. The results at very small $`k`$ are sensitive to the subtraction of power in the selection function, so we have constrained the smallest side of the triangle to have $`k_30.05h\mathrm{Mpc}^1`$. The scatter in the values of $`Q`$ is not small, but it is still apparent that there is a nonzero positive signal. Neighboring points are clearly correlated; for this reason we have not included error bars that describe the dispersion of values that enter into each point. A full treatment of triangle correlations, as described in § 4, is carried out in the analysis below. Dashed lines show the median and 68% width for the points in each panel. The average value of $`Q`$ in each window generally decreases as $`\theta _{23}`$ increases, as expected from gravitational instability, even for QDOT, where the scatter is largest.

Figures The Bispectrum of IRAS Redshift Catalogs, The Bispectrum of IRAS Redshift Catalogs, and The Bispectrum of IRAS Redshift Catalogs show the shape dependence $`Q(\theta )`$ for triangles with sides with two sides of ratio 0.4–0.6 separated by angle $`\theta `$. The solid curves show the tree-level (leading order) prediction from gravitational instability in perturbation theory with a spectral index $`n=1.4`$ (as in FKP), including perturbative redshift distortions (eq. ). The long-dashed curves show $`Q`$ averaged over a large number of synthetic 2LPT catalogs, and includes nonperturbative effects. Different symbols show results for $`k_1=0.1`$$`0.125h\mathrm{Mpc}^1`$ (triangles), $`k_1=0.125`$$`0.15h\mathrm{Mpc}^1`$ (squares), $`k_1=0.15`$$`0.175h\mathrm{Mpc}^1`$ (diamonds), $`k_1=0.175`$$`0.2h\mathrm{Mpc}^1`$ (circles). The signal is weak in the sparse and noisy QDOT survey, but in both the 2 Jy and 1.2 Jy surveys a shape dependence is apparent, although there is significant dispersion, especially at nearly colinear configurations. For 2 Jy and 1.2 Jy, open symbols show the direct result, and filled symbols are corrected for the average finite volume bias. The short-dashed curves show the local bias model of equation (16) using the results of the likelihood analysis below, including non-linear redshift distortions, estimator bias, correlations between triangles, and the non-Gaussianity of the bispectrum PDF, applied to the 2LPT result. While the finite sample correction removes the bias in the mean, the distribution is still skewed to low values, and so the filled symbols cannot be directly compared to theoretical predictions. It is visually apparent that the skewness is larger in the 2 Jy catalog, consistent with the findings in Paper I.

To assess to what extent the data follow the shape dependence predicted by gravitational instability and to extract quantitative bias parameters, we have carried out the likelihood procedure described in Section 4. There are potentially 4741 closed triangles binned with integer sides $`k_i/k_f=10,\mathrm{},40`$ (although the galaxy overdensity is not mapped into a grid, we use wavevectors integer multiples of $`k_f=0.005h\mathrm{Mpc}^1`$). Since the window function of the survey has a width of order $`\mathrm{\Delta }k/k_f5`$, the number of “independent” triangles is much smaller than 4741. We therefore use coarser bins, defined in terms of shape $`s`$, ratio $`r`$, and scale $`k`$ parameters (Peebles, 1980),

$$s\frac{k_1k_2}{k_3},r\frac{k_2}{k_3},kk_3,$$
(34)
where the “shape” parameter $`s`$ obeys $`0s1`$ ($`k_1k_2k_3`$), the ratio $`1r4`$, and the overall scale of the triangle satisfies $`10k/k_f40`$. The shape parameter is $`s=0`$ for isosceles triangles (and if $`r=1`$ these are equilateral triangles); for $`s=1`$ we have colinear triangles. Using 10 bins in each variable plus the closed triangle constraint yields 203 triangles. We find that using a finer binning leads to a more unstable behavior of numerical computation of the $`Q`$ eigenvalues and eigenvectors, while making the binning too coarse loses information on the shape and scale dependence of the bispectrum.

Figure The Bispectrum of IRAS Redshift Catalogs shows the results of the likelihood analysis for the 2 Jy and 1.2 Jy catalogs. Due to the sparse sampling (and thus large shot noise errors), the bias parameters cannot be significantly constrained from the QDOT survey, so we do not show results for it. The bottom panel in Figure The Bispectrum of IRAS Redshift Catalogs shows 68% and 90% likelihood contours as a function of $`1/b_1`$ and $`b_2/b_1^2`$, assuming Gaussian initial conditions and matter density $`\mathrm{\Omega }_m=0.3`$. The maximum likelihood analysis yields the results in Table 1, with error bars obtained from the 68% marginalized likelihoods shown in the upper panels of Figure The Bispectrum of IRAS Redshift Catalogs. Despite the apparently different amplitudes of $`Q`$ in Figures The Bispectrum of IRAS Redshift Catalogs and The Bispectrum of IRAS Redshift Catalogs, the bias parameters obtained for the two catalogs are remarkably similar. The differences between Figures The Bispectrum of IRAS Redshift Catalogs and The Bispectrum of IRAS Redshift Catalogs are largely due to different finite volume effects in the two surveys, which trace to their different selection functions. In particular, the shallower 2 Jy catalog is more strongly affected by finite volume corrections and skewness of the bispectrum PDF than the 1.2 Jy survey.

Figure The Bispectrum of IRAS Redshift Catalogs shows the effect of changing the matter density from $`\mathrm{\Omega }_m=0.3`$ to $`\mathrm{\Omega }_m=1`$ in the case of the 1.2Jy catalog (similar results hold for the 2 Jy catalog). As expected from the analytic results described in Section 3 and illustrated in Figure The Bispectrum of IRAS Redshift Catalogs, the results given in Table 1 change by about $`10\%`$ in the direction consistent with Figure The Bispectrum of IRAS Redshift Catalogs: for $`\mathrm{\Omega }_m=1`$ the expected configuration dependence is slightly stronger, and thus the bias must be slightly larger to account for the same measured value. Note however that the $`\mathrm{\Omega }_m`$ dependence is not nearly as strong as in redshift-space power spectrum measurements which determine the combination $`\beta \mathrm{\Omega }^{0.6}/b_1`$; to retain the same value of $`\beta `$, our result would have to come down to $`1/b_1=0.56`$ for $`\mathrm{\Omega }_m=1`$. Thus, to a large extent, measurements of the bispectrum can be used to lift the degeneracy present in the power spectrum.

We found it useful to track how the results change as we vary the number of eigenmodes kept in the analysis. For the two highest $`S/N`$ eigenmodes, which represent the overall amplitude and configuration dependence of the bispectrum, the maximum likelihood values for the bias parameters are very similar to those in Table 1, but with larger error bars. As we increase the number of eigenmodes from 2 to 203, the likelihood contours slowly contract (these additional eigenmodes all have signal to noise less than unity in our fiducial model); the maximum likelihood bias values for the 2 Jy catalog hardly change, whereas in the 1.2 Jy survey the maximum likelihood slowly oscillates in the linear bias direction within the range $`1/b_1=1.0`$–1.3.

Since these additional eigenmodes carry information on the scale dependence of the reduced bispectrum, they can help to constrain non-Gaussian initial conditions, as exemplified by the $`\chi ^2`$ model (Linde & Mukhanov, 1997; Antoniadis et al., 1997; Peebles, 1997, 1999a, 1999b), in which the linear density field is the square of a Gaussian field. As a result, the initial bispectrum is non-zero and has a different scaling, $`B_iP_i^{3/2}`$, and shape dependence from the gravitionally induced bispectrum in equation (4). However, in order to compare with observations, linear extrapolation of the initial conditions is insufficient, except at extremely large scales or high redshift, which are not probed by the IRAS surveys. Taking into account gravitational corrections, one finds that the overall scale-independence of $`Q`$ is partially restored by gravity over the scales we probe; however, the amplitude of $`Q`$ is roughly two times larger than in the Gaussian case. The configuration dependence resembles that of the Gaussian case, but with residual dependences on the ratios $`k_i/k_j`$, which would show up as a rather strong dependence on the variable $`r`$ in equation (34) (Scoccimarro, 2000a). In order to test this model against observations, we must take into account the effects of biasing as well, which are non-trivial for general non-Gaussian models. For the particular case of $`\chi ^2`$ initial conditions, however, a simple result follows,

$$Q_g=\frac{1}{b_1}Q+\frac{b_2^{\mathrm{eff}}}{b_1^2},$$
(35)
where $`b_2^{\mathrm{eff}}(38/\sqrt{3}\pi )b_21.53b_2`$ (Scoccimarro, 2000a). As we did in the Gaussian case, we shall assume that the same relation approximately holds in redshift-space as well. We have run 2LPT realizations and created mock 1.2 Jy catalogs for $`\chi ^2`$ initial conditions, as described in Paper I. We apply the same likelihood analysis as in the Gaussian case, with appropriate change for the mean, variance, and PDF of the bispectra in equation (32). The results of this analysis are shown in Figure The Bispectrum of IRAS Redshift Catalogs. If we consider only the two highest $`S/N`$ eigenmodes, which probe the overall amplitude and configuration dependence of the bispectrum, we obtain the $`68\%`$ contour shown in the bottom panel, which gives $`1/b_1=1.42_{0.49}^{+0.34}`$ and $`b_2/b_1^2=1.40\pm 0.58`$, consistent with the fact that the overall amplitude of $`Q`$ is larger than in the Gaussian case, and the configuration dependence is about the same. However, including the remaining eigenmodes, which are sensitive to $`r`$ and $`k`$ dependence, the maximum likelihood is driven to $`1/b_1=0`$, in an attempt to cancel the $`r`$ and $`k`$ dependence in the $`\chi ^2`$ model not seen in the data. If we consider that IRAS galaxies are not biased at large scales ($`1/b_11`$), then the 1.2 Jy bispectrum results rule out $`\chi ^2`$ initial conditions at the 95% confidence level. We note that this model may be considered ‘strongly’ non-Gaussian; models with substantially ‘weaker’ initial non-Gaussianity (e.g., with negligibly small initial three- and four-point correlations) are not excluded by this analysis.

We can also use this data to constrain the Lagrangian local bias model. This model also introduces additional dependence of the bispectrum on $`r`$, through the last term in equation (20), although with a much smaller amplitude than in the $`\chi ^2`$ model. We can therefore try to test whether the local Lagrangian or the Eulerian model is preferred. We have repeated the likelihood analysis using equation (20) in place of equation (16) to incorporate the change in the mean of $`Q`$ as a function of biasing parameters. We find that the goodness of fit, given by the value of the likelihood function at its maximum, does not change from one model to the other by more than 50%. The largest change is obtained for the 2 Jy survey, shown in Figure The Bispectrum of IRAS Redshift Catalogs , where we compare the Eulerian result (same as in Figure The Bispectrum of IRAS Redshift Catalogs) to the Lagrangian case (for the $`\mathrm{\Lambda }`$CDM parameters). A similar behavior as in the $`\chi ^2`$ model is observed: the contours shift towards higher bias, as listed in Table 1. These values of the bias parameters are consistent with those obtained in the Eulerian model. A larger dataset is required to distinguish between these two biasing schemes.

## 6 Conclusions

We have presented results for the reduced bispectrum $`Q`$ of three IRAS redshift catalogs, computed using the optimal weighting scheme of FKP, after subtracting contributions to the power from discreteness and from the shape of the selection function. That $`Q`$ is nonzero, positive, and generally of the shape predicted by gravitational instability is unmistakable. The ability to extract values for bias parameters is marginal; however, after taking into account non-linear redshift distortions and finite volume effects, likelihood analysis yields reasonable values for the parameters of a local Eulerian bias model, $`1/b_1=1.32_{0.58}^{+0.36}`$, $`1.15_{0.39}^{+0.39}`$ and $`b_2/b_1^2=0.57_{0.30}^{+0.45}`$, $`0.50_{0.51}^{+0.31}`$, for 2 Jy and 1.2 Jy respectively. This implies that IRAS galaxies trace dark matter with decreasing effective bias as a function of the mass density fluctuation amplitude, consistent with their being under-represented in the cores of clusters. We have also considered constraints on non-Gaussian initial conditions, using as an example the $`\chi ^2`$ model, in which the bispectrum develops a different amplitude and scale dependence than in the Gaussian case. If IRAS galaxies are not positively biased at large scales (i.e., if $`b_1^{IRAS}1`$), then $`\chi ^2`$ initial conditions are ruled out at the 95% confidence level.

The likelihood contours show that one combination of the bias parameters is known to better accuracy: the long axis of the likelihood contours in Figures The Bispectrum of IRAS Redshift CatalogsThe Bispectrum of IRAS Redshift Catalogs follows tracks of roughly constant skewness, $`S_3=3.25/b_1+3b_2/b^2`$ for $`n=1.4`$. Our values inferred for $`S_3`$ are consistent with previous work on the skewness of IRAS galaxies. For the 1.2 Jy sample, Bouchet et al. (1993) obtained $`S_3=1.5\pm 0.5`$; for the same catalog, Kim & Strauss (1988), who include a finite volume correction but use a different estimator, obtained $`S_3=2.83\pm 0.09`$. Our results translate via equation (19) (assuming $`n=1.4`$) into a finite-volume corrected value $`S_32.25\pm 0.9`$, comparable to that of Kim & Strauss (1988). To estimate the error bar in the skewness we assume that $`S_3`$ corresponds to sum of $`Q`$ over all triangles. The error is then roughly given by three times the width in the best determined direction in the ellipse in Figure The Bispectrum of IRAS Redshift Catalogs, $`\pm 0.3`$. An approximate estimate of $`S_3`$ that does not include the finite volume corrections can be obtained directly from Figure The Bispectrum of IRAS Redshift Catalogs, by using the fact that at the $`68\%`$ limit, $`Q/b_1+b_2/b_1^2`$ is in the range 0.2–0.8; this implies $`S_33Q_g1.5\pm 0.9`$, consistent with Bouchet et al. (1993). These are only rough estimates, as they ignore the skewness of the bispectrum (and $`S_3`$) PDF.

Our results have interesting implications for the relative clustering of IRAS and optically selected galaxies. From the ratio of power spectrum amplitudes between the Stromlo-APM survey and the combined QDOT-1.2Jy survey, Tadros & Efstathiou (1996) conclude that $`b_1^{opt}/b_1^{IRAS}1.2`$, while from the analysis of the APM projected three-point function, Frieman & Gaztañaga (1999) conclude that $`b_1^{opt}1`$. These results are consistent with our findings for $`b_1^{IRAS}`$. Future galaxy surveys will be able to quantify this issue with a greater level of accuracy.

We thank Stephane Colombi, Enrique Gaztañaga, Lam Hui, Simon Prunet, and Matias Zaldarriaga for useful discussions. Research supported in part by the NSF EPSCoR Grant at Kansas and the University of Kansas GRF, by the DOE and NASA grant NAG 5-7092 at Fermilab, and by NASA NAG5-2835 at the University of Florida. Parts of this work took place at the Aspen Center for Physics.