Are there arbitrarily large families of lines in $\Bbb R^3$ with average angle $\ge \pi/3$? 
Question: Can I have an arbitrarily large finite family of lines $\ell_1,\dotsc,\ell_n\subset\Bbb R^3$ so that the average angle between two (distinct) lines is $\ge \pi/3$?

We can assume that all lines pass through the origin, and by angle I mean the smaller one of the two angles defined by two intersecting lines.
If my integral calculus hasn't failed me, then the expected angle between two random lines with uniformly chosen directions is exactly $1\, \mathrm{rad}$, wheras $\pi/3\approx 1.0472 >1$.
So generic arrangements won't do it. Can I be more clever? Or is there perhaps some standard measure theory argument that tells me that in the limit, uniformly chosen lines is the best I can do? If so, can we put a bound on $n$ for where this stops being possible?
I know from random experiments that $n=17$ is possible, but already for $n=18$ I have no examples so far.

Update
The comment by alesia answers the question as posed (and I am looking forward to accept their answer): there are arbitrarily large families of lines with average angle $> \pi/3$.
Though one should note that the bound $\pi/3$ still seems kind of substantial: the more lines I have, the close the average will be to $\pi/3$.
I won't change the question (I might post a separate one), but the following seems to be the more fundamental question:

Is there a probability distribution on lines in $\Bbb R^3$ so that the expected value for the angle between two (independently sampled) lines is $>\pi/3$?

I think No and I am looking forward to ideas. The value $\pi/3$ is attained by alesia's construction (uniformly
choosing directions from an orthogonal basis).
 A: $\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.
A: Too long for a comment and please check my arguments. I think that there may be gaps there. This is however in the spirit of @DavidESpeyer's comment:
Given the Gram matrix formed by the inner products of $M$ (real or complex valued) vectors of dimension $n,$
$$
G=[\langle x_i,x_j\rangle]_{M\times M}
$$
the Welch bound (see Wikipedia) can be written as
$$
\sum_{1\leq i,j\leq M} \lvert\langle x_i,x_j\rangle\rvert^2 \geq \frac{\mathbb{tr}(G)^2}{n}
$$
which becomes
$$
\sum_{i} \langle x_i,x_i \rangle^2+
\sum_{i\neq j} \lvert\langle x_i, x_j \rangle\rvert^2 \geq
\frac{(\sum_{i} \langle x_i, x_i \rangle)^2}{n}.
$$
Keep in mind that we can consider unit nonzero vectors, and that a vector and its opposite represent the same line. Continuing:
$$
M(1)^2+\sum_{i\neq j} \lvert\langle x_i, x_j \rangle\rvert^2 \geq
\frac{M^2}{n},
$$
or
$$
\sum_{i\neq j} \lvert\langle x_i, x_j \rangle\rvert^2 \geq
\frac{M(M-1)}{n},
$$
which when divided by the number of pairs of vectors ($M(M-1)$) gives
$$
\mathbb{E}\{\lvert\langle x_i,x_j\rangle\rvert^2\}
\geq \frac{M(M-1)}{M(M-1)}=\frac{1}{n}=\frac{1}{3},
$$
since we have $n=3$, being in 3 dimensional real space.
This gives
$$
\cos^2(\theta)\geq \frac{1}{3}
$$
for the minimal angle, or $\theta\geq 54.735\cdots$ degrees somewhat under one radian.
Now, since all these angles are equal, it would seem that a perturbation of one line in this arrangement would actually increase the average angle. Assuming this, argument above is independent of $M$ it hasn't managed to rule out an arbitrarily large set of lines with this value. But for all angles below this value, it seems that one cannot have an equal angle arrangement.
Independently, I just found out (wikipedia again) that there are upper bounds to the maximum number of equiangular lines of $\binom{n+1}{2}$ in $n$-dimensional Euclidean space, which gives 6 lines for this case but 12 if we also take the opposite of each line. This might mean that if the perturbation argument above can be made rigorous then the vectors pointing to the centres of the 12 faces of a regular icosahedron may be the best possible for the angle corresponding to those vectors.
A: Start with $n$ times each coordinate axis, so $3n$ lines in total. The average angle between non necessarily distinct lines is $1/3*0 + 2/3*\pi/2 = \pi/3$. So the average angle between "distinct" lines (here, distinct lines may coincide geometrically due to multiplicities) is greater than $\pi/3$.
You can then perturb the lines to make them geometrically distinct while keeping an average angle larger than $\pi/3$.
Although that's not the question, I'm not sure what happens if one wants an infinite number of lines.
