The histogram of a sky too crowded to count

A survey reaches its confusion limit when the sources it is trying to count become so numerous that there is more than about one of them for every thirty beams. Past that point a peak in the map is a blend rather than a galaxy, the census acquires a knee it did not have at exactly the wrong flux, and no amount of further exposure separates the blends, because the thing blurring them is the telescope’s own diffraction and the sky itself, not noise.

That limit is a wall for anyone trying to make a catalogue. It is not a wall for anyone trying to measure how many sources there are. A confused map is still a map of a sum of discrete objects, and a sum of a random number of discrete objects has a distribution whose shape depends on how many there are and how bright. The histogram of the pixel values in a confused map is a measurement of the source counts below the confusion limit, and it was the first method ever devised for getting a number out of a sky that could not be counted.

It was devised by a graduate student, in 1957, to rescue a cosmological argument from the instrument that had produced it.

A student and a misleading map #

The second Cambridge radio survey had counted radio sources to a depth at which its beam was badly confused, and the counts it published rose far too steeply towards faint fluxes — steeply enough to argue that the universe had evolved and the steady-state theory was wrong. The Sydney group, surveying the same sky with a different instrument, could not find most of the faint Cambridge sources at all. The steep slope was a property of the beam: blends of faint sources had been counted as brighter ones.

Peter Scheuer’s thesis proposed not counting at all. Take the record of the telescope’s output — for a radio interferometer of that era, a chart-recorder trace, for a modern instrument a map — and form the histogram of its values. Call the value in one beam the deflection $D$. Every source on the sky contributes to $D$ in proportion to its flux times the beam’s response at its position, so $D$ is a sum over a random, Poisson-distributed set of sources with random fluxes and random positions. Its probability distribution, $P(D)$, can be computed from the counts $dN/dS$ and the beam shape alone. Fit the computed $P(D)$ to the measured histogram, and the counts come out.

Scheuer’s analysis gave a slope consistent with the Sydney counts and inconsistent with the steep Cambridge one, at fluxes where no source could be individually identified in either survey. The method sat largely unused for a decade, then became standard in radio astronomy, and it is now how the counts below the confusion limit are measured in every band where maps are confused, including the maps synthesised by radio interferometers whose beams are measured to a part in a thousand for exactly this purpose.

How a histogram is computed from a count #

The calculation Scheuer set up is worth stating, because its structure explains why the method works and where it will fail.

Each source contributes to a pixel its flux times the beam’s response at its offset, a number between zero and its full flux. Over the whole sky, what matters is therefore not the counts themselves but the rate $R(x)$ at which contributions of size $x$ occur in one pixel: the counts folded through the beam, so that a bright source far from the pixel’s centre contributes as a faint source at the centre would. The deflection $D$ is a sum of a Poisson-distributed number of such contributions, which makes it what probability theory calls a compound Poisson variable, and those have one clean property. The Fourier transform of $P(D)$ is

$$\tilde P(\omega) = \exp\!\left(\int_0^\infty R(x)\,\left(e^{i\omega x} - 1\right)dx\right).$$

Two consequences follow at once. The counts and the beam enter only through $R(x)$, so two different skies seen through two different beams can give the same histogram if their folded rates agree — which is why the beam must be measured independently, and why the Cambridge beam’s sidelobes were fatal. And the logarithm of the transform is linear in the counts: the histogram of a map with two populations is the convolution of the two histograms each would give alone. That is what makes it possible to ask, as the last sections here do, what a hypothetical extra population would do to a histogram already explained without it — it convolves, and a convolution can only widen.

The same structure is behind the error bar that comes from counting photons, which is the special case in which every contribution is the same size. Photon noise is the $P(D)$ of a sky of identical sources; a confused map is photon noise in which the photons come in a distribution of sizes, and the distribution is the thing being measured.

Why the shape carries the counts #

A pixel in a confused map receives light from many sources, most of them faint and a few of them near the confusion limit. Different parts of the count distribution leave different marks on the histogram, and seeing which marks is most of understanding the method.

The long tail to the right is made by the brightest sources in the map, the ones near and above the confusion limit. It is the most visible feature of the histogram and the least informative about the faint counts: it is shared by all three maps, because they share their bright counts, and it carries the sample variance of the few dozen bright objects that happen to fall in the field.

The informative part is the left side. A pixel’s lowest possible value is set by the beams that happen to contain the fewest sources — and since there are sources everywhere, the minimum is not zero but a floor set by the faint population that fills every beam. If faint sources are rare (a shallow slope, $b = 1.3$), some beams are nearly empty and the histogram rises abruptly from a sharp lower edge. If they are plentiful ($b = 1.7$), every beam holds several of them, the emptiest beams are less empty, and the lower edge is smeared. The faint counts are read off the side of the histogram where the sky is darkest, which is the side any glance at a map passes over.

