2 The method

The method we use to study regularity of the distribution of clusters was proposed by Dr. J. Pelt and was first applied to the cluster distribution by O. Toomet (1997).

2.1 The idea

As hints of regularity are seen in the correlation function and in the power spectrum, the first idea would be to use these statistics, but in their full three-dimensional form. For a general homogeneous density field the correlation function $\xi({\vec d})$ is a function of a spatial displacement vector ${\vec d}$and the power spectrum $P({\vec k})$ depends on the three-dimensional wave vector ${\vec k}$. The advantage of these statistics is that they describe directly the periodicity in the data. The disadvantages are larger, however. Firstly, these statistics are extremely noisy for real data samples, because we need much more individual amplitudes to populate a three-dimensional region than a line segment, as we have accustomed to do for isotropic $\xi(d)$ and P(k). Secondly, even simple real space patterns give rise to a number of interconnected amplitudes in the ${\vec d}$ and ${\vec k}$ space; they practically transform one 3-D distribution to another that has to be analyzed again. So, it will be easier (and more transparent) to study the spatial distribution directly.

For this purpose we shall use similar methods as are employed to find periods in observed time series, say, the brightness of a variable star in different moments. The most intuitive method is to use a trial period to fold the series into a "phase diagram'', and to find the value of a statistic describing the reliability of that trial period. An illustration of that is shown in Fig. 1.

The data in the phase diagram (middle and right panels of Fig. 1) are usually binned, and mainly two methods are used to find the best period. The first one, called "phase dispersion minimization'' (PDM), proposed by Stellingwerf (1978), does that by minimizing the sum of variances of data in individual phase bins (finding the narrowest light curve). Another method, called "epoch folding'' (Leahy et al. 1983), maximizes the variance of the means of the phase bins (looking for a maximum amplitude of the light curve). Though surprising at the first glance, these two methods are were shown to be equivalent by Schwarzenberg-Czerny (1989) and Davies (1990).

Phase diagrams already have been applied to the study of regularity in the galaxy distribution. Dekel et al. (1992) used phase diagrams to estimate the confidence levels of the regular signal, which was found in the pencil-beam surveys by Broadhurst et al. (1990), They selected as their statistics the maximum phase amplitudes and the asymmetry of the phase distribution. These statistics are more noisy and their sampling distributions can be found only by simulations. In the contrary, the PDM and epoch folding, which are based on analysis of variance, allow to derive exact sampling distributions.

If we search for a cubic regularity, we have to fold our cluster distribution into a phase cube, as shown in Fig. 2. However, we cannot use directly the methods described above, as we do not observe densities at certain points in space, but discrete objects. However, we could, in principle, calculate a density for any point. As suggested by Dr. J. Pelt, the closest analog to the analysis of variance in our case will be a search for the maximum of the total variance of density in the phase cube. By the way, in the case of time series this variance (the total variance of data) is fixed and does not depend on the trial period.

Figure 2: 2-D cubic folding. The sample volume (upper panel) is divided into smaller trial squares (dotted lines) and all trial squares are stacked together (lower panels). If the size of the trial squares is close to the real period, there will be a clear density enhancement (lower left panel) in the stacked distribution. Otherwise the stacked distribution is almost random (lower right panel). A kernel window is shown in the upper left corner of the stacked cell.

So, the statistic we shall use to estimate the cubic regularity of the point distribution can be written as

$$\kappa(d)=\frac{1}{\bar{n}^2_{\rm p}}\int_{C_{\rm p}}n^2_{\rm p} \left({\vec x}_{\rm p}\right)~{\rm d}{\vec x}_{\rm p},$$ (1)

where d is the trial period (cube size), $C_{\rm p}$ is the phase cube, ${\vec x}_{\rm p}$ are the phase coordinates, obtained from the data coordinates ${\vec x}$ by the folding and scaling transformation

$${\vec x}_{\rm p}={\vec x}/d-\lfloor{\vec x}/d\rfloor,$$ (2)

