text:
# Statistical Cosmology with Quadratic Density Fields

## 1 Introduction

In most popular versions of the gravitational instability model for the origin of cosmic structure, particularly those involving cosmic inflation (Guth 1981; Guth & Pi 1982), the initial fluctuations that seeded the structure formation process form a Gaussian random field (Adler 1981; Bardeen et al. 1986).

Gaussian random fields are the simplest fully-defined stochastic processes, which makes analysis of them relatively straightforward. Robust and powerful statistical descriptors can be constructed that have a firm mathematical underpinning and are relatively simple to implement. Second-order statistics such as the ubiquitous power-spectrum (e.g. Peacock & Dodds 1996) furnish a complete description of Gaussian fields. They have consequently yielded invaluable insights into the behaviour of large-scale structure in the latest generation of redshift surveys, such as the 2dFGRS (Percival et al. 2001). Important though these methods undoubtedly are, the era of precision cosmology we are now entering requires more thought to be given to methods for both detecting and exploiting departures from Gaussian behaviour.

The pressing need for statistics appropriate to the analysis of non-linear stochastic processes also suggests a need to revisit some of the fundamental properties cosmologists usually assume when studying samples of the Universe. Gaussian random fields have many useful properties. It is straightforward to impose constraints that result in statistically homogeneous fields, for example. Perhaps more relevantly one can understand the conditions under which averages over a single spatial domain are well-defined, the constraint of sample-homogeneity. The conditions under which such fields can be ergodic are also well established. It is known that smoothing Gaussian fields preserves Gaussianity, and so on. These properties are all somewhat related, but not identical. Indeed, as we shall see in Section 4, looking at the corresponding properties of non-linear fields turns up some interesting results and delivers warnings to be careful. Exploring these properties is the first aim of this paper.

Even if the primordial density fluctuations were indeed Gaussian, the later stages of gravitational clustering must induce some form of non-linearity. One particular way of looking at this issue is to study the behaviour of Fourier modes of the cosmological density field. If the hypothesis of primordial Gaussianity is correct then these modes began with random spatial phases. In the early stages of evolution, the plane-wave components of the density evolve independently like linear waves on the surface of deep water. As the structures grow in mass, they interact with other in non-linear ways, more like waves breaking in shallow water. These mode-mode interactions lead to the generation of coupled phases. While the Fourier phases of a Gaussian field contain no information (they are random), non-linearity generates non-random phases that contain much information about the spatial pattern of the fluctuations. Although the significance of phase information in cosmology is still not fully understood, there have been a number of attempts to gain quantitative insight into the behaviour of phases in gravitational systems. Ryden & Gramann (1991), Soda & Suto (1992) and Jain & Bertschinger (1998) concentrated on the evolution of phase shifts for individual modes using perturbation theory and numerical simulations. An alternative approach was adopted by Scherrer, Melott & Shandarin (1991), who developed a practical method for measuring the phase coupling in random fields that could be applied to real data. Most recently Chiang & Coles (2000), Coles & Chiang (2000), Chiang (2001) and Chiang, Naselsky & Coles (2002) have explored the evolution of phase information in some detail.

Despite this recent progress, there is still no clear understanding of how the behaviour of the Fourier phases manifests itself in more orthodox statistical descriptors. In particular there is much interest in the usefulness of the simplest possible generalisation of the (second-order) power-spectrum, i.e. the (third-order) bispectrum (Peebles 1980; Scoccimarro et al. 1998; Scoccimarro, Couchman & Frieman 1999; Verde et al. 2000; Verde et al. 2001; Verde et al. 2002). Since the bispectrum is identically zero for a Gaussian random field, it is generally accepted that the bispectrum encodes some form of phase information but it has never been elucidated exactly what form of correlation it measures. Further possible generalisations of the bispectrum are usually called polyspectra; they include the (fourth-order) trispectrum (Verde & Heavens 2001) or a related but simpler statistic called the second-spectrum (Stirling & Peacock 1996). Exploring the connection between polyspectra and non-linearly induced phase association is the second aim of this paper.

We investigate both these issues by reference to a particular form of non-Gaussian field that serves both as a kind of phenomenological paradigm and as a reasonably realistic model of non-linear evolution from Gaussian initial conditions. The model involves a field which is generated by a simple quadratic transformation of a Gaussian distribution, hence the term quadratic non-linearity. Quadratic fields have been discussed before from a number of contexts (e.g. Coles & Barrow 1987; Moscardini et al. 1991; Falk, Rangarajan & Srednicki 1993; Luo & Schramm 1993; Luo 1994; Gangui et al. 1994; Koyoma, Soda & Taruya 1999; Peebles 1999a,b; Matarrese, Verde & Jimenez 2000; Verde et al. 2000; Verde et al. 2001; Komatsu & Spergel 2001; Shandarin 2002; Bartolo, Matarrese & Riotto 2002); for further discussion see Section 3. Our motivation is therefore very similar to that of Coles & Jones (1991), which introduced the lognormal density field as an illustration of some of the consequences of a more extreme form of non-linearity involving an exponential transformation of the linear density field.

The plan is as follows. In the following section we introduce some fundamental concepts underlying statistical cosmology, more-or-less from first principles. We do this in order to allow the reader to see explicitly what assumptions underlie standard statistical practise. In Section 3 we then look at some of the contexts in which quadratic non-linearity may arise, either primordially or during the non-linear growth of structure from Gaussian fields. In Section 4 we revisit some of the basic properties used in Section 2 from the viewpoint of quadratic non-linearity and show how some basic implicit assumptions are violated. We then, in Section 5, explore how phase correlations arise in quadratic fields and relate these to higher-order statistics of quadratic fields. We discuss the lessons learned from this study in Section 6.

## 2 Basic Statistical Concepts

We start by giving some general definitions of concepts which we will later use in relation to the particular case of cosmological density fields. In order to put our results in a clear context, we develop the basic statistical description of cosmological density fields; see also, e.g., Peebles (1980).

### 2.1 Fourier Description

We follow standard practice and consider a region of the Universe having volume $`V_u`$, for convenience assumed to be a cube of side $`Ll_s`$, where $`l_s`$ is the maximum scale at which there is significant structure due to the perturbations. The region $`V_u`$ can be thought of as a “fair sample” of the Universe if this is the case. It is possible to construct, formally, a “realisation” of the Universe by dividing it into cells of volume $`V_u`$ with periodic boundary conditions at the faces of each cube. This device is often convenient, but in any case one often takes the limit $`V_u\mathrm{}`$. Let us denote by $`\overline{\rho }`$ the mean density in a volume $`V_u`$ and take $`\rho (𝐱)`$ to be the density at a point in this region specified by the position vector $`𝐱`$ with respect to some arbitrary origin. As usual, we define the fluctuation to be

$$\delta (𝐱)=[\rho (𝐱)\overline{\rho }]/\overline{\rho }.$$
(1)
We assume this to be expressible as a Fourier series:

$$\delta (𝐱)=\underset{𝐤}{}\delta _𝐤\mathrm{exp}(i𝐤𝐱)=\underset{𝐤}{}\delta _𝐤^{}\mathrm{exp}(i𝐤𝐱);\delta _𝐤=\frac{1}{V_u}_{V_u}\delta (𝐱)\mathrm{exp}(i𝐤𝐱)𝑑𝐱.$$
(2)
The Fourier coefficients $`\delta _𝐤`$ are complex quantities, $`\delta _𝐤=|\delta _𝐤|\mathrm{exp}(i\varphi _𝐤)`$ with amplitude $`|\delta _𝐤|`$ and phase $`\varphi _𝐤`$. Conservation of mass in $`V_u`$ implies $`\delta _{𝐤=0}=0`$ and the reality of $`\delta (𝐱)`$ requires $`\delta _𝐤^{}=\delta _𝐤`$.

If, instead of the volume $`V_u`$, we had chosen a different volume $`V_u^{}`$ the perturbation within the new volume would again be represented by a series of the form (2), but with different coefficients $`\delta _𝐤`$. Now consider a (large) number $`N`$ of realisations of our periodic volume and label these realisations by $`V_{u1}`$, $`V_{u2}`$, $`V_{u3}`$, …, $`V_{uN}`$. It is meaningful to consider the probability distribution $`𝒫(\delta _𝐤)`$ of the relevant coefficients $`\delta _𝐤`$ from realisation to realisation across this ensemble. One typically assumes that the distribution is statistically homogeneous and isotropic, in order to satisfy the Cosmological Principle, and that the real and imaginary parts of $`\delta _𝐤`$ have a Gaussian distribution and are mutually independent, so that

$$𝒫(w)=\frac{V_u^{1/2}}{(2\pi \alpha _k^2)^{1/2}}\mathrm{exp}\left(\frac{w^2V_u}{2\alpha _k^2}\right),$$
(3)
where $`w`$ stands for either $`\mathrm{Re}[\delta _𝐤]`$ or $`\mathrm{Im}[\delta _𝐤]`$ and $`\alpha _k^2=\sigma _k^2/2`$; $`\sigma _k^2`$ is the spectrum. This is the same as the assumption that the phases $`\varphi _𝐤`$ are mutually independent and randomly distributed over the interval between $`\varphi =0`$ and $`\varphi =2\pi `$. In this case the moduli of the Fourier amplitudes have a Rayleigh distribution:

$$𝒫(|\delta _𝐤|,\varphi _𝐤)d|\delta _𝐤|d\varphi _𝐤=\frac{|\delta _𝐤|V_u}{2\pi \sigma _k^2}\mathrm{exp}\left(\frac{|\delta _𝐤|^2V_u}{2\sigma _k^2}\right)d|\delta _𝐤|d\varphi _𝐤.$$
(4)
Because of the assumption of statistical homogeneity and isotropy, the quantity $`𝒫(\delta _𝐤)`$ depends only on the modulus of the wavevector $`𝐤`$ and not on its direction. It is fairly simple to show that, if the Fourier quantities $`|\delta _𝐤|`$ have the Rayleigh distribution, then the probability distribution $`𝒫(\delta )`$ of $`\delta =\delta (𝐱)`$ in real space is Gaussian, so that:

$$𝒫(\delta )d\delta =\frac{1}{(2\pi \sigma ^2)^{1/2}}\mathrm{exp}\left(\frac{\delta ^2}{2\sigma ^2}\right)d\delta ,$$
(5)
where $`\sigma ^2`$ is the variance of the density field $`\delta (𝐱)`$. This is a strict definition of Gaussianity. However, Gaussian statistics do not always require the distribution (4) for the Fourier component amplitudes. According to its Fourier expansion, $`\delta (𝐱)`$ is simply a sum over a large number of Fourier modes whose amplitudes are drawn from some distribution. If the phases of each of these modes are random, then the Central Limit Theorem will guarantee that the resulting superposition will be close to a Gaussian if the number of modes is large and the distribution of amplitudes has finite variance. Such fields are called weakly Gaussian.

### 2.2 Covariance Functions & Probability Densities

We now discuss the real-space statistical properties of spatial perturbations in $`\rho `$. The *covariance function* is defined in terms of the density fluctuation by

$$\xi (𝐫)=\frac{[\rho (𝐱)\overline{\rho }][\rho (𝐱+𝐫)\overline{\rho }]}{\overline{\rho }^2}=\delta (𝐱)\delta (𝐱+𝐫).$$
(6)
The angle brackets in this expression indicate two levels of averaging: first a volume average over a representative patch of the universe and second an average over different patches within the ensemble, in the manner of section section 2.1. Applying the Fourier machinery to equation (6) one arrives at the Wiener–Khintchin theorem, relating the covariance to the spectral density function or power spectrum, $`P(k)`$:

$$\xi (𝐫)=\underset{𝐤}{}|\delta _𝐤|^2\mathrm{exp}(i𝐤𝐫),$$
(7)
which, in passing to the limit $`V_u\mathrm{}`$, becomes

