text:
# Constructing, characterizing, and simulating Gaussian and higher–order point distributions

## I Introduction

Stochastic models are often used to describe physical phenomena. For spatial structures two broad classes of stochastic models have been established. One approach is based on random fields and the other one on random distributions of discrete objects, often only points, in space (see the contributions in  for recent applications and reviews). Stochastic point distributions are used to describe physical systems on vastly differing length–scales. The physical applications discussed in this article deal with the large–scale structures in the Universe formed by the distribution of galaxies. However, the methods are much more versatile.

Models for the dynamics of cosmic structures are often based on nonlinear partial differential equations for the mass density and velocity field. These models relate the initial mass density and velocity field, primarily modeled as Gaussian random fields, to the present day values of these fields. The nonlinear evolution leads to non–Gaussian features in the fields. However, observations supply us with the positions of galaxies in space. To compare theories with observations one has to relate fields with point distributions. Both deterministic or stochastic models have been used for this purpose so far (e.g. ).

Pursuing a direct approach, the observed spatial distribution of galaxies (galaxy clusters etc.) is compared to models for random point sets. Only a few attempts towards a dynamics of galaxies as discrete objects have been made (see e.g. ), however stochastic models are quite common. Following the works by , and  a purely stochastic description of the spatial distribution of galaxies as points in space is given in this article.

Models for stochastic point processes can be constructed using the physical interactions of the objects, typically leading to Gibbs processes (see e.g. , and the generalizations by ). Another approach to construct point processes is based on purely geometrical considerations, e.g. points are randomly distributed on randomly placed line–segment (see ). As a third possibility one can start from the characterization of point processes by the probability generating functional (p.g.fl.) and the expansion in terms of correlation functions. This is the way pursued in this work.

The simplest point process is a Poisson process showing no correlations at all. Since the galaxy distribution is highly clustered, a Poisson process is not a realistic model. The model with the next level of sophistication is a Gauss–Poisson process, the point distribution counterpart of a Gaussian random field. Whereas the properties of Gaussian random fields have been extensively studied, the Gauss–Poisson process has not been discussed in the cosmological literature in a systematic way. Some of the statistical properties of random fields directly translate to similar statistical properties of point distributions, but also important differences show up. The systematic inclusion of higher–order correlations, as well as the characterization, and the simulation algorithms for such point processes will be discussed.

Recently, a related class of stochastic models for the galaxy distribution, the halo–model, attracted some attention (see e.g. ). These models are based on the assumption that galaxies are distributed inside correlated dark matter halos. Using the probability generating functional, the two– and three–point correlation function will be calculated for this model, extending the results by  to include the effects of halo–substructure.