$n_{\rm p}$ is the number density in the phase cube and $\bar{n}_{\rm p}$ is the mean number density of the whole phase cube. The floor function $\lfloor\cdot\rfloor$ is defined to be the largest integer smaller or equal to x. As the phase coordinates (2) lie in the unit interval, the volume of the phase cube $V_{\rm p}=1$, and the regularity $\kappa=1$for a constant density distribution. Thus there is no regularity in the case $\kappa=1$. If there is a regular cubic signal in the point distribution, it is magnified by folding and by squaring its total amplitude. This makes the statistic extremely sensitive; we shall estimate its sensitivity below.

This statistic is also invariant with respect to translations in data space - the phases of the regular signal will change, but not its amplitude distribution. This is strictly true, however, only if we can neglect edge effects, i.e. if we were able to cut the data cube from an infinite volume.

For simple distributions this statistic can be found analytically. If we have a constant density data cube of a fixed size, the regularity signal will come from edge effects, since the partial cells will create a nonuniform phase cube density distribution. It is a bit tedious to count the different contributions, but it can be done. We give the formula in the Appendix (it is useful for checking the program, should the reader want to write one) and illustrate the result in Fig. 3. Such a figure, the dependence of a statistic on a trial period, is called a periodogram; in our case the regularity periodogram. The regularity according to the definition above is shown by a solid line; the lower dashed line shows a slightly modified statistic where the points that are nearby in real space are ignored when calculating densities. We shall argue below that the latter version works better.

Figure 3: The regularity periodogram for a constant density cube of size D (d is the test cube size, the period). The solid line is for the full density statistic, the dashed line - for the version without local correlations (see Sect. 2.4).

As we see, for a constant density cube the regularity $\kappa$ is unity for the periods where there are no edge effects (the data cube size is an integer multiple of the period). As the period grows, the number of foldings decreases and the relative strength of edge effects increases. A regularity amplitude of 1.1 is already fairly large.

2.2 The estimator

In order to get an estimator for the regularity $\kappa$we have to be able to calculate the density of points in the phase cube. Noting that the integral in the statistic can be written as

$$\begin{eqnarray*}\int_{C_{\rm p}}n^2_{\rm p}~{\rm d}{\vec x}_{\rm p}=\int_{C_{\rm p}}n_{\rm p}~{\rm d}N, \end{eqnarray*}$$

where N is the total number of points, and that

$$\begin{eqnarray*}\bar{n}_{\rm p}=\bar{n}\frac{V}{V_{\rm p}}=N, \end{eqnarray*}$$

($\bar{n}$ is the mean density of the sample, V - its total volume, $V_{\rm p}=1$ the phase cube volume), a natural estimator can be written as

$$\hat{\kappa}(d)=\frac{1}{N^2}\sum_i^N\sum_j^N K\left({\vec x}_i,{\vec x}_j\right),$$ (3)

where N is the total number of points in the sample and ${\vec x}_i$ are their phase coordinates. The kernel function $K({\vec x}_i,{\vec x}_j)$is used to estimate the phase cube density near data points:

$$n_{\rm d}(i)=\sum_j^N K\left({\vec x}_i,{\vec x}_j\right).$$ (4)

Any kernel with compact support can be used here. We used the simplest one

$$\begin{eqnarray*}K\left({\vec x}_i,{\vec x}_j\right)=\frac{1}{V_\varepsilon}\lef... ...}_i,\varepsilon),\\ 0,& \mbox{otherwise}. \end{array} \right. \end{eqnarray*}$$

Here $C({\vec x}_i,\varepsilon)$ is a cube of size $\varepsilon$, centered on ${\vec x}_i$ and with volume $V_\varepsilon=\varepsilon^3$. This definition is symmetric in i,j, although at first glance it does not look as though.

