At the resolution of most astronomical observations, stars appear as point-like sources, making them useful for studying how the telescope and, when applicable, the atmosphere shape their images. The resulting intensity distribution is described by the point spread function (PSF).

One important application requiring an accurate estimate of the PSF is weak gravitational lensing, a technique used extensively to study the distribution of dark matter and constrain the properties of dark energy. As light from distant galaxies travels toward us, its path is slightly deflected by the gravitational fields of intervening matter. The resulting distortions in galaxy shapes can be used to map dark matter halos and measure cosmic shear, the coherent deformation of galaxy images across large regions of the sky.

These gravitational distortions are extremely small, typically much weaker than the effects of atmospheric turbulence, optical aberrations, and other instrumental imperfections on the observed galaxy images. Recovering the gravitational lensing signal requires an accurate model of the PSF, allowing its effects to be accounted for when estimating galaxy shapes.

The PSF generally varies across the telescope’s field of view, but stars provide measurements at only a limited number of positions. At galaxy positions, the PSF must be estimated from these stellar measurements using a model of its spatial variation.

One approach is to represent the stellar PSFs using a small number of components, interpolate their coefficients across the field, and reconstruct the complete PSF at the required position.

From Stellar Measurements to Galaxy Positions

The brightness distribution of a galaxy is altered by the telescope’s optics and, for ground-based observations, by the atmosphere. Assuming the PSF is nearly constant across the galaxy image, the resulting image can be described approximately by a convolution

where G is the galaxy’s apparent brightness distribution before instrumental and atmospheric blurring, P_r is the PSF at position r in the focal plane, and N represents observational noise. The symbol * denotes convolution, with x representing the coordinates within the galaxy image.

Optical aberrations and atmospheric effects can change the size, ellipticity, and orientation of the PSF across the field of view, as illustrated below.

At the positions of stars, the PSF can be estimated directly from their observed images, but at galaxy positions, it must be reconstructed from measurements made at nearby stellar positions.

Simulated PSF Field

Comparing different reconstruction methods requires stellar PSF measurements at known positions and reference PSFs against which the reconstructions can be evaluated. Simulated observations provide both, with stellar images generated across the focal plane and exact PSFs available at positions where no star is observed.

The simulations presented here are inspired by the GREAT10 Star Challenge, a public competition organized in 2010 to assess methods for reconstructing spatially varying PSFs from sparsely sampled stellar images. The challenge considered different observing conditions, including optical effects and atmospheric turbulence, and evaluated the accuracy of the reconstructed PSF shapes and sizes. Like GREAT10, the simulations use simplified representations of astronomical PSFs, with controlled spatial variations designed to compare reconstruction methods under different observing conditions.

Three simulated PSF fields are considered, representing different observing conditions and levels of spatial complexity. The first describes smooth variations in atmospheric blurring, the second combines optical distortions with spatially correlated atmospheric fluctuations, and the third represents diffraction-based optics without atmospheric effects. Their principal characteristics are summarized below.

Simulating the stellar PSFs

The first scenario, G1, provides a baseline for comparison with the more complex G2 and S1 fields. It represents atmospheric blurring in ground-based observations using an elliptical Moffat profile:

where α controls the profile width, β determines the extent of its wings, and C is a normalization constant. The radial coordinate is modified to introduce ellipticity and orientation. In G1, the PSF width and ellipticity vary smoothly across the focal plane according to quadratic functions of position, providing a controlled mathematical benchmark.

Scenario G2 extends the Moffat-based model by introducing field-dependent optical aberrations (defocus, astigmatism, and coma) alongside spatially correlated atmospheric variations modeled with Gaussian random fields. Scenario S1 takes a different approach, calculating space telescope PSFs from the pupil geometry and wavefront aberrations through Fourier transforms. These aberrations are represented using Zernike polynomials, a set of mathematical functions commonly used to describe departures from an ideal optical wavefront.

The effects of these PSF models on astronomical images can be illustrated by convolving the same simulated galaxy with a representative PSF from each scenario, as shown below.

Across all three scenarios, the spatial sampling geometry remains identical. The field contains 320 stellar PSFs used to construct the PCA representation, 80 additional positions reserved for validation and hyperparameter selection, and 200 independent test positions reserved for final evaluation. Each PSF is represented by a 48 × 48 pixel image, normalized to unit total flux.

The simulated stellar images contain no observational noise, so reconstruction errors can be examined independently of measurement uncertainties, isolating the intrinsic representation and interpolation accuracy of each method.

Temporal variations of the PSF are not considered in these simulations. Each simulated field represents a fixed observing configuration, although atmospheric conditions and instrumental properties may change during real observations.

Reference PSFs and validation

Simulated data provide exact reference PSFs at every position, making it possible to evaluate reconstructions at locations where no stellar measurements are available. The measurements at the training positions are used to construct a PSF representation, whose parameters are then interpolated to these unseen locations. The complete PSFs are reconstructed from the interpolated parameters and compared with the corresponding reference PSFs to assess their accuracy.

The validation positions are used to select model parameters, while the test positions are reserved for the final assessment. All interpolation methods are evaluated using the same test positions, which were excluded from both the training of the models and the selection of their parameters.

Representing PSFs with Principal Component Analysis

To reconstruct the PSF at galaxy positions, we need a representation that describes its two-dimensional shape with a manageable number of parameters whose spatial variation can be modeled.

Each simulated stellar PSF is represented by a 48 × 48 pixel image, corresponding to 2,304 intensity values. Interpolating every pixel independently across the focal plane would require modeling thousands of spatially varying quantities, many of which are strongly correlated. A more compact description can be obtained by identifying the principal patterns of variation among the stellar PSFs.

Principal component analysis (PCA) provides such a representation. Applied to a collection of stellar images, PCA identifies a set of orthogonal components that describe the dominant variations around their mean. These components, often called eigenPSFs, are ordered according to the amount of variance they explain in the training images.

From stellar images to eigenPSFs

The stellar PSFs are first aligned, normalized, and arranged as vectors containing their pixel intensities. Their mean image is subtracted, and PCA is applied to the resulting collection of deviations. The principal components can then be reshaped into two-dimensional images, each describing a particular pattern of variation within the PSF field.

The mean PSF captures the average intensity distribution, while the first components describe the largest departures from that average. For the simulated Moffat field, these variations are associated with changes in width, ellipticity, and orientation, although an individual component need not correspond to a single physical property.

Any stellar PSF in the field can then be approximated by a linear combination of the mean image and a limited number of principal components:

where P_0 is the mean PSF, φ_k denotes the k-th eigenPSF, and a_k(r) is its coefficient at focal-plane position r. The number K specifies how many components are included in the representation.

The eigenPSFs describe patterns in the pixel intensities, rather than complete stellar images in their own right. Their individual pixel values may be positive or negative, and the coefficients determine how the corresponding patterns contribute to the reconstructed PSF.

Choosing the number of PCA components

Including more components generally improves the representation of the training images, but it also increases the number of coefficient fields that must be interpolated. Choosing K involves balancing the accuracy of the PCA representation against the complexity of the subsequent spatial reconstruction.

Truncating the PCA expansion introduces a representation error even when the coefficients are known exactly, distinct from the interpolation error introduced when they must be estimated at positions without observed stars. This representation error can be measured by projecting reference PSFs directly onto the selected PCA components, without spatial interpolation.

The effect of retaining additional components is illustrated below through the cumulative explained variance in the training images and the representation error measured at held-out test positions. In the smooth G1 field, the explained variance approaches 100% with relatively few components, whereas the diffraction-based S1 field requires a larger basis to capture its more varied image structures. Although the explained variance curves rapidly approach 100%, the representation errors continue to decrease as additional components are included.

PCA thus replaces the original collection of pixel intensities with a small set of coefficients for each stellar position. These coefficients are known for the training stars but must be estimated elsewhere in the field, where the corresponding PSFs are required for galaxy-shape measurements.

Related PCA-based methods have also been used to model spatial variations in the PSF across multiple exposures in weak-lensing surveys, as demonstrated by Jarvis and Jain (2004).

Spatial Interpolation of PCA Coefficients

Once the PCA basis has been constructed, each training star is described by a set of coefficients associated with its position in the focal plane. For each principal component, these coefficients form a spatially sampled field that must be interpolated to the positions of galaxies.

The interpolation methods considered here belong to three broad families: deterministic methods based on polynomial surfaces, distances, or radial basis functions; geostatistical methods that incorporate spatial correlations; and machine-learning approaches that learn relationships between position and PCA coefficients. The seven methods used in the reconstruction experiments are summarized below.

Polynomial and distance-based interpolation

Polynomial regression represents each coefficient field as a function of focal-plane position, using combinations of x, y, and their higher-order terms. A quadratic polynomial, for example, includes terms in x², xy, and y². Such models are well suited to smooth optical variations but may fail to capture smaller-scale fluctuations that depart from the assumed polynomial form.

Inverse-distance weighting (IDW) follows a different approach, estimating coefficients from nearby stellar measurements, with weights that decrease as the distance increases. It requires no explicit model of the spatial trend, although its predictions depend strongly on the distribution of neighboring stars and the chosen distance-weighting power.

Radial basis function (RBF) interpolation constructs a smooth surface from functions centered on the sampled positions. In the experiments considered here, thin-plate splines are used to model the coefficient fields, with a smoothing parameter controlling how closely the fitted surface follows the stellar measurements.

Kriging and machine learning

Kriging incorporates spatial correlation through a variogram, which describes how differences between measurements depend on their separation. Ordinary Kriging assumes an unknown constant mean, whereas Universal Kriging includes a spatially varying mean, or drift, represented here by linear or quadratic polynomials. The latter can describe broad optical trends while using the spatial correlation model to account for departures from those trends.

Machine-learning methods provide another possibility. A multilayer perceptron (MLP), a neural network, learns the relationship between focal-plane coordinates and PCA coefficients from the training stars. Its flexibility makes it suitable for nonlinear variations, although the accuracy of its predictions depends on the available training data, the network architecture, and regularization. A Random Forest provides a second machine-learning approach, combining predictions from multiple decision trees trained on the stellar positions and PCA coefficients. Unlike smooth interpolation models, its predictions are based on piecewise-constant regions of the input space, which can limit its ability to reproduce continuously varying PSF fields.

All methods are applied to the same stellar positions and PCA coefficients. Their parameters are selected using the validation dataset, while the independent test positions are reserved for comparing reconstruction accuracy. The three simulated PSF fields provide different spatial conditions under which to examine the relative strengths and limitations of these approaches.

For further reading on spatial interpolation methods and their validation, see the Spatial Interpolation & Geostatistics series, which covers the mathematical foundations of interpolation, geostatistics, machine learning, and spatial model evaluation.

Reconstructing the PSF

Once the PCA coefficients have been interpolated, the PSF can be reconstructed at any position within the modeled field, including the locations of galaxies where no stellar measurement is available. The reconstruction uses the same mean PSF and principal components derived from the training stars, replacing the measured coefficients with their interpolated values at the target position.

For a galaxy at position r_g, the reconstructed PSF is

where the estimated coefficients depend on the interpolation method, while the mean PSF and principal components remain unchanged.

The resulting image provides an estimate of the complete PSF intensity distribution, rather than only its width or ellipticity. In the simulated fields, this reconstruction can be compared directly with the reference PSF available at the same position, making it possible to examine differences in both the central profile and the outer structure.

The figure below illustrates this comparison using a test position from the S1 space-based simulation, whose diffraction rings and spikes provide a useful example for examining reconstruction differences beyond the central PSF profile. The reference PSF, its reconstruction, and the residual image are shown side by side.

Measuring Reconstruction Accuracy

The reconstructed PSFs can be compared with their exact counterparts at the 200 held-out test positions. This comparison must account for both the overall intensity distribution and the properties most relevant to galaxy-shape measurements, since a reconstruction may reproduce the central PSF profile closely while introducing small errors in its ellipticity or size. For this purpose, we use several complementary measures of accuracy. Other approaches to PSF model evaluation, including diagnostics based on spatial correlations in reconstruction residuals, are discussed by Rowe (2010).

Image reconstruction error

To assess how closely the reconstructed PSFs reproduce their reference images, we first examine the pixel-by-pixel differences, illustrated by the residual image in the preceding S1 example. These differences are then summarized across the full test set using the relative reconstruction error

where P_i denotes the reference PSF image at test position i, and the same symbol with a hat denotes its reconstruction. The Frobenius norm sums the squared pixel differences over the image, while the normalization expresses the total reconstruction error relative to the reference PSFs.

This measure includes errors in the PSF core, wings, and diffraction structure, but does not indicate which aspects of the image are most affected. It combines two contributions: the representation error caused by truncating the PCA expansion and the additional error introduced by spatially interpolating its coefficients.

To isolate the representation error, each reference PSF is first decomposed directly into PCA coefficients using the selected principal components. A PSF reconstructed from these coefficients contains no interpolation error, since the coefficients are obtained from the image itself rather than estimated from neighboring stars. The remaining difference between this truncated reconstruction and the reference PSF then measures the error introduced by retaining only a limited number of components.

The two reconstruction paths are illustrated below, using the same reference PSF to separate the effect of PCA truncation from the additional error introduced by interpolating its coefficients.

PSF ellipticity and size

For weak-lensing applications, the shape and size of the PSF are particularly important because inaccuracies in these properties can affect the inferred shapes of galaxies. Following the usual description in terms of second-order (quadrupole) moments, we characterize each PSF by a matrix Q_ij, which measures how its intensity is distributed around the image center. In the present analysis, the moments are calculated with a Gaussian weight to limit the influence of the faint outer regions.

The two ellipticity components are defined as

while the squared PSF size is

The first ellipticity component describes elongation along the image coordinate axes, while the second describes elongation along directions rotated by 45°. Together, they characterize the anisotropy of the PSF, whereas R² measures the overall spatial extent of its weighted intensity distribution.

For each quantity, reconstruction accuracy is evaluated using the root-mean-square error (RMSE) across the test positions. The ellipticity errors are reported separately for e₁ and e₂, and combined into a single measure

Errors in squared PSF size are expressed in squared pixels. These measurements complement the full-image error by identifying differences that are especially relevant when correcting observed galaxy shapes.

Additive and multiplicative biases

The RMSE measures the magnitude of reconstruction errors but does not distinguish random scatter from systematic effects. An interpolation method may, for example, reproduce the mean ellipticity while systematically underestimating its variations across the field.

To examine such effects, we fit a linear relation between reconstructed and reference ellipticity at the test positions:

Here m_j represents a multiplicative bias, c_j an additive offset, and ϵ_j the remaining residual. A negative multiplicative bias indicates that the reconstructed ellipticity varies less strongly than the reference values, while an additive bias shifts the reconstructed values by an approximately constant amount.

For the PSF size, we measure the mean fractional bias in R²:

This quantity indicates whether a method tends to overestimate or underestimate the squared PSF size. It is considered alongside the size RMSE, since positive and negative errors can cancel in the mean even when individual reconstructions remain inaccurate.

These measures provide a basis for comparing the interpolation methods across the three simulated PSF fields. The fitted additive and multiplicative biases describe systematic errors in the reconstructed PSF properties, rather than in the gravitational lensing signal itself. Determining how these errors affect cosmic shear measurements would require using the reconstructed PSFs to correct simulated galaxy images and evaluating the resulting biases in the recovered galaxy shapes.

