Closed-Form Exclusion for the Rayleigh-Benard Critical Constants

Two pure numbers fix the moment a heated fluid layer starts to move. Here they are to eighty certified digits, and here is the evidence that neither is algebraic of low degree.

A layer of fluid heated from below stays still until a single dimensionless group, the Rayleigh number, crosses a threshold. The threshold is a pure number. It does not depend on the fluid, the layer depth, gravity, or the temperature difference separately, only on the boundary conditions at the two plates. For the physically standard case of two rigid, perfectly conducting plates, that number is $1707.76\ldots$, and the horizontal wave number selected at onset is $3.1163\ldots$. Both have been in textbooks since Chandrasekhar (Chandrasekhar, 1961).

Three boundary-condition pairs are conventionally distinguished: free-free, rigid-rigid, and the asymmetric rigid-free case. Only the first has a closed form. Rayleigh solved it in 1916 and got $\mathrm{Ra}_c = 27\pi^4/4$ and $a_c = \pi/\sqrt{2}$ (Rayleigh, 1916). The other two have resisted, and the reason is structural rather than incidental: the free-free boundary conditions are satisfied by a single sine mode, which collapses the eigenvalue problem to an algebraic one, while rigid boundaries force a superposition of three exponential branches and a transcendental secular determinant.

This note does four things. It computes the rigid-rigid and rigid-free constants to eighty certified digits, extending the fifty-digit interval-analysis results of Jeng and Hassard (Jeng & Hassard, 1999) and Glomski and Johnson (Glomski & Johnson, 2012). It then asks a question that appears not to have been asked of these constants before: do they satisfy any polynomial with integer coefficients, or any linear relation over a basis of classical constants? Using integer-relation detection at eighty digits, the answer is no, within explicit and stated height bounds. It asks the same question of the pair jointly, since two constants can satisfy an exact relation even when neither is individually algebraic, and finds none. Finally it applies the method, validated against two known closed forms, to four further constants, from lattice statistics and radiative transport, with the same outcome.

I want to be exact about what that second result is. It is an exclusion, not an impossibility proof. Integer-relation search can rule a relation out up to a coefficient size; it cannot rule one out entirely. What follows is stated in those terms throughout.

The eigenvalue problem

Under the Boussinesq approximation, linearising about the conduction state and separating horizontal Fourier modes of wave number $a$ reduces marginal stability to a sixth-order ordinary differential equation for the vertical velocity amplitude $W(z)$ on $z \in [-1/2, 1/2]$:

$$ (D^2 - a^2)^3 W = -a^2 \mathrm{Ra} W, \qquad D = \frac{d}{dz}. $$

At each plate three conditions apply. All cases here take perfectly conducting plates, which gives $(D^2-a^2)^2 W = 0$. A rigid plate adds $W = 0$ and $DW = 0$; a free plate adds $W = 0$ and $D^2 W = 0$.

Substituting $W \sim e^{\lambda z}$ gives $(\lambda^2 - a^2)^3 = -a^2\mathrm{Ra}$. Writing $\tau = (a^2 \mathrm{Ra})^{1/3}$, the six roots are $\lambda = \pm\lambda_j$ with

$$ \lambda_j^2 = a^2 - \tau \omega_j, \qquad \omega_j \in \lbrace 1, e^{2\pi i/3}, e^{-2\pi i/3} \rbrace. $$

One branch is real in $\omega$ and two are a complex-conjugate pair, so the general solution mixes hyperbolic and trigonometric behaviour. This is the source of the transcendence.

Why free-free closes and the others do not

For two free plates, $W = \cos(\pi z)$ satisfies every boundary condition exactly. Then $D^2 W = -\pi^2 W$, the differential equation becomes the algebraic relation $\mathrm{Ra} = (\pi^2+a^2)^3/a^2$, and minimising over $a$ gives $a_c^2 = \pi^2/2$ and $\mathrm{Ra}_c = 27\pi^4/4$. The closed form exists because a single mode diagonalises the operator.

Rigid plates destroy this. No single exponential satisfies $W = DW = 0$ simultaneously, so all three $\lambda_j$ branches must be superposed and their amplitudes eliminated through a determinant. For the symmetric rigid-rigid case the even modes decouple and the condition reduces to a $3\times3$ determinant,