As the phase distribution is periodic by definition, any kernel we use has also to take into account that periodicity, when calculating densities for points near the boundaries of the phase cube. In practice we do this by creating a padding for the phase cube that extends from $-\varepsilon/2$ to $d+\varepsilon/2$, and letting the index j in Eq. (3) to run over the padded phase cube. This speeds up calculations considerably, since we have to create this padding once for every trial period, without that the periodicity conditions should have been checked for every point i in Eq. (4). The program to calculate the periodogram is written in C and it can be obtained by anonymous ftp from ftp.aai.ee:/pub/saar/regularity.tar.gz.

2.3 Error estimates

The probability distribution for the estimator (3) is easy to derive for Poisson-distributed data. As we use the regularity periodogram to search for small regular signals, this case will serve as a good point of reference.

Although we calculate the phase cube density by means of a kernel estimator, this is statistically equivalent to binning the phase cube into $\nu$ bins, where $\nu=1/V_\varepsilon$. In this binning approximation we can write

$$\begin{eqnarray*}\hat\kappa=\frac{1}{N^2\nu}\sum^\nu_{i=1}(n_p)^2_i= \frac{\nu}{N^2}\sum^\nu_{i=1}N^2_i, \end{eqnarray*}$$

where N_i are the occupation numbers of the phase bins and we have used the fact that $\Delta V_{\rm p}=V_\varepsilon=1/\nu$.

Figure 4: A regularity periodogram (left panel) for a Poisson cube with a 2% mixture of points (80 versus 4000) forming a regular distribution near a grid with a spacing of 100 (the cube size is 700). This periodogram is shown by a solid line, the two dotted lines show the range spanned by 100 regularity histograms for pure Poisson data. The right panel shows a slice of a thickness of 150 of the data cube.

In practice, the occupation numbers N_i should be at least a few tens, to ensure proper estimates of the phase cube density. Thus, we can approximate the Poisson distribution of intensity $\lambda=N/\nu$ by the Gaussian one with the mean and variance $\lambda$.

Let us now define a sum

$$X_\nu=\sum^\nu_{i=1}\left(\frac{N_i-\lambda}{\sqrt{\lambda}}\right)^2.$$ (5)

In the Gaussian approximation made above the quantity $X_\nu$ obeys a $\chi^2_\nu$ distribution.

Expanding the sum in (5) and using the expression for $\lambda$, we get

$$\hat{\kappa}=\frac{1}{N}X_\nu+1.$$ (6)

This allows us to estimate the mean and variance of $\hat{\kappa}$ as

$$\begin{eqnarray*}{E}[\hat{\kappa}]=\frac{1}{N}{E}[X_\nu]+1=1+\frac{\nu}{N} \end{eqnarray*}$$

and

$$\mbox{Var}[\hat{\kappa}]=\frac{1}{N^2}\mbox{Var}[X_\nu]=\frac{2\nu}{N^2}\cdot$$ (7)

As we see, the estimate is effective, but biased. The bias may be removed easily, using N_i(N_i-1) in the sum instead of (N_i)². Then

$$\begin{eqnarray*}E[\hat{\kappa}]=\frac{1}{N}E\left[X_\nu\right]+1-\frac{\nu}{N}=1. \end{eqnarray*}$$

The expression for the variance (7) shows that it depends only on the total number of data points and the number of bins used; the variance is uniform over all the period range. This is a very useful property. We tested the above formulae by calculating the regularity periodogram for a large number of simulations of Poisson cubes.

Now we are able to estimate the sensitivity of the estimator. Let us suppose that there is a small cubic signal, consisting of $N_{\rm s}$ additional points each in $\mu$bins. Expanding the sum of squares of bin population numbers for this case and taking into account the change in the normalization factor, we can write for the estimator

$$\hat{\kappa}=1+\frac{\nu}{\mu}\left(\frac{S}{N}\right)^2,$$ (8)

where $S=\mu N_{\rm s}$ is the total number of points in the signal, and $S\ll N$. Hence the amplitude of the signal is proportional to the square of the "signal-to-noise'' ratio S/N and is inversely proportional to the "filling factor'' of the structure $\mu/\nu$. Cubic arrangements that fill only grid vertices are easier to detect than those populating edges or faces of the grid.