Evaluation of PCA-Based PSF Modeling

The seven interpolation methods were evaluated across the three simulated PSF fields, giving 21 method–scenario combinations. The PCA representations use five components for G1, eight for G2, and nine for S1, with reconstruction accuracy measured at the same 200 held-out test positions in each scenario.

The comparison examines how the methods recover the complete PSF intensity distribution, its ellipticity components, and its squared size. Although these quantities are related, a method that minimizes the overall image error does not necessarily provide the most accurate estimates of PSF shape.

Reconstruction accuracy across observing scenarios

The first comparison concerns the ellipticity components, whose reconstruction errors are shown below for all seven methods and three observing scenarios.

In the smooth Moffat field G1, polynomial regression achieves the lowest combined ellipticity RMSE, 0.063 × 10⁻³, followed by Universal Kriging and RBF splines at 0.218 × 10⁻³ and 0.231 × 10⁻³, respectively. The relatively simple spatial dependence of the PSF parameters favors methods that capture broad, smoothly varying trends. More flexible models offer little advantage when those trends can already be described accurately by low-degree polynomials.

The relative performance of the methods changes in G2, where optical distortions are combined with spatially correlated atmospheric variations. Polynomial regression performs poorly, with a combined ellipticity RMSE of 15.292 × 10⁻³, while the RBF spline reduces this error to 2.487 × 10⁻³. Ordinary and Universal Kriging provide intermediate results, with combined errors of approximately 3.2 × 10⁻³. The additional fluctuations in G2 cannot be represented adequately by a single low-degree polynomial surface, whereas RBF and Kriging models can accommodate departures from a broad spatial trend.

The space-based scenario S1 presents a different situation. Although its individual PSFs contain diffraction rings and structures associated with the telescope pupil, their spatial variations remain sufficiently regular for polynomial regression and Universal Kriging to perform well. Universal Kriging achieves the lowest combined ellipticity RMSE, 0.359 × 10⁻³, closely followed by polynomial regression at 0.368 × 10⁻³.

The same comparison can be extended to squared PSF size and the complete reconstructed images.

For G1, Universal Kriging gives the smallest relative image reconstruction error, 0.016%, while polynomial regression remains the most accurate method for ellipticity. In S1, polynomial regression achieves the smallest image error, 0.111%, whereas Universal Kriging gives slightly better ellipticity and squared-size estimates. These differences show that the method providing the closest pixel-level reconstruction is not always the one that best preserves the particular PSF properties used in galaxy-shape measurements.

G2 remains the most challenging scenario across the principal accuracy measures. Even its leading RBF reconstruction has a relative image error of 0.657%, compared with 0.019% in G1 and 0.135% in S1 for the same method. The contrast is significant because the diffraction-based S1 profiles are visually more structured than the ground-based G2 profiles. The accuracy of the complete reconstruction depends not only on the structure of an individual PSF image, but also on how effectively it can be represented by the selected PCA components and how accurately their coefficients can be interpolated across the field.

The results do not identify a single interpolation method that is optimal for all three PSF fields and all accuracy measures. Instead, they indicate that the choice of model must account for the spatial behavior of the coefficient fields and the PSF properties whose accuracy is most important.

Systematic biases in PSF properties

The RMSE comparisons describe the magnitude of reconstruction errors, but they do not show whether a method systematically changes the range of ellipticities or sizes present across the field. This distinction is particularly relevant in G2, where several methods produce reconstructed PSFs whose properties vary less strongly than those of the reference images.

The effect is illustrated below by comparing reference and reconstructed ellipticity components and squared sizes for three interpolation methods in G2.

The RBF reconstruction closely follows the reference values, with multiplicative ellipticity biases of -0.77% for e₁ and -1.94% for e₂. Universal Kriging exhibits stronger attenuation, with biases of -3.93% and -3.60%, while polynomial regression reaches -38.69% and -32.59%. For the latter method, the reconstructed values are compressed toward the mean, consistent with its inability to reproduce the correlated spatial fluctuations present in the simulated coefficient fields.

The squared-size panel shows a related effect. The fitted slopes for RBF and Universal Kriging remain close to unity, whereas polynomial regression produces a substantially flatter relation and greater scatter. A fitted slope, however, describes how reconstructed sizes vary with the reference values; it is not equivalent to the mean fractional size bias introduced earlier.

The systematic effects across all three scenarios are summarized in the following comparison of multiplicative ellipticity biases, additive offsets, and mean fractional squared-size biases.

In G1, polynomial regression, RBF splines, Universal Kriging, and the multilayer perceptron produce relatively small ellipticity biases, whereas inverse-distance weighting and Random Forest show more pronounced attenuation. G2 exhibits the strongest negative multiplicative biases, particularly for polynomial regression and Random Forest, reflecting the difficulty of recovering its smaller-scale spatial variations. In S1, the leading methods generally preserve ellipticity variations more accurately, although their fitted biases are not uniformly zero.

The additive offsets and mean fractional size biases are generally smaller than the corresponding variation measured through RMSE, especially for the better-performing methods. This does not imply that their individual size estimates are equally accurate, since positive and negative residuals can compensate when averaged across the field.

Taken together, the benchmark results show that spatial interpolation can introduce systematic changes in PSF properties even when the reconstructed images appear close to their reference counterparts. The choice of interpolation method must therefore consider both the overall reconstruction error and the accuracy of the shape and size measurements derived from the PSF.

Discussion and Limitations

The numerical results illustrate how the accuracy of PSF modeling depends on both the representation of stellar images and the spatial variation of their PCA coefficients. In the smooth G1 field, low-degree polynomial regression reproduces the variations in PSF shape with high accuracy, while RBF splines and Kriging are better suited to the correlated fluctuations introduced in G2. Despite the more complex optical structure of individual PSFs in the diffraction-based S1 field, their spatial variations can be modeled accurately with relatively simple interpolation methods. This shows that the complexity of a PSF image is not, by itself, a reliable indicator of the difficulty of estimating that PSF at unsampled positions.

These differences in reconstruction accuracy must also be considered in relation to the two stages of the modeling pipeline. The reported errors combine the effects of PCA truncation and spatial interpolation. Although these contributions can be distinguished by projecting the reference PSFs directly onto the selected principal components, the numerical comparisons presented here evaluate the complete modeling pipeline rather than assigning separate error budgets to the two stages. Within each scenario, the interpolation methods share the same PCA representation, allowing their relative performance to be compared under consistent conditions. Comparisons between scenarios, however, also reflect differences in the underlying PSF models and in the number of components used to describe them.

Beyond the distinction between PCA representation and interpolation errors, the evaluation is limited by the assumptions of the simulations. In particular, the stellar PSFs contain no observational noise, whereas real measurements are affected by photon noise, detector noise, and uncertainties in the determination of stellar profiles. Such effects can influence both the PCA representation and the estimated coefficient fields, potentially changing the relative performance of the interpolation methods. The simulations also describe fixed observing configurations and do not account for temporal or wavelength-dependent PSF variations, which can be important in astronomical surveys.

In astronomical surveys, dedicated software packages such as PSFEx (Bertin 2011) and Piff (Jarvis et al. 2021) provide more comprehensive frameworks for fitting stellar PSFs and modeling their spatial variations, whereas the present experiments focus on PCA-based reconstruction and interpolation under controlled conditions.

The spatial sampling strategy imposes further limits on the interpretation of the results. The evaluation uses the same stellar sampling geometry across the three scenarios, with validation and test positions drawn from the modeled fields. This provides a consistent basis for comparing the methods, but does not establish how their accuracy would change with different stellar densities, irregular spatial coverage, or predictions outside the well-sampled regions. The confidence intervals reported for systematic biases also assume that reconstruction residuals at different test positions are statistically independent. However, spatially correlated PSF variations may produce correlated residuals at nearby positions, reducing the effective number of independent measurements and potentially leading to underestimated uncertainties.

Finally, the ellipticity and size biases measured here characterize errors in the reconstructed PSFs, not biases in the gravitational lensing signal itself. Their effect on cosmic shear measurements depends on how the PSF model is used to correct galaxy images, as well as on the galaxy population and shape-measurement procedure. Establishing that connection would require applying the reconstructed PSFs to simulated galaxy observations and measuring the resulting errors in recovered galaxy shapes.

Conclusion

Estimating the point spread function at galaxy positions is an essential step in accounting for instrumental and atmospheric effects in astronomical images. The approach examined here uses principal component analysis to represent stellar PSFs through a limited number of components, whose coefficients are interpolated across the field of view to reconstruct the PSF at positions where no direct stellar measurement is available.

The comparison of seven interpolation methods across three simulated PSF fields shows that reconstruction accuracy depends on the spatial variation of these coefficients and on the ability of the PCA representation to describe the underlying PSF structure. No single interpolation method performs best under all conditions, and the method providing the lowest full-image reconstruction error does not necessarily yield the most accurate estimates of PSF ellipticity or size.

For applications such as weak gravitational lensing, these results emphasize the importance of evaluating PSF models not only through their agreement with reference images, but also through the accuracy of the shape and size properties used in galaxy measurements.

You can explore the implementation, evaluation pipeline, and figure generation scripts in the GitHub companion repository.

References & Further Reading

1. Matheron, G. (1963). Principles of geostatistics. Economic Geology, 58(8), 1246–1266. DOI: 10.2113/gsecongeo.58.8.1246.
2. Chorley, R. J., & Haggett, P. (1965). Trend-surface mapping in geographical research. Transactions of the Institute of British Geographers, (37), 47–67. DOI: 10.2307/621689.
3. Shepard, D. (1968). A two-dimensional interpolation function for irregularly-spaced data. Proceedings of the 1968 23rd ACM National Conference, 517–524. DOI: 10.1145/800186.810616.
4. Moffat, A. F. J. (1969). A Theoretical Investigation of Focal Stellar Images in the Photographic Emulsion and Application to Photographic Photometry. Astronomy and Astrophysics, 3, 455–461. ADS.
5. Harder, R. L., & Desmarais, R. N. (1972). Interpolation using surface splines. Journal of Aircraft, 9(2), 189–191. DOI: 10.2514/3.44330.
6. Cressie, N. (1993). Statistics for Spatial Data. John Wiley & Sons, New York.
7. Kaiser, N., Squires, G., & Broadhurst, T. (1995). A Method for Weak Lensing Observations. The Astrophysical Journal, 449, 460–480. DOI: 10.1086/176071 · ADS · arXiv:astro-ph/9411005.
8. Breiman, L. (2001). Random Forests. Machine Learning, 45, 5–32. DOI: 10.1023/A:1010933404324 · Springer · PDF (UC Berkeley).
9. Jarvis, M., & Jain, B. (2004). Principal Component Analysis of PSF Variation in Weak Lensing Surveys. ADS · arXiv:astro-ph/0412234.
10. Goodman, J. W. (2005). Introduction to Fourier Optics (3rd ed.). Roberts & Company Publishers. ISBN: 978-0974707723. Google Books · Open Library. (A standard reference for diffraction theory, pupil functions, and image formation).
11. Bishop, C. M. (2006). Pattern Recognition and Machine Learning. Springer. ISBN: 978-0387310732.
12. Rowe, B. (2010). Improving PSF modelling for weak gravitational lensing using new methods in model selection. Monthly Notices of the Royal Astronomical Society, 404(1), 350–366. DOI: 10.1111/j.1365–2966.2010.16277.x · Oxford Academic · ADS · arXiv:0904.3056.
13. Bertin, E. (2011). Automated Morphometry with SExtractor and PSFEx. Astronomical Data Analysis Software and Systems XX, ASP Conference Series, 442, 435. ASP Conference Series · PDF. (PSFEx: a package for modeling spatial variations in point spread functions).
14. Kitching, T., Balan, S. T., Bernstein, G., et al. (2011). Gravitational Lensing Accuracy Testing 2010 (GREAT10) Challenge Handbook. Annals of Applied Statistics, Vol. 5, Issue. 3, pp. 2231–2263. 10.1214/11-AOAS484 · ADS · arXiv:1009.0779.
15. Kitching, T. D., Rowe, B., Gill, M., et al. (2013). Image Analysis for Cosmology: Results from the GREAT10 Star Challenge. The Astrophysical Journal Supplement Series, 205(2), 12. 10.1088/0067–0049/205/2/12 · ADS · arXiv:1210.1979.
16. Mandelbaum, R. (2018). Weak Gravitational Lensing for Precision Cosmology. Annual Review of Astronomy and Astrophysics, 56, 393–433. DOI: 10.1146/annurev-astro-081817–051928 · ADS · arXiv:1710.03235.
17. Jarvis, M., Bernstein, G. M., Amon, A., et al. (2021). Dark Energy Survey Year 3 Results: Point Spread Function Modelling. Monthly Notices of the Royal Astronomical Society, 501(1), 1282–1299. DOI: 10.1093/mnras/staa3679 · Oxford Academic. (Piff: a framework supporting different PSF model representations and spatial interpolation methods).

Interpolating Stars for Galaxy Shape Measurement was originally published in Stellar Priors on Medium, where people are continuing the conversation by highlighting and responding to this story.

What does it mean to understand a differential equation?

Is finding an explicit solution always necessary? And when no simple formula exists, what can the geometry of the possible motions still tell us?

These questions lead to a different way of thinking about dynamics. Instead of following one solution through time, we can study how entire families of trajectories are arranged, how they behave near special states, and how their long-term evolution changes when the system is perturbed.

The nonlinear pendulum provides a simple place to begin. Its equation is easy to write down, but it already shows why solving an equation and understanding a dynamical system are not always the same thing.

A deceptively simple equation

A pendulum is one of the simplest systems in classical mechanics.

A mass hangs from a pivot and moves under the action of gravity. If θ denotes its angular displacement from the downward vertical and ℓ the pendulum length, its motion is described by

For small θ,

so the equation becomes

We recognize the harmonic oscillator, with its well-known sinusoidal solution

In this regime, the motion can be written explicitly and is easy to interpret

However, once the small-angle approximation is removed, the sin θ term makes the system nonlinear and the equation is no longer solved by elementary trigonometric functions. Its exact solution instead involves elliptic functions, a class of special functions, so the explicit time evolution is less straightforward to interpret.

The nonlinear pendulum can behave in several qualitatively different ways. It may oscillate back and forth around its downward equilibrium, or, with sufficient initial velocity, rotate continuously around the pivot. The downward and upward equilibrium positions have different stability properties, and a special trajectory separates oscillatory motion from continuous rotation.

By the late nineteenth century, differential equations had become fundamental to mechanics and mathematical physics, but many nonlinear equations resisted simple analytical solutions, and their possible behaviors were not always apparent from the equations themselves. The difficulty was particularly important in celestial mechanics, where questions about long-term motion and stability remained even when explicit solutions were unavailable.

Henri Poincaré began developing a qualitative theory of differential equations in the 1880s, asking what could be learned about solutions even when they could not be written explicitly. Rather than following only individual solutions, he studied the geometry formed by entire families of trajectories.

If we take the example of the pendulum, the idea is to describe its instantaneous state not by its angle θ alone, but by a point in a plane whose coordinates are its angle and angular velocity,

