Solutions for 19 October 2026

Call the integral \(I\). It is the average of \(\ln|P-Q|^{2}\) over two independent uniform points \(P=(x,y)\), \(Q=(u,v)\) of the unit square. The logarithmic singularity on the diagonal \(P=Q\) is integrable (in two dimensions \(\int r\ln r\,dr\) converges), so \(I\) is finite. Attacking the four integrals one at a time is hopeless: the first two already produce a tangle of arctangents and logarithms in the remaining variables. The way in is to notice that the integrand only depends on the differences.

Step 1: collapse to a double integral. Put \(s=x-u\) and \(t=y-v\). If \(x,u\) are independent and uniform on \([0,1]\), then \(s\) has the tent density \(1-|s|\) on \([-1,1]\) (the length of \(\{u\in[0,1]:u+s\in[0,1]\}\)); likewise for \(t\), and \(s,t\) are independent because they are built from disjoint pairs of variables. Hence

using evenness in \(s\) and in \(t\).

Step 2: polar coordinates on a triangle. The last integrand is symmetric in \(s\) and \(t\), so \(I\) is \(8\) times the integral over the triangle \(0\le t\le s\le 1\). There put \(s=r\cos\theta\), \(t=r\sin\theta\) with \(0\le\theta\le\pi/4\) and \(0\le r\le\sec\theta\) (the side \(s=1\)). Since \(\ln(s^{2}+t^{2})=2\ln r\),

Step 3: the radial integral. For \(k\ge0\) and \(R>0\), \(\int_{0}^{R}r^{k}\ln r\,dr=R^{k+1}\Bigl(\dfrac{\ln R}{k+1}-\dfrac{1}{(k+1)^{2}}\Bigr)\). Take \(R=\sec\theta\) and write \(\tau=\tan\theta\); then \((\cos\theta+\sin\theta)R^{3}=R^{2}(1+\tau)\) and \(\cos\theta\sin\theta\,R^{4}=\tau R^{2}\). Collecting terms, the inner integral equals

because \(\frac12-\frac{1+\tau}{3}+\frac{\tau}{4}=\frac16-\frac{\tau}{12}\) and \(-\frac14+\frac{1+\tau}{9}-\frac{\tau}{16}=-\frac{5}{36}+\frac{7\tau}{144}\).

Step 4: the angular integral. Now \(w=\tan\theta\) gives \(dw=\sec^{2}\theta\,d\theta\) and \(\ln\sec\theta=\frac12\ln(1+w^{2})\), so

Integration by parts gives \(\int_{0}^{1}\ln(1+w^{2})\,dw=\bigl[w\ln(1+w^{2})-2w+2\arctan w\bigr]_{0}^{1}=\ln2-2+\frac{\pi}{2}\); this arctangent is where \(\pi\) enters. The substitution \(q=1+w^{2}\) gives \(\int_{0}^{1}w\ln(1+w^{2})\,dw=\frac12\int_{1}^{2}\ln q\,dq=\ln2-\frac12\). Therefore the bracket integrates to

and

Equivalently, the mean of \(\ln|P-Q|\) for two random points of the unit square is \(\frac{\pi+\ln 2}{3}-\frac{25}{12}\).

Check. The double integral of Step 1 was computed by composite Gauss-Legendre quadrature (60 and 80 nodes) on dyadically graded panels refined towards the corner \(s=t=0\) (40 and 50 levels); both agree with the closed form to within \(10^{-15}\). A Monte Carlo estimate of the original fourfold integral with \(10^{7}\) samples gave \(-1.61057\pm0.00040\), within one standard error.

We show that the solutions are exactly the functions

with real constants \(a,b,c,d,e,g,h,k,m\): a nine-dimensional space, which we call \(V\). The analogy with squares (where the alternating vertex sum vanishes only for \(a+bx+cy+d(x^{2}+y^{2})\)) suggests a four-dimensional answer, and testing a few special hexagons on polynomials of low degree reinforces that guess; both miss the cubic and quartic members such as \((x^{2}+y^{2})^{2}\). The right viewpoint is complex.

Step 1: two triangles. Identify the plane with \(\mathbb{C}\) and let \(\zeta=e^{2\pi i/3}\). A regular hexagon \(P_{1}\dots P_{6}\) with centre \(c\) has \(P_{k+3}-c=-(P_{k}-c)\), and \(P_{1},P_{3},P_{5}\) form an equilateral triangle, so \(\{P_{1},P_{3},P_{5}\}=\{c+\zeta^{j}u\}\) and \(\{P_{2},P_{4},P_{6}\}=\{c-\zeta^{j}u\}\) \((j=0,1,2)\), where \(u=P_{1}-c\ne0\). Conversely every \(c\) and every \(u\ne0\) arise. Writing

