Exact Slip Corrections to the Onset of Convection

Convection between slippery walls has only ever been solved numerically. Perturbing about the one case that closes exactly gives the correction as a clean pure number, and for insulating plates the whole family closes in closed form.

The onset of convection in a fluid layer heated from below is fixed by a pure number, the critical Rayleigh number. Which number depends on what the fluid does at the two plates. For stress-free boundaries Rayleigh solved the problem in 1916 and got $\mathrm{Ra}_c = 27\pi^4/4$ with squared wavenumber $a_c^2 = \pi^2/2$ (Rayleigh, 1916). For rigid no-slip boundaries there is no closed form; the value $1707.7617771\ldots$ has been computed numerically ever since Jeffreys and Chandrasekhar (Chandrasekhar, 1961).

Real surfaces are neither. A fluid at a solid wall slips a little, by an amount set by a slip length $\ell$, in Navier’s original boundary condition from 1823 (Navier, 1823). On ordinary surfaces $\ell$ is nanometres and the no-slip idealisation is excellent, but on superhydrophobic and liquid-infused surfaces it reaches tens or hundreds of micrometres (Lauga et al., 2007; Rothstein, 2010), which is not small compared with a thin convection cell. Convection between Navier-slip boundaries has accordingly become a subject in its own right, mostly through rigorous bounds on heat transport (Drivas et al., 2022) and through numerics, which report that the critical Rayleigh number falls and the critical wavenumber rises as slip increases.

I have not found that trend given a closed form in the accessible literature, though see the caveat on priority below. The stress-free case is an exactly solvable point of the slip family, sitting at $\ell = \infty$, and one can perturb about it in $\varepsilon = 1/\ell$. The result is that the leading corrections are exact rational multiples of powers of $\pi$:

$$ \mathrm{Ra}_c(\ell) = \frac{27\pi^4}{4} + \frac{12\pi^2}{\ell} + O(\ell^{-2}), \qquad a_c^2(\ell) = \frac{\pi^2}{2} + \frac{8}{9\ell} + O(\ell^{-2}). $$

The constant is $12\pi^2 = 118.435252813072303426\ldots$. Both results are confirmed against independent high-precision numerics to eighteen digits or better, and the derivation is checked in a stronger form than the constants themselves: the full wavenumber-dependent correction function is reproduced to nineteen digits at six separate wavenumbers.

There is also a case where the perturbation is unnecessary because the whole family closes. If the plates are insulating rather than conducting, the marginal mode is long-wavelength, the sixth-order problem degenerates to a polynomial one, and the critical Rayleigh number is exact at every slip length:

$$ \mathrm{Ra}_c(\ell) = \frac{720(2\ell+1)}{12\ell+1} . $$

Its endpoints are the two classical insulating constants, $120$ as $\ell \to \infty$ and $720$ as $\ell \to 0$, and it agrees with a full numerical solution to 33 digits at every slip length tested. The same velocity condition that leaves the conducting problem transcendental leaves the insulating one rational.

The slip family

Lengths are measured in units of the layer depth, so $z \in [-1/2, 1/2]$ and $\ell$ is a dimensionless slip length. Marginal stability reduces, as usual, to a sixth-order problem for the vertical velocity amplitude $W(z)$ at horizontal wavenumber $a$:

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

Both plates are perfectly conducting, giving $W = 0$ and $(D^2-a^2)^2 W = 0$ at $z = \pm 1/2$. The third condition carries the slip. Navier’s condition sets the tangential velocity proportional to its own normal derivative, $u = \ell \partial_n u$ with $n$ the inward normal. For a mode $w = W(z)e^{iax}$ continuity gives $u = (i/a) DW$, so $\partial_z u = (i/a) D^2 W$, and at the upper plate, where the inward normal points along $-z$,

$$ D^2 W + \varepsilon DW = 0 \quad \text{at } z = +\tfrac12, \qquad D^2 W - \varepsilon DW = 0 \quad \text{at } z = -\tfrac12, $$

with $\varepsilon = 1/\ell$. Setting $\varepsilon = 0$ recovers the stress-free conditions $D^2W = 0$; letting $\varepsilon \to \infty$ recovers no-slip, $DW = 0$. The family therefore interpolates between the two classical cases, and the solvable one sits at $\varepsilon = 0$, which is where a perturbation expansion can start.

