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$$?

$$\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 where $$\begin{equation*} Q_{m,n}:=P\Big(\bigcup_{1\le i\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\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. 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 $$$$Q_{m,n}=Me^{-nc/(2+o(1))}.$$$$ Setting now $$\de=Q_{m,n}\to0$$, we find the asymptotics of the needed $$n$$: $$$$n\sim2\frac{1-\ep^2}{\ep^2}\,\ln\frac{m(m-1)}{2\de}.$$$$