$$ \det M = 0, \qquad M_{1j} = \cosh(\lambda_j/2), \quad M_{2j} = \lambda_j \sinh(\lambda_j/2), \quad M_{3j} = \omega_j^2 \cosh(\lambda_j/2), $$

with $j = 1,2,3$ indexing the columns. The three rows are, in order, the conditions $W = 0$, $DW = 0$, and $(D^2-a^2)^2 W = 0$ at the plate.

The asymmetric rigid-free case has no midplane symmetry to exploit and requires the full $6\times6$ determinant over $\lbrace \cosh \lambda_j z, \sinh \lambda_j z \rbrace$. In both cases $\mathrm{Ra}_c$ is obtained by solving the determinant for $\mathrm{Ra}(a)$ and then minimising over $a$. This is the classical route of Pellew and Southwell (Pellew & Southwell, 1940) and Reid and Harris (Reid & Harris, 1958).

Certified values

Working in arbitrary precision, the determinant condition was solved for $\mathrm{Ra}(a)$ by Newton iteration and minimised over $a$ by a central-difference stationarity condition. Digit stability was established by repeating the whole computation at two working precisions and two step sizes: 120 digits with $h = 10^{-40}$, and 160 digits with $h = 10^{-55}$. The two runs agree to 81 significant digits in $a_c$ and to more than 120 in $\mathrm{Ra}_c$. The asymmetry is expected and is a useful internal check: $\mathrm{Ra}_c$ sits at a minimum, so an error $\delta a$ in the wave number perturbs it only at order $(\delta a)^2$.

For two rigid plates, to 80 digits:

Ra_c = 1707.
       7617771047 1835827375 5982297253 3944615499
       6330430342 7792799043 1141645326 1720675937

 a_c = 3.
       1163235548 2112203579 6173386585 7178246752
       2390095943 6186742173 1360101248 4416990067

For a rigid lower plate and a free upper surface, to 70 digits:

Ra_c = 1100.
       6496068876 7678462749 1888847271 1957506710
       5545241401 2234397846 5420212483

 a_c = 2.
       6823217576 9341424484 3892930110 7059823340
       0144144019 9325699178 0297015430

The rigid-free pair is quoted to 70 digits rather than 80. Every digit shown is certified: the wave number is independently confirmed to 74 significant digits by the cross-check in the next section, and the truncation to 70 leaves margin rather than reporting to the edge of what two methods agree on.

Validation

A high-precision claim is only as good as its independent checks, so there are four.

The free-free case was run through the identical pipeline as a control. It returns $657.5113644795164513459722456487595009357$ against $27\pi^4/4$, and $2.221441469079183123507940495030$ against $\pi/\sqrt{2}$, agreeing to every digit computed. A pipeline that reproduces the one exactly solvable case to 40 digits is not silently mis-stating boundary conditions.

The rigid-rigid case was recomputed by Chebyshev collocation (Trefethen, 2000), a completely different discretisation that never forms the characteristic roots $\lambda_j$. At 32 collocation points it returns 1707.76177710471835827375597, agreeing with the secular-determinant value to 26 significant digits, with the departure at exactly the level expected from spectral truncation at that resolution.

The rigid-free case was recomputed from the independent analytic secular equation of Glomski and Johnson (Glomski & Johnson, 2012), which is expressed in terms of $\phi\cot(\phi/2)$ and a separate hyperbolic-trigonometric ratio rather than a determinant. Solving their equation reproduces the values above to 74 significant digits in $a_c$ and 84 in $\mathrm{Ra}_c$.

Finally, the values agree with the published fifty-digit interval-analysis results (Jeng & Hassard, 1999; Glomski & Johnson, 2012), which are rigorous error-bounded enclosures rather than convergence estimates. One caveat is worth recording. The rigid-free Rayleigh number agrees with the printed value of Glomski and Johnson to 54 digits, but the wave number as I transcribed it from their Theorem 4.1 diverges from mine near the 45th place. Solving their own equation independently returns my digits to 74 places, so the difference is a transcription or typesetting matter in the printed digit string rather than a disagreement between the computations. I note it only so the reader who compares the two papers is not puzzled.

Exact-structure analysis

The constants are now known far past the point where integer-relation detection becomes informative. The PSLQ algorithm (Ferguson et al., 1999) takes a vector of real numbers and returns a small integer vector annihilating it, or fails; a failure at a given precision and coefficient bound is a genuine exclusion within those bounds (Bailey & Borwein, 2015).