The noise that exposure cannot remove #

The scatter between the pixels of a confused map has two parts, and only one of them responds to observing time.

The flat part of that curve is what an observer of a confused field eventually sees in their data: pixel scatter that stops improving. It used to be called confusion noise, and the name is misleading in a way that matters here. Noise is something to be reduced. The confusion “noise” is the signal the $P(D)$ method measures, and a survey that has integrated well past the crossing point has not wasted its time — it has measured the histogram’s shape precisely enough to fit.

That is the one sense in which a longer exposure of a confused field does buy depth. It does not buy fainter sources; it buys a sharper histogram, and the histogram’s shape is what carries the faint counts. The trade is explicit in how the Herschel deep fields were used. The instrument’s noise in the deepest maps was pushed well below the sky’s scatter, and the resulting histograms were fitted for counts down to a few millijanskys at 250 microns — roughly a decade below the flux at which a source could be individually trusted.

What the far-infrared histograms found #

The deepest far-infrared histograms were used to answer a question the individual detections could not: how much of the infrared background is made by galaxies of the kind already known, and at what flux.

The Herschel surveys detected individually only about a sixth of the 250-micron background. Fitting the histograms of the same maps extended the counts about a decade fainter and accounted for well over half of it — and the counts recovered from the histograms turned over below the confusion limit in the way the shallow slopes above predict, which is what makes the background converge rather than diverge at the faint end. What remains unaccounted for after that is within the uncertainty of the background’s own absolute level, which was measured against a foreground a hundred times brighter and is not known to better than some tens of per cent.

The physical content of that result is a statement about where the universe’s dust-obscured star formation happened. The sources making up the fitted counts are galaxies at redshifts of one to three, forming stars at tens to hundreds of solar masses a year behind enough dust that the ultraviolet indicators of star formation miss most of it. The far-infrared background carries about as much energy as all the optical and ultraviolet starlight put together, so roughly half of all the star formation there has ever been is recorded in it — and most of that half was measured, not by counting galaxies, but by fitting the histogram of the pixels between them.

Where the instrument’s noise does hide something #

Instrumental noise is Gaussian and it is added to every pixel, which means it smears the histogram’s shape, and the part of the shape it smears first is exactly the informative part.

That is why the method needs instrumental noise well below the confusion width, and why it works best on the deepest maps rather than on the widest. It is also why the skewness, rather than the width, is the statistic that carries the faint counts through a noisy measurement: a Gaussian has zero skewness whatever its width, so every bit of asymmetry in a histogram is a property of the sources.

Masking the sources that can be found #

The bright sources are both the least informative part of the histogram and the part with the most sample variance, and a real analysis removes them before fitting.

The iteration is a small version of a genuine circularity in defining the confusion limit. The limit is usually quoted as some multiple of the confusion noise, but the confusion noise depends on which sources have been counted as detected and removed, which depends on the limit. For steep counts the iteration converges quickly and the answer is stable; for shallow ones it can move by a factor of two, and published confusion limits for the same instrument differ for exactly that reason. The number quoted as “the” confusion noise of a far-infrared camera is a number that already has a detection threshold inside it.

A slope recovered from nothing that was counted #

With the bright sources masked and the noise known, the fit itself is a comparison between a measured histogram and histograms computed — or, as here, simulated — for a family of possible counts.

The asymmetry in that curve is the most important thing the figure shows, and it is not an accident of the simulation. Shallow faint counts are easy to rule out and steep ones are hard. A shallow slope means few faint sources, and few sources per beam make a histogram with distinctive, lumpy, non-Gaussian structure that a wrong model cannot imitate. A steep slope means many faint sources per beam, and a sum of many independent contributions tends to a Gaussian whatever the individual contributions are — the central limit theorem. Past a certain depth, the faint population contributes only a mean (which a real map does not measure, having had its mean removed) and a variance (which trades against every other contribution to the width). Its shape has been erased by its own abundance.

So the $P(D)$ method has a floor of its own, and the central limit theorem sets it. The histogram carries the counts down to the depth at which there are a few sources per beam, and not much further. For the Herschel maps that depth was a few millijanskys; for the deepest radio maps it is below a microjansky. In both cases it is about a decade below the conventional confusion limit — a real gain, and a gain whose end is set by arithmetic rather than by any instrument.

A background the counts do not explain #

The method’s sharpest use is not measuring counts that were expected but ruling out counts that would have been needed.