As the pendulum evolves, this point traces a curve. Different initial conditions generate different curves, and taken together these trajectories describe the geometry of the motion.

So, instead of asking only for the solution θ(t), we can ask how all possible trajectories are arranged in state space. This change of viewpoint is one of the central ideas of dynamical systems theory.

It allows us to identify equilibrium points, determine their stability, distinguish oscillation from rotation, and find the trajectories that separate different regimes of motion, often without requiring an explicit formula for the full time evolution.

The pendulum is simple enough for these ideas to be seen directly, but the same viewpoint extends much further, from nonlinear mechanics and celestial dynamics to stability theory and chaos.

From an Equation of Motion to a Dynamical System

The pendulum equation

is a second-order differential equation: determining its future motion requires both the initial angle and the initial angular velocity.

If we introduce

then the second-order equation can be written as a pair of first-order equations,

If the state of the pendulum is specified by the two-component vector

then its evolution can be written compactly as

The change of viewpoint, from the physical pendulum to a trajectory in state space, is summarized below.

This is a first-order autonomous dynamical system. It is first-order because only first derivatives appear, and autonomous because the right-hand side contains no explicit dependence on time. The evolution of the system at any instant is determined entirely by its current state.

The same form describes a much wider class of systems:

where the components of x specify the variables needed to describe the state of the system. A second-order equation for one variable therefore becomes a first-order system in two variables. More generally, higher-order equations and coupled equations can be rewritten in a similar form by enlarging the state vector.

At every point (θ, ω), the vector f(θ, ω) then specifies the instantaneous direction in which the state will move. Taken over the whole plane, these vectors define a vector field. Starting from an initial state (θ_0, ω_0), the evolving pendulum follows one curve through this field: its trajectory.

The collection of possible states forms the state space of the system. For mechanical systems, this is commonly called phase space. When the relevant phase space is two-dimensional, as for the pendulum considered here, it is called the phase plane. Each point represents a complete instantaneous state of the pendulum, while each trajectory represents one possible evolution in time.

💡 Note on phase space coordinates

In Hamiltonian mechanics, phase space uses (θ, p_θ) with p_θ = mℓ²ω. For fixed m and ℓ, this is just a constant rescaling of (θ, ω), so the qualitative phase-space geometry is unchanged.

The phase-plane representation is different from a conventional plot of θ(t). A time-series plot shows how one variable evolves for one chosen initial condition. The phase plane instead places all possible states in the same geometric space, allowing different kinds of motion to be compared directly.

For the pendulum, that geometry is organized into distinct dynamical regimes. Some trajectories form closed curves, others correspond to continuous rotations, and certain points remain fixed altogether.

To see why these different regimes appear, we can now examine the full set of trajectories in the phase plane.

Phase Plane and the Geometry of Motion

The phase plane makes it possible to examine all possible motions of the pendulum at once.

For the undamped pendulum, there is an additional quantity that organizes these trajectories: the total mechanical energy. Taking the lowest position of the pendulum, θ = 0, as the zero of gravitational potential energy, the total mechanical energy is

The first term is the kinetic energy and the second is the gravitational potential energy. In the absence of friction, this energy is conserved. Differentiating with respect to time gives

Using

the two terms cancel,

The pendulum thus remains on the same energy level throughout its motion. For a given value of E, its trajectory in the (θ, ω) plane must satisfy

Each trajectory is then a curve of constant energy. Plotting these trajectories for different initial conditions gives the phase portrait of the pendulum.

The phase portrait separates the motion into three qualitatively different regimes.

For sufficiently low energy, the pendulum cannot reach the vertically upward position. It moves back and forth around θ = 0, reversing direction whenever ω = 0. These bounded oscillations, usually called librations, appear in the phase plane as closed curves surrounding the downward equilibrium.

The boundary between libration and rotation occurs at a critical energy. At the vertically upward position,

and the energy is

At exactly this energy value, the trajectory forms the separatrix, the boundary dividing libration from rotation.

At higher energy, the pendulum has enough kinetic energy to pass over the top and continue moving in the same direction. This regime is called rotation; when the evolving variable is an angle, it is also commonly described as circulation. In the phase plane, rotational trajectories extend across successive values of θ rather than closing around the downward equilibrium.

During libration, ω changes sign at each turning point, whereas during rotation it retains the same sign. Since p_θ = mℓ²ω, the conjugate momentum behaves in the same way.

The qualitative behavior of the pendulum can thus be read directly from the geometry of its phase portrait, without solving explicitly for θ(t).

There is one further detail in this picture. The angle θ is periodic: values differing by 2π represent the same physical orientation. The periodicity of θ means that the pendulum's phase space can also be represented as a cylinder, obtained by identifying

Under this identification, librations remain closed loops around the stable equilibrium, while rotational trajectories wind continuously around the cylinder.

The phase portrait already distinguishes the main types of pendulum motion and shows how much information can be extracted from the geometry of trajectories alone. But mechanical motion and energy conservation are special features of this example. Before turning to fixed points and stability, it is useful to separate what belongs specifically to the pendulum from what applies more generally to dynamical systems.

Beyond the Pendulum: General Dynamical Systems

The pendulum provides a convenient example because its physical interpretation is immediate and its phase portrait already displays equilibria, closed trajectories, and a separatrix. But the same mathematical framework applies much more broadly.

A continuous-time autonomous dynamical system can be written in the general form

where

collects the variables required to specify the state of the system. The function f determines how that state changes from one instant to the next. In one dimension,

the state is described by a single variable. In two dimensions,

the state moves through a plane. More generally, an n-dimensional system evolves through an n-dimensional state space.

The state variables are simply the quantities needed to specify the system. They may represent interacting populations, chemical concentrations, electrical currents and voltages, or the coordinates and momenta of many bodies. They need not describe mechanical motion at all.

The pendulum has another property that is not shared by dynamical systems in general: in the absence of friction, it conserves mechanical energy. Energy conservation is specific to this example, not a defining property of dynamical systems. Other systems may dissipate energy, exchange it with their surroundings, or possess no meaningful mechanical energy at all.

Time can also enter the equations explicitly,

in which case the system is called non-autonomous. Here we will continue to focus mainly on autonomous systems, for which the current state alone determines the local direction of evolution.

Dynamical systems can also evolve in discrete steps rather than continuously, but the present discussion will remain with systems described by differential equations.

Despite these differences, the central geometric question remains the same: how does the system evolve through its state space?

We can begin by asking the opposite question: are there states that do not evolve at all? These are the fixed points of the system, and their response to small perturbations leads directly to the concept of stability.

Fixed Points and Local Stability

For an autonomous dynamical system

a fixed point is a state x_* for which

If the system starts exactly at this point, it remains there forever, but this alone tells us nothing about what happens to nearby states. A fixed point is called stable if sufficiently small perturbations remain small. If nearby trajectories not only remain close but converge toward the fixed point, it is asymptotically stable. On the other hand, if arbitrarily small perturbations can grow and carry the system away, the fixed point is unstable.

This distinction was placed on a systematic mathematical basis by Aleksandr Lyapunov in his work on the stability of motion.

For nonlinear systems, the dynamics close to a fixed point can often be understood by examining a linear approximation. The idea is to start close to the equilibrium and ask how a small displacement evolves: does it grow, decay, or oscillate?

Suppose the state is displaced slightly from x_*,

where δx is small. Expanding the vector field around the fixed point gives

where

is the Jacobian matrix of the system. The Jacobian describes how the vector field changes locally when the state is displaced from the fixed point.

Since f(x_*) = 0, the local dynamics reduce to

So we can use this linearized system to study how small perturbations evolve. The relevant information is contained in the eigenvalues and eigenvectors of the Jacobian. The eigenvalues are obtained from the characteristic equation

If an eigenvalue has a negative real part, the associated perturbation decays; if it has a positive real part, the perturbation grows. A nonzero imaginary part produces oscillations. The associated eigenvectors determine the corresponding directions in state space.

For a two-dimensional system, these different eigenvalue combinations correspond to recognizable local geometries, summarized below.

The last case requires particular care. If the eigenvalues have zero real part, linearization by itself may not determine the stability of the full nonlinear system. Additional information about the nonlinear terms, conserved quantities, or other properties of the system may be required.

We can now apply this classification to the two equilibrium positions of the pendulum.
Fixed points require ω = 0 and sin θ = 0. Since the equilibrium conditions repeat whenever θ is shifted by 2π, there are two physically distinct configurations,

corresponding respectively to the downward and inverted positions.

The Jacobian of the pendulum system is

At the downward equilibrium, θ = 0, substituting cos 0 = 1 into the Jacobian gives the characteristic equation and resulting eigenvalues

The eigenvalues are imaginary, meaning the perturbations are oscillatory. In the phase portrait, nearby trajectories form closed curves around the equilibrium, which is called a center.

As noted above, purely imaginary eigenvalues do not by themselves establish nonlinear stability. For the pendulum, however, the conserved energy provides the additional information: trajectories sufficiently close to the downward equilibrium remain on nearby closed energy curves. The equilibrium is therefore stable, but not asymptotically stable.

At the inverted equilibrium, θ = π, the result changes to

One eigenvalue is negative and the other positive. Perturbations contract along one eigendirection and grow along the other, so the inverted equilibrium is an unstable saddle point.

This also explains a feature already visible in the phase portrait. The stable and unstable directions of the saddle extend into the surrounding phase plane, forming the trajectories associated with the separatrix between libration and continuous rotation.

For the ideal pendulum, however, the local analysis also shows a limitation: the stable downward equilibrium does not attract nearby trajectories. With no mechanism for losing energy, a displaced pendulum continues to oscillate indefinitely. Real systems behave differently when dissipation is present.

Dissipation and Attractors

The ideal pendulum conserves energy, so each trajectory remains on a constant-energy curve of the phase portrait. A real pendulum, however, loses energy through effects such as friction and air resistance.

A simple model adds a damping term to the equation of motion,

where γ > 0 measures the strength of the damping.

The mechanical energy is still

but it is no longer conserved. Differentiating with respect to time and using the damped equations of motion gives

meaning energy decreases whenever the pendulum is moving. In the phase plane, trajectories no longer remain on curves of constant energy. Instead, they cross successive energy levels as the oscillations lose amplitude.

The effect can also be understood from the local stability analysis developed in the previous section. At the downward equilibrium (θ, ω) = (0, 0), the Jacobian becomes

Its characteristic equation is

with eigenvalues

For weak damping,

the eigenvalues are complex conjugates with a negative real part. The center of the undamped pendulum is therefore replaced by a stable spiral: nearby trajectories spiral inward as the oscillations decay.

For sufficiently strong damping, the eigenvalues become real and negative, and the equilibrium instead behaves as a stable node. In both cases, the downward equilibrium is now asymptotically stable.

Without damping, nearby trajectories remain on closed energy curves around the center. With damping, they cross successive energy levels and approach the downward equilibrium as time increases, replacing the closed curves of the undamped system with inward trajectories.

The equilibrium therefore acts as a local attractor: states starting sufficiently close to it converge toward the same long-term state. Dissipation changes not only the energy of the system, but also the geometry of its flow, turning a center into an attractor.

Poincaré and the Qualitative View of Dynamics

Fixed points, linearization, and attractors describe important features of a dynamical system, particularly near an equilibrium. But they do not by themselves describe how trajectories are arranged across the full state space.

This broader geometric question was central to the work of Henri Poincaré in the late nineteenth century. Instead of asking only whether a differential equation could be solved explicitly, Poincaré studied the geometry and long-term evolution of its trajectories: where they could go, whether they returned, and how different kinds of motion were related.

One important possibility is periodic motion. The closed libration trajectories of the undamped pendulum provide a simple example.

A trajectory is periodic if there exists a period T > 0 such that

In state space, this translates into a closed orbit.

For systems with several degrees of freedom, however, the geometry of a full trajectory can become difficult to visualize. Poincaré introduced a way of reducing this complexity by recording only the points where a trajectory crosses a chosen lower-dimensional surface.

Let this surface be denoted by Σ. If a trajectory intersects it successively at

the evolution can be represented by a Poincaré return map

The continuous flow is thereby represented by a discrete sequence of section crossings.

This reduction preserves important information about the dynamics. If a periodic orbit returns to the same point on the section after each cycle, that point satisfies

The periodic orbit of the continuous system has become a fixed point of the return map. Questions about the stability of the orbit can then be studied through the evolution of nearby points under repeated applications of P.

The same idea extends beyond periodic motion. A complicated continuous flow may produce a much simpler geometric pattern when viewed only through its successive intersections with a section. The return map can therefore make recurring patterns and stability properties easier to identify than in the original differential equations.

Poincaré’s aim was not necessarily to obtain an explicit expression for x(t), but to determine the geometry of the motion: equilibria, periodic orbits, invariant trajectories, and the relationships between them.

Poincaré later applied this viewpoint to celestial mechanics, where the motion of several gravitating bodies provided precisely the kind of nonlinear system for which explicit solutions were generally unavailable. Return maps and related geometric ideas became part of the foundation of modern dynamical systems theory and, later, the mathematical study of chaotic motion.

From Equations to Geometry

We began with a question: to understand a dynamical system, is it always necessary to solve its equations of motion explicitly?

The nonlinear pendulum shows that the answer is no. Without writing its full solution in terms of elliptic functions, the phase portrait already distinguishes libration from rotation, identifies the separatrix between them, and shows that the downward equilibrium is stable while the inverted equilibrium is unstable.

Local analysis adds another level of understanding. Linearization around a fixed point, together with the eigenvalues and eigenvectors of the Jacobian, determines whether nearby trajectories remain close, move away, or converge toward the equilibrium. Adding dissipation changes this geometry again, replacing the closed trajectories of the ideal pendulum by trajectories that approach an attractor.

Poincaré’s qualitative approach extends the same idea beyond the pendulum. Periodic orbits, return maps, and trajectories across state space can be studied without requiring an explicit expression for every possible motion.

Explicit solutions remain valuable when they are available. But dynamical systems theory shows that they are not the only way to understand a nonlinear system.

The Python code used to generate all figures and diagrams in this article is open-source and available in the GitHub companion repository.

References & Further Reading

1. Poincaré, H. (1881–1886). Mémoire sur les courbes définies par une équation différentielle / Sur les courbes définies par les équations différentielles. Journal de Mathématiques Pures et Appliquées. Four parts: Part I (1881), Part II (1882), Part III (1885), Part IV (1886).
2. Poincaré, H. (1890). Sur le problème des trois corps et les équations de la dynamique. Acta Mathematica, 13, 1–270. Henri Poincaré Papers.
3. Poincaré, H. (1892–1899). Les Méthodes Nouvelles de la Mécanique Céleste. Gauthier-Villars, Paris.
4. Lyapunov, A. M. (1892; French trans. 1907). Problème général de la stabilité du mouvement. Annales de la Faculté des Sciences de Toulouse, 9, 203–474. Numdam.
5. Birkhoff, G. D. (1927). Dynamical Systems. American Mathematical Society Colloquium Publications, Vol. 9.
6. Guckenheimer, J., & Holmes, P. (1983). Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields. Springer-Verlag.
7. Arnold, V. I. (1992). Ordinary Differential Equations. Springer-Verlag, Berlin Heidelberg. Springer.
8. Perko, L. (2001). Differential Equations and Dynamical Systems. 3rd ed., Springer, Texts in Applied Mathematics, Vol. 7. 10.1007/978–1–4613–0003–8 · Springer.
9. Hirsch, M. W., Smale, S., & Devaney, R. L. (2013). Differential Equations, Dynamical Systems, and an Introduction to Chaos. 3rd ed., Academic Press. Elsevier.
10. Strogatz, S. H. (2015). Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry, and Engineering. 2nd ed., Westview Press.