Thus we are able to detect a signal, if

$$\begin{eqnarray*}\frac{\nu}{\mu}\left(\frac{S}{N}\right)^2\geq \gamma\sigma(\hat{\kappa}) =\gamma\frac{\sqrt{2\nu}}{N}, \end{eqnarray*}$$

where $\sigma(\hat{\kappa})=\sqrt{\mbox{Var}[\hat{\kappa}]}$ is the rms error of $\hat{\kappa}$, and $\gamma$ is a coefficient that defines the confidence level required. We can write this as

$$S^2\geq\gamma\mu N\sqrt{\frac{2}{\nu}}\cdot$$ (9)

The regularity amplitude for another extreme case, when we are observing a purely cubic structure, is also easy to derive. In this case the amplitude depends only on the "filling factor'',

$$\begin{eqnarray*}\hat{\kappa}=\frac{\nu}{\mu}\cdot \end{eqnarray*}$$

This formula differs from the linear approximation above (8) only by the lack of the additive constant 1.

As an example, we take a typical application to galaxy clusters with $N\approx 1500$, $\nu=N/30\approx 50$, and choose $\gamma=3$. If we suppose that the number of signal bins is around 10, then a cubic signal with a total number of $S\geq 100$ clusters should be detectable already. If we know the amplitude of a periodogram and can estimate the filling factor of the structure, then we can use this formula also to estimate the number of points in the cubic arrangement.

An illustration of the sensitivity of the regularity statistic is shown in Fig. 4, where a tiny regular signal of a total of 80 points is mixed into a Poisson sample of 4000 points (S/N=0.02). These points are distributed near a grid with a spacing of 100, while the total size of the cube is 700. We show also the range spanned by the regularity periodograms of 100 pure Poisson realizations of the data. We see that such a small signal is readily extracted from data, there is a peak with confidence higher than 99% at the period of 100. The filling factor of the cubic structure $\mu/\nu$ is rather small for this arrangement, of course. A clear period will also generate both harmonics and sub-harmonics of smaller amplitude. The d=50 harmonic is well seen in the Figure, while sub-harmonics are hidden by edge effects.

We also show a slice of the data cube in this figure. Note that the regular signal can not be traced there by eye.

There are two possible strategies to determine the size of the bin (the width of the density kernel). The first is to keep this width $\varepsilon$ constant in data space. In that case the effective number of bins $\nu=(\varepsilon/d)^3$(d is the period) and

$$\begin{eqnarray*}\sigma(\hat{\kappa})=\frac{\sqrt{2}}{N}\left(\frac{l}{d}\right)^{3/2} \end{eqnarray*}$$

will depend strongly on the period. Also, a particular choice of $\varepsilon$ will restrict the usable period range. For periods only a few times larger than $\varepsilon$ the density smoothing would be excessive, erasing most of the signal. For larger periods, when the number of points per bin $N_{\rm B}=N(\varepsilon/d)^3$ becomes too small, the estimates may become very noisy.

Another strategy is to keep the phase resolution constant for all periods (the same number of bins, $\nu$). In this case we get similar smoothing for all periods, i.e. the mean number of points per bin $N_{\rm B}$ remains the same. As $\nu=N/N_{\rm B}$, the rms error can be written as

$$\begin{eqnarray*}\sigma(\hat{\kappa})=\frac{2}{N N_{\rm B}}; \end{eqnarray*}$$

the scatter of the estimate remains constant for all periods.

We have used both strategies, starting with that of the constant kernel width in data space, but we realized later that the second strategy is much better.

   
2.4 Effect of local correlations

The error analysis above was made for Poisson data samples. In real applications, our samples have substantial short-range correlations. This will increase both the amplitude and the variance of the estimate, and the clustered nature of the data could generate false regularity signals.