$$\xi (𝐫)=\frac{1}{(2\pi )^3}P(k)\mathrm{exp}(i𝐤𝐫)𝑑𝐤.$$
(8)
Averaging equation (7) over $`𝐫`$ gives

$$\xi (𝐫)_𝐫=\frac{1}{V_u}\underset{𝐤}{}|\delta _𝐤|^2\mathrm{exp}(i𝐤𝐫)𝑑𝐫=0.$$
(9)
The function $`\xi (𝐫)`$ is the *two–point* covariance function. In an analogous manner it is possible to define spatial covariance functions for $`N>2`$ points. For example, the three–point covariance function is

$$\zeta (r,s)=\frac{[\rho (𝐱)\overline{\rho }][\rho (𝐱+𝐫)\overline{\rho }][\rho (𝐱+𝐬)\overline{\rho }]}{\overline{\rho }^3}$$
(10)
which gives

$$\zeta (𝐫,𝐬)=\delta (𝐱)\delta (𝐱+𝐫)\delta (𝐱+𝐬),$$
(11)
where the spatial average is taken over all the points $`𝐱`$ and over all directions of $`𝐫`$ and $`𝐬`$ such that $`|𝐫𝐬|=t`$: in other words, over all points defining a triangle with sides $`r`$, $`s`$ and $`t`$. The generalisation of (10) to $`N>3`$ is obvious.

The covariance functions are related to the moments of the probability distributions of $`\delta (𝐱)`$. If the fluctuations form a Gaussian random field then the N-variate distributions of the set $`\delta _i\delta (𝐱_i)`$ are just multivariate Gaussians of the form

$$𝒫_N(\delta _1,\mathrm{},\delta _N)=\frac{1}{(2\pi )^{N/2}(\mathrm{det}C)^{1/2}}\mathrm{exp}\left(\frac{1}{2}\underset{i,j}{}\delta _iC_{ij}^1\delta _j\right).$$
(12)
The correlation matrix $`C_{ij}`$ can be expressed in terms of the covariance function

$$C_{ij}=\delta _i\delta _j=\xi (𝐫_{ij}).$$
(13)
It is convenient to go a stage further and define the N-point connected covariance functions as the part of the average $`\delta _i\mathrm{}\delta _N`$ that is not expressible in terms of lower order functions e.g.

$$\delta _1\delta _2\delta _3=\delta _1_c\delta _2\delta _3_c+\delta _2_c\delta _1\delta _3_c+\delta _3_c\delta _1\delta _2_c+\delta _1_c\delta _2_c\delta _3_c+\delta _1\delta _2\delta _3_c,$$
(14)
where the connected parts are $`\delta _1\delta _2\delta _3_c`$, $`\delta _1\delta _2_c`$, etc. Since $`\delta =0`$ by construction, $`\delta _1_c=\delta _1=0`$. Moreover, $`\delta _1\delta _2_c=\delta _1\delta _2`$ and $`\delta _1\delta _2\delta _3_c=\delta _1\delta _2\delta _3`$. The second and third order connected parts are simply the same as the covariance functions. Fourth and higher order quantities are different, however. The connected functions are just the multivariate generalisation of the cumulants $`\kappa _N`$ (Kendall & Stewart 1977). One of the most important properties of Gaussian fields is that all of their N-point connected covariances are zero beyond N=2, so that their statistical properties are fixed once the set of two–point covariances (13) is determined. All large-scale statistical properties are therefore determined by the asymptotic behaviour of $`\xi (r)`$ as $`r\mathrm{}`$. This simplifying property is not shared by non–Gaussian fields, a fact we shall explore in Section 4.

## 3 Quadratic Non-Linearity

In this section we discuss some of the circumstances wherein quadratic non-linearity may arise. At the outset we should admit that this study is primarily phenomenological and our intention is largely to use this as a model that displays some consequences of non-linearity.

### 3.1 Phenomenology

As far as we are aware, the first application of quadratic density fields in a cosmological setting was in Coles & Barrow (1987) who were studying the possible sample properties of non-Gaussian temperature fluctuations on the cosmic microwave background sky. They in fact studied a series of models called $`\chi _n^2`$ models obtained via the transformation

$$Y=X_1^2+X_2^2+\mathrm{}X_n^2,$$
(15)
where $`n`$ is the order and the $`X_i`$ are independent Gaussian random fields with identical covariance functions. Similar models were also explored by Moscardini et al. (1991) using $`Y`$ to model either the primordial density field or the primordial gravitation potential; see also (Koyoma, Soda & Taruya 1999; Verde et al. 2000; Matarrese, Verde & Jimenez 2000; Verde et al. 2001; Komatsu & Spergel 2001). The case $`n=1`$ is the quadratic model of the present paper. There was no physical motivation for this as a model of CMB temperature fluctuations; it was used simply because one could calculate analytic results for such a field to compare with similar results for a Gaussian. Indeed the $`\chi _n^2`$ random field has been used as a model for non-Gaussian phenomena in a wide range of fields, including surface physics and geology (e.g. Adler 1981).

### 3.2 Inflation

Since the early 1980s (e.g. Guth & Pi 1982) it has been commonly believed that the inflationary scenario of the very early Universe results in the imprint of primordial Gaussian fluctuations. Since then, and with the invention of increasingly complicated models of the inflationary process, it has become clearer that inflation can produce significant levels of non-Gaussianity. The simplest “slow-roll” models of inflation involving a single dominant scalar field do indeed produce Gaussian fluctuations, but including non-linear terms and back-reaction does produce some element of non-Gaussianity as do extra degrees of freedom (Salopek & Bond 1990; Salopek 1992; Falk, Rangarajan & Srednicki 1993; Gangui et al. 1994; Wang & Kamionkowski 2000; Bartolo, Matarrese & Riotto 2002).

Typically the non-linear contributions arising during inflation manifest themselves as higher-order contributions to the effective Newtonian gravitational potential, $`\mathrm{\Phi }`$ i.e.

$$\mathrm{\Phi }=\varphi +\alpha (\varphi ^2\varphi ^2)+\mathrm{},$$
(16)
where $`\varphi `$ is a Gaussian field (not the phase) and $`\alpha `$ is a constant which is vanishingly small in most models. The term in $`\varphi ^2`$ is needed to ensure $`\mathrm{\Phi }`$ has zero mean. Since the mean value of the Newtonian potential is not physically meaningful anyway this term is not really important. A more radical suggestion for an inflation-induced quadratic model is offered by Peebles (1999a,b). In this model the density field is given by

$$\rho (𝐱)=\frac{1}{2}m^2\psi (𝐱)^2,$$
(17)
where $`\psi `$ is a scalar field and $`m`$ is an effective mass. We return to the statistical consequences of this particular model in Section 3.

### 3.3 Gravitational Non-linearity

In a simple perturbative model, the non-linear density contrast at a point $`𝐫`$ can be modelled by the relation

$$\delta (𝐱)=\delta _1(𝐱)+ϵ\delta _2(𝐱)$$
(18)
where $`\delta _1(𝐱)`$ is a Gaussian random field, $`\delta _2(𝐱)`$ is a quadratic random field derived by squaring $`\delta _1`$ and $`ϵ`$ is a small factor that controls the degree of non-linearity. To be precise we should also include a constant term in this expression in order to ensure that $`\delta (𝐱)=0`$, but this does not play any role in the following for reasons mentioned above so we ignore it. Using a constant $`ϵ`$ is not rigorous but at least qualitatively it shows how the lowest non-linear corrections come into play. For a detailed discussion of perturbation theory done properly see Bernardeau et al. (2002).

### 3.4 Bias

Attempts to confront theories of cosmological structure formation with observations of galaxy clustering are complicated by the uncertain and possibly biased relationship between galaxies and the distribution of gravitating matter. A particular simple and useful way of modelling this relationship is through the idea of a local bias. In such models, the propensity of a galaxy to form at a point where the total (local) density of matter is $`\rho `$ is taken to be some function $`f(\rho )`$ (Coles 1993; Fry & Gaztanaga 1993). This boils down to a statement of the form

$$\delta _g(𝐱)=f[\delta (𝐱)],$$
(19)
where $`\delta _g`$ is the density contrast inferred from galaxy counts or other clustering statistics. In the simplest local bias models, $`f(\delta )`$ is a constant usually called $`b`$. Clearly a linear bias of this form simply scales the variance, covariance functions and power-spectra of the underlying field but has no effect on the detailed form of the statistical distribution. Models where the bias is non-linear (but still local) are useful as they subject constraints on the effect that the bias may have on galaxy clustering statistics, without making any particular assumption about the form of $`f`$ (Coles 1993). Fry & Gaztanaga (1993) discussed the implications of bias with the form

$$f(\delta )=\underset{n=0}{\overset{\mathrm{}}{}}b_n\delta ^n,$$
(20)
in which the $`b_n`$ cannot all be chosen independently because the mean of $`\delta _g`$ must again be zero. On scale where the density field is linear, one can therefore see that a non-linear bias with $`b_20`$ will result in quadratic contributions to $`\delta _g`$ even if they do not contribute significantly to $`\delta `$.

## 4 Asymptotic Properties

In developing the statistical background in Section 2, especially the Fourier description of random fluctuating fields, we made a number of assumptions along the way that were necessary in order for the resulting descriptors to be well-defined. The terminology relating to these assumptions is often used very loosely in cosmology, at least partly because they bear a relatively simple relationship to each other under when the fields one is dealing with are Gaussian. However, in the general case of non-Gaussian fields many subtleties arise relating to the presence of higher-order statistical correlations. It is especially important in the current era of high-precision developments in statistical cosmology also to be precise about the foundations.

In the following we discuss some of the large-scale properties of random fields in a formal fashion, with particular reference to the quadratic model. The behaviour of covariance functions on large-scales will turn out to be very important so, to take the simplest example, consider the Peebles model from Section 3.2. In this case we basically have a quadratic density field $`\delta =\psi ^2`$ where $`\psi `$ is assumed a statistically homogeneous random field. Let us suppose that $`\psi `$ has a well-behaved covariance function $`\mathrm{\Gamma }(𝐫)`$ and that the covariance function of $`\delta `$ is, as usual, $`\xi (𝐫)`$. It is trivial to show that that

$$\xi (𝐫)=\mathrm{\Gamma }^2(𝐫)$$
(21)
so that $`\xi (𝐫)`$ must be positive for all $`𝐫`$. In the more general case of a field of the form $`\delta =\psi +\alpha \psi ^2`$, such as the examples given in equations (15) & (18), the resulting covariance function has a behaviour of the form

$$\xi (𝐫)=\mathrm{\Gamma }(𝐫)+\alpha ^2\mathrm{\Gamma }^2(𝐫).$$
(22)
In this case the covariance function of the resulting field would contain terms of order $`\mathrm{\Gamma }(𝐫)`$, the corresponding covariance function of the underlying Gaussian field. If $`\mathrm{\Gamma }(𝐫)<0`$ on some scale then as long as $`\alpha `$ is small, the resulting $`\xi (𝐫)`$ need not be positive in this case.

### 4.1 Statistical Homogeneity

The formal definition of strict statistical homogeneity for a random field (also called stationarity) is that the set of finite-dimensional joint probability distributions, which we called $`𝒫_N(\delta _1,\mathrm{},\delta _N)`$ in Section 2.3, must be invariant under spatial translations, i.e.

