text:
# On the trispectrum as a gaussian test for cosmology

## 1 Introduction

Models of primordial fluctuations that generated the large-scale cosmological structures can be divided into two broad classes: gaussian and non-gaussian. The simplest versions of inflationary cosmology predict almost perfectly gaussian initial fluctuations, but such fluctuations can arise under more general conditions from the central limit theorem. The popularity of gaussian models is also due to their mathematical simplicity and the possibility they present for analytical calculations. On the other hand, density fields generated by topological defects (Vilenkin, 1985; Vachaspati, 1986; Hill et al., 1989; Turok, 1989; Albrecht & Stebbins, 1992) have non-gaussian initial conditions; the same can be said about fields generated by some versions of inflation (Allen et al., 1987; Kofman & Pogosyan, 1988; Salopek et al., 1989). Convincing evidence against gaussian initial conditions (GIC) would rule out many scenarios, and/or point us towards a physical theory for the origin of primordial fluctuations. Thus a test of the gaussian nature of the initial conditions is of great interest for cosmology.

Microwave background anisotropies probe cosmic fluctuations at a time when their statistical distribution should be close to its primeval form. To date, due to the limited signal-to-noise ratio of existing data, no conclusive evidence about the gaussianity of the initial conditions has been reached (e.g., Heavens (1998); Ferreira et al. (1998); Hinshaw et al. (1994); Falk et al. (1993); Luo & Schramm (1993); Gangui et al. (1994); Luo (1994); Smoot et al. (1994); Kogut et al. (1996); Banday et al. (1999); Tegmark & Bromley (1999)). An alternative is to analyze the present day density field of the galaxy distribution. This approach is complicated by the fact that the density field we observe today has already undergone non-linear gravitational evolution, and the observed non-gaussian nature of the galaxy distribution on small scales may be entirely the result of non-linear gravitational clustering combined with biasing effects. This non-linear growth makes any intrinsic non-gaussian signal more difficult to detect, and leads one to the conclusion that generally study of the microwave background is likely to be more profitable in detecting non-gaussian features (Verde et al., 2000c). However, one cannot exclude the possibility that non-gaussian features arise on physical scales which are difficult to probe with the microwave background, or non-GIC might produce a gaussian Sachs-Wolfe effect (e.g., Scherrer & Shaffer (1995)), so large-scale structure studies may still have a role to play in the test of the GIC hypothesis. An alternative method is to study number densities of rare/high-redshift objects (e.g., Chiu, Ostriker & Strauss (1998); Robinson, Gawiser & Silk (2000); Willick (2000); Matarrese et al. (2000); Verde et al. (2000a, b)).

There have been numerous studies to calculate the departures from gaussianity induced by gravity, and to set up a test of the gaussian nature of the initial conditions. There are essentially three different approaches: a) using N-body simulations, starting from GIC and several non-GIC, compare the resulting clustering properties of the evolved fields (e.g., Moscardini et al. (1991), Matarrese et al. (1991), Weinberg & Cole (1992)) b) investigate the topological properties of fields generated from gaussian and non-gaussian conditions and set up a comparison (e.g., Moscardini et al. (1990); Coles et al. (1993); Vogeley et al. (1994); Avelino (1997)) c) measure the moments of the density or the velocity fields and compare them with those predicted from a gaussian distribution or different non-gaussian models (e.g., Fry (1984); Fry & Scherrer (1994); Scherrer (1992); Catelan & Moscardini (1994); Catelan & Sherrer (1995), Gaztanaga & Fossalba (1998a,b), Kim & Strauss (1998); Bernardeau (1994)).

In this paper, we take the last approach, but work in Fourier space. We concentrate on the 4-point function, for the following reasons. Under rather general conditions, the two-point function carries no model-independent information about the gaussian nature of the field (Fan & Bardeen, 1994). The bispectrum (3-point function) is more promising, as it is zero for a gaussian field, but generally non-zero. However, it grows at second-order in perturbation theory, and this signal will place limits on the accuracy with which one can identify a primordial component (Verde et al., 2000c). The trispectrum (4-point function) has the advantage that it grows linearly, with contributions from non-linear growth appearing only weakly in perturbation theory. This will be quantified in Sec. 2. One might hope, therefore, that the linear growth extends to relatively small scales. In addition, the possibility is open that non-gaussian models have no bispectrum, or even negative initial bispectrum (e.g., Moscardini et al. (1991)).

One can also probe the 4-point function in real space (e.g., Luo & Schramm (1993); Lokas et al. (1995); Chodorowski & Bouchet (1996)), or look at a subset of configurations of the four-point function, through power correlations (Feldman et al., 1994; Stirling & Peacock, 1996). An alternative way to approach mildly non-gaussian fields is based on the Edgeworth expansion (Amendola, 1996).

The main difficulty with these approaches appears if the smoothed galaxy distribution is related to the underlying mass distribution by a non-linear transformation. In this circumstance, it is probably very difficult to distinguish non-GIC from a non-linear bias, so one has to assume that, at least on large scales, the bias is linear and deterministic (cf., Taruya et al. (2000)). In any eventuality, there is a realistic possibility that current large galaxy surveys such as the Anglo-Australian 2-degree Field (hereafter 2dF; Colless (1996)) and the Sloan Digital Sky Survey (York et al., 2000) will place some constraints on any initial departures from gaussian behavior.

This paper is organized as follows: In section 2 we review the relevant statistical properties of the density field through its $`n`$-point functions, and introduce the $`\tau `$ statistic as a non-gaussian discriminant. Section 3 contains the practical implementation of the method on numerical simulations. Redshift-space distortions are considered in Section 4 and conclusions are presented in Section 5.

## 2 Statistics of the density field in the linear and weakly non-linear regime

The statistical properties of the fractional overdensity field $`\delta (𝐱)=[\rho (𝐱)\overline{\rho }]/\overline{\rho }`$, can be characterized by the $`n`$-point correlation functions or, in Fourier space, by the $`n`$-point spectra. If the fluctuation field is gaussian, the connected part of the $`n`$-point function vanishes for $`n3`$ (e.g., Bertschinger (1992)). Thus the two-point function, or alternatively the power spectrum, completely specifies a gaussian distribution.

To linear order in perturbation theory an initially gaussian distribution remains gaussian; in particular the the connected $`n`$-point functions and $`n`$-point spectra of an initially gaussian distribution are zero as long as linear perturbation theory holds. On the other hand, if the initial conditions are non-gaussian, the $`n`$-point spectra in the linear regime are the primordial ones, scaled by the $`n^{th}`$-power of the linear growth factor. Therefore it would be possible to detect primordial non-gaussianity with a measurement of non-zero connected linear $`n`$-spectra for $`n3`$.

We concentrate here on the four-point function in Fourier space (the trispectrum), because, as we demonstrate below, second-order contributions vanish for gaussian primordial fluctuations, so we hope to be able to use linear theory to smaller scales. This is tested numerically. Moreover this quantity can be made independent of the cosmological model, the degree of linear bias, and the power spectrum at least on large scales as we will show below.

The data which we will work with are products of four Fourier coefficients (in practice we use the real part of this):

$$D_\alpha =\delta _{𝐤_1}\delta _{𝐤_2}\delta _{𝐤_3}\delta _{𝐤_4}.$$
(1)
The mean value of this is zero (by homogeneity), unless the four wavevectors form a quadrilateral. $`\alpha `$ labels the set of 4 wavevectors.

Note that the trispectrum itself is strictly the Fourier counterpart of the connected part of the four-point correlation function only. For a zero-mean field with non-zero connected (subscript $`c`$) four-point function,

$`\delta _{𝐤_1}\delta _{𝐤_2}\delta _{𝐤_3}\delta _{𝐤_4}`$ $`=`$
$`\delta _{𝐤_1}\delta _{𝐤_2}\delta _{𝐤_3}\delta _{𝐤_4}`$ $`+`$ $`\text{(2 perms.)}+\delta _{𝐤_1}\delta _{𝐤_2}\delta _{𝐤_3}\delta _{𝐤_4}_c`$ (2)
where

$`\delta _{𝐤_i}\delta _{𝐤_j}`$ $`=`$ $`(2\pi )^3P(k_i)\delta ^D(𝐤_i+𝐤_j)`$ (3)
$`\delta _{𝐤_1}\delta _{𝐤_2}\delta _{𝐤_3}\delta _{𝐤_4}_c`$ $`=`$ $`(2\pi )^3T(𝐤_i)\delta ^D(𝐤_1+𝐤_2+𝐤_3+𝐤_4).`$
Here $`P(k)`$ is the power spectrum, $`T`$ is the trispectrum, $`\delta ^D(𝐤)`$ is the three-dimensional Dirac delta function, the angle brackets indicate the ensemble average and $`(𝐤_i)`$ is a shorthand for $`(𝐤_1,𝐤_2,𝐤_3,𝐤_4)`$.

Note that $`\delta _{𝐤_1}\delta _{𝐤_2}\delta _{𝐤_3}\delta _{𝐤_4}`$ has a specific volume dependence when one works with a discrete rather than continuous Fourier transform. The gaussian part involves products of two Dirac delta functions, whereas the connected part has only one. For a volume-limited survey of volume $`V`$, the discrete Fourier transform changes these delta functions to $`\delta ^DV/(2\pi )^3`$ multiplied by a Kronecker delta, so any comparison between the size of non-gaussian and gaussian parts of $`D_\alpha `$ is volume-dependent. Any statement about the relative importance of the terms must therefore be made with some caution. We have to stress here that this is a general feature of $`\delta _{𝐤_1}\mathrm{}\delta _{𝐤_n}`$ that are evaluated in a finite volume: the relative importance of the connected and non-connected parts is volume dependent.

In order to demonstrate the lack of contribution to $`D_\alpha `$ from second-order perturbation theory, we expand the density field to second order as:

$$\delta (𝐱)=\delta ^{(1)}(𝐱)+\delta ^{(2)}(𝐱)$$
(4)
where $`\delta ^{(1)}`$ is the linearly-evolved overdensity. $`\delta ^{(2)}`$ is $`O(\delta ^{(1)2})`$ and, if initial conditions are gaussian, represents departures from gaussian behavior due to gravitational evolution (e.g., Goroff et al. (1986), Bouchet et al. (1992)). Strictly speaking the galaxy distribution is not a mildly non-linear field, but a highly non-linear field filtered on some smoothing scale. The operations of smoothing and evolution do not commute, but experiments (e.g., Matarrese et al. (1997), hereafter MVH97) show that perturbation theory works well if the smoothed field is not too non-linear.

To second order in perturbation theory the Fourier counterpart of the connected four-point correlation function can be obtained by applying to eq. (2) the substitution

$$\delta _𝐤=\delta _𝐤^{(1)}+ϵ\delta _𝐤^{(2)}+ϵ^2\delta _𝐤^{(3)}+\mathrm{}$$
(5)
and retaining only terms $`ϵ`$. Here $`ϵ`$ is simply a bookkeeping parameter that can then be set to unity and the expression for $`\delta _𝐤^{(2)}`$ can be found e.g., in Fry (1984) and Catelan et al. (1995). Hence

$`\delta _{𝐤_1}\delta _{𝐤_2}\delta _{𝐤_3}\delta _{𝐤_4}`$ $`=`$ $`\delta _{𝐤_1}^{(1)}\delta _{𝐤_2}^{(1)}\delta _{𝐤_3}^{(1)}\delta _{𝐤_4}^{(1)}+`$ (6)
$`\delta _{𝐤_1}^{(1)}\delta _{𝐤_2}^{(1)}\delta _{𝐤_3}^{(1)}\delta _{𝐤_4}^{(2)}+cyc.\text{ (4 terms)}.`$
There are 4 cyclical terms involving the second order $`\delta _𝐤^{(2)}`$ whose expression involves products of two linear coefficients, so for an initially gaussian field, the second-order contributions to the trispectrum are products of five coefficients and hence vanish. The first non-vanishing contribution to the connected trispectrum due to gravitational instability is $`ϵ^2`$ as is the next (third order) term in (5). In general, the trispectrum will depend on all orders of correlations present initially, but in the linear regime, the only contribution is from the initial trispectrum ($`T_{initial}`$). Specifically, if $`𝒟(t)`$ is the growth factor:

$$T=T_{initial}𝒟(t)^4.$$
(7)
For GIC, the trispectrum is zero, but $`D_\alpha _{GIC}`$ has a disconnected part

$$D_\alpha _{GIC}=(2\pi )^6P(k_1)P(k_3)\delta ^D(𝐤_1+𝐤_2)\delta ^D(𝐤_3+𝐤_4)+cyc.$$
(8)
If primordial fluctuations are non-gaussian, then $`D_\alpha D_\alpha _{GIC}`$. For $`k`$-vector configurations where all $`|𝐤|`$ are different, $`D_\alpha _{GIC}`$ vanishes, leaving $`D_\alpha =(2\pi )^3T(𝐤_1,𝐤_2,𝐤_3,𝐤_4)\delta ^D`$. We can thus characterize the departure from gaussian behavior by

$`\delta _{𝐤_1}\delta _{𝐤_2}\delta _{𝐤_3}\delta _{𝐤_4}`$ (9)
$`=`$ $`(2\pi )^6P(k_1)P(k_2)\delta ^D(𝐤_1+𝐤_3)\delta ^D(𝐤_2+𝐤_4)+cyc.`$
$`+`$ $`\tau \sqrt{P(k_1)P(k_2)P(k_3)P(k_4)}\delta ^D(𝐤_1+𝐤_2+𝐤_3+𝐤_4),`$
where we have introduced the quantity

$$\tau \frac{T(𝐤_1,𝐤_2,𝐤_3,𝐤_4)}{\sqrt{P(k_1)P(k_2)P(k_3)P(k_4)}}.$$
(10)
$`\tau `$ is independent of the volume, the cosmological model, linear bias and redshift (or alternatively the power spectrum amplitude), but in general depends on the $`𝐤_i`$ (i.e. the scale and the configuration of the quadrilateral formed by the 4 $`k`$-vectors). Here for simplicity we will use equilateral k-vector configurations; with this choice isotropy demands that $`\tau `$ depends only on $`|𝐤|`$, for a given shape. We can then use a likelihood method to estimate $`\tau (k)`$, and deduce whether it is consistent with the gaussian result ($`\tau =0`$).

### 2.1 Likelihood