Dynamical Systems: The Geometry Behind Differential Equations was originally published in Stellar Priors on Medium, where people are continuing the conversation by highlighting and responding to this story.

On September 8, 2026, OpenAI announced a proposed solution to the Navier–Stokes existence and smoothness problem, one of the seven Millennium Prize Problems. According to OpenAI, the successful search involved roughly 10,000 concurrent agents and reached its result after about 88 hours. A formal version was subsequently encoded and verified in Lean.

Three days later, the Clay Mathematics Institute issued an unusually strong statement: the Navier–Stokes problem had “apparently been settled”.

The result, if it survives the full process of mathematical scrutiny, addresses a question that has resisted analysis since the foundational work of Jean Leray in the 1930s.

The Navier–Stokes equations express Newton’s law of motion for a fluid, describing how velocity evolves under inertia, pressure, viscosity and external forces. Their form is compact, yet their three-dimensional behavior has proved remarkably difficult to control mathematically.

For three-dimensional flow, one fundamental question remained unresolved:

If the fluid starts smooth, must it remain smooth forever?

Viscosity dissipates kinetic energy and smooths variations in velocity. But three-dimensional fluid motion can also stretch and intensify vortices, generating increasingly fine spatial structure. The central difficulty is whether this nonlinear amplification can become strong enough to overcome viscous smoothing and produce a finite-time singularity.

Understanding that question requires looking at the mathematical structure of the equations: advection, vorticity and vortex stretching, energy dissipation, and the way Navier–Stokes behaves under changes of spatial scale.

The story therefore begins not with artificial intelligence, but with Newton’s second law and the dynamics of viscous fluid motion.

From Newton’s Law to Navier–Stokes

To understand where the Navier–Stokes equations come from, consider a small parcel of fluid moving through a flow. There are two natural ways to describe its motion.

In the Eulerian description, we stand at a fixed position x and describe the velocity field

This tells us the velocity of whatever fluid happens to pass through that point at time t.

In the Lagrangian description, we instead follow one particular parcel. If its position is X(t), then its trajectory satisfies

Newton’s second law applies to this moving parcel. What we need is not simply the change of velocity at a fixed point, but the change in velocity experienced by the parcel as it moves through the flow.

Along its trajectory, the parcel has velocity

Differentiating with respect to time and applying the chain rule gives

Since

the acceleration becomes

This is the material derivative.

Its two terms have different meanings. The first measures how the velocity changes with time at a fixed location. The second, the advection term (u · ∇)u (also called the convective term), measures the change experienced because the parcel moves through a velocity field that varies from place to place.

With the acceleration identified, Newton’s law can now be applied. For a fluid of density ρ,

Pressure contributes the force per unit volume

directed toward decreasing pressure.

For a Newtonian fluid with constant dynamic viscosity μ, viscous stresses contribute

The Laplacian compares the velocity at a point with the velocity in its surroundings. In this sense, viscosity tends to smooth sharp spatial variations in the flow.

Finally, let f denote an external force per unit mass. Newton’s second law becomes

Dividing by the constant density and introducing the kinematic viscosity

gives the incompressible Navier–Stokes equation:

A second equation expresses incompressibility:

It states that the flow cannot locally compress or expand: for any sufficiently small region, fluid enters and leaves at the same net rate. Pressure then does more than simply push the fluid from high-pressure toward low-pressure regions; it adjusts throughout the flow so that this incompressibility constraint remains satisfied.

The resulting system contains the competition that will matter throughout the rest of the article. Advection transports and deforms the velocity field, while viscosity smooths spatial variations.

The advection term (u · ∇)u also introduces the central mathematical complication: it makes the equation nonlinear. The velocity field appears both as the quantity being differentiated and as the field that determines how that quantity is transported. In other words, the velocity field is itself responsible for transporting the velocity field.

From Velocity to Vorticity

The velocity formulation tells us how the fluid moves, but the mechanism most closely associated with possible singular behavior becomes clearer when we look at how the fluid rotates.

The local rotation of the flow is described by the vorticity

In three dimensions, ω is a vector. Its direction gives the local axis of rotation, while its magnitude measures the strength of that rotation.

To find how vorticity evolves, we take the curl of the Navier–Stokes equation. The pressure term disappears because the curl of a gradient is zero:

For an incompressible flow, the resulting equation is

Using the material derivative introduced in the previous section, this can be written more compactly as

Each term has a distinct physical interpretation.

The material derivative on the left follows the vorticity as it moves with the fluid. The viscous term

acts to diffuse sharp variations in vorticity, just as viscosity smooths variations in velocity.

The crucial new term is

This is the vortex-stretching term.

To see what it means, we can imagine a thin tube of rotating fluid. If the surrounding flow stretches the tube along its axis, conservation of volume forces its cross-section to shrink. At the same time, the rotation can intensify: the vorticity becomes larger as the structure becomes thinner and more elongated.

This creates a possible feedback mechanism. As vorticity intensifies, the associated velocity gradients can stretch the vortex further; that stretching can then increase the vorticity again.

Viscosity acts in the opposite direction, diffusing these increasingly fine structures.

In two-dimensional incompressible flow, the vortex-stretching term disappears, removing an important mechanism for amplifying vorticity. But in three dimensions, that term remains active.

Why Two and Three Dimensions Behave Differently

For a two-dimensional incompressible flow, the velocity has the form

and the vorticity points perpendicular to the plane:

Because the flow has no variation in the third direction,

The vorticity equation therefore reduces to

What has disappeared is the vortex-stretching term capable of amplifying the vorticity through the flow’s own velocity gradients.

The equation is still nonlinear because the velocity field continues to advect the vorticity. What is missing is the self-amplification mechanism. In the unforced case, advection transports vorticity without creating larger extrema, while viscosity diffuses it. As a result, the maximum vorticity can be controlled for all time.

In three dimensions, by contrast, the term

need not vanish. Vortex stretching can increase the magnitude of the vorticity through the flow’s own velocity gradients. Unlike in two dimensions, there is now a mechanism that can create larger vorticity even as viscosity acts to smooth it.

Smooth two-dimensional incompressible Navier–Stokes solutions remain smooth for all time, whereas whether smooth three-dimensional solutions must remain smooth was the main question posed by the Millennium problem.

The mathematical term for this persistence of smoothness is regularity. A solution is regular when it remains sufficiently smooth, with the relevant quantities staying finite rather than developing a singularity.

How Viscosity Dissipates Kinetic Energy

Viscosity smooths spatial variations in the velocity field and continuously dissipates kinetic energy. The question is whether this combination of smoothing and energy dissipation is sufficient to preserve regularity in three dimensions.

This dissipative effect appears directly in the energy balance. Consider first the unforced Navier–Stokes equations, meaning that the external force is set to

Assume also that the velocity field either decays sufficiently rapidly at infinity or satisfies boundary conditions that produce no boundary contribution. Multiplying the Navier–Stokes equation by u and integrating over space shows how the total kinetic energy changes with time.

The time-derivative term gives

Up to the constant fluid density, this integral is the total kinetic energy of the flow:

Two other terms simplify because the flow is incompressible. The nonlinear advection term redistributes kinetic energy through the fluid but, after integration over the whole domain, does not change the total amount. The pressure term likewise gives no net contribution under the assumed boundary conditions.

Unlike advection and pressure, viscosity changes the total kinetic energy of the flow. After integration by parts,

The minus sign is essential: viscosity can only remove kinetic energy from an unforced flow.

Combining the terms gives the classical energy identity

Equivalently,

Hence, the total kinetic energy cannot grow. Moreover, viscosity continuously dissipates energy at a rate determined by the spatial gradients of the velocity field.

For a forced flow, the same calculation adds an energy-input term:

As a result, the forcing term can inject kinetic energy into the flow, while viscosity continues to dissipate it. However, the forced-flow energy balance derived above does not by itself rule out a singularity.

Why Finite Energy Does Not Guarantee Regularity

A bounded L² norm constrains the total kinetic energy, but not the magnitude of the velocity or its gradients at every point in space. Thus, an integral can remain finite even while a function becomes extremely large over an increasingly small region.

Consider a velocity field with a typical magnitude A inside a region of volume V. Its contribution to the kinetic energy is roughly

Suppose that the region containing the large velocity becomes smaller and smaller. If the amplitude grows while the occupied volume shrinks sufficiently quickly, the product A²V can remain bounded.

For example, in three dimensions a structure of characteristic size ℓ occupies a volume of order

If its velocity amplitude grows like

then

So the energy contribution remains bounded, even as the characteristic size of the structure becomes smaller and smaller.

The corresponding velocity gradient behaves differently: if the velocity changes by an amount A across a distance ℓ, then

With the same scaling A ∼ ℓ^{-3/2},

which becomes unbounded as ℓ → 0.

In other words, bounded energy does not prevent increasingly large local gradients: the energy stays of order one, while the gradient grows without bound as ℓ → 0.

Navier–Stokes Scaling and the Regularity Problem

The previous section showed that finite kinetic energy does not by itself prevent large velocities or gradients from developing at increasingly small spatial scales. But how do the Navier–Stokes equations behave under changes of scale?

The equations possess a particular scaling symmetry. Consider again the unforced equations. If u(x,t) and p(x,t) form a solution, then for any λ > 0,

and

also satisfy the Navier–Stokes equations with the same viscosity.

This means that the equations retain the same form when space, time, velocity and pressure are rescaled together in this particular way. This scaling shows how quantities such as velocity, kinetic energy, and Lᵖ norms change as progressively smaller spatial scales are examined.

Increasing λ corresponds to a characteristic length

Under this scaling, time contracts as

while the velocity amplitude grows as

And the L² norm, which determines the kinetic energy, transforms as

So, as the spatial scale becomes smaller, with λ → ∞,

The natural small-scale scaling of the Navier–Stokes equations makes the kinetic energy smaller, not larger.

More generally, the Lᵖ norm transforms as

However, one value of p is special:

Indeed,

Such a quantity is called scale-critical: it remains unchanged under the natural scaling of the equation. By contrast, the L² energy norm becomes smaller as the flow is rescaled toward finer structures. For p > 3, the exponent 1 - 3/p is positive, so the corresponding Lᵖ norms increase under the same small-scale rescaling. In particular,

This is why three-dimensional Navier–Stokes is often described as energy-supercritical. The basic energy estimate bounds the total kinetic energy, but under rescaling toward smaller spatial scales the L² energy decreases even while the characteristic velocity and its gradients increase. Thus, the energy bound does not exclude increasingly intense structures at smaller spatial scales.

The regularity problem is not simply a matter of keeping the total energy finite. Stronger bounds are needed on the velocity and its gradients at the increasingly small scales where a singularity could form.

What the Millennium Problem Actually Asked

The question of whether a smooth three-dimensional flow must remain smooth forever can be stated more precisely in terms of global smooth solutions.

Starting from smooth, divergence-free initial data, a smooth solution is known to exist for a short time. The unresolved question was whether it must remain smooth for all t ≥ 0, or whether smoothness can break down after a finite time.

The official problem statement by Charles Fefferman sets out two possible ways to resolve the problem: prove global smoothness, or exhibit finite-time breakdown. Statements A and B ask for global existence and smoothness for the unforced equations, respectively in the whole space and in a periodic domain. Statements C and D ask instead for smooth initial data and smooth external forcing for which no global smooth solution exists.

The logical structure of the Clay formulation can be summarized as follows.

The distinction between the unforced global-smoothness alternatives A/B and the forced blow-up alternatives C/D matters: a finite-time blow-up result does not need to involve an unforced fluid to satisfy the official Millennium problem formulation. A smooth force producing finite-time loss of regularity is sufficient.

Leray had already established the existence of global weak solutions in 1934. The unresolved issue was whether solutions starting from smooth initial data must remain smooth for all time.

A breakdown of regularity means that the solution cannot remain smooth for all time: as some finite time T is approached, a quantity required for smoothness becomes unbounded.

What the New Blow-Up Result Shows

The announced result takes the second route allowed by the Millennium problem: it exhibits a smooth three-dimensional Navier–Stokes flow whose velocity becomes unbounded after a finite time.

The flow begins from the simplest possible initial state,

The fluid is initially at rest. A smooth external force is then applied, and for every time before the singularity the resulting velocity and pressure fields remain smooth. The force itself is also smooth and compactly supported in space and time. Yet the velocity eventually becomes unbounded while the total kinetic energy remains bounded.

More precisely, the results yields a solution which is smooth for 0 ≤ t < 1, with

but

The two norms describe very different aspects of the flow. The bounded L² norm says that the total kinetic energy remains finite. The diverging L^∞ norm says that the maximum velocity becomes arbitrarily large somewhere in the fluid. Finite total energy therefore coexists with finite-time loss of regularity.

The singular behavior is produced by a shrinking vortex centered near the origin. Writing

for the remaining time before the singularity, the characteristic radial and axial dimensions of the vortex core behave approximately as

where h is a small positive constant. Both dimensions therefore shrink to zero, but the radial scale contracts slightly faster than the axial one. The core becomes progressively thinner and more elongated as t → 1.

At the same time, the dominant velocity components grow approximately as

The velocity diverges as the singular time is approached. Yet the volume occupied by the intense part of the flow shrinks rapidly enough that its kinetic-energy contribution scales like

which tends to zero for the range of h used in the construction. The region in which the velocity becomes arbitrarily large can therefore contribute progressively less to the total kinetic energy.

This geometry is characteristic of three-dimensional vortex dynamics. Fluid spirals inward toward the axis while being carried axially away from the center. Radial inflow increases the azimuthal velocity, while incompressibility drives axial motion away from the contracting core. Viscosity remains active throughout but does not prevent the characteristic velocity from increasing.

The most delicate part of the construction is not simply making such a velocity field large. For any chosen velocity and pressure fields one could formally define an external force that makes the Navier–Stokes equation hold. That would not solve the Millennium problem if the required force itself became singular.

What must be shown is stronger: the velocity becomes singular while the external force remains perfectly smooth.

To achieve this, the construction introduces carefully designed oscillatory perturbations around the shrinking vortex. Their nonlinear interaction compensates for the otherwise singular residual left by the leading vortex flow. Further corrections remove the remaining singular terms, so that the final force extends smoothly through the blow-up time even though the velocity does not.

The construction establishes statement C of the official Clay formulation in ℝ³. Because the relevant fields can also be localized appropriately in the periodic setting, it yields statement D as well.

The proposed result exhibits finite-time blow-up in a smooth forced Navier–Stokes flow while the initial data, external forcing and total kinetic energy remain within the admissible framework of the Millennium problem.

What the Result Does — and Does Not — Settle

The constructed flow is forced. It does not show that a smooth solution of the unforced equations, f = 0, can develop a finite-time singularity.

This still falls within the Millennium problem, since the official statement includes the forced case. But it leaves a mathematically distinct question untouched: whether nonlinear three-dimensional Navier–Stokes dynamics can produce finite-time blow-up without an externally applied force.

