What is the minimum over all distributions of unit vectors of the variance of the dot product of the vectors?
I will present an equivalent but simpler-looking formulation of the problem, and show a lower bound of (n/k − 1) / (n−1).  I also show a connection to an open problem in quantum information.  [Edit in revision 3: In earlier revisions, I claimed that an exact characterization of the cases in which the lower bound shown below is attained is likely to be difficult because an analogous question in the complex case includes an open problem about SIC-POVMs in quantum information.  However, this connection to SIC-POVMs was incorrect.  For details, see the section “Incorrect connection to SIC-POVMs in quantum information” below.]

Equivalent formulation

First, as was already pointed out in daniello’s answer, note that Var(xiTxj) = E[(xiTxj)2] − E[xiTxj]2 = E[(xiTxj)2].  So in the rest of the answer, we forget about variance and instead minimize maxi≠jE[(xiTxj)2].

Next, once we decide our goal is to minimize maxi≠jE[(xiTxj)2], we can ignore the constraint that E[xiTxj] = 0.  This is because if we have random unit vectors x1, …, xn, then we can negate each of them independently with probability 1/2 to satisfy E[xiTxj] = 0 without changing the value of the objective function maxi≠jE[(xiTxj)2].

Moreover, changing the objective function from maxi≠jE[(xiTxj)2] to (1/(n(n−1)))∑i≠jE[(xiTxj)2] does not change the optimum value.  The latter is at most the former because the average is at most the maximum.  However, we can always make the values of E[(xiTxj)2] for different choices of (i, j) (i≠j) equal by permuting the n vectors x1, …, xn randomly.

So for any n and k, the optimum value of the problem in question is equal to the minimum of (1/(n(n−1)))∑i≠jE[(xiTxj)2] where x1, …, xn are random variables which take unit vectors in ℝk as values.

However, by linearity of expectation, this objective function is equal to the expected value E[(1/(n(n−1)))∑i≠j(xiTxj)2].  Because the minimum is at most the average, there is no need to consider probability distributions anymore.  That is, the optimum value of the problem above is equal to the optimum value of the following:


  Choose unit vectors x1, …, xn ∈ ℝk to minimize (1/(n(n−1)))∑i≠j(xiTxj)2.


Lower bound

Using this equivalent formulation, we will prove that the optimal value is at least (n/k − 1) / (n−1).

For 1≤i≤n, let Xi=xixiT be the rank-1 projector corresponding to the unit vector xi.  Then, it holds that (xiTxj)2 = Tr(XiXj).

Let Y = ∑iXi.  Then, it holds that ∑i≠jTr(XiXj) = ∑i,jTr(XiXj) − n = Tr(Y2) − n.

The Cauchy–Schwarz inequality implies that Tr(Y2) ≥ (Tr Y)2/k = n2/k, and therefore ∑i≠jTr(XiXj) = Tr(Y2) − n ≥ n2/k − n.  By dividing by n(n−1), we obtain that the objective value is at least (n/k − 1) / (n−1).

In particular, when n=k+1, daniello’s answer is within a factor of 2 from the optimal value.

When is this lower bound attainable?

Attaining this lower bound (n/k − 1) / (n−1) is equivalent to making Y=(n/k)I.  I do not know the exact characterization when it is attainable, but following sufficient conditions exist:


When n=k+1, it is attainable by considering k+1 unit vectors which form a regular k-simplex centered at the origin, improving from 2/(k(k+1)) in daniello’s answer to the optimal 1/k2.
When n is a multiple of k, it is clearly attainable by fixing an orthonormal basis of ℝk and assigning each of the basis vectors to n/k of v1, …, vn.
More generally than the last bullet point, if it is attainable with some choice of k and both n=n1 and n=n2, then it is also attainable for the same k and n=n1+n2.  In particular, it is attainable if n=ak+b where a and b are integers satisfying a≥b≥0.


Although I have not checked details, it seems that any spherical 2-design gives a solution attaining this lower bound.

Incorrect connection to SIC-POVMs in quantum information

In earlier revisions, I stated:


  I suspect that answering this completely is a difficult question.  The reason is that if we instead consider the complex vector space ℂk, this question is related to an open problem in quantum information.


But this relation was incorrect.  I will explain why.


  More precisely, consider the following problem:
  
  
    Choose unit vectors x1, …, xn ∈ ℂk to minimize (1/(n(n−1)))∑i≠j|xi*xj|2.
  
  
  The lower bound above equally holds in this complex version.   Consider the case where n=k2 in the complex version.  Then the lower bound is equal to 1/(k+1).


So far, it was correct.


  A set of k2 unit vectors x1, …, xk2 ∈ ℂk attaining the lower bound is called a SIC-POVM in dimension k,


This part was incorrect.  A SIC-POVM is a set of k2 unit vectors x1, …, xn ∈ ℂk for which |xi*xj|2 = 1/(k+1) for all i≠j.  Note that here the requirement must hold for all pairs i≠j, not just the average over all pairs i≠j.  In the “Equivalent formulation” section, we showed the equivalence between minimizing the maximum and minimizing the average, but this was possible because x1, …, xn were random variables taking unit vectors there.  Here x1, …, xn are just unit vectors, so we cannot use the same trick.