The outline of this paper is as follows:
In Sect. II the properties of the probability generating functional (p.g.fl.) of a point process and the expansions in several types of correlation functions are briefly reviewed. The characterization of the Gauss–Poisson process is given in Sect. III, and the physical consequences of the constraints are discussed. The close relation to Poisson cluster) processes allows us to simulate a Gauss–Poisson process (see Appendix B 1. In Sect. IV simulations of the Gauss–Poisson processes and the line–segment process are used to show how the Gaussian approximation influences the $`J(r)`$-function, a statistic sensitive to higher–order correlations. A comparison of the galaxy distribution within the PSCz survey with a Gauss–Poisson processes illustrates the importance of higher–order correlations. In Sect. V the extension of the Gauss–Poisson point process to the $`n`$–point Poisson cluster process is discussed. Detailed results are derived for the three–point Poisson cluster process (the simulation recipe is give in Appendix B 2). The characterization of the general $`n`$–point process is discussed which is again detailed for the three–point case. Differences between a point process and a random field are highlighted in Sect. VI. Models for strongly correlated systems are mentioned in Sect. VII. The focus will be on the “halo model”. Using the formalism based on the p.g.fl., the correlation functions of the “halo model”, including the effects of halo–substructure, are calculated in Sect. VII B. In Sect. VIII some open problems are mentioned. An outlook is provided in Sect. IX. As an example the probability generating function (p.g.f.) of a random variable and its expansions in several kinds of moments is reviewed in Appendix A.

## II Product densities, factorial cumulants, and the probability generating functional

Probability generating functionals (p.g.fl.’s), and their expansions in different kinds of correlation measures have been used to describe noise in time series (e.g. ) and the electro–magnetic cascades occurring in air–showers (e.g. ). They have been employed in the theory of liquids (e.g. ) and other branches of many–particle physics (e.g. ). The mathematical theory of p.g.fl.’s for point processes is nicely reviewed in the book of . Stochastic methods based on p.g.f.’s have been introduced to cosmology by  (the p.g.fl. was presented by  in the discussion of this article), and became well–known following the work of  and . Focusing on the factorial moments (the volume averaged $`n`$–point densities) and on count–in–cells,  discussed several expansions of the p.g.f.’s. In the following only “simple” point processes will be considered: at each position in space at most one object is allowed. This assumption is physically well justified for galaxies. Also, for quantum systems the methods should be refined (see e.g.  for fermionic (determinantal) point processes).

An intuitive way to characterize a point process is to use $`n`$th–product densities: $`\varrho _n(𝐱_1,\mathrm{},𝐱_n)\mathrm{d}V_1\mathrm{}\mathrm{d}V_n`$ is the probability of finding a point in each of the volume elements $`\mathrm{d}V_1`$ to $`\mathrm{d}V_n`$. For stationary and isotropic point fields $`\varrho _1(𝐱)=\varrho `$ is the mean number density, and the product density (with a slight abuse of notation) is $`\varrho _2(𝐱_1,𝐱_2)=\varrho _2(r)`$ with $`r=|𝐱_1𝐱_2|`$ being the separation of the two points. The factorial cumulants $`c_{[n]}(𝐱_1,\mathrm{},𝐱_n)`$ are the irreducible or connected parts of the $`n`$th–product densities. E.g. for $`n=2`$

$$\varrho _2(r)=\varrho ^2+c_{[2]}(r)=\varrho ^2\left(1+\xi _2(r)\right),$$
(1)
and the second factorial cumulant $`c_{[2]}(r)`$ and the two–point correlation function $`\xi _2(r)`$ quantify the two–point correlations in excess of Poisson distributed points.

A systematic characterization of a point process is provided by the probability generating functional or a series of probability generating functions (see Appendix A). Equivalent to a random distribution of points in space, one considers a point process as a random counting measure. A realization is then a counting measure $`N`$, which assigns to each suitable set $`A^d`$ the number of points $`N(A)`$ inside. For suitable functions $`h(𝐱)`$ one defines the probability generating functional of a point process via

$$G[h]=𝔼\left[\mathrm{exp}\left(_^dN(\mathrm{d}𝐱)\mathrm{log}h(𝐱)\right)\right],$$
(2)
where $`^d`$ is the $`d`$–dimensional Euclidean space, and $`𝔼`$ is the expectation value, the ensemble average over realizations of the point process. Equivalently,

$$G[h]=𝔼\left[_ih(𝐱_i)\right],$$
(3)
where $`𝐱_i`$ are the particle positions in a realization. Consider $`k`$ compact disjoint sets $`A_j`$, and let $`n_j=N(A_j)`$ be the number of points inside $`A_j`$. The p.g.f. of the $`k`$-dimensional random vector $`(n_1,\mathrm{},n_k)`$ is then

$$P_k(𝐳)=P_k(z_1,\mathrm{},z_k)=𝔼\left[\underset{j=1}{\overset{k}{}}z_j^{n_j}\right].$$
(4)
Together with a continuity requirement the knowledge of all finite dimensional p.g.f.’s $`P_k`$ determines the p.g.fl. $`G`$ and the point process completely (e.g. ). One obtains the p.g.f. $`P_k(𝐳)=G[h^{}]`$ of the random vector $`𝐳`$ using

$$h^{}(𝐱)=1\underset{j=1}{\overset{k}{}}(1z_j)\text{1l}_{A_j}(𝐱),$$
(5)
where $`\text{1l}_A(𝐱)`$ is the indicator–function of the set $`A`$, with $`\text{1l}_A(𝐱)=1`$ for $`𝐱A`$ and zero otherwise. Several expansions of the p.g.fl. $`G[h]`$ are possible . The expansion in terms of the product densities $`\varrho _n`$ (the Lebesgue densities of the factorial moment measures) reads

$$\begin{array}{c}G[h+1]=1+\underset{n=1}{\overset{\mathrm{}}{}}\frac{1}{n!}_^dd𝐱_1\mathrm{}_^dd𝐱_n\hfill \\ \hfill \varrho _n(𝐱_1,\mathrm{},𝐱_n)h(𝐱_1)\mathrm{}h(𝐱_n).\end{array}$$
(6)
For the factorial cumulants $`c_{[n]}`$ or correlation functions $`\xi _n`$ one obtains ***The relations to the generating functionals $``$, $``$ and $`𝒢`$ defined by  are $`G[h]=[h]`$, $`G[h]=[h1]`$ and $`G[h]=\mathrm{exp}𝒢[h1]`$.

$`\mathrm{log}G[h+1]`$ $`={\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{n!}}{\displaystyle _^d}d𝐱_1\mathrm{}{\displaystyle _^d}d𝐱_n`$ (7)
$`c_{[n]}(𝐱_1,\mathrm{},𝐱_n)h(𝐱_1)\mathrm{}h(𝐱_n)`$ (8)
$`={\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{\varrho ^n}{n!}}{\displaystyle _^d}d𝐱_1\mathrm{}{\displaystyle _^d}d𝐱_n`$ (9)
$`\xi _n(𝐱_1,\mathrm{},𝐱_n)h(𝐱_1)\mathrm{}h(𝐱_n).`$ (10)
As a third possibility the p.g.fl. can be expanded around the origin:

$$\begin{array}{c}G[h]=J_0+\underset{n=1}{\overset{\mathrm{}}{}}\frac{1}{n!}_^dd𝐱_1\mathrm{}_^dd𝐱_n\hfill \\ \hfill j_n(𝐱_1,\mathrm{},𝐱_n)h(𝐱_1)\mathrm{}h(𝐱_n).\end{array}$$
(11)
The Janossy densities $`j_n(𝐱_1,\mathrm{},𝐱_n)\mathrm{d}V_1\mathrm{}\mathrm{d}V_n`$ are the probability that there are exactly $`n`$ points, each in one of the volume elements $`\mathrm{d}V_1`$ to $`\mathrm{d}V_n`$. Convergence issues of these expansions are discussed in Sect. VIII.

## III The Gauss–Poisson point process

For a stationary Poisson process with mean number density $`\varrho `$ the p.g.fl. is

$$\mathrm{log}G[h+1]=\varrho _^dd𝐱h(𝐱),$$
(12)
corresponding to a truncation of the expansion (7) after the first term. Truncating after the second term, one obtains the p.g.fl. for the Gauss–Poisson process

$$\begin{array}{c}\mathrm{log}G[h+1]=\varrho _^dd𝐱h(𝐱)+\hfill \\ \hfill +\frac{\varrho ^2}{2}_^dd𝐱_^dd𝐲\xi _2(|𝐱𝐲|)h(𝐱)h(𝐲),\end{array}$$
(13)
completely specified by its mean number density $`\varrho `$ and the two-point correlation function $`\xi _2(r)`$.

There is a close resemblance to random fields. For a homogeneous random field $`\rho (𝐱)`$ with mean $`𝔼[\rho ]=\overline{\rho }`$ the density contrast is defined by $`\delta (𝐱)=(\rho (𝐱)\overline{\rho })/\overline{\rho }`$. A homogeneous and isotropic Gaussian random field is stochastically fully specified by its mean $`\overline{\rho }`$ and its correlation function $`\xi _2^\delta (|𝐱𝐲|)=𝔼[\delta (𝐱)\delta (𝐲)]`$ . Here, $`𝔼`$ is the average over realizations of the random field. The higher (connected) correlation functions $`\xi _n^\delta =0`$ with $`n>2`$ vanish. Similar the correlation functions $`\xi _n`$ for $`n>2`$ vanish in a Gauss–Poisson process. However, also important differences between a Gaussian random field and a Gauss–Poisson point process show up.

### A Constraints on $`\xi _2(r)`$ and $`\varrho `$

A functional $`G[h]`$ defined by Eq. (13) is a p.g.fl. of a point process if and only if the $`P_k(𝐳)`$ as given in Eq. (4) are probability generating functions (p.g.f.’s). This will lead to restrictions on the two–point correlation function $`\xi _2`$ and the number density $`\varrho `$ as discussed by  and . A $`P_k(𝐳)`$ given by Eq. (4) always has to be positive and monotonic increasing with each component $`z_i`$ of $`𝐳`$, and hence $`\mathrm{log}P_k(𝐳)`$ is non–decreasing in each component of $`𝐳`$. With Eqs. (4), (5) and (13) one gets

$$\begin{array}{c}\frac{\mathrm{log}P_k(𝐳)}{z_l}=\varrho |A_l|+\hfill \\ \hfill +\varrho ^2\underset{j=1}{\overset{k}{}}_{A_l}d𝐱_{A_j}d𝐲\xi _2(|𝐱𝐲|)(z_j1)0\end{array}$$
(14)
for any $`z_j0`$, where $`|A_l|`$ is the volume of the set $`A_l`$. The rather obvious constraint $`\varrho 0`$ can be derived by setting $`z_j=1`$. For $`z_j=1,ji`$, and either $`z_i=0`$ or $`z_i1`$ the following two non–trivial constraints emerge:

$$\frac{\varrho }{|A_l|}_{A_l}d𝐱_{A_i}d𝐲\xi _2(|𝐱𝐲|)1$$
(15)
$$_{A_l}d𝐱_{A_i}d𝐲\xi _2(|𝐱𝐲|)\frac{|A_l|}{\varrho (z_i1)}\stackrel{z_i\mathrm{}}{}0.$$
(16)
One can show that these two conditions provide a necessary and sufficient characterization of $`\xi _2(r)`$ and $`\varrho `$, to assure that $`G[h]`$, as given in Eq. (13), is a p.g.fl. .

Eq. (15) constrains the shape and normalization of the two–point correlation functions admissible in a Gauss–Poisson process. For $`A_l=A=A_i`$

$$\frac{\varrho }{|A|}_Ad𝐱_Ad𝐲\xi _2(|𝐱𝐲|)=\overline{N}\sigma ^2(A)1,$$
(17)
where $`\sigma ^2(A)`$ are the fluctuations of count–in–cells in excess of a Poisson process, and $`\overline{N}=\varrho |A|`$ is the mean number of points inside the cell $`A`$. Hence, the total fluctuations of the number of points $`N`$ inside $`A`$ for a Gauss–Poisson process are

$$𝔼[(N\overline{N})^2]=\overline{N}+\overline{N}^2\sigma ^2(A)2\overline{N},$$
(18)
and must not exceed twice the value of the fluctuations in a Poisson process ($`\sigma ^2=0`$) for any domain $`A`$ considered. Another way of looking at constraint (15) is by taking $`A_l`$ as an infinitesimal volume element centered on the origin and $`A_j`$ equal to some volume $`A`$:

$$\varrho _Ad𝐲\xi _2(|𝐲|)1.$$
(19)
Consistent with Sect. III B this tells us that sitting on a point of the process on average at most one other point in excess of Poisson distributed points can be present.

Now consider two volume elements $`A_i=\mathrm{d}V_i`$ and $`A_l=\mathrm{d}V_l`$ separated by a distance of $`r`$, then Eq. (16) implies

$$\xi _2(r)0$$
(20)
Hence, only clustered point distributions can be modeled by a Gauss–Poisson process. Any zero crossing in $`\xi _2(r)`$ already indicates the presence of higher–order correlations.

### B A Gauss–Poisson process as a Poisson cluster process

A Gauss–Poisson process can be interpreted as a simple Poisson cluster process. This is important for simulations (see Appendix B 1).

A Poisson cluster process is a two–stage point process. First one chooses Poisson distributed cluster centers, the “parents”, with number density $`\varrho _p`$ and then attaches a second point process – the cluster to each cluster center (the cluster center is not necessarily part of the point process). The p.g.fl. of a Poisson cluster process is then given by

$$\mathrm{log}G[h]=_^dd𝐱\varrho _p\left(G_c[h|𝐱]1\right),$$
(21)
with $`G_c[h|𝐱]`$ being the p.g.fl. of the point process forming the cluster at center $`𝐱`$. Now consider the p.g.fl. of a cluster with at maximum two points (compare with Eq. (3)),

$$G_c[h|𝐱]=q_1(𝐱)h(𝐱)+q_2(𝐱)h(𝐱)_^dd𝐲f(|𝐱𝐲|)h(𝐲),$$
(22)
where $`q_1(𝐱)`$ is the probability that only one point, the cluster center at $`𝐱`$, is entering the cluster, whereas $`q_2(𝐱)`$ is the probability that a second point is added. Clearly, $`q_1(𝐱)+q_2(𝐱)=1`$. The probability density $`f(|𝐱𝐲|)`$ determines the distribution of the distance $`|𝐱𝐲|`$ of the second point $`𝐲`$ to the cluster center, and is normalized according to $`d𝐳f(|𝐳|)=1`$. By writing $`f(|𝐱𝐲|)`$ one assumes that the probability density $`f`$ is symmetric in $`𝐱`$ and $`𝐲`$. Indeed, the p.g.fl. Eq. (21) is invariant under interchanging $`𝐱`$ and $`𝐲`$, and this assumption does not impose any restrictions. Using this expression and Eq. (21) one obtains

$$\begin{array}{c}\mathrm{log}G[h+1]=_^dd𝐱\varrho _p(1+q_2(𝐱))h(𝐱)+\hfill \\ \hfill +_^dd𝐱_^dd𝐲\varrho _pq_2(𝐱)f(|𝐱𝐲|)h(𝐱)h(𝐲),\end{array}$$
(23)
which equals the p.g.fl. for the Gauss–Poisson process (13) for $`\varrho =\varrho _p(1+q_2)`$ and $`\xi _2(r)=2\frac{\varrho \varrho _p}{\varrho ^2}f(r)`$. Hence every Gauss–Poisson process is a Poisson cluster process of the above type, and vice versa.

### C Physical implications

From the preceding section one concludes that at maximum two points form a cluster in a Gauss–Poisson process. Therefore, no point distribution with large–scale structures can be modeled reliably with this kind of process. This has physical implications both for the galaxy distribution and percolating/critical systems.

More specific, from the observed galaxy distribution a scale–invariant two–point correlation function $`\xi _2(r)=Ar^\gamma `$ with $`\gamma 1.8`$ is deduced. Clearly such a correlation function does not satisfy the constraint (15). For $`0<\gamma 3`$ a cut–off at large scales has to be introduced. For the galaxy distribution a cut–off at approximately 20$`h^1`$Mpc is the lowest value which is roughly compatible with the observed two–point correlation function. Taking into account the observed number density of the galaxies, a cut–off even on this small scale does not help. Still the constraint (15) is strongly violated, indicating non negligible higher–order correlation functions (see also Sect. IV C). Similarly, a zero crossing or a negative $`\xi _2(r)`$ is violating the constraint (16) and also indicates that higher–order correlations are present. There are indications that the distribution of galaxies shows a negative $`\xi _2(r)`$ on some scale larger than 20$`h^1`$Mpc, followed by a positive peak at approximately 120$`h^1`$Mpc . A Gauss–Poisson process is not able to describe these features in the distribution of galaxies and galaxy clusters.

Also a percolating cluster shows scale–invariant correlations. The correlation length, specifying e.g. the exponential cut–off of the two–point correlation function, is going to infinity near the percolation threshold. Therefore, the geometry of the largest cluster cannot be modeled with a Gauss–Poisson processes. Higher–order correlations are essential to describe the morphology of such a system. This again illustrates that the tails of the distributions, in this case the asymptotic behavior of the two–point correlation function is essential.

To summarize these results: already by looking at the two–point correlation function and the density one is able to exclude a Gauss–Poisson process as a model. However one cannot turn the argument around and show that a given point distribution is compatible with a Gauss–Poisson process using the two–point correlation function alone. There are point processes with higher–order correlations satisfying the constraints (15,16) as discussed in Sect. IV A.

## IV Detecting deviations from a Gauss–Poisson process

After having outlined the basic theory of a Gauss–Poisson process, we discuss in this section how one can detect non–Gaussian features in a given point set.

### A The line–segment process

First a two dimensional analytic example is studied. In the line–segment process points are randomly distributed on line segments which are themselves uniformly distributed in space and direction. The number of points per line segment is a Poisson random variable. According to , p. 286

$$\xi _2(r)=\{\begin{array}{cc}\frac{1}{\pi l\varrho _s}\left(\frac{1}{r}\frac{1}{l}\right)\hfill & \text{ for }r<l\hfill \\ 0\hfill & \text{ for }rl;\hfill \end{array}$$
(24)
$`l`$ is the length of the line segments and $`\varrho _s`$ is the mean number density of line segments; $`l\varrho _s`$, $`\varrho /\varrho _s`$, $`\varrho `$ denote the mean length density, the mean number of points per line segment (which can be smaller than one), and the mean number density in space, respectively. A similar model for the distribution of galaxies was discussed by . On small scales $`rl`$, $`\xi _2(r)r^1`$, qualitatively similar to the observed two–point correlation function in the galaxy distribution.

This structured point process incorporates higher–order correlations. In Fig. 1 the line–segment process is shown in comparison to a Gauss–Poisson process with the same two–point correlation function for the parameters $`l=0.1`$, $`\varrho _s=201`$, $`\varrho =200`$, and $`\varrho =500`$. A number density $`\varrho >\varrho _s`$ violates the constraint Eq. (15) and no Gauss–Poisson process equivalent on the two–point level to such a line–segment process exists.

### B Detecting higher–order correlations

As can be seen from Fig. 1, the point processes are indistinguishable on the two–point level. For another example see . The differences between these point distributions can be investigated with statistical methods sensitive to higher–order correlations. One may use Minkowski functionals (, for reviews see ), percolation techniques , the minimum spanning tree , a method sensitive to three–point correlations , or directly calculate the higher moments . In the following the $`J`$–function is used to quantify the higher–order clustering .

To define $`J`$-function the spherical contact distribution $`F(r)`$ is needed, i.e. the distribution function of the distances $`r`$ between an arbitrary point and the nearest object in the point set. $`F(r)`$ is equal to one minus the void–probability function: $`F(r)=1P_0(r)`$. Another ingredient is the nearest neighbor distance distribution $`G(r)`$, defined as the distribution function of distances $`r`$ of an object in the point set to the nearest other point . For a Poisson process the probability to find a point only depends on the mean number density $`\varrho `$, leading to the well–known result

$$G_\mathrm{P}(r)=1\mathrm{exp}\left(\varrho |B_r|\right)=F_\mathrm{P}(r),$$
(25)
where $`|B_r|`$ is the volume of a $`d`$–dimensional sphere with radius $`r`$. The ratio

$$J(r)=\frac{1G(r)}{1F(r)}$$
(26)
was suggested by  as a probe for clustering of a point distribution. For a Poisson distribution $`J(r)=1`$ follows directly from Eq. (25). A clustered point distribution implies $`J(r)1`$, whereas regular structures are indicated by $`J(r)1`$. As discussed in  one can express the $`J(r)`$ function in terms of the $`n`$–point correlation functions $`\xi _n`$:

$$J(r)=1+\underset{l=1}{\overset{\mathrm{}}{}}\frac{(\varrho )^l}{l!}_{B_r}d𝐱_1\mathrm{}_{B_r}d𝐱_l\xi _{l+1}(0,𝐱_1,\mathrm{},𝐱_l).$$
(27)
$`B_r`$ is a $`d`$–dimensional sphere with radius $`r`$ centered on the origin. For a Gauss–Poisson process in two dimensions, i.e. $`\xi _n=0`$ for $`n>2`$, the above expression simplifiesUnfortunately  discussed this Gaussian approximation with examples of two–point correlation functions, which are not admissible in a Gauss–Poisson process.:

$$J(r)=1\varrho 2\pi _0^rdss\xi _2(s).$$
(28)
In Fig. 2 the results for $`F(r)`$, $`G(r)`$, and $`J(r)`$, estimated from several line–segment processes, and the Gauss–Poisson process are shown; all the processes investigated had the same two–point correlation function $`\xi _2(r)`$ given in Eq. (24). The line–segment process allows for larger voids than the Gauss–Poisson process, as seen from $`F_{\mathrm{line}}>F_{\mathrm{GP}}`$. On small scales the $`J(r)`$ of the line–segment process is well approximated by the $`J(r)`$ for the Gauss–Poisson process. However on large scales the Gauss–Poisson process shows significantly smaller $`J(r)`$ function than the line–segment process. The $`J(r)`$ function is known analytically for several point process models . In any of these cases a smaller $`J(r)`$ is an indication for stronger (positive) interaction between the points (see also ). Specifically for Gibbs–processes (see e.g. ) an attractive interaction leads to a monotonically decreasing $`J(r)`$ and a stronger interaction leads to smaller values of $`J(r)`$. Hence, the presence of higher–order correlation functions in the line–segment process gives rise to a reduced clustering strength, in the sense discussed above. Clearly, the signal of $`J(r)`$ also depends on the number density.

### C The non–Gaussian galaxy distribution

As already mentioned, the three–dimensional distribution of galaxies cannot be modeled in terms of a Gauss–Poisson process: the constraints on the density and two–point correlation function are violated. In the following this is illustrated with a volume–limited sample of 100$`h^1`$Mpc depth, extracted from the PSCz galaxy catalogue . The volume–limited sample incorporates 2232 galaxies with galactic latitude $`|b|>5^{}`$. A detailed description of the sample considered here may be found in . Estimators for the two–point correlation function are quite abundant (see  and references therein). The results presented here do neither depend on the estimator, nor on the exact sample geometry, which is indeed more complicated (see ). For the $`J(r)`$–function the minus estimator is used .

In Fig. 3 the estimated two–point correlation function is shown. The integral

$$\varrho _^3d𝐲\xi _2(|𝐲|)4.4$$
(29)
is violating the constraint (19), and the corresponding Gauss-Poisson process does not exist. Indeed higher–order correlations functions have been detected by  using factorial moments. By thinning the galaxy distribution (i.e. randomly sub–sampling), one generates a point set with the same correlation functions $`\xi _n`$ as the observed galaxy distribution, however with a reduced number of points. Since the number density enters linearly in the constraint (19), a comparison of the thinned galaxy distribution with a Gauss–Poisson process becomes feasible. The strongly interacting galaxy distribution, as indicated by the small values of $`J(r)`$, shows increasingly weaker interaction (higher values of $`J(r)`$) for the diluted subsamples (Fig. 3).

Now consider a sample with only 20% of the actual observed galaxies, where the constraint (19) is satisfied (compare with (29)). This diluted sample is compared to a Gauss–Poisson process with the same two–point correlation function. The $`\xi _2(r)`$ determined from the simulated Gauss–Poisson process matches perfectly with the observed correlation function (Fig. 3). On small scales the $`J(r)`$-function of the thinned PSCz is reasonably modeled by the Gauss–Poisson process. However, on large scales the Gauss–Poisson process shows stronger interactions, whereas the thinned galaxy sample, with its higher–order correlation functions, shows weaker interactions in the sense discussed in Sect. IV B.

## V Point processes with higher–order clustering

As already mentioned, the measured two–point correlation function of the galaxy distribution together with the observed density of galaxies violates the constraints Eqs. (15) and (16). Consequently the distribution of galaxies cannot be modeled with a Gauss–Poisson process. Even more compelling, there is a clear detection of higher–order correlations in the galaxy distribution (e.g. ). Hence, one is interested in analytical tractable approximations of the cumulant expansion (7). Hierarchical closure relations have been extensively studied (see Sect. VII A). In the following a truncation of the expansion (7) beyond the Gaussian term and the $`n`$–point Poisson cluster processes will be used.

Such a truncation may serve as a model for the galaxy distribution in the weakly nonlinear regime. Using perturbation theory one may show that $`\xi _3(\xi _2)^2`$ (see ). For large separations $`r`$ the correlation function $`\xi _2(r)`$ is smaller than unity, and consequently a truncation of the expansion (7) at $`n>2`$ provides a viable model for the large–scale distribution of galaxies.

The general Poisson cluster process is the starting point: consider the expansion of the cluster p.g.fl. $`G_c[h|𝐱]`$ in terms of Janossy densities conditional on the cluster center $`𝐱`$ (see Eq. (11)):

$$\begin{array}{c}G_c[h|𝐱]=J_0+\underset{n=1}{\overset{\mathrm{}}{}}\frac{1}{n!}_^dd𝐱_1\mathrm{}_^dd𝐱_n\hfill \\ \hfill j_n(𝐱_1,\mathrm{}𝐱_n|𝐱)h(𝐱_1)\mathrm{}h(𝐱_n).\end{array}$$
(30)
Explicit expression for the Janossy densities are given below. The p.g.fl. of a Poisson cluster process is then given by

$`G[h]`$ $`=\mathrm{exp}\left({\displaystyle _^d}d𝐱\varrho _c(G_c[h|𝐱]1)\right)`$ (31)
$`=\mathrm{exp}({\displaystyle _^d}\mathrm{d}𝐱\varrho _c({\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{n!}}{\displaystyle _^d}\mathrm{d}𝐱_1\mathrm{}{\displaystyle _^d}\mathrm{d}𝐱_n`$ (32)
$`j_n(𝐱_1,\mathrm{}𝐱_n|𝐱)h(𝐱_1)\mathrm{}h(𝐱_n)1)).`$ (33)
Here the probability $`q_0`$ of having no point in the cluster at $`𝐱`$ is assumed to be zero, i.e. $`J_0=0`$. This does not impose any additional constraints, it only leads to a redefinition of the number density of cluster centers $`\varrho _c^{}=\varrho _c(1+q_0)`$.

Using this more formal approach the p.g.fl. of the Gauss–Poisson process can be written in terms of the Janossy densities with $`j_n=0`$ for $`n>2`$:

$$\begin{array}{cc}\hfill j_1(𝐱_1|𝐱)& =q_1(𝐱)\delta ^D(𝐱𝐱_1),\hfill \\ \hfill j_2(𝐱_1,𝐱_2|𝐱)& =q_2(𝐱)\delta ^D(𝐱𝐱_1)f_2(𝐱_2|𝐱_1)2!.\hfill \end{array}$$
(34)
Here $`j_1`$ ($`j_2`$) are the probability densities for the spatial distribution of one (two) points in the cluster, multiplied by the probability $`q_1`$ ($`q_2`$) that there are exactly one (two) points in the cluster at $`𝐱`$. $`\delta ^D`$ is the $`d`$–dimensional Dirac distribution. $`f_2(𝐱_2|𝐱_1)`$ is the probability density of the second point $`𝐱_2`$ under the condition that there is a point at $`𝐱_1`$, normalized by $`d𝐱_2f_2(𝐱_2|𝐱_1)=1`$.

The p.g.fl. (31) is invariant under changes of the order of integration, implying that one can use the $`j_n(𝐱_1,\mathrm{}𝐱_n|𝐱)`$ symmetrically defined in all coordinates (including $`𝐱`$). With the additional assumption of homogeneity and isotropy one gets $`f_2(𝐱_2|𝐱_1)=f(|𝐱_1𝐱_2|)`$, as already used in Sect. III B for the construction of the Gauss–Poisson process.

### A The three–point Poisson cluster process

In a three–point Poisson cluster process Eq. (31) is truncated at the third order and at most three points per cluster are allowed. Additional to Eq. (34)

$$j_3(𝐱_1,𝐱_2,𝐱_3|𝐱)=q_3(𝐱)\delta ^D(𝐱𝐱_1)f_3(𝐱_2,𝐱_3|𝐱_1)3!$$
(35)
appears, with the probability $`q_3(𝐱)`$ that the cluster consists out of three points, and $`q_1+q_2+q_3=1`$ with $`q_i0`$. $`f_3(𝐱_2,𝐱_3|𝐱_1)`$ is the probability density that there are two points at $`𝐱_2`$, and $`𝐱_3`$, under the condition that one point is at $`𝐱_1`$, with the normalization $`d𝐱_2d𝐱_3f_3(𝐱_2,𝐱_3|𝐱_1)=1`$. Inserting these definitions one obtains

$$\begin{array}{c}\mathrm{log}G[h]=_^dd𝐱_1\varrho _pq_1(𝐱_1)(h(𝐱_1)1)+\hfill \\ \hfill +_^dd𝐱_1\varrho _p_^dd𝐱_2q_2(𝐱_1)f_2(𝐱_2|𝐱_1)(h(𝐱_1)h(𝐱_2)1)+\\ \hfill +_^d\mathrm{d}𝐱_1\varrho _p_^d\mathrm{d}𝐱_2_^d\mathrm{d}𝐱_3q_3(𝐱_1)f_3(𝐱_2,𝐱_3|𝐱_1)\times \\ \hfill \times \left(h(𝐱_1)h(𝐱_2)h(𝐱_3)1\right).\end{array}$$
(36)
As already mentioned, $`f_3(𝐱_2,𝐱_3|𝐱_1)`$ can be assumed to be symmetric in its three arguments. Slightly abusing notation, let $`f_2(𝐱_1,𝐱_2)`$ and $`f_3(𝐱_1,𝐱_2,𝐱_3)`$ be the symmetrically defined densities corresponding to $`f_2(𝐱_2|𝐱_1)`$ and $`f_3(𝐱_2,𝐱_3|𝐱_1)`$, and define

$$f_2^{(3)}(𝐱_1,𝐱_2)=_^dd𝐱_3f_3(𝐱_1,𝐱_2,𝐱_3).$$
(37)
Replacing $`h`$ by $`h+1`$ and rearranging the terms the factorial cumulant expansion of the three–point cluster process reads

$$\begin{array}{c}\mathrm{log}G[h+1]=_^dd𝐱_1h(𝐱_1)\varrho _p\left(1+q_2(𝐱_1)+2q_3(𝐱_1)\right)+\hfill \\ \hfill \begin{array}{cc}\hfill +& _^d\mathrm{d}𝐱_1_^d\mathrm{d}𝐱_2h(𝐱_1)h(𝐱_2)\times \hfill \\ & \times \varrho _p\left(q_2(𝐱_1)f_2(𝐱_1,𝐱_2)+3q_3(𝐱_1)f_2^{(3)}(𝐱_1,𝐱_2)\right)+\hfill \\ \hfill +& _^d\mathrm{d}𝐱_1_^d\mathrm{d}𝐱_2_^d\mathrm{d}𝐱_3h(𝐱_1)h(𝐱_2)h(𝐱_3)\times \hfill \\ & \times \varrho _pq_3(𝐱_1)f_3(𝐱_1,𝐱_2,𝐱_3).\hfill \end{array}\end{array}$$
(38)
Comparing Eq. (38) with the expansion (7) one arrives at

$`\varrho `$ $`=(q_1+2q_2+3q_3)\varrho _p=(1+q_2+2q_3)\varrho _p`$ (39)
$`\xi _2(𝐱_1,𝐱_2)`$ $`={\displaystyle \frac{2!}{(1+q_2+2q_3)^2\varrho _p}}\times `$ (40)
$`\times \left(q_2f_2(𝐱_1,𝐱_2)+3q_3f_2^{(3)}(𝐱_1,𝐱_2)\right)`$ (41)
$`\xi _3(𝐱_1,𝐱_2,𝐱_3)`$ $`={\displaystyle \frac{3!}{(1+q_2+2q_3)^3\varrho _p^2}}q_3f_3(𝐱_1,𝐱_2,𝐱_3),`$ (42)
and the correlation functions $`\xi _n`$ equal zero for $`n4`$. The simulation procedure for the three–point Poisson cluster process is described in Appendix B 2.

### B Constraints on $`\varrho `$, $`\xi _2`$ and $`\xi _3`$

By the definition of the three–point Poisson cluster process, the probability densities $`f_20`$, $`f_30`$, and $`f_2^{(3)}0`$ and consequently $`\xi _2(r)0`$ for all $`r`$, as well as $`\xi _30`$. This is a generic feature of Poisson cluster processes.

The Gauss–Poisson process, defined through the truncation of the cumulant expansion after the second term, is equivalent to the two–point Poisson cluster process (see Sect. III B). Unfortunately, this equivalence does not hold for the higher $`n`$–point processes anymore. The general three–point process is defined as point process with a factorial cumulant expansion truncated after the third term. Proceeding similar to Sect. III A, necessary conditions for the existence of such a point process can be derived (compare with Eq. (14)):

$$\begin{array}{c}0\varrho |A_l|+\underset{i=1}{\overset{k}{}}\varrho ^2|A_l||A_i|\overline{\xi }_2(A_l,A_i)(z_i1)+\hfill \\ \hfill +\underset{i=1}{\overset{k}{}}\underset{j=1}{\overset{k}{}}\frac{\varrho ^3}{2}|A_l||A_i||A_j|\overline{\xi }_3(A_l,A_i,A_j)(z_i1)(z_j1),\end{array}$$
(43)
with the volume–averaged correlation functions

$$\begin{array}{c}\overline{\xi }_n(A_1,\mathrm{},A_n)=\frac{1}{|A_1|\mathrm{}|A_n|}\times \hfill \\ \hfill \times _{A_1}\mathrm{d}𝐱_1\mathrm{}_{A_n}\mathrm{d}𝐱_n\xi _n(𝐱_1,\mathrm{},𝐱_n),\end{array}$$
(44)
and for consistency $`\overline{\xi }_1(A)=1`$. Again, for $`z_i=1`$ one obtains $`\varrho 0`$. The non–trivial constraints read:

$`1`$ $`\varrho |A_s|\overline{\xi }_2(A_l,A_s)\varrho ^2|A_s|^2\overline{\xi }_3(A_l,A_s,A_s),`$ (45)
$`0`$ $`\overline{\xi }_3(A_l,A_s,A_s),`$ (46)
$`1`$ $`\varrho |A_s|\overline{\xi }_2(A_l,A_s)+\varrho |A_r|\overline{\xi }_2(A_l,A_r)`$ (47)
$`{\displaystyle \frac{\varrho ^2}{2}}|A_s|^2\overline{\xi }_3(A_l,A_s,A_s){\displaystyle \frac{\varrho ^2}{2}}|A_r|^2\overline{\xi }_3(A_l,A_r,A_r)`$ (48)
$`\varrho ^2|A_s||A_r|\overline{\xi }_3(A_l,A_s,A_r).`$ (49)
Eq. (45) can be derived by setting $`z_s=0`$ and $`z_i=1`$ for all $`is`$, Eq. (46) follows from $`z_s\mathrm{}`$ and $`z_i=1`$ for all $`is`$. Using $`z_r,z_s\mathrm{}`$ does not lead to new constraints. With $`z_r=0=z_s`$, $`rs`$, and $`z_i=1`$ for all $`ir,s`$ one obtains Eq. (47). No additional constraint arises by setting $`z=2`$.

Eq. (46) implies $`\xi _30`$. Eq. (45) and Eq. (47) are the extension of the constraint (17). The terms proportional to $`\overline{\xi }_3`$ can balance the terms with $`\overline{\xi }_2`$, and a clustering point processes with a number density higher than in a Gauss–Poisson process is possible. Moreover, $`\xi _2`$ is not constrained to positive values anymore. Hence, already by including three–point correlations, a point process model with a two–point correlation function $`\xi _2`$ having a zero crossing becomes admissible. This answers affirmatively the question by , whether there exists a general three–point cluster processes with a negative second moment. However, in the three–point Poisson cluster process discussed in the preceding section a $`\xi _20`$ is required illustrating that the three–point Poisson cluster processes form only a subset of all possible three–point processes.

### C The $`n`$–point Poisson cluster process

It is now clear how to construct the $`n`$–point Poisson cluster process. Let $`q_m`$ be the probability of having $`m`$ points per cluster with $`_{m=1}^nq_m=1`$.

$$j_n(𝐱_1,\mathrm{},𝐱_n|𝐱)=q_n(𝐱)\delta ^D(𝐱𝐱_1)f_n(𝐱_2,\mathrm{},𝐱_n|𝐱_1)n!,$$
(50)
determines the distribution of the $`n`$ points inside the cluster ($`f_1=1`$). As above the $`f_n`$ are assumed to be symmetric in all their arguments, and for $`n>m`$

$$\begin{array}{c}f_m^{(n)}(𝐱_1,\mathrm{},𝐱_m)=_^dd𝐱_{m+1}\mathrm{}_^dd𝐱_n\hfill \\ \hfill f_n(𝐱_1,\mathrm{},𝐱_m,𝐱_{m+1},\mathrm{},𝐱_n),\end{array}$$
(51)
and $`f_m=f_m^{(m)}`$. Inserting Eq. (50) into Eq. (31) and after some algebraic manipulations one can compare term by term with the factorial cumulant expansion (7) of the p.g.fl.:

$`\varrho `$ $`=\varrho _p{\displaystyle \underset{m=1}{\overset{n}{}}}mq_m`$ (52)
$`\xi _k(𝐱_1,\mathrm{},𝐱_k)`$ $`={\displaystyle \frac{k!\varrho _p}{\varrho ^k}}{\displaystyle \underset{m=k}{\overset{n}{}}}\left({\displaystyle \genfrac{}{}{0pt}{}{m}{k}}\right)q_mf_k^{(m)}(𝐱_1,\mathrm{},𝐱_k),`$ (53)
with $`kn`$, and $`\xi _k=0`$ for $`k>n`$. The statistical properties of this $`n`$–point Poisson cluster process are now completely specified by the correlation functions $`\xi _k`$ with $`kn`$ and the mean density $`\varrho `$. Eqs. (52) and the normalization of the $`f_m`$ can be used to determine the $`f_m`$ as well as $`\varrho _p`$ and $`q_m`$ from given correlation functions $`\xi _n`$ and the number density $`\varrho `$. A simulation algorithm similar to the one described in Appendix B 2 can be constructed.

### D The general $`n`$–point process

The general $`n`$–point process is defined as the point process resulting from the factorial cumulant expansion truncated after the $`n`$th term. Proceeding similarly to Sect. III A one arrives at the constraint equations

$$\begin{array}{c}0\varrho |A_l|+\varrho ^2\underset{i=1}{\overset{k}{}}|A_l||A_i|\overline{\xi }_2(A_l,A_i)(z_i1)+\mathrm{}+\hfill \\ \hfill +\frac{\varrho ^nn}{n!}\underset{i_1,\mathrm{},i_{n1}=1}{\overset{k}{}}|A_l||A_{i_1}|\mathrm{}|A_{i_{n1}}|\times \\ \hfill \times \overline{\xi }_n(A_l,A_{i_1},\mathrm{},A_{i_{n1}})(z_{i_1}1)\mathrm{}(z_{i_{n1}}1).\end{array}$$
(54)
It is now possible to compute the constraints for the $`n`$–point process, in close analogy to the three–point process in Sect. V B.  gave necessary and sufficient conditions for the existence of a generalized Hermite distribution (closely related to this $`n`$–point process). They discuss the constraints for a slightly different expansion of the p.g.f. Unfortunately, the transformation of their expansion to the expansion in terms of correlation functions is as tedious as the direct calculation of the constraints.

## VI Random fields vs. point processes

A random field $`u(𝐱)`$ is in the simplest case a real–valued function on $`^d`$ . In cosmology the initial mass–density field is often modeled as a Gaussian random field (see e.g. ). The nonlinear evolution of the density field unavoidably introduces higher–order correlations. A random field $`u(x)`$ is stochastically characterized by its characteristic functional (e.g. )

$$\mathrm{\Phi }^u[v]=𝔼^u\left[\mathrm{exp}\left(i_^dd𝐱v(𝐱)u(𝐱)\right)\right],$$
(55)
where $`𝔼^u`$ denotes the expectation value over realizations of the random field $`u`$. In close analogy to the expansion (A5) of the characteristic function of a random variable in terms of cumulants, one obtains the expansion of the characteristic functional

$$\begin{array}{c}\mathrm{ln}\mathrm{\Phi }^u[v]=\underset{n=1}{\overset{\mathrm{}}{}}\frac{i^n}{n!}_^dd𝐱_1\mathrm{}_^dd𝐱_n\hfill \\ \hfill c_n^u(𝐱_1,\mathrm{},𝐱_n)v(𝐱_1)\mathrm{}v(𝐱_n)\end{array}$$
(56)
in terms of $`n`$-point cumulants $`c_n^u(𝐱_1,\mathrm{},𝐱_n)`$. Here $`c_1^u(𝐱)=𝔼^u[u(𝐱)]=\overline{u}`$ is the mean value. The correlation function of the field is

$$\xi _2^u(𝐱_1,𝐱_2)=\frac{c_2^u(𝐱_1,𝐱_2)}{\overline{u}^2}=\frac{𝔼^u[u(𝐱_1)u(𝐱_2)]}{\overline{u}^2}1,$$
(57)
and similar for higher–order correlation functions (see e.g. ). The well known characteristic functional of the Gaussian random field reads

$$\begin{array}{c}\mathrm{ln}\mathrm{\Phi }^u[v]=i_^dd𝐱_1\overline{u}v(𝐱_1)\hfill \\ \hfill \frac{1}{2}_^dd𝐱_1_^dd𝐱_2c_2^u(𝐱_1,𝐱_2)v(𝐱_1)v(𝐱_2),\end{array}$$
(58)
with the covariance function $`c_2^u(𝐱_1,𝐱_2)`$.

The characteristic functional of a point process is defined by

$$\mathrm{\Phi }[h]=𝔼\left[\mathrm{exp}\left(i_^dN(\mathrm{d}𝐱)h(𝐱)\right)\right],$$
(59)
and the relation to the p.g.fl. is $`\mathrm{\Phi }[h]=G[\mathrm{e}^{ih}]`$. An expansion into cumulants $`c_k(𝐱_1,\mathrm{},𝐱_k)`$ is also possible:

$$\begin{array}{c}\mathrm{log}\mathrm{\Phi }[h]=\underset{k=1}{\overset{\mathrm{}}{}}\frac{i^k}{k!}_^dd𝐱_1\mathrm{}_^dd𝐱_k\hfill \\ \hfill c_k(𝐱_1,\mathrm{},𝐱_k)h(𝐱_1)\mathrm{}h(𝐱_k).\end{array}$$
(60)
The cumulants $`c_k()`$ should not be confused with the factorial cumulants $`c_{[k]}()`$.

### A A theorem of Marcinkiewicz

A theorem of Marcinkiewicz  states that if the characteristic function $`\phi (s)`$ of a random variable (see Appendix A) is the exponential of a polynomial with finite degree larger than two, then the positive definiteness of the probability distribution is violated (see e.g. ). The generalized Marcinkiewicz theorem for characteristic functionals  tells us that this expansion has to be either infinite or a polynomial in $`h(𝐱)`$ (or $`v(𝐱)`$) of degree less than or equal to two. This directly applies to the expansion of the characteristic functionals of a random field $`\mathrm{\Phi }^u[v]`$ and a point process $`\mathrm{\Phi }[h]`$ in terms of cumulants.

However, for a point process one can see that the expansion of the p.g.fl. in terms of factorial cumulants $`c_{[k]}()`$ (or correlation functions $`\xi _k()`$) allows a truncation at a finite $`k>2`$. As long as constraints on the density and the correlation functions are fulfilled, the point process is well defined. Although the p.g.fl. was used mainly in the context of discrete events, it seems worthwhile to consider the characterizations of random fields with factorial cumulants.

Another systematic expansion is provided by the Edgeworth series. It was successfully applied in cosmology to quantify the one point probability distribution function for the smoothed density field on large scales . Recently,  showed how to use the truncated Edgeworth series to generate realizations of non–Gaussian random fields with predefined correlation properties. The truncated Edgeworth series also violates the positive definiteness of the probability distribution, but  restore the positivity, reintroducing higher correlations, leading to a “leaking” into higher correlations.

The cumulants $`c_k`$ and the factorial cumulants $`c_{[k]}`$ of a random variable are related by $`c_k=_{l=1}^ks_2(k,l)c_{[l]}`$ (see Eq. (A16)). Looking at the Poisson cluster processes discussed in the preceding sections, one observes that such a relation must not hold between the cumulants and the factorial cumulants of a point process. As an example consider the three–point Poisson cluster process with $`c_{[n]}()=0`$ for all $`n>3`$. A finite $`c_{[3]}()`$ leads to non–zero $`c_n()`$ for all $`n`$ (see Appendix C for details.)

### B The Poisson model

In cosmology the point distribution is often related to the mass density field assuming the “Poisson model”. The value of the mass density field is assumed to be proportional to the local number density, and the point distribution is constructed by “Poisson sampling” the correlated mass density field. If the mass density field is itself a realization of a random field, the resulting point process is called a Cox process, or doubly stochastic process. Within this model one may show that the cumulants $`c_n^u`$ of the density field are proportional to the factorial cumulants $`c_{[n]}`$ of the point distribution : $`c_{[n]}=c_n^u\rho ^2/\overline{u}^2`$. It is important to notice that this relates the characteristic functional of the random field $`\mathrm{\Phi }^u[v]`$ with the p.g.fl of the point process $`G[h]`$. For the correlation functions one obtains

$$\xi _n(𝐱_1,\mathrm{},𝐱_n)=\xi _n^u(𝐱_1,\mathrm{},𝐱_n).$$
(61)
Hence, this model allows the direct comparison of predictions from analytical calculations with the observed correlation functions in the galaxy distribution.

Clearly the question arises, what is wrong with the simple picture that one starts with a Gaussian random field and “Poisson sample” it to obtain the desired point distribution. The answer is that a Gaussian random field is an approximate model for a mass density field only if the fluctuations are significantly smaller than the mean mass density. Otherwise negative mass densities (i.e. negative “probabilities” for the Poisson sampling) would occur. Only in the limit of vanishing fluctuations a Poisson sampled Gaussian random field becomes a permissible model. However, in this limit one is left with a pure Poisson process.

## VII Models for strongly correlated systems

In the Sects. III and V several types of point processes were discussed, all featuring a truncated factorial cumulant expansion. As argued at the beginning of Sect. V, such a truncation is feasible for the matter distribution in the Universe, as long as $`\xi _2(r)<1`$, i.e. for points with large separations. Mainly two approaches have been followed to model the galaxy distribution also on small scales with $`\xi _2(r)1`$. The hierarchical models are briefly discussed in the next section and in Sect. VII B an extension of the halo–model is presented.

### A Hierarchical models

In cosmology one often starts with a scale–invariant correlation function $`\xi _2(r)r^\gamma `$ and assumes some closure relations for the $`\xi _n`$. Especially the hierarchical ansatz $`\xi _n=Q_n_{\text{trees}}^{n1}\xi _2`$ was extensively studied (e.g. , and more recent ).  discuss conditions for the coefficients $`Q_n`$ such that the expansion of the p.g.f.’s in terms of the count–in–cells converges. In this case the count–in–cells uniquely determine the point process. As illustrated in Sect. VIII with the log–normal distribution, a non–converging expansion does not necessarily imply that the stochastic model is not well–defined. It only implies that such a point process model is not completely specified by its correlation functions. For critical systems similar expansion in terms of correlation functions are typically divergent (see e.g. , chapt. 41).

As another closure relation Kirkwood  employed the following approximation

$$\begin{array}{c}\varrho _3(𝐱_1,𝐱_2,𝐱_3)=\hfill \\ \hfill \varrho ^3(1+\xi _2(𝐱_1,𝐱_2))(1+\xi _2(𝐱_2,𝐱_3))(1+\xi _2(𝐱_3,𝐱_1))\end{array}$$
(62)
to calculate thermodynamic properties of fluids using the BBGKY hierarchy. This closure relation is exact for the log–normal distribution (e.g. ). Empirically however one finds that this ansatz is disfavored as a model for the galaxy distribution .

### B The generalized halo model

In Sect. V several types of Poisson cluster processes were constructed by starting with Poisson distributed centers and attaching a secondary point process, the cluster, to each point. One can generalize this procedure by considering cluster centers given by already correlated points. One possibility, is to iterate the construction principle of the simple Poisson cluster process leading to the $`m`$–th order Neyman–Scott processes . If one is only interested in the first few correlation functions, the full specification of the point process is not necessary. Within the halo model (see e.g. ) it is specifically easy to calculate the correlation functions. The difference to the Poisson cluster processes discussed previously is that the cluster centers now may be correlated themselves. The major physical assumption entering is that the properties of the clusters (halos) are independent from the positions and correlations of the cluster centers.

Consider a point process for the cluster centers, the parents, specified by the p.g.fl. $`G_p[h]`$. Independent from the distribution of the centers, a cluster with a p.g.fl. $`G_c[h|𝐲]`$ is attached to each center $`𝐲`$. Then the p.g.fl. of this cluster process is given by the “folding” of the two p.g.fl.’s :

$$G[h]=G_p\left[G_c[h|]\right].$$
(63)
Using the expansion (7), these p.g.fl.’s are given by

$`\mathrm{log}G_p[h+1]`$ $`={\displaystyle \underset{m=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{\varrho _p^m}{m!}}{\displaystyle _^d}d𝐲_1\mathrm{}{\displaystyle _^d}d𝐲_m`$ (64)
$`\xi _m^{(p)}(𝐲_1,\mathrm{},𝐲_m)h(𝐲_1)\mathrm{}h(𝐲_m),`$ (65)
$`\mathrm{log}G_c[h+1|𝐲]`$ $`={\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{n!}}{\displaystyle _^d}d𝐱_1\mathrm{}{\displaystyle _^d}d𝐱_n`$ (66)
$`c_{[n]}(𝐱_1,\mathrm{},𝐱_n|𝐲)h(𝐱_1)\mathrm{}h(𝐱_n),`$ (67)
where $`\varrho _p`$ is the number density and the $`\xi _n^{(p)}`$ are the correlation functions of the parent process. The $`c_{[n]}(\mathrm{}|𝐲)`$ are the factorial cumulants specifying the point distribution in a cluster, conditional on the cluster center $`𝐲`$. $`c_{[1]}(𝐱|𝐲)`$ is the halo profile, with the mean number of points per halo $`\mu =d𝐱c_{[1]}(𝐱|𝐲)`$. The $`c_{[n]}`$, $`n2`$ quantify the halo substructure. A halo without substructure is an inhomogeneous Poisson process, and completely characterized by $`c_{[1]}(𝐱|𝐲)`$ and $`c_{[n]}=0`$, $`n2`$. Combining Eqs. (63) and (64)

$$\begin{array}{c}\mathrm{log}G[h+1]=\underset{m=1}{\overset{\mathrm{}}{}}\frac{\varrho _p^m}{m!}_^dd𝐲_1\mathrm{}_^dd𝐲_m\hfill \\ \hfill \xi _m^{(p)}(𝐲_1,\mathrm{},𝐲_m)\underset{i=1}{\overset{m}{}}\left(G_c[h+1|𝐲_i]1\right),\end{array}$$
(68)
one immediately recovers the p.g.fl. of the Poisson cluster process Eq. (21) by setting $`\xi _m^{(p)}=0`$ for $`m2`$ ($`\xi _1^{(p)}=1`$).

### C $`\xi _2`$ and $`\xi _3`$ in the generalized halo model

In the standard halo model the clusters are simply modeled by an inhomogeneous Poisson process, whereas the centers are given by a correlated point process, typically determined from the evolved density distribution. Based on these assumptions one can calculate the correlation functions $`\xi _n`$ for the halo model . Both theoretical models as well as observations suggest that dark matter caustics lead to substructure inside halos . Also recent high–resolution $`N`$–body simulations suggest that 15%–40% of the simulated halos show a significant amount of substructure (see  and references therein). To generalize the halo model, the correlations inside the halo are taken into account.

Consider the expansion of $`G_c[h|𝐲_i]1`$ in $`h`$:

$`G_c[h|𝐲_i]1={\displaystyle _^d}d𝐱_1h(𝐱_1)c_{[1]}(𝐱_1|𝐲_i)+`$ (69)
$`+{\displaystyle \frac{1}{2!}}{\displaystyle _^d}\mathrm{d}𝐱_1{\displaystyle _^d}\mathrm{d}𝐱_2h(𝐱_1)h(𝐱_2)\times `$ (70)
$`\times \left(c_{[1]}(𝐱_1|𝐲_i)c_{[1]}(𝐱_2|𝐲_i)+c_{[2]}(𝐱_1,𝐱_2|𝐲_i)\right)+`$ (71)
$`+{\displaystyle \frac{1}{3!}}{\displaystyle _^d}\mathrm{d}𝐱_1{\displaystyle _^d}\mathrm{d}𝐱_2{\displaystyle _^d}\mathrm{d}𝐱_3h(𝐱_1)h(𝐱_2)h(𝐱_3)\times `$ (72)
$`\times (c_{[1]}(𝐱_1|𝐲_i)c_{[1]}(𝐱_2|𝐲_i)c_{[1]}(𝐱_3|𝐲_i)+`$ (73)
$`+3c_{[2]}(𝐱_1,𝐱_2|𝐲_i)c_{[1]}(𝐱_3|𝐲_i)+c_3(𝐱_1,𝐱_2,𝐱_3|𝐲_i))+`$ (74)
$`+O[h^4].`$ (75)
After inserting this expansion into Eq. (68) and collecting terms proportional to powers of $`\varrho h()`$, with the mean number density $`\varrho =\varrho _p\mu `$, one can directly compare with the expansion (7) and read off the correlation functions:

$$\begin{array}{c}\xi _2(𝐱_1,𝐱_2)=\hfill \\ \hfill =\frac{1}{\varrho _p\mu ^2}_^dd𝐲_1\left(c_{[1]}(𝐱_1|𝐲_1)c_{[1]}(𝐱_2|𝐲_1)+c_{[2]}(𝐱_1,𝐱_2|𝐲_1)\right)+\\ \hfill +\frac{1}{\mu ^2}_^dd𝐲_1_^dd𝐲_2\xi _2^{(p)}(𝐲_1,𝐲_2)c_{[1]}(𝐱_1|𝐲_1)c_{[1]}(𝐱_2|𝐲_2),\end{array}$$
(76)
$`\xi _3(𝐱_1,𝐱_2,𝐱_3)=`$ (77)
$`=`$ $`{\displaystyle \frac{1}{\varrho _p^2\mu ^3}}{\displaystyle _^d}\mathrm{d}𝐲_1(c_{[1]}(𝐱_1|𝐲_1)c_{[1]}(𝐱_2|𝐲_1)c_{[1]}(𝐱_3|𝐲_1)+`$ (78)
$`+3c_{[1]}(𝐱_1|𝐲_1)c_{[2]}(𝐱_2,𝐱_3|𝐲_1)+c_3(𝐱_1,𝐱_2,𝐱_3|𝐲_1))+`$ (79)
$`+{\displaystyle \frac{3}{\varrho _p\mu ^3}}{\displaystyle _^d}\mathrm{d}𝐲_1{\displaystyle _^d}\mathrm{d}𝐲_2\xi _2^{(p)}(𝐲_1,𝐲_2)\times `$ (80)
$`\times (c_{[1]}(𝐱_1|𝐲_1)c_{[1]}(𝐱_2|𝐲_2)c_{[1]}(𝐱_3|𝐲_2)+`$ (81)
$`+c_{[1]}(𝐱_1|𝐲_1)c_{[2]}(𝐱_2,𝐱_3|𝐲_2))+`$ (82)
$`+{\displaystyle \frac{1}{\mu ^3}}{\displaystyle _^d}\mathrm{d}𝐲_1{\displaystyle _^d}\mathrm{d}𝐲_2{\displaystyle _^d}\mathrm{d}𝐲_3\xi _3^{(p)}(𝐲_1,𝐲_2,𝐲_3)\times `$ (83)
$`\times c_{[1]}(𝐱_1|𝐲_1)c_{[1]}(𝐱_2|𝐲_2)c_{[1]}(𝐱_3|𝐲_3).`$ (84)
Similarly the higher $`n`$–point correlation functions can be calculated.

In current calculations of the two– and three–point functions for the halo model  the galaxies inside the halos are modeled as an inhomogeneous (finite) Poisson process. The halo profile $`c_{[1]}(𝐱|𝐲)`$ is conditional on the cluster center $`𝐲`$, but no substructure inside halos is present, i.e. $`c_{[n]}=0`$ for $`n2`$. In this case the above expressions simplify to the result of .

The simulation of such a point distribution can be carried out in a multi–step approach similar to the simulation of the Gauss–Poisson process (Appendix B 1). First generate the correlated cluster centers, e.g. by using a Gauss–Poisson process or a low–resolution simulation, and then attach a secondary point process either modeled as an inhomogeneous Poisson or $`n`$–point Poisson cluster process.

### D Halo substructure

The following discussion shall serve mainly as an illustration of how to incorporate halo substructure in calculations of the correlation function. To keep things simple the following assumptions are made: the halo profile $`c_{[1]}(𝐱|𝐲)=c_{[1]}(|𝐱𝐲|)`$ is independent from the mass of the halo, and $`c_{[2]}(𝐱_1,𝐱_2|𝐲)`$ factors into $`c_{[1]}(𝐱_1|𝐲)c_{[1]}(𝐱_2|𝐲)\gamma (|𝐱_1𝐱_2|)`$, as expected for locally isotropic substructures. Let $`P^{(p)}(k)=\frac{\mathrm{d}𝐱}{(2\pi )^3}\xi _2^{(p)}(|𝐱|)\mathrm{e}^{i𝐤𝐱}`$ be the power spectrum of the spatial distribution of the halo centers, and let $`\stackrel{~}{c}_{[1]}(k)`$ and $`\stackrel{~}{\gamma }(k)`$ be the Fourier transform of $`c_{[1]}`$ and $`\gamma `$ respectively. The power spectrum of the galaxy distribution in the generalized halo model is then

$$\begin{array}{c}P(k)=\frac{(2\pi )^3}{\varrho _p\mu ^2}|\stackrel{~}{c}_1(k)|^2+\frac{(2\pi )^6}{\mu ^2}|\stackrel{~}{c}_{[1]}(k)|^2P^{(p)}(k)+\hfill \\ \hfill +\frac{(2\pi )^3}{\varrho _p\mu ^2}d𝐤^{}|\stackrel{~}{c}_{[1]}(k^{})|^2\stackrel{~}{\gamma }(|𝐤𝐤^{}|).\end{array}$$
(85)
This first two terms are the result of , the additional term accounts for halo–substructure and involves a folding of $`\stackrel{~}{c}_{[1]}`$ with $`\stackrel{~}{\gamma }`$ in Fourier–space. Similar expressions can be derived from Eq. (77) for the bispectrum. Quantitative predictions for the galaxy distribution, similar to the investigations by , will be the topic of future work.

## VIII Some open problems

Our investigations rested on the assumption that the correlation functions exist and that the expansions of the p.g.fl. converge. In this case the p.g.fl., and consequently the point process, is determined completely by the correlation functions. The first assumption, the existence of the correlation functions (the factorial cumulants), does not impose dramatic restrictions for the models. In classical systems the mean number of points $`𝔼[N]`$ as well as the factorial moments $`𝔼[N(N1)\mathrm{}(Nn+1)]`$ should be finite in any bounded domain. For the $`n`$–point Poisson cluster processes, discussed in the preceding sections, at maximum $`n`$ points reside in a cluster, which are themselves distributed according to a Poisson process with constant number density. Clearly in such a simple situation both assumptions are satisfied. However, even for physically well motivated models, the convergence of the expansion of the p.g.fl. may not be guaranteed, although the point process itself and the correlation functions are well–defined.

Perhaps the best known example of a probability distribution which is not fully specified by its moments is the log–normal distribution. The probability density of a log–normal random variable is given by

$$p(x)=\frac{1}{\sigma x\sqrt{2\pi }}\mathrm{exp}\left(\frac{(\mathrm{log}x\overline{x})^2}{2\sigma ^2}\right)$$
(86)
with parameters $`\overline{x}`$ and $`\sigma ^2`$, the mean and variance of $`\mathrm{log}x`$. The moments (see Eq. (A1)) $`m_n=\mathrm{exp}(n\overline{x}+\frac{1}{2}n^2\sigma ^2)`$ are well–defined, however the expansion (A5) of the characteristic function is not convergent. And indeed  showed that the probability density

$$p^{}(x)=p(x)\left(1+ϵ\mathrm{sin}\left(\frac{2\pi k}{\sigma ^2}(\mathrm{log}x\overline{x})\right)\right)$$
(87)
where $`0<ϵ<1`$ and $`k`$ is a positive integer, has moments identical to the moments of the log–normal distribution. A comparison of $`p(x)`$ and $`p^{}(x)`$ is shown in Fig. 4.

A log–normal random field (an “exponentiated” Gaussian random field) is positive at any point in space, and a point process can be constructed using the value of the field as the local number density. The multivariate log–normal distribution, and the log–normal random field inherit the behaviour of the moments of the simple log–normal distribution. The point distribution obtained from the “Poisson sampled” log–normal random field is not characterized completely by its correlation functions as already discussed by . See also  for a similar approach towards this “log–Gaussian Cox process”.

In a Poisson cluster process (and also in the halo–model) the point distribution inside the cluster is specified independently from the distribution of the centers. This constructive approach, and the truncation of the moment expansion, guarantee the existence of these processes. A characterization result for the generalized Hermite distribution, closely related to the general $`n`$–point process considered in Sect. V D, is discussed by . Also well–defined point processes which do not impose such a truncation of the moment expansion are possible. A simple model is the line–segment process used in Sect. IV A, where the number of points per cluster is a Poisson random variable. Attempts towards a general characterization of point processes were conducted by , and partially succeeded for the case of infinitely divisible point processes.

One can show that any regular infinitely divisible point processes is a Poisson cluster process (e.g. , regular means that a cluster with an infinite number of points has probability zero). An infinitely divisible point processes may be constructed as a superposition of any number of independent point processes. It is interesting to note that the log–normal distribution is infinitely divisible , although the expansion of the characteristic function in terms of moments (A5) does not converge.

On small scales the galaxy correlation function is scale invariant: $`\xi _2r^\gamma `$. If a cut–off at some large scales is present, and the constraints for the density and the correlation functions are satisfied, a model based on a Poisson cluster process becomes feasible. Unfortunately, the superposition of independent point processes, as implied by the infinite divisibility of a Poisson cluster process, does not seem to be a good model assumption for the interconnected network of correlated walls and filaments, as observed in the galaxy distribution (e.g. ). The correlation functions for the galaxy distribution are close to zero for large separations, but from current observations one can not infer a definite cut–off. As discussed in Sect. III C for the Gauss–Poisson process, the large–scale behaviour of the correlation functions plays an important role in the construction of the Poisson cluster processes. Moreover, the dynamical equations governing the evolution of large–scale structures are non–local (see  and references therein). Therefore it seems worthwhile to consider also point process models which are not infinitely divisible. Unfortunately, beyond infinitely divisible point processes it is not clear what kind of properties the correlation functions have and especially what kind of additional constraints arise.

## IX Summary and conclusion

The Gaussian random field, fully specified by the mean and its correlation function, is one of the reference models employed in cosmology. Typical inflationary scenarios suggest that the primordial mass–density field is a realization of a Gaussian random field. Non–Gaussian features in the present day distribution of mass may be attributed either to the non–linear process of structure formation, or to a non–Gaussian primordial density field. Observations of the large–scale distribution of galaxies however provide us with a distribution of points in space. The process of galaxy formation may introduce further non–Gaussian features in the galaxy point distribution. In this paper a direct approach towards the characterization of this point set was pursued. The statistical properties of the point distribution can be specified by the sequence of correlation functions $`\xi _n`$. In close analogy to the Gaussian random field, a Gaussian point distribution, the Gauss–Poisson point process, was constructed. This random point set is fully specified by its mean number density $`\varrho `$, the two–point correlation function $`\xi _2(r)`$, and $`\xi _n=0`$ for $`n>2`$. Important constraints on $`\varrho `$ and $`\xi _2(r)`$, not present for the Gaussian random field, show up. Namely, $`\xi _2(r)0`$ for all $`r`$, and the variance of the number of points must not exceed twice the value of a Poisson process. The violation of these constraints indicates non–Gaussian features in the galaxy distribution. The equivalence of the Gauss–Poisson point process with a Poisson cluster point process leads to a simple simulation algorithm for such a point distribution. Using the $`J`$–function, higher–order correlations were detected in both a two–dimensional example and the galaxy distribution. The comparison with the Gauss–Poisson point process allows us to quantify the level of significance of these non–Gaussian features. Using these methods  could show that the distribution of galaxy clusters may not be modeled by a Gauss–Poisson process at a significance level of 95%.

The formal approach based on the probability generating functional (p.g.fl.) facilitated the definition, the characterization, and the simulation of the Gauss–Poisson point process. The inclusion of higher–order correlation functions was straightforward, leading to the $`n`$–point Poisson cluster processes. Both the definition and the simulation algorithm were detailed for the three–point Poisson cluster process. The Gauss–Poisson point process and the two–point Poisson cluster process are equivalent. However, this is not true for the $`n`$–point case $`n>2`$ anymore. The set of general $`n`$–point processes, resulting from a truncation of the cumulant expansion of the p.g.fl. after the $`n`$–th order, contains all $`n`$–point Poisson cluster processes as a true subset. This was discussed for the three–point case explicitly. Although models based on the $`n`$–point Poisson cluster process are not the most general ones, they cover a broad range of clustering point distributions. A Poisson cluster process can be simulated easily and is especially helpful for comparing statistical methods and estimators.

The inclusion of more and more points in the randomly placed clusters is only one way to extend the Gauss–Poisson point process into the strongly–correlated regime. In the halo model one allows for correlations between the cluster centers. Typically the halo (i.e. the galaxy cluster) is modeled without substructure. Again using the p.g.fl., the influence of correlations inside a halo on the $`n`$–point correlation functions of the resulting point distribution could be calculated.

All the models discussed above offer some insight into certain aspects of the clustering of the galaxy distribution. As argued in the preceding section, point process models which are not decomposable into independent point processes seem more appropriate. Unfortunately, even basic mathematical questions concerning the (complete) characterization of these models in terms of moments and beyond are still open.

### Acknowledgments

I would like to thank Lothar Heinrich, Klaus Mecke, Katja Schladitz, Peter Schücker, Alex Szalay, and Herbert Wagner for suggestions and helpful comments. Especially I would like to thank Istvan Szapudi for several interesting and beneficial discussions and Claus Beisbart and Alvaro Dominguez for discussions and extensive comments on the manuscript. I acknowledge support from the NSF grant AST 9802980 and from the Sonderforschungsbereich 375 für Astroteilchenphysik der DFG.

## A Characteristic and generating functions of random variables

A short review dealing with the characteristic and probability generating function (p.g.f.) of a random variable and their expansions in terms of moments, cumulants, factorial moments, and factorial cumulants is given. This Appendix is meant to serve as an illustration highlighting the analogies between expansions of the probability generating functional (p.g.fl.) and the probability generating function (see also ). To keep this summary simple it is assumed that the moments etc. exist and the expansions converge. For a more thorough treatment of characteristic and generating functions see e.g. .

The moments of a random variable $`x`$ with probability distribution $`F`$ are defined by

$$m_k=𝔼\left[x^k\right]=dF(x)x^k.$$
(A1)
If $`n`$ is a discrete random variable, especially if $`n`$ is integer–valued and greater equal zero, it is often more convenient to work with the factorial moments:

$$m_{[k]}=𝔼\left[n^{[k]}\right]=dF(n)n^{[k]}=\underset{n=0}{\overset{\mathrm{}}{}}p_nn^{[k]},$$
(A2)
where $`n^{[k]}=n(n1)\mathrm{}(nk+1)`$, and $`p_n`$ is the probability that the random variable takes the value $`n`$. Similarly, for point processes it is more convenient to work with product densities (or factorial moment measures), instead of moment measures.

The characteristic function $`\phi (t)`$, $`t`$ of a distribution $`F`$ is defined as

$$\phi (s)=_{\mathrm{}}^{\mathrm{}}dF(x)\mathrm{e}^{isx},$$
(A3)
and serves as a generating function for the moments. Expanding the exponential one can easily verify that

$$m_k=(i)^k\frac{\mathrm{d}^k\phi (s)}{\mathrm{d}s^k}|_{s=0}.$$
(A4)
By inverting one obtains the expansion of $`\phi (s)`$ in terms of moments:

$$\phi (s)=\underset{k=0}{\overset{\mathrm{}}{}}\frac{(is)^k}{k!}m_k.$$
(A5)
The probability generating function (p.g.f.) $`P(z)`$ of a random variable $`n`$ is defined as

$$P(z)=𝔼[z^n].$$
(A6)
Note that $`P(\mathrm{e}^{it})=\phi (t)`$. For a nonnegative integer–valued random variable one obtains the expansions

$`P(z)`$ $`={\displaystyle \underset{k=0}{\overset{\mathrm{}}{}}}p_kz^k,`$ (A7)
$`P(1+t)`$ $`=1+{\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{t^k}{k!}}m_{[k]},`$ (A8)
in terms of the probabilities $`p_n`$. $`P(1+t)`$ serves as the generating function for the factorial moments $`m_{[k]}`$. Similarly, the product densities (factorial moment measures) for a point process can be derived as functional derivatives (Frechet derivatives) of the probability generating functional.

Using $`P(\mathrm{e}^{it})=\phi (t)`$ one can derive the relation between moments and factorial moments:

$`\phi (s)`$ $`={\displaystyle \underset{k=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(is)^k}{k!}}m_k=P(\mathrm{e}^{is})`$ (A9)
$`=1+{\displaystyle \underset{t=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(is)^t}{t!}}{\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}m_{[k]}s_2(t,k),`$ (A10)
where

$$s_2(t,k)=\underset{l=1}{\overset{k}{}}(1)^{kl}\frac{l^t}{(kl)!l!}$$
(A11)
are the Stirling numbers of the second kind, the number of partitions which split the set $`\{1,\mathrm{},t\}`$ into $`k`$ pairwise disjoint nonempty sets (see e.g. ). Since $`s_2(t,k)=0`$ for $`k>t`$, which is also respected by the expression (A11), one finally arrives at

$$m_t=\underset{k=1}{\overset{t}{}}s_2(t,k)m_{[k]}.$$
(A12)
One considers not only the moments $`m_k`$ of a random variable, but also the cumulants $`c_k`$ defined by

$$\mathrm{exp}\left(\underset{k=1}{\overset{\mathrm{}}{}}\frac{(is)^k}{k!}c_k\right)=1+\underset{k=1}{\overset{\mathrm{}}{}}\frac{(is)^k}{k!}m_k=\phi (s).$$
(A13)
Clearly, $`\mathrm{log}\phi (s)`$ serves as a generating function for the cumulants $`c_k`$. Perhaps the best known cumulant is the variance

$$c_2=\sigma ^2=m_2m_1^2.$$
(A14)
Similarly, for nonnegative integer valued random variable the factorial cumulants $`c_{[k]}`$ are defined by

$$\mathrm{exp}\left(\underset{k=1}{\overset{\mathrm{}}{}}\frac{t^k}{k!}c_{[k]}\right)=1+\underset{k=1}{\overset{\mathrm{}}{}}\frac{t^k}{k!}m_{[k]}=P(1+t).$$
(A15)
Hence, $`\mathrm{log}P(1+t)`$ serves as the generating function of the factorial cumulants $`c_{[k]}`$. The correlation functions $`\xi _n`$ used in cosmology, are the densities of the normalized factorial cumulant measures of a point processes, corresponding to the factorial cumulants $`c_{[k]}`$ of a discrete random variable. Setting $`t=\mathrm{e}^{is}1`$ in Eq. (A15) and comparing term by term with Eq. (A13) one obtains the same relation between cumulants and factorial cumulants, as between moments an factorial moments:

$$c_k=\underset{l=1}{\overset{k}{}}s_2(k,l)c_{[l]}.$$
(A16)
## B Simulation algorithms

### 1 The Gauss–Poisson point process

As discussed in Sect. III B every Gauss–Poisson process is a Poisson cluster process and therefore can be simulated easily. For a given number density $`\varrho `$ and a two-point correlation function $`\xi _2(r)`$ fulfilling the constraints (15) and (16), realizations of the Gauss–Poisson process can be generated straightforwardly. With the normalization $`C_2=_^dd𝐱\xi _2(|𝐱|)`$ and $`_^dd𝐱f(|𝐱|)=1`$ one calculates the quantities needed for the simulation: $`f(r)=\xi _2(r)/C_2`$, $`q_2=\frac{\varrho C_2}{2\varrho C_2}`$, $`q_1=1q_2`$, and $`\varrho _p=\varrho (1\varrho C_2/2)`$. The constraint (19) now can be written as $`C_2\varrho 1`$. The simulation is carried out in two steps:

* First generate the cluster centers according to a Poisson distribution with number density $`\varrho _p`$.
* For each cluster center $`𝐱`$ draw a uniform random number $`z`$ in $`[0,1]`$. If $`z<q_1`$, then keep only the point $`𝐱`$. If $`zq_1`$ then keep the point $`𝐱`$ and additionally chose a random direction on the unit sphere and a distance $`d`$ with the probability density $`f`$ and place the second point according to them.

To get the correct point pattern inside a given window, one also has to use cluster centers outside the window to ensure that any possible secondary point inside the window is included.

### 2 The three–point Poisson cluster process

In the following the algorithm for the simulations of the three–point cluster process is given. The expressions (39) together with the normalization conditions for $`f_2`$ and $`f_3`$ serve as a starting point. An algorithm similar to the one for the Gauss–Poisson process described in Appendix. B 1 can be constructed:

Given the number density $`\varrho `$, and the two– and three–point correlation functions, $`\xi _2`$ and $`\xi _3`$, one defines

$$C_2=_^dd𝐱\xi _2(|𝐱|)\text{ and }C_3=_^dd𝐱_^dd𝐲\xi _3(0,𝐱,𝐲).$$
(B1)
Using the normalization of $`f_2`$ and $`f_3`$ one obtains

$`q_2`$ $`={\displaystyle \frac{3(C_2\varrho C_3\varrho ^2)}{(63C_2\varrho +C_3\varrho ^2)}},`$ (B2)
$`q_3`$ $`={\displaystyle \frac{C_3\varrho ^2}{(63C_2\varrho +C_3\varrho ^2)}},`$ (B3)
$`\varrho _p`$ $`={\displaystyle \frac{\varrho }{6}}(63C_2\varrho +C_3\varrho ^2),`$ (B4)
$`q_1`$ $`=1q_2q_3,`$ (B5)
resulting in

$`f_2(𝐱_1,𝐱_2)`$ $`={\displaystyle \frac{1}{C_2C_3\varrho }}(\xi _2(𝐱_1,𝐱_2)`$ (B6)
$`\varrho {\displaystyle _^d}\mathrm{d}𝐱_3\xi _3(𝐱_1,𝐱_2,𝐱_3)),`$ (B7)
$`f_2^{(3)}(𝐱_1,𝐱_2)`$ $`={\displaystyle \frac{1}{C_3}}{\displaystyle _^d}d𝐱_3\xi _3(𝐱_1,𝐱_2,𝐱_3),`$ (B8)
$`f_3(𝐱_1,𝐱_2,𝐱_3)`$ $`={\displaystyle \frac{1}{C_3}}\xi _3(𝐱_1,𝐱_2,𝐱_3).`$ (B9)
Since $`\varrho _c`$ and the $`q_n`$ are positive numbers, the constraints $`C_3\varrho ^2C_2\varrho 2+C_3\varrho ^2/3`$ must be satisfied and the relation $`\xi _2=C_2f_2+C_3\varrho (f_2f_2^{(3)})`$ holds. The algorithm now reads:

* First generate the cluster centers according to a Poisson distribution with number density $`\varrho _p`$.
* For each cluster center $`𝐱`$ draw a uniform random number $`z`$ in $`[0,1]`$. If $`z<q_1`$, then keep only the point $`𝐱`$. If $`q_1z<q_1+q_2`$ then keep the point $`𝐱`$ and additionally chose a random point according to the probability density $`f_2`$. If $`q_1+q_2z`$ then keep the point $`𝐱`$, chose a second point according to the probability density $`f_2^{(3)}`$, and a third point according to $`f_3`$.

## C Cumulants and factorial cumulants

Consider the expansion (7) of the p.g.fl. in terms of the factorial cumulants $`c_{[k]}(𝐱_1,\mathrm{},𝐱_k)`$:

$$\begin{array}{cc}\hfill \mathrm{log}G[h]& =\underset{k=1}{\overset{\mathrm{}}{}}\frac{1}{k!}_^dd𝐱_1\mathrm{}_^dd𝐱_k\hfill \\ & c_{[k]}(𝐱_1,\mathrm{},𝐱_k)(h_11)\mathrm{}(h_k1)\hfill \\ & =\underset{k=1}{\overset{\mathrm{}}{}}\frac{1}{k!}_^dd𝐱_1\mathrm{}_^dd𝐱_kc_{[k]}(𝐱_1,\mathrm{},𝐱_k)\hfill \\ & \left\{(1)^k+\underset{n=1}{\overset{k}{}}(1)^{kn}\underset{IJ_n^k}{}h_{I_1}\mathrm{}h_{I_n}\right\},\hfill \end{array}$$
(C1)
with $`h_1=h(𝐱_1)`$ and $`J_n^k`$ is formed by the ordered subsets of $`\{1,\mathrm{},k\}`$ with $`nk`$ distinct entries. Hence, a subset $`IJ_n^k`$ consists out of $`n`$ distinct numbers $`\{I_1,\mathrm{},I_n\}`$ with $`I_1<\mathrm{}<I_n`$, e.g $`J_2^3=\{\{1,2\},\{1,3\},\{2,3\}\}`$. Using $`\mathrm{\Phi }[h]=G[\mathrm{e}^{ih}]`$ one obtain

$$\begin{array}{ccc}\multicolumn{3}{c}{\mathrm{log}\mathrm{\Phi }[h]=\underset{k=1}{\overset{\mathrm{}}{}}\frac{1}{k!}_^d\mathrm{d}𝐱_1\mathrm{}_^d\mathrm{d}𝐱_kc_{[k]}(𝐱_1,\mathrm{},𝐱_k)\times }\\ & \times \left\{(1)^k+\underset{n=1}{\overset{k}{}}(1)^{kn}\underset{IJ_n^k}{}\mathrm{exp}\left(i(h_{I_1}+\mathrm{}+h_{I_n})\right)\right\},\hfill \\ \hfill =& \underset{k=1}{\overset{\mathrm{}}{}}\frac{1}{k!}_^d\mathrm{d}𝐱_1\mathrm{}_^d\mathrm{d}𝐱_kc_{[k]}(𝐱_1,\mathrm{},𝐱_k)\times \hfill \\ & \times \left\{(1)^k+\underset{m=0}{\overset{\mathrm{}}{}}\frac{i^m}{m!}\underset{n=1}{\overset{k}{}}(1)^{kn}\underset{IJ_n^k}{}(h_{I_1}+\mathrm{}+h_{I_n})^m\right\}.\hfill \end{array}$$
(C2)
In the $`m`$ sum the first term with $`m=0`$ equals $`(1)^k`$, canceling with the $`(1)^k`$ inside the braces:

$$\begin{array}{ccc}\multicolumn{3}{c}{\mathrm{log}\mathrm{\Phi }[h]=\underset{k=1}{\overset{\mathrm{}}{}}\frac{1}{k!}_^d\mathrm{d}𝐱_1\mathrm{}_^d\mathrm{d}𝐱_kc_{[k]}(𝐱_1,\mathrm{},𝐱_k)\times }\\ & \times \left\{\underset{m=1}{\overset{\mathrm{}}{}}\frac{i^m}{m!}\underset{n=1}{\overset{k}{}}(1)^{kn}\underset{IJ_n^k}{}(h_{I_1}+\mathrm{}+h_{I_n})^m\right\}\hfill \\ \hfill =& \underset{m=1}{\overset{\mathrm{}}{}}\frac{i^m}{m!}[\underset{k=1}{\overset{\mathrm{}}{}}\frac{1}{k!}_^d\mathrm{d}𝐱_1\mathrm{}_^d\mathrm{d}𝐱_kc_{[k]}(𝐱_1,\mathrm{},𝐱_k)\times \hfill \\ & \times \left\{\underset{n=1}{\overset{k}{}}(1)^{kn}\underset{IJ_n^k}{}(h_{I_1}+\mathrm{}+h_{I_n})^m\right\}].\hfill \end{array}$$
(C3)
The expression inside the braces $`[\mathrm{}]`$ in (C3) equals

$$\begin{array}{c}_^dd𝐱_1\mathrm{}_^dd𝐱_mc_k(𝐱_1,\mathrm{},𝐱_m)h_1\mathrm{}h_m.\hfill \end{array}$$
(C4)
This fixes the relation between the cumulants and factorial cumulants. However, there is no straightforward way to simplify this expression. Above all the theorem of Marcinkiewicz demands that as soon as $`k>2`$, the $`k`$–sum always has to be an infinite sum (see Sect. VI A). Hence, the cumulants $`c_k()`$ depend on an infinite alternating sum of the factorial cumulants $`c_{[k]}()`$, and vice versa.


format: text

text:
# Decoherence in circuits of small Josephson junctions

## Abstract

We discuss dephasing by the dissipative electromagnetic environment and by measurement in circuits consisting of small Josephson junctions. We present quantitative estimates and determine in which case the circuit might qualify as a quantum bit. Specifically, we analyse a three junction Cooper pair pump and propose a measurement to determine the decoherence time $`\tau _\phi `$.

>

The rate at which phase coherence is lost is an important factor in deciding whether a system qualifies as a quantum bit. Until now there exist very few experimental tests  or theoretical arguments  to determine this time in proposed quantum bits based on small Josephson junctions (squbits), and moreover only a lower bound of $`5`$ ns has been determined so far. In this letter we estimate quantitatively the dephasing time $`\tau _\phi `$ caused by the electromagnetic environment in single- and multijunction Josephson circuits both at zero and at non-zero temperatures. This dephasing is induced by the same fundamental processes which cause the decoherence in a squbit, thus yielding a direct measure of the decoherence time of a squbit . Based on our analysis we determine certain limiting factors for the realisation of a successful squbit experiment, and as a concrete example we investigate the three junction Cooper pair pump in which the coherent nature of the charge transport induces deviations from the accurate quantized transfer . We demonstrate that the crossover at $`f=f_C`$ as a function of the pumping frequency $`f`$ from incoherent to coherent charge transport yields a direct measure of $`\tau _\phi 1/f_C`$. We also discuss the limitations in implementing the so-called quasiparticle traps often used to suppress single electron effects in Cooper pair transistors and Cooper pair boxes .

We start by relating the dephasing rate of the Josephson junction circuit under consideration to the impedance of the electromagnetic environment that it is imbedded in. By dephasing we mean the deviation of the Josephson phase $`\phi (t)`$ across the junction circuit from its initial value $`\phi (0)`$. Specifically we are interested in the rms value of the phase deviation $`\sqrt{(\mathrm{\Delta }\phi )^2}\sqrt{[\phi (t)\phi (0)]^2}`$. Voltage fluctuations induced by dissipative circuit elements result in the phase-phase correlation function $`J(t)=\left[\phi (t)\phi (0)\right]\phi (0)`$ , which, based on the fluctuation-dissipation theorem, can be expressed in the form

$`J(t)`$ $`=`$ $`2{\displaystyle \underset{0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{d\omega }{\omega }}{\displaystyle \frac{\mathrm{Re}Z_t(\omega )}{R_K}}\times `$ (2)
$`\left\{\mathrm{coth}({\displaystyle \frac{\mathrm{}\omega }{2k_BT}})\left[\mathrm{cos}(\omega t)1\right]i\mathrm{sin}(\omega t)\right\},`$
where $`Z_t(\omega )`$ is the impedance seen by the circuit/junction whose phase fluctuations we want to determine, and $`R_K=h/e^225.8`$ k$`\mathrm{\Omega }`$ is the resistance quantum. It is straightforward to see that $`(\mathrm{\Delta }\phi )^2`$ and $`J(t)`$ are related by

$$\left(\mathrm{\Delta }\phi \right)^2=2\mathrm{R}\mathrm{e}J(t).$$
(3)
Other possible sources of decoherence  besides the dissipative environment will not be discussed here.

First we consider the case of an array of $`N`$ superconducting tunnel junctions. In the following analysis the array is assumed to be homogeneous, $`C_1=C_2=\mathrm{}=C_NC`$, and the electromagnetic environment to be purely resistive $`Z(\omega )=R_e`$ as shown in Fig. 1a. (Figure 1b refers to the three junction Cooper pair pump to be discussed later.) These simplifications make results more transparent, but it is straightforward to generalise the present method also to other circuits with an arbitrary environment, $`Z(\omega )`$. With the assumptions mentioned we obtain for the total resistance seen by the array  $`Z_t(\omega )=R_e/(1+i\omega \tau )`$, where $`\tau =R_eC/N`$. The real part of this can be written in form

$$\mathrm{Re}Z_t(\omega )=\frac{R_e}{1+\omega ^2\tau ^2}.$$
(4)
Inserting this into Eq. (2) and using the result for $`J(t)`$ when Re$`Z_t(\omega )`$ assumes the Lorentzian form of Eq. (4)  we obtain an expression for $`(\mathrm{\Delta }\phi )^2`$. In the limit of zero temperature ($`T0`$) we immediately get for $`t\tau `$

$$\left(\mathrm{\Delta }\phi \right)^2=4\frac{R_e}{R_K}\left[\mathrm{ln}\left(t/\tau \right)+\gamma \right],$$
(5)
where $`\gamma 0.57721`$ is Euler’s constant. In the case of non-zero temperature we consider only the long time ($`\pi k_BTt/\mathrm{}1`$) limit which is relevant in most cases, except in the limit of large $`R_e`$. At a realistic measurement temperature, e.g. $`T=50`$ mK, the result is valid in the range $`t50`$ ps, which is the region we are interested in. Long time expansion yields

$$\left(\mathrm{\Delta }\phi \right)^24\frac{R_e}{R_K}\left[\frac{\pi k_BT}{\mathrm{}}t\mathrm{ln}\left(\frac{2\pi k_BT\tau }{\mathrm{}}\right)+\gamma \right].$$
(6)
The long time expansion is valid only at non-zero temperatures and therefore Eq. (5) cannot be recovered from Eq. (6) in the limit of $`T0`$.

We can also apply the same method to an individual (the $`i`$th) junction to find the phase fluctuations $`(\mathrm{\Delta }\phi _i)^2`$ across it. With the same assumptions as before, we get for the single junction $`\mathrm{Re}Z_{t,i}(\omega )=\mathrm{Re}Z_t(\omega )/N^2`$ . This immediately yields the relation $`_i\sqrt{(\mathrm{\Delta }\phi _i)^2}=\sqrt{(\mathrm{\Delta }\phi )^2}`$, which can be shown to hold also with an arbitrary electromagnetic environment $`Z(\omega )`$ in series with the array.

If we define the dephasing time $`\tau _\phi `$ as the value of $`t`$ for which $`(\mathrm{\Delta }\phi )^2=(\pi /2)^2`$, we obtain for zero temperature:

$$\tau _\phi =\tau \mathrm{exp}\left(\frac{\pi ^2}{16}\frac{R_K}{R_e}\gamma \right),(T=0).$$
(7)
Dropping out the small constant terms in Eq. (6) we can write the result at finite temperature in the form:

$$\tau _\phi \frac{\pi }{16}\frac{\mathrm{}}{k_BT}\frac{R_K}{R_e},(T>0).$$
(8)
For an array with $`N=3`$ and $`C=10^{15}`$ F, and for the environment of $`R_e=1`$ k$`\mathrm{\Omega }`$ we obtain $`\tau _\phi 1.5`$ $`\mu `$s at zero temperature. A resistance of the environment of the order of the free space impedance $`R_e=Z_0377`$ $`\mathrm{\Omega }`$, yields $`\tau _\phi 1.510^5`$ s. With the same parameters at $`T=50`$ mK the decoherence is very fast: $`\tau _\phi 0.77`$ ns and $`\tau _\phi 2.1`$ ns for $`R_e=1`$ k$`\mathrm{\Omega }`$ and $`R_e=Z_0`$, respectively. Figure 2 shows the dependence of $`\tau _\phi `$ on the resistance of the environment, $`R_e`$, for a homogeneous three junction array at several different temperatures. Also the zero temperature limit (solid line) corresponding to Eq. (7), is shown. It forms an envelope curve for the finite temperature curves calculated from Eq. (6). It is also seen that only in the limit of low environmental resistance, $`R_e1`$ $`\mathrm{\Omega }`$, one can obtain long decoherence times $`\tau _\phi >1`$ $`\mu `$s at realistic measurement temperatures. Such values of $`\tau _\phi ^1`$ would possibly allow practical quantum logic operations and measurements to be performed by fast RF-gate lines and by using an RF-SET (radio frequency single electron transistor) as an electrometer .

In the previous analysis we considered an array of Josephson junctions without gates connected capacitively to the islands. These capacitances together with the impedance of the gate lines, $`Z_{gi}`$, also affect decoherence and should be taken into account. Therefore a quantitative analysis was also applied to the symmetric ($`C_iC`$, $`C_{gi}C_g`$ and $`Z_{gi}R_g`$ with all $`i`$) three junction Cooper pair pump which includes gate lines connected to the islands, as shown in Fig. 1b. The total impedance seen by the array is

$$Z_t(\omega )=\frac{R_e}{1+i\omega \tau \frac{6i\omega C+3g}{6i\omega C+2g}},$$
(9)
where $`g^1=R_g+1/(i\omega C_g)`$ is the impedance of the gate line and $`\tau =R_eC/3`$. The explicit expression for $`\tau _\phi `$ from Re$`Z_t(\omega )`$ does not assume a simple form but can be calculated numerically. Figure 3 shows the influence of the gate lines to the dephasing time as a function of $`R_g`$ with several different values of $`R_e`$. In the case of disconnected gate lines ($`R_g\mathrm{}`$) $`\tau _\phi `$ naturally approaches the value of the array without gates, as seen from the figure. In the limit of vanishing $`R_g`$ we also recover the result of an array by replacing $`C`$ by the effective capacitance $`C_{\mathrm{eff}}=C(C+C_g/2)/(C+C_g/3)`$ in Eq. (4).

The influence of dissipation in the gate lines on the dephasing rate is counterintuitive at the first sight: gate lines, even resistive ones, make the dephasing time longer than in an array without gates (Fig. 3). The reason behind this is twofold. Firstly, in our estimates we are interested in the fluctuations of the total phase difference across the array $`\phi `$, not in those of the individual phases $`\phi _i`$. Because of the series connection with additive phase differences, each gate line induces an exactly opposite, i.e. a cancelling fluctuation in the neighbouring junctions. Thus the noise of the gate resistors does not contribute to the noise in $`\phi `$. On the other hand, the gate lines decrease the impedance seen by the whole array, and this way the noise in the total phase also decreases. The longest dephasing time is therefore obtained with non-resistive gate lines. Due to the anticorrelated fluctuations in $`\phi _i`$ the sum rule for $`\sqrt{(\mathrm{\Delta }\phi _i)^2}`$, verified earlier for an $`N`$ junction array without gates, does not hold anymore.

A multijunction Josephson pump provides an interesting testground for quantum coherence  and, on the other hand, it may eventually qualify as a metrologically accurate current standard when coherence is suppressed by a very dissipative environment . Here we discuss the three junction pump (Fig. 1b) whose characteristics are determined by the two energies $`E_J`$, the Josephson coupling energy, and $`E_C`$, the charging energy. In the ideally coherent adiabatic regime, the phase across the array, $`\phi `$, is constant and no Landau-Zener band crossing occurs , and the optimum charge transfer through the array attains an approximate value (in the lowest order in $`E_J/E_C`$)

$$\frac{I}{2ef}19\frac{E_J}{E_C}\mathrm{cos}(\phi ).$$
(10)
Here $`I`$ is the current induced by operating the gates and $`f`$ is the frequency at which the system makes a wind around a degeneracy node of the charging energy, i.e. the frequency of the harmonic gate voltages $`V_{g1}`$. (The two gate voltages are phase-shifted by $`\pi /2`$.) The pumped charge per cycle, $`Q=I/f`$, is related but not equal to the geometric phase (Berry’s phase)  accumulated during one cycle along the closed path on the gate plane ($`V_{g1},V_{g2}`$). Contrary to the pump in the normal state , the coherent adiabatic Josephson pump lacks the ability to pump single charges virtually free of errors, and the relative deviations, $`9E_J/E_C\mathrm{cos}(\phi )`$, from the quantized transport are large even for very small values of $`E_J/E_C`$: for example, they can be as large as 9% for $`E_J/E_C=0.01`$. Yet, if the gates are operated slow enough, which means $`f1/\tau _\phi `$, the $`\mathrm{cos}(\phi )`$ term averages to $`\mathrm{cos}(\phi )=0`$ during one cycle or during the integration time of the measurement, and the pumping becomes accurate. Another limit comes from the Landau-Zener band crossing, which sets an upper limit for the operation frequency $`f_{LZ}E_J^2/(\mathrm{}E_C)`$. For typical $`E_J=0.1`$ meV and $`E_C=1`$ meV we obtain $`f_{LZ}10`$ GHz.

Based on these limitations we expect the following dependence of the pump performance at different frequencies (Fig. 4).

1. $`f<\tau _\phi ^1f_C`$: $`I/2ef1`$ because pumping is adiabatic but the phase is undetermined $`(\mathrm{cos}(\phi )=0)`$. Yet at the lowest frequencies the current becomes very small to measure and the accuracy will be lost in practise.
2. $`f_C<f<f_{LZ}`$: $`I/2ef19E_J/E_C\mathrm{cos}(\phi )`$, pumping is adiabatic and coherent.
3. $`f>f_{LZ}`$: $`I/2ef`$ decays because the condition for no band crossing is lost and charge transport does not follow the gating sequence adiabatically.

Since $`\tau _\phi `$ is presently expected to fall in the range $`\tau _\phi 5`$ ns, in a carefully designed experiment, we would have $`f_C200`$ MHz, yielding a clear separation of the three pumping regimes. In particular, if the decoherence time of a squbit and thus also $`\tau _\phi `$ turns out to be of the order of $`1`$ $`\mu `$s, which would allow quantum computation by Josephson qubits in this respect, $`\tau _\phi ^1`$ would give an experimentally convenient crossover frequency in the MHz range.

Our estimations of the decoherence time also bring up the issue of using so-called quasiparticle traps in single Cooper pair boxes and transistors. The parity effect has been observed over the years in several experiments . It is the manifestation of pure Cooper pair effects (not necessarily coherent) in small superconducting islands. Up to now, using aluminium structures, parity effect manifested by a 2e-periodic gate modulation, can be observed reliably only in Josephson junction circuits which are embedded in a resistive environment. Our estimates of $`\tau _\phi `$ now set the limit of how dissipative the quasiparticle traps can be, in order not to destroy the coherent state of the qubits too fast.

To perform an experiment of Fig. 4 with distinct regimes one has to have a setup with long enough dephasing time $`\tau _\phi 10`$ ns. This means that the on-chip resistances should be very low, which limits the use of quasiparticle traps. It has already been shown by experiment that with high enough resistances in the biasing circuit the pumped current becomes accurate . Combined with the fact that the parity effect, i.e. the quasiparticle free Cooper pair effect is very difficult to observe without quasiparticle traps, one needs to seek alternative measurement schemes of the pump. One way to avoid decoherence induced by a quasiparticle current is to fabricate an on-chip loop of a gated array terminated by an on-chip SET-electrometer . The electrometer could then be used to measure the number of Cooper pairs pumped through the array into it. However, it turns out that in this case the pumped current is accurate ($`I=2ef`$) due to charge conservation implied by the terminating classical capacitance.

Our suggestion is to realise the experiment shown in Fig. 4 by using a closed superconducting (phase) biasing circuit on the chip with an inductance in series. This way the real part of the impedance of the series environment, $`\mathrm{Re}Z_t(\omega )`$, vanishes. Further, gate lines do not induce any extra decoherence to the system, as shown before, and can therefore be as resistive as needed to filter the feed lines. Thus the major source of decoherence is the resistive impedance of the quasiparticle traps if needed. The pumped current can be measured by a SQUID ammeter inductively connected to the coil in the biasing circuit. This kind of a setup might give a low enough decoherence rate.

In conclusion, we have presented a method to quantitatively estimate the decoherence time, $`\tau _\phi `$, due to dissipative electromagnetic environment in circuits consisting of small Josephson junctions. This method allows us, among other things, to discuss the suitability of the system in consideration as a quantum bit. We also suggest a direct measurement of $`\tau _\phi ^1`$ as a crossover pumping frequency between coherent and incoherent pumping in the single Cooper pair pump.

Acknowledgments: This work has been supported by the Academy of Finland under the Finnish Centre of Excellence Programme 2000-2005 (Project No. 44875, Nuclear and Condensed Matter Programme at JYFL) and EU (contract IST-1999-10673). We thank Klavs Hansen and Matias Aunola for many useful discussions.