In order to demonstrate this we generated a segment Cox process with the same correlation properties as our cluster sample. The use of segment Cox processes in correlation analysis has been advocated by Pons-Bordería et al. (1999). A segment Cox process places random segments of length l in space and then generates Poisson-distributed points along these segments. An important property of a segment Cox process is that while its long-range distribution is strictly homogeneous ($\xi(r)=0$, if $r\geq l$), its short-range correlation function $\xi(r)\sim r^{-2}$is similar to the correlations we find in observations.

We describe segment Cox processes and discuss how to specify the parameters of the process in the appendix. Here we list those for our sample: the size of the cube is 700, the segment length l=60, the density of segments $\lambda_{\rm S}=1.65\times 10^{-6}$ (the mean number of segments for the sample cube is 566), the line density of points $\lambda_{\rm L}=0.12$ (the mean number of points per segment is 7). The total number of points in the sample is 2057. As this is only an example, we use dimensionless units here, but if we would multiply these by 1 h^-1 Mpc, all these sizes would be comparable to those of cluster samples.

Figure 5: The correlation function (lower panel) for a segment Cox process simulating the spatial distribution of rich clusters of galaxies. The straight line shows the approximation $\xi (r)=(r/20)^{-2}$. The upper panel shows a slice of the simulation cube with a thickness of 1/7 of the cube size.

We show the correlation function of our simulated segment Cox sample and plot a slice of the sample cube in Fig. 5. The spatial distribution of points looks slightly peculiar (like a segment Cox process should), but the correlation function that we show is close to that of rich clusters of galaxies. We plot the function $\xi (r)=(r/20)^{-2}$ there as an example.

Now that we have our simulated Cox sample we can calculate its regularity periodogram. In order to avoid edge effects we choose the periods this time by dividing the size of the data cube by successive integers, so that all cells fit exactly in the cube. This would not be a good recipe for observations, as we loose a lot of interesting periods this way, but here we know that there are no real periods in the data.

We show this periodogram in Fig. 6 with a dashed line. As expected, the estimate is biased, exceeding the expected value of unity for all periods, and the bias grows with period. The latter effect is also easy to understand, as for larger periods the number of stacked cells is smaller and the influence of clusters in individual cells is larger.

We could try to correct for this effect, determining the short-range correlation in advance, and comparing our periodograms with those built for Cox samples. This is not easy, however; while one can usually construct a segment Cox process that describes the two-point correlations of a sample, Cox samples often look rather different from observational ones. Hence one could never be sure if the comparison is a fair one.

Fortunately, there is an easier way to handle this, and again it is a trick known to the periodicity community. As Dr. J. Pelt suggested, we could modify our estimator, discarding any local contributions to the phase cube density. In other words, when calculating the sum in formula (4), we do not count these data points that lie in the same real space cell (of size d) as the point i itself.

Figure 6: Regularity histograms for a segment Cox process simulating short-range correlations. The dotted line shows the periodogram as defined above by formula (3), the solid line - for the estimator that discards local correlations.

This effectively eliminates local correlations and leaves only long-range correlations we are interested in. We show this periodogram as the solid line in Fig. 6. This line wiggles nicely around the expected value in the range predicted by the estimate of the variance, and it does not "feel'' short-range correlations at all.

We have used this estimator in all the rest of the paper. We have to note that in this case the normalization is strictly correct for very short periods only. The sum in the estimator (3) has to be changed to

$$\begin{eqnarray*}\sum_\nu n^2_i \rightarrow \sum_\nu n_i n'_i, \end{eqnarray*}$$

where

$$\begin{eqnarray*}n'_i=n_i\left(1-(d/L)^3\right) \end{eqnarray*}$$

(here d is the period and L is the size of the sample). Thus the expected value of $\hat{\kappa}$will change to

$$\begin{eqnarray*}E[\hat{\kappa}]=1-(d/L)^3=f(d/L). \end{eqnarray*}$$

However, as we have to use the reduced estimator for observed samples, as explained below, we do not have to worry about normalization.