In 2009 the ARCADE 2 balloon experiment measured the absolute temperature of the sky at centimetre wavelengths and found, above the microwave background’s blackbody, a radio background several times brighter than the sum of all the radio sources ever counted. That sum is well determined: it converges, because the counts are shallow at faint fluxes, and the counted sources account for perhaps a fifth or a quarter of what ARCADE 2 saw after the Milky Way’s own emission was modelled and removed. The excess agreed with older absolute measurements at lower frequencies, which had been harder to interpret.

If the excess is made of discrete sources, they must lie below the flux at which the radio surveys can count them. The histogram can say how far below.

This is the argument Condon and colleagues made in 2012 with a deep map at 3 GHz, whose histogram was fully explained by extrapolating the known counts. Any population making the excess had to be faint enough to vanish into the Gaussian core of that histogram, which forced it below about a microjansky, and numerous enough to supply the excess from there, which forced it to outnumber by a wide margin every galaxy known to exist in the same patch of sky. No known class of galaxy is that numerous and that faint in the radio, and a new one would have to be invisible at every other wavelength too.

The conclusion drawn from this is not that the excess is wrong. It is that if it is real it is probably not made of sources at all — or not of sources in the sense that the counts describe. The candidates then are diffuse: synchrotron emission from a large, faint halo around the Milky Way that the foreground models mis-subtract, emission from the turbulent gas in galaxy clusters and the filaments between them, or something more exotic. None has been established, and the radio synchrotron background remains the one background in the sky whose level and whose counted sources disagree by a factor of several.

What a histogram can and cannot know #

The simulated maps have no clustering. Real galaxies cluster, so the number of sources in a beam is not Poisson-distributed but over-dispersed, and the histogram of a clustered sky is wider and more skewed than a Poisson one with the same counts. A $P(D)$ fit that ignores clustering reads the extra width as extra faint sources. Every serious far-infrared analysis has to model the clustering alongside the counts, and at 500 microns, where the beam is widest and the sources most strongly clustered, the correction is not small.

The sources are points. At the fluxes where confusion bites in the far infrared and the radio this is a good approximation, since the beams are much larger than the galaxies; in the optical it is not, and the optical “confusion” of the deepest images is a problem of overlapping galaxy profiles that the $P(D)$ formalism does not describe.

The mock fit uses the same kind of model that made the mock. Its success says that the histogram contains enough information to recover a slope when the model family is right; it does not say the family is right. Real counts are not two power laws, and a fit with the wrong family returns a wrong slope with a small error bar — the error bar measures how well the histogram constrains the parameters of the model, not whether the model describes the sky.

And the excess figures treat the hidden population as single-flux. A population with a spread of fluxes adds more variance for the same total, because its brightest members contribute disproportionately, so the single-flux case is the most forgiving one and the constraint drawn is the loosest.

An old paradox, turned into a test #

The histogram gives the confusion limit a second meaning, and it is the one that makes the connection back to Olbers’ paradox exact.

Olbers’ sky, where every sight line ends on a star, would have no histogram worth the name. Every beam would be filled edge to edge with stellar surface, every pixel would have the same value, and $P(D)$ would be a spike. A saturated sky carries no information about the number of sources in it, because the sources have become a wall. The real sky is the opposite limit: sparse enough that the number of sources in a beam fluctuates, and those fluctuations are why its counts can be measured below the depth at which it can be counted.

The central limit theorem is the bridge between the two. As the number of sources per beam rises, the histogram narrows relative to its mean and turns Gaussian, and the information about the individual sources drains out of it; at infinite depth and uniform density it would become Olbers’ spike. The depth at which a histogram stops carrying the counts is the depth at which a sky begins, statistically, to look like a wall — and it is set by exactly the quantity that resolved the paradox in the first place: how many sources a sight line meets before it runs out of universe.

Still open: a sky measured with no sources in it at all #

The histogram measures the counts of sources that are too crowded to separate. The obvious extension, now a whole programme, drops the separation entirely.

The fluctuations of a confused map have spatial structure as well as a histogram, and the structure carries the clustering of the sources — the angular power spectrum of the far-infrared background has been measured this way by Planck and Herschel, as a measurement of how the dusty galaxies that make it trace the dark matter. If the map is made in a narrow spectral line rather than in broad continuum, the fluctuations also carry a redshift, and the result is a three-dimensional map of line emission — carbon monoxide, ionised carbon, neutral hydrogen — made without detecting a single galaxy. Whether line-intensity mapping of this kind can measure the history of star formation and the growth of structure at redshifts beyond any galaxy survey, and how to separate the faint cosmological lines from foregrounds a thousand times brighter, is where the measurement of what the sky’s darkness is made of is currently being taken.