The result also does not imply that ordinary physical fluids are expected to reach infinite velocity. A mathematical singularity marks a breakdown of the smooth solution of the idealized equations. At sufficiently small physical scales, assumptions underlying the continuum model may themselves cease to provide an adequate description of a real fluid. The theorem concerns what the Navier–Stokes equations permit mathematically, not a prediction that an experimentally observable fluid must literally attain infinite speed.

OpenAI has released both an analytical proof and a Lean formalization. The Clay Mathematics Institute has stated that the problem has “apparently been settled,” while emphasizing that the process of evaluating the result and assigning credit will be deliberately unhurried.

Conclusion

The Navier–Stokes equations describe a balance between inertia, pressure and viscosity: inertia transports the flow, pressure enforces incompressibility, and viscosity smooths velocity variations and dissipates kinetic energy.

The difficulty appears in three dimensions. Vortex stretching can amplify vorticity and generate increasingly fine spatial structure. The classical energy estimate shows that the total kinetic energy cannot increase in an unforced flow, but it does not prevent the velocity or its gradients from becoming very large in a sufficiently small region. The scaling of the equations makes the limitation sharper: the L² energy decreases under small-scale rescaling even while characteristic velocities and gradients increase.

In the proposed solution, a smooth, initially stationary fluid subjected to smooth forcing develops a shrinking vortex whose velocity becomes unbounded in finite time while its total energy remains finite. If the proof is confirmed, it shows that viscous dissipation does not guarantee global smoothness: the three-dimensional Navier–Stokes equations permit finite-time breakdown under the forced conditions included in the official Millennium problem.

The Python code used to generate all figures and diagrams in this article is open-source and available in the GitHub companion repository.

References & Further Reading

1. Leray, J. (1934). Sur le mouvement d’un liquide visqueux emplissant l’espace. Acta Mathematica, 63, 193–248. 10.1007/BF02547354 · Springer.
2. Caffarelli, L., Kohn, R., & Nirenberg, L. (1982). Partial regularity of suitable weak solutions of the Navier–Stokes equations. Communications on Pure and Applied Mathematics, 35(6), 771–831. 10.1002/cpa.3160350604 · Wiley.
3. Escauriaza, L., Seregin, G. A., & Šverák, V. (2003). L³,∞-solutions of the Navier–Stokes equations and backward uniqueness. Russian Mathematical Surveys, 58(2), 211–250. 10.1070/RM2003v058n02ABEH000609 · IOPscience.
4. Fefferman, C. L. (2006). Existence and Smoothness of the Navier–Stokes Equation. In J. Carlson, A. Jaffe & A. Wiles (Eds.), The Millennium Prize Problems, 57–67. Clay Mathematics Institute / American Mathematical Society. Official Clay PDF.
5. Tao, T. (2016). Finite time blowup for an averaged three-dimensional Navier–Stokes equation. Journal of the American Mathematical Society, 29(3), 601–674. 10.1090/jams/838 · arXiv:1402.0290.
6. OpenAI. (2026). On the Navier–Stokes Millennium Prize Problem. OpenAI Research, 8 September 2026. OpenAI Research Announcement.
7. OpenAI. (2026). Finite Time Blowup for Navier–Stokes. Research preprint, September 2026. OpenAI Proof PDF.
8. Clay Mathematics Institute. (2026). Navier–Stokes Announcement. 11 September 2026. Clay Mathematics Institute Announcement.

Navier–Stokes Blow-Up: The Mathematics Behind the Millennium Problem was originally published in Stellar Priors on Medium, where people are continuing the conversation by highlighting and responding to this story.

By the second half of the eighteenth century, Newtonian gravity had provided a unified description of planetary motion. The motions of the planets could be understood as variations on Kepler’s ellipses, all governed by the same inverse-square law.

The two-body picture, however, omits the mutual gravitational influence of the planets. The Solar System is a many-body system. Jupiter pulls on Saturn, Saturn pulls on Jupiter, and both perturb the other planets. However small those effects may be over one orbit, they have centuries in which to accumulate.

Astronomers had observed that the planets did not return exactly to the positions predicted by simple Keplerian motion, raising the question of what caused these discrepancies. Were the orbital elements merely oscillating around stable values, or could the perturbations slowly build up until the long-term configuration of the Solar System itself changed?

The mathematician Joseph-Louis Lagrange came to the question through his earlier work in celestial mechanics. During the 1760s he had already studied the motion of the Moon, the satellites of Jupiter, and the perturbations of Jupiter and Saturn, problems in which small gravitational interactions produced observable departures from simple Keplerian motion.

In the 1770s, Lagrange began giving the long-term evolution of planetary orbits a mathematical form. He studied how the inclinations, nodes, eccentricities and other elements of planetary orbits change under their mutual gravitational attraction. Rather than attempting to calculate every small displacement directly, he followed the slow evolution of the orbits themselves. In 1774, his work on the secular variations of planetary inclinations and nodes reached the Paris Academy.

Pierre-Simon Laplace extended this approach to other orbital elements and to the coupled evolution of the planetary system as a whole. His work increasingly focused on the secular effects produced by the mutual gravitational interactions of the planets.

Observations of the motions of Jupiter and Saturn showed slow discrepancies that had eluded earlier calculations. Lagrange’s work suggested that a simple secular drift could not account for the observed behavior.

Laplace would eventually show that the strange long-term behavior of Jupiter and Saturn was tied to their near 5:2 resonance: Jupiter completes almost five orbits in the time Saturn completes two. What looked like a persistent drift was instead part of a very long oscillation.

Laplace’s work resolved one of the most prominent anomalies in planetary motion, but it did not settle the more general question of whether the mutual perturbations of all the planets could remain bounded over very long times.

If every planet continuously perturbs every other planet, and those perturbations accumulate over centuries and millennia, is the Solar System actually stable?

From Keplerian Orbits to Perturbed Motion

For an isolated planet orbiting the Sun, Newtonian gravity reduces the motion to the two-body problem. Its geometry can be described by a small set of orbital elements. The semimajor axis a sets the size of the orbit, the eccentricity e its shape, and the inclination i the tilt of its orbital plane. Additional angular elements specify the orientation of that plane and the direction of perihelion.

A real planetary trajectory is not one fixed ellipse repeated forever. However, at any instant, one can define the Keplerian orbit that the planet would follow if all perturbing forces were suddenly removed while its position and velocity remained unchanged. This instantaneous ellipse is called the osculating orbit.

As the perturbing forces gradually alter the planet’s motion, the corresponding osculating ellipse changes with it. The orbit can therefore still be described in terms of Keplerian elements, but those elements are no longer constants

while the angular elements specifying the orientation of the orbit also evolve with time.

This description separates the motion into two distinct parts. The planet moves rapidly around the Sun on an approximately Keplerian orbit, while the parameters describing that orbit evolve much more slowly.

A planet may complete a revolution in months or years, whereas the orientation of its perihelion, the tilt of its orbital plane, or its eccentricity may vary over thousands or even hundreds of thousands of years. The main problem of perturbation theory is not to replace Keplerian motion with an entirely different description, but to determine how the Keplerian orbit itself changes under the weak influence of the other planets.

What mattered for the long-term motion was whether these small perturbations cancelled over repeated orbits or accumulated over time. A perturbation that reverses sign as the planets move around their orbits may largely cancel over many revolutions, producing only a bounded oscillation. Another contribution may survive the averaging over orbital motion and act coherently over much longer timescales. Such slowly accumulating effects came to be known as secular perturbations.

This time-scale distinction became closely tied to the question of stability. If a planet’s eccentricity merely oscillates between nearby values, its orbit remains confined to a limited range of shapes. If, however, a secular term were to drive the eccentricity steadily upward, an initially nearly circular orbit could eventually become highly elongated, and neighboring planetary orbits might begin to intersect. The deviations from Keplerian motion were therefore not merely corrections required for accurate astronomical tables.

These deviations from Keplerian motion cast doubt on the long-term stability of the planetary system itself: were the observed variations bounded, or could the cumulative effect of mutual gravitational perturbations eventually drive the planets far from their present orbital configuration?

Fortunately, the Solar System contains a natural hierarchy. The Sun overwhelmingly dominates the motion of each planet, while the planets exert much weaker forces on one another. This separation of scales is what makes a perturbative treatment possible: the dominant Keplerian motion can be solved first, while the weaker planetary interactions are introduced as corrections to it.

Perturbing a Keplerian Orbit

To make the separation between the dominant Keplerian motion and the weaker planetary interactions explicit, we first consider only the Sun and a single planet, and neglect the planet’s mass compared with the Sun’s. The planet’s heliocentric position vector r obeys

where r=∣r∣ is the distance from the Sun and M⊙ is the solar mass. This is the exactly solvable Keplerian problem that provides the unperturbed motion.

Once the other planets are included, additional gravitational terms appear, but for the major planets these remain small compared with the central attraction of the Sun. So, the full problem can be viewed schematically as

In Hamiltonian mechanics, the same separation can be written more compactly as

The Hamiltonian H is the function that describes the dynamics of the system in terms of its positions and momenta. Here H₀ represents the exactly solvable Keplerian two-body problem, while H₁ contains the weaker gravitational interactions between the planets. The parameter ε measures the relative strength of those perturbations.

In celestial mechanics, the planetary perturbations are represented by the disturbing function R. With a common sign convention,

so R represents the same perturbing interactions as εH₁ (up to the sign convention), written in a form suited to studying how the orbital elements evolve. In heliocentric coordinates, it can be viewed schematically as

where the direct term describes the gravitational attraction of the perturbing planet, while the indirect term accounts for the acceleration of the Sun itself.

The indirect term arises because heliocentric coordinates use the Sun as their origin, even though the Sun itself is accelerated by the other planets. Expressed in terms of the orbital elements and relative geometry of the planets, the disturbing function provides a systematic way to determine how the semimajor axis, eccentricity, inclination and orbital orientation evolve.

The disturbing function rewrites this additional gravitational interaction in terms of the orbital elements and relative geometry of the interacting planets, allowing their effects on semimajor axis, eccentricity, inclination and orbital orientation to be calculated systematically.

When expanded, the disturbing function contains many terms involving combinations of orbital angles, eccentricities and inclinations, reflecting the different geometrical configurations that recur as the planets move around the Sun.

These terms do not all matter in the same way over long timescales. Many oscillate rapidly as the relative positions of the planets change and largely cancel over repeated revolutions; others vary much more slowly and can shape the long-term evolution of the orbits.

From Periodic Perturbations to Secular Motion

To separate the rapidly oscillating perturbations from the slower variations that persist over many orbital periods, one can average over the rapidly changing orbital phases. Terms that alternate in sign as the planets move around the Sun largely cancel over many revolutions, while more slowly varying contributions remain.

What survives this averaging is the secular evolution: the much slower change in the shapes and orientations of the planetary orbits.

The separation of timescales is large: the mean longitude of a planet advances by 2π during every revolution, while its perihelion and orbital plane may evolve over tens or hundreds of thousands of years. This gap is precisely what makes orbital averaging useful: the fast variables describe where the planets are on their orbits, while the slow variables describe how the orbits themselves change.

The simplest secular approximation removes the dependence on the rapidly varying orbital phases. What remains is an averaged interaction between the planetary orbits themselves. In this picture, a planet is no longer treated primarily as a point moving around its ellipse, but almost as if its mass were spread continuously along that orbit. The rapid motion has been averaged away, leaving the much slower gravitational coupling between the ellipses.

The averaged equations no longer describe where each planet is along its orbit, but how the orbit itself evolves over long timescales. The semimajor axes are approximately constant at this level of approximation, while the eccentricities, inclinations and orbital orientations continue to evolve. The quantities of interest are therefore not primarily the orbital periods themselves, but the way the shapes and planes of the orbits precess and exchange angular momentum.

For the eccentricity dynamics, a convenient pair of variables is

and

where e is the eccentricity and ϖ is the longitude of perihelion.

The quantities h and k then represent the eccentricity and perihelion orientation as a single vector in the (k,h) plane.

A similar pair of variables can be introduced for the inclination and the orientation of the orbital plane. At this linear order, the inclination dynamics can be treated in much the same way, so the discussion below can be restricted to the eccentricity variables.

After averaging and retaining only the lowest relevant orders in eccentricity and inclination, the secular equations become linear. What originally looked like a complicated nonlinear gravitational interaction between several moving planets is reduced to the evolution of a finite set of slowly varying orbital elements.

If the point (k,h) associated with a planet moves around the origin, the eccentricity and longitude of perihelion both change. If its distance from the origin remains bounded, so does the eccentricity. The question of whether an orbit becomes progressively more elongated can therefore be translated into a question about the solutions of the averaged secular equations.

A secular term proportional to time would represent growth without bound within the approximation, whereas a sum of oscillatory terms would keep the motion confined to a finite range.

However, considering the secular evolution of one planet in isolation would not be realistic. Jupiter affects Saturn, Saturn feeds back on Jupiter, and both interact more weakly with the rest of the planetary system. The eccentricity variables of the planets are therefore dynamically coupled.

Lagrange–Laplace Theory: The Solar System as a Set of Coupled Modes

The theory now known as Lagrange–Laplace secular theory emerged from a sequence of closely connected investigations rather than a single joint calculation. Building on Lagrange’s work on the secular variation of orbital elements, Laplace extended the analysis and helped develop a general description of the slow, coupled evolution of the planetary system.

After averaging over the rapid orbital phases and retaining the lowest relevant orders in eccentricity and inclination, the calculation leads to a set of linear coupled differential equations. For the eccentricity variables introduced above, these can be written in modern notation as

and

where j labels the planet and the matrix A describes the secular coupling between the planetary orbits. Its coefficients depend principally on the planetary masses and semimajor axes. The important point is not their detailed form, but the fact that the eccentricity of each planet is coupled to those of the others.

The structure is analogous to a system of coupled oscillators. Pendulums connected by weak springs do not generally oscillate independently; the collective system possesses normal modes, characteristic patterns of motion in which all parts of the system oscillate together with a single frequency. An arbitrary motion is then a superposition of those modes.

The coupled equations take a more compact form when the two real eccentricity variables are combined into the complex variable

The two real equations then combine into

or, in matrix notation,

In modern linear-algebra language, the problem is solved by diagonalizing A. If its eigenvalues are g_1, g_2, ..., each eigenvector represents a collective secular mode evolving as

and the motion of any individual planet is a superposition of these modes,

with coefficients determined by the initial state of the system.

This result shows a planet does not possess an independent secular oscillation of its own; its long-term eccentricity and perihelion motion contain contributions from collective modes of the planetary system. The same construction applies to inclinations and nodal longitudes, producing another family of secular frequencies.

For stability, the relevant feature is that the linear solutions are oscillatory. In the Laplace–Lagrange approximation, conservation of angular momentum implies that the secular frequencies are real, provided the planets revolve around the Sun in the same direction. The factors exp(i g_p t) remain bounded: eccentricities and orbital orientations can vary continually without containing a term that simply grows without limit.

A striking historical example of how a very slow oscillation could masquerade as a secular drift was already visible in the motions of Jupiter and Saturn.

Jupiter and Saturn: Explaining the Great Inequality