The regularity periodogram for this estimator for a constant density cube can also be found analytically; the formula is given in the appendix and demonstrated in Fig. 3.

2.5 Border and selection effects

It is good, of course, that our estimator is sensitive to extremely small cubic signals. However, the same sensitivity makes it also prone for edge effects as we see in Fig. 4; incomplete boundary cells may hide the signal we are seeking for.

Thus, first of all we have to account for edge effects very carefully. The case illustrated in Fig. 4 is a little extreme, as the incomplete cells form a clean cubic signal, but it shows clearly the possible amplitude of edge effects. One possibility would be to use a minus-estimator, discarding all those cells that intersect the boundary of the sample. This is a mathematically elegant possibility. However, it is of little use in practice, since boundaries of real samples are usually so complex that we would have to discard most of the data.

To show the effect of selection, we generated a spherical sample of the same linear dimensions, with a diameter of 700 and the same number of points, 2000. The density of this sample falls linearly with the radius, going to zero at the boundary of the sphere. This mimics observational samples rather well. Ten periodograms for this sample are shown in Fig. 7. The main effect is a radical change in the shape of the periodogram. Secondly, the boundary-generated oscillations are gone, as the boundaries are fuzzy. And thirdly, the scatter of the curves remains constant. This is expected, as both the total number of points (2000) and the number of points per density bin (30) are the same.

Selection effects could be accounted for by introducing weights for data points that are inversely proportional to the normalized density of the sample (the selection function) at that point $\bar{n}_i$, defining:

$$\begin{eqnarray*}\hat{\kappa}'(d)=\frac{1}{N(N-1)}\sum_i^N\sum_j^N\frac{w_i}{\ba... ... \frac{w_j}{\bar{n}_j}K\left({\mathbf x}_i,{\mathbf x}_j\right), \end{eqnarray*}$$

where the additional weights w_i can be chosen to minimize the sampling variance of $\hat{\kappa}'(d)$. As these weights have to be calculated only once, it also would speed up the calculation of periodograms. The trouble with this approach is that the choice of the weights is rather difficult. This approach also mimics a volume-limited constant mean density sample, leading to abrupt boundaries and thus amplifying edge effects.

An easier approach is to build reference samples using the selection functions determined for data samples, and to account for edge effects by comparing the real sample and a reference sample with the same geometry, using a reduced regularity periodogram k:

$$k(d)=\frac{\hat{\kappa}(d)_{{\rm obs}}}{\hat{\kappa}(d)_{{\rm ref}}},$$ (10)

where the subscripts "obs'' and "ref'' denote the observed data and simulated reference samples, respectively. In order not to introduce discreteness noise by this operation, the number of points in the reference sample should be larger than in the data sample, but not too much (the variance of the estimate is proportional to N^-2). We have used in this paper reference samples with ten times more points than data samples have. We shall use this reduced periodogram throughout the rest of the paper, calling it simply the regularity periodogram. The reduced regularity k(d) is an unbiased estimator of the regularity $\kappa(d)$, but its variance will grow as $\hat{\kappa}^{-2}(d)_{{\rm ref}}$. This limits the range of the usable periods.

Figure 7: Regularity periodograms for ten realizations of a sphere of diameter 700, with 2000 particles generated to represent a density falling linearly with radius; the limiting density at the surface of the sphere is zero. Notice the absence of edge effects.

The reference periodogram will take account of both the edge and selection effects. Interestingly, the reduced estimator (10) might be the only one used in astronomy that actually gains from selection effects - the fact that the density of data points gets smaller at the sample boundaries helps to reduce significantly the amplitude of edge effects.

Apart from the spurious signal the boundaries create, they have another effect - they destroy the translational invariance of the statistic. This can be most easily seen when folding a one-dimensional analog of a data cube, a constant density line segment. This means that we should create our reference samples very carefully, as any small shift in the shape between the reference sample and the observed sample would give a spurious signal. Fortunately, fuzzy boundaries of observed samples, due to selection effects, practically eliminate this problem.