Perturbing about the solvable point

At $\varepsilon = 0$ the even mode is exactly $W_0 = \cos(\pi z)$. It satisfies every boundary condition identically, because $\cos(\pi/2) = 0$ kills $W_0$, $D^2W_0$ and $D^4W_0$ at the plates simultaneously. Substituting gives the algebraic marginal curve

$$ \mathrm{Ra}_0(a) = \frac{(\pi^2+a^2)^3}{a^2}, $$

whose minimum over $a$ is Rayleigh’s result. The one boundary quantity that does not vanish is the slope, $DW_0(\tfrac12) = -\pi$, and this is exactly what the slip condition couples to. The perturbation therefore enters at first order.

Write $W = W_0 + \varepsilon W_1 + \ldots$ and $\mathrm{Ra} = \mathrm{Ra}_0 + \varepsilon \mathrm{Ra}_1 + \ldots$, and let $L = (D^2-a^2)^3 + a^2\mathrm{Ra}_0$, so that $LW_0 = 0$. At first order,

$$ L W_1 = -a^2 \mathrm{Ra}_1 W_0 , $$

with $W_1 = 0$ and $(D^2-a^2)^2W_1 = 0$ at the plates as before, but now with an inhomogeneous third condition inherited from the slip term:

$$ D^2 W_1 = -DW_0 = \pi \quad \text{at } z = \tfrac12 . $$

The operator $L$ contains only even derivatives and is formally self-adjoint, so $\int (W_0 L W_1 - W_1 L W_0) dz$ collapses to a boundary concomitant. For $(D^2-a^2)^3$ that concomitant is

$$ J = \left[ uv^{(5)} - u’v^{(4)} + u’‘v’‘’ - u’’‘v’’ + u^{(4)}v’ - u^{(5)}v \right] - 3a^2\left[ uv’‘’ - u’v’’ + u’‘v’ - u’’‘v \right] + 3a^4\left[ uv’ - u’v \right]. $$

Evaluating at $z = \tfrac12$ with $u = W_0$ and $v = W_1$, most terms die. From $W_0$ we have $W_0 = 0$, $W_0’’ = 0$, $W_0^{(4)} = 0$, together with $W_0’ = -\pi$ and $W_0’‘’ = \pi^3$. From $W_1$ we have $W_1 = 0$, $W_1’’ = \pi$, and, from the conducting condition, $W_1^{(4)} = 2a^2W_1’’ = 2a^2\pi$. Only three products survive, coming from $-u’v^{(4)}$, from $-u’’‘v’‘$, and from the $-3a^2$ bracket via $-u’v’’$:

$$ J(\tfrac12) = 2a^2\pi^2 - \pi^4 - 3a^2\pi^2 = -\pi^2(\pi^2+a^2). $$

Both $W_0$ and $W_1$ are even, and every term in $J$ pairs derivatives of orders summing to five, so $J$ is odd and the two plates contribute equally: the total is $2J(\tfrac12)$. Since $\int_{-1/2}^{1/2}\cos^2(\pi z) dz = \tfrac12$, the solvability condition reads $-\tfrac12 a^2\mathrm{Ra}_1 = -2\pi^2(\pi^2+a^2)$, giving the first-order correction as an explicit function of wavenumber:

$$ \mathrm{Ra}_1(a) = \frac{4\pi^2(\pi^2+a^2)}{a^2}. $$

The two constants

To first order the minimising wavenumber may be held fixed, because $\mathrm{Ra}_0$ is stationary there and the shift contributes only at second order. Evaluating at $a_0^2 = \pi^2/2$,

$$ \mathrm{Ra}_1(a_0) = \frac{4\pi^2\left(\pi^2 + \tfrac{\pi^2}{2}\right)}{\tfrac{\pi^2}{2}} = 12\pi^2 , $$

and the $\pi$-dependence of the ratio cancels exactly, leaving a rational multiple of $\pi^2$.

The wavenumber shift follows from the same two functions. Writing $s = a^2$, stationarity of $\mathrm{Ra}_0 + \varepsilon \mathrm{Ra}_1$ gives $s - s_0 = -\varepsilon \mathrm{Ra}_1’(s_0)/\mathrm{Ra}_0’’(s_0)$. Both derivatives are elementary at $s_0 = \pi^2/2$:

$$ \mathrm{Ra}_0’‘(s_0) = 18, \qquad \mathrm{Ra}_1’(s_0) = -\frac{4\pi^4}{s_0^2} = -16 , $$

so $a_c^2 = \pi^2/2 + (8/9)\varepsilon$. Both are pure numbers with no $\pi$ left in them, which is a small surprise: the $18$ and the $-16$ each carry $\pi^4$ that cancels in the ratio.

The same two derivatives fix part of the second order. The envelope contribution from the wavenumber shift is $-\mathrm{Ra}_1’^2/(2\mathrm{Ra}_0’’) = -256/36$, so

$$ c_2 = \mathrm{Ra}_2(a_0) - \frac{64}{9}, $$

where $\mathrm{Ra}_2(a_0)$ requires second-order perturbation theory that I have not carried out.

There is a structural reason to expect second order to be qualitatively worse, and it explains why first order came out so clean. At the unperturbed solution the cube-root parameter degenerates exactly:

$$ \tau = (a^2\mathrm{Ra}_0)^{1/3} = \left( (\pi^2+a^2)^3 \right)^{1/3} = \pi^2 + a^2 , $$

so the first characteristic root is $\lambda_1^2 = a^2 - \tau = -\pi^2$, which is precisely the $\cos(\pi z)$ mode. That is why the unperturbed problem closes: one of the three branches lands exactly on a half-integer number of wavelengths across the layer. The other two branches do not. They sit at $\lambda^2 = a^2 + \tfrac12(\pi^2+a^2) \mp \tfrac{i\sqrt3}{2}(\pi^2+a^2)$, which are irrational and complex, and $W_1$ is built from hyperbolic and circular functions of them. Any second-order coefficient must therefore involve transcendental values at those arguments, not just powers of $\pi$.

That expectation is borne out. Fitting $c_2$ to 25 digits and testing it against $\lbrace 1, \pi^2, \pi^4, \pi^6 \rbrace$, $\lbrace 1, \pi^2, 1/\pi^2 \rbrace$ and $\lbrace 1, \pi^2, \pi^4, 1/\pi^2 \rbrace$ returns only vectors at the spurious height threshold. So the clean rational-times-$\pi$ structure is a first-order phenomenon, not a property of the whole expansion.

This also closes off a tempting idea. Since $\varepsilon$ runs from $0$ at the solvable end to $\infty$ at no-slip, an expansion with closed-form coefficients at every order would express $1707.7617771\ldots$ as the endpoint of an exactly known series. The second-order obstruction above, together with the fact that the no-slip point sits at infinite $\varepsilon$ rather than inside any disc of convergence, makes that route unavailable. Numerically $c_2 = -17.69212625334931\ldots$, hence $\mathrm{Ra}_2(a_0) = -10.58101514223820\ldots$, which I could not identify in closed form against bases built from $1, \pi, \pi^2, \pi^4$. The $-64/9$ is exact; the rest of $c_2$ is not claimed.

Verification

The derivation was checked against an independent numerical solution of the full slip problem, using the even-mode $3\times3$ secular determinant with the slip row $\lambda_j^2\cosh(\lambda_j/2) + \varepsilon\lambda_j\sinh(\lambda_j/2)$, solved in 50-digit arithmetic.

Three checks matter, in increasing order of strength.

At $\varepsilon = 0$ the solver returns $657.511364479516451345972245649$, which is $27\pi^4/4$ to every digit computed. As $\varepsilon \to \infty$ it flows to the no-slip values: at $\varepsilon = 10^5$ it gives $\mathrm{Ra}_c = 1707.6525\ldots$ and $a_c = 3.11628\ldots$, against the known limits $1707.76177710\ldots$ and $3.11632355\ldots$. The family is the right family.

Extracting the slope at the solvable end, a Vandermonde fit over $\varepsilon = 10^{-6}$ to $10^{-10}$ gives $c_1 = 118.435252813072303426012841788$, against $12\pi^2 = 118.435252813072303426013891999$: agreement to 23 significant digits. The corresponding fit for the wavenumber gives $0.888888888888888888$, against $8/9$, to 18 digits.

The strongest check is not of either constant but of the function behind them. Fixing $a$ and not minimising, the predicted correction $\mathrm{Ra}_1(a) = 4\pi^2(\pi^2+a^2)/a^2$ can be compared wavenumber by wavenumber:

Agreement to nineteen digits at six independent wavenumbers is not a coincidence at a stationary point; it confirms the concomitant calculation as an identity in $a$.

The insulating family closes exactly

Everything above is perturbative, because the fixed-temperature problem is transcendental away from $\varepsilon = 0$. Change the thermal boundary condition and that stops being true.

If the plates are insulating rather than perfectly conducting, so that the heat flux rather than the temperature is held fixed, the marginal mode is long-wavelength: $a_c \to 0$, and the critical Rayleigh number is the $a \to 0$ limit of the marginal curve. This is the standard fixed-flux situation, and it is what makes the problem tractable, because at $a = 0$ the eigenvalue problem degenerates to a polynomial one.

Write $W = a^2\hat{W}$ and expand in $a^2$. At leading order the temperature perturbation is constant, $\Theta_0 = 1$, and the velocity satisfies

$$ D^4 \hat{W}_0 = \mathrm{Ra}, $$

subject to the velocity conditions at the plates. The next order carries the eigenvalue: $D^2\Theta_1 = \Theta_0 - \hat{W}_0$, and insulating plates force $D\Theta_1 = 0$ at both ends, so integrating across the layer gives the solvability condition

$$ \int_{-1/2}^{1/2} \hat{W}_0 , dz = 1 . $$

That is now pure algebra. The even solution is $\hat{W}_0 = \mathrm{Ra}z^4/24 + c_2z^2 + c_0$. The slip condition $D^2\hat{W}_0 + \varepsilon D\hat{W}_0 = 0$ at $z = 1/2$ gives

$$ c_2 = -\frac{\mathrm{Ra}(6+\varepsilon)}{48(2+\varepsilon)}, $$

and $\hat{W}_0(1/2) = 0$ fixes $c_0$. Evaluating the integral leaves $-\mathrm{Ra}/480 - c_2/6 = 1$, and substituting $c_2$ collapses everything to

$$ \mathrm{Ra}_c(\varepsilon) = \frac{720(2+\varepsilon)}{12+\varepsilon}, \qquad\text{equivalently}\qquad \mathrm{Ra}_c(\ell) = \frac{720(2\ell+1)}{12\ell+1}. $$

This is not an asymptotic statement. It is exact for every slip length, and it contains both classical insulating constants as its endpoints: $\ell \to \infty$ gives $\mathrm{Ra}_c = 120$, the free-insulating value, and $\ell \to 0$ gives $720$, the rigid-insulating one. The whole one-parameter family of critical Rayleigh numbers between them is a rational function of the slip length.

Checking it against a full numerical solution of the sixth-order problem, with the insulating condition imposed as $D(D^2-a^2)^2W = 0$ and the limit $a \to 0$ taken by fitting the marginal curve in powers of $a^2$:

Thirty-three digits at five independent slip lengths, with rational values reproduced exactly.

The long-wavelength framework here is not mine. It is the expansion of Chapman and Proctor (Chapman & Proctor, 1980), following Childress and Spiegel, for convection between poorly conducting boundaries. Their equation (3.3) is the statement that the leading velocity satisfies $D^4P(z) = -1$ with $P$ fixed by the velocity conditions at the plates, which is the equation solved above up to normalisation, and their $O(\varepsilon^2)$ balance supplies the same solvability condition. What is added here is only the polynomial itself. They state explicitly, in their equation (2.7), that they treat three cases: two stress-free boundaries, two rigid boundaries, and rigid-below with stress-free-above, remarking that the velocity conditions “arise only in the details of certain polynomials” relegated to their appendix. Navier slip is not among them, and a slip length does not appear in their analysis.

That is a more useful statement about priority than a failed search. The method is forty-five years old and completely standard; the endpoints $120$ and $720$ are older still. Carrying the same polynomial through with a slip parameter left in is the step I have not found taken, and it is a small step.

The contrast with the conducting case is the interesting part. The same velocity boundary condition, the same operator, the same one-parameter family; changing only the thermal condition at the plates turns a problem whose endpoint is the stubbornly transcendental $1707.7617771\ldots$ into one whose entire family is a ratio of two linear polynomials. What makes the difference is that insulating plates push the critical wavenumber to zero, and at zero wavenumber the sixth-order problem degenerates to $D^4\hat{W} = \mathrm{Ra}$, which polynomials solve.

