Astronomical Data & Computation: Counting Photons, Finding Rhythms, Building Pictures

The night mode on a phone takes a dozen dim, grainy frames and stacks them into one clean picture. Astronomers have been doing the same thing since the 1970s, when the charge-coupled device (invented in 1969 by Willard Boyle and George Smith, Nobel Prize 2009) replaced photographic plates. A faint galaxy is not in any single exposure. It emerges from the arithmetic, and the arithmetic tells you how sure to be.

Now scale it up. The European Space Agency's Gaia mission released a catalogue of about 1.8 billion stars and other sources in June 2022. The Vera C. Rubin Observatory in Chile began its ten-year survey in 2026 and can issue up to about 10 million alerts about changing objects in a single night. The 2019 picture of the black hole in M87 was assembled from roughly five petabytes of recordings that were flown to two supercomputers on hard drives. Nobody reads that by eye. Astronomy has become a data science, and this piece is what an astronomer's data education still looks like years later.

Every pixel is a count: noise, signal and stacking

A detector does not measure brightness. It counts electrons freed by photons. Photons arrive at random, so a source that gives on average 100 counts gives 100 plus or minus 10: a spread equal to the square root of the count (Poisson statistics). That is the physical floor, called shot noise. On top come the night sky's own glow, which also arrives as random photons, and the detector's read noise, a fixed few electrons added every time the chip is read out. Put together, for a star measured inside an aperture of n pixels:

Read it in two regimes. Bright star, dark sky: the S under the root dominates, and the signal-to-noise ratio grows as the square root of S, so four times the exposure gives twice the SNR. Faint star in a bright sky: the sky term dominates, and the result is the same square-root law. Either way, the lesson is that doubling the SNR costs four times the telescope time, which is why big mirrors and long exposures matter. The twist is read noise. It is paid per frame, so one 60-second exposure beats sixty 1-second ones, unless the sky is bright enough that sky noise swamps it anyway. Try both.

Without JavaScript: in a dark sky, one 60 s frame detects the faint star at SNR of about 8; sixty 1 s frames give about 1.5. In a city sky one 60 s frame gives about 3, ten of them stacked about 9.

This is also why astronomers calibrate. Before the counts mean anything they subtract the bias and dark current (the signal in a closed shutter) and divide by a flat field (a picture of an evenly lit screen, which maps the sensitivity of each pixel and the shadows of dust). Skip those and a smooth sky becomes a landscape of fake structure.

The periodogram: finding a rhythm in data with gaps

Stars pulse, planets pass in front of them, binaries eclipse. All are periodic, and the tool is the Fourier transform. But a telescope cannot watch for a month without stopping. The sun comes up, clouds roll in, the Moon is bright. The observations fall on unevenly spaced times, and the ordinary Fourier transform assumes an even grid. The answer is the Lomb–Scargle periodogram (Lomb 1976, Scargle 1982): it fits a sine wave at each trial frequency by least squares and plots how much of the variation each frequency explains, with no need for an even grid. VanderPlas (2018) gives a modern account.

The gaps leave a signature. If you observe at about the same time every night, you sample the star once a day, and a true period of P gives ghost peaks at the frequencies 1/P ± 1 per day, called aliases. An observer at one site cannot tell a 0.7-day rhythm from a 2.3-day one, except by other clues. The cure is to sample at other times, and that is why astronomers spread networks of telescopes around the world: when the star sets in Chile it is rising over Australia.

Without JavaScript: with one site, the periodogram of a 1.7-day signal (0.59 cycles/day) shows a second peak almost as high at 0.41 cycles/day (the alias at 1 − 0.59). With three sites the alias drops to under 4% of the main peak.

Fitting with error bars: the posterior is the result