$$𝒫_N(\delta (𝐱_1),\mathrm{},\delta (𝐱_N))=𝒫_N(\delta (𝐱_1+𝐱),\mathrm{},\delta (𝐱_N+𝐱))$$
(23)
for any $`𝐱`$. This must be true for all orders $`N`$. For a Gaussian random field in which the form of $`𝒫_N(\delta _1,\mathrm{},\delta _N)`$ is given by equation (12), necessary and sufficient conditions for $`\delta (𝐱)`$ to be strictly homogeneous is that the covariance function $`\delta (𝐱_1)\delta (𝐱_2)`$ is a function of $`𝐱_1𝐱_2`$ only (Adler 1981). Statistical isotropy can be added by requiring rotation-invariance. One can define weaker versions of homogeneity and isotropy according to which only the moments of the distribution be translation invariant. For example, second-order homogeneity and isotropy (all that is required for the analysis of power-spectra or two-point covariance functions) basically means that the function $`\xi (𝐫)`$ does not depend on either the origin or the direction of $`𝐫`$, but only on its modulus. Since the properties of a Gaussian random field depend only on second-order properties, this weaker condition is sufficient to require condition (23) in this case.

This does not mean that any function satisfying translation and rotation invariance is necessarily the covariance function of a homogeneous random field (in either the strict or second-order sense). For one thing the power spectrum must be positive (or zero) for all $`k`$, which places a constraint on the shape of any possible $`\xi (𝐫)`$ – which must be convex. The result (9) also implies that

$$\underset{r\mathrm{}}{lim}\frac{1}{r^3}_0^r\xi (r^{})r^2𝑑r^{}=0$$
(24)
for such fields. A perfectly homogeneous distribution would have $`P(k)0`$ and $`\xi (r)`$ would be identically zero for all $`r`$. Note, however, that it is possible for fields obeying either (21) or (22) to be statistically homogeneous.

### 4.2 Sample Homogeneity

Statistical homogeneity plays a vital role in both the analysis of clustering and the formal development of the theory of cosmological perturbation growth. Unfortunately the use of the word “homogeneity” in this context leads to a confusion regarding the more fundamental use of this word in cosmology. Standard cosmologies are based on the Cosmological Principle, which requires our Universe to be homogeneous and isotropic on large scales. More loosely, it needs to be sufficiently homogeneous and isotropic that the Robertson-Walker metric and the Friedmann equations furnish an adequate approximation to the evolution of the Universe. Statistical homogeneity as described above is a much weaker requirement than the requirement that one realisation from a probability ensemble (our Universe) has asymptotically small fluctuations when smoothed on a sufficiently large scale. In fact, the analogous relation to (24) in the case of sample homogeneity is far stronger:

$$\underset{r\mathrm{}}{lim}_0^r\xi (r^{})r^2𝑑r^{}=0$$
(25)
as this requires the real fluctuations in density to be asymptotically small within a single realisation. Notice that the requirement for sample homogeneity is such that in general the covariance function must changes sign, from positive at the origin, where $`\xi (0)=\sigma ^20`$, to negative at some $`r`$ to make the overall integral (24) converge in the correct way.

It is clear from this discussion that statistical homogeneity does not require sample homogeneity. Revisiting the quadratic model now reveals another interesting point: the Peebles model (21) can not be sample homogeneous, even if it is statistically homogeneous. If we want $`\delta `$ to have a covariance function matching observations, say $`\xi (r)r^2`$, then the underlying Gaussian field must have $`\mathrm{\Gamma }(r)r^1`$ which violates the constraint for it to be sample homogeneous.

This model behaves in a similar way to fractal models of the type discussed, for example, by Coleman & Pietronero (1992). Mention of fractal models tends to send mainstream cosmologists screaming into the hills, but the lack of sample homogeneity they describe need not be damaging to standard methods. To give a prominent illustration, consider the behaviour of the gravitational potential field defined by a Gaussian random field with the Harrison-Zel’dovich spectrum. In such a case, the fluctuations will always be small (of order $`10^5`$ to be consistent with observations) but they are independent of scale, and thus there is never a scale at which sample homogeneity is reached. It is not particularly important for the purposes of galaxy clustering studies that the universe obeys the property of sample homogeneity. What is more important is estimates of statistical properties obtained from different samples vary with respect to the ensemble-averaged property in a fashion which is under control for large samples. This does not require asymptotic convergence to homogeneity.

Finally, note that the general quadratic model (22) can be sample homogeneous if $`\mathrm{\Gamma }(r)`$ obeys condition (25) and $`\alpha `$ is sufficiently small. Perturbative corrections, such as those described in Section 3.3, do not therefore necessarily induce sample inhomogeneity.

### 4.3 Ergodicity and Fair Samples

We have already introduced the idea of a “fair sample hypothesis”, which is basically that averages over finite patches of the Universe can be treated as averages over some probability ensemble. Peebles (1980), for example, gives the definition of a fair sample hypothesis in a number of ways. First, he states

> “..the fair sample hypothesis is taken to mean that the universe is statistically homogeneous and isotropic.”

Later we find

> “Samples from well separated spots are uncorrelated, and the collection of such samples is a statistical ensemble generated by many independent applications…”

This second definition is close to the one we used in Section 2, but it is clear that it is stronger than the first one. Related to the fair sample hypothesis, but not identical to it, is the so-called ergodic property, which is that averages over an infinite domain within a single realization can be treated as averages over the probability ensemble. The second definition of a fair sample is a stronger statement than the ergodic property, since it involves the properties of finite patches rather than an infinite domain within a single realisation from the probability ensemble.

Ergodic properties are extremely difficult to prove, but results do exist for Gaussian random fields (Adler 1981). Intriguingly, in this case the result is extremely simple. The necessary and sufficient condition for a statistically homogeneous Gaussian random field to be ergodic is that its power-spectrum (defined above) should be continuous. Continuity of the power-spectrum leads, by standard Fourier analysis, to the result that

$$\underset{r\mathrm{}}{lim}\frac{1}{r^3}_0^r[\xi (r^{})]^2r^2𝑑r^{}=0.$$
(26)
This requires the covariance function to be decreasing. In fact, any statistically homogeneous Gaussian field will be ergodic if $`\xi (𝐫)0`$ as $`𝐫\mathrm{}`$. Notice then that a Gaussian random field can be ergodic without being sample homogeneous.

A general form of this ergodic theorem does not exist for arbitrary non-Gaussian random fields, but fortunately this does not matter. What is needed for statistical cosmology is not an ergodic property but something closer to a version of the fair sample hypothesis. The ergodic property does not require the fair-sample property, as it involves averages over infinite domains of a single realisation. The fair sample hypothesis does not require ergodicity, either. From this we can conclude that the ergodic property is irrelevant and use of this term should be avoided.

There is, however, one particularly neat relationship between ergodicity and statistical homogeneity for the Peebles model. Notice if we take $`\delta \psi ^2`$ and require $`\psi `$ to be a Gaussian random field with covariance function $`\xi (𝐫)`$ then the covariance function of $`\delta `$ is just $`\xi ^2(𝐫)`$, exactly the form that appears in equation (26). If we take $`\psi `$ to be an ergodic Gaussian random field then this guarantees the resulting quadratic field must be at least second-order statistically homogeneous.

### 4.4 Asymptotic Independence and Smoothing

The considerations we have discussed above generally lead to requirements that the correlations between points become small as the separation between the points grows large. For a Gaussian field, if $`\xi (𝐫)0`$ as $`r\mathrm{}`$ then the probability distributions tend to an asymptotically independent form. This must be the case because such a field contains only second-order correlations. For example, in the limit that the correlation matrix $`C_{ij}`$ becomes diagonal, the 2-point Gaussian probability density $`𝒫_2(\delta _1,\delta _2)𝒫_1(\delta _1)𝒫_1(\delta _2)`$, consistent with the requirement of independence i.e. $`P(A,B)=P(A)P(B)`$. Similar results will hold for higher order $`N`$–point distributions. This result means that for Gaussian fields absence of (second-order) correlation, i.e.

$$X_1X_2=X_1X_2,$$
(27)
means independence. Lack of correlation only requires full independence for Gaussian fields. Independence always implies lack of correlation whether the field is Gaussian or not.

We can see then that even though a non-Gaussian field, such as a quadratic model, may be uncorrelated on large scales, consistent with the requirements above, this does not necessarily mean that points are asymptotically independent.

The reason for discussing this is that it is very relevant to what happens to a random field as it is smoothed on successively larger scales. This smoothing is equivalent to filtering the field with a low pass filter. The filtered field, $`\delta (𝐱;R_f)`$, may be obtained by convolution of the “raw” density field with some function $`F`$ having a characteristic scale $`R_f`$:

$$\delta (𝐱;R_f)=\delta (𝐱^{})F(|𝐱𝐱^{}|;R_f)𝑑𝐱^{}.$$
(28)
The filter $`F`$ has the following properties: $`F=\mathrm{constant}R_f^3`$ if $`|𝐱𝐱^{}|R_f`$, $`F0`$ if $`|𝐱𝐱^{}|R_f`$, $`F(𝐲;R_f)𝑑𝐲=1`$.

If the underlying density field is Gaussian then the filtered field will also be Gaussian. This is a result of the fact that filtering essentially constructs a weighted average of the underlying field and any sum of Gaussian variates is itself a Gaussian variate (e.g. Kendall & Stuart 1977). According to the Central Limit Theorem, the sum of a large number of independent variates drawn from a distribution with finite variance also tends to a Gaussian distribution. One would imagine, therefore, that if distant points were asymptotically independent then the effect of filtering on a non-Gaussian field is to “gaussianize” it. In fact this is assumed in standard statistical cosmology. We know the small-scale distribution is non-Gaussian, but averaging over sufficiently large smoothing windows is assumed to recover the linear field, or something close to it. But for a general non–Gaussian field how quickly do points have to become independent in order for filters to Gaussianize the distribution?

The answer to this question is given by Fan & Bardeen (1995): it depends on a function called the Rosenblatt dependence rate which governs the rate at which the maximum value of $`|P(AB)P(A)P(B)|`$ tends to zero at large separations ($`A`$ and $`B`$ are any combination of values of $`\delta _i`$ in different locations). This is quite a technical issue, and we leave the details to Fan & Bardeen (1995). On the other hand when the field in question is a local transformation of a Gaussian random field (such as in the quadratic model) then there is a simple result that the covariance function of the underlying field must fall at least as quickly as $`1/r^3`$ as $`r\mathrm{}`$. This is sometimes referred to as the pseudo-Markov condition (Adler 1981). If this criterion is satisfied then all local transformations of the underlying field satisfy the mixing rate condition and consequently become Gaussian if smoothed on sufficiently large scales.

Let us look at this issue in the light of the Peebles model. Peebles (1999b) shows that this model is fully self-similar. When smoothed on larger and larger scales the distribution function does not tend to a Gaussian but retains the same ($`\chi ^2`$) form.

At first sight this looks extremely surprising. Suppose we imagine a simple version of a $`\chi ^2`$ random field built upon discrete cells of some size $`R`$. Suppose the field were uniform within each cell but that adjacent cells were generated independently from the $`\chi ^2`$ distribution. Filtering on a scale $`R_f`$ in this case simply corresponds to adding neighbouring (uncorrelated) cells, producing something like a sum of $`N=(R_f/R)^3`$ independent values from a quadratic model. This would give rise to a resulting field of the form (15). As $`R_f`$ becomes larger, $`N`$ increases and the resulting distribution changes to a distribution of $`\chi _n^2`$ of larger and larger $`n`$. I=t is well known that, as $`n\mathrm{}`$, the distribution of $`Y_n`$ is asymptotically close to the Gaussian as expected from the central limit theorem.

In the case of the Peebles model, however, the mixing rate condition is not satisfied. Although distant points are asymptotically independent, the rate at which they tend to independence is not sufficient to produce a Gaussian distribution. Regardless of the scale of smoothing, as long as the covariance function is chosen to be scale-free, the density field retains a distribution having the same shape. This demonstrates that this model is, in fact, a kind of fractal as we discussed above.

