It is well-known that the expected distance between two uniformly random points on (the perimeter of) a unit circle is $4/\pi$. Here I would like to ask about the three-dimensional case:
What is the expected area of a (planar) circle through three uniformly random points on the surface of a unit sphere?
I conjecture that the answer is $\sqrt{2\pi}$ (see the "UPDATE" below for an explanation).
My attempt
Without loss of generality, assume that the first random point is $(1,0,0)$.
From sphere point picking, the other two random points' coordinates are, for $i=1,2$:
$x_i=\sqrt{1-(2u_i-1)^2}\cos{(2\pi\theta_i)}$
$y_i=\sqrt{1-(2u_i-1)^2}\sin{(2\pi\theta_i)}$
$z_i=2u_i-1$
where $u_i$ and $\theta_i$ are independent uniformly random real numbers between $0$ and $1$.
The normal to the plane that passes through the three points, is
$\newcommand\mycolv[1]{\begin{bmatrix}#1\end{bmatrix}}$
$\textbf{n}=\mycolv{y_1z_2-y_2z_1\\z_1(x_2-1)-z_2(x_1-1)\\y_2(x_1-1)-y_1(x_2-1)}$
The distance from the origin to the plane that passes through the three points, is
$h=\dfrac{|y_1z_2-y_2z_1|}{\left|\textbf{n}\right|}$
The area of the circle is $\pi(1-h^2)$.
If I'm not mistaken, the expected area of the circle is $\int_0^1 \int_0^1 \int_0^1 \int_0^1 \pi(1-h^2) du_1 du_2 d\theta_1 d\theta_2$. But my computer will not calculate this.
UPDATE
Using Desmos, if I let $u_1$ and $u_2$ equal specific numbers from $0$ to $1$, Desmos calculates$\int_0^1 \int_0^1 \pi(1-h^2)d\theta_1 d\theta_2$.
So I methodically calculated this integral for all the combinations of $u_1=0,0.1,0.2,…,1$ and $u_2=0,0.1,0.2,…,1$, and took the average of the values of the integral. The result is approximately $2.507$.
Could the true answer be as elegant as $\sqrt{2\pi}\approx 2.507$ ?
Best Answer
The answer turns out to be
$$ \frac{4\pi}{5} \approx 2.513274123. $$
Below is a simulation of $10^7$ trials using Mathematica 13, demonstrating a sample mean that is quite close to the exact value above:
Proof. Let $A, B, C$ be independent, uniformly sampled points on the unit sphere $\mathbb{S}^2 $ centered at the origin.
Step 1. It is easy to check that the distance $L$ between the origin and the line $\overline{AB}$ has the PDF of the form
$$ f_L(l) = \begin{cases} 2l, & 0 < l < 1, \\ 0, & \text{elsewhere}. \end{cases} $$
Indeed, assuming WLOG that $A = (0, 0, 1)$ is the north pole,
$$ \mathbf{P}(L \leq l) = \mathbf{P}( \text{[$z$-coordinate of $B$]} \leq 2l^2 - 1) = l^2 $$
and hence the desired claim follows.
Step 2. Now we condition on the event $\{L = l\}$. In doing so, we may assume that the line $\overline{AB}$ lies in the $xy$-plane, passes through the point $(l, 0, 0)$, and is perpendicular to the $x$-axis. In other words, we may assume that
$$ A = (l, \sqrt{1-l^2}, 0) \qquad\text{and}\qquad B = (l, -\sqrt{1-l^2}, 0). $$
Now, given the value of $l$, we introduce the following parametrization of $\mathbb{S}^2$:
$$ \Psi_l(r, \theta) = r \left(\sin\alpha, 0, \cos\alpha\right) + \sqrt{1-r^2} \bigl[ (\cos\alpha, 0, -\sin\alpha) \cos\theta + (0, 1, 0) \sin\theta \bigr] $$
Here, $-l \leq r \leq l$ and $-\pi \leq \theta \leq \pi$, and $\alpha = \arcsin(r/l)$. This parametrization is designed so that the following property holds:
Below is an animated example of this circle where $l = 2/3$ and $r$ varies from $-l$ to $l$:
Then a bit of tedious computation shows that the area element is computed as
$$ \mathrm{d}A = \left\| \frac{\partial \Psi_l}{\partial r} \times \frac{\partial \Psi_l}{\partial \theta} \right\| \, \mathrm{d}r\mathrm{d}\theta = \left|1 - \sqrt{\frac{1-r^2}{l^2-r^2}}\cos\theta\right| \, \mathrm{d}r\mathrm{d}\theta. $$
So, if $R$ denotes the distance between the origin and the plane spanned by the points $A, B, C$, then
\begin{align*} \mathbf{P}(R \leq \rho \mid L = l) &= \frac{1}{4\pi} \int_{\operatorname{Dom}(\Phi_l)} \mathbf{1}_{\{|r| \leq \rho\}} \, \mathrm{d}A \\ &= \frac{1}{4\pi} \int_{-l}^{l} \int_{-\pi}^{\pi} \left|1 - \sqrt{\frac{1-r^2}{l^2-r^2}}\cos\theta\right| \mathbf{1}_{\{|r| \leq \rho\}} \, \mathrm{d}\theta\mathrm{d}r. \end{align*}
Computing the inner integral by invoking the identity
$$ \int_{-\pi}^{\pi} |1 - k \cos\theta| \, \mathrm{d}\theta = 4\left( \frac{\pi}{2} - \arccos\left(\frac{1}{k}\right) + \sqrt{k^2 - 1} \right), $$
which holds for $k \geq 1$, the conditional CDF is recast as
\begin{align*} \mathbf{P}(R \leq \rho \mid L = l) &= \int_{-l}^{l} \frac{1}{\pi} \left( \frac{\pi}{2} - \arccos\sqrt{\frac{l^2-r^2}{1-r^2}} + \sqrt{\frac{1-l^2}{l^2-r^2}} \right) \mathbf{1}_{\{|r| \leq \rho\}} \, \mathrm{d}r. \end{align*}
From this, the conditional PDF of $R$ given $L$ is computed as
$$ f_{R \mid L} (r \mid l) = \frac{2}{\pi} \left( \frac{\pi}{2} - \arccos\sqrt{\frac{l^2-r^2}{1-r^2}} + \sqrt{\frac{1-l^2}{l^2-r^2}} \right) \mathbf{1}_{\{r \leq l\}}. $$
Then by a tedious calculation (combined with integration by parts), we can check that the distribution of $R$ is fairly elementary, with the PDF given by
$$ f_{R}(r) = \frac{3}{2}( 1 - r^2 ). $$
Therefore the desired answer is
$$ \mathbf{E}[\pi(1-R^2)] = \frac{3\pi}{2} \int_{0}^{1} (1 - r^2)^2 \, \mathrm{d}r = \boxed{\frac{4\pi}{5}} \approx 2.513274123. $$