Robbins’s constant

Solution. The answer is Robbins’s constant: writing $L_1 = \ln(1 + \sqrt 2)$ and $L_2 = \ln(2 + \sqrt 3)$, $$\boxed{\mathbb{E}[|P_1 - P_2|] = \frac{4 + 17\sqrt{2} - 6\sqrt{3} - 7\pi + 21 L_1 + 42 L_2}{105} \approx 0.6617.}$$ This was first computed in closed form by David P. Robbins in 1978. The derivation uses $$\mathbb{E}[|P_1 - P_2|] = \int_{[0,1]^6} \!\! \sqrt{(x_1\!-\!x_2)^2 + (y_1\!-\!y_2)^2 + (z_1\!-\!z_2)^2} \, dV,$$ which reduces (after change of variable to $(u, v, w) = (|x_1-x_2|, \ldots)$ with triangular marginals) to $$\iiint_{[0,1]^3} \!\! \sqrt{u^2+v^2+w^2} \cdot 8(1-u)(1-v)(1-w) \, du \, dv \, dw,$$ a 3D integral. Its evaluation involves careful integration by parts and yields the explicit combination of $\sqrt 2$, $\sqrt 3$, $\pi$, and two logarithms above.

Compare:

• 1D (unit segment): $\mathbb{E}[|X_1 - X_2|] = \tfrac13$ (clean).

• 2D (unit square, interior): $\mathbb{E}[|P_1 - P_2|] = \tfrac{2 + \sqrt 2 + 5 \ln(1+\sqrt 2)}{15} \approx 0.5214$. (The distinct boundary-of-square constant, for two points on the perimeter, is $\tfrac{3 + \sqrt 2 + 5\ln(1+\sqrt2)}{12} \approx 0.7351$; do not conflate the two.)

• 3D (unit cube): Robbins’s constant above.

Each successive dimension introduces new transcendental terms, and no closed form is known for general $n$-dimensional cubes. ■

import numpy as np
from math import sqrt, log, pi
p = np.random.rand(3, 2, 10**6)
d = np.linalg.norm(p[:,0] - p[:,1], axis=0)
exact = (4 + 17*sqrt(2) - 6*sqrt(3) + 21*log(1+sqrt(2))
         + 42*log(2+sqrt(3)) - 7*pi)/105
print(f"sim: {d.mean():.5f}   exact: {exact:.5f}")
# sim: 0.66156   exact: 0.66171