Notice that in the more normal case where we take the quadratic contribution to represent only the non-linear effect on initially Gaussian fluctuations then it is guaranteed to satisfy the mixing rate condition as long as the Gaussian field does that generates it. It is therefore justified to assume that the model described in Section 3.3 does become Gaussian when smoothed on large scales.

## 5 Quadratic Phase Coupling

In section 2 we pointed out that a convenient definition of a Gaussian field could be made in terms of its Fourier phases, which should by independent and uniformly distributed on the interval $`[0,2\pi ]`$. A breakdown of these conditions, such as the correlation of phases of different wavemodes, is a signature that the field has become non-Gaussian. In terms of cosmic large-scale structure formation, non-Gaussian evolution of the density field is symptomatic of the onset of non-linearity in the gravitational collapse process, suggesting that phase evolution and non-linear evolution are closely linked. A relatively simple picture emerges for models where the primordial density fluctuations are Gaussian and the initial phase distribution is uniform. When perturbations remain small evolution proceeds linearly, individual modes grow independently and the original random phase distribution is preserved. However, as perturbations grow large their evolution becomes non-linear and Fourier modes of different wavenumber begin to couple together. This gives rise to phase association and consequently to non-Gaussianity. It is clear that phase associations of this type should be related in some way to the existence of the higher order connected covariance functions, which are traditionally associated with non-linearity and are non-zero only for non-Gaussian fields. In this sections such a relationship is explored in detail using an analytical model for the non-linearly evolving density fluctuation field. Phase correlations of a particular form are identified and their connection to the covariance functions is established.

### 5.1 A simple non-linear model

We adopt the simple perturbative expansion of equation (18) in order to model the non-linear evolution of the density field. Although the equivalent transformation in formal Eulerian perturbation theory is a good deal more complicated, the kind of phase associations that we will deal with here are precisely the same in either case. In terms of the Fourier modes, in the continuum limit, we have for the first order Gaussian term

$$\delta _1(𝐱)=d^3k|\delta _𝐤|\mathrm{exp}[i\varphi _𝐤]\mathrm{exp}[i𝐤𝐱]$$
(29)
and for the second-order perturbation

$$\delta _2(𝐱)=\left[\delta _1(𝐱)\right]^2=d^3kd^3k^{}|\delta _𝐤||\delta _𝐤^{}|\mathrm{exp}[i(\varphi _𝐤+\varphi _𝐤^{})]\mathrm{exp}[i(𝐤+𝐤^{})𝐫].$$
(30)
The quadratic field, $`\delta _2`$, illustrates the idea of mode coupling associated with non-linear evolution. The non-linear field depends on a specific harmonic relationship between the wavenumber and phase of the modes at $`𝐤`$ and $`𝐤^{}`$. This relationship between the phases in the non-linear field, i.e.

$$\varphi _𝐤+\varphi _𝐤^{}=\varphi _{𝐤+𝐤^{}},$$
(31)
where the RHS represents the phase of the non-linear field, is termed quadratic phase coupling.

### 5.2 The two-point covariance function

The two-point covariance function can be calculated using the definitions of section 2, namely

$$\xi (r)=\delta (𝐱)\delta (𝐱+𝐫).$$
(32)
Substituting the non-linear transform for $`\delta (𝐱)`$ (equation 18) into this expression gives four terms

$$\xi (r)=\delta _1(𝐱)\delta _1(𝐱+𝐫)+ϵ\delta _1(𝐱)\delta _2(𝐱+𝐫)+ϵ\delta _2(𝐱)\delta _1(𝐱+𝐫)+ϵ^2\delta _2(𝐱)\delta _2(𝐱+𝐫).$$
(33)
The first of these terms is the linear contribution to the covariance function whereas the remaining three give the non-linear corrections. We shall focus on the lowest order term for now.

As we outlined in Section 2, the angle brackets $``$ in these expressions are expectation values, formally denoting an average over the probability distribution of $`\delta (𝐱)`$. Under the fair sample hypothesis we replace the expectation values in equation (32) with averages over a selection of independent volumes so that $`_{\text{vol, real}}`$. The first average is simply a volume integral over a sufficiently large patch of the universe. The second average is over various realisations of the $`\delta _k`$ and $`\varphi _k`$ in the different patches. Applying these rules to the first term of equation (33) and performing the volume integration gives

$$\xi _{11}(r)=d^3kd^3k^{}|\delta _𝐤||\delta _𝐤^{}|\mathrm{exp}[i(\varphi _𝐤+\varphi _𝐤^{})]_{\text{real}}\delta _D(𝐤+𝐤^{})\mathrm{exp}[i𝐤^{}𝐬],$$
(34)
where $`\delta _D`$ is the Dirac delta function. The above expression is simplified given the reality condition

$$\delta _𝐤=\delta _𝐤^{},$$
(35)
from which it is evident that the phases obey

$$\varphi _𝐤+\varphi _𝐤=0\text{mod}[2\pi ].$$
(36)
Integrating equation (34) one therefore finds that

$$\xi _{11}(r)=d^3k|\delta _𝐤|^2_{\text{real}}\mathrm{exp}[i𝐤𝐬].$$
(37)
so that the final result is independent of the phases. Indeed this is just the Fourier transform relation between the two-point covariance function and the power spectrum we derived in section 2.1.

### 5.3 The three-point covariance function

Using the same arguments outlined above it is possible to calculate the 3-point connected covariance function, which is defined as

$$\zeta (𝐫,𝐬)=\delta (𝐱)\delta (𝐱+𝐫)\delta (𝐱+𝐬)_c.$$
(38)
Making the non-linear transform of equation (18) one finds the following contributions

$`\zeta (𝐫,𝐬)`$ $`=`$ $`\delta _1(𝐱)\delta _1(𝐱+𝐫)\delta _1(𝐱+𝐬)_c+ϵ\delta _1(𝐱)\delta _1(𝐱+𝐫)\delta _2(𝐱+𝐬)_c+\text{perms}(121,211)`$ (39)
$`+ϵ^2\delta _1(𝐱)\delta _2(𝐱+𝐫)\delta _2(𝐱+𝐬)_c+\text{perms}(212,221)`$
$`+ϵ^3\delta _2(𝐱)\delta _2(𝐱+𝐫)\delta (𝐱+𝐬)_c.`$
Again we consider first the lowest order term. Expanding in terms of the Fourier modes and once again replacing averages as prescribed by the fair sample hypothesis gives

$$\zeta _{111}(𝐫,𝐬)=d^3kd^3k^{}d^3k^{\prime \prime }|\delta _𝐤||\delta _𝐤^{}||\delta _{𝐤^{\prime \prime }}|\mathrm{exp}[i(\varphi _𝐤+\varphi _𝐤^{}+\varphi _{𝐤^{\prime \prime }})]_{\text{real}}\delta _D(𝐤+𝐤^{}+𝐤^{\prime \prime })\mathrm{exp}[i𝐤^{}𝐫]\mathrm{exp}[i𝐤^{\prime \prime }𝐬].$$
(40)
Recall that $`\delta _1`$ is a Gaussian field so that $`\varphi _𝐤`$, $`\varphi _𝐤^{}`$ and $`\varphi _{𝐤^{\prime \prime }}`$ are independent and uniformly random on the interval $`[0,2\pi ]`$. Upon integration over one of the wavevectors the phase terms is modified so that its argument contains the sum $`(\varphi _𝐤+\varphi _𝐤^{}+\varphi _{𝐤𝐤^{\prime \prime }})`$, or a permutation thereof. Whereas the reality condition of equation (35) implies a relationship between phases of anti-parallel wavevectors, no such conditions hold for modes linked by the triangular constraint imposed by the Dirac delta function. In other words, except for serendipity,

$$\varphi _𝐤+\varphi _𝐤^{}+\varphi _{𝐤𝐤^{\prime \prime }}0.$$
(41)
In fact due to the circularity of phases, the resulting sum is still just uniformly random on the interval $`[0,2\pi ]`$ if the phases are random. Upon averaging over sufficient realisations, the phase term will therefore cancel to zero so that the lowest order contribution to the 3-point function vanishes, i.e. $`\zeta _{111}(𝐫,𝐬)=0`$. This is not a new result, but it does explicitly illustrate how the vanishing of the three-point connected covariance function arises in terms of the Fourier phases.

Next we consider the first non-linear contribution to the 3-point function given by

$$\zeta _{112}(𝐫,𝐬)=ϵ\delta _1(𝐱)\delta _1(𝐱+𝐫)\delta _2(𝐱+𝐬),$$
(42)
or one of its permutations. In this case one of the arguments in the average is the field $`\delta _2(𝐱)`$, which exhibits quadratic phase coupling of the form (31). Expanding this term to the point of equation (40) using the definition (30) one obtains

$`\zeta _{112}(𝐫,𝐬)`$ $`=`$ $`{\displaystyle d^3kd^3k^{}d^3k^{\prime \prime }d^3k^{\prime \prime \prime }|\delta _𝐤||\delta _𝐤^{}||\delta _{𝐤^{\prime \prime }}||\delta _{𝐤^{\prime \prime \prime }}|\mathrm{exp}[i(\varphi _𝐤+\varphi _𝐤^{}+\varphi _{𝐤^{\prime \prime }}+\varphi _{𝐤^{\prime \prime \prime }})]_{\text{real}}}`$ (43)
$`\times \delta _D(𝐤+𝐤^{}+𝐤^{\prime \prime }+𝐤^{\prime \prime \prime })\mathrm{exp}[i𝐤^{}𝐫]\mathrm{exp}[i(𝐤^{\prime \prime }+𝐤^{\prime \prime \prime })𝐬].`$
Once again the Dirac delta function imposes a general constraint upon the configuration of wavevectors. Integrating over one of the $`𝐤`$ gives $`𝐤^{\prime \prime \prime }=𝐤𝐤^{}𝐤^{\prime \prime }`$ for example, so that the wavevectors must form a closed loop. This general constraint however, does not specify a precise shape of loop, instead the remaining integrals run over all of the different possibilities. At this point we may constrain the problem more tightly by noting that most combinations of the $`𝐤`$ will contribute zero to $`\zeta _{112}`$. This is because of the circularity property of the phases and equation (41). Indeed, the only nonzero contributions arise where we are able to apply the phase relation obtained from the reality constraint, equation (36). In other words the properties of the phases dictate that the wavevectors must align in anti-parallel pairs: $`𝐤=𝐤^{}`$, $`𝐤^{\prime \prime }=𝐤^{\prime \prime \prime }`$ and so forth.

There is a final constraint that must be imposed upon the $`𝐤`$ if $`\zeta `$ is the connected $`3`$-point covariance function. In a graph theoretic sense, the general (unconnected) $`N`$-point function $`\delta _{l_1}(𝐱_1)\delta _{l_2}(𝐱_2)\mathrm{}\delta _{l_N}(𝐱_N)`$ can be represented geometrically by a sum of tree diagrams. Each diagram consists of $`N`$ nodes of order $`l_i`$, representing the $`\delta _{l_i}(𝐱_i)`$, and a number of linking lines denoting their correlations; see Fry (1984) or Bernardeau (1992) for more detailed accounts. Every node is made up of $`l_i`$ internal points, which represent a factor $`\delta _𝐤=|\delta _𝐤|\mathrm{exp}(i\varphi _𝐤)`$ in the Fourier expansion. According to the rules for constructing diagrams, linking lines may join one internal point to a single other, either within the same node or in an external node (see Figure 1). The connected covariance functions are represented specifically by the subset of diagrams for which every node is linked to at least one other, leaving none completely isolated (Figure 1; right). This constraint implies that certain pairings of wavevectors do not contribute to the connected covariance function (Figure 2).

The above constraints may be inserted into equation (43) by re-writing the Dirac delta function as a product over Delta functions of two arguments, appropriately normalised. There are only two allowed combinations of wavevectors so we have