We have tested this by shifting the reference samples generated to use with the observed galaxy cluster samples. The geometry of these samples is rather complicated and slight shifts between the observed and simulated shapes could easily happen. The test showed that shifts up to 20 h^-1 Mpc in every coordinate direction (the samples consist of two cones with characteristic sizes of 350 h^-1 Mpc) practically do not change regularity periodograms.

2.6 Voronoi samples

The regularity statistic is a typical universal frequentist statistic that can be applied to all point configurations. If we were using the Bayesian approach, we also should have specified a statistical model. However, this universality has also a drawback - without a statistical model we do not know if there exist non-cubic arrangements which also will produce a cubic signal. Thus, in order to understand the proposed statistic, we have to test it on various spatial arrangements.

The first class of arrangements we use are the Voronoi models of large-scale structure, introduced by van de Weygaert & Icke (1989) and van de Weygaert (1991).

These models are based on the Voronoi tessellation of space. In a Poisson-Voronoi point process the centers of voids serve as seeds of structure, and have a Poisson distribution, while the points (clusters of galaxies in our case) are located at the vertices of the structure formed by the expanding voids. Although the model starts from a Poisson distribution, it is well clustered and has a power-law correlation function $\xi\sim r^{-2}$, close to that of the galaxy clusters. This amazing fact was discovered by van de Weygaert & Icke (1989) and verified exactly by Heinrich et al. (1998). In this model clusters of galaxies form superclusters, which together with voids form a cellular supercluster-void network, similar to that observed. The only free parameter of the model is the mean density of void centers that determines the mean diameter of voids.

We used a program by R. van de Weygaert to generate the cluster sample. The size of the simulation box was 700 h^-1 Mpc; the number of seeds was 431, chosen to obtain a mean diameter of voids 115 h^-1 Mpc in accordance with the distribution of Abell clusters (Paper I). We generated ten realizations of the model; the number of clusters in the samples varied between 2884 and 2948 (for details see Einasto et al. 1997c, hereafter Paper III).

To characterize the regularity of a model we calculated the reduced regularity periodograms k(d) and the correlation functions $\xi(r)$ of clusters. The results for ten realizations are shown in Fig. 8. The correlation function has a deep minimum around 80 h^-1 Mpc and a secondary maximum at 140 h^-1 Mpc; on larger scales it flattens out to a mean value around zero. The geometric interpretation of this behaviour was discussed in Paper III: the minimum corresponds to the mean distance between superclusters and surrounding voids; the secondary maximum can be identified with the mean distance between superclusters across voids. In a Voronoi model the mean size of voids is well fixed, thus the presence of a secondary maximum is expected. On still larger scales the behavior of the correlation function depends on the regularity of the distribution of rich superclusters. In a Voronoi model void centers as well as rich superclusters are randomly distributed, thus we expect no correlation on very large scales. As we see from Fig. 8, on very large scales the correlation function is indeed close to zero which corresponds to a Poisson distribution. In contrast to the correlation function the features in the regularity periodograms k(d)have a rather low amplitude (there are suspicious signals at 190 h^-1 Mpc and 220 h^-1 Mpc), but overall the periodograms are featureless, resembling the regularity periodograms of pure Poisson samples. Although the Voronoi model is cellular, the cells are not cubic or regular.

Figure 8: Regularity periodograms (upper panel) and correlation functions (lower panel) for ten realizations of Voronoi models.

2.7 Quasi-regular models

Figure 9: Regularity periodograms (upper panel) and correlation functions (lower panel) for ten realizations of quasi-regular line models.

Next we shall choose a quasi-regular model distribution. We use mixed models, i.e. samples with two populations of clusters. The clusters of one population are randomly distributed, while in the second population they are located in superclusters which form a regular rectangular network with a period of 130 h^-1 Mpc. In our model superclusters are randomly located along rods which form a regular rectangular network with a period of 130 h^-1 Mpc. The positions of the rods have been randomly shifted with the rms deviation of $\pm 20$ h^-1 Mpc (for details see Paper III).