the condition reads \(T_{c}(-u)-T_{c}(u)=0\): for every \(c\), the function \(T_{c}\) is even. The condition is linear in \(f\) and invariant under translations.

Step 2: every member of \(V\) works. For integers \(p,q\ge0\),

and the inner sum is \(3\) if \(i\equiv l\pmod 3\) and \(0\) otherwise. If \(p,q\le2\) then \(0\le i,l\le2\), so only \(i=l\) survives, every surviving term has even degree \(2i\) in \(u\), and the sum is even in \(u\). The real and imaginary parts of \(z^{p}\bar z^{q}\) with \(p,q\le2\) are \(1,x,y,x^{2}-y^{2},xy,x^{2}+y^{2},x(x^{2}+y^{2}),y(x^{2}+y^{2}),(x^{2}+y^{2})^{2}\), which span \(V\). For instance \(\sum_{j}|c+\zeta^{j}u|^{4}=3(|c|^{2}+|u|^{2})^{2}+6|c|^{2}|u|^{2}\).

Step 3: smooth solutions lie in \(V\). Let \(f\) be a \(C^{\infty}\) solution, and use the Wirtinger derivatives \(\partial=\frac12(\partial_{x}-i\partial_{y})\), \(\bar\partial=\frac12(\partial_{x}+i\partial_{y})\); since \(f\) is real, \(\bar\partial^{\,3}f=\overline{\partial^{3}f}\). Fix \(c,u\) and put \(h(t)=T_{c}(tu)\) for real \(t\). Then \(h(t)=h(-t)\) for \(t\ne0\), hence for all \(t\) by continuity, so \(h'''(0)=0\). For a direction \(v\in\mathbb{C}\),

Summing over \(v=\zeta^{j}u\) multiplies the \((i,l)\) term by \(\sum_{j}\zeta^{j(i-l)}\); as \(i-l\in\{3,1,-1,-3\}\), only \((3,0)\) and \((0,3)\) survive, and

Since \(u^{3}\) runs over all of \(\mathbb{C}\), \(\partial^{3}f\equiv0\) and hence \(\bar\partial^{\,3}f\equiv0\). Every operator \(\partial^{a}\bar\partial^{\,b}\) with \(a+b=6\) has \(a\ge3\) or \(b\ge3\), so it kills \(f\); as \(\partial_{x}=\partial+\bar\partial\) and \(\partial_{y}=i(\partial-\bar\partial)\), all sixth-order partial derivatives of \(f\) vanish and \(f\) is a polynomial of degree at most \(5\). Write it uniquely as \(f=\sum a_{pq}z^{p}\bar z^{q}\). Since \(\partial^{3}(z^{p}\bar z^{q})=p(p-1)(p-2)z^{p-3}\bar z^{q}\) sends distinct monomials to distinct monomials, \(\partial^{3}f=0\) forces \(a_{pq}=0\) for \(p\ge3\), and likewise \(\bar\partial^{\,3}f=0\) forces \(a_{pq}=0\) for \(q\ge3\). So \(f\in V\).

Step 4: from continuous to smooth. Let \(f\) be a continuous solution and \(\varphi_{\varepsilon}\ge0\) a smooth bump supported in the disc of radius \(\varepsilon\) with integral \(1\). The function \(f_{\varepsilon}(z)=\int f(z-y)\varphi_{\varepsilon}(y)\,dy\) is smooth, and since every translate \(P_{k}-y\) of a regular hexagon is a regular hexagon,

By Step 3, \(f_{\varepsilon}\in V\). As \(\varepsilon\to0\), \(f_{\varepsilon}\to f\) uniformly on every closed disc \(K\) (uniform continuity). Restriction to \(K\) is injective on \(V\) (a polynomial vanishing on an open set is zero), so \(V|_{K}\) is a finite-dimensional, hence closed, subspace of \(C(K)\), and \(f|_{K}\in V|_{K}\). The representing element of \(V\) is unique, so it is the same for all \(K\), and \(f\in V\).

Conclusion. The continuous solutions are exactly the nine-parameter family above. Beyond the affine functions and \(x^{2},xy,y^{2}\), the functions \(x(x^{2}+y^{2})\), \(y(x^{2}+y^{2})\) and \((x^{2}+y^{2})^{2}\) also satisfy the condition.

Check. By computer, the condition was imposed at 400 random regular hexagons (both orientations) on the 36-dimensional space of real polynomials of degree at most 7: the solution space has dimension exactly 9, the nine basis functions above satisfy the condition at 2000 further random hexagons to within \(2\times10^{-11}\) (rounding), and test functions such as \(x^{3}\), \(x^{2}y^{2}\) and \((x^{2}+y^{2})x^{2}\) fail.