The discipline matters more than the search. With $n$ basis elements and $D$ correct digits, PSLQ will return a spurious relation with coefficients of size roughly $10^{D/n}$ purely by counting, so a returned relation is meaningful only if its coefficients are far below that threshold, and an exclusion is meaningful only if the height bound is stated. A first pass at 40 digits duly produced apparent minimal polynomials of degree 7 to 9 for all four constants, with coefficients of order $10^4$ to $10^5$ against a spurious threshold of $10^{40/8} \approx 10^5$. Re-evaluating those candidate polynomials on the 80-digit values gives relative residuals of $2\times10^{-40}$ to $4\times10^{-42}$, precisely the noise floor of the 40-digit input that generated them, where a genuine relation would give $10^{-75}$. All four were artefacts. This is the failure mode that makes casual constant-hunting unreliable, and it is why the bounds below are quoted explicitly.

Repeating the search on the 80-digit values, with the coefficient bound for degree $d$ set to $10^{70/(d+1)}$:

No integer polynomial annihilates any of the four constants at these degrees and heights. A separate search tested each constant, and its cube, logarithm, and quotient by $\pi^4$, for a linear relation over the basis

$$ \lbrace 1, \pi, \pi^2, \pi^3, \pi^4, \pi^6, \sqrt{2}, \sqrt{3}, \ln 2, \ln 3, \zeta(3), G, \gamma \rbrace $$

with $G$ Catalan’s constant and $\gamma$ the Euler-Mascheroni constant, at coefficient bound $10^5$ and tolerance $10^{-70}$. Every combination returned nothing. The cube and the $\pi^4$ quotient were included because the free-free solution has the shape $27\pi^4/4$ and the characteristic equation enters through $\tau^3$, so those are the forms a surviving closed relative of the solvable case would most plausibly take.

Are the two constants jointly algebraic?

Asking whether each constant separately has a closed form is the obvious question, but it is not the only one. In the solvable case the two are bound together: eliminating nothing at all, $\mathrm{Ra}_c a_c^2 = (\pi^2 + a_c^2)^3$, and separately $\pi^6 = 8 a_c^6$. A pair can satisfy an exact algebraic relation even when neither member is individually algebraic, and if the rigid cases did so it would be a genuine closed-form result about the onset of convection, just not of the expected shape.

That relation fixes a natural weighting: $\mathrm{Ra}$ counts as four powers of $a$, and $\pi$ as one. Searching the span of monomials $\mathrm{Ra}^i a^j \pi^k$ of fixed total weight $W$ therefore covers the free-free relation at $W = 6$, and the search does recover it, returning $\pi^6 - 8a^6$ with height 8 against a spurious threshold of $10^{7.2}$.

Run on the rigid cases at $W = 6, 7, 8, 9, 10$, it returns nothing. Every vector has height at or above its threshold:

Neither pair is jointly algebraic over $\pi$ at these weights and heights. The same holds for the natural derived quantities at criticality: $\tau_c = (a_c^2\mathrm{Ra}_c)^{1/3} = 25.5017974035\ldots$ and the real characteristic root $q_0 = \sqrt{\tau_c - a_c^2} = 3.9737041793\ldots$ are neither algebraic of degree at most four nor expressible over $\lbrace 1, \pi, \pi^2, \sqrt{2}, \sqrt{3} \rbrace$.

The same test on five other universal constants

An exclusion is only worth reporting if the method that produced it can be shown to find closed forms when they exist. The free-free control does that for the Rayleigh-Benard pipeline. The integer-relation step deserves the same treatment, and it also generalises, so it is worth applying to constants of a different type.

One caveat first, and it is the reason the search above nearly missed its own target. A multiplicative closed form built from gamma values is invisible to an integer-relation search on a linear basis; it appears only on a logarithmic one. Watson’s constant, the return integral for the simple cubic lattice, is exactly of this kind (Watson, 1939):

$$ W_3 = \int_0^\infty e^{-t} I_0(t/3)^3 dt = \frac{\sqrt{6}}{32\pi^3} \Gamma(1/24)\Gamma(5/24)\Gamma(7/24)\Gamma(11/24). $$