Our treatment here follows that of MVH97. For simplicity, we would like to be able to adopt a gaussian likelihood (italics will refer to the form of the likelihood function, not the initial statistics), which will be the a posteriori probability distribution for $`\tau `$ if we assume a uniform prior:

$`(\tau )`$ $`=`$ $`{\displaystyle \frac{1}{(2\pi )^{\frac{M}{2}}(\mathrm{det}C)^{\frac{1}{2}}}}\times `$
$`\mathrm{exp}\left[{\displaystyle \frac{1}{2}}{\displaystyle \underset{\alpha \beta }{}}(D_\alpha \mu _\alpha )C_{\alpha \beta }^1(D_\beta \mu _\beta )\right].`$
Here $`M`$ is the number of data, which have means $`\mu _\alpha (\tau )D_\alpha `$, and that can be evaluated using equation 9. The covariance matrix is $`C_{\alpha \beta }(D_\alpha \mu _\alpha )(D_\beta \mu _\beta )`$. If $`\tau =0`$ the likelihood analysis is effectively a $`\chi ^2`$ analysis.

In the general (i.e. non-gaussian) case we cannot justify the applicability of the gaussian likelihood mathematically from the central limit theorem, but numerical experiments (section 3) support its use when many modes are used.

The covariance matrix involves the eight-point correlation function, computable from Wick’s theorem and, for continuous and discrete fields by the methods of MVH97:

$`\delta _1\mathrm{}..\delta _8_{\mathrm{total}}=`$ (12)
$`\delta _1\delta _2\delta _3\delta _4\delta _5\delta _6\delta _7\delta _8+\mathrm{}`$ 105 terms
$`+\delta _1\delta _2\delta _3\delta _4\delta _5\delta _6\delta _7\delta _8+\mathrm{}`$ 280 terms
$`+\delta _1\delta _2\delta _3\delta _4\delta _5\delta _6\delta _7\delta _8+\mathrm{}`$ 210 terms
$`+\delta _1\delta _2\delta _3\delta _4\delta _5\delta _6\delta _7\delta _8+\mathrm{}`$ 35 terms
$`+\delta _1\delta _2\delta _3\delta _4\delta _5\delta _6\delta _7\delta _8+\mathrm{}`$ 28 terms
$`+\delta _1\delta _2\delta _3\delta _4\delta _5\delta _6\delta _7\delta _8+\mathrm{}`$ 56 terms
$`+\delta _1\mathrm{}\delta _8.`$
These expressions apply both to continuous fields and discrete point processes, and the averages on the r.h.s. are all connected parts only. Under the hypothesis of GIC, as long as linear theory holds111As shown above, second order contributions vanish for gaussian primordial fluctuations, therefore linear theory breaks down when third-order contributions become important., all terms involving the connected $`m`$-point correlation with $`m`$ greater or equal than two are non-zero only in the presence of shot noise. Moreover for equal-side configurations the terms involving $`\delta _l\delta _m\delta _n`$ are also zero, since one cannot make a triangle from any 3 of the wavevectors. For clarity we neglect shot noise in the main text, but the calculations have been performed including it.

For GIC, the only non-zero terms come from the two-point function, i.e. the power spectrum. We can test the gaussian hypothesis by including only these gaussian terms in the covariance matrix, reducing the process to $`\chi ^2`$ analysis. This should also be a good approximation for mildly non-gaussian fields (cf. Heavens (1998) for the microwave background). With this covariance matrix, a measurement of $`\tau `$ which is significantly non-zero would only rule out gaussian statistics, but would not necessarily recover the correct value of $`\tau `$, which would require use of the correct non-gaussian covariance matrix.

### 2.2 A priori estimation of the error on $`\tau `$

We can readily calculate the expected variance of our estimator for $`\tau `$. The calculation follows MVH97 so we only sketch the calculation here. Ignoring shot noise, for square configurations we have

$$\mu _\alpha =D_\alpha =(2\pi )^6P^2(\delta ^D)^2+(2\pi )^3\tau P^2\delta ^D$$
(13)
and the variance is:

$$\sigma _\alpha ^2=D_\alpha ^2=4P^4\left[(2\pi )^3\delta ^D\right]^4$$
(14)
since only 4 terms are non-zero out of the 105. Hence

$$\sigma _\tau ^2=\frac{^2\mathrm{ln}}{\tau ^2}_{\tau =0}=\underset{\alpha }{}\frac{\mu _\alpha (\tau =0)P^2}{\sigma _\alpha ^2}.$$
(15)
The number of uncorrelated squares in a thin shell in $`k`$-space of width $`\delta (\mathrm{ln}k)`$ is $`\pi k^3g\delta (\mathrm{ln}k)`$ where $`g=V/(2\pi )^3`$ is the density of states. Considering contributions from all the shells to the continuum limit we obtain

$$\sigma _\tau ^2=\frac{1}{2^3\pi ^2}\frac{1}{4V}_{k_{min}}^{k_{max}}k^3d(\mathrm{ln}k).$$
(16)
$`k_{max}`$ is set by the breakdown of second-order perturbation theory and $`k_{min}`$ by the size of the sample. Therefore the error on $`\tau `$ scales as $`\sqrt{V}`$ and the maximum signal-to-noise is obtained by splitting the volume into subunits. This procedure to reduce the error when dealing with higher-order spectral statistics is widely used in other fields such as signal processing222Where it is referred to as segment averaging. and fluid mechanics see e.g., Brillinger (1975), Brillinger & Rosenblatt (1967), Lii et al. (1976) and has been already applied in cosmology by Matarrese et al. (1997) and Verde et al. (1998). This appears to give something for nothing, but of course this is not true. The resolution of the paradox is given in section 2.3.

The minimum size of the samples may be set by the breakdown of perturbation theory, or possibly by aliasing difficulties if there is significant power on scales larger than the box.

Note that the analysis in this section, and therefore the scaling (16), holds only as long as the field is close to gaussian so that non-gaussian elements of the covariance matrix can be neglected. As we said before, if the field is close to gaussian we are justified in adopting a gaussian likelihood, however if the distribution appears to be highly non-gaussian this error calculation and volume dependence is no longer valid. Since effectively $`P(k)`$ is estimated from the data themselves, in principle one should propagate any error on $`P(k)`$ through into the error on $`\tau `$. However, since the cosmic variance for the trispectrum is much larger than that on the power spectrum, this effect should be negligible. We confirmed this with numerical experiments.

The square configuration we considered above is not the only possible one. We could also consider other configurations in which all the $`k`$-vectors in the quadruplet have the same modulus. These configurations can be parametrized by the angle $`\theta `$ between the first two vectors and $`90^{}\theta 0^{}`$: in the square case $`\theta =90^{}`$, in the ‘degenerate’ case $`\theta =0^{}`$. The above error analysis holds for all $`\theta >0^{}`$. For the ‘degenerate’ case of $`\theta =0^{}`$, there are modifications to the number of cyclic terms contributing to the means and covariance matrix, and also to the number of independent quadrilaterals. In this case, we have $`\mu _\alpha =2P^2[(2\pi )^3\delta ^D]^2+(2\pi )^3\tau P^2\delta ^D`$, $`\sigma _\alpha ^2=24P^4[(2\pi )^3\delta ^D]^4`$, and the number of uncorrelated quadruplets is $`2\pi k^3g\delta (\mathrm{ln}k)`$. These reduce the error by a factor $`\sqrt{2\times 2^4/24}=1.15`$, if $`k_{max}`$ remains the same.

### 2.3 Volume-splitting procedure

Volume-splitting appears to improve the signal-to-noise. We show here in outline that it gives no more information than relaxing the fixed quadrilateral shape and including more modes.

Consider the $`n`$ point spectrum and imagine that we wish to compute the expected error on the amplitude $`A`$: $`\delta _{k_1}\mathrm{}\delta _{k_n}=AS\delta ^D`$ where $`\delta ^DV`$ and $`S`$ is some fixed function. The leading term in the covariance, neglecting shot noise, will be $`P^n(\delta ^D)^n`$. As already seen before the density of states is $`gV`$. If we fix the configuration of the $`k`$-vectors, then the expected error on $`A`$, for a thin shell in $`k`$-space, will be given by:

$$\sigma _A^2\frac{gS^2(\delta ^D)^2}{P^n(\delta ^D)^n}$$
(17)
which implies the curious result that $`\sigma _AV^{(n3)/2}`$. However, if we allow the shape of the $`k`$-vector configuration to change we obtain:

$$\sigma _A^2\frac{g^{(n1)}S^2(\delta ^D)^2}{P^n(\delta ^D)^n}$$
(18)
where the exponent $`n1`$ of the density of states comes from the presence of the Dirac delta function, i.e. the fact that the polygon has to close and therefore has $`n1`$ degrees of freedom. In this case we obtain the (more intuitive result) that: $`\sigma _A1/\sqrt{V}`$ when all shapes are allowed for.

In the case of the power spectrum (as in the example of Press et al. (1992)) the two possibilities coincide ($`n`$=1), for this reason splitting the volume does not improve the errors. Note also that in the bispectrum case ($`n`$=3) of Matarrese et al. (1997) and Verde et al. (1998), the error is independent of the size of the sample because the triangle shape is fixed. For the trispectrum ($`n`$=4) $`\sigma _A\sqrt{V}`$. Here, as in all cases, one has the choice between considering all shapes, or keeping the shape fixed and subdividing.

## 3 N-body tests

In this section, we perform tests with real-space results from N-body simulations. We analyze square configurations, and ensure that no wavevector occupies more than one square, so the covariance matrix is approximately diagonal (cf. MVH97 for the bispectrum). The simulation is a $`128^3`$ particle, 100 $`h^1`$ Mpc side box, CDM-like simulation from the Hydra consortium (Couchman et al., 1995), with parameters $`\mathrm{\Omega }=1.0`$, $`\mathrm{\Lambda }=0.0`$, $`\sigma _8=0.64`$, $`\mathrm{\Gamma }=0.25`$, and GIC. Shot noise is completely negligible in all scales of interest ($`k<1`$; here and hereafter $`k`$ is in units of h/Mpc).

The breakdown of linear theory for the trispectrum is not known in advance, but from MVH97 one can expect that the leading-order corrections for the trispectrum (i.e. contributions from the third order in perturbation theory) should be small at least up to $`k=0.55`$ h Mpc-1. A likelihood computation (cf., section 2.1) shows that for this square configuration, linear theory for the 4-point function should be valid up to $`k_{\mathrm{max}}=0.67`$, which sets the upper limit for our further analysis, giving 1250 squares in total. This limit is also apparent in Fig. 1, which plots 333The notation here follows MVJ97: $`D_\alpha `$ is a single measurement of a statistical quantity i.e. the trispectrum for a given (e.g., square) $`k`$ configuration on a given scale (e.g., $`k`$=0.5). $`D_\alpha `$ is the ensemble average of this quantity. Since we assume that the volume under consideration is a “fair sample”, then the average of all $`D_\alpha `$ in the simulation on the same scale and with the same $`k`$ configuration is an estimate of $`D_\alpha `$.$`D_\alpha `$. For this wavenumber limit, $`\chi ^2`$ is shown as a function of $`\tau `$ in Fig. 2: the minimum reduced $`\chi ^2`$ is 0.8 and $`\tau =(1.5\pm 5.5)\times 10^4`$ is consistent with zero. The a priori estimation for the error of Eq.(16) yields $`\sigma _\tau =5.6\times 10^4`$. Note that the estimate of $`\tau `$ is independent of linear biasing because of the presence of the same number of $`\delta `$ factors in the numerator and denominator of Eq. 10.

We have also analyzed another N-body simulation with very non-GIC. The initial conditions were set by applying the following mapping to a gaussian field:

$$\delta \delta ^2\delta ^2$$
(19)
(the so-called $`\chi ^2`$ model for initial conditions). The linear power spectrum is a power law with spectral index $`n=1`$. The simulation we used has $`80^3`$ particles, in a 200 $`h^1`$ Mpc side box, with $`\sigma _8=0.5`$. In this simulation the shot noise is not negligible, but it can be included with the methods of MVH97. This family of initial conditions could arise from evolution through inflation of an isocurvature CDM model for structure formation and is motivated by the fact that it accommodates galaxy formation at high redshift (e.g., Peebles (1999a, b)). Despite the low normalization, in this case we find that for all $`k`$-space shells with $`k<0.67`$ the minimum $`\chi ^2`$ value for $`\tau `$ is nonzero at high significance: considering all $`k`$-vectors in the range $`0.15<k<0.67`$ for equilateral configurations and in the range $`0.15<k<0.75`$ for degenerate configurations, we obtain $`\tau _{eq}=(9.2\pm 1.6)\times 10^5`$ and $`\tau _{deg}=(6.4\pm 1.3)\times 10^5`$ respectively. Thus we can reject the GIC hypothesis with confidence, but note that for this highly non-gaussian model the covariance matrix is not accurate, so we can not rely on the actual measurement of $`\tau `$ being reliable.

Incidentally, for this simulation, the bispectrum also differs in the mildly non-linear regime from the one that would have originated by gravitational instability from GIC. However, at the bispectrum level, there is degeneracy between biasing and initial non gaussianity: in practice it would not be possible to assess if there is substantial bias/anti-bias or if the initial conditions were truly non-gaussian. As we have mentioned, the trispectrum method is bias-independent, at least as long as the bias is not strongly non-linear.

### 3.1 Subdividing the Volume

We have tested the subdivision procedure with static simulations of a non-gaussian density field derived from a gaussian $`\delta (x)`$

$$\delta (x)\delta (x)+ϵ\delta ^2(x)$$
(20)
where $`ϵ`$ is some parameter. For simplicity, we take a white noise power spectrum between $`k_{\mathrm{min}}`$ and $`k_{\mathrm{max}}`$. Details of the resulting 2-point and 4-point functions are given in the Appendix, but they are summarized as

$$P(k)=P_g(k)+2ϵ^2P_g(k)^2V_k/(2\pi )^3$$
(21)
where $`P_g(k)`$ is the power spectrum of the underlying gaussian $`\delta `$ field, and $`V_k`$ is a $`k`$-space volume defined in the appendix. The 4-point function is