The absence of a secular drift in the planetary semimajor axes found in the calculations of Lagrange and Laplace seemed to conflict with a problem astronomers had known for generations: the slow, unexplained variations in the motions of Jupiter and Saturn that would later become known as the Great Inequality of Jupiter and Saturn.

By comparing ancient observations with those of his own time, Edmond Halley concluded that Jupiter’s mean motion appeared to be gradually increasing while Saturn’s appeared to be decreasing. He incorporated empirical secular corrections into his planetary tables. Later calculations by Leonhard Euler and Lagrange attempted to derive these variations from the mutual gravitational interaction of the two planets, but their results did not resolve the discrepancy. In 1773, Johann Heinrich Lambert also found that more recent observations seemed to show the trend reversing, strengthening the possibility that the apparent drift belonged to a much longer periodic variation rather than an indefinitely growing secular one.

This was the problem Laplace returned to in the mid-1780s. His work with the secular perturbation equations had given him a reason to suspect that the observed drift was not secular at all. If the known first-order secular terms could not explain it, perhaps an extremely slow periodic perturbation could.

The clue was the orbital periods of Jupiter and Saturn. They lie close to a 5:2 commensurability: Jupiter completes almost five revolutions while Saturn completes two. Let λ_J and λ_S denote their mean longitudes, angular variables that track the planets’ average progress around their orbits. The combination

is therefore a slow angular variable measuring the gradual phase mismatch associated with the near 5:2 commensurability.

If n_J and n_Sdenote the planets’ mean motions (their mean orbital angular frequencies), then, to leading Keplerian order,

and hence

Had the orbital periods been in an exact 5:2 ratio, this combination would vanish. Because the commensurability is only approximate, 2n_J − 5n_S is small but non-zero. The angle φ therefore evolves on a timescale far longer than either orbital period, allowing associated perturbative terms to remain coherent for centuries before reversing.

The result resolved the apparent contradiction. Observations separated by centuries had sampled only portions of a very long oscillation, making a bounded periodic variation look like a persistent acceleration or deceleration. Laplace could account for ancient and modern observations within Newtonian gravity without introducing an irreversible secular drift.

This analysis also exposed an important limitation of simple orbital averaging. Most rapidly varying angular combinations disappear when averaged over many revolutions, but near a resonance some combinations evolve slowly enough to survive. A perturbation need not be large to have an important long-term effect; it may instead vary extremely slowly.

The Great Inequality strengthened the case for long-term order in the Solar System, while also showing that resonances could complicate the simple secular picture. Explaining the apparent drift of Jupiter and Saturn was not the same as proving that the planetary system as a whole was stable.

Had Lagrange and Laplace Proved the Solar System Stable?

The linear secular approximation yielded a decomposition of the slow evolution of eccentricities, inclinations and orbital orientations into regular modes after averaging over the rapid orbital motion. The semimajor axes showed no secular drift, while the remaining elements underwent bounded oscillations rather than simple runaway growth.

But this did not prove that the real Solar System would remain stable indefinitely. It showed stability only within the simplified, lowest-order model used by Lagrange and Laplace.

Linear Lagrange–Laplace theory assumes that the planetary masses are small compared with the Sun and retains only the lowest relevant orders in eccentricity and inclination. Higher-order terms introduce additional couplings, while near-resonant combinations of frequencies can generate small denominators and effects absent from the linearized system. Nor does the success of a finite-order approximation guarantee that the full perturbation series converges indefinitely.

The word stable also hides several different questions. Do the semimajor axes remain close to their present values? Do the eccentricities remain bounded? Do the planets avoid close encounters? Or do nearby initial conditions continue to produce nearby trajectories? These properties are related, but they are not equivalent. A planetary system may remain physically intact while its detailed evolution is chaotic.

Lagrange–Laplace theory established an important but limited form of stability. Within its approximation, planetary perturbations do not simply accumulate without bound; much of the long-term evolution can instead be understood as a superposition of regular secular modes. That theory remains fundamental to modern celestial mechanics, even though accurate descriptions of the Solar System must also account for nonlinear interactions and resonances.

During the nineteenth century, attempts to push perturbation theory to higher orders made those limitations increasingly difficult to ignore. The calculations became more intricate, and the question shifted from finding the next correction to asking whether the method itself could establish the long-term stability of the planetary system.

By the end of the century, Henri Poincaré would show that the problem was more fundamental. In studying the three-body problem, he found that nonlinear gravitational dynamics could possess an intrinsic complexity and sensitivity that no finite collection of corrected Keplerian ellipses could fully capture.

Lagrange–Laplace was only an approximation to the full nonlinear dynamics of the planetary system.

Further Reading on Related Topics
From Newton to Hamilton: The Road to Modern Dynamics
Dynamical Systems: The Geometry Behind Differential Equations

The Python code used to generate all figures and diagrams in this article is open-source and available in the GitHub companion repository.

References & Further Reading

1. Halley, E. (1749). Tabulae astronomicae. Accedunt de usu tabularum praecepta. Edited posthumously by J. Bevis. London: William Innys. Google Books · Wikimedia Commons, 1752 English/Latin edition.
2. Euler, L. (1749). Recherches sur la question des inégalités du mouvement de Saturne et de Jupiter. Pièces qui ont remporté le prix de l’Académie royale des Sciences, Prix de 1748, 1–123. Euler Archive.
3. Lagrange, J.-L. (1766). Solution de différents problèmes de calcul intégral. Mémoires de l’Académie royale des Sciences de Turin pour les années 1762–1765, III, 179–380. Reprinted in Œuvres de Lagrange, Tome I, 471–668. Archive.org.
4. Lagrange, J.-L. (1774). Recherches sur les équations séculaires des mouvemens des nœuds, et des inclinaisons des orbites des planètes. Académie Royale des Sciences, 78 pp. Google Books.
5. Lambert, J. H. (1775). Résultat des recherches sur les irrégularités du mouvement de Saturne et de Jupiter. Nouveaux Mémoires de l’Académie royale des Sciences et Belles-Lettres de Berlin, année 1773, 216–221. Google Books · BBAW.
6. Lagrange, J.-L. (1781). Théorie des variations séculaires des éléments des Planètes. Première Partie. Nouveaux Mémoires de l’Académie royale des Sciences et Belles-Lettres de Berlin, 125–207. Reprinted in Œuvres de Lagrange, Tome V. Wikisource.
7. Laplace, P.-S. (1787). Mémoire sur les inégalités séculaires des planètes et des satellites. Mémoires de l’Académie royale des Sciences de Paris, année 1784, 1–50. Reprinted in Œuvres complètes de Laplace, Tome XI, 49–92. Wikisource · BnF.
8. Laplace, P.-S. (1788). Théorie de Jupiter et de Saturne. Mémoires de l’Académie royale des Sciences de Paris, année 1785, 33–160; continued in année 1786, 201–234. Reprinted in Œuvres complètes de Laplace, Tome XI, 95–239. Wikisource.
9. Laplace, P.-S. (1798). Exposition du système du monde (2nd ed.). Paris: Duprat. See especially Livre IV, De la théorie de la pesanteur universelle, and the discussion of planetary perturbations. Wikisource · Gallica.
10. Brouwer, D., & Clemence, G. M. (1961). Methods of Celestial Mechanics. Academic Press. 10.1016/C2013–0–08160–1 · ScienceDirect.
11. Murray, C. D., & Dermott, S. F. (1999). Solar System Dynamics. Cambridge University Press. 10.1017/CBO9781139174817 · Cambridge University Press · Companion site.
12. Laskar, J. (1989). A numerical experiment on the chaotic behaviour of the Solar System. Nature, 338, 237–238. 10.1038/338237a0 · ADS.
13. Laskar, J. (1990). The chaotic motion of the Solar System: A numerical estimate of the size of the chaotic zones. Icarus, 88(2), 266–291. 10.1016/0019–1035(90)90084-M · ADS.
14. Sussman, G. J., & Wisdom, J. (1992). Chaotic Evolution of the Solar System. Science, 257(5066), 56–62. 10.1126/science.257.5066.56 · ADS.
15. Laskar, J. (2008). Chaotic diffusion in the Solar System. Icarus, 196(1), 1–15. 10.1016/j.icarus.2008.02.017 · ADS · arXiv:0802.3371.
16. Laskar, J., & Gastineau, M. (2009). Existence of collisional trajectories of Mercury, Mars and Venus with the Earth. Nature, 459, 817–819. 10.1038/nature08096 · ADS.
17. Laskar, J. (2013). Is the Solar System Stable? In B. Duplantier, S. Nonnenmacher & V. Rivasseau (eds.), Chaos: Poincaré Seminar 2010, Progress in Mathematical Physics, 66, 239–270. Birkhäuser. 10.1007/978–3–0348–0697–8_7 · arXiv:1209.5996.
18. Laskar, J. (2024). Stability of the Solar System. Scholarpedia, 19(4), 5216. 10.4249/scholarpedia.5216 · Scholarpedia.

Is the Solar System Really Stable? Lagrange and Laplace’s Discovery was originally published in Stellar Priors on Medium, where people are continuing the conversation by highlighting and responding to this story.

In 1974, Russell Hulse discovered a pulsar with some peculiar behavior. Its apparently regular pulse period was changing in a systematic but initially puzzling way.

Over the following years, PSR B1913+16 would become one of the most important systems for testing General Relativity. Its importance came not only from the results themselves, but from how much physics could be extracted from the arrival times of its radio pulses.

A Puzzling Pulsar

In the summer of 1974, Russell Hulse was searching for pulsars with the 305-meter radio telescope at Arecibo Observatory in Puerto Rico. The survey formed part of his PhD work at the University of Massachusetts Amherst under the supervision of Prof. Joseph Taylor.

Hulse spent long, irregular shifts inside the frigid hum of the computer control room, sheltered from the tropical heat. The purpose of the survey was to search the Galactic plane for previously unknown pulsars and enlarge the small population then available for study.

On July 2, Hulse spotted the initial weak signal of PSR B1913+16 during a routine survey run. The source produced a radio pulse roughly every 59 ms.

Such regularity is characteristic of a pulsar. A neutron star can contain more mass than the Sun compressed into an object only tens of kilometers across. If its magnetic axis is tilted relative to its rotation axis, radiation emitted near the magnetic poles sweeps through space like the beam of a lighthouse. Whenever that beam crosses Earth, a radio telescope records a pulse. This makes pulsars behave like cosmic clocks.

On August 25, Hulse attempted to obtain a more accurate period for PSR B1913+16. Instead, two measurements taken during the same observing session differed by 27 μs, far too much to dismiss as experimental uncertainty.

His first reaction was to suspect a bug in his software or an equipment error. A second attempt two days later produced an even worse disagreement. Hulse then checked his code, probed the equipment for possible faults, and developed a special observing procedure for PSR B1913+16.

But the problem persisted. Sometimes the pulses arrived a little sooner than expected; at other times they arrived a little later. The changes were too large and too systematic to dismiss as irregularities in the neutron star’s rotation.

Comparing observations from September 1 and 2, Hulse noticed a pattern. When he plotted the measured pulse period as a function of time, the two curves nearly overlapped if one was shifted by about 45 minutes.
And then the explanation came to him: the variations in the observed pulse period were Doppler shifts produced by the pulsar’s orbital motion around a companion star.

When PSR B1913+16 was moving toward Earth, successive pulses arrived slightly closer together, producing a shorter observed period. When it was moving away, they arrived farther apart, producing a longer observed period. And because the Doppler pattern repeated within only a few hours, the pulsar had to be moving in a very tight orbit.

The decisive confirmation that the system was a binary came on September 16: the period reached a clear minimum and then began increasing again, as the binary interpretation predicted.

By September 18, Hulse had completed the analysis and wrote Joseph Taylor with the news that PSR B1913+16 was in a high-velocity binary orbit with a period of about eight hours.

The discovery of the companion solved the immediate puzzle. But the changing pulse period contained much more than evidence for a binary: it carried information about the motion and geometry of the orbit itself.

Reconstructing the Binary Orbit from Pulsar Timing

Once the Doppler origin of the changing pulse period was understood, the timing variations could, in principle, be used to reconstruct the binary orbit.

As a first approximation, the system PSR B1913+16 can be treated as a Keplerian binary: two bodies moving under their mutual gravitational attraction according to Newtonian gravity. Each star follows an elliptical orbit around their common center of mass. Relativistic effects introduce much smaller corrections to this underlying motion, which we will consider later.

For such a Keplerian binary, the radial velocity of the observed star along our line of sight can be written as

Here, e is the orbital eccentricity, with e = 0 for a circular orbit and larger values corresponding to more elongated ellipses.

The angle ν is the true anomaly, the instantaneous angular position of the pulsar measured from periastron, the point of closest approach to its companion. The angle ω, called the argument of periastron, specifies the orientation of periastron within the orbital plane, while K sets the overall scale of the radial velocity.

The first quantity to emerge directly from the repeating Doppler pattern is the orbital period:

The shape of the observed Doppler modulation already showed that the orbit was not circular. Fitting it with the Keplerian model gives an eccentricity of

With this eccentricity, the pulsar spends much of each orbit relatively far from its companion, moving more slowly, before accelerating as it approaches periastron. This is Kepler’s second law made visible through radio timing: the line joining the two stars sweeps out equal areas in equal times, so the orbital motion becomes faster as the separation decreases.

Consequently, the true anomaly ν does not increase uniformly with time. Near periastron it changes rapidly, while farther away it changes more slowly. This produces the characteristic asymmetry of the observed Doppler curve, with line-of-sight velocities reaching roughly 300 km/s.

The amplitude of the Doppler curve contains information about the size of the pulsar’s orbit. The velocity scale K is

where P_b is the orbital period, a₁ the semimajor axis of the pulsar’s orbit around the center of mass, and i the inclination of the orbital plane to our line of sight.

The Doppler effect measures only motion along the line of sight, thus a₁ and i cannot be determined separately. Timing measures their combination

known as the projected semimajor axis. It is conventionally expressed as a light-travel time,

and for PSR B1913+16,

Thus the projected semimajor axis corresponds to about 2.342 light-seconds, or roughly 700,000 km. The changing orbital position shifts the pulse arrival times by amounts of order seconds, making the orbital motion readily measurable.

One final parameter, T₀, specifies the reference time at which the pulsar passes through periastron.

The Keplerian timing orbit can be summarized by five measured parameters:

Together, P_b, e, x, ω, and T₀ reconstruct the ordinary Keplerian orbit from the arrival times of the pulsar’s radio pulses, even though the two stars cannot be spatially resolved.

Relativistic Timing and the Neutron-Star Masses

The Keplerian solution described the dominant orbital motion of PSR B1913+16 very well. But as the pulse arrival times were measured with increasing precision, small systematic departures from that model became detectable.

The Keplerian description remained the underlying orbit, with General Relativity predicting additional, much smaller effects superimposed on that motion.

In pulsar timing, these are described by Post-Keplerian parameters, directly measurable quantities that characterize departures from the ordinary Keplerian orbit. In General Relativity, once the Keplerian parameters are known, their values depend on the masses of the two stars, m₁ and m₂.

The relativistic information did not emerge from a single observing run. After the 1974 discovery, Hulse and Taylor continued timing PSR B1913+16 at Arecibo. By 1976, after more than a year of observations spanning about 390 days, Taylor, Hulse and collaborators had already detected a slow change in the orientation of the orbit.

