text:
# Quantifying angular clustering in wide-area radio surveys

## 1 Introduction

**footnotetext: E-mail: cab@astro.ox.ac.uk

Describing the large-scale structure of the Universe is of fundamental importance for testing theories of galaxy and structure formation and for measuring the cosmological parameters. The largest structures require delineation by the deepest, widest surveys, currently represented by surveys for radio AGN. These contain objects to redshifts of at least $`z4`$: the radio emission marking these objects is not affected by dust obscuration, large-scale calibration effects should be minimal, and the number of objects in the current generation of radio surveys such as WENSS, FIRST and NVSS reach $`10^6`$ over substantial fractions of the sky.

We see the distribution of galaxies projected on the sky, but this is still useful to quantify: it is easy to assemble a large sample of objects and the angular clustering can be de-projected (in a global statistical manner) to measure the spatial clustering, conclusions being reached in the absence of complete redshift information. There are many sophisticated methods for quantifying the angular distribution of galaxies. These include spherical harmonic analysis (e.g. Baleisis et al. 1998), percolation analysis (Bhavsar & Barrow 1983) and minimal spanning trees (Krzewina & Saslaw 1996). Chiang & Coles (2000) emphasize the importance of maintaining the phase information of the clustering for describing morphology. In contrast, here we use two of the crudest statistics for describing angular structure: the two-point angular correlation function and the moments of counts-in-cells. It is well-known that these methods lose much of the clustering information: two very different distributions can have the same two-point correlation function. However, these statistics are simple to interpret and hence reveal the fundamental observational problems and survey limitations. They provide simple points of contact with prediction, have well-understood statistical errors and together provide a consistency check. They must be understood and must give consistent results before application of more powerful techniques can be considered.

Correlation function analyses (Peebles 1980), widely used since the early days of clustering investigations, have been extensively applied in the optical regime, for example to the APM survey (Maddox et al. 1996). Here, the correlation function typically shows a power-law behaviour $`w(\theta )\theta ^{1\gamma }`$ with $`\gamma 1.8`$ on small scales, with a steepening break to larger scales. The key difference between angular correlation function analyses in the optical and radio regimes is in the latter, the wide redshift range of radio sources washes out much of the clustering amplitude through the superposition of unrelated redshift slices. Hence an angular clustering signal has only been measurable in the most recent radio surveys, initially with marginal detections in the Green Bank 87GB survey (Kooiman et al. 1995) and the Parkes-MIT-NRAO (PMN) survey (Loan et al. 1997).

The latest generation of deep radio surveys – FIRST (Becker et al. 1995), WENSS (Rengelink et al. 1998) and NVSS (Condon et al. 1998) – reveal the imprint of structure more clearly. The correlation function has been measured for WENSS by Rengelink et al. (1998), and for FIRST by Cress et al. (1996) and, in a pioneering and innovative series of papers, by Magliocchetti et al. (1998). These studies concluded that the slope of the correlation function for radio galaxies was steep ($`\gamma >2`$). Cress & Kamionkowski (1998) and Magliocchetti et al. (1999) modelled the 3D clustering from these analyses, including the behaviour of bias with epoch. Magliocchetti et al. (1998) also carried out a counts-in-cells analysis of the FIRST survey, detecting significant skewness.

These results motivated our present investigation. The NVSS had never been investigated for large-scale structure effects, and we wished to determine if the more extensive sky coverage and source list of $`2\times 10^6`$ objects led to conclusions compatible with FIRST, with higher signal-to-noise ratio affording further insight. We wished to understand how robust the conclusions from FIRST were, given the issue of over-resolution and the consequent need to “combine” catalogue sources from multiple-component radio galaxies. We wished to examine the compatibility of results from counts-in-cells and correlation-function analyses. With these aims in mind, NVSS and FIRST, at the same frequency but at resolutions differing by a factor of 9, suggest an ideal comparative study. Our initial results (measurement of the NVSS angular correlation function) were presented in Blake & Wall (2002).

To proceed we first describe the two surveys, NVSS and FIRST. Section 3 summarizes the clustering statistics we use and Section 4 discusses the observational issues bound to impact upon large-scale structure analyses. Sections 5 and 6 derive angular correlation functions and counts-in-cells for each of NVSS and FIRST, and Section 7 compiles the conclusions.

## 2 The radio surveys: NVSS and FIRST

### 2.1 NVSS

The NVSS (NRAO VLA Sky Survey, Condon et al. 1998) was carried out with the VLA at an observing frequency of 1.4 GHz over the period 1993 – 1996 and covers the whole sky north of declination $`40^{}`$ (33,884 square degrees or 82 per cent of the celestial sphere). The source catalogue contains $`1.8\times 10^6`$ sources and is claimed to be 99 per cent complete at integrated flux density $`S_{1.4\mathrm{GHz}}=3.5`$ mJy and 50 per cent complete at 2.5 mJy. These figures are differential completenesses, i.e. 99 per cent of all sources with $`S_{1.4\mathrm{GHz}}=3.5`$ mJy appear in the NVSS catalogue. The survey was performed with the VLA in D configuration, with DnC configuration used for fields at high zenith angles ($`\delta <10^{},\delta >78^{}`$), and the FWHM of the synthesized beam is about 45 arc-seconds. The raw fitted source parameters are processed by a computer program provided by the survey team called NVSSlist, which performs the deconvolution and corrects for known biases to produce source diameters and integrated flux densities. Details are given by Condon et al. (1998); NVSSlist version 2.16 (March 2001) was used for this investigation. Figure 1 is a plot of NVSS catalogue entries with $`S_{1.4\mathrm{GHz}}>200`$ mJy.

### 2.2 FIRST

The FIRST (Faint Images of the Radio Sky at Twenty centimetres) survey, also carried out at an observing frequency of 1.4 GHz but with the VLA in B configuration, began in 1993 and continues. It covers the North and South Galactic caps and the B configuration of the VLA yields an angular resolution of about 5 arcsec. Details of the survey design, analysis and catalogue generation are given in Becker, White & Helfand (1995) and White et al. (1997). These papers claim that the survey is 95 per cent complete at $`S_{1.4\mathrm{GHz}}=2`$ mJy and 80 per cent complete at 1 mJy. These figures are cumulative completenesses, i.e. 95 per cent of sources with $`S_{1.4\mathrm{GHz}}>2`$ mJy appear in the FIRST catalogue.

The latest publicly-available catalogue dates from 15 October 2001. It contains 771,076 sources and covers a total of 8565 square degrees (7954 in the north galactic gap and 611 in the south), or 21 per cent of the celestial sphere. The raw catalogue contains a number of spurious entries representing sidelobe responses from nearby brighter sources. As described in White et al. (1997), the FIRST survey team developed an oblique decision-tree program to identify and flag these sidelobes. In the 15 October 2001 catalogue, 28,017 sources are flagged and were excluded from our analysis.

In Figure 2 we plot FIRST catalogue entries with $`S_{1.4\mathrm{GHz}}>50`$ mJy in the northern sky, and we outline the contiguous region for which our large-scale structure analysis was performed ($`123^{}<\alpha <247^{}`$, $`5^{}<\delta <58^{}`$).

## 3 Quantifying the angular clustering of galaxies

A common way of quantifying the clustering in an angular distribution of galaxies is by using the angular correlation function, $`w(\theta )`$. This compares the observed (clustered) distribution to a random (unclustered) distribution of points across the same survey area, by simply measuring the fractional increase in the number of close pairs separated by angle $`\theta `$. Specifically, if we let $`DD(\theta )`$ be the number of unique pairs of galaxies with separations $`\theta \theta +\delta \theta `$, and $`RR(\theta )`$ be the number of random pairs in the same separation range, then $`w(\theta )`$ can be estimated as $`w(\theta )=(DDRR)/RR`$. Further investigation reveals that the statistical error on $`w`$ may be minimized by averaging over a large number of random sets ($`RR\overline{RR}`$) and by using a different estimator for $`w`$ (see Landy & Szalay 1993). We adopt the Landy-Szalay estimator for our investigations, as this has minimal statistical bias and variance.

Another simple way to quantify the galaxy distribution is to grid the sky into cells of fixed area and shape, and count the number of sources that lie in each cell. Counts-in-cells yields the probability distribution $`P(N)`$ of finding $`N`$ sources in a cell, and the moments of the distribution such as the variance $`\mu _2=\overline{(N\overline{N})^2}`$ and skewness $`\mu _3=\overline{(N\overline{N})^3}`$ (the horizontal bar indicates an average over cells). A clustered distribution produces a higher variance than a random distribution because cells may lie in clusters or voids, broadening the probability distribution $`P(N)`$. Skewness is important because, assuming Gaussian primordial perturbations and linear theory, the skewness of counts-in-cells remains zero (Peebles 1980). Measurement of a non-zero skewness therefore indicates either non-linear gravitational clustering or non-Gaussian initial conditions.