$$\delta _D(𝐤+𝐤^{}+𝐤^{\prime \prime }+𝐤^{\prime \prime \prime })\frac{1}{2V_u}[\delta _D(𝐤+𝐤^{\prime \prime })\delta _D(𝐤^{\prime \prime \prime }+𝐤^{\prime \prime \prime })+\delta _D(𝐤+𝐤^{\prime \prime \prime })\delta _D(𝐤^{}+𝐤^{\prime \prime })].$$
(44)
Integrating over two of the $`𝐤`$ and using equation (36) eliminates the phase terms and leaves the final result

$$\zeta _{112}(𝐫,𝐬)=\frac{1}{V_u}d^3kd^3k^{}|\delta _𝐤|^2|\delta _𝐤^{}|^2_{\text{real}}\mathrm{exp}[i𝐤^{}𝐫]\mathrm{exp}[i(𝐤+𝐤^{})𝐬].$$
(45)
The existence of this quantity has therefore been shown to depend on the quadratic phase coupling of Fourier modes. The relationship between modes and the interpretation of the tree diagrams is also dictated by the properties of the phases.

One may apply the same rules to the higher order terms in equation (39). It is immediately clear that the $`\zeta _{122}`$ terms are zero because there is no way to eliminate the phase term $`\mathrm{exp}[i(\varphi _𝐤+\varphi _𝐤^{}+\varphi _{𝐤^{\prime \prime }}+\varphi _{𝐤^{\prime \prime \prime }}+\varphi _{𝐤^{\prime \prime \prime \prime }})]`$, a consequence of the property equation (41). Diagrammatically this corresponds to an unpaired internal point within one of the nodes of the tree. The final, highest order contribution to the 3-point function is found to be

$$\zeta _{222}(𝐫,𝐬)=\frac{1}{V_u^2}d^3kd^3k^{}d^3k^{\prime \prime }|\delta _𝐤|^2|\delta _𝐤^{}|^2|\delta _{𝐤^{\prime \prime }}|^2_{\text{real}}\mathrm{exp}[i(𝐤𝐤^{})𝐫]\mathrm{exp}[i(𝐤^{}𝐤^{\prime \prime })𝐬],$$
(46)
where the phase and geometric constraints allow 12 possible combinations of wavevectors.

### 5.4 Cubic non-linearity and higher order

The above ideas extend simply to higher order where the non-linear field is represented by a perturbation series that does not truncate at the quadratic term. At the next highest order for example, the series includes $`\delta _3=\delta _1^3`$, which introduces a cubic phase coupling in the Fourier expansion. Although quadratic phase coupling is essential as a minimum requirement for the three point covariance function, cubic phase coupling is not the minimum requirement for the next highest order covariance function. Indeed, quadratic coupling is sufficient to provide contributions to all of the n-point covariance functions due to the way the phases dictate that wave-vectors must arrange themselves into antiparallel pairs. For the 4-point covariance function the cubic term allows for a different diagrammatic representation: a star as opposed to a snake topology. However, in terms of the constituent wave-vectors, loops (as in Figure 2) contributing to the star topologies are a symmetric subset of those contributing to the snake topologies.

### 5.5 Power-spectrum and Bispectrum

The formal development of the relationship between covariance functions and power-spectra developed above suggests the usefulness of higher–order versions of $`P(k)`$. It is clear from the arguments of Section 5.2 that a more convenient notation for the power-spectrum than that introduced in section 2.1 is

$$\delta _𝐤\delta _𝐤^{}=(2\pi )^3P(k)\delta _D(𝐤+𝐤^{}).$$
(47)
The connection between phases and higher-order covariance functions obtained in Section 5.3 also suggests defining higher-order polyspectra of the form

$$\delta _𝐤\delta _𝐤^{}\mathrm{}\delta _{𝐤^{(n)}}=(2\pi )^3P_n(𝐤,𝐤^{},\mathrm{}𝐤^{(n)})\delta _D(𝐤+𝐤^{}+\mathrm{}𝐤^{(n)})$$
(48)
where the occurrence of the delta-function in this expression arises from a generalisation of the reality constraint given in equation (36); see, e.g., Peebles (1980). Conventionally the version of this with $`n=3`$ produces the bispectrum, usually called $`B(𝐤,𝐤^{},𝐤^{\prime \prime })`$ which has found much effective use in recent studies of large-scale structure (Peebles 1980; Scoccimarro et al. 1998; Scoccimarro, Couchman & Frieman 1999; Verde et al. 2000; Verde et al. 2001; Verde et al. 2002). It is straightforward to show that the bispectrum is the Fourier-transform of the (reduced) three–point covariance function by following similar arguments as in Section 5.2; see, e.g., Peebles (1980).

Note that the delta-function constraint requires the bispectrum to be zero except for $`k`$-vectors ($`𝐤`$, $`𝐤^{}`$, $`𝐤^{\prime \prime }`$) that form a triangle in $`k`$-space. From Section 5.3 it is clear that the bispectrum can only be non-zero when there is a definite relationship between the phases accompanying the modes whose wave-vectors form a triangle. Moreover the pattern of phase association necessary to produce a real and non-zero bispectrum is precisely that which is generated by quadratic phase association. This shows, in terms of phases, why it is that the leading order contributions to the bispectrum emerge from second-order fluctuations of a Gaussian random field. The bispectrum measures quadratic phase coupling.

Three-point phase correlations have another interesting property. While the bispectrum is usually taken to be an ensemble-averaged quantity, as defined in equation (48), it is interesting to consider products of terms $`\delta _𝐤\delta _𝐤^{}\delta _{𝐤^{\prime \prime }}`$ obtained from an individual realisation. According to the fair sample hypothesis discussed above we would hope appropriate averages of such quantities would yield an estimate of the bispectrum. Note that

$$\delta _𝐤\delta _𝐤^{}\delta _{𝐤^{\prime \prime }}=\delta _𝐤\delta _𝐤^{}\delta _{𝐤𝐤^{}}=\delta _𝐤\delta _𝐤^{}\delta _{𝐤+𝐤^{}}^{}\beta (𝐤,𝐤^{}),$$
(49)
using the requirement (36) and the triangular constraint we discussed above. Each $`\beta (𝐤,𝐤^{})`$ will carry its own phase, say $`\varphi _{𝐤,𝐤^{}}`$, which obeys

$$\varphi _{𝐤,𝐤^{}}=\varphi _𝐤+\varphi _𝐤^{}\varphi _{𝐤+𝐤^{}}.$$
(50)
It is evident from this that it is possible to recover the complete set of phases $`\varphi _𝐤`$ from the bispectral phases $`\varphi _{𝐤,𝐤^{}}`$, up to a constant phase offset corresponding to a global translation of the entire structure (Chiang & Coles 2000). This furnishes a conceptually simple method of recovering missing or contaminated phase information in a consistent way, an idea which has been exploited, for example, in speckle interferometry (Lohmann, Weigelt & Wirnitzer 1983). In the case of quadratic phase coupling, described by equation (31), the left-hand-side of equation (50) is identically zero leading to a particularly simple approach to this problem, which we shall explore in future work.

## 6 Discussion

In this paper we have addressed two main issues, using the quadratic model as an illustrative example. First we showed explicitly how this non-Gaussian model has properties that contradict standard folklore based on the assumption of Gaussian fluctuations. We used this model to distinguish carefully between various inter-related concepts such as sample homogeneity, statistical homogeneity, asymptotic independence, ergodicity, and so on. We showed the conditions under which each of these is relevant and deployed the quadratic model for particular examples in which they are violated.

We then used the quadratic model to show, for the first time, how phase association arises in non-linear processes which has exactly the correct form to generate non-zero bispectra and three–point covariance functions. The magnitude of these statistical descriptors is of course related to the magnitude of the Fourier modes, but the factor that determines whether they are zero or non-zero is the arrangement of the phases of these modes.

The connection between polyspectra and phase information is an important one and it opens up many lines of future research, such as how phase correlations relate to redshift distortion and bias. Also, we assumed throughout this study that we could straightforwardly take averages over a large spatial domain to be equal to ensemble averages. Using small volumes of course leads to sampling uncertainties which are quite straightforward to deal with in the case of the power-spectra but more problematic for higher-order spectra like the bispectrum. Understanding the fluctuations about ensemble averages in terms of phases could also lead to important insights.

## Acknowledgements

PC acknowledges useful discussions with John Peacock, Simon White and Sabino Matarrese about matters related to Section 4 of this paper. This work was supported by PPARC grant number PPA/G/S/1999/00660.

## References

Adler R.J., 1981, The Geometry of Random Fields. John Wiley & Sons, New York.

Bardeen J.M., Bond J.R., Kaiser N., Szalay A.S., 1986, ApJ, 304,

Bartolo N., Matarrese S., Riotto A., 2002, Phys. Rev. D., 65, 103505

Bernardeau F., 1992, ApJ, 392, 1

Bernardeau F., Colombi S., Gaztanaga E., Scoccimarro R., 2002, Phys. Rep., in press, astro-ph/0112551

Chiang L.-Y., 2001, MNRAS, 325, 405

Chiang L.-Y., Coles P., 2000, MNRAS, 311, 809

Chiang L.-Y., Coles P., Naselsky P., 2002, MNRAS, in press, astro-ph/0207584

Coleman P.H., Pietronero L., 1992, Phys. Rep. 231, 311

Coles P., 1993, MNRAS, 262, 1065

Coles P., Barrow J.D., 1987, MNRAS, 228, 407

Coles P., Chiang L.-Y., 2000, Nature, 406, 376

Coles P., Jones B.J.T., 1991, MNRAS, 248, 1

Falk T., Rangarajan R., Srednicki M., 1993, ApJ, 403, L1

Fan Z., Bardeen J.M., 1995, Phys. Rev. D., 51, 6714

Fry J.N., 1984, ApJ, 279, 499

Fry J.N., Gaztanaga E., 1993, ApJ, 413, 447

Gangui A., Lucchin F., Matarrese S., Mollerach S., 1994, ApJ, 430, 447

Guth A.H., 1981, Phys. Rev. D., 23, 347

Guth A.H., Pi S.-Y., 1982, Phys. Rev. Lett., 49, 1110

Hobson M.P., Jones A.W., Lasenby A., 1999, MNRAS, 309, 125

Jain B., Bertschinger E., 1998, ApJ, 509, 517

Kendall M., Stuart A., 1977, The Advanced Theory of Statistics, Vol. 1. Griffin & Co, London

Komatsu E., Spergel D.N., 2001, Phys. Rev. D., 63, 063002

Koyama K., Soda J., Taruya A., 1999, MNRAS, 310, 1111

Lohmann A.W., Weigelt G., Wirnitzer B., 1983, Appl. Optics, 22, 4028

Luo X., 1994, ApJ, 427, L71

Luo X., Schramm D.N., 1993, ApJ, 408, 33

Matarrese S., Verde L., Jimenez R., 2000, ApJ, 541, 10

Moscardini L., Matarrese S., Lucchin F., Messina A., 1991, MNRAS, 248, 424

Peacock J.A., Dodds S., 1996, MNRAS, 267, 1020

Peebles P.J.E., 1980, The Large-Scale Structure of the Universe. Princeton University Press, Princeton NJ.

Peebles P.J.E., 1999a, ApJ, 510, 523

Peebles P.J.E., 1999b, ApJ, 510, 531

Percival W.J. et al. 2001, MNRAS, 327, 1297

Ryden B.S., Gramann M., 1991, ApJ, 383, L33

Salopek D.S., 1992, Phys. Rev. D., 45, 1139

Salopek D.S., Bond J.R., Phys. Rev. D., 42, 3936

Scherrer R.J., Melott A.L., Shandarin S.F., 1991, ApJ, 377, 29