The box size was taken as 690 h^-1 Mpc, the period was set to 115 h^-1 Mpc, the number of randomly located poor superclusters was chosen to be 3000, and the minimal number of rich superclusters on rods was 4. We generated ten realizations of the model; the total number of clusters varies between 9218 and 9519, and the total number of clusters in quasi-regular population in rods varied between 3328 and 4216. The regularity periodograms and the correlation functions are shown in Fig. 9. The correlation functions are oscillating with a period of 115 h^-1 Mpc; the amplitude of oscillations decreases very slowly with distance r. The regularity periodograms have two well-pronounced maxima at d=115 h^-1 Mpc and d=230 h^-1 Mpc. The scatter of the regularity periodograms is very small.

2.8 Angular sensitivity

Obviously, the method is sensitive to the direction of the axes of the trial cubes. The orientation of a cube can be described by three Euler angles, so there is a considerable amount of freedom here. If the trial cubes are rotated $45^{\circ }$ along one axis of the cubic alignment, there will be a mixture of cubic signals - one with the previous period from the direction perpendicular to that face, and another with the period of $\sqrt{2}$ times larger. If the trial cubes are oriented along the long diagonals, from a vertex to the opposing vertex of the original cube, the cubic network will be cubic again with a single period of $\sqrt{3}$times the original one. For other orientations the signal will be weak or absent. This property is illustrated in Fig. 10. The reason for this behaviour is clear: the value of k(d) depends on the number of clusters which coincide in the stacked cell. If clusters located in superclusters of different original cells lie in the stacked cell in different locations (which happens when orientations differ considerably), then k(d)cannot exceed by much the expected value of unity for random samples.

Figure 10: Angular sensitivity of regularity periodograms (for a quasi-regular model). The solid, dashed and dot-dashed lines are for trial cubes oriented at $\theta =0^{\circ }$, $22.5^{\circ }$ and $45^{\circ }$ with respect to the grid.

2.9 $\mathsfsl{N}$-body cluster samples

We have compared the regularity periodograms and correlation functions also for several cluster samples found in numerical simulations of the large-scale structure of the Universe. We used four models analyzed by E99. These models were calculated using a PM code with 128³ particles and a 256³ mesh in a cube of size 720 h^-1 Mpc. The models include the standard CDM model with the density of matter $\Omega_{\rm M}=1$ (model CDM61), and a LCDM model with the vacuum energy term, $\Omega_{\rm M}=0.2$ and $\Omega_{\Lambda}=0.8$ (model CDM62). To obtain a model which better represents the observed power spectrum (see E99) we used two initial power-law spectra with the indices n=1 on large scales, n=-1.5 on small scales, and a sharp transition on a scale of k=0.05 h Mpc^-1 (models DPS6 and DPS3, the latter with 256³ particles). The last model we used (DPP6) has an additional high-amplitude peak in the power spectrum near the maximum. In most of the models clusters of galaxies were selected using the friends-of-friends algorithm. Only in the high-resolution model DPS3, computed by Gramann & Suhhonenko (1999), maxima of the smoothed density field were used for cluster identification.

Figure 11: Regularity periodograms (upper panel) and correlation functions (lower panel) for simulated CDM and double-power-law models.

Regularity periodograms and correlation functions for simulated cluster samples are shown in Fig. 11. These functions have been calculated for the whole cluster sample in the simulation box. The correlation functions of these model samples oscillate slightly, except for the model DPP6 with a high peak in the power spectrum that causes high-amplitude oscillations. The regularity periodograms of all models are practically featureless, even that for the model DPP6; small peaks and valleys are due to shot noise. This is understandable, as the initial density fluctuations are, by definition, isotropic and cannot develop into a cubic structure.

Thus, application of the regularity statistic to these three different, but typical classes of spatial arrangement of objects in cosmology, shows that it well detects the cubic signal, if it is present, and does not "feel'' other regularities.