Whereas the angular correlation function bins pair separations into small intervals, a counts-in-cells analysis combines information from a range of angular scales up to the cell size, effectively measuring an average $`w(\theta )`$ (see equation 1). By avoiding the binning of angular separations, the counts-in-cells is less affected by “shot noise” and is a more sensitive probe of long-range correlations. A simple relation exists between $`w(\theta )`$ and $`\mu _2`$ ($`\mathrm{\S }6.1`$), hence we can verify that they are consistent for a given distribution. Neither quantity provides a complete statistical description of the galaxy distribution: this can be encompassed by the hierarchy of correlation functions or the full probability distribution $`P(N)`$; both are harder to interpret and to compare with theory. Nor are $`w(\theta )`$ and $`\mu _2`$ well-suited for describing the morphology of the distribution (do galaxies cluster in filaments, sheets or clumps?) or its topology (how do the filaments or sheets join up to form the global pattern?). However, these statistics are easy to measure, provide a simple point of contact with prediction, and can reveal observational problems with the survey data.

## 4 Observational effects impacting on large-scale structure studies

### 4.1 Resolution effects

The NVSS and FIRST surveys are at the same observing frequency (1.4 GHz) but at resolutions differing by a factor of 9, and therefore provide an excellent comparative study of survey resolution effects on the observed properties of radio galaxies. In Figure 3 we plot NVSS and FIRST catalogue entries with integrated flux densities $`S_{1.4\mathrm{GHz}}>3.5`$ mJy (at which threshold both surveys are claimed to be complete) in a randomly chosen $`3^{}\times 2^{}`$ patch of sky. This region contains 230 FIRST entries and 228 NVSS entries. We see that:

* Most radio sources are detected in both surveys.
* Of the objects appearing in just one survey, many more appear in NVSS than FIRST.
* The surface densities are almost identical because many FIRST objects are multiple components of the same galaxy.

These facts can be explained by the superior angular resolution of FIRST, which picks up multiple radio components unresolved by NVSS. But this high resolution provides FIRST with much poorer sensitivity to surface brightness, losing flux from extended sources even well above the survey limit. Thus the objects appearing solely in NVSS in Figure 3 have fallen below the 3.5 mJy threshold in FIRST. The small number of sources which only appear in FIRST are due to either statistical flux-density errors, source variability, unrecognised sidelobes or noise spikes.

The fact that FIRST under-estimates flux densities is confirmed by Figure 4, which plots the average ratio of FIRST to NVSS integrated flux density for sources matched between the surveys for different NVSS flux-density bands. The matched sources are chosen to be isolated in both surveys (having nearest neighbours more distant than 2 arcmin) to prevent confusion from galaxies resolved as multiple radio components. Figure 4 shows that on average FIRST loses 10 per cent of the flux density from 3.5 mJy sources, a fraction which decreases with increasing flux density.

In Figure 5 we overplot differential source counts for FIRST and NVSS for 1 mJy $`<S_{1.4\mathrm{GHz}}<`$ 1 Jy. The curves are in rough agreement. However, the NVSS source count has an unphysical shape below $`S_{1.4\mathrm{GHz}}10`$ mJy. This distortion arises because NVSSlist uses a different algorithm to convert raw fitted peak amplitudes to integrated flux densities, depending on whether a source is classified as extended or not (Condon et al. 1998). In addition, a larger number of bright ($`>100`$ mJy) sources are detected in NVSS than in FIRST. This is due to the finer angular resolution of FIRST breaking up bright sources into more radio components.

These comparisons demonstrate the significant impact of survey resolution and surface brightness sensitivity on the observed properties of radio galaxies. However, the effect of these flux biases on the deduced large-scale clustering should be minimal: the broadness of the radio galaxy luminosity function ensures that the observed clustering is not a strong function of flux density above 3.5 mJy.

### 4.2 Gradients in source surface density

Both the NVSS and FIRST surveys suffer from systematic fluctuations in source surface density across the sky above flux-density thresholds at which they are complete, affecting any attempt to quantify the large-scale structure present. Figures 6 and 7 illustrate these variations. In both surveys the magnitude of the effect depends on the flux-density threshold. The FIRST survey contains 10 per cent fluctuations in surface density over the sky at the stated completeness limit $`S_{1.4\mathrm{GHz}}=2`$ mJy, dropping to below 5 per cent at 10 mJy. The variations appear to correlate with the different observing periods. The NVSS suffers from 2 per cent fluctuations at the completeness limit $`S_{1.4\mathrm{GHz}}=3.5`$ mJy, with significant density steps at the declinations where the array configuration is changed from D to DnC, but these fluctuations have become insignificant by 15 mJy. These NVSS effects originate from difficulties in compensating for the sparse $`uv`$-plane coverage of the survey, constructed from “snapshot” radio observations, and are purely declination-dependent: the projection of array baselines on the sky changes with declination, whereas the data acquisition and analysis software should be insensitive to the right ascension of sources. Note that the NVSS flux-density errors are dominated by a constant additive bias that affects weak sources much more strongly than powerful sources (unlike a multiplicative calibration error).

A varying source density $`\sigma `$ will spuriously enhance the measured value of the angular correlation function $`w(\theta )`$. This is because the number of close pairs of galaxies in any region depends on the local surface density ($`DD\overline{\sigma ^2}`$), but the number of close pairs in the comparison random distribution depends on the global average surface density ($`RR(\overline{\sigma })^2`$). Systematic fluctuations mean that $`\overline{\sigma ^2}>(\overline{\sigma })^2`$, thus $`w(\theta )`$ is increased. The variance of counts-in-cells, as quantified by the statistic $`y(L)`$6.1), will also be increased: a spread in the mean surface density across the cells will inevitably broaden the overall probability distribution $`P(N)`$, which is constructed from fluctuations about those means. We can show that on angular scales less than those on which $`\sigma `$ is varying, both $`w(\theta )`$ and $`y(L)`$ are subject to the same spurious constant offset $`\overline{\delta ^2}`$, where $`\delta =(\sigma \overline{\sigma })/\overline{\sigma }`$ is the surface overdensity.

To estimate the magnitude of this effect, take a simple toy model in which a survey is divided into two equal areas between which there is a fractional surface density shift $`ϵ`$. For this model, $`\overline{\delta ^2}=ϵ^2/4`$. Thus for example the expected offsets on $`w(\theta )`$ are $`\overline{\delta ^2}2.5\times 10^3`$ for FIRST sources above 2 mJy ($`ϵ0.1`$) and $`\overline{\delta ^2}2\times 10^4`$ for NVSS sources above 3.5 mJy ($`ϵ0.03`$).

In this study, analysis of the radio surveys was restricted to flux-density ranges for which surface gradients were negligible. The alternative approach is to modulate the random comparison sets with the same surface gradients as contained in the data. The gradients are not known in advance and must be measured from the data itself; it was found that this could not be done sufficiently accurately to subtract the offset completely.

### 4.3 Multiple-component sources

The large linear sizes and complex morphologies of radio sources mean that a single radio galaxy can be resolved in a radio survey (and appear in a survey catalogue) as two or more closely-separated components of radio emission. When investigating the clustering of individual galaxies, these multiple-component sources will produce spurious clustering at small separations. There is no set of criteria that will reliably distinguish multiple-component sources from closely-separated independent galaxies; we instead choose to incorporate their effect into the fitted clustering model. As it turns out, this effect may be successfully disentangled from that of galaxy-galaxy clustering.

The influence of multiple-component sources on the angular correlation function is quantified in detail in our analysis of the NVSS $`w(\theta )`$ (Blake & Wall 2002), in which it is stressed that even a tiny fraction of giant radio galaxies can substantially affect clustering measurements at small angles. This is because a very small number of close pairs determine the value of $`w(\theta )`$: we find that we cannot neglect the effect of radio galaxies of size $`\theta `$ until $`\theta 0.1^{}`$. At angles $`\theta <0.1^{}`$, $`w(\theta )`$ effectively measures the size distribution of radio galaxies; at $`\theta >0.1^{}`$, the clustering of individual galaxies dominates the pair count. This reasoning is evidenced by the observed numbers and sizes of giant radio sources (Lara et al. 2001) as well as a clear break in the measured $`w(\theta )`$ at $`\theta 0.1^{}`$ (see Figure 8). Hence the influence of multiple-component sources may be disentangled from that of galaxy clustering.

At angles where the clustering of individual galaxies dominates the pair count, the fact that these galaxies may be split into multiple radio sources is unimportant. For if the mean number of radio components per galaxy is $`\overline{n}`$, then the number of pair separations $`N_p`$ at any angle is increased by a factor $`(\overline{n})^2`$, and the source surface density $`\sigma `$ is increased by a factor $`\overline{n}`$. As $`N_p\sigma ^2\times (1+w)`$, the measured correlation function is unaffected.

The effect of multiple-component sources on the moments of counts-in-cells is modelled in Section 6.2.

## 5 The radio galaxy angular correlation function

### 5.1 Measurement of the NVSS angular correlation function