Computing $W_3$ to 62 digits by quadrature with an analytic tail correction, and handing the logarithm to PSLQ against $\lbrace \ln 2, \ln 3, \ln \pi, \ln\Gamma(k/24) \rbrace$, returns the coefficient vector $(2, 9, -1, 6, -2, -2, -2, -2)$ with residual $1.6\times10^{-61}$. That is Watson’s formula, recovered from the number alone, with coefficients of order ten against a spurious threshold of $10^{7.8}$. The method demonstrably works.

Applied to the constants where no closed form is known, it returns nothing:

• $W_4$ and $W_5$, the four- and five-dimensional analogues, tested on logarithmic gamma bases with denominators $N \in \lbrace 3,4,6,8,12,16,20,24,30,36,48,60 \rbrace$. Every returned vector had coefficients at the spurious threshold. This is consistent with the four-dimensional case belonging to the Calabi-Yau and L-value setting rather than admitting gamma products (Guttmann, 2010).
• The simple cubic spanning-tree entropy, $1.6733893029701967322834306216555980752577$, against $\lbrace 1, \ln 2, \ln 3, \ln 5, G/\pi, \zeta(3)/\pi^2, \pi, \sqrt{2}, \sqrt{3}, 1/\pi \rbrace$. No relation below the spurious threshold of $10^{4.7}$. Here the control is the two-dimensional case, whose entropy is $4G/\pi$; the same quadrature returns that value with zero error at full working precision.

One more constant is worth adding, from a third field and of a more favourable type. Hopf’s constant $q_\infty$ is the extrapolation length of the Schwarzschild-Milne problem, the universal distance beyond a half-space boundary at which the asymptotic diffusion solution extrapolates to zero. It governs neutron diffusion and stellar atmospheres, and unlike the stability constants it is not a root of a transcendental determinant but an elementary definite integral (Finch, 2008):

$$ q_\infty = \frac{1}{\pi}\int_0^{\pi/2} \left( \frac{3}{\sin^2\theta} - \frac{1}{1 - \theta\cot\theta} \right) d\theta . $$

The two terms diverge separately at the origin and their difference tends to $6/5$, so the singularity is removable. Evaluating it by quadrature with the series taken analytically near zero, and independently by Hopf’s Bernoulli-coefficient series, gives agreement to 45 digits:

$$ q_\infty = 0.710446089598763072732524141699153671993201334 . $$

That is the most promising shape in this collection, since an integrand built from $\sin$, $\cot$ and rational functions is exactly the sort that sometimes evaluates. It does not. The constant is not algebraic of degree at most eight, and no relation survives over $\lbrace 1, \pi, \pi^2, 1/\pi, 1/\pi^2, 1/\pi^3 \rbrace$, over $\lbrace 1, 1/\pi^2, \ln 2, \ln 3, \zeta(3)/\pi^3, \zeta(3)/\pi^2, G/\pi \rbrace$, or over a wider basis including $\gamma$ and $\sqrt{3}$. The same holds for the bare integral $\pi q_\infty - 6/\pi$.

One further case is worth recording with an explicit caveat about its weight. The Thomas-Fermi constant $B = -y’(0) = 1.5880710226113753127186845\ldots$, the initial slope of the solution to $y’’ = y^{3/2}x^{-1/2}$ with $y(0)=1$ and $y(\infty)=0$, has been without a closed form since 1927, and it is a universal constant in a strong sense: Lieb and Simon proved that it fixes the leading coefficient of the ground-state energy of an atom as $Z \to \infty$, through $E \sim -(3/7)BZ^{7/3}$. It is the most physically consequential constant in this collection.

It is also the one that resists a naive approach. Shooting cannot reach high precision, because the mode that must be tuned away grows only as $x^{0.772}$, so resolving even 40 digits would require integrating to $x \sim 10^{52}$. Direct shooting confirms the published value to nine significant digits and no further.

The way through is to integrate the other way. Inward from large $x$ the offending mode decays rather than grows, so the integration is stable, and the starting condition is supplied by the Coulson-March expansion $y = (144/x^3)\sum_k c_k(Ax^{-\lambda})^k$ with $\lambda = (\sqrt{73}-7)/2$. Two facts make this practical. The amplitude $A$ has two sign branches, and because the scaling symmetry $y \mapsto \alpha^3 y(\alpha x)$ preserves $\operatorname{sign}(A)$, they are genuinely distinct families rather than a normalisation; the physical solution needs $A < 0$, since $y(0)$ is finite while $144/x^3$ is not. And that same scaling means any single inward integration suffices: whatever $y(0)$ comes out, the physical constant is recovered as $B = y(0)^{-4/3},|y’(0)|$. No shooting is required at all.