Schmalzing J., Gorski K.M., 1998, MNRAS, 297, 355

Scoccimarro R., Colombi S., Fry J.N., Frieman J.A., Hivon E., Melott A.L., 1998, ApJ, 496, 586

Scoccimarro R., Couchman H.M.P., Frieman J.A., 1999, ApJ, 517, 531

Shandarin S.F., 2002, MNRAS, 331, 865

Soda J., Suto Y., 1992, ApJ, 396, 379

Stirling A.J., Peacock J.A., 1996, MNRAS, 283, L99

Verde L., Heavens A.F., 2001, ApJ, 553, 14

Verde L., Jimenez R., Kamionkowski M., Matarrese S., 2001, MNRAS, 325, 412

Verde L., Wang L., Heavens A.F., Kamionkowski M., 2000, MNRAS, 313, 141

Verde L. et al. 2002, MNRAS, in press, astro-ph/0112161

Wang L., Kamionkowski M., 2000, Phys. Rev. D., 61, 063504


format: text

text:
# Theoretical Analysis and Simulations of the Generalized Lotka-Volterra Model

## I Introduction

In recent years there has been considerable interest in the collection and analysis of large volumes of economic data. Such data includes the distributions of the income and wealth of individuals Pareto1897 ; Mandelbrot1961 ; Mandelbrot1951 ; Mandelbrot1963 ; Atkinson1978 ; Levy1997 ; Takayasu1997 ; Jogi1998 , the market values of publicly traded companies as well as their short and long term fluctuations Mantegna1995 ; Gopikrishnan1998 ; Plerou1999 ; Gopikrishnan1999 . A common observation is that distributions of economic data exhibit a power-law behavior of the form

$$P(w)w^{1\alpha }.$$
(1)
where the variable $`w`$ represents the wealth of an individual or market value of a company and $`\alpha `$ is the exponent that provides the best fit to the empirical data. Empirical studies show that the distribution the wealth of individuals in different countries follows the power-law behavior described by Eq. (1), with $`1<\alpha <2`$ Pareto1897 ; Mandelbrot1951 ; Mandelbrot1961 ; Mandelbrot1963 ; Atkinson1978 ; Levy1997 . These results stimulated theoretical studies in attempt to construct models that reproduce the power-law behavior and predict the value of $`\alpha `$ Levy1995 ; Levy1996b ; Solomon1996 ; Sornette1997 ; Zanette1997 ; Biham1998 ; Malcai1999 .

In this paper we study a stochastic dynamical model, based on the Lotka-Volterra system that gives rise to the power-law distribution of Eq. (1). The model consists of coupled dynamic equations which describe the discrete time evolution of the basic system components $`w_i`$, $`i=1,\mathrm{},N`$. The structure of these equations resembles the logistic map and they are coupled through the average value $`\overline{w}(t)`$. The dynamics includes autocatalysis both at the individual level and at the community level as well as a saturation term. To model the non-stationary conditions we introduce a time dependent parameter into the saturation term in each of these equations. We find that the system components spontaneously evolve into a power-law distribution of the form of Eq. (1), even in the presence of non-stationary external conditions. Furthermore, it is shown analytically that the exponent $`\alpha `$ depends only on the ratio of the constant drift component (social security) and the fluctuating component (investments). It is found to be insensitive to variations in the saturation term that describes the level of economic activity and varies between periods of prosperity and depression.

The paper is organized as follows. In Sec. II we present the generalized Lotka-Volterra model under non-stationary conditions. Analytical results and predictions are presented in Sec. III and compared with the results of numerical simulations in Sec. IV. A summary is presented in Sec. V.

## II The Model

The generalized Lotka-Volterra system Levy1995 ; Levy1996b ; Solomon1996 ; Biham1998 ; Malcai1999 ; Biham2001  describes the evolution in discrete time of $`N`$ dynamic variables $`w_i`$, $`i=1,\mathrm{},N`$. In ecological systems, $`w_i`$ represents the population size of the $`i`$th specie, while in economic systems it may represent the wealth of an individual investor or the market value of a publicly traded firm. At each time step $`t`$, an integer $`i`$ is chosen randomly in the range $`1iN`$, which is the index of the dynamic variable $`w_i`$ to be updated at that time step. A random multiplicative factor $`\lambda (t)`$ is then drawn from a given distribution $`\mathrm{\Pi }(\lambda )`$, which is independent of $`i`$ and $`t`$. It will later be convenient to express this multiplicative factor by

$$\lambda (t)=\lambda +\eta (t)$$
(2)
where $`\lambda `$ is the average value of $`\mathrm{\Pi }(\lambda )`$ and

$$D=\lambda ^2\lambda ^2$$
(3)
is its standard deviation. The system is then updated according to

$`w_i(t+1)`$ $`=`$ $`[1+\lambda (t)]w_i(t)+{\displaystyle \underset{j=1}{\overset{N}{}}}a_{i,j}w_j(t){\displaystyle \underset{j=1}{\overset{N}{}}}a_{j,i}w_i(t){\displaystyle \underset{j=1}{\overset{N}{}}}c_{i,j}w_i(t)w_j(t)`$
$`w_j(t+1)`$ $`=`$ $`w_j(t),j=1,\mathrm{},N;ji.`$ (4)
where $`a_{i,j}`$ and $`c_{i,j}`$ are constants. This is an asynchronous update mechanism. The first term on the right hand side of Eq. (4), describes the effect of stochastic auto-catalysis at the individual level. In an ecological system this term represents variations in the population of a given specie, including births and deaths that may be affected by external conditions but are not affected by the interaction with other species. In a stock-market system it represents the increase (or decrease) by a random factor $`\lambda (t)`$ of the capital of the investor $`i`$ between time $`t`$ and $`t+1`$. The second and third terms in Eq. (4), describe the interaction between different dynamical variables. In an ecological system, the second term represents the dependence of population $`i`$ on the availability of food, in the form of population $`j`$. The third term represents the fact that population $`i`$ itself may be the food of some other species. In an economic system the second and third terms represent trade between investors or firms $`i`$ and $`j`$, such as buying and selling, respectively. The fourth term in Eq. (4), describes saturation effects due to the competition for limited resources. In an ecological model, this term implies that large populations tend to exhaust the available resources on which they depend. The saturation parameters $`c_{i,j}`$ are large for populations $`i`$ and $`j`$ that consume the same type of food. In an economic system this term has to do with the saturation due to the finite size of the economy.

To simplify the analysis we will consider in this paper a simple case in which the $`w_i`$’s interact in a uniform fashion with each other. This case is obtained by choosing $`a_{i,j}=a/N`$ and $`c_{i,j}=c/N`$. With this choice Eq. (4) will be reduced to

$`w_i(t+1)`$ $`=`$ $`[1+\lambda (t)]w_i(t)+a\overline{w}(t)cw_i(t)\overline{w}(t)`$
$`w_j(t+1)`$ $`=`$ $`w_j(t),j=1,\mathrm{},N;ji.`$ (5)
where

$$\overline{w}(t)=\frac{1}{N}\underset{i=1}{\overset{N}{}}w_i(t).$$
(6)
is the average value of the dynamical variables at time t. Here the random term $`\lambda (t)`$ was shifted to $`\lambda (t)a`$ but its distribution around the average (and thus the value of the standard deviation $`D`$) remain unchanged. The second term in Eq. (5), may now describe the effect of auto-catalysis at the community level. In an economic model, this term can be related to the social security policy or to general publicly funded services which every individual receives. It prevents an individual $`w_i`$ from falling below a certain fraction of the average $`\overline{w}`$. The third term in Eq. (5), describes saturation or the competition for limited resources. It has the effect of limiting the growth to values sustainable for the current conditions and resources. Within the ecological context, here the interactions between populations are uniform, describing the case in which all of them consume the same type of food. We refer to Eq. (5) as the generalized Lotka-Volterra system because when averaged over $`i`$ and over $`\lambda (t)`$, this system tends to approach a Lotka-Volterra-like equation Lotka1925 ; Volterra1926

$$w(t+1)=(1+\lambda +a)w(t)cw^2(t).$$
(7)
where $`w(t)\overline{w}(t)`$. Computer simulations show that after some equilibration time the system described by Eq. (5) approaches steady-state conditions. Even at steady state, $`\overline{w}`$ exhibits fluctuations. However, its average over long time scales approaches a constant value given by

$$\overline{w}_t=(\lambda +a)/c.$$
(8)
In previous studies of the system described by Eq. (5) the parameters $`a`$ and $`c`$ were considered as constants, corresponding to steady conditions of the market. In fact, the typical dynamics of microscopic market models Levy1994 ; Levy1995 ; Solomon1996 ; Levy1997  is generically not in a steady state. The effect of varying market conditions can be studied by considering the parameters $`a`$ and $`c`$ and the distribution $`\mathrm{\Pi }(\lambda )`$ as slowly varying functions of time. We will show below that systems described by Eq. (5) lead, under very general conditions, to a power-law distribution of the $`w_i`$’s of the form of Eq. (1). Moreover, it will be shown that the exponent $`\alpha `$ is insensitive to variations in the parameter $`c`$, namely it depends only on $`a`$ and $`D`$. In order to examine the effect of variations in the economic conditions we will now introduce an explicit time dependence into the third term, as well as a more general dependence on the $`w_j`$’s. The dynamic equation will now take the form

$`w_i(t+1)`$ $`=`$ $`[1+\lambda (t)]w_i(t)+a\overline{w}(t)C(w_1,\mathrm{},w_n,t)w_i(t)`$
$`w_j(t+1)`$ $`=`$ $`w_j(t),j=1,\mathrm{},N;ji,`$ (9)
where $`C(w_1,\mathrm{},w_n,t)`$ is a general function of the $`w_j`$’s that includes an explicit time dependence.

## III Theoretical Analysis

In order to study the dynamics of the generalized Lotka-Volterra model, it will be convenient to denote the change of $`w_i`$ in a single time step by $`\mathrm{\Delta }w_i(t)=w_i(t+1)w_i(t)`$. We introduce a set of normalized variables

$$x_i=\frac{w_i}{\overline{w}},i=1,\mathrm{},N.$$
(10)
The change $`\mathrm{\Delta }x_i(t)=x_i(t+1)x_i(t)`$ in a single time step is given by

$$\mathrm{\Delta }x_i\frac{1}{\overline{w}}\mathrm{\Delta }w_i\frac{w_i}{\overline{w}^2}\mathrm{\Delta }\overline{w}$$
(11)
up to first order in powers of the $`\mathrm{\Delta }w_i`$’s. Considering the time dependence of the average $`\overline{w}`$ one should remember that at any time step $`t`$ of Eq. (9) only one of the $`w_i`$’s is chosen randomly and updated. Moreover, there is no correlation between the chosen $`w_i`$ and $`\lambda (t)`$. Thus, the time evolution of $`\overline{w}`$ should be considered on a longer time scale of the order of $`N`$ moves. However, for simplicity we evaluate $`\mathrm{\Delta }\overline{w}`$ by averaging Eq. (9) for $`\mathrm{\Delta }w_i`$ over $`i=1,\mathrm{},N`$ at a given time $`t`$. We make an independent random choice of $`\lambda (t)`$ for each $`i`$, which we denote by $`\lambda _i(t)=\lambda +\eta _i(t)`$. The time dependence of $`\overline{w}`$ is given by

$$\mathrm{\Delta }\overline{w}=\frac{1}{N}\underset{j=1}{\overset{N}{}}\eta _j(t)w_j(t)+\left[\lambda +a\right]\overline{w}(t)C(w_1,\mathrm{},w_n,t)\overline{w}(t).$$
(12)
The dynamics of the $`x_i`$’s is thus given by