Our measurements of the NVSS $`w(\theta )`$ for different flux-density thresholds are described in Blake & Wall (2002); the results for 10 mJy and 20 mJy are compared in Figure 8. The angular correlation function at all thresholds can be fit by a sum of two power-laws. A convincing interpretation is that the steep small-angle power-law is created by multiple-component radio sources and hence indicates their size distribution, whereas the shallow large-angle power-law describes the clustering between different radio galaxies.

Table 1 displays the results of fitting the function $`w(\theta )=A\theta ^\alpha +B\theta ^\beta `$ to the measurements for flux-density thresholds $`S_{1.4\mathrm{GHz}}=10,15,20`$ mJy. The fits were performed to angles $`\theta >0.02^{}`$, safely above the resolution limit of the NVSS ($`\theta _{\mathrm{res}}=0.0125^{}`$). Figure 8 reveals a fall-off in $`w(\theta )`$ with decreasing $`\theta `$ in the lowest separation bins. This is not real, but is a signal-to-noise problem caused by the failure of the survey to resolve weak double sources with separations slightly greater than the beam-width. As $`\theta `$ increases, a rapidly increasing fraction of doubles of size $`\theta `$ can be successfully resolved.

For all flux-density thresholds, the slope of the clustering power-law is consistent with $`\alpha =0.8`$ (in agreement with other classes of objects) with an amplitude $`A1\times 10^3`$ (with $`\theta `$ in degrees). Blake & Wall (2002) describe a preliminary analysis of the implications for spatial clustering.

### 5.2 Comparison with the FIRST angular correlation function

For comparison, we measured $`w(\theta )`$ from the FIRST catalogue over the region $`123^{}<\alpha <247^{}`$, $`5^{}<\delta <58^{}`$ for all objects above flux-density thresholds $`S_{1.4\mathrm{GHz}}=10`$ mJy and 2 mJy. As a precautionary measure, we placed circular masks of radius $`0.5^{}`$ around all sources with $`S_{1.4\mathrm{GHz}}>1`$ Jy; this left respectively 88,873 and 305,872 objects at the two thresholds. Figure 9 displays the results, compared to the best-fitting NVSS sum of two power-laws at 10 mJy. The FIRST measurement at $`\theta 0.1^{}`$ is known to be contaminated with sidelobes as described by Cress et al. (1996).

The FIRST 10 mJy and NVSS 10 mJy measurements agree well at small angles, in the regime of multiple-component sources. At bigger angles, the large FIRST error bars mean that $`w(\theta )`$ is poorly constrained, although the measurement is consistent with the NVSS result ($`\chi _{\mathrm{red}}^2=1.84`$ for $`\theta >0.02^{}`$ excluding the two points at $`\theta 0.1^{}`$). The FIRST 2 mJy $`w(\theta )`$ measurement has much smaller errors but is offset by the source surface density gradients described in Section 4.2. To verify this, we re-ran the 2 mJy analysis for the more uniform declination region $`42^{}<\delta <57^{}`$, prompted by Figure 7, using right ascension range $`107^{}<\alpha <263^{}`$. The re-determined FIRST 2 mJy $`w(\theta )`$ is in much better agreement with the NVSS result at large angles (see Figure 9). The amplitude of the small-angle FIRST $`w(\theta )`$ drops between 10 mJy and 2 mJy due to the increased surface density, consistent with the hypothesis that $`w(\theta )`$ in this regime is entirely governed by multiple-component sources (see Blake & Wall 2002).

## 6 Comparison with the radio galaxy counts-in-cells

We also performed a counts-in-cells analysis of the radio surveys. Our motivation was twofold: to quantify the clustering imprint with an independent statistic to compare with $`w(\theta )`$, and to make contact with previous work (Magliocchetti et al. 1998, 1999).

### 6.1 Relation of counts-in-cells variance to $`w(\theta )`$

Consider an unclustered distribution of sources, distributed randomly and independently with surface density $`\sigma `$. The expected number of sources in a cell of area $`S`$ is $`<N>=\sigma \times S`$. The expected probability distribution $`P(N)`$ is the Poisson distribution with mean $`<N>`$ and variance $`<N>`$. We define the following statistic to quantify the increased variance of a clustered distribution:

$$y=\frac{\mu _2\overline{N}}{(\overline{N})^2}$$

Hence $`<y>=0`$ for no clustering, as $`<\mu _2>=<\overline{N}>=<N>`$. (Actually, as there is statistical error in the denominator, $`y`$ contains a slight bias which may be neglected). We can show (Peebles 1980, equation 36.6) that for a given $`w(\theta )`$, the expected value of $`y`$ is

$$<y>=\frac{_{\mathrm{cell}}_{\mathrm{cell}}w(\theta )𝑑S_1𝑑S_2}{S^2}=w(\theta )\frac{dG_p}{d\theta }𝑑\theta $$
(1)
where in the final expression, $`dG_p`$ is the fraction of all pairs of area elements within the cell lying in the separation range $`\theta \theta +d\theta `$ (using the notation of Landy & Szalay 1993). $`dG_p`$ may be calculated analytically for simple geometries:

* For a square of side $`L`$, putting $`x=\theta /L`$,

