Are almost all $k$-tuples of vectors in high dimensional space almost orthogonal? I would like to know, in high dimensions, if you pick $k$ vectors (not necessarily distinct) are they likely to be $\textit{almost}$ orthogonal. I make this precise below but I am also open to other interpretation of the question. 
Since we only care about orthogonality, we can assume the vectors are unit vectors. Thus I want to consider a set $\textbf{V}=\{x_1,...,x_k\}\subset S^n$, where $S^n$ is the unit sphere. We can thus view $\textbf{V}\in \underbrace{S^n\times\cdots\times S^n}_{k-times}=\mathcal{S}_{n,k}$. Let $\mathcal{O}\subset\mathcal{S}_{n,k}$ denote the subset of $k$-tuples of orthogonal vectors. Let $\mathcal{U}_\epsilon=\{x\in\mathcal{S}_{n,k} : d(x,\mathcal{O}) <\epsilon\}$ be the $\epsilon$ neighborhood of $\mathcal{O}$ in $\mathcal{S}_{n,k}$ (using, say, the $l^2$ product metric, I don't think this matters). Let $\mu$ be the obvious probability product measure on $\mathcal{S}_{n,k}$. My questions is then:
Is it true that $$\forall \epsilon > 0\textrm{ and } \forall k\textrm{ }\lim_{n\rightarrow\infty}\mu(\mathcal{U}_\epsilon) = 1$$
If not what can we say about this limit?
 A: We will write $P:= \mu$ (for more probabilistic notation) and we will show that $\lim_{n \to \infty} P(U_{\epsilon}) = 1$.
Lemma 1: For all $\epsilon>0$ there exists some $\delta>0$ such that if $(x^1,...,x^k) \in \prod_1^k S^n$ then $\sum_{i<j} |\langle x^i,x^j \rangle|< \delta$ implies $d(x, \mathcal O)<\epsilon$.
Proof: This is a consequence of the Gram-Schmidt orthogonalization process. Specifically, we define $u^1:=x^1$, then $u^2:=x^2-\frac{\langle x^2,u^1\rangle}{|u^1|^2} u^1$, then $u^3=x^3 - \frac{\langle x^3,u^2\rangle}{|u^2|^2} u^2- \frac{\langle x^3, u^1 \rangle}{|u^1|^2} u^1$, and so on. Then let $e^i:=u^i/|u^i|$, so the $e^i$ form an orthonormal set which was constructed from the original $x^i$ using formulas which depend only on the values of the inner products $\langle x^i,x^j \rangle$. By carefully keeping track and using induction, one may show that $|x^j - u^j|$ (and thus $|x^j-e^j|$) can be made arbitrarily small by making the $|\langle x^i,x^j \rangle|$ very small. $\Box$
Lemma 2: If $U,V$ are independent and uniformly distributed on $S^{n-1}$ then $P(|\langle U, V \rangle|>\epsilon) \leq \frac{1}{n\epsilon^2}$
Proof: By conditioning on $V$ and using independence we may write $$P(|\langle U, V \rangle|>\epsilon) = \Bbb E\bigg[ P\big(|\langle U, V \rangle|>\epsilon \;\big| \;V\big)\bigg] = \int_{S^n} P(|\langle U, v \rangle|>\epsilon)\;\sigma(dv) = P(|\langle U, e_1 \rangle|>\epsilon)$$ where $\sigma$ is the uniform measure on $S^n$, and $e_1=(1,0,...,0)$. The final equality follows by rotational invariance of $\sigma$. Writing $U=(U_1,...,U_n)$ we see that the $U_j$ are identically distributed and $\sum U_j^2=1$. Thus $1 = E[\sum U_j^2] = nE[U_1^2]$, so that $E[U_1^2]=1/n$. Consequently, $$P(|\langle U, e_1 \rangle|>\epsilon) = P(|U_1|>\epsilon) \leq \frac{E[U_1^2]}{\epsilon^2} = \frac{1}{n\epsilon^2}$$ At the end we have used Chebyshev's inequality with the second moment. $\Box$
Proof of the original claim: By combining the two lemmas, we see that $$P(U_{\epsilon}^c) = P\big( d(\{X^1,...,X^k\}, \mathcal O)>\epsilon\big) \leq P \bigg( \sum_{i<j} |\langle X^i,X^j \rangle |> \delta \bigg) $$$$\leq \sum_{i<j} P\bigg(|\langle X^i,X^j \rangle| > \frac{\delta}{k(k-1)}\bigg) \leq\frac{k^3(k-1)^3}{n \delta^2} \stackrel{n\to \infty}{\longrightarrow} 0$$
Remarks: it was essential that $k$ is fixed here (though the final line shows that $k$ is allowed to depend on $n$, so long as it grows slowly enough). Also note that we only used pairwise independence in this argument (full mutual independence was not needed).
A: It is true, and it can be seen as a consequence of the concentration of measure on the sphere.
Let $P$ be the uniform probability measure on $\mathbb S^{n}$. Given a subset $A\subset \mathbb S^n$, call $A_t:=\{y\in\mathbb S^n:dist(y,A)<t\}$ its $t$-neighbourhood. Then:
Concentration of measure: If $P(A)\geq\frac12$, then $P(A_t^c)\leq 2e^{-nt^2/64}$. (the specific constants are not important here)
Apply this to two complementary emispheres to obtain that for any $\nu\in\mathbb S^n$ we have $P(|x\cdot \nu|\geq t)\leq 4e^{-nt^2/64}$, i.e. the measure concentrates near any equator when the dimension increases.
In particular it follows immediately that the result you ask for is true if $k=2$: once you have chosen the first vector $x_1$, the probability that the second you choose will be near the hyperplane $x_1^\perp$ goes to $1$ as $n\to\infty$ (by the rotational invariance of the measure on the sphere it doesn't matter which specific $x_1$ you chose in the first step).
To prove it for general $k$ use induction. I'll sketch an argument. The case $k=2$ (or $k=1$ if you want) is the base case and it's proved. Suppose now the statement is true for $k$. With probability close to $1$ the first $k$ vectors are $\epsilon$-orthogonal. Move them one by one until $x_i$ is close to $e_i$, the $i^{th}$ standard unit base vector (by the rotation invariance we can compute probabilities in this new position). From the concentration of measure
$$ P(\{|x\cdot e_1|\geq t\}\cup \ldots\cup\{|x\cdot e_k|\geq t\})\leq 4ke^{-nt^2/64}$$
Therefore with probability close to $1$ the $(k+1)^{th}$ vector has small scalar products with the previous $x_i$'s. This implies that with probability close to $1$ the family is $\epsilon$-orthogonal.
