Precise probability that $m$ random vectors in $n$ dimensional space are nearly-orthogonal Consider $m$ vectors $v_1,\dots,v_m$ in $\mathbb R^n$, drawn uniformly and independetly from unit sphere. It is pretty straightforward from Chebyshev inequality that 
$$
\mathrm P (\forall i\ne j \ |v_i \cdot v_j|\leq \varepsilon) \to 1\  \text{as} \ n \to \infty.
$$ 
But what about quantitative version of this limit, i.e. if we define 
$$
f(m, \varepsilon, \delta) = \min\{n : \mathrm P(\forall i\ne j\ |v_i \cdot v_j| \leq \varepsilon)\geq 1 - \delta\} 
$$
what can we say about asymptotic behavior of $f$?
 A: $\newcommand{\ep}{\varepsilon}
\newcommand{\Ga}{\Gamma}
\newcommand{\de}{\delta}$
For $\ep\in(0,1)$, let
\begin{equation*}
 P_{m,n}:=P\Big(\bigcap_{1\le i<j\le m}\{|v_i\cdot v_j|\le\ep\}\Big)
 =1-Q_{m,n},
\end{equation*}
where 
\begin{equation*}
 Q_{m,n}:=P\Big(\bigcup_{1\le i<j\le m}\{|v_i\cdot v_j|>\ep\}\Big). 
\end{equation*}
By Bonferroni inequalities,
\begin{equation*}
 Mp\ge Q_{m,n}\ge Mp-R/2,
\end{equation*}
where 
\begin{equation*}
 p:=P(|v_1\cdot v_2|>\ep), 
\end{equation*}
\begin{equation*}
 M:=m(m-1)/2,
\end{equation*}
\begin{equation*}
 R:=\sum_{1\le i<j\le m}\;\sum_{1\le k<l\le m,\,(k,l)\ne(i,j)}
 P(|v_i\cdot v_j|>\ep,|v_k\cdot v_l|>\ep). 
\end{equation*}
If $\{i,j\}\cap\{k,l\}=\emptyset$, then $P(|v_i\cdot v_j|>\ep,|v_k\cdot v_l|>\ep)=P(|v_i\cdot v_j|>\ep)\,P(|v_k\cdot v_l|>\ep)=p^2$. 
If $\{i,j\}\cap\{k,l\}\ne\emptyset$ but $\{i,j\}\ne\{k,l\}$, then, using the iid condition on the $u_i$'s and the spherical symmetry, for (say) the unit vector $e_1$ of the standard orthonormal basis of $\mathbb R^n$, we have 
\begin{multline}
 P(|v_i\cdot v_j|>\ep,|v_k\cdot v_l|>\ep)
 =P(|v_1\cdot v_2|>\ep,|v_1\cdot v_3|>\ep) \\ 
 =P(|e_1\cdot v_2|>\ep,|e_1\cdot v_3|>\ep)
 =P(|e_1\cdot v_2|>\ep)\,P(|e_1\cdot v_3|>\ep)=p^2. 
\end{multline}
So, $P(|v_i\cdot v_j|>\ep,|v_k\cdot v_l|>\ep)=p^2$ for any $i,j,k,l$ such that $1\le i<j\le m,\,1\le k<l\le m,\,(k,l)\ne(i,j)$. So, 
\begin{equation*}
R=M(M-1)p^2.  
\end{equation*}
Next, 
\begin{equation*}
 p=P(|e_1\cdot v_1|>\ep)=K_nI_n,
\end{equation*}
where, with $n\to\infty$, 
\begin{equation*}
 K_n:=\frac{\Ga(n/2)}{\Ga(1/2)\Ga((n-1)/2)(n-1)^{1/2}}\to1/\sqrt\pi, 
\end{equation*}
\begin{equation*}
 I_n:=(n-1)^{1/2}\int_{\sqrt c}^\infty(1+t^2)^{-n/2}\,dt=e^{-nc/(2+o(1))}, 
\end{equation*}
\begin{equation*}
 c:=\frac{\ep^2}{1-\ep^2}. 
\end{equation*}
Collecting the pieces, we see that 
\begin{equation}
 Q_{m,n}=Me^{-nc/(2+o(1))}. 
\end{equation}
Setting now $\de=Q_{m,n}\to0$, we find the asymptotics of the needed $n$:
\begin{equation}
 n\sim2\frac{1-\ep^2}{\ep^2}\,\ln\frac{m(m-1)}{2\de}. 
\end{equation}