$$\frac{dG_p}{dx}=\{\begin{array}{cc}2x(\pi 4x+x^2)\hfill & 0x1\hfill \\ 2x[\pi 2+4\sqrt{x^21}4\mathrm{arccos}(1/x)x^2]\hfill & 1x\sqrt{2}\hfill \end{array}$$
(2)
* For a circular cell of diameter $`L`$, putting $`x=\theta /L`$,

$$\frac{dG_p}{dx}=\frac{16x}{\pi }\left(\mathrm{arccos}(x)x\sqrt{1x^2}\right)$$
(3)

We now calculate $`<y>`$ for a power-law angular correlation function. If the survey has angular resolution $`\theta _{\mathrm{res}}`$ then

$$w(\theta )=\{\begin{array}{cc}1\hfill & \theta <\theta _{\mathrm{res}}\hfill \\ (\theta /\theta _0)^\alpha \hfill & \theta >\theta _{\mathrm{res}}\hfill \end{array}$$

Hence from equation 1,

$$<y>=_0^{\theta _{\mathrm{res}}}\frac{dG_p}{d\theta }𝑑\theta +_{\theta _{\mathrm{res}}}^{\theta _{\mathrm{max}}}\left(\frac{\theta }{\theta _0}\right)^\alpha \frac{dG_p}{d\theta }𝑑\theta $$

where $`\theta _{\mathrm{max}}=L`$ for a circular cell of diameter $`L`$ and $`\theta _{\mathrm{max}}=L\sqrt{2}`$ for a square cell of side $`L`$. To unveil the clustering pattern we measure the variation of $`<y>`$ with cell size $`L`$, for fixed cell shape. This is elegantly done by the substitution $`x=\theta /L`$:

$$<y(L)>=_0^{\theta _{\mathrm{res}}/L}\frac{dG_p}{dx}𝑑x+\left(\frac{\theta _0}{L}\right)^\alpha _{\theta _{\mathrm{res}}/L}^{x_{\mathrm{max}}}x^\alpha \frac{dG_p}{dx}𝑑x$$
(4)
The quantity $`dG_p/dx`$ is determined purely by the cell shape, not size (see equations 2 and 3). For small $`x`$, $`dG_p/dx=kx`$ independently of cell shape ($`k=8`$ for circular cells and $`k=2\pi `$ for square cells). Assuming that $`\theta _{\mathrm{res}}L`$, this solves the first integral of equation 4:

$$<y(L)>=\frac{k}{2}\left(\frac{\theta _{\mathrm{res}}}{L}\right)^2+\left(\frac{\theta _0}{L}\right)^\alpha _{\theta _{\mathrm{res}}/L}^{x_{\mathrm{max}}}x^\alpha \frac{dG_p}{dx}𝑑x$$

* If the slope of $`w(\theta )`$ is steep enough ($`\alpha >2`$), the non-Poisson clustering will be dominated by close pairs ($`x1`$) and for all cell shapes we can neglect edge effects by assuming that $`dG_p/dx=kx`$ and $`x_{\mathrm{max}}=\mathrm{}`$. This allows us to solve the second integral of equation 4 and obtain:

$$<y(L)>=k\left(\frac{\theta _{\mathrm{res}}}{L}\right)^2\left[\frac{1}{(\alpha 2)}\left(\frac{\theta _0}{\theta _{\mathrm{res}}}\right)^\alpha \frac{1}{2}\right]$$
(5)
* Otherwise, cell shape is important and the exact solution is

$`<y(L)>`$ $`=`$ $`k\left({\displaystyle \frac{\theta _{\mathrm{res}}}{L}}\right)^2\left[{\displaystyle \frac{1}{(\alpha 2)}}\left({\displaystyle \frac{\theta _0}{\theta _{\mathrm{res}}}}\right)^\alpha {\displaystyle \frac{1}{2}}\right]`$ (6)
$`+`$ $`\left({\displaystyle \frac{\theta _0}{L}}\right)^\alpha \left[{\displaystyle _{x_{\mathrm{min}}}^{\mathrm{}}}x^\alpha {\displaystyle \frac{dG_p}{dx}}𝑑x{\displaystyle \frac{k}{(\alpha 2)x_{\mathrm{min}}^{\alpha 2}}}\right]`$
where $`x_{\mathrm{min}}=\theta _{\mathrm{lim}}/L_{\mathrm{min}}`$, and the integral must be solved numerically.

Thus equation 6 shows that in general, $`<y>`$ has a variation with cell size of the form

$$<y(L)>=aL^2+bL^\alpha $$
(7)
where $`a`$ and $`b`$ are constants. The limited angular resolution $`\theta _{\mathrm{res}}`$ reduces the variance: the existence of an object in a cell limits the available space in which other objects can appear. This effect varies with cell size because it depends on the scale of the resolution relative to the cell size. If a survey has sharp enough angular resolution then the first term of equation 7 may be neglected (if $`\alpha <2`$), and $`<y>L^\alpha `$. This is the case for the FIRST survey but not for the NVSS; thus equation 9 of Magliocchetti et al. (1998) only contains the second term of our equation 7.

The following statistic is commonly defined to characterize the departure of the skewness from a Poisson distribution:

$$z=\frac{\mu _33\mu _2+2\overline{N}}{(\overline{N})^3}$$
(8)
This has the expectation value

$$<z>=\frac{_{\mathrm{cell}}_{\mathrm{cell}}_{\mathrm{cell}}W(\theta _{12},\theta _{13},\theta _{23})𝑑S_1𝑑S_2𝑑S_3}{S^3}$$

where $`W(\theta _{12},\theta _{13},\theta _{23})`$ is the three-point angular correlation function. The significance of the skewness is summarized in Section 3 above.

### 6.2 Effect of multiple-component sources on the counts-in-cells moments

The existence of multiple-component sources increases the moments of counts-in-cells. This is because the fraction of radio sources within a cell that are split into multiple components varies from cell to cell, which acts to broaden the probability distribution of counts-in-cells. The simplest model is to suppose that a fraction $`e`$ of the galaxies are double sources (i.e. two sources appearing at the same location in space) and a fraction $`f`$ are triple sources. It can be shown (Blake 2002) that the expected offsets in the variance and skewness statistics are

$$\mathrm{\Delta }y=\frac{1}{\sigma S}\left(\frac{2e+6f}{1+e+2f}\right)$$
(9)
$$\mathrm{\Delta }z=\frac{1}{(\sigma S)^2}\left(\frac{6f}{1+e+2f}\right)$$
(10)
where $`S`$ is the cell area and $`\sigma `$ is the surface density of all components (i.e. catalogue entries). Thus skewness is sensitive only to triple sources.

This simple model neglects the fact that the components of a radio galaxy have a range of non-zero separations, but this should only matter if the cell size is not much greater than the maximum component separation ($`0.1^{}`$). A more sophisticated treatment is to model the separation distribution by the effective (small-angle) $`w(\theta )`$ of Figure 8. We can compute the effect on the variance statistic $`y`$ using equation 7, thus the general expression for $`y(L)`$ can be modified to

$$<y(L)>=aL^2+bL^\alpha +cL^\beta $$

where $`c`$ is a constant and $`\alpha `$ and $`\beta `$ are respectively the slopes of the shallow (galaxy clustering) and steep (multiple component) $`w(\theta )`$ power-laws.

It is easy to show that this more sophisticated model reduces in the appropriate limit to the more simple treatment initially outlined. As component separations tend to zero, the slope $`\beta `$ of the effective $`w(\theta )`$ becomes large and we can use the “steep clustering” approximation of equation 5, which reproduces the dependence $`yL^2`$ of equation 9.

### 6.3 Error on the counts-in-cells moments

Variance and skewness measurements are subject to statistical error due to averaging over a finite number of cells $`N_c`$. Calculating the standard error on the statistics $`y`$ and $`z`$ in the case of a random (unclustered) distribution yields

$$\sigma _y=\sqrt{\frac{2}{N_c(\overline{N})^2}}$$
(11)
$$\sigma _z=\sqrt{\frac{6}{N_c(\overline{N})^3}}$$
(12)
The probability distribution of the clustered data does not depart greatly from a Poisson distribution ($`y1`$, $`z1`$), so that these expressions are very good approximations to the actual statistical errors.

### 6.4 Measuring the counts-in-cells moments from a real survey

To measure the counts-in-cells we used the simple technique of defining a grid of touching circular cells of diameter $`L`$ on the sky and counting the number of sources that fall inside each cell. This only utilises a fraction $`\pi /4`$ of the available area, but the cell shape is constant over the sky. We note that the methodology of counts-in-cells was revolutionized by Szapudi (1998) who showed that it was valid to throw a very large number of randomly-placed cells over the sky, heavily oversampling the survey area. We prefer the former, simpler approach for a first investigation aiming to show consistency with the $`w(\theta )`$ analysis.

Surveys do not encompass the whole sky: there are boundaries and masked regions. Hence some cells in the grid are partially filled, the $`i`$th cell having fraction of useful area $`f_i`$ (say). To determine $`f_i`$ for each cell we populated the sky with random points subject to the same boundaries and masks as the real survey. Counting the number of random points that fall in each cell accurately measures the useful area. We then boosted the data count in the $`i`$th cell by a factor $`1/f_i`$, unless $`f_i`$ was less than a threshold $`f_{\mathrm{rej}}=0.75`$ in which case we rejected the cell. To measure the factors $`f_i`$ accurately enough it is essential to average over a sufficiently large number of random realizations that statistical noise does not dominate. Let there be $`m`$ random sets, each with the same surface density as the survey. The lower limit on $`m`$ is determined by the following considerations:

* The correction of cell counts by factors $`1/f_i`$ creates extra variance in the counts-in-cells (because all cells are corrected, whether they are partially filled or not). This extra systematic variance must be much smaller than the statistical error of equation 11. This condition is equivalent to

$$mm_1=\sqrt{\frac{N_c}{2}}\left(1+\frac{1}{\overline{N}}\right)$$
* The cell areas must be determined precisely enough that a negligible fraction of “unspoilt” cells are rejected with $`f_i<f_{\mathrm{rej}}`$. This condition is equivalent to

$$mm_2=\frac{1}{\overline{N}(1f_{\mathrm{rej}})^2}$$

When evaluating the moments of the counts-in-cells distribution we assume that all cells are populated independently. This is not strictly true given that clustered sources have correlated positions, but the assumption should be a very good approximation if the cells are large enough. With this consideration in mind, we adopted a minimum cell size $`L_{\mathrm{min}}=0.3^{}`$ for our analysis; below this cell size the average number of sources contained in a cell drops below $`\overline{N}=1`$ at the relevant flux-density thresholds.

### 6.5 Measurement of the NVSS counts-in-cells moments

Using these methods, we measured the counts-in-cells variance statistic $`y(L)`$ from the NVSS for circular cells of diameters $`0.3^{}<L<10^{}`$ for flux-density thresholds $`S_{1.4\mathrm{GHz}}=10`$ mJy and 20 mJy. As described in Blake & Wall (2002), we masked out all NVSS catalogue entries within $`5^{}`$ of the Galactic plane and also within 22 additional masked regions around radio sources that appear in the NVSS catalogue as a large number of separate elliptical Gaussians. As discussed in Section 5.1, the difficulty in resolving faint doubles with separations just greater than the beam-width causes some variation in the effective survey angular resolution, as evidenced by the fall-off in $`w(\theta )`$ with decreasing $`\theta `$ in the lowest separation bins of Figures 8 and 9. To ensure a consistent value of $`\theta _{\mathrm{res}}`$ for the prediction of the counts-in-cells variance from $`w(\theta )`$, we ran a linking algorithm that combined together all pairs of sources with angular separations less than 1 arcmin, thus artificially establishing $`\theta _{\mathrm{res}}=1`$ arcmin. This left 289,981 NVSS catalogue entries above 20 mJy and 516,782 entries above 10 mJy.

The variance measurements at the two flux-density thresholds are plotted in Figure 10, with error bars from equation 11. We also plot the predictions generated from the double power-law model for $`w(\theta )`$ with best-fit coefficients at the appropriate flux-density thresholds as listed in Table 1. The counts-in-cells variance is visually consistent with the prediction of the measured angular correlation function, verifying the agreement of these two independent methods of quantifying angular structure in the NVSS. The $`\chi ^2`$ statistics between the data and the prediction are $`\chi _{\mathrm{red}}^2=0.69`$ at 20 mJy and $`\chi _{\mathrm{red}}^2=1.64`$ at 10 mJy. However, the value of $`\chi ^2`$ is not strictly meaningful because the variances for different cell sizes are not independent. Note that multiple-component sources produce an approximately constant offset in $`yL^2`$ in Figure 10. The variance statistic for the 10 mJy threshold falls below that for the 20 mJy measurement because this offset is proportional to $`1/\sigma `$, where $`\sigma `$ is the source surface density (see equation 9).

The angular correlation function model used to generate the predictions is a sum of two power-laws, representing respectively multiple-component sources ($`w(\theta )\theta ^{3.1}`$) and galaxy clustering ($`w(\theta )\theta ^{0.8}`$). To assess the contribution of each effect to the overall counts-in-cells variance we separately converted each of the two 20 mJy power-laws to a variance and plotted the result on Figure 10. For small cell sizes ($`L<1^{}`$) the extra variance produced by multiple-component sources dominates. For larger cell sizes the extra variance produced by galaxy clustering becomes increasingly important. This type of transition is expected, as the number of extra pairs at angle $`\theta `$ scales as $`w(\theta )\times 2\pi \theta d\theta `$, which varies as $`\theta ^{2.1}`$ for multiple-component sources and $`\theta ^{+0.2}`$ for pairs of galaxies.

By converting the measured $`w(\theta )`$ into an equivalent variance we have shown that the angular correlation function and counts-in-cells analyses are entirely consistent. Alternatively, we can make an independent determination of the angular correlation function parameters by finding the best fit to $`y(L)`$. We varied the amplitude and slope of the clustering power-law $`w(\theta )=A\theta ^\alpha `$, whilst keeping the multiple-component $`w(\theta )`$ parameters fixed at their best-fitting values from Table 1. Thus for each point in the $`(A,\alpha )`$ grid we obtained a model $`y(L)`$, which we compared with the measured $`y(L)`$ using the $`\chi ^2`$ statistic. The best-fitting clustering parameters were in good agreement with those derived from the original $`w(\theta )`$ analysis (see Figure 11).

Figure 12 displays measurements of the skewness statistic $`z(L)`$ (equation 8) from the NVSS at 20 mJy and 10 mJy. The result is a significant detection of skewness. However, the fact that $`z`$ scales roughly as $`L^4`$ suggests that the skewness is dominated by multiple-component sources, especially as the amplitude of the skewness scales with surface density in accordance with equation 10. On Figure 12 we plot the multiple-component skewness predictions at 10 mJy and 20 mJy assuming that 7 per cent of radio galaxies are doubles (see Blake & Wall 2002) and 1 per cent are triples (i.e. $`e=0.07`$, $`f=0.01`$). This fraction of triple sources is reasonable; analyzing NVSS sources into groups using link-length $`\theta _{\mathrm{link}}=5`$ arcmin (a reasonable approximation to the range of dominance of multiple-component sources, see Figure 8) produces groups of which 9.8 per cent are doubles and 1.5 per cent are triples. These skewness predictions due to multiple components produce a fairly good fit to the data. The NVSS hence provides no convincing evidence for cosmological skewness.

### 6.6 Comparison with the FIRST counts-in-cells variance

We also derived counts-in-cells for the FIRST radio survey, previously studied by Magliocchetti et al. (1998). Our results are not directly comparable to those of Magliocchetti et al. (1998) as we make no attempt to combine multiple-component sources. We analyzed the three FIRST sub-samples used in the angular correlation function study of Section 5.2. To ensure a consistent survey angular resolution we first ran the source-combining algorithm with link-length $`\theta _{\mathrm{res}}=0.003^{}`$.

The variance results for the three FIRST sub-samples are plotted in Figure 13. In all samples there is a strong contribution from multiple components ($`yL^2=`$ constant) which dominates at small cell sizes. The amplitude of this contribution varies as $`1/\sigma `$ where $`\sigma `$ is the source surface density (equation 9), which accounts for the overall shift between the 10 mJy and 2 mJy samples (between which $`\sigma `$ varies by a factor $`3.4`$). Angular resolution effects are unimportant for FIRST; hence there is no dip in the variance at small $`L`$. Galaxy clustering becomes important at higher $`L`$ ($`yL^\alpha `$). The difference between the two samples at 2 mJy arises from the surface density gradients present in the whole-sky sample. Gradients offset the counts-in-cells variance by $`\mathrm{\Delta }y=`$ constant (Section 4.2), as is apparent from the increasing difference between the plotted triangles and open circles in Figure 13 (which plots $`yL^2`$ against $`L`$). Comparing Figures 13 and 10, the greater angular resolution of FIRST with respect to NVSS leads to an increased abundance of multiple-component sources and thus a greater variance signal for a given flux-density threshold.

### 6.7 Comparison with previous work

Our results differ from those of Magliocchetti et al. (1998) in two respects:

1. Magliocchetti et al. (1998) reported a much steeper slope for the correlation function, $`\gamma =2.5\pm 0.1`$ (where $`w(\theta )\theta ^{1\gamma }`$) compared to our NVSS measurement $`\gamma =1.83\pm 0.05`$ (Table 1). Magliocchetti et al. adopted a combining algorithm for multiple components and assumed that, after the operation of this algorithm, any remaining pairs were independent radio galaxies. However, no set of criteria can unequivocally distinguish multiple components from independent galaxies, and a small number of residual multiple-component sources can have a dramatic effect on the clustering statistics on angular scales up to several arc-minutes (see Blake & Wall 2002). As demonstrated by Figure 8, multiple components are the dominant provider of close pairs in the NVSS up to $`\theta 0.1^{}`$, whereas the fraction of closely-separated pairs combined by Magliocchetti et al. becomes negligible by $`0.02^{}`$. Thus the steep slope found by Magliocchetti et al. may be a manifestation of the remaining multiple-component sources, and not of galaxy clustering. Furthermore, it is difficult to understand a steep slope $`\gamma =2.5`$ persisting to small angles, given that the $`\gamma 1.8`$ clustering law is obeyed by all other studied classes of galaxy including local optically-selected galaxies (both spirals and ellipticals, e.g. Loveday et al. 1995) and high-redshift QSOs (Croom et al. 2001).
2. Magliocchetti et al. (1998) reported a cosmological skewness. We suggest that this may also be produced by uncombined doubles. The relation $`zy^2`$ used by Magliocchetti et al. as evidence for the non-linear gravitational growth of perturbations (their equation 20) is also naturally produced in a model where multiple components dominate – it results from combining our equations 9 and 10. An alternative explanation lies in the different flux-density limits (10 mJy for NVSS versus 3 mJy for FIRST). The FIRST sample may contain a non-negligible fraction of low-redshift starburst galaxies which trace non-linear clustering and hence are a source of skewness.

## 7 Conclusions

We have quantified the angular clustering in the NVSS and FIRST radio surveys using two independent methods: the two-point angular correlation function and the variance on counts-in-cells. Our results may be summarized as follows:

1. The results of angular correlation function and counts-in-cells analyses of the surveys are entirely consistent.
2. The larger area and greater number of sources in the NVSS yield a much clearer description of the clustering imprint. The correlation function has two contributions: that due to multiple components of the same galaxy, dominant at $`\theta <0.1^{}`$, and that due to clustering between galaxies, which dominates at larger angles. A clear break in $`w(\theta )`$ is evident between these scales. Both of these contributions are needed to explain the observed variance on counts-in-cells.
3. The clustering part of the correlation function has a slope consistent with that measured in the optical regime, $`w(\theta )\theta ^{0.8}`$; this is confirmed by our counts-in-cells measurements.
4. Both the NVSS and FIRST surveys suffer from systematic fluctuations in source surface density at flux-density thresholds at which they purport to be complete.

Our work disagrees with some previous conclusions drawn from the FIRST survey:

1. We find a galaxy correlation slope consistent with that measured in the optical, $`\gamma 1.8`$, in contrast to somewhat steeper slopes reported in previous analyses. These steeper slopes may have been produced by residual multiple component radio sources.
2. The skewness reported by Magliocchetti et al. (1998) may also be due to these residual double sources.

This investigation has improved our understanding of the methodology of angular clustering analyses for large-scale radio surveys, of relevant observational effects present in such surveys, and of the derived structural parameters. There is now the basis to use these surveys to derive three-dimensional information on the very largest structural scales, adopting more powerful statistical methods in conjunction with the redshift databases to be provided by surveys such as 2dF and SDSS.

## acknowledgments

We thank Lance Miller and Steve Rawlings for helpful comments on earlier drafts of this paper. We especially acknowledge Manuela Magliocchetti for extensive and helpful discussions.


format: text

text:
# Double Bubbles in the 3-Torus

## 1. Introduction

Our Central Conjecture 2.1 states that the ten different types of two-volume enclosures pictured in Figure 1 comprise the complete set of surface area minimizing *double bubbles* in the flat cubic three-torus T3. Our numerical results, summarized in Figure 2, indicate the volumes for which we conjecture that each type of double bubble minimizes surface area. Our main theorem, Theorem 4.1, states that given any fixed ratio of volumes, for small volumes, the minimizer is the standard double bubble. This result applies to any smooth flat Riemannian manifold of dimension three or four with compact quotient by its isometry group.

The *double bubble problem* is a two-volume generalization of the famous isoperimetric problem. The isoperimetric problem seeks the least-area way to enclose a single region of prescribed volume. About 200 BC, Zenodorus argued that a circle is the least-perimeter enclosure of prescribed area in the plane (see ). In 1884, Schwarz  proved by symmetrization that a round sphere minimizes perimeter for a given volume in R3. Isoperimetric problems arise naturally in many areas of modern mathematics. Ros  provides a beautiful survey.

Soap bubble clusters seek the least-area (least-energy) way to enclose and separate several given volumes. Bubble clusters have served as models for engineers, architects, and material scientists (the chapter by Emmer in  is a nice survey of architectual applications, and the text by Weaire and Hutzler , an introduction to the physics of space-filling bubble clusters or *foams*, discusses numerous other applications).

### Existence and regularity

The surface area minimizing property of bubble clusters can be codified mathematically in various useful ways, using the rectifiable currents, varifolds, or $`(𝐌,ϵ,\delta )`$-minimal sets of geometric measure theory (see ). In three dimensions, mathematically idealized bubble clusters consist of constant-mean-curvature surfaces meeting smoothly in threes at $`120^{}`$ along smooth curves, which meet in fours at a fixed angle of approximately $`109^{}`$ (\[24, Theorems II.4, IV.5, IV.8\], or see \[13, Section 13.9\]).

### Recent results

The existence and regularity of solutions to the double bubble problem played a key role in the proof by Hutchings, Morgan, Ritoré, and Ros (, ; see \[13, Chapter 14\]) that the standard double bubble, the familiar shape consisting of three spherical caps meeting one another at 120 degree angles, provides the least-area way to separate two volumes in R3. The proof relied on a component bound that had been developed by Hutchings \[9, Theorem 4.2\], and a detailed stability argument to rule out the possible remaining candidates. Reichardt *et al.*  extended these results to R4.

Contemporaneously with the development of results presented here, Corneli, Holt, Leger, and Schoenfeld  produced a more or less complete solution to the T2 version of the problem. Their proofs also rely on regularity theory, which in two dimensions implies that bubbles are bounded by circular arcs. A variational bound on the number of components of bubble clusters in surfaces due to Wichiramala  provided considerable added simplifications. A sequence of proofs using the techniques of plane geometry then eliminated all but five candidates (of which only four are expected to appear). Examining the T2 candidates was helpful to us in our work on the T3 project.

### Bubbles in the three-torus

In comparison with R3 and T2, the double bubble problem in T3 appears to be more difficult. In the torus it is not possible to push through a component bound like Hutchings’, since a key step in his proof is to show that the double bubble has an axial symmetry. Nor is a variational bound after Wichiramala forthcoming, due to the additional topological complications in three dimensions. Until some new approach provides a component bound, there probably will be no definitive results. Indeed, the single bubble for the three-torus is not yet completely understood, although there are partial results. The smallest enclosure of half of the volume of the torus was shown by Barthe and Maurey \[2, Section 3\] to be given by two parallel two-tori. Morgan and Johnson \[15, Theorem 4.4\] show that the least-area enclosure of a small volume is a sphere. Spheres, tubes around geodesics, and pairs of parallel two-tori are shown to be the only types of area minimizing enclosures for most tori by the work of Ritoré and Ros (\[19, Theorem 4.2\], ).

### Theorem about small volumes

Theorem 4.1 states that any sequence of area minimizing double bubbles of decreasing volume and fixed volume ratio has a tail consisting of standard double bubbles. The central difficulty is to bound the curvature. This accomplished, we show that the bubble lies inside some small ball that lifts to R3, where a minimizer is known to be standard . The result extends to any flat 3- or 4-manifold with compact quotient by the isometry group.

Our proof goes roughly as follows. From the original sequence of double bubbles we generate a new sequence by rescaling the manifold at each stage so that one of the volumes is always equal to one. We can then apply compactness arguments and area estimates to the rescaled sequence to show that certain subsequences of translates have non-trivial limits. These limits are used to obtain a curvature bound on the original sequence. With such a bound, we can apply monotonicity to conclude that if the volumes are small, the double bubble is contained in a small ball. We conclude that it must be the same as the minimizer in R3 or R4, *i.e.* that it must be the standard double bubble by  or ).

### Plan of the Paper.

Section 2 reviews the methods leading to our Central Conjecture 2.1 and to Figures 1 and 2. Section 3 surveys some subconjectures. Section 4 focuses on the proof of Theorem 1 on small volumes. Section 5 shifts from the cubic to other tori and discusses other conjectures and candidates, including a “Hexagonal Honeycomb.”

### Acknowledgments.

This paper has its origins in a problem given by Frank Morgan to the audience of his two week course on Geometric Measure Theory and the Proof of the Double Bubble Conjecture, at the Clay Mathematics Institute (CMI) Summer School on the Global Theory of Minimal Surfaces, held at the Mathematical Sciences Research Institute (MSRI) in June and July of 2001.

We owe our deepest gratitude to Frank Morgan for introducing us to this problem and taking the time to guide us along as we have worked on it, by many helpful discussions and constructive comments. CMI and MSRI are to be thanked for making the Summer School possible, and for providing research and travel support to the authors. We should like to extend our thanks to Joel Hass, David Hoffman, Arthur Jaffe, Antonio Ros, Harold Rosenberg, Richard Schoen, and Michael Wolf, who organized the Summer School. We benefited tremendously from the outstanding lectures and informal discussions, for which more people are owed thanks than we can name here. Three faculty participants, Michael Dorff, Denise Halverson, and Gary Lawlor, of Brigham Young University, are to be singled out and thanked for helping us get our computations started. We thank Manuel Ritoré and John Sullivan for helpful conversations. Other participants who made notable contributions include Baris Coskunuzer, Brian Dean, Dave Futer, Tom Fleming, Jim Hoffman, Jon Hofmann, Paul Holt, Matt Kudzin, Jesse Ratzkin, Eric Schoenfeld, and Jean Steiner.

Yair Kagan of New College of Florida created the picture of the Hexagonal Honeycomb in Figure 4. The authors would like thank him for this contribution.

Corneli participated in the 2001 “SMALL” undergraduate summer research Geometry Group, supported by Williams College and by the National Science Foundation through an REU grant to Williams and an individual grant to Morgan.

## 2. The Conjecture

### Generating Candidates.

Many possibilities for double bubbles in the three-torus T3 were proposed in brainstorming sessions by participants in the Clay/MSRI Summer School. In order to classify the candidates we used the following method. Starting with a standard double bubble, we imagined one of the two volumes growing until the bubble enclosing it wrapped around the torus and encountered an obstruction. Following the principals of regularity for bubble clusters, if a bubble collided with itself we opened the walls up, whereas if two different bubbles collided we allowed them to stick together. We then repeated this procedure for these new double bubbles, sometimes changing the perspective slightly slightly or becoming a bit more fanciful (e.g. the Center Cylinder of Figure 3 or Gary Lawlor’s Fire Hydrant of Figure 3, also known as “Scary Gary”).

### Producing the Phase Diagram.

Brakke’s *Surface Evolver*  was used to closely approximate the minimal area that a double bubble of each type needs to enclose specified volumes. Our initial simulations gave us the approximate surface area for each candidate double bubble, tested on partitions $`v_1:v_2:v_3`$ of a unit volume taken in increments of $`0.01`$. From the data obtained in these simulations, we found the least-area competitor for each volume triple. Figure 1 shows the candidates we found to be minimizing for some set of volumes. The phase diagram appearing in Figure 2 is the result of refining our initial computations along the boundaries with a 0.005 increment.

###### Central Conjecture 2.1.

The ten double bubbles pictured in Figure 1 represent each type of surface area minimizing two-volume enclosure in a flat, cubic three-torus, and these types are minimizing for the volumes illustrated in Figure 2.

### Comments.

One might expect that minimizers would be found among the various regularity-satisfying conglomerations of topological spheres and products of spheres and homotopically non-trivial tori, other possibilities being excessively complex. This was borne out in our computations. It is interesting to note, however, that not all of the simple possibilities along appeared as minimizers, for example the Transverse Cylinders pictured in Figure 3. The various double bubbles of Figure 3, while stable for a certain range of volumes, are never area minimizing. A challenging unsolved problem is to find all of the stable non-minimizing bubbles. Note that a given type from Figure 1 might be stable for a much wider range of volumes than those for which it actually minimizes surface area.

It is worth observing that all of the conjectured minimizers for the double bubble problem on T2 (a Standard Double Bubble, Band with Lens, Symmetric Chain, and Double Band) are echoed here in at least two ways. The Double Cylinder, the Slab Cylinder, the Double Slab, and the Cylinder String are T2 minimizers $`\times `$ T1. There are also more direct analogues, as is seen by comparing, for example, the three- and two-dimensional Standard Double Bubbles, or the Délauney and Symmetric Chains. See Corneli *et al.*  for more on the $`T^2`$ minimizers.

## 3. Subconjectures

We now present a list of natural subconjectures, suggested either by the phase diagram or by examination of the pictures made with *Surface Evolver*.

One immediate observation is that the edges of our phase diagram appear to characterize single bubbles in T3.

###### Conjecture 3.1 (Ritoré and Ros , , ).

The optima for the isoperimetric problem in a cubic T3 are the sphere, cylinder, and slab.

We still do not know proofs for the following intuitive conjectures about the double bubble problem (although Theorem 4.1 gives partial results on Conjecture 3.3):

###### Conjecture 3.2.

An area-minimizing double bubble in T3 has connected regions and complement.

###### Conjecture 3.3.

For two small volumes the standard double bubble is optimal.

###### Conjecture 3.4.

For one very small volume and two moderate volumes, the Slab Lens is optimal.

###### Conjecture 3.5.

The first phase transition as small equal, or close to equal, volumes grow is from the standard double bubble to a chain of two bubbles bounded by Délaunay surfaces.

Délaunay surfaces are constant-curvature surfaces of rotation, and as such have full rotational symmetry. From the *Surface Evolver* pictures, it appears that the conjectured surface area minimizers always have the maximial symmetry, given the constraints , a fact which leads to the following natural conjecture.

###### Conjecture 3.6.

A minimizer is as symmetric as possible, given its topological type.

We will conclude this section with a proposition which represents the first step towards establishing a symmetry property. Specifically, we prove that for any double bubble, there is a pair of parallel planes (actually two-tori) that cut both regions in half. Hutchings \[9, Theorem 2.6\] used the fact that in Rn a double bubble has two perpendicular planes that divide both regions in half. This is an easy step in a difficult proof that the function that gives the least-area to enclose two given volumes is concave. Before proving a basic result of a similarly elementary flavor for the torus, we mention that we conjecture that the much deeper concavity result also holds:

###### Conjecture 3.7.

The least area to enclose and separate two given volumes in the three-torus is a concave function of the volumes.

If one could prove Conjecture 3.7, one would then be able to apply other ideas in Hutchings’ paper \[9, Section 4\] to obtain a functional bound on the number of components of a minimizing bubble.

###### Proposition 3.8 (Deluxe Ham Sandwich Theorem***Name suggested by Eric Schoenfeld, who gave an independent proof for $`T^2`$.).

In a rectangular torus, if a double bubble lies inside a cylinder $`S^1\times D^2`$, then there is a plane that cuts both volumes in half.

###### Proof.

This is a generalization of the standard argument for the ham sandwich theorem in Euclidean space. Assume that the identification occurs in the vertical direction. Take a circle in the $`xy`$ plane indexed by $`\alpha [0,\pi ]`$. Then for each $`\alpha `$, rotate the surface by $`\alpha `$ in the $`xy`$ plane and consider the family of pairs of flat tori parallel to the $`xz`$ plane that are at distance 1/2 from each other. Then there is at least one such pair of planes that cuts the volume $`V_1`$ in half, for each $`\alpha `$. There may be an interval of such pairs for a given $`\alpha `$. However, some one parameter family of these planes may be chosen which varies continuously as a function of $`\alpha `$, and one such plane cuts $`V_2`$ in half. ∎

## 4. Small Volumes

Conjecture 3.3 stated that small volumes are best enclosed by a standard double bubble. Theorem 4.1 proves that for any fixed volume ratio, the standard double becomes optimal when the volumes are sufficiently small. This result holds for many flat three- and four-dimensional manifolds (and see Remark 4.4).

The standard double bubble consists of three spherical caps meeting at 120 degrees. (If the volumes are equal, the middle surface is planar.) It is known to be minimizing in R3  and R4 .

###### Theorem 4.1.

Let $`M`$ be a smooth flat Riemannian manifold of dimension three or four, such that $`M`$ has compact quotient by the isometry group. Fix $`\lambda (0,1]`$. Then there is an $`ϵ>0`$, such that if $`0<v<ϵ`$, a perimeter-minimizing double bubble of volumes $`v,\lambda v`$ is standard.

###### Remark 4.2.

Solutions to the double bubble problem exist for all volume pairs in manifolds with compact quotient by their isometry group. The proof is the same as in Morgan \[13, Section 13.7\].

For the proof we will regard a double bubble as a pair of $`3`$\- or $`4`$dimensional rectifiable currents, $`R_1`$ and $`R_2`$, each of multiplicity one, of volumes $`V_1=𝐌(R_1)`$ and $`V_2=𝐌(R_2)`$. The total perimeter of such a double bubble is $`\frac{1}{2}(𝐌(R_1)+𝐌(R_2)+𝐌((R_1+R_2)))`$. Here $`𝐌`$ denotes the mass of the current, which can be thought of as the Hausdorff measure of the associated rectifiable set (counting multiplicities). For a review of the pertinent definitions, see notes from Morgan’s course at the Summer School  or the texts by Morgan  or Federer .

Before proceeding it is helpful to fix some notation. By the Nash and subsequent embedding theorems, we may assume $`M`$ is a submanifold of some fixed RN. We will consider a sequence of perimeter-minimizing double bubbles in $`M`$, containing the volumes $`v`$ and $`\lambda v`$, as $`v0`$. For each $`v`$, $`M_v`$ will denote $`s_v(M)`$ in $`s_v(𝐑^N)`$, where $`s_v`$ is the scaling map that takes regions with volume $`v`$ to similar regions with volume $`1`$. In particular, $`s_v`$ maps our perimeter-minimizing double bubble containing volumes $`v`$ and $`\lambda v`$ to a double bubble which we call $`S_v`$ which contains volumes $`1`$ and $`\lambda `$, and which is of course perimeter-minimizing for these volumes.

We focus on the case of dimension three; the proof for dimension four is essentially identical.

###### Lemma 4.3.

There is a $`\gamma >0`$ such that if $`R`$ is a region in an open Euclidean 3-cube $`K`$ and $`\mathrm{vol}(R)\mathrm{vol}(K)/2`$, then

$$\mathrm{area}(R)\gamma (\mathrm{vol}(R))^{2/3}.$$

###### Proof.

Let $`\gamma _0`$ be such an isoperimetric constant for a cubical 3-torus, so that $`\mathrm{area}(P)\gamma _0(\mathrm{vol}(P))^{2/3}`$ for all regions $`P𝐓^3`$. Such a $`\gamma _0`$ exists by the isoperimetric inequality for compact manifolds \[13, Section 12.3\]. Make the necessary reflections and identifications of the cube $`K`$ to obtain a torus containing a region $`R^{}`$ with eight times the volume and eight times the boundary area of $`R`$. The claim follows, with $`\gamma =\gamma _0/2`$. ∎

###### Proof of Theorem 4.1.

The first step is to show that the sequence $`S_v`$, suitably translated and rotated, has a subsequence that converges as $`v0`$ and has $`V_10`$ in the limit. Our argument also shows that there is a subsequence that converges to a limit with $`V_20`$, but does not show that there is a subsequence where both volumes are non-zero in the limit. This is because while we are exerting ourselves trapping the first volume in a ball, the second one may wander off to infinity.

We first show the existence of a covering $`𝒦_v`$ of $`M_v`$ with bounded multiplicity, consisting of cubes contained in $`M_v`$, each of side-length $`L`$. Lemma 1 will give us a positive lower bound on the volume of the part of $`R_1`$ that is inside one of these cubes for each $`S_v`$. We will then apply a standard compactness theorem to show that a subsequence of the $`S_v`$’s converges.

Take a maximal packing by balls of radius $`\frac{1}{4}L`$. Enlargements of radius $`\frac{1}{2}L`$ cover $`M_v`$. Circumscribed cubes of edge-length $`L`$ provide the desired covering $`𝒦_v`$. To see that the multiplicity of this covering is bounded, consider a point $`pM_v`$. The ball centered at $`p`$ with radius 2$`L`$ contains all the cubes that might cover it, and the number of balls of radius $`\frac{1}{4}L`$ that can pack into this ball is bounded, implying that the multiplicity of $`𝒦_v`$ is also bounded by some $`m>0`$.

Now consider some such covering, with $`L=2`$. By Lemma 1, there is an isoperimetric constant $`\gamma `$ such that

$$\mathrm{area}((R_1K_i))\gamma (\mathrm{vol}(R_1K_i))^{2/3},$$

and therefore, since $`\mathrm{max}_k\mathrm{vol}(R_1K_k)\mathrm{vol}(R_1K_i)`$ for any $`i`$,

(4.1)
$$\mathrm{area}((R_1K_i))\gamma \frac{\mathrm{vol}(R_1K_i)}{(\mathrm{max}_k\mathrm{vol}(R_1K_k))^{1/3}}.$$

Note that the total area of the surface is greater than $`1/m`$ times the sum of the areas in each cube, and the total volume enclosed is less than the sum of the volumes, so summing Equation 4.1 over all the cubes $`K_i`$ in the covering $`𝒦_v`$ yields

$$\mathrm{area}(S_v)\mathrm{area}(R_1)m\gamma \frac{V_1}{(\mathrm{max}_k\mathrm{vol}(R_1K_k))^{1/3}}$$

and

$$(\underset{k}{\mathrm{max}}(\mathrm{vol}(R_1K_k)))^{1/3}m\gamma \frac{V_1}{\mathrm{area}(S_v)}\delta >0,$$

because $`V_1=1`$ and it is easy to show that $`\mathrm{area}(S_v)`$ is bounded (since there is a bounded way of enclosing the volumes). Translate each $`M_v`$ so that the cube where the maximum occurs is centered at the origin of RN, and rotate so that the tangent space of each $`M_v`$ at the origin is equal to a fixed R3 in RN. The limit of the $`M_v`$ will be equal to this R3. Since a cube with edge-length $`L`$ centered at the origin fits inside a ball of radius $`2L`$ centered at the origin, we have

$$\mathrm{vol}(R_1B(0,2L))\delta ^3$$

for every $`S_v`$. By the compactness theorem for locally integral currents (\[13, pp. 64,88\]\[23, Section 27.3, 31.2, 31.3\]) we know that a subsequence of the $`S_v`$ has a limit, which we will call $`D`$, with the property that $`\mathrm{vol}(R_1)\delta ^3`$. This completes the first step.

Since $`D`$ is contained in the limit of the $`M_v`$, namely, the copy of R3 chosen above, and each $`S_v`$ is minimizing for its volumes, a standard argument shows that the limit $`D`$ is the perimeter-minimizing way to enclose and separate the given volumes in R3 (cf. \[13, 13.7\]). In the limit, $`V_2`$ could be zero, in which case D is a round sphere. If both volumes are non-zero, D is the standard double bubble ( and , or see \[13, Chapter 14\]).

Our goal is now to prove that the double bubble of volumes $`v`$, $`\lambda v`$ lies inside a trivial ball, when $`v`$ is small enough. This can be accomplished using the monotonicity theorem for mass ratio \[1, Section 5.1(1)\], which implies that for a perimeter-minimizing bubble cluster, a small ball around any point on the surface contains some substantial amount of area. This will limit the number of disjoint balls we can place on the surface. The monotonicity theorem applies only to surfaces for which the mean curvature is bounded, *i.e.*, for which there is a $`C`$ such that for smooth variations

$$\frac{dA}{dV}C.$$

Accordingly, the second major step in the proof is to obtain such bound on the curvature, for $`v`$ small. It suffices to show that all smooth variation vector fields have the property that changes in the volume of $`S_v`$ and in the area of $`S_v`$ are controlled. Take a smooth variation vector field $`F`$ in Rn such that for $`D`$,

$$dV_1/dt=_{R1}(Fn)𝑑A=c>0$$

and

$$dV_2/dt=0.$$

(Note that we need the first volume to be non-zero, or its variation could be zero.) For $`v`$ small enough the subsequence of $`S_v`$ headed towards the limit $`D`$ has the property that $`dV_1/dt`$ is approximately $`c`$ and $`dV_2/dt`$ is approximately 0.

By the argument in the first step, we can translate each $`S_v`$ similarly, so that a subsequence of the subsequence above converges to a minimizer $`D^{}`$ in R3 where the second volume is non-trivial. This time we take a smooth variation vector field $`F^{}`$, such that for $`v`$ small enough the subsequence of $`S_v`$ headed to $`D^{}`$ has the property that $`dV_1/dt`$ is approximately 0 and $`dV_2/dt`$ is approximately $`c^{}>0`$. This proves that for this subsequence the change in volume is bounded below.

Now we need to show that the change in area is bounded above. This follows from the fact that every rectifiable set can be thought of as a varifold \[13, Section 11.2\]. By compactness for varifolds \[1, Section 6\], the $`S_v`$, situated so that the first volume does not disappear, converge as varifolds to some varifold, $`J`$. The first variation of the varifolds also converge, *i.e.*, $`\delta S_v\delta J`$, see . The first variation of a varifold is a function representing the change in area. Therefore, far enough out in the sequence the change in area of the $`S_v`$ under $`F`$ is bounded close to the change in area of $`J`$ under $`F`$, which is finite. Similarly, the change in area of the $`S_v`$ under $`F^{}`$ is bounded.

We conclude that the mean curvature of $`S_v`$ is bounded for two independent directions in the two-dimensional space of volume variations, and hence for all variations. This completes the second step.

The third and final step is to show that all of the surface area is contained in some ball in $`M_v`$, of fixed radius for all $`v`$. Eventually, as $`v`$ shrinks and $`M_v`$ grows, this ball will have to be trivial in $`M_v`$. We will then use the result that the optimal double bubble in R3 is standard to show that our double bubble is standard as well.

As the $`S_v`$ are approaching the minimizer in R3, there must be a bound $`A`$ on the perimeters of the $`S_v`$. By monotonicity of mass ratio \[1, Section 5.1(1)\], every unit ball centered at a point of $`S_v`$ contains perimeter $`\delta >0`$. Therefore there are at most $`A/\delta `$ such disjoint balls.

We claim that the $`S_v`$ are eventually connected. There is an upper bound on the diameter of any component, $`2A/\delta `$. Since we are controlling curvature, our components cannot become too small. Since every unit ball contains at least $`\delta `$ area, we also have a lower bound on the area of each component, when unit balls are trivial. Unless eventually the $`S_v`$ are connected, you can arrange to get in the limit a disconnected minimizer in R3, a contradiction.

Hence $`S_v`$ is contained in a ball of radius $`2A/\delta `$ for all $`v`$. Since our original manifold has compact quotient by its isometry group, there is a radius such that balls in the original manifold of that radius or smaller are topologically trivial. Hence, as we expand the manifold, eventually balls of radius $`2A/\delta `$ can be lifted to R3, which means that they are Euclidean. Hence, $`S_v`$ is eventually contained in a Euclidean ball, and is therefore the standard double bubble (, , \[13, Chapter 14\]).

Finally, since $`S_vM_v`$ is simply a scaled version of the original double bubble in $`M`$, we conclude that the original double bubble is standard as desired. ∎

###### Remark 4.4.

Given $`n`$, $`m`$, similar arguments show that for any smooth $`n`$-dimensional Riemannian manifold with compact quotient by the isometry group, given $`0\lambda 1`$, there are $`C`$, $`ϵ>0`$, such that for any $`0<v<ϵ`$, a minimizing cluster with $`m`$ prescribed volumes between $`\lambda v`$ and $`v`$ lies inside a ball of diameter at most $`Cv^{1/n}`$.

To further deduce that the cluster smoothly approximates a Euclidean minimizer would require knowing that convergence weakly and in measure, under bounded mean curvature, implies $`C^1`$ convergence, as is known for hypersurfaces without singularities (\[1, Section 8\], see \[14, Section 1.2\])

## 5. Special Tori

Changing the shape of the torus, by stretching it or by skewing some or all of its angles, would certainly change the phase diagram of Figure 2:

###### Conjecture 5.1.

In the special case of a very long $`T^3`$ the Double Slab is optimal for most volumes.

Special tori may have special minimizers:

###### Conjecture 5.2.

For the special case of a torus based on a relatively short $`\frac{\pi }{3}`$-rhombic right prism, the Hexagonal Honeycomb prism of Figure 4 is a perimeter-minimizing double bubble for which both regions and the exterior each have one third the volume, or when two volumes are equal and the third is close.

Indeed, for such volumes the Hexagonal Honeycomb ties the Double Slab, just as in the $`\frac{\pi }{3}`$-rhombic two-torus a Hexagonal Tiling ties the Double Band .

In the *triple bubble problem* for the Face Centered Cubic (FCC) and Body Centered Cubic (BCC) tori we would expect to find minimizers that lift to R3 as periodic foams with cells of finite volume. The Weaire-Phelan foam, a counterexample to Kelvin’s conjectured best way to divide space into unit volumes (, ; see \[13, Chapter 15\], ), would be expected to appear as lifts of some solutions to the triple bubble problem in the BCC torus. Kelvin’s foam might appear as a solution to the triple bubble problem in the FCC torus. One might also see *eight* Kelvin cells in the BCC torus, or *sixteen* Kelvin cells in the standard cubic torus. Eight Weaire-Phelan cells also fit in the standard cubic torus. Furthermore, by scaling, higher numbers of cells fit into these tori.

In contrast with the triple bubble problem and problems with more volume constraints, it is extremely unlikely that a solution to the double bubble problem in any torus would lift to a division of R3 into finite volumes, since by regularity, singular curves meet in fours. This means that singular points look locally like the cone over a tetrahedral frame, and hence have four volumes coming together. Since a region will never be adjacent to another component of the same region (because the dividing wall could be removed to decrease area and maintain volumes), a foam generated using a fundamental domain coming from a double double in the three-torus would have to exhibit the strange property that the singular curves never meet.

###### Conjecture 5.3.

There are no least-area divisions of T3 into three volumes that lift to a foam in R3.

This conjecture suggests that it is not likely that there will be are any other special minimizers for the double bubble problem.

###### Conjecture 5.4.

The double bubbles of Figure 1 together with the Hexagonal Honeycomb of Figure 4 comprise the complete set of area minimizing double bubbles for all three tori.

As a final remark, in light of the fact that the triple bubble problem in the torus seems likely to produce so many interesting candidates, we would like to mention one final conjecture.

###### Conjecture 5.5.

For the triple bubble problem in a cubic T3 in the case where one of the volumes is small, the minimizers will look like the double bubbles of Figure 2 with a small ball attached. The phase diagram will look just like our Figure 3.