$$\delta _{𝐤_1}\delta _{𝐤_2}\delta _{𝐤_3}\delta _{𝐤_4}_{\mathrm{conn}.}(48ϵ^2P_g^3+ϵ^4P_g^4V_k/(2\pi )^3)V.$$
(22)
We choose $`P_g=2(2\pi )^3`$ in the range $`0.1<k<0.837`$ and zero elsewhere, and $`ϵ=0.4`$. In this case, the non-gaussian terms in the covariance matrix are small in comparison with the gaussian one and we ignore them. We split the volumes into cubes of side $`20`$ and $`24h^1`$ Mpc, and make repeated simulations to reduce errors. In all, we analyzed 460 of the smaller boxes and 266 of the larger boxes, so the total volume is the same in the two cases. The results are shown in Fig. 3. The expected value of $`\tau `$ is 1040, and is recovered within the errors. Note also that the errors (160 and 250 respectively) scale as expected $`N`$ where $`N`$ is the number of boxes, and the signal-to-noise is indeed higher for the ensemble of smaller boxes, as expected.

## 4 Complications with a real survey

In section 3 we showed that the trispectrum is a useful discriminant between GIC and non-GIC in a very idealized case where the field is unbiased, and the positions of the particles are known in real space. Note however that galaxy catalogues have an average number of galaxies per unit volume that is about two orders of magnitude smaller than the one of the simulation used, and use the redshift as a third spatial coordinate. The resulting redshift-space map of the galaxy distribution is therefore distorted and shot noise is significant. This is not a trivial issue since we are pushing into the mildly non-linear regime, which is significantly affected by non-perturbative contributions.

In this section, we show how to deal with these effects.

### 4.1 Redshift-space distortions

It is a convenient approximation to split redshift-space distortions into two components (Kaiser, 1987): a large scale distortion responsible of the squashing known also as the ‘bull’s eye’ effect and a small-scale radial smearing responsible for ‘Fingers-of-God’. The large-scale effect on individual Fourier components of the density fluctuation is well described by multiplication by the Kaiser factor $`(1+\beta \mu ^2)`$, where $`\mu `$ is the cosine of the angle between the $`k`$-vector and the line of sight, and $`\beta \mathrm{\Omega }_0^{0.6}/b`$ with $`b`$ the linear bias parameter.

The small-scale effect is hard to treat exactly and could potentially erase the signal we are trying to detect; in fact virialized motions on small scales produce a radial smearing and the associated Finger-of-God effect contaminates the wavelengths we are interested in. A successful model that fits the power spectrum in numerical simulation reasonably well (e.g., Hatton & Cole (1998); Peacock & Dodds (1994)) assumes that the small-scale velocity field is uncorrelated with density and has an exponential velocity distribution. Although not rigorously theoretically-motivated, it has been shown in VHMM98 that this modeling of the small-scale velocity dispersion works well as an addition to perturbation theory. In the Fourier domain, and in the distant observer approximation, the exponential velocity dispersion gives a damping factor, that, combined with the boosting factor of the large-scale effect (Kaiser, 1987) gives (Peacock & Dodds, 1994):

$$\delta _𝐤\delta _𝐤\frac{(1+\beta \mu ^2)}{\sqrt{1+k^2\sigma ^2\mu ^2/2}},$$
(23)
where $`\sigma `$ is the pairwise velocity dispersion of galaxies. With this model for redshift-space distortions:

$$D_\alpha (k)D_\alpha (k,\mu )=D_\alpha (k)\frac{(1+\beta \mu ^2)^4}{\sqrt{_{i=1,..,4}(1+k_i^2\sigma ^2\mu _i^2/2)}}.$$
(24)
As in VHMM98 we allow $`\sigma `$ to be scale-dependent (to fit the power spectrum – there is some observational evidence for this; see Hamilton & Tegmark (2000)) and we reject the $`k`$-vectors aligned too closely with the line-of-sight (see Scoccimarro et al. (1999) for an alternative model).

We have tested the model (24) for square configurations by performing a $`\chi ^2`$ analysis for the parameter $`\tau `$ on the unbiased redshift-space catalogue created from the N-body simulation with GIC.

Following VHMM98 the covariance matrix is modified as follows. The power spectrum is first replaced by $`P(𝐤)P(k)(1+\beta \mu ^2)^2`$ and the resulting covariance matrix is then divided by $`\sqrt{_{i=1}^8(1+k_i^2\sigma ^2\mu _i^2/2)}`$, where the index $`i`$ runs over the 8 $`k`$-vectors that form the two squares.

The knowledge of the real-space power spectrum is required in this model: in a realistic application (like for example Sloan and 2dF) an accurate fit for the galaxy real-space power spectrum will be known. Thus for the present work we shall assume that the real-space power spectrum is known: only the velocity dispersion $`\sigma `$ needs to be determined, and this can be done by fitting the power spectrum. The limit of validity of the small-scale redshift-distortion model for the trispectrum can also in principle be determined, although in this case the limiting factor is the breakdown of perturbation theory. A likelihood analysis of the redshift-space Fourier modes can give $`\sigma `$ in a similar manner to VHMM98.

The result for the redshift-space analysis is shown in Fig. 4. Again the true value for $`\tau `$ is well recovered within the errors: $`\tau =(1.0\pm 5.5)\times 10^4`$ (reduced $`\chi ^2`$ value is 0.95) for the square configuration and $`\tau =(0.8\pm 4.5)\times 10^4`$, for the degenerate configuration in good agreement with the predicted error of $`4.3\times 10^4`$ obtained as outlined in section 2.2 .

As a test of the velocity dispersion, we consider the quadrupole-to-monopole ratio of $`D_\alpha `$, following its use with the power spectrum Hamilton (1992) and the bispectrum (VHMM98, Scoccimarro et al. (1998)). For the degenerate configuration all vectors have the same $`\mu ^2`$, so the quadrupole-to-monopole ratio can be simply defined by $`R_TD_\alpha ^{(2)}/D_\alpha ^{(0)}`$, where the quadrupole and monopole moments of $`T`$ are