Stepping inward with a Taylor method and matching to the series in $\sqrt{x}$ near the origin, two runs at deliberately different settings, differing in Taylor order, step ratio, starting radius, matching point and matching-series length, agree to 114 significant digits. Conservatively,

B = 1.
    5880710226 1137531271 8684509423 9501094527 4662167482
    5616765677 4181665519 6115430926 2332033970 1384286652
    6981402442

Against this, the search at 105 digits excludes far more than before: no integer polynomial of degree at most 16 with height below $10^{105/2(d+1)}$, which is a bound of $10^{21}$ at degree four and $10^{11.7}$ at degree eight, and no relation over $\lbrace 1,\pi,\pi^2,\pi^3 \rbrace$, over $\lbrace 1,\pi,\pi^2,\ln 2,\ln 3 \rbrace$, over $\lbrace 1, 2^{1/3}, 3^{1/3}, 9^{1/3}, 3^{1/7} \rbrace$, or over $\lbrace 1,\pi,\zeta(3),G,\gamma \rbrace$. The Thomas-Fermi constant now carries an exclusion comparable in strength to the hydrodynamic ones rather than the token bound of a first pass.

So four constants from hydrodynamic stability, three from lattice statistics, and one from radiative transport all decline to have a closed form, under a procedure that recovers Watson’s gamma product and Rayleigh’s $27\pi^4/4$ when pointed at them.

What this establishes

The positive content is a set of reference values at eighty digits, cross-validated by four independent routes, for constants whose best published precision was fifty digits and whose textbook precision is seven.

The negative content is narrower than it may look, and I would rather understate it. These searches show that if $\mathrm{Ra}_c$ for rigid boundaries is algebraic, its minimal polynomial has degree above 12 or coefficients larger than the stated bounds; and that it admits no simple linear expression over the classical basis tested. That is real evidence against a Rayleigh-style closed form of the kind the free-free case enjoys, and it is consistent with the structural argument in the third section. It is not a proof of transcendence, and nothing here constitutes one. Proving that a root of a transcendental determinant is not algebraic is a different and much harder problem, of the kind that remains open even for far more studied constants.

One qualification belongs here, because it cuts against the pessimistic reading. The absence of a closed form is a property of these particular points in parameter space, not of the underlying problem. Replacing the rigid condition by Navier slip embeds the no-slip case in a one-parameter family whose other end is Rayleigh’s solvable case, and near that end the problem is not opaque at all: the leading corrections come out as $12\pi^2$ and $8/9$ exactly, derived in a companion note. The rigid point is hard; the family containing it is not uniformly so.

The broader point is methodological. Constants like these sit in a gap. The experimental-mathematics community has combed the lattice sums and Watson integrals for closed forms for decades (Watson, 1939), but constants that arise in fluid mechanics and are quoted numerically in engineering references have largely escaped that treatment. The Rayleigh-Benard constants turn out to have no low-height algebraic structure, which is a mildly disappointing but genuinely informative answer, and one worth recording so the question is not repeatedly reopened. Other classical stability thresholds, the critical Reynolds number for plane Poiseuille flow among them, are natural candidates for the same treatment.

Reproducibility

Every number above follows from the equations as stated. The determinant condition, the two-precision convergence protocol, the Chebyshev cross-check, and the PSLQ bounds are each specified completely enough to be re-derived independently, and the free-free control gives any reimplementation an immediate correctness test: if it does not return $27\pi^4/4$, it is wrong before it reaches the cases that matter.

References

Letchford, B. (2026, July 27). Closed-Form Exclusion for the Rayleigh-Benard Critical Constants. benletchford.com. https://benletchford.com/writing/closed-form-exclusion-rayleigh-benard/

@misc{letchford2026closed,
  author = {Letchford, Ben},
  title = {Closed-Form Exclusion for the Rayleigh-Benard Critical Constants},
  year = {2026},
  month = jul,
  url = {https://benletchford.com/writing/closed-form-exclusion-rayleigh-benard/}
}