A measurement without an uncertainty is not a measurement. Fitting a model to data should give a best value and how far the data let you trust it. The cleanest picture is Bayes's rule (from Probability & Statistics, 34/05): the probability of the parameters given the data is proportional to the likelihood of the data given the parameters, times what you believed before. With Gaussian noise the likelihood is exp(−χ²/2), where χ² is the sum of squared residuals measured in units of the noise.

The posterior for several parameters rarely has a formula. So astronomers draw samples from it with Markov chain Monte Carlo (Metropolis and colleagues 1953; Hastings 1970). A walker proposes a small random step in parameter space. If the new point fits better, it moves. If it fits worse, it moves anyway with a probability equal to the ratio of the likelihoods. Run it long enough and the places it visits are samples of the posterior, and the cloud they make is the honest answer. Software such as emcee (Foreman-Mackey et al. 2013) does this for most published exoplanet transits.

Without JavaScript: with 200 points and noise of 0.3%, about 40 points fall in the transit, so the depth is known to about 0.05% (σ/√40), and the radius ratio to about 2.5%. Noise of 1% widens that to 0.16%, so the dip is only about six times its own error.

Notice what the cloud shows that a single number cannot: the tilt. Depth and duration trade off against each other, so the cloud leans. Different fields argue about the details. Bayesians treat the posterior as the result; many frequentist astronomers prefer χ² minimisation with bootstrap or profile-likelihood errors. For well-measured data with Gaussian noise the answers agree, and they part only when data are poor or the priors matter.

Selection effects: the catalogue is not the sky

Every survey sees only what is bright enough to detect. Suppose stars scatter in true brightness. At a given distance, the dim ones fall below the survey's limit, and what remains has an average brightness above the real average. The farther away you look, the more extreme the cut: far away you catalogue only the rare luminous ones. This is Malmquist bias (named for Gunnar Malmquist's work in the 1920s), and it affects every distance measure that relies on average brightness. For a population with mean absolute magnitude μ and spread σ, a flux limit sets a cut M_lim = m_lim − 5 log₁₀(d/10 pc), and the mean of the stars that pass is shifted brighter:

The exoplanet census shows it on a grand scale. The NASA Exoplanet Archive lists more than 6,000 confirmed planets as of late 2026. The first found around a Sun-like star, 51 Pegasi b in 1995 (Mayor and Queloz, Nobel Prize 2019), was a "hot Jupiter", a giant planet that orbits in four days. That was not because hot Jupiters are common. They are easy: a planet at 0.05 AU crosses its star as seen from Earth with a probability of about 9% (the star's radius over the orbit's radius), whereas for Earth at 1 AU it is about 0.5%, and its dip is about a hundred times shallower too. The planets we count are the ones our methods find. Astronomers therefore measure completeness (what fraction of real planets of this kind would we have found?) and correct for it, rather than report raw counts.

Without JavaScript: with a limit of 17, spread 1.5 mag, at 2,000 pc, about 56% of the stars are detected and the detected ones average 0.9 mag brighter than the truth.

Interferometry: telescopes that measure Fourier components

The sharpness of a telescope is set by its diameter in wavelengths. To see something 40 millionths of an arcsecond wide, like the ring around the black hole in M87, at a wavelength of 1.3 mm you would need a dish the size of the Earth. You cannot build one, so you build a pair of antennas, then more. Each pair measures one number: the visibility, which is one Fourier component of the sky picture. Its position (u, v) in the Fourier plane is the baseline between the two antennas, in wavelengths:

This is the van Cittert–Zernike theorem (1934, 1938). Eight stations give 28 pairs. As the Earth turns, each baseline sweeps an arc, adding more points, but most of the Fourier plane stays empty. The "image" you get by transforming the points back is the true sky convolved with the pattern of gaps: the dirty image. The art, since Jan Högbom's CLEAN algorithm (1974) and its modern successors, is deciding what the gaps hide. The Event Horizon Telescope's team used several independent imaging methods precisely because different assumptions fill in the empty plane differently. The 2017 campaign used eight telescopes at six sites, at 1.3 mm. Earth's diameter gives θ ≈ 1.3 mm / 12,742 km = 21 microarcseconds, close to the ring size of about 42 µas.