Hulse then moved on to other work, while Taylor and collaborators continued regular precision timing of the system at Arecibo. Over the following several years, the longer observational baseline made progressively weaker deviations from the Keplerian orbit measurable. By the early 1980s, the accumulated timing data were precise enough to use these relativistic effects to constrain the masses of the two neutron stars.

The first of these effects to become clearly measurable was a slow change in the orientation of the orbit.

Periastron Advance

In Newtonian two-body gravity, an isolated elliptical orbit retains the same orientation from one revolution to the next.

For PSR B1913+16, however, the point of closest approach, the periastron, slowly moves forward around the orbit. This relativistic effect is known as the advance of periastron. The measured rate is

The expression obtained from General Relativity is,

where

is the solar mass expressed in convenient time units.

With the orbital period P_b and eccentricity e already known from the Keplerian timing solution, the remaining mass dependence is essentially through the sum m₁ + m₂.

The measured periastron advance thus determines the total mass of the binary:

But many different pairs of values for m₁ and m₂ could produce the same total mass. Therefore, a second relativistic effect is needed to separate them, one that depends on m₁ and m₂ differently.

The Einstein Delay

In the highly eccentric orbit of PSR B1913+16, both the pulsar’s speed and its distance from the companion vary significantly during each revolution. Near periastron, it moves faster and lies deeper in the companion’s gravitational field.

Both effects alter the rate at which the pulsar’s clock appears to run. Special-relativistic time dilation increases with orbital speed, while gravitational redshift changes as the pulsar moves through the companion’s gravitational potential.

In pulsar timing, these two periodic contributions are combined into the Einstein delay, whose amplitude is denoted by γ. Here, γ is a timing parameter measured in seconds; it should not be confused with the Lorentz factor of special relativity.

In General Relativity,

which gives for PSR B1913+16

An amplitude of `4.29 ms` is readily measurable against a pulsar period of only `59 ms`, particularly when the signal is accumulated over many years of observations.

The different dependence of the Einstein delay γ on m₁ and m₂ supplies the second independent mass constraint.

Determining the Component Masses

The resulting mass constraints can be represented geometrically in the m₁ — m₂ plane: the measured periastron advance defines one curve containing all pairs of masses compatible with ω̇. The Einstein delay defines another.

If General Relativity describes the system correctly, the actual neutron-star masses must lie at their intersection.

By 1982, using pulse-arrival measurements extending from 1974 to 1981, Joseph Taylor and Joel Weisberg had accumulated enough data to solve for the two component masses, obtaining values of roughly 1.42 M☉ and 1.41 M☉ with the precision then available. Continued Arecibo observations over the following decades refined those values to much higher precision.

A later thirty-year timing solution gives

That analysis used arrival times accumulated over decades; by then individual five-minute observations could yield pulse arrival times at the level of tens of microseconds.

The first two relativistic measurements had already determined m₁ and m₂. This meant that any additional observable became a genuine test of General Relativity: its value had to be predicted from those already measured masses.

Orbital Decay as a Test of General Relativity

Over successive observations, a small but systematic trend emerged. The interval between successive orbits was becoming shorter: the 7.75-hour orbital period was gradually decreasing.

Measuring a Changing Orbital Period

If the orbital period were perfectly constant, then after measuring one periastron passage at time T₀, the subsequent passages would occur at

After N revolutions, the pulsar would therefore be expected to return to periastron at these predicted times.

Over time, however, the observed periastron passages began to occur slightly earlier than this constant-period model predicted.

As years of timing measurements accumulated, PSR B1913+16 was found to reach periastron progressively earlier than predicted by a constant-period orbit. The effect is extremely small from one revolution to the next, but the binary completes more than a thousand orbits every year. In timing terms, this systematic drift corresponds to a negative orbital-period derivative,

meaning that the orbital period is gradually becoming shorter.

This tiny change becomes measurable because it accumulates. To leading order, the difference between the observed periastron time and that expected for a constant-period orbit grows approximately as

The displacement scales as N², making even a very small change from one orbit to the next measurable after thousands of revolutions.

A Prediction with the Masses Already Fixed

General Relativity predicts a specific value of dP_b/dt once the properties of the binary are known.

The Keplerian timing solution has already provided P_b and e, while the periastron advance and Einstein delay have independently fixed m₁ and m₂. General Relativity therefore makes a definite prediction:

The measured orbital decay can then be compared directly with this prediction.
Unlike the first two Post-Keplerian parameters, dP_b/dt is no longer needed to determine the masses: it provides an independent test of the theory.

The observed decrease in orbital period thus requires a mechanism for removing energy from the orbit.

Gravitational Radiation and Orbital Energy Loss

In Newtonian gravity, an isolated two-body orbit can continue indefinitely while conserving its mechanical energy. General Relativity changes this picture because gravity is not an instantaneous force acting at a distance: it is a dynamical field whose disturbances propagate through spacetime.

Far from the binary, where the gravitational field is weak, the spacetime geometry can be written as a small perturbation of flat spacetime,

where η_μν is the Minkowski metric of flat spacetime and h_μν is a small perturbation produced by the gravitational field.

A time-dependent perturbation h_μν can propagate outward as a gravitational wave.

Why a Binary Radiates

Gravitational radiation does not begin with a changing mass monopole or dipole. Conservation of mass-energy and momentum removes those possibilities, so the first radiative contribution is the mass quadrupole moment.

For a binary, the quadrupole tensor I_ij describes how the two masses are distributed around their common centre of mass. As the stars orbit one another, this mass distribution changes with time.

Far from the system, the gravitational-wave perturbation is proportional to the second time derivative of that quadrupole:

Here, I_ij denotes the mass quadrupole tensor in a normalization that absorbs the conventional factor of 2, R is the distance to the source, and t − R/c is the retarded time. The superscript TT denotes the transverse-traceless part of the perturbation.

As the stars orbit, the quadrupole moment changes continuously, generating gravitational waves that carry energy away from the system.

The Role of the Eccentric Orbit

The highly eccentric orbit of PSR B1913+16 significantly enhances its gravitational-wave emission.

During each revolution, the stars spend much of their time relatively far apart before sweeping rapidly through periastron. The rapid motion near periastron produces the strongest time variation of the quadrupole moment, concentrating much of the gravitational-wave emission into that part of the orbit.

For an eccentric binary, the Peters–Mathews calculation introduces the enhancement factor

For PSR B1913+16, with e ≃ 0.617,

This corresponds to a gravitational-wave energy loss about 12 times greater than for a circular binary with otherwise comparable orbital parameters.

Combining the gravitational radiation calculation with the orbital dynamics gives the predicted secular period change,

The period change measured from Earth also includes small contributions from the binary’s motion within the Galaxy. Correcting for these gives the intrinsic orbital decay.
Arecibo did not directly detect the gravitational waves from PSR B1913+16 passing through Earth. Instead, it measured what their emission was doing to the binary that produced them.

As gravitational waves carried energy away, the orbit gradually changed, yielding the cumulative shifts in pulse arrival times measured at Arecibo. That indirect chain provided the first compelling observational evidence that gravitational radiation carries energy in accordance with General Relativity.

From Indirect Evidence to Gravitational-Wave Astronomy

The Hulse–Taylor observations inferred gravitational radiation from its effect on the binary orbit. Arecibo did not detect a gravitational wave passing through Earth. Instead, it measured the gradual orbital changes produced as the system lost energy:

gravitational radiation → orbital energy loss → orbital decay → pulse timing

After correcting for Galactic and kinematic effects, decades of observations showed that the measured decay followed the relativistic prediction to high precision. Hulse and Taylor were awarded the 1993 Nobel Prize in Physics for the discovery of the binary pulsar, which opened this new way of testing relativistic gravity.

From the Source to the Detector

More than forty years after the discovery of PSR B1913+16, gravitational-wave astronomy entered a different regime.

In 2015, LIGO observed GW150914, produced by the merger of two black holes. This time, the gravitational wave itself was detected directly, through the extremely small differential change in the lengths of the interferometer arms as the wave passed through the detector.

The two experiments probed different parts of the same process: Hulse–Taylor measured the effect of gravitational radiation on the source, whereas LIGO measured the effect of the propagating wave on a distant detector.

binary source → gravitational waves → distant detector

Beyond Hulse–Taylor

Binary-pulsar timing also continued to develop.

The Double Pulsar PSR J0737−3039A/B, discovered in 2003, provided a still richer relativistic laboratory. Its orbital period is only about 2.45 h, both neutron stars were detectable as pulsars, and several independent Post-Keplerian effects could be measured with substantially greater precision.

Pulsar timing has since been extended on an entirely different scale. Pulsar Timing Arrays monitor many stable millisecond pulsars spread across the sky, looking for correlated variations in their pulse arrival times produced by very-low-frequency gravitational waves.

In 2023, several collaborations reported evidence for a common nanohertz gravitational-wave signal consistent with the expected gravitational-wave background.

The techniques are very different, but they share the same basic resource that made PSR B1913+16 so powerful: extremely regular astronomical clocks whose arrival times can be measured over long periods.

The Hulse–Taylor binary therefore occupies an important place between classical pulsar astronomy and modern gravitational-wave astronomy. Long before an interferometer directly measured spacetime strain, the changing orbit of two neutron stars had already shown that gravitational radiation was carrying energy away from an astrophysical system.

The Python code used to generate all figures and diagrams in this article is open-source and available in the GitHub companion repository.

References & Further Reading

1. Peters, P. C., & Mathews, J. (1963). Gravitational Radiation from Point Masses in a Keplerian Orbit. Physical Review, 131(1), 435–440. 10.1103/PhysRev.131.435 · ADS.
2. Peters, P. C. (1964). Gravitational Radiation and the Motion of Two Point Masses. Physical Review, 136(4B), B1224–B1232. 10.1103/PhysRev.136.B1224 · ADS · CaltechTHESIS (PDF).
3. Hulse, R. A., & Taylor, J. H. (1975). Discovery of a pulsar in a binary system. The Astrophysical Journal, 195, L51–L53. 10.1086/181708 · ADS.
4. Taylor, J. H., Hulse, R. A., Fowler, L. A., Gullahorn, G. E., & Rankin, J. M. (1976). Further observations of the binary pulsar PSR 1913+16. The Astrophysical Journal, 206, L53–L58. 10.1086/182131 · ADS.
5. Taylor, J. H., & Weisberg, J. M. (1982). A new test of general relativity — Gravitational radiation and the binary pulsar PSR 1913+16. The Astrophysical Journal, 253, 908–920. 10.1086/159690 · ADS.
6. Damour, T., & Deruelle, N. (1985). General relativistic celestial mechanics of binary systems. I. The post-Newtonian motion. Annales de l’I.H.P. Physique théorique, 43(1), 107–132. ADS · Numdam.
7. Hulse, R. A. (1994). The discovery of the binary pulsar. Reviews of Modern Physics, 66(3), 699–710. 10.1103/RevModPhys.66.699 · ADS · Nobel Lecture (PDF).
8. Taylor, J. H. (1994). Binary Pulsars and Relativistic Gravity. Reviews of Modern Physics, 66(3), 711–719. 10.1103/RevModPhys.66.711 · ADS · Nobel Lecture (PDF).
9. Weisberg, J. M., & Taylor, J. H. (2005). The Relativistic Binary Pulsar B1913+16: Thirty Years of Observations and Analysis. Binary Radio Pulsars, ASP Conference Series, 328, 25–31. ADS · arXiv:astro-ph/0407149.
10. Kramer, M., et al. (2006). Tests of General Relativity from Timing the Double Pulsar PSR J0737–3039A/B. Science, 314(5796), 97–102. 10.1126/science.1132305 · ADS · arXiv:astro-ph/0609417.
11. Damour, T. (2007). Binary Systems as Test-Beds of Gravity Theories. IHES Lecture Series. ADS · arXiv:0704.0749.
12. Abbott, B. P., et al. (LIGO Scientific Collaboration and Virgo Collaboration) (2016). Observation of Gravitational Waves from a Binary Black Hole Merger. Physical Review Letters, 116(6), 061102. 10.1103/PhysRevLett.116.061102 · ADS · arXiv:1602.03837.
13. Agazie, G., et al. (NANOGrav Collaboration) (2023). The NANOGrav 15 yr Data Release: Evidence for a Gravitational-wave Background. The Astrophysical Journal Letters, 951(1), L8. 10.3847/2041–8213/acdac6 · ADS · arXiv:2306.16213.

How Two Neutron Stars Revealed Gravitational Waves Before LIGO was originally published in Stellar Priors on Medium, where people are continuing the conversation by highlighting and responding to this story.

📌 Spatial Interpolation & Geostatistics Series:
• Part 1 (Here): Polynomials & B-Splines
• Part 2: Local Interpolation Models (NN, IDW & RBF)
• Part 3: Geostatistics & Kriging
• Part 4: Spatial Machine Learning

It’s mid-summer and a heatwave has struck Switzerland.

Weather sensors in Geneva record significant spatial variation. Official regional forecasts report a uniform 32°C (90°F), but local microclimates vary substantially across the basin. Urban center stations register 33°C (91°F) from urban heat retention. Atmospheric lapse rates drive lower temperatures across Mont Salève and the Jura foothills, and shoreline sensors along Lac Léman record 29°C (84°F).

So, imagine you live in Geneva and want to know the temperature at your home. How would you proceed?

You might start by looking at the weather map, which shows temperature readings from various weather stations across the region. Then, you would estimate the temperature at your home based on the surrounding weather stations.

This is an example of spatial interpolation, the process of estimating the values of a variable at unsampled locations based on the values of the variable at sampled locations. From everyday tasks like mapping regional weather, geological layers, and hydrological flows to specialized fields like optical engineering and astronomy, the underlying mathematical challenge remains the same.

In this 4-part series, we rigorously benchmark classical spatial interpolation methods alongside modern spatial machine learning models using a single, consistent test case: a simulated spatial model of Geneva, Switzerland.

All algorithms in this article are implemented from first principles using NumPy and benchmarked against production routines in SciPy. All Python code, datasets, and visualizations are open-source and available in the GitHub companion repository.

Note on Notation: To ensure clean, responsive rendering across all mobile and web browsers, this series represents inline variables and multi-indices using clean subscript typography (e.g., z_obs, z_pred) instead of heavy equation blocks.

A Spatial Interpolation Taxonomy

With decades of development spanning diverse mathematical paradigms, spatial interpolation offers a vast array of methods. To avoid getting lost, we classify them across three key properties:

1. Range (Global vs. Local): Do all data points influence every prediction (Global), or only the nearest neighbors (Local)?
2. Exactness (Exact vs. Approximate): Does the surface pass exactly through every observed sensor point (Exact), or does it smooth out noise (Approximate)?
3. Paradigm (Deterministic, Stochastic, or Data-Driven ML): Is the surface built on pure geometric calculations (Deterministic), statistical uncertainty and spatial random processes (Stochastic), or empirical feature mappings (Data-Driven ML)?

In fact, every spatial interpolation method is defined by a specific combination of these three properties. In this series, we cover eleven of the most widely used algorithms, grouping them by shared properties across each part.