$$\mathrm{\Delta }x_i=x_i(t)\left[\eta (t)a\frac{1}{N}\underset{j}{}\eta _j(t)x_j(t)\right]+a.$$
(13)
Consider the sum

$$r(N)=\underset{j=1}{\overset{N}{}}\eta _j(t)x_j(t).$$
(14)
If the $`x_j`$’s exhibit a distribution of the form

$$P(x)x^{1\alpha },$$
(15)
then the second moment of the distribution of $`r(N)`$ satisfies Shlesinger1993 ; Shlesinger1994

$$r^2(N)^{1/2}=\{\begin{array}{ccc}N^{1/2},\hfill & 2<\alpha \hfill & \\ N^{(3\alpha )/2},\hfill & 1<\alpha <2\hfill & \\ N,\hfill & 0<\alpha <1.\hfill & \end{array}$$
(16)
In the first case, the distribution of the $`x_j`$’s exhibits a finite second moment and $`r(N)/N0`$ in the limit $`N\mathrm{}`$. In the second case the second moment of $`P(x)`$ diverges and $`r(N)`$ follows a Lévy distribution.

In both cases, namely for $`\alpha >1`$, we obtain that in the (thermodynamic) limit $`N\mathrm{}`$:

$$\frac{1}{N}x_i\underset{j}{}\eta _j(t)x_j(t)0.$$
(17)
Thus, under the assumption that $`P(x)`$ follows Eq. (15) with $`\alpha >1`$, we obtain to a good approximation that for large values of $`N`$

$$\mathrm{\Delta }x_i=\left[\eta (t)a\right]x_i+a,i=1,\mathrm{},N.$$
(18)
We see that the dynamics of the normalized variable $`x_i`$ is reduced to a set of identical decoupled linear Langevin equations, which do not depend on the function $`C(w_1,\mathrm{},w_n,t)`$ or on the mean value $`\lambda `$ of the multiplicative noise. These equations can be cast into a general framework of multiplicative processes of the form

$$\mathrm{\Delta }x(t)=\eta (t)G(x(t))+F(x(t)).$$
(19)
Eq. (18) can then be recovered by taking $`F(x_i)=a(1x_i)`$ and $`G(x_i)=x_i`$. By using a suitable change of variables to $`y=y(x)`$ that satisfies

$$\frac{dy}{dx}=\frac{1}{G(x)}$$
(20)
one can reduce Eq. (19) to a Langevin equation in which the term $`\eta (t)`$ appears as an additive noise, rather than a multiplicative noise such as $`\eta (t)G(x(t))`$ Richmond2001 . The time evolution of $`y(t)`$ is obtained from Eq. (19) by using the chain differential rule up to second order in $`\mathrm{\Delta }x`$ (and first order in $`D`$)

$$\mathrm{\Delta }y\frac{dy}{dx}\mathrm{\Delta }x+\frac{1}{2}\frac{d^2y}{dx^2}\mathrm{\Delta }x^2.$$
(21)
Inserting $`\mathrm{\Delta }x(t)`$ from Eq. (19) and using the change of variables described in Eq. (20) we obtain

$$\mathrm{\Delta }y\frac{1}{G}(F+\eta G)+\frac{1}{2}\frac{d}{dx}(\frac{1}{G})(F+\eta G)^2.$$
(22)
We now approximate the second order term by averaging over the noise term $`\eta `$, that satisfies $`\eta =0`$ and $`\eta ^2=D`$. We obtain

$$\mathrm{\Delta }y\eta +\frac{F}{G}\frac{1}{2}\frac{dG}{dx}(D+\frac{F^2}{G^2}).$$
(23)
Assuming that $`F^2/G^2D`$, Eq. (23) is reduced to a discrete-time Langevin equation

$$\mathrm{\Delta }y\eta +J(y)$$
(24)
where the drift force $`J(y)`$ takes the form

$$J(y)=\frac{F}{G}\frac{D}{2}\frac{dG}{dx}.$$
(25)
The Fokker-Planck equation corresponding to Eq. (24) is vanKampen1981

$$\frac{P(y,t)}{t}=\frac{}{y}\left[J(y,t)P(y,t)\right]+\frac{D}{2}\frac{^2P(y,t)}{y^2}$$
(26)
where $`P(y,t)`$ is the probability distribution of $`y`$ at time $`t`$. The solution of this equation under the stationary condition $`P(y,t)/t=0`$ is

$$P(y)=\mathrm{exp}\left[\frac{2}{D}^yJ(y^{})𝑑y^{}\right].$$
(27)
Thus, the distribution $`P(x)=P(y)dy/dx`$ of the original variables $`x_i`$ is

$$P(x)=\frac{1}{G^2(x)}\mathrm{exp}\left[\frac{2}{D}^x\frac{F(x^{})}{G^2(x^{})}𝑑x^{}\right].$$
(28)
By taking $`F(x)=a(1x)`$ and $`G(x)=x`$ we obtain a power-law distribution of the form Solomon2001

$$P(x)=x^{1\alpha }\mathrm{exp}\left[\frac{1\alpha }{x}\right].$$
(29)
where

$$\alpha =1+\frac{2a}{D}.$$
(30)
The exponent $`\alpha `$ thus depends on a single parameter $`a/D`$, namely on the ratio between the global drift coefficient $`a`$ and the fluctuations measured by $`D`$.

Another way to derive Eq. (30) from dynamical models of the form (18) was shown in Refs. Takayasu1997 ; Sornette1997 . It is based on the fact that, under steady-state conditions, linear Langevin equations of the form (18) satisfy Takayasu1997 ; Sornette1997 ; Kesten1973

$$(\eta a+1)^\alpha =1.$$
(31)
Considering $`\eta a`$ as a small parameter and expanding Eq. (31) in a power series up to second order we obtain

$$\alpha =1+\frac{2a}{D+a^2}.$$
(32)
Assuming that $`a^2D`$ we reproduce Eq. (30).

## IV Numerical Simulations and Results

To examine the theoretical predictions presented in Sec. III we have performed computer simulations of the generalized Lotka-Volterra system described by Eq. (9) with different choices of $`C(w_1,\mathrm{},w_n,t)`$. It was found that after some equilibration time the distribution $`P(x)`$ reaches a steady state, and exhibits a power-law behavior. For $`C(w_1,\mathrm{},w_n,t)=c\overline{w}`$ \[as in Eq. (5)\], $`\overline{w}(t)`$ fluctuates around some average value, given by Eq. (8). For a general function $`C(w_1,\mathrm{},w_n,t)`$, that exhibits an explicit time dependence, $`\overline{w}(t)`$ continues to vary according to this function and its temporal average does not reach a steady state.

To examine how robust the power-law distribution is under varying conditions we have simulated Eq. (9) with

$$C(w_1,\mathrm{},w_n,t)=c_o\left(1+\mathrm{sin}\frac{2\pi t}{T}\right)\underset{j=1}{\overset{N}{}}w_j^2$$
(33)
(for $`c_0=0.001`$, $`T=2\times 10^5`$ and $`\lambda =0.002`$) and compared the results with the case $`C(w_1,\mathrm{},w_n,t)=c\overline{w}`$ (for $`c=1`$ and $`\lambda =0.01`$). The time dependence of $`\overline{w}`$ in both cases is shown in Fig. 1(a). The distributions $`P(x)`$ obtained from the simulations in these two cases are shown in Fig. 1(b). The distributions are found to be nearly identical and exhibit a power-law behavior characterized by the same exponent $`\alpha `$. The exponent $`\alpha `$ is also found to be independent of $`\lambda `$. Note that the power-law behavior is maintained even for $`C(w_1,\mathrm{},w_n,t)0`$, where $`\overline{w}(t)`$ does not reach a steady state and diverges to infinity (for $`\lambda >0`$) or collapses to $`0`$ (for $`\lambda <0`$) Malcai1999 . This can lead to changes by orders of magnitude in the total wealth or the population size without affecting the exponent $`\alpha `$.

To examine the theoretical prediction for the distribution, given by Eq. (29), and the exponent $`\alpha `$, given by Eq. (30), we have compared these predictions to the results of numerical simulations. This comparison for the distribution $`P(x)`$ is shown in Fig. 2. The simulations were done for $`N=1000`$, $`a=0.00023`$, $`C(w_1,\mathrm{},w_n,t)=0.01\overline{w}`$, and $`\lambda (t)`$ uniformly distributed in the range $`0.0\lambda (t)0.1`$ (namely, $`D=0.00083`$). We found that the simulation results are in very good agreement with the theoretical predictions. A power law distribution is found for a range of almost three orders of magnitude, with the exponent $`\alpha =1.52`$. This is close to the theoretical prediction of $`\alpha =1+2a/D=1.55`$. The distribution has a peak at $`x_0=(\alpha 1)/(\alpha +1)`$, that using Eq. (30) can be expressed by $`x_0=a/(a+D)`$. Above $`x_0`$ the distribution $`P(x)`$ behaves like power-law while below it $`P(x)`$ decays exponentially. This provides an effective lower cutoff for the range of $`x`$ in which a power-law behavior is observed. This result can be compared to a somewhat simpler model studied earlier, in which the value of the lower cutoff $`x_{min}`$ is imposed as a constraint Malcai1999 . In this model, using the sum rules for the probability and the total wealth, it was found that $`x_{min}=11/\alpha `$. Using Eq. (30) it can be expressed as $`x_{min}=2a/(2a+D)`$. These predictions for the lower bounds in the two models satisfy $`x_0<x_{min}<2x_0`$, namely, they are in good agreement in light of the broad distribution of $`x`$.

To examine the prediction given by Eq. (30) for the exponent $`\alpha `$, we present in Fig. 3 a comparison between this prediction and the numerical results for $`\alpha `$ as a function of $`a/D`$. The numerical results are presented for $`N=1000`$ ($``$). The prediction of Eq. (30) (solid line), shows a good agreement with the numerical results for $`a/D>0.2`$. The numerical results for the range of small $`a/D`$ converge to the theoretical prediction as the value of $`N`$ is increased. As shown in Ref. Malcai1999  the infinite system limit, $`N\mathrm{}`$, and the vanishing coupling limit, $`a/D0`$, do not commute. On the one hand, for any finite $`N`$ and $`a/D0`$ the exponent $`\alpha 0`$. On the other hand, for any fixed positive value of $`a/D`$ (no matter how small) and $`N\mathrm{}`$ the exponent $`\alpha 1`$. The majority of empirical results in which $`1<\alpha <2`$, indicate that the second case is highly relevant and that the theoretical predictions of Eqs. (29) and (30) broadly apply. Thus, for $`N\mathrm{}`$, both the exponent $`\alpha `$ of the power-law decay, and the lower bound $`x_0`$ depend only on a single parameter $`a/D`$. In the economic context, this parameter represents the ratio between the fixed income of minimal wage jobs or social security payments and the level of fluctuations of the speculative income/loss.

## V Summary

We have studied the dynamics of stochastic Lotka-Volterra systems under non-stationary conditions using both analytical and numerical techniques. For this class of models, we found that in order to obtain a power-law distribution, it is sufficient that relative returns of the agents are stochastically equivalent. The assumption that the distribution $`\mathrm{\Pi }(\lambda )`$ of the multiplicative noise, is independent of $`i`$, means that there are no investors or strategies that can obtain ’abnormal’ returns. This can be related to to the ’efficient market hypothesis’, which assumes that the market pricing mechanism is so efficient that it reaches the ’right price’ before any of the agents can take systematic advantage. Therefore, the presence of a power-law distribution may be a sign of ’market efficiency’, by analogy with Boltzmann distributions in statistical mechanics systems, which characterize thermal equilibrium. Here we have shown that the power-law distribution is stable even under non-stationary economic conditions, that are represented by the time dependence of the saturation term $`C(w_1,\mathrm{},w_n,t)`$. We found that even under such conditions the distribution of the (normalized) dynamical variables $`x_i`$ follow a power-law distribution with an exponent $`\alpha `$. An expression for $`\alpha `$ in terms of ratio of the parameters $`a`$ and $`D`$ was obtained \[Eq. (30)\]. In the economic context, the parameter $`a`$ represents the minimal wage or social security payments, while $`D`$ represents the level of fluctuations in speculative income/loss. These results provide the distribution of wealth in a society in terms of the social security policy and the volatility of the stock market. They also provide a connection between the incomes/wealths of the poorest and the richest sectors of the society as a function of a single parameter.


text:
# How the BAL quasars are quiet.

## 1 Introduction

Recent observations show two discrepancies in the standard interpretation of the spectra of quasars:

i) the relative Doppler or expansion frequency shift of spectral lines is strictly constant, while the observations show slight variations ; the explanation by a variation of the fine structure constant is hardly credible;

