$\def\RR{\mathbb{R}}$I have been working on and off on the problem where we sample ordered pairs $(x,y)$ of lines, with duplicates allowed. To make sure the problem is clearly understood, I'll restate it: For $\ell_1$ and $\ell_2$ two lines through the origin in $\RR^3$, define $d(\ell_1, \ell_2)$ to be the acute angle that they make. In other words, if $v_1$ and $v_2$ are unit vectors on $\ell_1$ and $\ell_2$, then $d(\ell_1, \ell_2)$ is $\cos^{-1} |v_1 \cdot v_2|$.
The problem now is that we have $n$ lines $\ell_1$, $\ell_2$, ..., $\ell_n$. We sample two lines $\ell_i$ and $\ell_j$ independently at random from these $n$. The problem is to show that the expected value of $d(\ell_i, \ell_j)$ is at most $\pi/3$. As has been observed many times, equality occurs when we have three mutually orthogonal lines: $\tfrac{3}{9} \cdot 0 + \tfrac{6}{9} \cdot \tfrac{\pi}{2} = \tfrac{\pi}{3}$. If we can handle the discrete case, we can also handle more general probability distributions, but I find the discrete case easier to think about.
Here are the results I have gotten:
Proposition 1 For $n \geq 3$, in any local maximum, any line in the configuration is perpendicular to at least two other lines in the configuration.
Proposition 2 The result is true for $n \leq 4$ lines.
Proposition 3 The result is true for families of $3n$ lines $(u_1, v_1, w_1, u_2, v_2, w_2, \ldots, u_n, v_n, w_n)$ where each $(u_i, v_i, w_i)$ is orthogonal.
It looks like it might be true that Proposition 1 could be strengthened to
Conjecture 4 For $n \geq 3$, any local maximum contains three orthogonal lines.
Proposition 5 If we can show Conjecture 4, the result follows.
Here are the proofs. We first need:
Lemma 6 Consider $d(x,y)$ as a function of $x \in S^2$, with $y$ fixed. Away from $\pm y$ and $y^{\perp}$, this is a smooth function and its Hessian (using the standard Riemannian metric on $S^2$) is positive semidefinite.
Proof Normalize $y$ to be the north pole and use polar coordinates $(\theta, \phi)$ for $x$, so $d(x,y) = \min(\phi, \pi - \phi)$. This is clearly smooth away from the equator $\phi = \pi/2$ and from the poles (where spherical coordinates aren't defined). I'll work in the northern hemisphere, where $\phi < \pi/2$, so $d(x,y) = \phi$.
Thanks to user3658307, I know the formula for the spherical Hessian. All the second derivatives vanish, as does the derivative with respect to $\theta$, and the Hessian is just
$$\begin{bmatrix} \cos \phi \sin \phi \tfrac{\partial \phi}{\partial \phi} &0 \\ 0&0 \end{bmatrix} =
\begin{bmatrix} \cos \phi \sin \phi & 0 \\ 0&0 \\ \end{bmatrix}.$$
Since $0 < \phi < \pi/2$, we have $\cos \phi \sin \phi>0$ and this is positive semidefinite. $\square$
Corollary 7 Let $x_1$, $x_2$, ..., $x_n$ be points on $S^2$. Suppose that $x_n$ is not orthogonal to and of the $x_i$ or equal to $\pm x_i$. Then the Hessian of $\sum_{i,j} d(x_i, x_j)$ with respect to $x_n$ is positive semidefinite.
Proof The sum of positive semidefinite matrices is positive semidefinite. $\square$
Proof sketch for Proposition 1 Using Corollary 7 and the second derivative test on Riemannian manifolds, we see that a local maximum cannot have any point $x_n$ as in the hypotheses of Corollary 7. Furthermore, taking some of the $x_i$ equal to $\pm x_n$ only makes things worse, since $d(x, y)$ is convex up for $x$ near $\pm y$. (I don't know the right language to say "convex up" for a non-smooth function on a Riemannian manifold, but I'm pretty confident that I'm right.) Even taking $x_n$ orthogonal to one $x_i$ doesn't help, since we could then wiggle $x_n$ along $x_i^{\perp}$ and $\sum_{j \neq i,n} d(x_j, x_n)$ will be convex up there by the Corollary. So, in any local maximum, each $x_n$ must be orthogonal to at least two distinct $x_i$. $\square$
Proof of Proposition 2: For $n \leq 3$ lines, the result is clear.
So we turn to the case of $4$ lines. Let $\Gamma$ be the graph with vertices $1$, $2$, $3$, $4$ and an edge $(i,j)$ if $\ell_i \perp \ell_j$. Then Proposition 1 says that every vertex of $\Gamma$ has degree $\geq 2$. Every such graph on $4$ vertices contains a $4$-cycle; assume WLOG that $\ell_1 \perp \ell_2$, $\ell_2 \perp \ell_3$, $\ell_3 \perp \ell_4$ and $\ell_4 \perp \ell_1$.
Then $\text{Span}(\ell_1, \ell_3) \perp \text{Span}(\ell_2, \ell_4)$. Since we are in $3$-dimensional space, one of these two spaces must be one-dimensional, so either $\ell_1 = \ell_3$ or $\ell_2 = \ell_4$, WLOG $\ell_2 = \ell_4$. Then either Prop 3 applied to $\ell_3$, or a direct computation, shows that the optimal choice is $\ell_1 \perp \ell_3$. $\square$
Lemma 8 Let $\ell_1$, $\ell_2$, $\ell_3$ be orthogonal and let $\ell$ be any other line. Then $\sum_{i=1}^3 d(\ell_i, \ell) \leq \pi$.
Proof This either follows directly from Prop 1 or from Prop 2. $\square$.
Proof of Prop 3 We break up our expected value into terms where both lines are chosen from the same triple $(u_i, v_i, w_i)$, and terms where the lines come from distinct triples.
When both lines come from the same triple, the expected value is $(3/9)(0)+(6/9)(\pi/2) = \pi/3$.
For two different triples $(u_i, v_i, w_i)$ and $(u_j, v_j, w_j)$, we have $9$ summands. Break them into three sums of $3$ summands each, and bound each of them above by Lemma 8. $\square$.
I note that every equality case that we know is of the form in Proposition 3 (which is why I started considering it).
Finally, I'll prove Proposition 5: Our proof is by induction on $n$, with the base cases $n \leq 2$ clear. If $n \geq 3$, let $(\ell_1, \ell_2, \ell_3)$ be mutually orthogonal. Then break our expected value of $d(\ell_i, \ell_j)$ into four sums, depending on whether $i$ and $j$ are $\leq 3$ or $\geq 4$.
If both are $\leq 3$, the expected value is $\pi/3$.
If $i\leq 3$ and $j \geq 4$, fix $j$ and bound the sum on $i$ by Lemma 8. If $i \geq 4$ and $j \leq 3$, do the same thing with the roles of $i$ and $j$ fixed.
Finally, if $i$, $j \geq 4$, we have the result by induction. $\square$
Some more thoughts: Define $\Gamma$ as in the proof of Prop 2. If $\Gamma$ has a triangle, then we can split it off as in the proof of Prop 5 and induct. If $\Gamma$ has a $4$-cycle, then we can use the $\text{Span}(\ell_1, \ell_3) \perp \text{Span}(\ell_2, \ell_4)$ trick in the proof of Prop 2, although I'm not actually sure how we finish the argument from there.
If we can partition the vertices of $\Gamma$ into two classes $I$ and $J$ with $\leq 2$ edges between vertices of different classes, then I think I have one more trick: Apply a rotation matrix $g \in SO(3)$ to $\ell_j$ for $j \in J$ while fixing $\ell_i$ for $i \in I$; since $SO(3)$ is $3$-dimensional and there are only two summands which are convex down, I think you can use a computation like the proof of Prop 1 to show that this cannot be a local maximum. But I am not good enough at Riemannian geometry to be confident in this.
The first graph I know for which none of these tricks apply (i.e. all vertices degree $\geq 2$, no triangles or quadrilaterals, and $3$-edge connected) is the Petersen graph.