In this article, we focus exclusively on the Part 1 interpolation methods. We evaluate how these global coordinate equations model the Geneva simulation. Our goal is to create a reference baseline before exploring local neighborhood methods in Part 2, geostatistics in Part 3, and spatial machine learning in Part 4.

Defining the Baseline: Geneva’s Geography

To get started, we generate a synthetic City Sensor Grid of 100 sensors normalized between coordinate boundaries -1 and 1. The geography contains:

• Lac Léman: A curved crescent-shaped lake where no land-based weather sensors are located. Nearby land experiences lake cooling.
• Jura Mountains & Mont Salève: High elevation zones where temperature decays linearly (the standard temperature lapse rate of -0.0065°C/meter).

📌 Why use Synthetic Data?
In a real-world application, temperatures are only measured at a limited number of scattered weather stations; the exact temperature at unmonitored locations remains unknown. By using a synthetic dataset, we have access to the exact, noiseless generating equation (the Ground Truth) at every point in space. This allows us to calculate absolute errors across the entire 10,000-point grid, verify our interpolation algorithms with complete mathematical accuracy, and visualize exactly where they underperform (e.g. at boundaries or sharp gradients).

Synthetic Field Formulation

We model the thermal structure of the Geneva basin via spatial superposition of three additive physical processes:

1. Regional Macroclimate Gradient (T_macro = 32.0 - 2.0x + 1.2y): Prevailing regional weather trend across the basin.
2. Environmental Lapse Rate (T_lapse = -0.0065 × (Elevation - 375)): Standard thermodynamic cooling of -6.5°C per 1,000m elevation gain over the Jura Mountains (1,175m) and Mont Salève (1,275m).
3. Lac Léman Thermal Buffer (T_lake = -2.5 exp(-d_lake / 0.15)): Microclimate lake breeze cooling at the shoreline that decays exponentially inland.

Synthetic sensor measurements are generated at N = 100 land stations by adding Gaussian noise to the true temperature: z_obs = T_true + ε, where ε ~ N(0, 0.25°C) represents measurement noise.

The map below illustrates the spatial setup, showing discrete noisy temperature readings across the 100 weather stations and the underlying ground truth temperature field:

# (Extracts land points excluding Lake Geneva)
# Computes elevation based on coordinates:
# elevation = 375m (base) + Jura Peak (NW) + Saleve Peak (SE)
# Temp = T_macro_gradient + T_elevation_lapse + T_lake_cooling + noise

(See code/data_generator.py in the GitHub repository).

Now that we have created the synthetic data and mapped the ground-truth temperatures across the Geneva basin, we are in a position to evaluate our first interpolation method.

Bivariate Polynomial Fitting (Trend Surface Analysis)

We start with Bivariate Polynomial Fitting (a process often called Trend Surface Analysis). Because a polynomial fits a single equation to the entire canvas, it is a Global, Approximate, Deterministic method.

We fit a model of the form:

where d is the polynomial degree. The implementation in Python using NumPy is shown below:

def fit_bivariate_polynomial(x, y, z, degree=2):
    columns = []
    for r in range(degree + 1):
        for s in range(degree + 1 - r):
            columns.append((x**r) * (y**s))
    A = np.column_stack(columns)
    coefficients, _, _, _ = np.linalg.lstsq(A, z, rcond=None)
    return coefficients

This method relies on solving the normal equations via Ordinary Least Squares (OLS) to find the coefficients c_{r,s} that minimize the sum of squared errors at our sensor locations.

Before evaluating the performance of global polynomial fitting on the Geneva grid, we consider the second of the global, approximate and deterministic interpolation methods mentioned earlier.

Bivariate B-Splines (Piecewise Polynomials)

Bivariate Basis Splines (commonly called B-Splines) constitute a more flexible alternative to global polynomials. They divide the coordinate grid into smaller rectangular patches (defined by a grid of knots) and fit local low-degree polynomials (typically cubic, p = 3) to each patch, maintaining smooth derivatives across patch borders.

In Python, we fit and predict bivariate splines using SciPy’s bisplrep and bisplev:

from scipy.interpolate import bisplrep, bisplev

# Fit a cubic bivariate B-spline (kx=ky=3)
# 's' is the smoothing factor. If s=0.0, it interpolates exactly (unstable in gaps).
tck = bisplrep(x_obs, y_obs, z_obs, kx=3, ky=3, s=1.0)

# Evaluate on a 1D vector grid
grid_z = bisplev(x_vec, y_vec, tck).T

Unlike global polynomials, B-splines possess local support, meaning that changes in one sensor’s value won’t drastically swing predictions on the opposite side of the map. However, they introduce two additional parameters (hyper-parameters): the knot grid placement and the smoothing factor s. The impact of these parameters will become clear in the next sections.

Visualizing Polynomial and B-Spline Fitted Surfaces

Let’s now look at the results. The side-by-side 2×3 grid below shows the reconstructed temperature fields across global polynomial and B-spline configurations:

If we compare the predicted surfaces above against the ground truth field, we notice several patterns:

• Global Polynomial Smoothness (Degrees 2 & 3): Quadratic and cubic surfaces produce overly simplified, parabolic temperature contours. They capture the broad regional decline from north to south, yet remain far too rigid to depict localized features like Lac Léman’s cool shoreline boundary or mountain peak microclimates.
• Degree 5 Border Oscillations: Increasing the polynomial degree to 5 induces boundary instability along the domain edges. Artificial temperature contours (reaching as low as 15.4°C or spiking unrealistically) wrap around outer coordinates. This illustrates Runge’s phenomenon, where fitting high-degree polynomials over sparse sample points leads to runaway oscillations near domain boundaries.
• B-Spline Local Flexibility & Tuning (s = None vs. s = 1.0 vs. s = 0.0): Piecewise B-splines adapt much better to local curvature. Auto-smoothing (s = None) introduces slight boundary artifacts in data gaps, whereas tuned smoothing (s = 1.0) produces a realistic, continuous temperature gradient. Unregularized exact splines (s = 0.0) become ill-conditioned, generating spike artifacts across unmonitored coordinate gaps.

Let’s now analyze the spatial error distribution of each estimator across the coordinate grid. To this end, we examine their pointwise residual error maps (ẑ - z) below:

• Degree 2 & 3 (Shoreline & Peak Underfitting): Displays dark blue patches over Lac Léman (lake cooling underpredicted by up to -1.4°C) and bright red patches over mountain peaks (overpredicted by +1.9°C), because low-degree polynomials flatten localized temperature features into broad regional averages.
• Degree 5 (Boundary Runge Oscillation): Reduces errors near central sensors close to zero, but introduces concentric residual ring oscillations along outer borders, reaching a -7.58°C underprediction at data-sparse grid corners.
• B-Spline Lake Gap Fluctuations (s = None): Unregularized splines exhibit spurious oscillations in data-sparse regions, yielding an unphysically cold 19.3°C lake prediction.
• B-Spline Regularized Surface (s = 1.0): Suppresses boundary oscillations and reduces domain-wide residual errors to a nearly uniform field (Grid RMSE = 0.322°C).

Performance Comparison: Global Polynomials vs. B-Splines

Looking at the previous maps, you likely guessed that B-splines have an edge over polynomials. Now, let’s confirm this quantitatively across two complementary error metrics:

1. Sensor RMSE (In-Sample): The root mean square error calculated only at the 100 weather station locations.
2. Grid RMSE/MAE (Ground Truth): The root mean square and mean absolute errors calculated across all 10,000 coordinates of the continuous map compared to the simulator’s true underlying temperature equations (where the actual temperature range is [22.95°C, 33.05°C

The results are shown below:

1. Degree 2 Bivariate Fit (Spatial Underfitting)

This model demonstrates spatial underfitting, a condition where a model is too structurally simple to capture the underlying spatial patterns of the data. The small number of parameters in a quadratic surface prevents it from representing the true complexity of the temperature field. It captures the broad regional trend but fails to model localized shoreline cooling and peak microclimates. Mathematically, this high model bias compresses the predicted temperature range ([22.36°C, 33.77°C]), under-representing spatial variance.

2. Degree 3 Bivariate Fit (Persistent Underfitting)

Even with a cubic degree, the model still suffers from underfitting. It yields a very minor improvement over the quadratic baseline, but because a cubic surface is still globally smooth, it cannot model sharp boundaries and localized thermal features (like the crescent lake cooling gradient).

3. Degree 5 Bivariate Fit (Spatial Overfitting & Runge’s Phenomenon)

At degree 5, the model enters a regime of spatial overfitting, where a model fits noisy sample points so closely that it fails to generalize to unsampled areas. While in-sample sensor error decreases from 0.43°C to 0.30°C, the reconstructed surface shows notable spatial degradation. To accommodate localized hotspots, the polynomial oscillates erratically across coordinate gaps. Along data-sparse domain borders, the maximum absolute error reaches 7.58°C, expanding the predicted range down to an unphysical minimum of 15.37°C (compared to true basin minimums near 22.1°C). As noted before, this is the 2D spatial manifestation of Runge’s phenomenon: reducing in-sample error comes at the cost of boundary instability.

4. Bivariate B-Splines (Smoothing Regularization & Tradeoffs)

Unlike rigid global polynomials, mi B-Splines can adapt to localized variations because they are piecewise cubic equations. However, their performance depends entirely on the smoothing factor s:

• Automatic Smoothing (s = None): SciPy’s default routine produces a conservative fit. The surface captures regional trends but introduces localized artifacts near data-sparse boundaries, where temperature estimates decrease to 19.28°C (Grid RMSE = 0.532°C).
• Tuned Smoothing (s = 1.0): Regularization with s = 1.0 provides optimal surface fidelity. Grid RMSE reduces to 0.322°C, accurately modeling lake cooling and terrain lapse rates alongside numerical stability.
• Exact Interpolation (s = 0.0): Zero-smoothing exact fitting (s = 0.0) leads to ill-conditioned system matrices. Across sensor-free regions, basis functions generate runaway oscillations (Grid RMSE > 5,000,000°C), exemplifying multi-dimensional Runge instability.

Bias-Variance Decomposition

We just saw that B-Splines with a s = 1.0 achieved the lowest total prediction error RMSE/MAE. But what actually causes those errors? To find out, we can decompose that error into directional bias and spatial error variance. This breakdown helps isolate the mechanisms governing polynomial and spline performance.

The Mean Squared Error (MSE) of any estimator ẑ(x, y) over a continuous domain decomposes into squared bias and residual spatial error variance:

where the bias is defined as:

which measures the overall mean directional offset, and the spatial variance of residual errors is defined as:

where e_i represents the error at sensor i and e_bar is the mean error.

The results of the bias-variance decomposition are summarized in the table below:

Insights from the Decomposition:

1. Low Global Squared Bias (Bias² ≤ 0.0095) vs. Spatial Underfitting: Global squared bias remains minimal for low-degree polynomials despite spatial underfitting because domain-wide errors cancel out. A Degree 2 surface captures Geneva’s mean temperature (≈ 27.5°C), but localized overpredictions over Lac Léman (+1.8°C) offset underpredictions on mountain slopes (-1.4°C), yielding a near-zero mean bias (MBE = +0.097°C). Structural underfitting is instead reflected by high Spatial Error Variance (Var(e) = 0.3475), as rigid equations generate large, spatially varying residuals across distinct terrain zones.
2. Spatial Error Variance Dominates Total Error: The primary contributor to grid error is spatial error variance (Var(e)). The inability of rigid coordinate functions to adapt to localized microclimates (such as shoreline breezes and mountain lapse rates) causes residual errors to vary substantially across space, elevating spatial error variance despite a small global mean bias.
3. Degree 5 Tradeoff: Increasing the polynomial to degree 5 reduces squared bias to its minimum value in Part 1 (0.0015) but spatial variance increases to 0.3565 due to Runge boundary oscillations.
4. Regularization via Spline Tuning (s = 1.0): Piecewise B-splines reduce squared bias by nearly 50% relative to quadratic polynomials (from 0.0094 to 0.0049) through localized knot adaptation. Setting the smoothing factor to s = 1.0 acts as a spatial regularizer, maintaining low squared bias (0.0049) while achieving a 71.5% reduction in spatial error variance (from 0.3475 to 0.0988) and establishing a balanced bias-variance tradeoff.

Summary: Pros & Cons of Bivariate Surface Models

Drawing from the analyses of the previous sections, we are finally in a position to synthesize their key architectural trade-offs below:

Conclusion: The Limits of Coordinate-Based Equations

Global and piecewise polynomial formulations have fundamental limitations in spatial interpolation applications:

1. Global Polynomials (High Bias vs. Boundary Instability): Low-order polynomials are too spatially rigid to capture localized microclimates (such as shoreline cooling breezes or mountain lapse rates). In contrast, high-order polynomials develop spurious oscillations near grid boundaries resulting from Runge’s phenomenon.
2. Piecewise Splines (Regularization Sensitivity): B-Splines address global rigidity by utilizing localized knots. But they exhibit significant sensitivity to the smoothing factor s. Unregularized exact fits (s = 0) induce numerical ill-conditioning across unmonitored areas, whereas default auto-smoothing produces generalized trends rather than precise local predictions.

Ultimately, the problem is that these coordinate-based mathematical functions f(x, y) treat space as a passive coordinate grid rather than a physical continuum of connected neighborhoods.

To map complex, real-world geographies like Geneva, we must shift our paradigm. Instead of fitting global or piecewise equations to the entire grid, we should look at distance-based local neighborhood methods.

In Part 2: Nearest Neighbor, IDW & RBF Splines

To map complex, real-world geographies like Geneva, we must shift our paradigm. Instead of fitting global or piecewise equations to the entire grid, we should look at distance-based local neighborhood methods.

In the next installment, we evaluate Nearest Neighbor, Inverse Distance Weighting (IDW), and Radial Basis Function (RBF) splines.

👉 Continue with Part 2: Nearest Neighbor, IDW & RBF.

References & Further Reading

1. Chorley, R. J., & Haggett, P. (1965). Trend-surface mapping in geographical research. Transactions of the Institute of British Geographers, (37), 47–67. DOI: 10.2307/621689
2. Unwin, D. J. (1975). An Introduction to Trend Surface Analysis. Concepts and Techniques in Modern Geography (CATMOG №5), Geo Abstracts, University of East Anglia. ISBN: 0-902246-51-8
3. Dierckx, P. (1993). Curve and Surface Fitting with Splines. Oxford University Press. ISBN: 978-0198534402. DOI: 10.1093/oso/9780198534419.001.0001
4. Hastie, T., Tibshirani, R., & Friedman, J. (2009). The Elements of Statistical Learning: Data Mining, Inference, and Prediction (2nd ed.). Springer. ISBN: 978-0387848570. Free Online PDF
5. Burden, R. L., Faires, J. D., & Burden, A. M. (2016). Numerical Analysis (10th ed.). Cengage Learning. ISBN: 978-1305253667.
6. Cressie, N. (1993). Statistics for Spatial Data (Revised ed.). John Wiley & Sons. ISBN: 978-0471002550. Wiley Library
7. NumPy & SciPy Documentation: np.linalg.lstsq & scipy.interpolate.bisplrep — Python implementations for OLS polynomial regression and bivariate B-splines (bisplev)

Mastering Spatial Interpolation: Polynomials & Splines was originally published in Stellar Priors on Medium, where people are continuing the conversation by highlighting and responding to this story.