ii) the concentration of iron in quasars whose spectrum is largely shifted requires a generation of this metal in old stars , implying a too large increase of the age of the Universe.

The standard model seems unable to explain that the loud quasars do not have broad lines.

Twenty years ago, it seemed difficult to interpret various observed redshifts by the Doppler effect , but the lack of a credible alternative led to find solutions which do not convince all astrophysicists.

In 1968, Giordmaine et al.  observed frequency shifts in experiments using high power ultrashort laser pulses; this observation is now common, used for analysis and called ”Impulsive Stimulated Raman Scattering” (ISRS). The theory of ISRS shows that this effect has no intensity threshold . The intensities of the lasers make ISRS non-linear, so that the frequency shifts depend on these intensities; at low intensities, the effect becomes linear, so that it is similar to a Doppler effect, except for a small dispersion. This qualitative difference justifies a new name: ”Coherent Raman Effect on time-Incoherent Light” (CREIL).

In the next section, the properties of the CREIL are recalled. In the following sections, the influence of each of these properties on the standard model of quasars is given, solving very simply some problems, in particular the noisy quasars cannot show broad spectral lines.

## 2 Properties of the ”Coherent Raman Effect on Incoherent Light” (CREIL).

Any alternative to the Doppler effect in the vacuum requires a time-incoherence of the light: If, during a certain time, the received frequency of a coherent light is lower than the emitted frequency, the number of cycles, that is the number of wavelengths between the source and the receiver is increased; as the wavelength is an unit of length, the distance is increased, it is a Doppler effect. Thus, we understand that the length of the pulses of light is critical, the pulses of light must verify the condition set by G. L. Lamb when the length of the pulses is a condition for the observation of an optical effect: the ”ultrashort light pulses”, must be ”shorter than all relevant time constants” .

In an homogeneous gas, two relevant time constants must be longer than the pulses which make the usual incoherent light, that is longer than about 10 nanoseconds:

i) the collisional time $`\tau `$ must be longer than 10 ns, for instance $`\tau =100ns`$. For hydrogen atoms at 10000K, the density of atoms $`N`$ is $`N1,7.10^{21}`$ atomes/m3; for H2 at 100K, $`N2,4.10^{22}`$; these values are orders of magnitude because the mean diameter of the atoms depends on the criterion used to define a collision. This condition allows that all molecules on a wave surface remain in phase and scatter coherent waves which rebuild wave surfaces identical to the incident wave surfaces: no blur of the images.

ii) The period corresponding to the Raman transition must be larger than 10 ns, that is in the radiofrequencies range. As the period of the beats of the incident and scattered beams is much longer than the duration of the pulses, the beats do not appear, the incident and scattered beams interfere into a single frequency beam; the frequency shift is proportional to the ratio of the amplitudes and to the Raman frequency.

Although the coherent scattering is much more powerful than an incoherent Raman scattering, the previous conditions make CREIL much weaker than ISRS which uses dense matter and infrared Raman frequencies; an other source of weakness is that the populations of the levels involved in the Raman transition are nearly the same because they are close, so that the opposed virtual transitions nearly cancel the shifts that they produce. Thus, it seems impossible to get in the labs long enough paths to observe the CREIL while the Universes provides long paths in low pressure gases.

Previous papers  show that the cosmological redshift may be produced by CREIL if there are, in the average, some active molecules per litre in the intergalactic space; it is a very rough evaluation because a precise computation would require for each molecule present in the space, the knowledge of a lot of traces of Raman polarisation tensors.

Suppose that a redshifting gas has absorption lines; the absorption by a line occurs from the redshifted frequency to the absolute frequency; if the redshift is large, the line is wide and weak. If the gas has many absorption lines, the lines are mixed, they cannot be observed.

If ISRS and CREIL were simple Raman effects, the active molecules would be excited, and, without collisions, their de-excitations would be slow; ISRS and CREIL are parametric effects which mix the described Raman effect with an other which de-excite the molecules and blue-shift (an) other beam(s); in ISRS, the other beam is a laser beam, in CREIL, its role may be plaid by the thermal radiation which is amplified. In this last case, all frequencies are low, resonant.

## 3 Distance (age) of the quasars

The distance and the age of the quasars is deduced from the redshift of the sharp emission or absorption lines, which are the most redshifted lines. It was remarked that the neighbourhood of the most redshifted quasars contains a lot of matter, in particular molecular hydrogen detected by the 21 cm nuclear spin coupling line. It is reasonable to think that where the pressure is low, the molecular hydrogen is partly ionised by the UV into H$`{}_{}{}^{+}{}_{2}{}^{}`$. This gas cannot be detected by its absorption or emission lines because it has many relatively weak lines which are shifted during their absorption or emission by the CREIL; if the pressure is high enough to produce collisions which forbid the CREIL, these collisions destroy H$`{}_{}{}^{+}{}_{2}{}^{}`$ too quickly to allow its observation.

The CREIL distorts the standard scale of distances deduced from the redshifts. This scale remains probably good in the intergalactic space, because there, the redshift is produced by expansion or by the CREIL in a nearly homogenous medium; but, where there is much molecular hydrogen, in particular in the neighbourhood of the quasars the distances must be strongly reduced.

## 4 Absence of broad lines in the loud quasars

Suppose that a quasar has an accretion disk larger, but similar to Saturn’s disk, except maybe for an inner part which protects the remainder from burning. The friction of the stones on the halo of the quasar charges them, just as the drops of rain are charged in our clouds. While the lightnings discharge a volume of our clouds, the discharges of regions of the disk are limited to its plane, so that the electromagnetic field radiated by sheet-shaped currents propagate more in the direction of the axis of the disk.

The disk is only partly and locally discharged by the lightnings; as the disk is inhomogeneous, the sign of its charges depend on the distance to the quasar; thus, the rotation of the disk produces circular currents which induce a magnetic field similar to the field introduced by de Kool and Begelman: near the disk, its orientation may change, so that it may have zero values, in particular if the disk has gaps.

Raman transitions inside a Zeeman structure of the hydrogen atom require a principal quantum number larger than one, that is the atom must have been excited by a Lyman absorption; then the excited atom is active for CREIL were the magnetic field is not null; for atoms pumped by a Lyman $`\alpha `$ transition, neglecting the electronic spin, and supposing that the de-exciting Raman component is very strong, the computation of the CREIL requires only the trace of a single, simple Raman polarisability tensor. If the Lyman lines are redshifted while they are absorbed, they appear so wide and weak that they cannot be seen.

If a line of sight makes large enough an angle with the axis of the quasar, it crosses places where the magnetic field is nearly zero, so that, as there is no CREIL in these places, the absorption of the Lyman lines is stable in the spectrum and the lines appear. The large width of the lines may be a consequence of a low CREIL (by a low magnetic field, or collisional decrease), or, as this process is close to the quasar, it may depend slightly on the line of sight. This description of the origin of the BAL applies to the BEL.

Thus, the loudness, the presence of broad lines is not a property of the quasars, but depends on the direction of observation: BAL quasars, observed in the direction nearly perpendicular to the axis of the disk, cannot radiate much electromagnetic field in this direction so that they are ”quiet”.

## 5 The Lyman forest

In the standard model, the Lyman forest corresponds to clouds, that is to variations of the density of hydrogen atoms. It may be the result of the absorption in a chemically homogenous medium and a space-variable redshifting power of the gas resulting of variations of the magnetic (or electric) field. But how could a space-variable magnetic field be generated in the far halo? Maybe by electric instabilities in the plasma, or it could be the magnetic part of an extremely low frequency electromagnetic field radiated by the quasar and propagating very slowly in the plasma, more probably by magnetised satellites of the quasar.

In one of these hypotheses, the modulation of the absorption which structures the Lyman forest is self-amplifying: Follow the light; if, in a region, the magnetic field is very low, there is no redshift, the Lyman $`\alpha `$ line saturates, so that, in the following region there is no Lyman $`\alpha `$ pumping, thus no CREIL even if a magnetic field appears; a high magnetic field is necessary to shift the spectrum by a small number of remaining CREIL-active atoms ( hydrogen pumped by other Lyman lines, and other atoms). When the line is shifted of a fraction of its linewidth, the Lyman $`\alpha `$, the CREIL and the redshift become strong, the absorption is spread and seems low in the spectrum; a return to the saturation requires a low magnetic field. This discussion is valid where the forest is saturated, but its results remain true elsewhere: only small fluctuations of the field may induce the forest.

Two effects explain that the density of lines in the high redshift region of the forest generally decreases: i) the density of hydrogen being higher and its illumination higher, the CREIL is larger, the spectrum is spread, so that the decrease is an artefact; ii) in this region, some small satellites may fall into the accretion disk.

Along a line of sight, the frequencies are shifted near the magnetic satellites while the Lyman lines are written far from these satellites. When an eclipse is observed by a variation of the intensity of the quasar, if small variations of the redshift can be detected, the gap corresponding to the satellite in the Lyman forest may be found. It could help a study of the satellite and the quasar.

## 6 High redshift quasars

The high redshift quasars appear in regions of the Universe in which clouds of molecular hydrogen are detected by the nuclear spin coupling transition at 21 cm. They seem dusty too.

If the pressure of molecular hydrogen allows a very long collisional time, the UV ionises a part of the hydrogen into H$`{}_{}{}^{+}{}_{2}{}^{}`$ which is active in CREIL, thus cannot be observed from its absorption or emission spectra. The CREIL increases the redshift of the quasar; it amplifies the thermal radiation, possibly up to 100K, so that this thermal radiation may be confused with radiation of hot dust.

If a galaxy appears next to a quasar, and has a lower redshift, it may, however, be further than the quasar. Thus its line of sight may cross the neighbourhood of the quasar, so that absorptions of metal lines may have the same origin, appear with the same redshift.

## 7 Conclusion

This draft intends only to show that taking the Coherent Raman Effect on time-Incoherent Light (CREIL) into account solves very easily, without any new physics or new matter, problems which appear in the study of the quasars: \- the absence of BAL lines in the loud quasars is deduced from the hypothesis of a disk similar to Saturn’s; \- the discrepancies in the observed spectra may be produced by the small dispersion of the CREIL; \- the quasars are younger than they seem, so that they can contain metals made in old stars.

The use of CREIL allows to find several explanations of a particular observation; the present qualitative study must be specified quantitatively to choose the best explanation, maybe the standard explanation: for instance the thermal radiation may be partly produced by hot dust, partly by CREIL. Having several possible solutions is better than having none !

The model of quasars which takes CREIL into account uses only commonly observed types of astrophysical objects and simple regular physics; with much less work, it seems to give interpretations of observations on the quasars better than the standard model; it does not require dark matter, but if this matter is necessary for the gravitational stability, it proposes partially ionised H2. Logically, is it acceptable to neglect CREIL a priori?

The present demonstration of the necessity to take CREIL into account in astrophysics weakens the two main proofs of the big bang, expansion and 2.7K radiation.