Without JavaScript: with two antennas the dirty image is a set of stripes; with all eight and 12 hours of rotation the bright ring emerges, with ripples around it from the gaps.

Scale: when the catalogue has billions of rows

A routine job: match one catalogue to another, to ask whether this source in an X-ray survey is that star in an optical one. The naive way compares every row to every row. For N objects that is N(N−1)/2 comparisons. At 10,000 objects, 50 million: instant. At 1.8 billion (Gaia) it is about 1.6 × 10¹⁸: at 100 million comparisons a second, a single computer would need roughly 500 years. The fix is an index. Sort the sky into cells (HEALPix, equal-area pixels on the sphere, Górski et al. 2005; or a k-d tree, Bentley 1975) and compare only an object's own cell and its neighbours. The work is then proportional to N, and the same match takes seconds to minutes.

The same thinking runs the alert streams. Rubin's pipeline compares each new image to a template within minutes, and sends the changes to "brokers" that sort them, because no human could look at ten million a night. A graduate is expected to ask of every task: how does this grow with N, and what would I have to change if N were a thousand times bigger? The community uses shared tools for this: FITS files (designed around 1979–81 and still the standard image format), the Astropy library for Python, and public archives such as ESA's Gaia archive, NASA's MAST and the NASA Exoplanet Archive. That openness is why a student in Nairobi or Santiago can analyse the same data as a professor at a large observatory.

Without JavaScript: for N = 1.8 billion, all-pairs matching needs about 1.6 × 10¹⁸ tests (centuries on one machine), while a grid index needs a few billion (seconds to minutes).

What a graduate just sees

• A pixel is a count with a square-root scatter, and doubling the SNR costs four times the exposure.
• Calibration comes first: bias, dark and flat, or the sky is full of ghosts.
• Uneven sampling makes aliases; observe at other times (or other sites) and they go away.
• A fit is a cloud, not a point: the error bar comes from the width of the posterior, and parameters can tilt together.
• Every catalogue is a flux-limited sample: ask what is missing before you ask what is there.
• An interferometer measures Fourier components, and the picture depends on what you assume about the gaps.
• When N grows a thousand times, an N² method is a different project; ask how the work scales.

Where it shows up next

• Planetary Science & Exoplanets (36/40): transit and radial-velocity searches use exactly these periodograms, fits and completeness corrections.
• Galaxies & the Milky Way (36/42): galaxy surveys are flux-limited; Malmquist-type corrections run through every distance and luminosity function.
• Compact Objects & High-Energy Astrophysics (36/43): pulsar searches and the Event Horizon Telescope are signal processing before they are physics.
• Machine Learning (34/32): classifying transients and finding lenses in survey images uses neural networks, which have to cope with imbalanced classes and selection effects.
• In the world: the same skills lead to jobs in data science, medical imaging, remote sensing and software engineering.

If you want to go deeper

• 📘 Textbook. Željko Ivezić, Andrew Connolly, Jacob VanderPlas and Alexander Gray, Statistics, Data Mining, and Machine Learning in Astronomy (Princeton, updated edition 2019). The standard text, with Python code for every figure.
• 📘 Textbook. Steve Howell, Handbook of CCD Astronomy (Cambridge, 2nd edition 2006). Noise, calibration and photometry; and Ian Wall and Brian Jenkins, Practical Statistics for Astronomers (Cambridge, 3rd edition 2012).
• 📘 Gentler. Jacob VanderPlas, Python Data Science Handbook (O'Reilly, 2016), with a free online version, for the tools themselves.
• 🌐 Free text. David Hogg, Jo Bovy and Dustin Lang, "Data analysis recipes: Fitting a model to data" (2010), a free tutorial on doing fits honestly; and the tutorials at learn.astropy.org.