On priority, and on how hard this was

I should be candid about what kind of result this is, because it would be easy to oversell.

Both derivations are elementary. First-order perturbation theory about a solvable eigenvalue problem, and a long-wavelength expansion closed by a solvability condition, are standard methods, and any reader who wanted $12\pi^2$ or the rational formula could obtain either in an afternoon once the question is posed. Nothing here resisted derivation. The constant was not previously out of reach; on the evidence I have, it simply had not been written down in this form, and that is a much weaker claim.

It is also a claim I cannot substantiate, and it is worth saying precisely what was checked. The recent preprint literature on Navier-slip Rayleigh-Benard is concentrated on rigorous bounds for the Nusselt number and on scaling laws interpolating between the free-slip $\mathrm{Ra}^{5/12}$ and no-slip $\mathrm{Ra}^{1/3}$ regimes (Drivas et al., 2022); searches there return no explicit large-slip expansion of the critical Rayleigh number, and no coefficient of this kind. That is a genuine negative result for the modern literature, and it is also the part of the literature least likely to contain it.

The risk lies earlier, and I could not check it. Sparrow, Goldstein and Jonsson treated mixed and imperfectly conducting boundary conditions for this problem in 1965, and the more recent slip literature includes work I could not access behind paywalls. For the insulating result the position is better, because the relevant paper was accessible: Chapman and Proctor set up exactly this expansion and restrict it to three velocity cases, none of them Navier slip (Chapman & Proctor, 1980), and the endpoints trace back to Hurle, Jakeman and Pike (Hurle et al., 1967). That is a checked statement rather than a failed search, though it still does not exclude the intervening literature. A perturbative coefficient like this is exactly the sort of thing that appears as an unremarked intermediate step in a longer paper, in a thesis, or in a textbook exercise. Anyone intending to cite this should check that first. If it is already known, the value here is a clean derivation and an unusually stringent numerical verification, not priority.

What I am confident of is correctness rather than novelty. The identity $\mathrm{Ra}_1(a) = 4\pi^2(\pi^2+a^2)/a^2$ holds to nineteen digits at six wavenumbers, the constants follow from it by elementary algebra, and the family reproduces both classical limits.

What it is good for, and what it is not

The practical statement is that a slippery cell convects sooner, by a definite amount:

$$ \mathrm{Ra}_c \simeq 657.5114 + 118.4353 \cdot \frac{d}{\ell}, $$

with $d$ the layer depth. A cell with a slip length ten times its depth sits about 1.8 per cent above the free-slip threshold rather than the 160 per cent that no-slip walls would impose. This gives the numerically observed trend a closed form and a coefficient, and it provides an exact benchmark for codes that treat slip boundaries, in the same way that $27\pi^4/4$ benchmarks the stress-free case.

The honest limits are worth stating plainly. This is an asymptotic result about the large-slip end, $\ell \gg d$. It says nothing useful about ordinary no-slip walls, which sit at the opposite, non-perturbative end of the same family, and it is the second order that would set the width of its window of validity, which I have only pinned down in part. The expansion is in the boundary condition, not in the amplitude, so it is a statement about linear onset alone and carries no information about heat transport above threshold. And while the free-free control and the no-slip limit together make it very unlikely that the underlying problem has been set up wrongly, the derivation is a perturbation argument checked numerically rather than a theorem with error bounds.

What it does establish is that the slip family, unlike the no-slip point it flows to, is not opaque to exact methods near its solvable end. The companion note on closed-form exclusion for the classical onset constants argues that $1707.76\ldots$ itself is not algebraic of low degree and admits no simple closed form. The contrast is the point: the difficulty there is genuine and structural, but it is a property of that particular point in the family, not of the problem as a whole.

References

Letchford, B. (2026, July 27). Exact Slip Corrections to the Onset of Convection. benletchford.com. https://benletchford.com/writing/exact-slip-correction-convection-onset/

@misc{letchford2026exact,
  author = {Letchford, Ben},
  title = {Exact Slip Corrections to the Onset of Convection},
  year = {2026},
  month = jul,
  url = {https://benletchford.com/writing/exact-slip-correction-convection-onset/}
}