$`D_\alpha ^{(2)}`$ $`=`$ $`{\displaystyle \frac{5}{2}}{\displaystyle _1^1}D_\alpha (k,\mu )(3\mu ^21)𝑑\mu `$
$`D_\alpha ^{(0)}`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle _1^1}D_\alpha (k,\mu )𝑑\mu .`$ (25)
In Fig. 5 we show how the model (24) reproduces the observed $`R_T`$ in the GIC simulation. We follow the same procedure as VHMM98, where $`\sigma `$ is allowed to be slightly scale-dependent, constrained to fit the power spectrum, and $`k^2\mu ^20.3`$. It would be possible to define a similar quantity also for non-degenerate configurations, but in this case $`D_\alpha `$ in redshift space would depend on two different angles to the line-of-sight which leads to more complexity.

For the non-gaussian isocurvature field of Section 3.1, we find a $`\tau `$ which is inconsistent with zero: e.g., for degenerate configurations $`\tau _{deg}=(6.4\pm 1.3)\times 10^5`$. For non-gaussian models, in general, the signal $`\tau `$ will be shape-dependent, (as for example in the particular case illustrated in section 5). This dependence on the $`k`$-vector configuration might therefore hold some extra information or allow stronger detection of a non-gaussian signal.

### 4.2 Expected performance from PSCz and Sloan

The PSCz survey (Saunders et al., 1998) is the largest nearly-all-sky survey, containing around 15000 galaxies. The space density of galaxies is not very high, so shot noise is important beyond about 40 $`h^1`$ Mpc for $`k=0.67`$. Assuming the initial field is close to gaussian, we use the gaussian covariance matrix, but with shot noise included.

Before being able to perform the trispectrum analysis on this survey however still there are a number of unresolved issues. The most important are the practical effects of the rapidly-varying selection function and the radial nature of the redshift-space distortions. In fact the technique described in section 4.1 assumes that the spherical nature of the distortion can be neglected and the sky can be considered flat (the so-called ‘distant-observer approximation’, or ‘plane-parallel approximation’). We must emphasize that this approximation does not hold for the PSCz, and one really has to go beyond the plane-parallel approximation, by using spherical harmonics for example (e.g., Heavens & Taylor (1995)).

In the case of the 2dF and Sloan surveys, due to the bigger volume available, the problem of the radial nature of the redshift-space distortions is much reduced, and the selection function issue can be tackled as illustrated in the appendix. The 1$`\sigma `$ error achievable on $`\tau `$ is very encouraging: $`1\times 10^4`$ for Sloan, if we restrict the analysis to the linear regime ($`k<0.3h`$) and divide the survey conservatively into 100 $`h^1`$ Mpc side cubes. This reduces to $`4\times 10^3`$ if wavenumbers up to $`k=0.7h`$ Mpc-1 are included. The error can be reduced even further by considering smaller boxes.

#### 4.2.1 Meaning of non zero value for $`\tau `$

Most inflationary models predict deviations from gaussianity of the form described by equation (26)–applied to the potential rather than the density fluctuation field; this kind of non-gaussianity affects mainly the bispectrum (e.g., Verde et al. (2000c)). An order of magnitude calculation shows that for the non-gaussian models analyzed in Verde et al. (2000c), the minimum deviation from gaussianity that can be detected with the large-scale-structure bispectrum is at least a factor of a few smaller than the one detectable with the large-scale-structure trispectrum; in these cases an analysis of the cosmic microwave background bispectrum is thus more promising. However, as already mentioned, primordial non-gaussianity might arise on physical scales which are difficult to probe with the microwave background, the large-scale-structure bispectrum might deviate from the perturbation theory prediction due to the presence of bias or antibias, and non-gaussian models might have no primordial bispectrum, or even negative initial bispectrum. It is nevertheless clear that other plausible non-gaussian models such as the $`\chi ^2`$ model described in section 3 can easily be ruled out with a trispectrum analysis of Sloan/2dF data sets.

## 5 Conclusions

We have presented a method for using the trispectrum, the four-point function in Fourier space, to discriminate between GIC and non-GIC from linear large-scale structure data. The advantage of this method is that for mildly non-gaussian fields the analysis in real space is independent of the power spectrum normalization, linear biasing, and cosmology, and, as a spectral method, the covariance matrix is more simply computed, in contrast to $`n`$point correlation functions. In redshift space, cosmology and bias enter only through the measurable quantity $`\beta `$, we show how to deal with redshift-space distortions (in the distant-observer approximation), we include in the calculations effects of shot noise following the method presented in MVH97 (although for clarity the shot noise contribution is ignored in the text and reported in the appendix) and it is straightforward to model the effects of varying selection function using the method of MVH97 (see the Appendix). The non-linear evolution of the trispectrum is expected to be rather weak in perturbation theory, so our expectation is that linear theory should provide a good description up to scales where the field becomes significantly non-linear. This is born out (and quantified) in simulations. Fig. (1) illustrates this point. The equilateral configurations bispectrum for the same simulation considered here, agrees with second-order perturbation theory up to $`k0.55`$ (MVH97). Fig (1) shows linear perturbation theory is adequate for square trispectrum configurations up to $`k0.67`$. We parameterize the departures from gaussian statistics by introducing the related quantity $`\tau `$, which is the 4-point correlation function in units of $`\sqrt{P(k_1)P(k_2)P(k_3)P(k_4)}`$. This quantity, in specific cases, can give us a meaningful measure of ‘non-gaussianity’. We compute the gaussian variance of $`\tau `$, which allows a test of the gaussian hypothesis. The error on $`\tau `$ can be computed straightforwardly, and is found to be in good agreement with internal error from N-body simulations. For mildly non-gaussian fields, we can expect that the use of a gaussian covariance matrix is still quite adequate, and in such cases a measurement of $`\tau `$ can reliably be made, subject to the requirement that the bias is linear on large scales. For highly non-gaussian fields, the gaussian hypothesis can be rejected, but the measurement of $`\tau `$ will be unreliable as the covariance matrix we use may neglect important terms.

Application to data will probably have to wait for completion of the 2dF and Sloan surveys, as the distant-observer approximation is not a good one for the nearly all-sky PSCz survey. In addition, the high shot noise in the PSCz survey means that the expected errors on $`\tau `$ are large enough to make it hard to distinguish between GIC and some significantly non-gaussian fields. The 2dF and Sloan surveys should do much better. One could compute $`\tau `$ from volumes of side 100 $`h^1`$ Mpc with a random error of less than one per cent for a $`\chi ^2`$ model as the one considered here, and could therefore put tight constraints on non-gaussian models. It would be helpful to apply this method to defect models, but unfortunately none is available at the required stage of evolution.

The major uncertainty in this analysis is the effect of bias. A non-linear bias term acts like a non-gaussian initial field, and, as we are using linear perturbation theory here, we cannot use polygons of different shape to lift the degeneracy, as done by MVH97 to lift the degeneracy between bias and gravitational evolution in the bispectrum. The best hope is to constrain a combination of initial non-gaussianity and quadratic bias, and to argue that if the measured value of $`\tau `$ is zero, then it would require a conspiracy unless both effects were absent. In principle it is possible to distinguish these two effects by considering differences at higher order in perturbation theory, as an Eulerian non-linear bias is applied to the evolved field, whereas a primordial non-gaussian field of the same mathematical form is applied to the initial field, and thus the evolution is different. It is an open question as to whether any realistic survey would have the signal-to-noise to do this.

## APPENDIX

We consider the transformation of a gaussian field as follows:

$$\delta (x)\delta (x)+ϵ\delta ^2(x)$$
(26)
This can be seen as the first two terms of a Taylor expansion of any non-gaussian field originating as a local mapping from a underlying gaussian one.

We choose a toy model for the power spectrum of the underlying gaussian field $`P_g`$: a top hat function between some $`k`$ minimum and maximum. The power spectrum for the resulting field will then be:

$`P(k)=P_g(k)+{\displaystyle \frac{2}{(2\pi )^3}}ϵ^2{\displaystyle P_g(k^{})P_g(𝐤𝐤^{})d^3k^{}}`$
$`P_g(k)+2ϵ^2P_g(k)^2V_k/(2\pi )^3`$ (27)
where $`V_k`$ is the volume where the integrand in the previous equation is non-zero. The relevant quantities for our analysis, substituting the Dirac delta function by $`V/(2\pi )^3`$ are the following:

$$\delta _{𝐤_1}\delta _{𝐤_2}\delta _{𝐤_3}\delta _{𝐤_4}_{\mathrm{conn}.}(48ϵ^2P_g^3+ϵ^4P_g^4V_k/(2\pi )^3)V,$$
(28)
$`\mathrm{SIGNAL}\tau {\displaystyle \frac{\delta _{𝐤_1}\delta _{𝐤_2}\delta _{𝐤_3}\delta _{𝐤_4}_{\mathrm{conn}.}V}{\delta _{𝐤_1}\delta _{𝐤_2}\delta _{𝐤_3}\delta _{𝐤_4}}}`$
$`{\displaystyle \frac{48ϵ^2P_g^3+ϵ^4P_g^4V_k/(2\pi )^3}{(P_g^2+4ϵ^2P_g^3V_k/(2\pi )^3+4ϵ^4P_g^4[V_k/(2\pi ^3)]^2)}}`$ (29)
$`\mathrm{NOISE}\sqrt{(2\pi )^3V}{\displaystyle \frac{2}{\sqrt{k_{max}^3k_{min}^3}}}`$ (30)
The covariance matrix is dominated by the gaussian terms if $`\delta _{𝐤_1}\delta _{𝐤_2}\delta _{𝐤_3}\delta _{𝐤_4}_{\mathrm{conn}.}\delta _{𝐤_1}\delta _{𝐤_2}\delta _{𝐤_3}\delta _{𝐤_4}`$, but at the same time we want to maximize the signal-to-noise ratio.

This considerations suggest that the best choice is: $`P(k)=2(2\pi )^3`$ if $`0.1<k<0.837`$ and zero elsewhere and $`ϵ=0.4`$. This allows us to use the gaussian covariance matrix, but also to have a reasonable signal-to-noise for the non-gaussian component from a single box of size of a few tens of Mpc.

### Shot Noise

In the application of the method to a real survey the shot noise might be important. Here we report how to modify the relevant expressions in the presence of shot noise. In what follows $`\overline{n}`$ denotes the average number of particles per unit volume, and the superscript $`d`$ indicates that these results apply to discrete, point processes. We follow the generating functional approach of MVH97 to find the following:

$$\delta _i\delta _j_c^d(2\pi )^3\left[P(k_i)+\frac{1}{\overline{n}}\right]\delta ^D(𝐤_i+𝐤_j)$$
(31)
$`\delta _l\delta _m\delta _n_c^d`$ $`(2\pi )^3\delta ^D(𝐤_l+𝐤_m+𝐤_n)\times `$ (32)
$`\left\{B_{lmn}+{\displaystyle \frac{1}{\overline{n}}}\left[P(k_l)+P(k_m)+P(k_n)\right]+{\displaystyle \frac{1}{\overline{n}^2}}\right\}`$ ,
where $`B_{lmn}`$ denotes the Bispectrum,

$`\delta _o\delta _p\delta _q\delta _r_c^d`$ $``$ $`(2\pi )^3\delta ^D(𝐤_o+𝐤_p+𝐤_q+𝐤_r)\times `$ (33)
$`\{{\displaystyle \frac{1}{\overline{n}}}[B_{(o+p)qr}+\text{perm. (6 terms)}]+`$
$`{\displaystyle \frac{1}{\overline{n}^2}}[(P_{o+p+q}+\text{cyc. (4 terms)})+`$
$`P_{o+p}+P_{o+q}+P_{o+r}]+`$
$`{\displaystyle \frac{1}{\overline{n}^3}}\},`$
$`\delta _1\mathrm{}\delta _5_c^d`$ $``$ $`(2\pi )^3\delta ^D(𝐤_1+𝐤_2+𝐤_3+𝐤_4+𝐤_5)\times `$ (34)
$`\{{\displaystyle \frac{1}{\overline{n}^2}}[B_{(1+2)(3+4)5}+\text{perm. (15 terms)}]+`$
$`{\displaystyle \frac{1}{\overline{n}^2}}[B_{(1+2+3)45}+\text{perm. (10 terms)}]\}+`$
$`{\displaystyle \frac{1}{\overline{n}^3}}[(P_{1+2}+\text{cyc. (10 terms)})+`$
$`+P_1+P_2+P_3+P_4+P_5]+`$
$`{\displaystyle \frac{1}{\overline{n}^4}}\},`$
$`\delta _1\mathrm{}\delta _6_c^d`$ $``$ $`(2\pi )^3\delta ^D(𝐤_1+\mathrm{}+𝐤_6)\times `$ (35)
$`\{{\displaystyle \frac{1}{\overline{n}^3}}[B_{12(3+\mathrm{}+6)}+\text{perm. (15 terms)}+`$
$`B_{(1+2)(3+4)(5+6)}+\text{ perm. (15 terms.)}+`$
$`B_{1(2+3)(4+5+6)}+\text{perm. (60 terms)}]+`$
$`{\displaystyle \frac{1}{\overline{n}^4}}[P_1+\mathrm{}+P_6+`$
$`P_{1+2}+\text{perm. (15 terms)}+`$
$`P_{1+2+3}+\text{ perm. (10 terms)}]+`$
$`{\displaystyle \frac{1}{\overline{n}^5}}\}.`$
$`\delta _1\mathrm{}\delta _8_c^d`$ $``$ $`(2\pi )^3\delta ^D(𝐤_1+\mathrm{}+𝐤_8)\times `$ (36)
$`\{{\displaystyle \frac{1}{\overline{n}^5}}[B_{12(3+\mathrm{}+8)}+\text{perm. (28 terms)}+`$
$`B_{1(2+3)(4+\mathrm{}+8)}+\text{perm.(168 terms)}+`$
$`B_{1(2+3+4)(5+\mathrm{}+8)}+\text{perm.(280 terms)}+`$
$`B_{(1+2)(3+4+5)(6+7+8)}+\text{perm.(280 terms)}+`$
$`B_{(1+2)(3+4+5+6)(7+8)}+\text{perm.(210 terms)}]+`$
$`{\displaystyle \frac{1}{\overline{n}^6}}[P_1+\mathrm{}+P_8+`$
$`P_{1+2}+\text{perm. (28 terms)}+`$
$`P_{1+2+3}+\text{ perm. (56 terms)}+`$
$`P_{1+2+3+4}+\text{ perm. (35 terms)}+`$
$`{\displaystyle \frac{1}{\overline{n}^7}}\}.`$
Neglecting shot noise when it is not negligible has two effects. The main effect is to overestimate $`H`$. E.g., the generalization of (8) is $`TT_{SN}`$, where

$`(2\pi )^3T_{SN}=`$
$`(2\pi )^3(P(k_1)`$ $`+{\displaystyle \frac{1}{\overline{n}}})(P(k_3)+{\displaystyle \frac{1}{\overline{n}}})\delta ^D(𝐤_1+𝐤_2)\delta ^D(𝐤_3+𝐤_4)`$
$`+\delta _1\mathrm{}\delta _4_c^d`$ (37)
The second effect is that the errors are underestimated. In particular, even in the gaussian case, there are additional connected terms in the correlations arising from shot noise.

### Selection function

The above analysis is valid for volume-limited samples, where the mean number density of galaxies is independent of position. In a realistic catalogue the presence of a varying selection function makes the mean density position-dependent and might even induce spurious detection of non-gaussianity. However it is possible to treat this effect accurately by the generating functional approach as in MVH97. If $`n(𝐱)`$ is the mean density of such a catalogue, following Feldman et al. (1994), we can define a working fluctuation field composed by subtracting from the real catalogue a synthetic catalogue with no clustering, but with the same selection function, and then weighting suitably the combination:

$$F(𝐱)=\gamma w(𝐱)[n(𝐱)\alpha n_s(𝐱)]$$
(38)
where $`w(𝐱)`$ is an arbitrary weighting function, $`\gamma `$ is a normalization factor, $`n_s(𝐱)=n(𝐱)/\alpha `$, and $`\alpha `$ is the “dilution” of the synthetic catalogue. We will then consider the limit $`\alpha 0`$ to avoid shot noise in the synthetic catalogue. The $`n`$-point correlation functions in Fourier space are easily calculated by considering $`F`$ as the superposition of the process $`f=\gamma w(𝐱)n(𝐱)`$ and $`f_s=\alpha \gamma w(𝐱)n_s(𝐱)`$. The generating functional for $`F`$ will then be: $`𝒵_F(𝒥)=𝒵_f(𝒥)𝒵_{f_s}(\alpha 𝒥)`$. Ignoring shot noise,

$$𝒵_{f_s}(\alpha 𝒥)=\alpha id^3x𝒥_s$$
(39)
where $`𝒥_s=\gamma w(𝐱)n_s(𝐱)_ms_m\mathrm{exp}(i𝐤_m𝐱)`$, and the ansatz  for the generating functional for the field $`f`$ is:

$$\begin{array}{c}𝒵_f[𝒥]=\\ \mathrm{exp}[id^3x𝒥(\stackrel{}{x})\frac{1}{2}d^3xd^3x^{}𝒥(\stackrel{}{x})𝒥(\stackrel{}{x}^{})\xi _{\mathrm{conn}.}^{(2)}(\stackrel{}{x},\stackrel{}{x}^{})\\ \frac{i}{6}d^3xd^3x^{}d^3x^{\prime \prime }𝒥(\stackrel{}{x})𝒥(\stackrel{}{x}^{})𝒥(\stackrel{}{x}^{\prime \prime })\xi _{\mathrm{conn}.}^{(3)}(\stackrel{}{x},\stackrel{}{x}^{},\stackrel{}{x}^{\prime \prime })\\ +\frac{1}{24}d^3xd^3x^{}d^3x^{\prime \prime }d^3x^{}𝒥(\stackrel{}{x})𝒥(\stackrel{}{x}^{})𝒥(\stackrel{}{x}^{\prime \prime })\mathrm{}𝒥(\stackrel{}{x}^{})\xi _{\mathrm{conn}.}^{(4)}(\stackrel{}{x},\mathrm{},\stackrel{}{x}^{})]\end{array}$$
(40)
Following the same procedure as in Matarrese et al. (1997) it is possible to obtain the $`n`$-point function in Fourier space by differentiating the generating functional; shot noise can be easily included by modifying the generating functional .

The $`\tau `$ in presence of a spatially-varying selection function therefore becomes:

$$\tau \tau \frac{I_{44}}{I_{22}^2}$$
(41)
where

$$I_{ij}d^3𝐱w^i(𝐱)n^j(𝐱)$$
(42)
In general the effect of a spatially-varying selection function can be summarized as follows:

$`\gamma ^2`$ $`=`$ $`1/I_{22}`$ (43)
$`V`$ $``$ $`{\displaystyle \frac{I_{NN}}{I_{22}^{N/2}}};\delta ^D{\displaystyle \frac{I_{NN}}{(2\pi )^3I_{22}^{N/2}}}`$ (44)
$`{\displaystyle \frac{1}{\overline{n}^q}}`$ $``$ $`{\displaystyle \frac{I_{N(Nq)}}{I_{NN}}}`$ (45)
## Acknowledgments

LV acknowledges NASA grant NAG5-7154 and TMR fellowship. We are grateful to Sabino Matarrese, Avery Meiksin, M. Strauss and the anonymous referee for helpful comments. The simulation with gaussian initial conditions was obtained from the data bank of cosmological N-body simulations provided by the Hydra consortium (http://coho.astro.uwo.ca/pub/data.html) and produced using the Hydra N-body code (Couchman et al., 1995). LV thanks Alison Stirling for providing the N-body simulation with non-gaussian initial conditions.

Fig. 1.–The 4-point correlation function in Fourier space (disconnected and connected) $`D_\alpha \delta _{𝐤_1}\mathrm{}\delta _{𝐤_4}`$ is approximated by averaging square configurations with binned wavevectors from a CDM-like N-body simulation (solid line). Also shown is the linear perturbation theory prediction for gaussian initial conditions i.e. $`\mu _\alpha (\tau =0)`$ from equation (9), for a square configuration of wavevectors. Errors are errors in the mean for each bin. For this square configuration, linear perturbation theory breaks down around $`k=0.67h`$ Mpc-1.

Fig. 2.–Minimum $`\chi ^2`$ analysis for the parameter $`\tau `$ for the gaussian CDM-like N-body simulation. Only the ‘square’ configuration for the trispectrum has been considered here and $`0.15k0.67`$ (in units of h Mpc-1). $`\tau `$ is a measure of the connected part of the trispectrum; for gaussian initial conditions $`\tau =0`$ (vertical line). To produce this graph, we use the gaussian initial conditions covariance matrix to compute the likelihood (equation 11). See sections 2.1 and 3 for further details.

Fig. 3.–The $`\chi ^2`$ analysis for 460 volumes of $`20h^1`$ Mpc side (dot-dashed line) and 266 volumes of $`24h^1`$ Mpc side (dashed line) for $`0.1<k<0.837`$ (in units of h Mpc-1). The $`y`$ axis is normalized so that 1-$`\sigma `$ limits are at $`y=1`$. The predicted value for $`\tau `$ is $`\tau 1040`$ (vertical lines). This illustrates the scaling with the volume of the error on $`\tau `$ (see Section 3.1 for more details).

Fig. 4.–Minimum $`\chi ^2`$ analysis for the parameter $`\tau `$ from the redshift-space unbiased CDM-like N-body simulation (see text for details). Only the square configuration has been considered here. The value for the velocity dispersion parameter $`\sigma 980`$km/s although it is slightly scale-dependent. The range of $`k`$-vectors considered is $`k0.67h`$ Mpc-1.

Fig. 5.– Quadrupole-to-monopole ratio for the trispectrum $`R_T`$ for degenerate square configuration of a GIC simulation. The continuous line is the theoretical $`R_T`$ using the redshift-space distortions as in VHMM98 (see section 4.1 for details).


format: text

text:
# Virial coefficients and equations of state for mixtures of hard discs, hard spheres and hard hyperspheres

## I Introduction

Depending on the nature of the independent variables, it is well known that the full thermodynamic description of a given system requires the availability of either one of the thermodynamic potentials or of the corresponding equations of state. In the case of classical fluids, the usual expression for the equation of state (EOS) is a relationship among the pressure $`p`$, the density $`\rho `$ and the temperature $`T`$ of the fluid. However, except for very few ideal systems, the explicit (exact) form of this relationship is not known in general. Therefore, given the fact that an experimental determination of the EOS for every particular fluid is of course not practical, researchers have concentrated their efforts in proposing “reasonable” empirical or semi-empirical approximations. Perhaps the neatest example of such an approach, already produced over a century ago and providing an essentially correct qualitative picture of the thermodynamic properties of gases and liquids, is the celebrated van der Waals EOS. An alternative proposal, first introduced by Kammerlingh Onnes as a mathematical representation of experimental results on the equation of state of gases and liquids, is the so called virial equation of state. Written for the compressibility factor $`Zp/\rho k_BT`$ (where $`k_B`$ is the Boltzmann constant) as a power series in $`\rho `$, it provides an expansion that gives the deviation from ideal gas behaviour in ascending powers of the density, namely

$$Z=1+\underset{n=2}{\overset{\mathrm{}}{}}B_n\rho ^{n1}.$$
(1)
Here, the $`B_n`$ are the virial coefficients. For a simple fluid, the virial coefficients are in general functions of $`T`$ alone, while in the case of mixtures they also depend on composition. An interesting aspect of the virial equation of state, discovered many years after it was introduced, is the fact that it can be derived rigorously using statistical mechanics. This in turn implies that the virial coefficients are not merely empirical constants but are rather related to intermolecular interactions in a well-defined manner. It is unfortunate, however, that in general the actual computation of the virial coefficients is a formidable task and that the radius of convergence of the series in Eq. (1) is not known.

In the case of hard-core fluids (rods, discs, spheres and hyperspheres) the virial coefficients are easier to compute than those corresponding to any other intermolecular interactions. Even in these model systems and with the exception of hard rods, where the exact EOS and corresponding virial coefficients are known, only the first few of such coefficients are available so far because the number of cluster integrals involved increases very rapidly with the order of the coefficient. The availability is scarcer for mixtures than for simple fluids, although recently there has been renewed interest in partially remedying this deficiency . The limited knowledge of the virial coefficients has proved useful, however, in the sense that proposals for the EOS of these systems may be judged amongst other things by how well the virial coefficients arising in such proposals compare with the exact values. The rationale here is that a “good” theoretical EOS should lead to an accurate prediction of the value of the virial coefficients. Alternatively, the few available virial coefficients may also be profitably used to construct a rational (Padé or Levin) aproximation to the EOS or to guide in the construction of a theoretical EOS, as was the case with the Carnahan-Starling (CS) EOS for a simple hard-sphere fluid.

A variety of (approximate) EOS for simple fluids composed of hard discs and hard spheres are available in the literature and some work has also been reported on simple fluids of hard hyperspheres . Very recently we introduced a simple recipe to derive the compressibility factor of a multicomponent mixture of $`d`$-dimensional hard hyperspheres in terms of that of the corresponding single component system. A straightforward consequence of such a recipe is that one can easily derive the explicit expressions for the virial coefficients of the mixture in terms of those of the single component fluid for all $`d`$. The major aim of this paper is to further assess the usefulness of our recipe by comparing the predictions of the values of the virial coefficients of binary mixtures of hard discs, hard spheres and hard hyperspheres in $`d=4`$ and $`d=5`$ with those obtained through Monte Carlo integration. Along the way we will introduce slight modifications to our original proposal and consider an alternative EOS to cope with mixtures very disparate in size. This new EOS is consistent with the forms suggested by Wheatley for the virial coefficients of hard discs  and hard spheres .

The paper is organized as follows. In the next Section and in order to make the paper self-contained, we will briefly recall the main results of our previous work and write down explicitly the composition-independent virial coefficients of a $`d`$-dimensional binary mixture of (additive) hard hyperspheres. Section III deals with a new proposal for the EOS of the binary mixture that attempts to incorporate additional features not present in the original recipe. In Section IV we focus on the case of hard spheres because of its intrinsic interest. We close the paper in Section V with an assessment of the different results. Some mathematical details of our derivations are presented in two Appendices.

## II Approximate equation of state and virial coefficients for a binary mixture of $`d`$-dimensional hard spheres

In this Section we first consider a multicomponent mixture of $`d`$-dimensional (additive) hard hyperspheres and then restrict ourselves to binary mixtures of hard hyperspheres in $`d`$ dimensions.

Let the number of components be $`N`$, the total number density of the mixture be $`\rho `$, the set of mole fractions be $`\{x_i\}`$, and the set of diameters be $`\{\sigma _i\}`$, $`(i=1,2,\mathrm{},N)`$. The packing fraction is $`\eta =_{i=1}^N\eta _i=v_d\rho \sigma ^d`$, where $`\eta _i=`$ $`v_d\rho _i\sigma _i^d`$ is the partial packing fraction due to species $`i`$, $`\rho _i=\rho x_i`$ is the partial number density corresponding to species $`i`$, $`v_d=(\pi /4)^{d/2}/\mathrm{\Gamma }(1+d/2)`$ is the volume of a $`d`$-dimensional sphere of unit diameter and $`\sigma ^n_{i=1}^Nx_i\sigma _i^n`$. In previous work  we proposed a simple equation of state for the mixture, $`Z_\text{m}(\eta )`$, consistent with a given EOS for a single component system, $`Z_\text{s}(\eta )`$, at the same packing fraction $`\eta `$, namely

$$Z_\text{m}(\eta )=1+\left[Z_\text{s}(\eta )1\right]2^{1d}\mathrm{\Delta }_0+\frac{\eta }{1\eta }\left(1\mathrm{\Delta }_0+\frac{1}{2}\mathrm{\Delta }_1\right),$$
(2)
with

$$\mathrm{\Delta }_p\frac{\sigma ^{d+p1}}{\sigma ^d^2}\underset{m=p}{\overset{d1}{}}\frac{(d+p1)!}{m!(d+p1m)!}\sigma ^{mp+1}\sigma ^{dm}.$$
(3)
In the one-dimensional case, Eq. (2) yields the exact result $`Z_\text{m}(\eta )=Z_\text{s}(\eta )`$. Further, for binary mixtures in $`d=2`$, $`3`$, $`4`$ and $`5`$, it proved to be very satisfactory when a reasonably accurate $`Z_\text{s}(\eta )`$ was taken .

From the virial expansion $`Z_\text{s}(\eta )=1+_{n=2}^{\mathrm{}}b_n\eta ^{n1}`$, where $`b_n`$ are (reduced) virial coefficients, and according to Eqs. (1)–(3), the $`n`$-th virial coefficient of the mixture is given by

$$B_n=v_d^{n1}\sigma ^d^{n1}\left(2^{1d}\mathrm{\Delta }_0b_n+1\mathrm{\Delta }_0+\frac{1}{2}\mathrm{\Delta }_1\right),$$
(4)
and explicitly for $`n2`$

$$B_n=v_d^{n1}\sigma ^d^{n3}\underset{m=0}{\overset{d}{}}\frac{(d1)!}{m!(dm)!}\sigma ^{dm}\left[(2^{1d}b_n1)(dm)\sigma ^{m+1}\sigma ^{d1}+\frac{d}{2}\sigma ^m\sigma ^d\right].$$
(5)
Now we particularize to binary mixtures. In that case, it is useful to define composition-independent coefficients $`B_{n_1,n_2}`$ as

$$B_n=\underset{n_1=0}{\overset{n}{}}B_{n_1,nn_1}\frac{n!}{n_1!(nn_1)!}x_1^{n_1}x_2^{nn_1}.$$
(6)
Our objective is to use our model to get approximate expressions of $`B_{n_1,n_2}`$ as explicit functions of $`\sigma _1`$, $`\sigma _2`$, and $`d`$. In Appendix A it is shown that the result is

$`B_{n_1,n_2}`$ $`=`$ $`v_d^{n1}\sigma _1^{d(n1)}\alpha ^{d(n_21)}[C_{n_1,n_2}^{(1)}\alpha ^d+C_{n_1,n_2}^{(2)}\alpha ^{d1}+C_{n_1,n_2}^{(3)}(1+\alpha )^{d1}`$ (8)
$`+C_{n_2,n_1}^{(3)}\alpha (1+\alpha )^{d1}+C_{n_2,n_1}^{(2)}\alpha +C_{n_2,n_1}^{(1)}],`$
where $`nn_1+n_2`$, $`\alpha \sigma _2/\sigma _1`$ and

$$C_{n_1,n_2}^{(1)}=\frac{n_1(n_11)\left[(n_12)b_n+2^{d1}n_2\right]}{n(n1)(n2)},$$
(9)
$$C_{n_1,n_2}^{(2)}=\frac{n_1n_2(n_11)\left(b_n2^{d1}\right)}{n(n1)(n2)},$$
(10)
$$C_{n_1,n_2}^{(3)}=\frac{n_1n_2\left[(n_21)2^{2d}b_n+n_1n_2\right]}{n(n1)(n2)}.$$
(11)
The special case of $`n=2`$ is obtained by first setting $`b_n=2^{d1}`$ and then setting $`n=2`$. The result is

$$C_{n_1,n_2}^{(1)}=2^{d2}n_1(n_11),C_{n_1,n_2}^{(2)}=0,C_{n_1,n_2}^{(3)}=\frac{1}{2}n_1n_2.$$
(12)
This gives the exact second virial coefficient.

Equation (8) expresses the reduced coefficient

$$B_{n_1,n_2}^{}(\alpha )B_{n_1,n_2}(\sigma _1,\sigma _2)\sigma _1^{d(n1)}\alpha ^{d(n_21)}$$
(13)
as a polynomial in $`\alpha `$ of degree $`d`$. Although it is known that the exact coefficients do not have in general such a polynomial structure , most of the proposals to date for $`d=2`$ and $`d=3`$ are polynomials . Let us see which consistency conditions the approximation (8)–(11) fulfills. First, it verifies the obvious symmetry property $`B_{n_1,n_2}(\sigma _1,\sigma _2)=B_{n_2,n_1}(\sigma _2,\sigma _1)`$, which implies that $`B_{n_1,n_2}^{}(\alpha )=\alpha ^dB_{n_2,n_1}^{}(1/\alpha )`$. In addition, the coefficients (9)–(11) satisfy the properties

$$C_{n_1,n_2}^{(1)}+C_{n_2,n_1}^{(2)}+2^{d1}C_{n_2,n_1}^{(3)}=\frac{n_1}{n}b_n,$$
(14)
$$dC_{n_1,n_2}^{(1)}+(d2)C_{n_1,n_2}^{(2)}+2^{d1}C_{n_2,n_1}^{(3)}=\frac{n_1}{n(n1)}\left[2^{d1}n_2+d(n_11)b_n\right],$$
(15)
$$C_{n,0}^{(1)}=b_n,C_{0,n}^{(1)}=C_{n,0}^{(2)}=C_{0,n}^{(2)}=C_{n,0}^{(3)}=C_{0,n}^{(3)}=0.$$
(16)
Equation (14) implies that

$$C_{n_1,n_2}^{(1)}+C_{n_2,n_1}^{(1)}+C_{n_1,n_2}^{(2)}+C_{n_2,n_1}^{(2)}+2^{d1}\left[C_{n_1,n_2}^{(3)}+C_{n_2,n_1}^{(3)}\right]=b_n,$$
(17)
which guarantees that

$$B_{n_1,n_2}^{}(\alpha =1)=v_d^{n1}b_n,$$
(18)
i.e., if both species have the same size we recover the one-component case. The same situation appears if $`x_1=0`$ or $`x_2=0`$, which means that $`B_{n,0}^{}(\alpha )=v_d^{n1}b_n\alpha ^d`$ and $`B_{0,n}^{}(\alpha )=v_d^{n1}b_n`$. This is verified as a consequence of Eq. (16). A subtler consistency condition is

$$\frac{B_{n_1,n_2}^{}}{\alpha }|_{\alpha =1}=d\frac{n_1}{n}v_d^{n1}b_n,$$
(19)
which ensures that the derivative of the excess free energy of mixing with respect to $`\sigma _2`$ is zero when the spheres are of the same size . In the approximation (8), Eq. (19) requires

$$dC_{n_1,n_2}^{(1)}+(d1)C_{n_1,n_2}^{(2)}+2^{d2}(d1)C_{n_1,n_2}^{(3)}+2^{d2}(d+1)C_{n_2,n_1}^{(3)}+C_{n_2,n_1}^{(2)}=d\frac{n_1}{n}b_n.$$
(20)
If we interchange the roles of $`n_1`$ and $`n_2`$ in Eq. (20) and add the original and the transformed equations, we get Eq. (17). If, instead, we subtract both equations, we obtain

$$d\left[C_{n_1,n_2}^{(1)}C_{n_2,n_1}^{(1)}\right]+(d2)\left[C_{n_1,n_2}^{(2)}C_{n_2,n_1}^{(2)}\right]2^{d1}\left[C_{n_1,n_2}^{(3)}C_{n_2,n_1}^{(3)}\right]=d\frac{n_1n_2}{n}b_n.$$
(21)
Thus, enforcement of Eq. (20) is equivalent to enforcement of Eq. (21). The property (15) guarantees that the condition (21) is fulfilled by the coefficients (9)–(11).

The most stringent conditions appear in the limits $`\alpha 0`$ and $`\alpha \mathrm{}`$. In that case, Eq. (8) yields

$$B_{n_1,n_2}^{}(\alpha )=\{\begin{array}{cc}v_d^{n1}\left[C_{n_1,n_2}^{(3)}+C_{n_2,n_1}^{(1)}\right]\hfill & (\alpha 0),\hfill \\ v_d^{n1}\alpha ^d\left[C_{n_1,n_2}^{(1)}+C_{n_2,n_1}^{(3)}\right]\hfill & (\alpha \mathrm{}).\hfill \end{array}$$
(22)
On the other hand, the exact result is

$$B_{n_1,n_2}^{}(\alpha )=\{\begin{array}{cc}v_d^{n1}\frac{n_2}{n}b_{n_2}\hfill & (\alpha 0),\hfill \\ v_d^{n1}\alpha ^d\frac{n_1}{n}b_{n_1}\hfill & (\alpha \mathrm{}).\hfill \end{array}$$
(23)
The above condition means that in the limit of infinite size asymmetry, the smaller spheres contribute to the total pressure as if they were a one-component system in a free volume equal to the total volume minus the volume occupied by the larger spheres. Our approximation (8)–(11) gives the correct forms for the asymptotic behaviours of $`B_{n_1,n_2}^{}`$, but not the correct coefficients (except for $`n_2=1`$ in the limit $`\alpha 0`$ and for $`n_1=1`$ in the limit $`\alpha \mathrm{}`$, apart from the trivial cases $`n_1=0,n`$ and $`n_2=0,n`$). In fact, in our approximation the virial coefficient $`B_n`$ depends parametrically on $`b_n`$ only, and not on the previous coefficients $`b_{n_1},n_1n`$, as is needed in Eq. (23). This is a consequence of the form of our approximation in which the compressibility factor of the mixture at a given packing fraction $`\eta `$ is expressed in terms of the compressibility factor of the monodisperse system at the same packing fraction.

Another consistency condition is

$$\frac{B_{n1,1}^{}}{\alpha }|_{\alpha =0}=v_d^{n1}d\frac{n1}{n}.$$
(24)
This comes from the exact value of the derivative of the excess free energy with respect to molecular size in the limit of infinite size asymmetry . In general, this condition is not verified by our approximation either. We consider the condition (23) to be more important than (24) because it is not restricted to $`n_2=1`$ (or $`n_1=1`$) and also because it refers to the leading term in the limit of infinite size asymmetry, while the condition (24) refers to the sub-leading term.

Table I gives the known values of the coefficients $`b_n`$ for $`2d5`$. From them, Eqs. (8)–(12) can be used to estimate the coefficients $`B_{n_1,n_2}^{}`$ for any value of $`\alpha `$ in these dimensions. We will come back to this point later on. The values of Table I can also be used to check to which degree Eq. (23) is violated by the approximation (8)–(11). The less favorable case corresponds to $`n_1=1`$ in the limit $`\alpha 0`$ (or, equivalently, $`n_2=1`$ in the limit $`\alpha \mathrm{}`$). Let us then define

$$R_n\frac{C_{1,n1}^{(3)}+C_{n1,1}^{(1)}}{[(n1)/n]b_{n1}}.$$
(25)
According to Eqs. (9) and (11), we have

$$R_n=\frac{\left(2^{2d}+n3\right)b_n+2^{d1}1}{(n1)b_{n1}}.$$
(26)
Table II shows the values of $`R_n`$ that can be obtained from the known virial coefficients. In the cases $`d=2`$ and $`d=3`$ the known values of the ratio $`R_n`$ differ from 1 less than 5%. In fact, by assuming $`R_n1`$, Eq. (26) can be used to estimate $`b_9`$ and $`b_{10}`$. This gives $`b_99.39`$ ($`d=2`$), $`b_986.7`$ ($`d=3`$), $`b_{10}10.44`$ ($`d=2`$) and $`b_{10}103.6`$ ($`d=3`$). Estimates based on Padé approximants are  $`b_99.37`$ ($`d=2`$), $`b_993.1`$ ($`d=3`$), $`b_{10}10.55`$ ($`d=2`$) and $`b_{10}123.2`$ ($`d=3`$). On the other hand, the deviations of $`R_n`$ from 1 become more important as the dimensionality increases and are already relatively large for $`d=5`$. Thus, at least for those high dimensionalities, it appears to be convenient to modify the prescription (8) as to satisfy the requirement (23). This is the subject of the next Section, where we will also discuss in more detail how good the estimates of $`B_{n_1,n_2}^{}`$, as obtained after using Eqs. (8)–(12) with the values of the coefficients $`b_n`$ given in Table I, are for any value of $`\alpha `$ and $`2d5`$.

## III A new proposal for the equation of state of a binary mixture

In order to account for the conditions imposed by Eq. (23) without greatly sacrificing the simplicity of our original recipe, let us now explore the possibility of keeping the structure (8), but with expressions for the coefficients different from (9)–(11). The requirement (23) implies that

$$C_{n_1,n_2}^{(3)}+C_{n_2,n_1}^{(1)}=\frac{n_2}{n}b_{n_2}.$$
(27)
In addition, the condition (18) yields Eq. (17). Note that Eq. (14) is more restrictive than Eq. (17). Since we want our modified expressions to remain as close as possible to the original ones, we will take (14) as a second condition. From Eqs. (27) and (14) we have

$$C_{n_1,n_2}^{(1)}=\left(2^{d1}1\right)^1\left[\frac{n_1}{n}\left(2^{d1}b_{n_1}b_n\right)+C_{n_2,n_1}^{(2)}\right],$$
(28)
$$C_{n_1,n_2}^{(3)}=\left(2^{d1}1\right)^1\left[\frac{n_2}{n}\left(b_nb_{n_2}\right)C_{n_1,n_2}^{(2)}\right].$$
(29)
We need a third condition to close the problem. Since $`C_{n_1,n_2}^{(2)}`$ does not enter into Eq. (22) and in agreement with our philosophy of departing from the original formulation as less as possible, one possibility is to keep Eq. (10); however, the resulting approximation is consistent with the conditions (16) but does not satisfy the requirement (20). If we enforce the fulfillment of Eq. (20) \[or, equivalently, of Eq. (21)\], then $`C_{n_1,n_2}^{(2)}`$ must have the following form

$$C_{n_1,n_2}^{(2)}=D_{n_1,n_2}\frac{n_1}{n}\frac{b_nb_{n_1}}{2^{2d}1},$$
(30)
where $`D_{n_1,n_2}`$ must be symmetric, but otherwise remains so far unknown. To close the problem, we impose the more restrictive property (15). This gives rise to

$`D_{n_1,n_2}`$ $`=`$ $`{\displaystyle \frac{1}{n\left[2+2^{d1}(d3)\right]}}[(2^{d1}1)(2^{d1}db_n){\displaystyle \frac{n_1n_2}{n1}}`$ (32)
$`+{\displaystyle \frac{2^{d1}d}{12^{2d}}}(nb_nn_1b_{n_1}n_2b_{n_2})].`$
Equations (28)–(32) close the modified version of our approximation.

It must be noted that Eqs. (30) and (32) are meaningless if $`d=2`$. This is related to the fact that if $`d=2`$ it is impossible to enforce conditions (18), (19) and (23) simultaneously on a polynomial of second degree for $`B_{n_1,n_2}^{}`$, as the number of conditions exceeds the number of unknowns. As a matter of fact, Eqs. (28) and (29) are enough to uniquely determine $`B_{n_1,n_2}^{}`$ in the case $`d=2`$, irrespective of the value of $`C_{n_1,n_2}^{(2)}`$. The result is

$$B_{n_1,n_2}^{}=v_2^{n1}\left[\frac{n_2}{n}b_{n_2}+\left(b_n\frac{n_1}{n}b_{n_1}\frac{n_2}{n}b_{n_2}\right)\alpha +\frac{n_1}{n}b_{n_1}\alpha ^2\right].$$
(33)
This was the form proposed by Wheatley  by imposing conditions (18) and (23). However, condition (19) is not verified unless $`n_1(b_nb_{n_1})=n_2(b_nb_{n_2})`$, which is only true if $`b_n=n`$, i.e., if one considers the scaled particle theory . For the sake of comparison, the form that stems from our original approximation, Eqs. (8)–(11), is (for $`d=2`$)

$$B_{n_1,n_2}^{}=v_2^{n1}\frac{n_1n_2}{n(n1)}\left[\frac{n_21}{n_1}b_n+1+2(b_n1)\alpha +\left(\frac{n_11}{n_2}b_n+1\right)\alpha ^2\right].$$
(34)
In summary, the modified expressions for the virial coefficients are given by Eq. (33) for $`d=2`$ and by Eqs. (8) and (28)–(32) for $`d3`$. We are now in a position to compare the theoretical predictions we have provided in this and the previous Section with the available Monte Carlo data for the virial coefficients. We note that the coefficients $`B_{n_1,n_2}`$ corresponding to $`n5`$ have been recently evaluated for $`d=2`$, $`d=3`$ , $`d=4`$ and $`d=5`$ . Further, the coefficients corresponding to $`n=6`$ for $`d=2`$ are also known . Because of its special physical interest and for reasons that will become apparent later, the case of hard spheres ($`d=3`$) will be addressed separately in Section IV.

We begin our assessment of our estimates of $`B_{n_1,n_2}`$ with the binary mixture of hard discs. This implies using Eq. (34) for the original recipe or using Eq. (33) in the modified version, which coincides with Wheatley’s proposal . In Figs. 1–4 we compare the results of taking these two routes with the ones of Refs. . As clearly seen in the figures, the theoretical prescriptions do a very good job for the whole range of $`\alpha `$, the modified version (33) being almost perfect.

In the case of hard hyperspheres in $`d=4`$ and $`d=5`$, Figs. 510 illustrate the performance of Eq. (8) together with (9)–(12) and of Eqs. (8) and (28)–(32) with respect to the “exact” data of Ref. . In these dimensionalities, it is clear that while the modified version of the estimates for the virial coefficients continues to exhibit a very good performance, the use of Eqs. (8)–(12) presents limitations for $`n=4`$ and especially for $`n=5`$. Nevertheless, the EOS for the binary mixture derived from Eq. (2) (which yields the original estimates of the coefficients $`B_{n_1,n_2}`$), once the proper $`Z_\text{s}(\eta )`$ is used in both $`d=4`$ and $`d=5`$, turns out to be very accurate.

Once we have analyzed the behaviour of the virial coefficients, the important question that now arises is, which is the EOS consistent with the modified virial coefficients? In the case $`d=2`$, from Eq. (33) one gets

$`Z_\text{m}(\eta )`$ $`=`$ $`{\displaystyle \frac{x_1}{1\eta _2}}Z_\text{s}\left({\displaystyle \frac{\eta _1}{1\eta _2}}\right){\displaystyle \frac{\sigma _2\sigma _1}{\sigma _2}}+{\displaystyle \frac{x_2}{1\eta _1}}Z_\text{s}\left({\displaystyle \frac{\eta _2}{1\eta _1}}\right){\displaystyle \frac{\sigma _1\sigma _2}{\sigma _1}}`$ (36)
$`+Z_\text{s}(\eta ){\displaystyle \frac{\sigma ^2}{\sigma _1\sigma _2}}.`$
The EOS consistent with the coefficients (28)–(32) can also be explicitly obtained for $`d3`$ by carrying out the resummation of the corresponding virial series. The result is (see Appendix B for details)

$`Z_\text{m}(\eta )`$ $`=`$ $`1+{\displaystyle \frac{K}{\sigma ^d}}[Z_\text{s}(\eta )1]+{\displaystyle \frac{K_0}{\sigma ^d}}x_1x_2{\displaystyle \frac{\eta }{1\eta }}`$ (38)
$`+{\displaystyle \frac{K_1}{\sigma _2^d}}x_1\left[{\displaystyle \frac{1}{1\eta _2}}Z_\text{s}\left({\displaystyle \frac{\eta _1}{1\eta _2}}\right)1\right]+{\displaystyle \frac{K_2}{\sigma _1^d}}x_2\left[{\displaystyle \frac{1}{1\eta _1}}Z_\text{s}\left({\displaystyle \frac{\eta _2}{1\eta _1}}\right)1\right],`$
where $`K`$, $`K_0`$, $`K_1`$ and $`K_2`$ are independent of the density and are given by

$$K=x_1^2\sigma _1^d\left(1\frac{K_1}{\sigma _2^d}\right)+x_2^2\sigma _2^d\left(1\frac{K_2}{\sigma _1^d}\right)+2x_1x_2\left(\sigma _{12}^d\frac{K_1+K_2+K_0}{2^d}\right),$$
(39)
$$K_0=\frac{1}{2^{2d}+d3}\left[\sigma _1^d+\sigma _2^d+\left(2^{d1}1\right)\sigma _1\sigma _2\left(\sigma _1^{d2}+\sigma _2^{d2}\right)2^d\sigma _{12}^d\right],$$
(40)
$`K_1`$ $`=`$ $`{\displaystyle \frac{1}{(2^d4)(2^{2d}+d3)\sigma _{12}}}\{\sigma _{12}^d[2^d(d2)\sigma _1+8(12^{d1}+d2^{d3})\sigma _2]`$ (43)
$`(2^{d1}2)\sigma _1\sigma _2\left(\sigma _1^{d1}+\sigma _2^{d1}\right)(d2)(2^{d1}1)\sigma _1^2\sigma _2^{d1}(2^{d1}d)\sigma _1^{d1}\sigma _2^2`$
$`(d2)\sigma _1^{d+1}[d42^{d1}(d3)]\sigma _2^{d+1}\},`$
$`K_2`$ $`=`$ $`{\displaystyle \frac{1}{(2^d4)(2^{2d}+d3)\sigma _{12}}}\{\sigma _{12}^d[2^d(d2)\sigma _28(12^{d1}+d2^{d3})\sigma _1]`$ (46)
$`(2^{d1}2)\sigma _1\sigma _2\left(\sigma _1^{d1}+\sigma _2^{d1}\right)(d2)(2^{d1}1)\sigma _2^2\sigma _1^{d1}(2^{d1}d)\sigma _2^{d1}\sigma _1^2`$
$`(d2)\sigma _2^{d+1}[d42^{d1}(d3)]\sigma _1^{d+1}\}.`$
Here, $`\sigma _{12}(\sigma _1+\sigma _2)/2`$. Note that Eq. (36) for hard discs is included in the structure of Eq. (38), except that then Eqs. (39)–(46) cannot be applied and we have instead $`K_0=0`$, $`K_1=\sigma _2(\sigma _2\sigma _1)`$, $`K_2=\sigma _1(\sigma _1\sigma _2)`$ and $`K=\sigma ^2\left(1x_1K_1/\sigma _2^2x_2K_2/\sigma _1^2\right)`$.

Equation (38) complies with the exact results given in Eq. (23). It is clearly more complicated than the original recipe, Eq. (2), and should in principle be useful in particular for mixtures involving components of disparate sizes. In contrast to Eq. (2), the proposal (38) expresses $`Z_\text{m}(\eta )`$ in terms not only of $`Z_\text{s}(\eta )`$ but also involves $`Z_\text{s}\left(\frac{\eta _1}{1\eta _2}\right)`$ and $`Z_\text{s}\left(\frac{\eta _2}{1\eta _1}\right)`$.

We note that if $`d=3`$, it follows from Eqs. (39)–(46) that $`K_0=0`$, $`K_1=\sigma _2(\sigma _1\sigma _2)^2`$, $`K_2=\sigma _1(\sigma _1\sigma _2)^2`$ and $`K=\sigma ^3\left(1x_1K_1/\sigma _2^3x_2K_2/\sigma _1^3\right)`$. Consequently, Eq. (38) for hard spheres becomes

$`Z_\text{m}(\eta )`$ $`=`$ $`{\displaystyle \frac{x_1}{1\eta _2}}Z_\text{s}\left({\displaystyle \frac{\eta _1}{1\eta _2}}\right){\displaystyle \frac{(\sigma _1\sigma _2)^2}{\sigma _2^2}}+{\displaystyle \frac{x_2}{1\eta _1}}Z_\text{s}\left({\displaystyle \frac{\eta _2}{1\eta _1}}\right){\displaystyle \frac{(\sigma _1\sigma _2)^2}{\sigma _1^2}}`$ (48)
$`+Z_\text{s}(\eta ){\displaystyle \frac{2\sigma ^2\sigma _1\sigma _2\sigma ^4}{\sigma _1^2\sigma _2^2}}.`$
It is interesting to note that Eqs. (36) and (48) can be written in the common form

$`Z_\text{m}(\eta )`$ $`=`$ $`Z_\text{s}(\eta )+x_1\left[{\displaystyle \frac{1}{1\eta _2}}Z_\text{s}\left({\displaystyle \frac{\eta _1}{1\eta _2}}\right)Z_\text{s}(\eta )\right]\left({\displaystyle \frac{\sigma _2\sigma _1}{\sigma _2}}\right)^{d1}`$ (50)
$`+x_2\left[{\displaystyle \frac{1}{1\eta _1}}Z_\text{s}\left({\displaystyle \frac{\eta _2}{1\eta _1}}\right)Z_\text{s}(\eta )\right]\left({\displaystyle \frac{\sigma _1\sigma _2}{\sigma _1}}\right)^{d1}.`$
Equation (50) reduces to Eq. (36) if $`d=2`$ and to Eq. (48) if $`d=3`$. In addition, it becomes exact for $`d=1`$, i.e., $`Z_\text{s}(\eta )=Z_\text{m}(\eta )=(1\eta )^1`$. It seems then very tempting to assume Eq. (50) for arbitrary $`d`$. With such an assumption, after using Eqs. (B7), (B8), (B12) and (B13), we simply get for the composition-independent virial coefficients

$$B_{n_1,n_2}^{}=v_d^{n1}\left[b_n\frac{n_1\alpha ^d+n_2}{n}+\left(b_{n_1}b_n\right)\frac{n_1}{n}\alpha (\alpha 1)^{d1}+\left(b_{n_2}b_n\right)\frac{n_2}{n}(1\alpha )^{d1}\right].$$
(51)
Of course, this includes the cases $`d=2`$ \[cf. Eq. (33)\] and $`d=3`$ \[cf. Eq. (60)\]. However, Eq. (51) differs from the prescription (28)–(32) if $`d4`$. Although Eq. (51) is consistent with conditions (18), (19) \[for $`d3`$\] and (23), it does not reproduce the exact second virial coefficient $`B_{1,1}^{}(\alpha )=v_d(1+\alpha )^d/2`$ (except, of course, for $`d3`$). In fact, we have checked that the performance of Eqs. (50) and (51) is very poor for $`d4`$, so that they will not be further considered here.

## IV The equation of state of binary hard-sphere mixtures

Up to here, although ocasionally we have mentioned explicit features of hard discs and hard spheres, the developments presented in the previous Sections apply for general dimensionality $`d`$. Given their intrinsic importance, in this Section we focus on the results for the EOS in three dimensions. In order to set the framework of our discussion, we shall first recall two other EOS in the literature which share with ours the idea of obtaining $`Z_\text{m}(\eta )`$ from the knowledge of $`Z_\text{s}(\eta )`$.

The first one is due to Hamad who in 1994  derived a consistency condition involving the derivatives of the contact values of the radial distribution functions $`g_{ij}`$ with respect to $`\sigma _k`$ . Imposing this condition, he proposed the following equation of state for hard-sphere mixtures

$$Z_\text{m}^\text{H}(\eta )=Z_\text{s}(\eta )+\frac{3\eta }{(1\eta )^3}\left[\frac{\sigma ^2^3}{\sigma ^3^2}\eta +\frac{\sigma \sigma ^2}{\sigma ^3}(1\eta )1\right],$$
(52)
where we have used the label H to refer to Hamad’s proposal. It must be noted that Eq. (52) was not explicitly written in Ref. , where the results are expressed in a rather more involved form. The second EOS for binary mixtures of hard spheres that we want to consider was proposed more recently by Barrio and Solana (label BS). It reads

$$Z_\text{m}^{\text{BS}}(\eta )=1+\frac{1}{4}(1+\beta \eta )\left(1+3\frac{\sigma \sigma ^2}{\sigma ^3}\right)\left[Z_\text{s}(\eta )1\right],$$
(53)
where $`\beta `$ is adjusted as to reproduce the exact third virial coefficient , namely

$$\begin{array}{c}B_{3,0}^{}=v_3^210\alpha ^3,B_{2,1}^{}=v_3^2\left(\frac{1}{3}+2\alpha +5\alpha ^2+\frac{8}{3}\alpha ^3\right),\hfill \\ B_{1,2}^{}=v_3^2\left(\frac{8}{3}+5\alpha +2\alpha ^2+\frac{1}{3}\alpha ^3\right),B_{0,3}^{}=v_3^210.\hfill \end{array}$$
(54)
More explicitly, in our notation $`\beta `$ is written as

$$\beta =\frac{B_3}{v_3^2\sigma ^3^2(1+3\sigma \sigma ^2/\sigma ^3)}\frac{10}{4}.$$
(55)
We now first concentrate on the composition-independent virial coefficients. Our original recipe, Eqs. (8)–(11), gives for $`d=3`$

$`B_{n_1,n_2}^{\text{SYH}}`$ $`=`$ $`v_3^{n1}{\displaystyle \frac{n_1n_2}{2n(n1)(n2)}}\{2(n_1+3n_24)+{\displaystyle \frac{n_21}{n_1}}(n_1+2n_24)b_n`$ (58)
$`+\left[2(n_15n_2+4)+(n_1+4n_25)b_n\right]\alpha +\left[2(n_25n_1+4)+(n_2+4n_15)b_n\right]\alpha ^2`$
$`+[2(n_2+3n_14)+{\displaystyle \frac{n_11}{n_2}}(n_2+2n_14)b_n]\alpha ^3\},`$
where we have incorporated the label SYH to distinguish it from the other results. On the other hand, in the three-dimensional case, Eqs. (28)–(32) yield

$`B_{n_1,n_2}^\text{W}`$ $`=`$ $`v_3^{n1}[{\displaystyle \frac{n_2}{n}}b_{n_2}+({\displaystyle \frac{2n_2n_1}{n}}b_n+{\displaystyle \frac{n_1}{n}}b_{n_1}2{\displaystyle \frac{n_2}{n}}b_{n_2})\alpha `$ (60)
$`+({\displaystyle \frac{2n_1n_2}{n}}b_n+{\displaystyle \frac{n_2}{n}}b_{n_2}2{\displaystyle \frac{n_1}{n}}b_{n_1})\alpha ^2+{\displaystyle \frac{n_1}{n}}b_{n_1}\alpha ^3],`$
which coincides with the form proposed by Wheatley  and so the label in this case is W. The virial coefficients associated with Eqs. (52) and (53) are also readily derived. In the former case, first it is easy to obtain

$$B_n^\text{H}=v_3^{n1}\sigma ^3^{n1}\left[b_n\frac{3}{2}n(n1)+3(n1)\frac{\sigma \sigma ^2}{\sigma ^3}+\frac{3}{2}(n1)(n2)\frac{\sigma ^2^3}{\sigma ^3^2}\right].$$
(61)
Now, by making $`m=0`$ and $`m=1`$ in Eq. (A9) and $`m=1`$ in Eq. (A6), and inserting the results into Eq. (61), we get the composition-independent coefficients arising from Hamad’s proposal,

$`B_{n_1,n_2}^\text{H}`$ $`=`$ $`v_3^{n1}\{{\displaystyle \frac{n_2}{n}}[b_n{\displaystyle \frac{3}{2}}n_1(n1+n_2)]+{\displaystyle \frac{3}{2}}{\displaystyle \frac{n_1n_2}{n}}(3n_21)\alpha `$ (63)
$`+{\displaystyle \frac{3}{2}}{\displaystyle \frac{n_1n_2}{n}}(3n_11)\alpha ^2+{\displaystyle \frac{n_1}{n}}[b_n{\displaystyle \frac{3}{2}}n_2(n1+n_1)]\alpha ^3\}.`$
As far as Eq. (53) is concerned, one readily derives (for $`n3`$)

$`B_n^{\text{BS}}`$ $`=`$ $`{\displaystyle \frac{v_3^{n1}}{4}}(b_n+\beta b_{n1})\sigma ^3^{n1}\left(1+3{\displaystyle \frac{\sigma \sigma ^2}{\sigma ^3}}\right)`$ (64)
$`=`$ $`{\displaystyle \frac{v_3^{n1}}{4}}\left(b_n{\displaystyle \frac{10}{4}}b_{n1}\right)\sigma ^3^{n1}\left(1+3{\displaystyle \frac{\sigma \sigma ^2}{\sigma ^3}}\right)+{\displaystyle \frac{v_3^{n3}}{4}}b_{n1}B_3\sigma ^3^{n3}.`$ (65)
Now, after some algebra also involving the substitution of the results obtained by taking $`m=0`$ and $`m=1`$ in Eq. (A9) into Eq. (65), the composition-independent virial coefficients for the Barrio-Solana EOS turn out to be given by

$`B_{n_1,n_2}^{\text{BS}}`$ $`=`$ $`v_3^{n1}{\displaystyle \frac{b_n}{4n(n1)}}\left[n_2(3n_2+n4)+3n_1n_2\alpha +3n_1n_2\alpha ^2+n_1(3n_1+n4)\alpha ^3\right]`$ (68)
$`v_3^{n1}{\displaystyle \frac{3b_{n1}n_1n_2}{8n(n1)(n2)}}[n+2n_24+(n6n_2+4)\alpha `$
$`+(n6n_1+4)\alpha ^2+(n+2n_14)\alpha ^3].`$
These four approximations for $`B_{n_1,n_2}^{}`$ satisfy the requirements (18) and (19). The condition (23) is only satisfied by the W approximation, while the condition (24) is only satisfied by the H approximation. By construction, the BS approximation, Eq. (68), reproduces the exact third virial coefficients (54). It is interesting to point out the approximations (58), (60) and (63) also yield the exact value of the third virial coefficients.

In order to assess the accuracy of the previous approximate expressions for $`B_{n_1,n_2}^{}`$ with $`n4`$, we will take the case $`n_1=3`$, $`n_2=1`$ as a benchmark for comparison since the coefficient $`B_{3,1}^{}`$ has been obtained exactly for the interval $`0\alpha 2/\sqrt{3}10.1547`$ . This exact result is

$`B_{3,1}^{}`$ $`=`$ $`{\displaystyle \frac{v_3^3}{4}}(1+9\alpha +36\alpha ^2+21\alpha ^3+{\displaystyle \frac{27}{2}}\alpha ^4+{\displaystyle \frac{27}{10}}\alpha ^5`$ (70)
$`{\displaystyle \frac{108}{5}}\alpha ^6{\displaystyle \frac{648}{35}}\alpha ^7{\displaystyle \frac{81}{14}}\alpha ^8{\displaystyle \frac{9}{14}}\alpha ^9).`$
Setting $`n_1=3`$ and $`n_2=1`$ in Eqs. (58), (60), (63) and (68), respectively, yields

$$B_{3,1}^{\text{SYH}}=\frac{v_3^3}{4}\left[1+\left(\frac{b_4}{2}+1\right)\alpha +\left(2b_45\right)\alpha ^2+3\left(\frac{b_4}{2}+1\right)\alpha ^3\right],$$
(71)
$$B_{3,1}^\text{W}=\frac{v_3^3}{4}\left[1+\left(28b_4\right)\alpha +\left(5b_459\right)\alpha ^2+30\alpha ^3\right],$$
(72)
$$B_{3,1}^\text{H}=\frac{v_3^3}{4}\left[b_418+9\alpha +36\alpha ^2+3(b_49)\alpha ^3\right]$$
(73)
and

$$B_{3,1}^{\text{BS}}=\frac{v_3^3}{16}\left[b_415+(3b_415)\alpha +(3b_4+75)\alpha ^2+(9b_445)\alpha ^3\right].$$
(74)
Note that Eqs. (71) and (72) would coincide if $`b_4`$ were equal to $`18`$, as happens with the CS equation of state. However, they differ if we take the exact value $`b_4=18.36477`$. The maximum deviation of (72) from the exact result (70) is about $`1.8\%`$ (at $`\alpha 0.08`$), while that of (71) is about $`3.7\%`$ (at $`\alpha 0.09`$). Also, if $`b_4`$ were equal to $`19`$, Eq. (73) would coincide with Eq. (72), while using the exact value for $`b_4`$ Eq. (73) implies that $`lim_{\alpha 0}B_{3,1}^\text{H}=0.091v_3^3`$ instead of the exact value $`lim_{\alpha 0}B_{3,1}^{}=\frac{1}{4}v_3^3`$, which, however, is satisfied by the approximations (71) and (72). Comparison between Hamad’s approximation (73) and the exact result (70) for $`\alpha 2/\sqrt{3}1`$ shows that the relative error of (73) monotonically decreases from about $`64\%`$ at $`\alpha =0`$ to about $`18\%`$ at $`\alpha =2/\sqrt{3}1`$. Finally, according to Eq. (74) with the exact value of $`b_4`$, $`lim_{\alpha 0}B_{3,1}^{\text{BS}}=0.210v_3^3`$, while in the interval $`\alpha 2/\sqrt{3}1`$ the relative error of (74) monotonically decreases from about $`16\%`$ at $`\alpha =0`$ to about $`1.7\%`$ at $`\alpha =2/\sqrt{3}1`$. This relative error is always larger than that of Eqs. (71) and (72), but smaller than that of Eq. (73).

Apart from the exact result for $`B_{3,1}^{}`$, as pointed out before, numerical values of the coefficients $`B_{n_1,n_2}^{}`$ corresponding to $`n5`$ are also available for hard spheres. In Figs. 11 and 12 we display the comparison between the theoretical estimates we have just discussed for the virial coefficients and the numerical results. It is clear that the estimates $`B_{n_1,n_2}^{\text{SYH}}`$, $`B_{n_1,n_2}^\text{W}`$ and $`B_{n_1,n_2}^{\text{BS}}`$ for $`n=4`$ and $`n=5`$ are remarkably close to one another and to the exact results, except perhaps in the case of $`B_{1,4}^{\text{BS}}`$. On the other hand, the prediction of $`B_{n_1,n_2}^\text{H}`$ is much poorer.

Once we have briefly discussed the performance of the various approximate expressions for the composition-independent virial coefficients, we may next wonder how well the approximate EOS for a binary mixture of hard spheres perform with respect to the simulation data. To address this issue, we have first to specify what $`Z_\text{s}(\eta )`$ is. Two choices will be made here. On the one hand, since perhaps the most popular and widely known EOS for a single component fluid of hard spheres is the CS EOS, we will also use it here. On the other hand, a more accurate EOS for the same system is obtained if a rescaled Padé approximation taking into account the known virial coefficients is used. In Tables III and IV we present a comparison of the results obtained using the various approximate EOS for the mixture with these two choices of $`Z_\text{s}(\eta )`$, two different mole fractions and a diameter ratio $`\alpha =0.3`$. In Table III, the values corresponding to the widely used Boublík-Mansoori-Carnahan-Starling-Leland (BMCSL) EOS  have also been included. One immediately notes that, as already stated in the literature, the EOS given in Eq. (2) using the CS EOS for the single component system provides a very accurate account of the simulation data, better than any of the other EOS for mixtures considered, including that of Eq. (48). If instead of the CS EOS the rescaled Padé approximation is taken for $`Z_\text{s}(\eta )`$, then it is Eq. (48) which gives in general the best results, although the original recipe (2) and the BS EOS (53) are not too far behind, the latter being especially accurate in the case $`x_1=0.75`$. It should be noted that, irrespective of the choice for $`Z_\text{s}(\eta )`$, both $`Z_\text{m}^\text{H}(\eta )`$ and $`Z_\text{m}^{\text{BS}}(\eta )`$ for these size ratio and packing fractions always underestimate the simulation results, a comment that also applies to the BMCSL EOS.

## V  Concluding remarks

In this paper we have further assessed the merits and limitations of a simple recipe recently introduced to derive the compressibility factor of a multicomponent mixture of $`d`$-dimensional hard spheres. We have now considered the case of binary mixtures and looked in particular at the composition-independent virial coefficients. The comparison with the available exact and simulation results confirmed the usefulness of our approach except perhaps when the mixture involves components of very disparate sizes (especially for high dimensionalities). Guided by the conditions that the virial coefficients must fulfill in certain limits, a slight modification of the form of the original coefficients was made, trying to sacrifice simplicity as less as possible while at the same time improving the accuracy of their numerical predictions. It is fair to state that the (modified) composition-independent virial coefficients are in excellent agreement with the reported values for all dimensions and all size ratios.

Due to its particular relevance, the case of the binary mixture of hard spheres $`(d=3)`$ was analyzed in some detail, both in connection with the virial coefficients and through the consideration of different EOS for this system that exist in the literature. As for the virial coefficients, the trends pointed out above about the reliability of the original and modified prescriptions hold also for $`d=3`$. The comparison with the simulation results for the compressibility factor of these mixtures given in Tables III and IV indicates the superiority of the original recipe for the EOS provided one uses the CS EOS for $`Z_\text{s}(\eta )`$. On the other hand, if the more accurate rescaled Padé approximation for $`Z_\text{s}(\eta )`$ is used, it is the new EOS (48) (i.e., the one obtained after resumming the virial series with the modified coefficients) which gives the best performance, although only slightly better than $`Z_\text{m}(\eta )`$ as given by Eq. (2) with $`d=3`$. Both equations are also of comparable accuracy to the one of the presumably best EOS available presently for hard-sphere mixtures, namely the rescaled Padé approximant for mixtures introduced by Malijevský and Veverka . This latter equation of state rests on a different philosophy, however, in the sense that it makes use of the known virial coefficients of the mixture rather than the $`Z_s(\eta )`$, as in our approach and those of Hamad  and Barrio and Solana .

In conclusion, we want to add that both the original proposal of the EOS of a binary mixture of $`d`$-dimensional hard spheres, Eq. (2), and the modified version, Eq. (38), provide two simple and yet accurate EOS at least with respect to the composition-independent virial coefficients and in comparison with the (relatively scarce) simulation data. Further confirmation of their usefulness depends on new simulations which we hope will be encouraged by the results presented in this paper.

###### Acknowledgements.

A.S. and S.B.Y. acknowledge partial support from the DGES (Spain) through grant No. PB97-1501 and from the Junta de Extremadura (Fondo Social Europeo) through grant No. IPR99C031. They are also grateful to the DGES (Spain) for sabbatical grants No. PR2000-0117 and No. PR2000-0116, respectively. The research of M.L.H. was supported in part by DGAPA-UNAM under Project IN103100.

## A Derivation of Eqs. (8)–(11)

In order to identify the composition-independent coefficients $`B_{n_1,n_2}`$ from Eq. (5), we need to expand the moments $`\sigma ^p`$. After some algebra, we get

$`\sigma ^d^{n3}\sigma ^{dm}\sigma ^{m+1}\sigma ^{d1}`$ $`=`$ $`\sigma _1^{d(n1)}{\displaystyle \underset{n_1=0}{\overset{n}{}}}\alpha ^{d(nn_11)}{\displaystyle \frac{(n3)!}{n_1!(nn_1)!}}x_1^{n_1}x_2^{nn_1}`$ (A6)
$`\times [n_1(n_11)(n_12)\alpha ^d+n_1(n_11)(nn_1)\alpha ^{d1}`$
$`+n_1(n_11)(nn_1)\alpha ^{dm}+n_1(nn_1)(nn_11)\alpha ^{dm1}`$
$`+n_1(n_11)(nn_1)\alpha ^{m+1}+n_1(nn_1)(nn_11)\alpha ^m`$
$`+n_1(nn_1)(nn_11)\alpha `$
$`+(nn_1)(nn_11)(nn_12)],`$
$`\sigma ^d^{n2}\sigma ^{dm}\sigma ^m`$ $`=`$ $`\sigma _1^{d(n1)}{\displaystyle \underset{n_1=0}{\overset{n}{}}}\alpha ^{d(nn_11)}{\displaystyle \frac{(n2)!}{n_1!(nn_1)!}}x_1^{n_1}x_2^{nn_1}`$ (A9)
$`\times [n_1(n_11)\alpha ^d+n_1(nn_1)\alpha ^{dm}`$
$`+n_1(nn_1)\alpha ^m+(nn_1)(nn_11)],`$
where $`\alpha \sigma _2/\sigma _1`$. Now we use the identities

$$\underset{m=0}{\overset{d}{}}\frac{d!}{m!(dm)!}\alpha ^m=(1+\alpha )^d,$$
(A10)
$$\underset{m=0}{\overset{d}{}}\frac{(d1)!}{m!(dm)!}(dm)\alpha ^m=(1+\alpha )^{d1},$$
(A11)
to get

$`{\displaystyle \underset{m=0}{\overset{d}{}}}{\displaystyle \frac{(d1)!}{m!(dm)!}}(dm)\sigma ^d^{n3}\sigma ^{dm}\sigma ^{m+1}\sigma ^{d1}=\sigma _1^{d(n1)}{\displaystyle \underset{n_1=0}{\overset{n}{}}}\alpha ^{d(nn_11)}{\displaystyle \frac{(n3)!}{n_1!(nn_1)!}}`$ (A12)
$`\times x_1^{n_1}x_2^{nn_1}[n_1(n_11)(n_12)2^{d1}\alpha ^d+n_1(n_11)(nn_1)2^{d1}\alpha ^{d1}`$ (A13)
$`+2n_1(n_11)(nn_1)\alpha (1+\alpha )^{d1}+2n_1(nn_1)(nn_11)(1+\alpha )^{d1}`$ (A14)
$`+n_1(nn_1)(nn_11)2^{d1}\alpha +(nn_1)(nn_11)(nn_12)2^{d1}],`$ (A15)
$`{\displaystyle \underset{m=0}{\overset{d}{}}}{\displaystyle \frac{d!}{m!(dm)!}}\sigma ^d^{n2}\sigma ^{dm}\sigma ^m=\sigma _1^{d(n1)}{\displaystyle \underset{n_1=0}{\overset{n}{}}}\alpha ^{d(nn_11)}{\displaystyle \frac{(n2)!}{n_1!(nn_1)!}}x_1^{n_1}x_2^{nn_1}`$ (A16)
$`\times \left[n_1(n_11)2^d\alpha ^d+2n_1(nn_1)(1+\alpha )^d+(nn_1)(nn_11)2^d\right].`$ (A17)
From here we finally obtain Eqs. (8)–(11).

## B Derivation of Eqs. (38)–(46)

Let us rewrite Eq. (8) as

$`B_{n_1,n_2}`$ $`=`$ $`v_d^{n1}[C_{n_1,n_2}^{(1)}\sigma _1^{d(n_11)}\sigma _2^{dn_2}+C_{n_1,n_2}^{(2)}\sigma _1^{d(n_11)+1}\sigma _2^{dn_21}+C_{n_1,n_2}^{(3)}2^{d1}\sigma _{12}^{d1}\sigma _1^{d(n_11)+1}\sigma _2^{d(n_21)}`$ (B2)
$`+12],`$
where $`\sigma _{12}(\sigma _1+\sigma _2)/2`$ and the notation $`12`$ indicates that a contribution similar to the three summands appearing in the brackets in which the roles of $`(n_1,n_2)`$ and of $`(\sigma _1,\sigma _2)`$ are interchanged (i.e., $`n_1n_2`$, $`\sigma _1\sigma _2`$ and vice versa) should be added. The compressibility factor of the mixture is then given by

$`Z_\text{m}(\eta )1`$ $`=`$ $`{\displaystyle \underset{n=2}{\overset{\mathrm{}}{}}}\rho ^{n1}{\displaystyle \underset{n_1=0}{\overset{n}{}}}B_{n_1,n_2}{\displaystyle \frac{n!}{n_1!n_2!}}x_1^{n_1}x_2^{n_2}`$ (B3)
$`=`$ $`{\displaystyle \underset{n=2}{\overset{\mathrm{}}{}}}(v_d\rho )^{n1}{\displaystyle \underset{n_1=0}{\overset{n}{}}}{\displaystyle \frac{n!}{n_1!n_2!}}x_1^{n_1}x_2^{n_2}[C_{n_1,n_2}^{(1)}\sigma _1^{d(n_11)}\sigma _2^{dn_2}+C_{n_1,n_2}^{(2)}\sigma _1^{d(n_11)+1}\sigma _2^{dn_21}`$ (B6)
$`+C_{n_1,n_2}^{(3)}2^{d1}\sigma _{12}^{d1}\sigma _1^{d(n_11)+1}\sigma _2^{d(n_21)}]`$
$`+12.`$
Now we insert Eqs. (28)–(32) and sum over $`n_1`$ and $`n`$. To do that, we use the identities

$$x_1\left[Z_\text{s}\left(\eta \right)1\right]=\underset{n=2}{\overset{\mathrm{}}{}}(v_d\rho )^{n1}\underset{n_1=0}{\overset{n}{}}\frac{n!}{n_1!n_2!}x_1^{n_1}x_2^{n_2}b_n\frac{n_1}{n}\sigma _1^{d(n_11)}\sigma _2^{dn_2},$$
(B7)
$$x_2\left[Z_\text{s}\left(\eta \right)1\right]=\underset{n=2}{\overset{\mathrm{}}{}}(v_d\rho )^{n1}\underset{n_1=0}{\overset{n}{}}\frac{n!}{n_1!n_2!}x_1^{n_1}x_2^{n_2}b_n\frac{n_2}{n}\sigma _1^{dn_1}\sigma _2^{d(n_21)},$$
(B8)
$$\sigma ^d\left[Z_\text{s}(\eta )1\right]=\underset{n=2}{\overset{\mathrm{}}{}}(v_d\rho )^{n1}\underset{n_1=0}{\overset{n}{}}\frac{n!}{n_1!n_2!}x_1^{n_1}x_2^{n_2}b_n\sigma _1^{dn_1}\sigma _2^{dn_2},$$
(B9)
$$x_1x_2\frac{Z_\text{s}(\eta )1}{\sigma ^d}=\underset{n=2}{\overset{\mathrm{}}{}}(v_d\rho )^{n1}\underset{n_1=0}{\overset{n}{}}\frac{n!}{n_1!n_2!}x_1^{n_1}x_2^{n_2}b_n\frac{n_1n_2}{n(n1)}\sigma _1^{d(n_11)}\sigma _2^{d(n_21)},$$
(B10)
$$x_1x_2\frac{\eta }{(1\eta )\sigma ^d}=\underset{n=2}{\overset{\mathrm{}}{}}(v_d\rho )^{n1}\underset{n_1=0}{\overset{n}{}}\frac{n!}{n_1!n_2!}x_1^{n_1}x_2^{n_2}\frac{n_1n_2}{n(n1)}\sigma _1^{d(n_11)}\sigma _2^{d(n_21)}.$$
(B11)
In addition,

$$\frac{x_1}{1\eta _2}Z_\text{s}\left(\frac{\eta _1}{1\eta _2}\right)x_1=\underset{n=2}{\overset{\mathrm{}}{}}(v_d\rho )^{n1}\underset{n_1=0}{\overset{n}{}}\frac{n!}{n_1!n_2!}x_1^{n_1}x_2^{n_2}b_{n_1}\frac{n_1}{n}\sigma _1^{d(n_11)}\sigma _2^{dn_2},$$
(B12)
$$\frac{x_2}{1\eta _1}Z_\text{s}\left(\frac{\eta _2}{1\eta _1}\right)x_2=\underset{n=2}{\overset{\mathrm{}}{}}(v_d\rho )^{n1}\underset{n_1=0}{\overset{n}{}}\frac{n!}{n_1!n_2!}x_1^{n_1}x_2^{n_2}b_{n_2}\frac{n_2}{n}\sigma _1^{dn_1}\sigma _2^{d(n_21)}.$$
(B13)
This allows one to obtain $`Z_\text{m}(\eta )`$ in terms of $`Z_\text{s}`$ evaluated at $`\eta `$, $`\eta _1/(1\eta _2)`$ and $`\eta _2/(1\eta _1)`$. After some algebra, we arrive at Eqs. (38)–(46).